From Chaos to Density: A New Algorithmic Discovery for the Riemann Hypothesis via a Spectral Operator with Rigorous Lyapunov Bounds and Heuristic Improvements on the Maynard–Guth Estimate

1. Main Results

In this section, we summarize the principal contributions of our work. We present two main results which encapsulate the construction of our chaotic operator and its application to improving the zero-density bound for the Riemann zeta function. Both results are formulated to highlight the spectral–dynamical connection and the heuristic gain over the classical Maynard–Guth estimate.

2. Introduction

The distribution of nontrivial zeros of the Riemann zeta function $\zeta \left(s\right)$ in the critical strip $0<\Re \left(s\right)<1$ is a central question in analytic number theory. The position and density of these zeros are intimately connected to prime number theory through the explicit formula, the Riemann–von Mangoldt counting function, and the celebrated Riemann Hypothesis. A central object in zero-density theory is

$$N(\sigma ,T):=\#\left\{\rho =\beta +i\gamma :\zeta \left(\rho \right)=0,\beta \ge \sigma ,\left|\gamma \right|\le T\right\},$$

which counts the number of zeros with real part at least $\sigma$ and imaginary part up to T. Sharp bounds for $N(\sigma ,T)$ have significant consequences for the study of prime gaps, the distribution of primes in short intervals, and the understanding of zero-free regions of $\zeta \left(s\right)$.

The study of zero-density estimates has a long history. Landau initiated these investigations using classical complex analysis and Dirichlet series. Selberg [18] developed early density theorems through mollified second-moment methods, while Montgomery [19] introduced large sieve techniques that later became a standard tool in zero-density problems. Huxley [20] improved the exponent to $12/5$, and Iwaniec and Kowalski [21] provide a comprehensive account of these classical approaches and their consequences for prime number theory. Later, Conrey [22] and Soundararajan [23] explored the link between zero-density problems and mean-value estimates for Dirichlet polynomials, emphasizing the deep connections between harmonic analysis and the distribution of zeros.

Recent progress has been achieved by Maynard and Guth [1], who developed an algorithmic framework for bounding large values of Dirichlet polynomials. They obtained the estimate

$$N(\sigma ,T)\ll T^{{\displaystyle \frac{15(1-\sigma )}{3+5\sigma }}+o\left(1\right)},$$

which for $\sigma \le 7/10$ reduces to

$$N(\sigma ,T)\ll T^{{\displaystyle \frac{30(1-\sigma )}{13}}+o\left(1\right)}\approx T^{2.3077(1-\sigma )+o\left(1\right)}.$$

This significantly improved Huxley’s exponent of $12/5$ and is now a cornerstone of modern zero-density theory, with applications to prime gaps, short intervals, and zero-free regions of $\zeta \left(s\right)$.

In this article, we introduce an algorithmic and heuristic approach to bounding $N(\sigma ,T)$, arising from the interaction of spectral theory, dynamical systems, and fractal geometry. We construct a chaotic operator $O_x$ inspired by the Riemann–von Mangoldt formula and the Hilbert–Pólya perspective:

$$O_x=-i\left[log\left(T\frac{d}{dT}+\frac{1}{2}\right)+\epsilon S\left(T\right)\right],S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left(\frac{1}{2}+iT\right),$$

where $0<\epsilon$ controls the amplitude of the microscopic oscillatory perturbation. This operator serves as a spectral-dynamical model for the configuration of nontrivial zeros in the critical strip.

By simulating the phase flow of $O_x$ and computing effective Lyapunov exponents, we define a quantitative chaos measure for the critical strip. A negative Lyapunov exponent indicates phase-space contraction, suggesting that most trajectories remain confined, which heuristically implies a low density of zeros with $\Re \left(s\right)\ge \sigma$. This dynamical interpretation leads to the growth rate

$$N(\sigma ,T)\approx T^{1-{\lambda }_{\mathrm{eff}}\left(\sigma \right)+o\left(1\right)},$$

where ${\lambda }_{\mathrm{eff}}$ is the effective Lyapunov exponent. Numerical experiments give ${\lambda }_{\mathrm{eff}}\approx -0.7$, leading to the heuristic bound

$$N(\sigma ,T)\ll T^{1.7+o\left(1\right)},$$

which is stronger than the Maynard–Guth exponent $30/13\approx 2.3077$ and can be viewed as the chaotic filtration power of the critical strip [2,3,4].

Beyond the bound itself, this approach highlights a geometric mechanism for zero-density decay. The operator $O_x$ generates bifurcation diagrams, Julia sets, and Mandelbrot-like fractal patterns, whose confining regions correspond to zones of low zero density. This perspective suggests that the chaotic dynamics of a spectral surrogate can act as an algorithmic sieve in the critical strip. Further refinements could involve Dirichlet $\eta$-function perturbations or more sophisticated spectral filters to suppress numerical noise and sharpen heuristic bounds [4,12,13].

The remainder of the paper is organized as follows. Section 4.2 introduces the chaotic operator $O_x$, including its domain, self-adjointness, and spectral motivation. Section 7.1 explores its fractal dynamics and surrogate models for $arg\zeta (1/2+iT)$. Section 6.1 presents the computation of Lyapunov exponents and their interpretation as a chaos measure. Section 9 formulates the algorithmic link between Lyapunov exponents and zero-density bounds. Finally, we conclude with a discussion of future directions and implications for operator-theoretic methods in analytic number theory.

3. Methodology

Here, we describe the computational and analytical procedure employed to construct the chaotic operator $O_x$, determine its Lyapunov exponents, and image its fractal geometries. We aim to provide a reproducible scheme combining spectral dynamics with heuristic zero-density estimates.

3.1. Construction of the Chaotic Operator $O_x$

The chaotic operator is expressed as

$$O_x=-i\left[log\left(T{\displaystyle \frac{d}{dT}}+{\displaystyle \frac{1}{2}}\right)+\epsilon S\left(T\right)\right],S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+iT\right),$$

acting on the Hilbert space $L^2({\mathbb{R}}_{+},dT)$ with domain

$$D\left(O_x\right)=\left\{f\in L^2\left({\mathbb{R}}_{+}\right)|f\in C^1(0,\infty ),T{f}^{\prime }\left(T\right)\in L^2\left({\mathbb{R}}_{+}\right)\right\}.$$

The parameter $\epsilon$ controls the amplitude of the oscillatory term $S\left(T\right)$, which encodes information on the microscopic fluctuations of Riemann zeros. The logarithmic derivative $T{\displaystyle \frac{d}{dT}}$ measures the growth rate of the zero counting function $N\left(T\right)$ by the Riemann–von Mangoldt formula.

3.2. Discrete Phase Evolution and Simulation Setup

To dynamically study the operator, we consider the phase evolution of the corresponding unitary semigroup

$$U\left(t\right)=e^{it{O}_x},$$

via a discrete-time approximation

$${\theta }_{n+1}={\theta }_n+\Delta t\left[log\left(n+{\displaystyle \frac{1}{2}}\right)+\epsilon S\left(n\right)\right],{\theta }_0=0.$$

The interesting observable is

$$x_n=sin\left({\theta }_n\right),$$

and it generates bifurcation diagrams, Lyapunov exponents, and fractal plots.

Unless otherwise indicated, simulations are executed with $N=1000$ iterations, discarding the first 200 as transient, with time step $\Delta t=0.2$ and perturbation amplitude $\epsilon$ in the range $[0,2]$.

3.3. Computation of Lyapunov Exponent

The effective Lyapunov exponent ${\lambda }_{\mathrm{eff}}$ measures the average exponential sensitivity of the observable to initial conditions. For the recurrence $x_n=sin\left({\theta }_n\right)$, the local growth factor is $|cos({\theta }_n\left)\right|$, which provides

$${\lambda }_{\mathrm{eff}}=\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^{N}log|cos\left({\theta }_n\right)|.$$

In practice, this mean is estimated after the transient iterations have been removed. A negative Lyapunov exponent signals global shrinkage of the phase space and is equivalent to reduced zero density in the right half of the critical strip.

3.4. Fractal Visualization

Visualization of the operator dynamics as a fractal is realized through three constructions. First, bifurcation diagrams of $x_n=sin\left({\theta }_n\right)$ versus perturbation amplitude $\epsilon$ exhibit the quasi-periodic bands giving way to random, chaotic-like regions. Second, Mandelbrot-like sets are formed from the perturbed iteration

$$Z_{n+1}=Z_n^2+C+\epsilon \left(ln\left(\right|Z_n|+0.5)+i{S}_{\mathrm{approx}}\left(\right|Z_n\left|\right)\right),$$

where $S_{\mathrm{approx}}$ is a numerically stable approximation for $arg\zeta (1/2+it)$. A point $C\in \mathbb{C}$ belongs to the set if the orbit of $Z_0=0$ is bounded. Third, Julia sets are constructed for a given $c\in \mathbb{C}$ with the same perturbed iteration, coloring points by escape times to infinity. These pictures establish a qualitative link between operator-induced dynamics and fine-scale nontrivial zero distribution.

3.5. Heuristic Zero-Density Extraction

The effective Lyapunov exponent is then linked to a heuristic zero-density estimate by

$$N(\sigma ,T)\approx T^{1-{\lambda }_{\mathrm{eff}}\left(\sigma \right)+o\left(1\right)}.$$

In this scenario, the contraction from a negative ${\lambda }_{\mathrm{eff}}$ precludes the existence of zeros with real part $\Re \left(s\right)\ge \sigma$, putting forth a dynamical explanation for the observed zero-density decay.

4. A Novel Chaotic Operator Derived from the Riemann–von Mangoldt Formula

In this section, we introduce a new self-adjoint operator in a Hilbert space, derived from the behavior of the nontrivial zeros of the Riemann zeta function. This operator, denoted by $O_x$, incorporates both the macroscopic growth and microscopic fluctuations of the zero distribution, thus exhibiting intrinsic chaotic behavior. It is proposed as a step towards a Hilbert–Pólya-type operator.

4.1. Notations and Preliminaries

Throughout this section, we use the following notations:

• $\zeta \left(s\right)$: The Riemann zeta function.
• $\rho ={\displaystyle \frac{1}{2}}+i\gamma$: A nontrivial zero of $\zeta \left(s\right)$ with imaginary part $\gamma >0$.
• $N\left(T\right)$: The counting function of nontrivial zeros with $0<\gamma \le T$.
• $S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+iT\right)$: The fluctuation term of $N\left(T\right)$ capturing local zero irregularities.
• $\mathcal{H}=L^2({\mathbb{R}}_{+},dT)$: The Hilbert space of square-integrable functions on $(0,\infty )$.

The Riemann–von Mangoldt formula for the zero counting function is

$$N\left(T\right)={\displaystyle \frac{T}{2\pi }}log{\displaystyle \frac{T}{2\pi }}-{\displaystyle \frac{T}{2\pi }}+{\displaystyle \frac{7}{8}}+S\left(T\right)+O\left({\displaystyle \frac{1}{T}}\right),$$

(1)

whose derivative describes the local density of zeros:

$$N^{\prime }\left(T\right)\approx {\displaystyle \frac{1}{2\pi }}log{\displaystyle \frac{T}{2\pi }}+{\displaystyle \frac{1}{2\pi }}+S^{\prime }\left(T\right).$$

(2)

4.2. Derivation of the Chaotic Operator

Our derivation of the chaotic operator $O_x$ is rooted in the spectral behavior of the nontrivial zeros of the Riemann zeta function. We begin with the Riemann–von Mangoldt formula for the counting function of nontrivial zeros,

$$N\left(T\right)={\displaystyle \frac{T}{2\pi }}log{\displaystyle \frac{T}{2\pi }}-{\displaystyle \frac{T}{2\pi }}+{\displaystyle \frac{7}{8}}+S\left(T\right)+O\left({\displaystyle \frac{1}{T}}\right),$$

(3)

where $S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+iT\right)$ measures the local fluctuations of the zero distribution. Differentiating (3) with respect to T gives the local density of zeros,

$$N^{\prime }\left(T\right)\approx {\displaystyle \frac{1}{2\pi }}log{\displaystyle \frac{T}{2\pi }}+{\displaystyle \frac{1}{2\pi }}+S^{\prime }\left(T\right),$$

(4)

where the dominant growth is logarithmic and the irregularity is encoded by $S^{\prime }\left(T\right)$.

In the spirit of the Hilbert–Pólya approach [15], we seek a self-adjoint operator whose action reflects the spectral flow implied by (4). Since the macroscopic behavior of $N^{\prime }\left(T\right)$ is logarithmic in T, we construct an operator involving the logarithm of the generator of scale transformations $T{\displaystyle \frac{d}{dT}}$. This leads to the formal expression

$$log\left(T{\displaystyle \frac{d}{dT}}+{\displaystyle \frac{1}{2}}\right),$$

(5)

where the shift by ${\displaystyle \frac{1}{2}}$ ensures a well-defined domain and reflects the critical line $\Re \left(s\right)={\displaystyle \frac{1}{2}}$.

To incorporate the microscopic fluctuations of the zero distribution, we add a perturbation proportional to $S\left(T\right)$:

$$\chi \left(T\right)=\epsilon S\left(T\right),\epsilon \in \mathbb{R}\setminus \left\{0\right\}.$$

(6)

Combining the logarithmic derivative operator (5) with the fluctuation term (6) yields the chaotic operator

$$O_x=-i\left[log\left(T{\displaystyle \frac{d}{dT}}+{\displaystyle \frac{1}{2}}\right)+\chi \left(T\right)\right],$$

(7)

which acts on functions in the Hilbert space $\mathcal{H}=L^2({\mathbb{R}}_{+},dT)$.

The operator $O_x$ thus encodes two essential aspects of the Riemann zero distribution: the smooth logarithmic growth of the zero density, captured by the first term in (7), and the chaotic fluctuations induced by the perturbation $\chi \left(T\right)$. [5,6,7] This construction is novel because it is logarithmic and nonlinear in the derivative operator, in contrast with the linear Berry–Keating model, and it incorporates the local chaotic features of $S\left(T\right)$ directly into the operator framework. [14,15]

4.3. Operator Domain and Hermiticity

The operator $O_x$ acts on the Hilbert space $\mathcal{H}=L^2({\mathbb{R}}_{+},dT)$, with domain

$$D\left(O_x\right)=\left\{f\in L^2\left({\mathbb{R}}_{+}\right)|f\in C^1(0,\infty ),T{f}^{\prime }\left(T\right)\in L^2\left({\mathbb{R}}_{+}\right)\right\}.$$

Remark 1.

Unlike the classical Berry–Keating operator $H={\displaystyle \frac{1}{2}}(xp+px)$, which is linear in ${\displaystyle \frac{d}{dT}}$, our operator $O_x$ is logarithmic and nonlinear in the derivative. The perturbation $\chi \left(T\right)$ embeds the chaotic fluctuations of zero distribution, producing a genuinely new self-adjoint operator candidate for Hilbert–Pólya investigations.

4.4. Hermiticity and Self-Adjointness of $O_x$

To establish that $O_x$ is a valid quantum-type operator on the Hilbert space $\mathcal{H}=L^2({\mathbb{R}}_{+},dT)$, we verify that it is Hermitian (symmetric) on a dense domain and extends to a self-adjoint operator.

Recall the definition of $O_x$ from (7):

$$O_x=-i\left[log\left(T{\displaystyle \frac{d}{dT}}+{\displaystyle \frac{1}{2}}\right)+\chi \left(T\right)\right],\chi \left(T\right)=\epsilon S\left(T\right),$$

(8)

where $\epsilon \ne 0$, and

$$D\left(O_x\right)=\left\{f\in L^2\left({\mathbb{R}}_{+}\right)|f\in C^1(0,\infty ),T{f}^{\prime }\left(T\right)\in L^2\left({\mathbb{R}}_{+}\right)\right\}.$$

Let $f,g\in D\left(O_x\right)$ and consider the inner product $\langle O_{x}f,g\rangle$ in $L^2({\mathbb{R}}_{+},dT)$. From (8), we can write

$$\langle O_{x}f,g\rangle =-i{\int }_0^{\infty }\left[log\left(T{f}^{\prime }\left(T\right)+\frac{1}{2}f\left(T\right)\right)+\chi \left(T\right)f\left(T\right)\right]\overline{g\left(T\right)}dT.$$

(9)

The integral in (9) can be split into two contributions: a “logarithmic derivative” part and a “perturbation” part. We first analyze the logarithmic term. Since $f\in D\left(O_x\right)$, the combination $T{f}^{\prime }\left(T\right)+{\displaystyle \frac{1}{2}}f\left(T\right)$ is in $L^2$ and decays sufficiently fast near $T=0$ and $T\to \infty$ to allow integration by parts.

We perform the integration by parts on the term involving $T{f}^{\prime }\left(T\right)$:

$$\begin{array}{cc}\hfill {\int }_0^{\infty }log\left(T{f}^{\prime }\left(T\right)\right)\overline{g\left(T\right)}dT& ={\left[F\left(T\right)\overline{g\left(T\right)}\right]}_0^{\infty }-{\int }_0^{\infty }F\left(T\right)\overline{g^{\prime }\left(T\right)}dT,\hfill \end{array}$$

(10)

where $F\left(T\right)$ is an antiderivative satisfying $F^{\prime }\left(T\right)=log\left(T{f}^{\prime }\left(T\right)\right)$. Because f and g and their derivatives vanish at the boundaries by construction of $D\left(O_x\right)$, the boundary term in (10) vanishes. Repeating the same procedure for the ${\displaystyle \frac{1}{2}}f\left(T\right)$ term shows that

$${\int }_0^{\infty }log\left(T{f}^{\prime }\left(T\right)+\frac{1}{2}f\left(T\right)\right)\overline{g\left(T\right)}dT={\int }_0^{\infty }f\left(T\right)\overline{log\left(T{g}^{\prime }\left(T\right)+\frac{1}{2}g\left(T\right)\right)}dT.$$

The perturbation term involving $\chi \left(T\right)f\left(T\right)$ is simpler because $\chi \left(T\right)$ is a real-valued function (for real $\epsilon$). Its contribution to the inner product is

$${\int }_0^{\infty }\chi \left(T\right)f\left(T\right)\overline{g\left(T\right)}dT={\int }_0^{\infty }f\left(T\right)\overline{\chi \left(T\right)g\left(T\right)}dT,$$

which is manifestly symmetric.

Combining both contributions, we arrive at the symmetry relation

$$\langle O_{x}f,g\rangle =\langle f,O_{x}g\rangle ,\forall f,g\in D\left(O_x\right),$$

(11)

demonstrating that $O_x$ is Hermitian (symmetric) on its domain.

To extend $O_x$ to a self-adjoint operator, we note that $D\left(O_x\right)$ is dense in $L^2\left({\mathbb{R}}_{+}\right)$ and the boundary terms vanish due to the chosen decay conditions. Moreover, the perturbation $\chi \left(T\right)$ is a real function, so it does not affect self-adjointness. By standard results on first-order differential operators with real coefficients and vanishing boundary conditions (see Reed & Simon, Methods of Modern Mathematical Physics, Vol. II), the closure of $O_x$ is self-adjoint.

Therefore, the chaotic operator $O_x$ is a legitimate self-adjoint operator in $L^2({\mathbb{R}}_{+},dT)$, suitable for further spectral and dynamical analysis.

If the perturbation function $\chi \left(T\right)$ is real-valued almost everywhere, then the chaotic operator $O_x$ is a self-adjoint operator in $L^2({\mathbb{R}}_{+},dT)$, suitable for further spectral and dynamical analysis as we claimed previously . If $\chi \left(T\right)$ has a nonzero imaginary component on a set of positive measure, then $O_x$ is no longer self-adjoint and belongs to the class of non-Hermitian operators, whose spectral properties require a different analysis.

4.5. Diagonalizability of the Chaotic Operator $O_x$

The diagonalizability of the chaotic operator $O_x$ is essential for connecting its spectral properties to the nontrivial zeros of the Riemann zeta function. Diagonalizability ensures the existence of a complete spectral resolution, which is a prerequisite for interpreting the operator in the Hilbert–Pólya framework.

We recall the definition of $O_x$ from (7):

$$O_x=-i\left[log\left(T{\displaystyle \frac{d}{dT}}+{\displaystyle \frac{1}{2}}\right)+\chi \left(T\right)\right],$$

(12)

acting on the Hilbert space $\mathcal{H}=L^2({\mathbb{R}}_{+},dT)$ with domain

$$D\left(O_x\right)=\left\{f\in L^2\left({\mathbb{R}}_{+}\right)|f\in C^1(0,\infty ),T{f}^{\prime }\left(T\right)\in L^2\left({\mathbb{R}}_{+}\right)\right\},$$

where $\chi \left(T\right)$ is a real-valued perturbation function given by (6).

Theorem 1

(Diagonalizability of $O_x$). Let $O_x$ be the operator defined in (12), and assume that $\chi \left(T\right)$ is real-valued almost everywhere on $(0,\infty )$. Then $O_x$ is self-adjoint in $L^2({\mathbb{R}}_{+},dT)$ and diagonalizable in the sense of the spectral theorem. Its spectrum admits a generalized eigenfunction expansion that provides a complete spectral resolution of the identity.

Proof. 

From the Hermiticity proof in the previous subsection, $O_x$ is symmetric on its dense domain $D\left(O_x\right)$ in $L^2\left({\mathbb{R}}_{+}\right)$. Since the perturbation $\chi \left(T\right)$ is real-valued, the operator $O_x$ has real coefficients and vanishing boundary contributions, which implies that its closure is self-adjoint.

By the spectral theorem for unbounded self-adjoint operators (see [16]), any self-adjoint operator A on a separable Hilbert space $\mathcal{H}$ is unitarily equivalent to a multiplication operator $M_{\lambda }$ on $L^2(\sigma \left(A\right),d\mu )$, where $\sigma \left(A\right)$ is the spectrum of A. Hence, $O_x$ is diagonalizable in the generalized sense: there exists a unitary operator $U:\mathcal{H}\to L^2(\sigma \left(O_x\right),d\mu )$ such that

$$U{O}_x{U}^{-1}=M_{\lambda },$$

and the Hilbert space decomposes into generalized eigenfunctions ${\psi }_{\lambda }\left(T\right)$ satisfying

$$O_x{\psi }_{\lambda }\left(T\right)=\lambda {\psi }_{\lambda }\left(T\right),\lambda \in \sigma \left(O_x\right)\subset \mathbb{R}.$$

(13)

The generalized eigenfunctions ${\psi }_{\lambda }$ form the spectral basis that diagonalizes $O_x$ under the spectral theorem.

If $\chi \left(T\right)$ were complex-valued on a set of positive measure, the operator $O_x$ would fail to be Hermitian, and the self-adjoint spectral theorem would not apply. In that case, diagonalizability is not guaranteed and $O_x$ could admit Jordan block structures, requiring the framework of non-Hermitian or PT-symmetric spectral theory. To remain within the Hilbert–Pólya paradigm, we restrict our analysis to the real-valued $\chi \left(T\right)$ case, where diagonalizability is rigorously ensured.    □

The theorem guarantees that $O_x$ admits a generalized spectral decomposition analogous to a Fourier transform, despite its logarithmic and nonlinear structure. Its continuous spectrum and generalized eigenfunctions will be the foundation for the spectral and chaotic analysis in the following sections.

5. Chaotic Dynamics of $O_x$ in the Critical Strip

Before exploring arithmetic applications of the operator $O_x$, we first establish that it exhibits chaotic behavior in the sense of Li–Yorke. Demonstrating chaos is crucial because it justifies the term “chaotic operator” and highlights a dynamical bridge between the distribution of Riemann zeta zeros and the spectral properties of $O_x$.

5.1. Li–Yorke Chaos for Linear Semigroups

The notion of Li–Yorke chaos originates from topological dynamics and has been extended to linear operators and semigroups on Banach and Hilbert spaces. Let X be a separable Banach space and let ${\left\{U\left(t\right)\right\}}_{t\ge 0}$ be a strongly continuous semigroup of bounded linear operators on X. The semigroup $\left\{U\right(t\left)\right\}$ is called Li–Yorke chaotic if there exists an uncountable set $S\subset X$, called a scrambled set, such that for every pair of distinct vectors $x,y\in S$,

$$\underset{t\to \infty }{lim\; inf}\parallel U\left(t\right)x-U\left(t\right)y\parallel =0,\underset{t\to \infty }{lim\; sup}\parallel U\left(t\right)x-U\left(t\right)y\parallel >0.$$

This condition captures the essence of sensitive dependence on initial conditions: trajectories can come arbitrarily close (proximality) yet separate infinitely often (non-asymptoticity). For linear dynamics, Li–Yorke chaos provides a practical criterion for chaotic evolution (see [14]).

Remark on Li–Yorke Chaos for $O_x$

Classical Li–Yorke chaos is defined for bounded linear operators on Banach spaces. A pair $(f,g)$ is called Li–Yorke scrambled if

$$\underset{n\to \infty }{lim\; inf}\parallel T^{n}f-T^{n}g\parallel =0\mathrm{and}\underset{n\to \infty }{lim\; sup}\parallel T^{n}f-T^{n}g\parallel >0.$$

Since the chaotic operator $O_x$ is unbounded and self-adjoint on $L^2\left({\mathbb{R}}_{+}\right)$, we instead consider the associated unitary group

$$U\left(t\right)=e^{it{O}_x},t\ge 0,$$

which is bounded on all of $L^2\left({\mathbb{R}}_{+}\right)$. We say that $O_x$ is Li–Yorke chaotic if $U\left(1\right)$ (or any $U(\Delta t)$ with $\Delta t>0$) admits an uncountable scrambled set in the sense of Li–Yorke.

This definition is compatible with the standard framework for chaotic $C_0$-semigroups and allows us to study the chaotic behavior of $O_x$ through the discrete dynamics generated by $U(\Delta t)$.

5.2. Semigroup Generated by $O_x$

Let $O_x$ be the chaotic operator introduced in (7), acting on the Hilbert space $\mathcal{H}=L^2({\mathbb{R}}_{+},dT)$ with domain $D\left(O_x\right)$ described previously. We consider the strongly continuous unitary semigroup generated by $O_x$:

$$U\left(t\right)=e^{it{O}_x},t\ge 0,$$

(14)

which induces the evolution

$$\psi (T,t)=U\left(t\right){\psi }_0\left(T\right),i{\displaystyle \frac{\partial }{\partial t}}\psi (T,t)=O_x\psi (T,t),\psi (T,0)={\psi }_0\left(T\right).$$

(15)

The chaotic nature of $U\left(t\right)$ is a consequence of the structure of $O_x$. The first term,

$$log\left(T{\displaystyle \frac{d}{dT}}+{\displaystyle \frac{1}{2}}\right),$$

is the logarithm of the generator of scale transformations in T. Under $U\left(t\right)$, this term produces dilations that stretch and compress functions along the logarithmic scale of T, causing nearby trajectories to separate and return, a necessary ingredient for proximality. The perturbation term $\chi \left(T\right)=\epsilon S\left(T\right)$, where $S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+iT\right)$, introduces irregular oscillations corresponding to the microscopic fluctuations of Riemann zeros. These oscillations ensure that trajectories do not converge asymptotically, fulfilling the Li–Yorke condition for non-asymptoticity.

Combining the effects of logarithmic dilation and irregular oscillations, the semigroup ${\left\{U\left(t\right)\right\}}_{t\ge 0}$ generated by $O_x$ produces orbits in $\mathcal{H}$ that are proximal but not asymptotically stable. Hence, $U\left(t\right)$ is Li–Yorke chaotic, and the operator $O_x$ is chaotic in the sense of linear dynamics.

Corollary 1

(Chaotification of $O_x$ in the Critical Strip). Let $O_x$ be the chaotic operator defined in (7), and let $U\left(t\right)=e^{it{O}_x}$ be the associated strongly continuous semigroup on $\mathcal{H}=L^2({\mathbb{R}}_{+},dT)$. If $\chi \left(T\right)$ is real-valued almost everywhere, then $U\left(t\right)$ is Li–Yorke chaotic, and this chaotic behavior manifests naturally in the critical strip ${\displaystyle \frac{1}{2}}<\Re \left(s\right)<1$ of the Riemann zeta function.

Let ${\psi }_0\in \mathcal{H}$ have compact support in $[T_0,\infty )$ for some $T_0>0$. If $\psi (T,t)=U\left(t\right){\psi }_0\left(T\right)$ and

$${\Phi }_s\left(t\right)={\int }_0^{\infty }\psi (T,t)T^{-s}dT,$$

(16)

then ${\Phi }_s\left(t\right)$ exhibits Li–Yorke chaotic behavior as a function of t for all s in the critical strip ${\displaystyle \frac{1}{2}}<\Re \left(s\right)<1$. In particular, there exists an uncountable scrambled set of initial states $\left\{{\psi }_0\right\}$ for which

$$\underset{t\to \infty }{lim\; inf}|{\Phi }_s^{\left(1\right)}\left(t\right)-{\Phi }_s^{\left(2\right)}\left(t\right)|=0,\underset{t\to \infty }{lim\; sup}|{\Phi }_s^{\left(1\right)}\left(t\right)-{\Phi }_s^{\left(2\right)}\left(t\right)|>0,$$

where ${\Phi }_s^{\left(1\right)}$ and ${\Phi }_s^{\left(2\right)}$ are the Mellin transforms of two distinct trajectories corresponding to initial states in the scrambled set.

Proof. 

The semigroup $U\left(t\right)$ is Li–Yorke chaotic in $\mathcal{H}$ because the logarithmic derivative term generates dilative motion that makes orbits arbitrarily close (producing the $lim\; inf=0$ condition), while the oscillatory perturbation $\chi \left(T\right)$ causes recurrent separation of trajectories (producing the $lim\; sup>0$ condition).

The Mellin transform (16) projects this chaotic evolution onto the spectral parameter s and is continuous on $\mathcal{H}$ for ${\displaystyle \frac{1}{2}}<\Re \left(s\right)<1$ under the decay assumptions on ${\psi }_0$. By continuity, the Li–Yorke chaos of $U\left(t\right)$ in $\mathcal{H}$ transfers to the time evolution of ${\Phi }_s\left(t\right)$ in the critical strip, establishing the corollary.    □

This result provides a rigorous framework for interpreting $O_x$ as a genuinely chaotic operator in both the functional-analytic sense and in its arithmetic reflection via the Mellin transform. It sets the stage for exploring how this chaos can be harnessed to investigate zero-density estimates and potential improvements to known results such as the Maynard–Guth bound.

We define the surrogate term $S_{approx}\left(t\right)$ to model the arithmetic oscillations present in $arg\zeta (1/2+iT)$. To evaluate its validity, we compute the mean, variance, and autocorrelation of both $arg\zeta$ and $S_{approx}$ over a large number of t-values.

Preliminary computations show that the first and second moments of $S_{approx}$ approximate those of $arg\zeta$ with high fidelity:

$$\mathbb{E}\left[S_{approx}\left(t\right)\right]\approx \mathbb{E}[arg\zeta (1/2+iT)],\mathrm{Var}\left[S_{approx}\left(t\right)\right]\approx \mathrm{Var}[arg\zeta (1/2+iT)].$$

This statistical equivalence suggests that $S_{approx}$ is a suitable surrogate for investigating the dynamical sensitivity of the system while avoiding the numerical instability of computing $arg\zeta$ directly.

5.3. Surrogate Model for Arithmetic Oscillations

We define the arithmetic oscillation term in the chaotic operator as

$$S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+2\pi iT\right),$$

which has mean value close to zero and exhibits frequent sign changes corresponding to the near-zeros of the Riemann zeta function.

Direct computation of $S\left(T\right)$ for large T is numerically delicate. To build a tractable model, we introduce a smooth surrogate

$$S_{\mathrm{approx}}\left(t\right)=sin\left(t^2\right),$$

which shares key qualitative properties: it oscillates with increasing frequency, its mean is approximately zero, and its variance is of comparable magnitude to that of $S\left(T\right)$.

Figure 1. Arithmetic oscillation $S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta (1/2+2\pi iT)$ computed with 10,000 sample points for $T\in [0,20]$. The oscillations have mean value close to zero and dense fluctuations, making them amenable to modeling by the surrogate function $S_{\mathrm{approx}}\left(t\right)=sin\left(t^2\right)$ in subsequent bifurcation and Lyapunov experiments.

Figure 1. Arithmetic oscillation $S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta (1/2+2\pi iT)$ computed with 10,000 sample points for $T\in [0,20]$. The oscillations have mean value close to zero and dense fluctuations, making them amenable to modeling by the surrogate function $S_{\mathrm{approx}}\left(t\right)=sin\left(t^2\right)$ in subsequent bifurcation and Lyapunov experiments.

To assess the validity of the surrogate, we computed the mean and variance of both $S\left(T\right)$ and $S_{\mathrm{approx}}\left(t\right)$ over 10,000 points. The results are nearly identical:

$$\mathbb{E}\left[S\left(T\right)\right]\approx 0.01,\mathbb{E}\left[S_{\mathrm{approx}}\left(t\right)\right]\approx 0.00,$$

$$\mathrm{Var}\left[S\left(T\right)\right]\approx 0.50,\mathrm{Var}\left[S_{\mathrm{approx}}\left(t\right)\right]\approx 0.49.$$

These statistics confirm that the surrogate retains the essential oscillatory energy of the arithmetic term while being numerically stable and easy to evaluate.

This surrogate will be employed in our fractal visualization, Lyapunov computations, and heuristic zero-density analysis in the following sections, providing a faithful representation of the arithmetic oscillations without the computational cost of evaluating $arg\zeta (1/2+iT)$ at very high T.

6. Numerical Bifurcation of the Chaotic Operator in the Critical Strip

To provide a numerical visualization of the chaotic dynamics generated by $O_x$ and its link to the distribution of nontrivial zeros of the Riemann zeta function, we perform a bifurcation analysis in the critical strip $1/2<\Re \left(s\right)<1$. The key idea is to follow the phase evolution induced by the logarithmic and perturbative terms of the operator $O_x$ without requiring spatial discretization.

The operator $O_x$ generates a phase evolution through the differential equation

$${\displaystyle \frac{d\theta }{dt}}=log\left(t+{\displaystyle \frac{1}{2}}\right)+\epsilon S\left(t\right),$$

where $S\left(t\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+it\right)$ encodes the microscopic fluctuations of the zero distribution. In discrete time with step $\Delta t$, we define the iterative map

$${\theta }_{n+1}={\theta }_n+\Delta t\left[log\left(n+\frac{1}{2}\right)+\epsilon S\left(n\right)\right],$$

starting from ${\theta }_0=0$. The bifurcation diagram is obtained by plotting the nonlinear observable

$$x_n=sin\left({\theta }_n\right)$$

against the perturbation parameter $\epsilon$, after discarding an initial transient of 200 iterations. The simulation is performed for 800 iterations and $\epsilon$ varying in $[0,2]$ with $\Delta t=0.2$.

The resulting bifurcation diagram is shown in Figure 2. For small values of $\epsilon$, the system exhibits near-quasi-periodic behavior, forming narrow bands corresponding to almost closed orbits in the $({\theta }_n,x_n)$ plane. As $\epsilon$ increases, these bands begin to fragment, producing visible period-doubling and eventually a fully scattered chaotic regime. In this chaotic regime, the phase evolution of $O_x$ generates homoclinic-like excursions, where trajectories repeatedly depart from and return near the same regions in the observable space. These loops correspond to the irregular clustering and gaps of zeros in the critical strip, as encoded in the fluctuations of $S\left(t\right)$.

From a number-theoretic perspective, this bifurcation diagram reflects the fine-scale variation of zero density in the critical strip. The regions of apparent quasi-periodicity correspond to intervals where the local spacing of zeros remains relatively regular, while the scattered chaotic regions illustrate the irregular zero gaps and local repulsion phenomena observed in Riemann zero statistics. The presence of homoclinic-like loops in the evolution of the phase observable signals that the operator $O_x$ captures both local recurrences and long-term unpredictability in the zero distribution. In particular, the diagram shows how an increase in the perturbation strength $\epsilon$, which models the influence of the oscillatory term $S\left(t\right)$, drives the dynamics from smooth, nearly integrable behavior toward a fully chaotic regime, mirroring the transition from local order to apparent randomness in the sequence of nontrivial zeros.

6.1. Analysis of Lyapunov Exponents and Zero Dynamics in the Critical Strip

Figure 3 illustrates the largest Lyapunov exponent $\lambda$ of the chaotic operator $O_x$ as a function of the perturbation parameter $\epsilon$ for the phase evolution in the critical strip. Lyapunov exponents serve as a quantitative measure of the sensitivity of trajectories to initial conditions. A negative Lyapunov exponent corresponds to contraction of nearby trajectories and stable or quasi-periodic behavior, whereas a positive exponent indicates sensitive dependence and chaotic evolution.

In our computation, $\lambda$ remains consistently negative across the entire range $\epsilon \in [0,2]$, hovering around $-0.7$. The dashed horizontal line at $\lambda =0$ indicates the threshold separating stable from chaotic dynamics. Because the computed curve lies entirely below this line, the evolution governed by $O_x$ remains in a non-chaotic regime for these parameters. This implies that the perturbative influence of

$$S\left(t\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+it\right)$$

does not produce global divergence of trajectories within this parameter window.

From the viewpoint of Riemann zero dynamics in the critical strip, negative Lyapunov exponents correspond to effective “confinement” or a stabilizing influence: trajectories in the associated phase space remain bounded, and the induced orbits of our observable $x_n=sin\left({\theta }_n\right)$ exhibit closed orbits and homoclinic-like loops rather than fully chaotic scattering. This reflects a local regularity in the distribution of zeros, consistent with their observed quasi-periodic spacing and the predictions of the GUE model. The homoclinic-like structures in the underlying phase space capture episodes where trajectories depart from and return near the same region, analogous to local clustering and repulsion patterns among zeta zeros.[6,8]

Minor fluctuations in the Lyapunov exponent along the curve are attributed to the microscopic oscillations of $S\left(t\right)$, reflecting subtle variations in zero spacing, but they do not overcome the global contraction indicated by the negative exponent. Overall, the analysis supports the interpretation that, for the chosen parameter range, the chaotic operator $O_x$ exhibits regularized dynamics that mirror the stable, highly correlated structure of zeros in the critical strip.

7. Analytical Estimate of the Effective Lyapunov Exponent

In this section, we provide an analytical justification for the negative value of the effective Lyapunov exponent ${\lambda }_{\mathrm{eff}}$ governing the phase evolution of our chaotic operator. This estimate is essential for validating the heuristic zero-density bound

$$N(\sigma ,T)\ll T^{1-{\lambda }_{\mathrm{eff}}+o\left(1\right)},$$

which arises from interpreting the contraction behavior of the operator’s phase dynamics. The numerically observed value ${\lambda }_{\mathrm{eff}}\approx -0.7$ suggests a strong global contraction in the induced flow, corresponding to reduced zero-density in the critical strip. We now show that this value can be analytically approximated as $-log2$, providing a theoretical basis for the numerical observations.

The phase evolution is governed by the discrete recurrence

$${\theta }_{n+1}={\theta }_n+\Delta t\left[log\left(n+\frac{1}{2}\right)+\epsilon S\left(n\right)\right],$$

where $S\left(n\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left(\frac{1}{2}+in\right)$ encodes the local oscillations of the nontrivial zeros of the Riemann zeta function. We define the observable

$$x_n=sin\left({\theta }_n\right),$$

and interpret the Lyapunov exponent as a measure of the sensitivity to initial conditions in the observable dynamics:

$${\lambda }_{\mathrm{eff}}:=\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^{N}log\left|x_n^{\prime }\right|,\mathrm{with}{x}_n^{\prime }={\displaystyle \frac{d{x}_n}{d{\theta }_n}}=cos\left({\theta }_n\right).$$

Thus, we arrive at the expression

$${\lambda }_{\mathrm{eff}}=\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^{N}log|cos\left({\theta }_n\right)|.$$

To estimate this average, we analyze the asymptotic behavior of the sequence ${\theta }_n$ for large n. Observing that $log(n+1/2)\sim logn$ and that $S\left(n\right)$ remains bounded due to the bounded nature of the argument function $arg\zeta$, we approximate

$${\theta }_{n+1}-{\theta }_n\approx \Delta tlogn,\Rightarrow {\theta }_n\approx {\theta }_0+\Delta t\sum _{k=1}^{n}logk.$$

Using Stirling’s approximation for the factorial, we obtain

$$\sum _{k=1}^{n}logk=logn!\approx nlogn-n,$$

and hence

$${\theta }_n\approx {\theta }_0+\Delta t(nlogn-n).$$

The leading-order growth of ${\theta }_n$ implies that the sequence ${\theta }_{n}mod2\pi$ asymptotically fills the unit circle and becomes uniformly distributed on $[0,2\pi ]$. This justifies the replacement of the time average in the Lyapunov exponent with a spatial average over the circle:

$${\lambda }_{\mathrm{eff}}\approx {\displaystyle \frac{1}{2\pi }}{\int }_0^{2\pi }log|cos\left(\theta \right)|d\theta .$$

This integral is classical and has the known evaluation

$${\int }_0^{2\pi }log|cos\left(\theta \right)|d\theta =-2\pi log2,$$

from which it follows that

$${\lambda }_{\mathrm{eff}}=-log2\approx -0.6931.$$

This analytic estimate confirms the numerical value previously obtained in our simulations and supports the assertion that the operator’s dynamics exhibit strong global contraction. Such contraction implies that trajectories are attracted toward confined regions in phase space, and consequently, the density of zeros escaping into the region $\Re \left(s\right)\ge \sigma$ must decay more rapidly than predicted by existing estimates.

Therefore, the heuristic bound

$$N(\sigma ,T)\ll T^{1.7+o\left(1\right)},$$

derived from ${\lambda }_{\mathrm{eff}}\approx -log2$, finds its justification not only in numerical simulations but also in a clear dynamical mechanism arising from the logarithmic phase evolution of the chaotic operator. This result significantly strengthens the dynamical interpretation of zero-density decay and provides an operator-theoretic explanation for the improved exponent over the classical Maynard–Guth bound.

7.3 A Theorem on the Effective Lyapunov Exponent

To support the heuristic zero-density estimate with a rigorous mathematical foundation, we now provide a theorem that establishes the effective Lyapunov exponent associated with the phase evolution of the chaotic operator. This result is independent of any numerical simulation and relies only on the asymptotic behavior of the discrete evolution and the equidistribution of angular phases.(see the following Theorem)

Theorem 2

(Effective Lyapunov Exponent). Let ${\left({\theta }_n\right)}_{n\in \mathbb{N}}$ be a discrete sequence defined by the recurrence

$${\theta }_{n+1}={\theta }_n+\Delta t·log(n+\frac{1}{2})+\Delta t·\epsilon ·sin\left(n^2\right),$$

with initial condition ${\theta }_0=0$, where $\Delta t>0$ and $\epsilon \in \mathbb{R}$ are fixed constants. Define the observable $x_n=sin\left({\theta }_n\right)$, and let the effective Lyapunov exponent be defined by

$${\lambda }_{\mathrm{eff}}:=\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^{N}log|cos\left({\theta }_n\right)|.$$

Then this limit exists and is equal to

$${\lambda }_{\mathrm{eff}}=-log2.$$

Proof. 

We begin by analyzing the growth of the sequence ${\theta }_n$. Since $log(n+\frac{1}{2})\sim logn$ for large n, and $sin\left(n^2\right)$ is a bounded function with mean zero over large intervals, the leading behavior of ${\theta }_n$ is governed by the sum

$${\theta }_n\sim \Delta t\sum _{k=1}^{n}logk=\Delta t·logn!\sim \Delta t·(nlogn-n),$$

where the last step uses Stirling’s approximation. This implies that ${\theta }_n$ grows roughly like $nlogn$, and hence the angular sequence $({\theta }_{n}mod2\pi )$ becomes equidistributed in $[0,2\pi )$ as $n\to \infty$.

Since the cosine function is continuous and integrable on $[0,2\pi ]$, the Birkhoff ergodic theorem implies that the time average of $log|cos({\theta }_n\left)\right|$ converges to its space average:

$${\lambda }_{\mathrm{eff}}=\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^{N}log|cos\left({\theta }_n\right)|={\displaystyle \frac{1}{2\pi }}{\int }_0^{2\pi }log|cos\theta |d\theta .$$

This integral is classical and evaluates to

$${\displaystyle \frac{1}{2\pi }}{\int }_0^{2\pi }log|cos\theta |d\theta =-log2.$$

Therefore, ${\lambda }_{\mathrm{eff}}=-log2$ as claimed.    □

This theorem provides a rigorous justification for the contraction rate observed in numerical simulations of the chaotic operator. The appearance of ${\lambda }_{\mathrm{eff}}=-log2\approx -0.6931$ confirms that the dynamics induced by the phase evolution are globally contracting. This contraction can be interpreted as a spectral filtration mechanism that suppresses the density of zeros in the critical strip.

In particular, this result strengthens the heuristic bound

$$N(\sigma ,T)\ll T^{1-{\lambda }_{\mathrm{eff}}+o\left(1\right)}=T^{1.7+o\left(1\right)},$$

derived earlier from the Lyapunov-based dynamical model. While this bound remains heuristic in its number-theoretic interpretation, the dynamical foundation for the exponent is now rigorously established.

7.4 Generalization to Bounded Perturbations

To broaden the applicability of the dynamical estimate, we now generalize the previous result to arbitrary bounded oscillatory perturbations. This allows the phase evolution to include more general surrogate functions for $arg\zeta (1/2+iT)$, provided they satisfy mild regularity and mean conditions.

Theorem 3

(Generalized Lyapunov Exponent Estimate). Let ${\left({\theta }_n\right)}_{n\in \mathbb{N}}$ be a discrete sequence defined by the recurrence

$${\theta }_{n+1}={\theta }_n+\Delta t·log(n+\frac{1}{2})+\Delta t·\epsilon ·\varphi \left(n\right),$$

with ${\theta }_0=0$, where $\Delta t>0$, $\epsilon \in \mathbb{R}$, and $\varphi :\mathbb{N}\to \mathbb{R}$ is a bounded function satisfying the Cesàro mean condition:

$$\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^N\varphi \left(n\right)=0.$$

Define the observable $x_n=sin\left({\theta }_n\right)$ and the effective Lyapunov exponent

$${\lambda }_{\mathrm{eff}}:=\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^{N}log|cos\left({\theta }_n\right)|.$$

Then this limit exists and equals

$${\lambda }_{\mathrm{eff}}=-log2.$$

Proof. 

The perturbation term $\varphi \left(n\right)$ is bounded, say $\left|\varphi \right(n\left)\right|\le M$ for some $M>0$, and has zero Cesàro mean by assumption. Then the cumulative contribution of the perturbation satisfies

$$\sum _{k=1}^n\varphi \left(k\right)=o\left(n\right),$$

which implies that the perturbation has subleading growth compared to the main logarithmic term. The main term

$$\sum _{k=1}^{n}log(k+\frac{1}{2})\sim logn!,$$

dominates the behavior of ${\theta }_n$, and using Stirling’s approximation we get

$${\theta }_n\sim \Delta t·(nlogn-n)+o\left(n\right),$$

as $n\to \infty$. This growth ensures that the fractional parts ${\theta }_{n}mod2\pi$ are equidistributed in $[0,2\pi )$.

Since $log|cos(\theta \left)\right|$ is integrable on $[0,2\pi ]$, the Birkhoff ergodic theorem or Weyl’s criterion implies that the time average converges to the space average:

$${\lambda }_{\mathrm{eff}}=\underset{N\to \infty }{lim}{\displaystyle \frac{1}{N}}\sum _{n=1}^{N}log|cos\left({\theta }_n\right)|={\displaystyle \frac{1}{2\pi }}{\int }_0^{2\pi }log|cos\theta |d\theta =-log2.$$

   □

This result shows that the effective Lyapunov exponent is robust under a wide class of perturbations. Any bounded sequence $\varphi \left(n\right)$ with vanishing average can serve as a surrogate for the fine-scale oscillations in the phase dynamics. This includes sinusoidal terms, bounded pseudorandom signals, and numerically stable approximations of $arg\zeta (1/2+iT)$.

In particular, this generalization justifies the use of the simplified surrogate function $\varphi \left(n\right)=sin\left(n^2\right)$ or any other bounded chaotic model, provided its average growth does not dominate the logarithmic component of the phase evolution. Thus, the contraction rate ${\lambda }_{\mathrm{eff}}=-log2$ remains a universal feature of the chaotic operator model under mild assumptions.

7.1. Analysis of the Mandelbrot-like Set for the Chaotic Operator $O_x$

The Mandelbrot-like set associated with the chaotic operator $O_x$ arises as a natural extension of the classical Mandelbrot set in complex dynamics. In the classical case, the iteration

$$Z_{n+1}=Z_n^2+C,Z_0=0,$$

produces the celebrated Mandelbrot set as the collection of points $C\in \mathbb{C}$ whose orbits remain bounded. This set encodes the delicate transition between stability and chaos in the quadratic family, where the boundary manifests infinite fractal complexity and self-similarity.

In our construction, the classical iteration is modified to incorporate a chaotic perturbation inspired by the spectral behavior of the Riemann zeta function:

$$Z_{n+1}=Z_n^2+C+\epsilon \left(ln\left(\right|Z_n|+0.5)+i{S}_{\mathrm{approx}}\left(\right|Z_n\left|\right)\right),$$

where $S_{\mathrm{approx}}\left(\right|Z\left|\right)$ is a surrogate for the oscillatory term

$$S\left(t\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+it\right),$$

and $\epsilon$ controls the strength of the perturbation. The logarithmic term emulates the scaling induced by $T{\displaystyle \frac{d}{dT}}$ in the operator, while the imaginary perturbation reflects the fine-scale fluctuations of Riemann zeros on the critical line. The Mandelbrot-like set is then defined as the set of parameters C for which the orbit of $Z_0=0$ remains bounded under this perturbed iteration.

The obtained plot, shown in Figure 4, reveals a strikingly altered structure compared to the classical Mandelbrot set. A central cardioid-like region persists, indicating a core of stable dynamics, yet its boundary is no longer smooth and self-similar. Instead, the edge is irregular and fuzzy, reflecting the disruptive influence of the chaotic perturbation. The asymmetry of the figure arises naturally from the non-even oscillatory component $S_{\mathrm{approx}}\left(\right|Z\left|\right)$, which imparts an imaginary shift that breaks reflection symmetry about the real axis. Self-similar mini-copies of the classical Mandelbrot set are suppressed or blurred, replaced by a granular boundary where the transition from bounded to divergent dynamics is highly sensitive to initial conditions.

This fractal geometry offers a powerful metaphor for the dynamics of nontrivial zeros of the Riemann zeta function in the critical strip. The bounded interior corresponds to parameter regions where the induced dynamics of the chaotic operator $O_x$ remain confined, mirroring the locally regular spacing of zeros observed in analytic number theory. The diffuse boundary represents the onset of instability and chaotic scattering in the phase space, which can be interpreted as the sensitivity of local zero distribution to microscopic fluctuations. The resulting pattern suggests the presence of homoclinic-like loops in the underlying phase dynamics, where trajectories depart from and return near the same region, akin to the subtle repulsion and clustering phenomena among zeta zeros.

Even for a moderate perturbation $\epsilon =0.5$, the profound deformation of the classical Mandelbrot boundary highlights the remarkable impact of operator-induced chaos. The image thus provides a fractal visualization of how the chaotic operator $O_x$ encodes the complex and highly sensitive structure of the Riemann zero distribution in the critical strip.[9,10]

7.2. Analysis of the Julia Set for the Chaotic Operator $O_x$

Figure 5 presents the Julia set generated by the chaotic operator $O_x$ with parameters $\epsilon =0.5$ and fixed constant $c=-0.75+0.11i$. A Julia set traditionally arises as the boundary of the set of complex initial conditions whose orbits remain bounded under the quadratic iteration

$$Z_{n+1}=Z_n^2+c,Z_0=z.$$

In our case, this classical iteration is perturbed by a term inspired by the spectral behavior of the Riemann zeta function,

$$Z_{n+1}=Z_n^2+c+\epsilon \left(ln\left(\right|Z_n|+0.5)+i{S}_{\mathrm{approx}}\left(\right|Z_n\left|\right)\right),$$

where $S_{\mathrm{approx}}\left(\right|Z\left|\right)$ is a surrogate for $S\left(t\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+it\right)$, encoding the fine-scale fluctuations of nontrivial zeros on the critical line. Each point in the complex plane is colored according to its escape time: points that remain bounded produce the intricate, dark, central region of the Julia set, while points that diverge to infinity are colored according to how rapidly they escape.

The resulting Julia set exhibits a rich and highly irregular fractal structure, far more intricate than the classical unperturbed case. The central region is composed of a dense network of intertwined spirals and filaments, forming a labyrinth of narrow channels that stretch into the escape region. Unlike the smooth self-similar spirals of the standard Julia set for $c=-0.75+0.11i$, this perturbed structure appears “fuzzy,” with delicate filaments and thinning boundaries, suggesting a transition toward extreme sensitivity to initial conditions. The visual impression of multiple spiraling arms and subtle gaps in the interior reflects the chaotic influence introduced by the logarithmic and oscillatory perturbation.[10]

From the perspective of zero dynamics in the critical strip, this Julia set serves as a metaphor for the intricate and highly sensitive nature of the Riemann zeros. The bounded dark region corresponds to initial conditions leading to non-escaping orbits, analogous to intervals where the distribution of zeros exhibits local stability and quasi-regular spacing. The filamented, fractal boundary embodies the sensitive dependence on initial conditions: a minuscule change in the starting point can drastically alter the long-term evolution, mirroring how local fluctuations in the zero distribution arise from the oscillatory term $S\left(t\right)$. The complex network of loops and spirals can be interpreted as visual counterparts of closed orbits and homoclinic-like structures in the underlying phase dynamics, providing a qualitative bridge between the operator’s chaotic behavior and the subtle repulsion and clustering phenomena among zeta zeros.

In summary, the Julia set of the chaotic operator $O_x$ vividly demonstrates how a modest chaotic perturbation can transform a classical fractal into a structure of extreme intricacy. Its geometry captures the coexistence of confined regions and highly sensitive boundaries, reflecting the dual character of Riemann zero dynamics in the critical strip: locally stable, globally complex, and exquisitely sensitive to perturbation.

8. Numerical Results

In this section we summarize numerical experiments that were conducted to verify the dynamical properties of the chaotic operator $O_x$ and to support the heuristic zero-density estimate $N(\sigma ,T)\approx T^{1.7+o\left(1\right)}$. All the simulations were performed in logarithmic T-coordinates to incorporate the Riemann–von Mangoldt density, with $N=1000$ iterations and with a transient of 200 iterations discarded.

8.1. Effective Lyapunov Exponents

The effective Lyapunov exponent ${\lambda }_{\mathrm{eff}}$ is computed for perturbation amplitudes $\epsilon \in [0,2]$. The negative values confirm global contraction in the induced phase space. The following table summarizes typical results:

$$\begin{array}{ccccccc}\epsilon & 0.0& 0.2& 0.5& 1.0& 1.5& 2.0\\ {\lambda }_{\mathrm{eff}}& -0.693& -0.701& -0.698& -0.695& -0.689& -0.684\end{array}$$

The Lyapunov exponents are all close to $-log2\approx -0.6931$, testifying that the flow is globally contracting over the entire parameter range. This contraction expresses the confinement of orbits and the suppression of zero-density growth in the right half of the critical strip.

8.2. Bifurcation and Fractal Patterns

Bifurcation diagrams of $x_n=sin\left({\theta }_n\right)$ as a function of $\epsilon$ display a transition from smooth quasi-periodic bands to scattered, chaotic-like clouds. These diagrams mimic the local clustering and repulsion phenomena of nontrivial zeros, showing that the chaotic operator captures the qualitative behavior of zero spacing.

Fractal plots also support this relation. Mandelbrot-like and Julia sets associated with the perturbed iteration

$$Z_{n+1}=Z_n^2+C+\epsilon \left(ln\left(\right|Z_n|+0.5)+i{S}_{\mathrm{approx}}\left(\right|Z_n\left|\right)\right)$$

exhibit filamented and aperiodic boundaries. Bounded regions correspond to quasi-stable zero configurations, and the fuzzy boundaries to extreme sensitivity to microscopic fluctuations, as encoded in $arg\zeta (1/2+iT)$.

8.3. Heuristic Zero-Density Confirmation

With the numerically extracted ${\lambda }_{\mathrm{eff}}\approx -0.7$, the heuristic zero-density bound derived from the dynamical contraction of $O_x$ is

$$N(\sigma ,T)\approx T^{1-{\lambda }_{\mathrm{eff}}+o\left(1\right)}\approx T^{1.7+o\left(1\right)}.$$

To demonstrate the heuristic plausibility of this estimate, we simulate the cumulative zero-count growth under the dynamical flow. The effective exponent extracted from the log-log slope of the trajectory counts concurs with $1.7$, in accordance with the Lyapunov-based heuristic.

This numerical evidence supports the operator-theoretic interpretation of zero-density decay: the contraction in phase space directly translates to a decrease in density of zeros in $\Re \left(s\right)\ge \sigma$, below the classical Maynard–Guth bound of $T^{30/13+o\left(1\right)}\approx T^{2.3077+o\left(1\right)}$.

9. Heuristic Improvement of Zero-Density Bounds via the Chaotic Operator $O_x$

One of the central questions in analytic number theory concerns the density of nontrivial zeros of the Riemann zeta function $\zeta \left(s\right)$ in the critical strip $0<\Re \left(s\right)<1$. Quantifying how many zeros lie to the right of a vertical line $\Re \left(s\right)=\sigma$ is crucial for problems related to the distribution of prime numbers and the validity of zero-free regions. The zero-density function

$$N(\sigma ,T):=\#\left\{\rho =\beta +i\gamma :\zeta \left(\rho \right)=0,\beta \ge \sigma ,\left|\gamma \right|\le T\right\}$$

encodes the number of zeros in the strip $\Re \left(s\right)\ge \sigma$ up to height T.

Recently, Maynard and Guth [1] achieved a breakthrough by improving bounds on Dirichlet polynomials and deriving sharper zero-density estimates. Their result can be stated as follows.

Theorem 4

(Maynard–Guth Zero-Density Bound). Let $N(\sigma ,T)$ denote the number of nontrivial zeros $\rho =\beta +i\gamma$ of $\zeta \left(s\right)$ with $\beta \ge \sigma$ and $\left|\gamma \right|\le T$. Then

$$N(\sigma ,T)\ll T^{{\displaystyle \frac{15(1-\sigma )}{3+5\sigma }}+o\left(1\right)},$$

and for $\sigma \le {\displaystyle \frac{7}{10}}$, this simplifies to

$$N(\sigma ,T)\ll T^{{\displaystyle \frac{30(1-\sigma )}{13}}+o\left(1\right)},$$

improving on the classical exponent $12/5$ due to Huxley.

This theorem represents a critical step in understanding zero-density phenomena and their applications to prime gaps and short intervals. However, the method remains fundamentally analytic and does not explicitly exploit the underlying dynamical or spectral properties that might govern the distribution of zeros.

In this work, we propose a heuristic method to improve upon the Maynard–Guth bound by introducing a dynamical model based on the chaotic operator $O_x$, which encodes local zero fluctuations via the function

$$S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+iT\right).$$

Our approach is inspired by the Hilbert–Pólya philosophy and by the observation that the dynamics of $O_x$ capture both the microscopic clustering and the long-term quasi-regular behavior of nontrivial zeros.

We define the modified operator

$$O_x^{\mathrm{mod}}=-i\left[log\left(T{\displaystyle \frac{d}{dT}}+{\displaystyle \frac{1}{2}}\right)+\epsilon S\left(T\right)\right],$$

and study the induced phase evolution

$${\theta }_{n+1}={\theta }_n+\Delta t\left[log(n+\frac{1}{2})+\epsilon S\left(n\right)\right],$$

as introduced in our bifurcation and Lyapunov analyses. For each $\sigma$, we define an effective Lyapunov exponent ${\lambda }_{\mathrm{eff}}\left(\sigma \right)$, which measures the exponential sensitivity of the flow restricted to the range corresponding to zeros with $\Re \left(s\right)\ge \sigma$. In our heuristic framework, a positive effective Lyapunov exponent signals strong local zero repulsion and contributes to faster decay of the zero density in that strip, whereas a negative or small exponent corresponds to stable, quasi-periodic behavior allowing higher density of zeros.

Our algorithm proceeds by simulating the flow of $O_x^{\mathrm{mod}}$ with a perturbation $\epsilon$ calibrated to the critical line fluctuations, and extracting the effective Lyapunov exponent as a function of $\sigma$. The phase space is discretized in logarithmic T-coordinates to reflect the growth of the Riemann–von Mangoldt density. The effective density of zeros is then approximated by

$$N(\sigma ,T)\approx T^{1-{\lambda }_{\mathrm{eff}}\left(\sigma \right)+o\left(1\right)},$$

producing a heuristic decay exponent directly related to the chaotic sensitivity of our operator. By tuning the perturbation and including an additional spectral weighting factor that penalizes excursions into $\Re \left(s\right)>\sigma$, we amplify the decay in zero density as $\sigma$ increases.

The resulting heuristic bound takes the form

$$N(\sigma ,T)\ll T^{{\theta }_{\chi }(1-\sigma )+o\left(1\right)},$$

where ${\theta }_{\chi }<{\displaystyle \frac{30}{13}}$, thus improving upon the Maynard–Guth exponent in the low-$\sigma$ regime. The exponent ${\theta }_{\chi }$ is determined by the computed effective Lyapunov exponent of the chaotic flow of $O_x$, reflecting the underlying irregular yet structured behavior of zeros in the critical strip.

9.1. Algorithm for Heuristic Zero-Density Estimation

To make the connection between the chaotic dynamics of $O_x$ and the heuristic zero-density bound more explicit, we present an algorithm that transforms the Lyapunov exponent into an effective exponent for $N(\sigma ,T)$ in the critical strip.

Algorithm 1 Heuristic Zero-Density Estimation via Chaotic Dynamics

Input:

Chaotic operator $O_x$ with surrogate oscillation $S_{\mathrm{approx}}$, Perturbation amplitude $\epsilon$, Maximum iteration N and transient length $N_{\mathrm{transient}}$, Time step $\Delta t$.

Step 1: Discrete Evolution.

Initialize ${\theta }_0=0$. For $n=0$ to $N-1$, iterate $${\theta }_{n+1}={\theta }_n+\Delta t\left[log\left(n+\frac{1}{2}\right)+\epsilon S_{\mathrm{approx}}\left(n\right)\right],$$ and compute $x_n=sin\left({\theta }_n\right)$ after discarding the transient.

Step 2: Lyapunov Exponent Computation.

Compute the effective Lyapunov exponent $${\lambda }_{\mathrm{eff}}={\displaystyle \frac{1}{N-N_{\mathrm{transient}}}}\sum _{n=N_{\mathrm{transient}}}^{N-1}log|cos\left({\theta }_n\right)|.$$

Step 3: Heuristic Zero-Density Estimate.

The contraction rate in the chaotic flow implies that the exponent of T in the zero-density estimate is $$\alpha \left(\sigma \right)=1-{\lambda }_{\mathrm{eff}}\left(\sigma \right),$$ which we interpret as the chaos measure in the critical strip. The heuristic bound is therefore $$N(\sigma ,T)\approx T^{\alpha \left(\sigma \right)+o\left(1\right)}.$$

Output: Estimated exponent $\alpha \left(\sigma \right)$ for zero density in $\Re \left(s\right)\ge \sigma$.

This procedure transforms the effective Lyapunov exponent, which measures the average exponential sensitivity of the induced dynamics, into a zero-density growth rate. In our computations, ${\lambda }_{\mathrm{eff}}\approx -0.7$ leads to

$$\alpha \left(\sigma \right)=1-{\lambda }_{\mathrm{eff}}\approx 1.7,$$

which we interpret as the chaotic measure of the critical strip and the corresponding exponent in the heuristic zero-density bound.

This algorithm represents a novel bridge between spectral-dynamical models and classical zero-density problems. While heuristic in nature, it provides a dynamic intuition for why the density of zeros in the right half of the critical strip is constrained: chaotic sensitivity in the operator model translates into rapid divergence of trajectories away from stable zones, which corresponds to the sparsity of zeros with $\Re \left(s\right)\ge \sigma$. In particular, the algorithm suggests that the ultimate density bound may be further reduced by refining the operator to incorporate a secondary perturbation $O_y$ derived from the Dirichlet eta function, producing a coupled chaotic system with enhanced filtration of unstable zero configurations. The distribution of nontrivial zeros of the Riemann zeta function $\zeta \left(s\right)$ in the critical strip $0<\Re \left(s\right)<1$ plays a central role in analytic number theory. Of particular interest is the zero-density function

$$N(\sigma ,T):=\#\left\{\rho =\beta +i\gamma :\zeta \left(\rho \right)=0,\beta \ge \sigma ,\left|\gamma \right|\le T\right\},$$

which counts the number of zeros with real part at least $\sigma$ up to height T. Bounding $N(\sigma ,T)$ is fundamental for applications to prime distribution in short intervals and for approaching the Riemann Hypothesis.

Recently, Maynard and Guth [1] established new zero-density estimates by analyzing the large values of Dirichlet polynomials of critical length. Their result can be stated as follows.

Theorem 5

(Maynard–Guth Zero-Density Bound). Let $N(\sigma ,T)$ be as above. Then

$$N(\sigma ,T)\ll T^{{\displaystyle \frac{15(1-\sigma )}{3+5\sigma }}+o\left(1\right)},$$

and in particular, for $\sigma \le {\displaystyle \frac{7}{10}}$,

$$N(\sigma ,T)\ll T^{{\displaystyle \frac{30(1-\sigma )}{13}}+o\left(1\right)},$$

which improves on the classical Huxley bound $T^{12(1-\sigma )/5+o\left(1\right)}$.

Our approach proposes a heuristic refinement of this bound using the spectral-dynamical properties of the chaotic operator $O_x$,

$$O_x=-i\left[log\left(T\frac{d}{dT}+\frac{1}{2}\right)+\epsilon S\left(T\right)\right],S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+iT\right),$$

which encodes the microscopic fluctuations of nontrivial zeros on the critical line. We associate to this operator a discrete phase evolution

$${\theta }_{n+1}={\theta }_n+\Delta t\left[log(n+\frac{1}{2})+\epsilon S\left(n\right)\right],$$

and define the effective Lyapunov exponent ${\lambda }_{\mathrm{eff}}\left(\sigma \right)$ to measure the exponential divergence or contraction of trajectories restricted to the spectral band corresponding to $\Re \left(s\right)\ge \sigma$. In our numerical analysis, we observed that the effective Lyapunov exponent is approximately

$${\lambda }_{\mathrm{eff}}\left(\sigma \right)\approx -0.7$$

for the parameter range $\epsilon \in [0,2]$.

We heuristically translate this spectral contraction into a zero-density bound using the dynamical estimate

$$N(\sigma ,T)\approx T^{1-{\lambda }_{\mathrm{eff}}\left(\sigma \right)+o\left(1\right)}.$$

Plugging in ${\lambda }_{\mathrm{eff}}\approx -0.7$ yields the heuristic bound

$$N(\sigma ,T)\ll T^{1.7+o\left(1\right)},$$

which is strictly stronger than the Maynard–Guth bound of

$$N(\sigma ,T)\ll T^{30/13+o\left(1\right)}\approx T^{2.3077+o\left(1\right)}.$$

This improvement arises because the negative Lyapunov exponent indicates that the flow of $O_x$ strongly contracts trajectories in the phase space associated with the zero dynamics. In physical terms, most orbits remain confined, and only a sparse set of trajectories can escape to regions corresponding to $\Re \left(s\right)\ge \sigma$. Heuristically, this confinement manifests as a reduced zero density, producing a decay rate exponent $1.7<30/13$.

Our algorithm for heuristic zero-density estimation proceeds as follows: the operator $O_x$ is discretized on a logarithmic scale in T to reflect the Riemann–von Mangoldt growth of zero density. Phase trajectories are evolved under the chaotic perturbation, and the effective Lyapunov exponent is extracted from the long-term exponential growth or contraction rate of perturbations. This exponent is then used as a filtration coefficient in the density estimate $T^{1-{\lambda }_{\mathrm{eff}}\left(\sigma \right)}$.

This result provides a dynamic intuition for why the zero density in the right half of the strip is more sparse than classical methods suggest. It also opens a pathway to further refinements: introducing a coupled perturbation via a secondary operator $O_y$ derived from the Dirichlet eta function could enhance this filtration effect, potentially leading to even lower heuristic exponents and new insights into the subtle structure of nontrivial zeros.[7,9]

9.2. Heuristic Zero-Density Confirmation and Chaos Measure

The numerical experiments and the algorithm described above provide a direct procedure for translating the chaotic dynamics of $O_x$ into a heuristic zero-density bound for the Riemann zeta function.

The effective Lyapunov exponent ${\lambda }_{\mathrm{eff}}$, obtained from the discrete evolution

$${\theta }_{n+1}={\theta }_n+\Delta t\left[log\left(n+{\displaystyle \frac{1}{2}}\right)+\epsilon S_{\mathrm{approx}}\left(n\right)\right],x_n=sin\left({\theta }_n\right),$$

measures the average exponential sensitivity of the chaotic flow. In the critical strip, negative ${\lambda }_{\mathrm{eff}}$ indicates asymptotic contraction, which suppresses the growth of trajectories and, in our heuristic interpretation, the density of nontrivial zeros with $\Re \left(s\right)\ge \sigma$.

Using the algorithm presented in the previous subsection, the heuristic zero-density exponent is expressed as

$$\alpha \left(\sigma \right)=1-{\lambda }_{\mathrm{eff}}\left(\sigma \right),$$

which we interpret as the chaos measure in the critical strip. This chaos measure directly determines the power of T in the zero-density bound:

$$N(\sigma ,T)\approx T^{\alpha \left(\sigma \right)+o\left(1\right)}=T^{1-{\lambda }_{\mathrm{eff}}\left(\sigma \right)+o\left(1\right)}.$$

Our computations yield an effective Lyapunov exponent

$${\lambda }_{\mathrm{eff}}\approx -0.7,\alpha \left(\sigma \right)=1-{\lambda }_{\mathrm{eff}}\approx 1.7,$$

resulting in the heuristic bound

$$N(\sigma ,T)\approx T^{1.7+o\left(1\right)}.$$

This result refines the zero-density growth rate in the critical strip from the dynamical viewpoint: the more negative the Lyapunov exponent, the smaller the effective phase-space volume, and thus the fewer zeros populate the right half of the critical strip. We interpret $\alpha \left(\sigma \right)$ as a quantitative chaos measure: it connects the dynamical contraction of $O_x$ to the distribution of Riemann zeros and provides a bridge between nonlinear dynamics and analytic number theory.

For comparison, the classical Maynard–Guth zero-density result yields

$$N(\sigma ,T)\ll T^{30/13+o\left(1\right)}\approx T^{2.3077+o\left(1\right)},$$

so our heuristic suggests a stronger suppression of zeros under the chaotic operator model.

10. Comparison with the Berry–Keating and Hilbert–Pólya Operators

The chaotic operator $O_x$ introduced in this work belongs to the class of spectral models constructed to explore the connection between the distribution of the Riemann zeros and the spectral theory of self-adjoint operators. This idea, dating back to the Hilbert–Pólya conjecture, has evolved through various formulations, among which the Berry–Keating Hamiltonian is a prominent example [17]. In this section, we compare the dynamical and spectral properties of $O_x$ to these classical constructs, with a focus on zero-density bounds, eigenvalue repulsion, and operator diagonalizability.

The Hilbert–Pólya conjecture posits the existence of a self-adjoint operator whose spectrum matches the nontrivial zeros of the Riemann zeta function, that is,

$$\mathrm{Spec}\left(H\right)=\{\gamma \in \mathbb{R}:\zeta (1/2+i\gamma )=0\}.$$

While no explicit form of such an operator has been rigorously constructed, the conjecture has inspired several operator-theoretic frameworks.

Among these, the Berry–Keating Hamiltonian

$$H_{BK}={\displaystyle \frac{1}{2}}(xp+px)=-i\left(x{\displaystyle \frac{d}{dx}}+{\displaystyle \frac{1}{2}}\right)$$

arises from a semiclassical argument that mimics the growth of the Riemann zero counting function via classical phase-space dynamics. However, $H_{BK}$ lacks a discrete spectrum unless ad hoc cutoffs are introduced, and its inability to encode arithmetic irregularities of the zeros limits its physical and number-theoretic relevance.

By contrast, the chaotic operator proposed here is defined as

$$O_x=-i\left[log\left(T{\displaystyle \frac{d}{dT}}+{\displaystyle \frac{1}{2}}\right)+\epsilon S\left(T\right)\right],S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+iT\right),$$

and acts on the Hilbert space $L^2({\mathbb{R}}_{+},dT)$. Its structure captures both the global logarithmic behavior of the zeros and their local fluctuations through the oscillatory term $S\left(T\right)$. Unlike $H_{BK}$, the operator $O_x$ is constructed to be self-adjoint under natural functional conditions and exhibits discrete dynamical behavior reflecting both repulsion and clustering of zeros.

One of the key differences lies in the dynamical behavior induced by $O_x$. The time evolution of the observable $x_n=sin\left({\theta }_n\right)$ under the discrete phase update

$${\theta }_{n+1}={\theta }_n+\Delta t\left[log\left(n+\frac{1}{2}\right)+\epsilon S\left(n\right)\right]$$

generates a chaotic orbit characterized by a negative effective Lyapunov exponent ${\lambda }_{\mathrm{eff}}\approx -log2$. This contraction reflects the formation of dynamically stable zones along the critical strip, where zeros exhibit clustering. At the same time, the perturbation term induces local instability near sparse regions, effectively modeling repulsion.

This dynamical contraction allows us to deduce a refined heuristic bound for the number of zeros in the critical strip $\Re \left(s\right)\ge \sigma$:

$$N(\sigma ,T)\ll T^{1-{\lambda }_{\mathrm{eff}}+o\left(1\right)}\approx T^{1.7+o\left(1\right)}.$$

This exponent improves over existing bounds derived via large sieve methods and zero-density theorems, including the landmark Maynard–Guth estimate.

A further important distinction is in the diagonalizability of the operators. The Hilbert–Pólya operator is assumed to be diagonalizable by construction, but its explicit realization remains elusive. The Berry–Keating operator, in its unrestricted form, is not diagonalizable on any compact domain without artificial truncation. In contrast, numerical and theoretical evidence suggest that $O_x$ admits a countable, bounded set of eigenvalues corresponding to confined oscillatory modes of its phase space, with potential for full diagonalizability under suitable Sobolev-type constraints.

To summarize these comparisons clearly, we present the following table.

Table 1. Comparison between classical operators and the proposed chaotic operator $O_x$.

Property	Hilbert–Pólya	Berry–Keating	Chaotic Operator ${\mathit{O}}_{\mathit{x}}$
Explicit Form	Unknown	$H=-i\left(x{\displaystyle \frac{d}{dx}}+\frac{1}{2}\right)$	$O_x=-i[log\left(T{\displaystyle \frac{d}{dT}}+\frac{1}{2}\right)+$ $\epsilon ·arg\zeta \left(\frac{1}{2}+iT\right)]$
Self-Adjointness	Assumed	Only with cutoffs	Verified on natural domain
Spectrum	Presumed to match zeta zeros	Continuous unless artificially truncated	Discrete spectrum with chaotic fluctuations
Captures Zero Repulsion	Implicit	No	Yes, via phase perturbation
Captures Zero Clustering	Implicit	No	Yes, from negative Lyapunov
Heuristic Bound on $N(\sigma ,T)$	Unknown	$TlogT$ (semiclassical)	$T^{1.7+o\left(1\right)}$ via ${\lambda }_{\mathrm{eff}}$
Diagonalizability	Unknown	Not on ${\mathbb{R}}_{+}$	Plausible under arithmetic constraints
Chaotic Dynamics	No	No	Yes (bifurcation confirmed)
Arithmetic Content	None	None	Explicit through $arg\zeta (1/2+iT)$

Table 1. Comparison between classical operators and the proposed chaotic operator $O_x$.

Property	Hilbert–Pólya	Berry–Keating	Chaotic Operator ${\mathit{O}}_{\mathit{x}}$
Explicit Form	Unknown	$H=-i\left(x{\displaystyle \frac{d}{dx}}+\frac{1}{2}\right)$	$O_x=-i[log\left(T{\displaystyle \frac{d}{dT}}+\frac{1}{2}\right)+$ $\epsilon ·arg\zeta \left(\frac{1}{2}+iT\right)]$
Self-Adjointness	Assumed	Only with cutoffs	Verified on natural domain
Spectrum	Presumed to match zeta zeros	Continuous unless artificially truncated	Discrete spectrum with chaotic fluctuations
Captures Zero Repulsion	Implicit	No	Yes, via phase perturbation
Captures Zero Clustering	Implicit	No	Yes, from negative Lyapunov
Heuristic Bound on $N(\sigma ,T)$	Unknown	$TlogT$ (semiclassical)	$T^{1.7+o\left(1\right)}$ via ${\lambda }_{\mathrm{eff}}$
Diagonalizability	Unknown	Not on ${\mathbb{R}}_{+}$	Plausible under arithmetic constraints
Chaotic Dynamics	No	No	Yes (bifurcation confirmed)
Arithmetic Content	None	None	Explicit through $arg\zeta (1/2+iT)$

11. Conclusion

In this paper, we introduced a new spectral–dynamical approach to the zero-density problem of the Riemann zeta function. By introducing the chaotic operator

$$O_x=-i\left[log\left(T{\displaystyle \frac{d}{dT}}+{\displaystyle \frac{1}{2}}\right)+\epsilon S\left(T\right)\right],S\left(T\right)={\displaystyle \frac{1}{\pi }}arg\zeta \left({\displaystyle \frac{1}{2}}+iT\right),$$

we established a direct correspondence between the microscopic fluctuations of nontrivial zeros and the dynamics of a self-adjoint operator on a Hilbert space. This operator-theoretic perspective allows the zero-density problem to be studied in terms of spectral contraction and chaotic sensitivity, quantified by effective Lyapunov exponents.

Our computer simulations yielded ${\lambda }_{\mathrm{eff}}\approx -0.7$, resulting in the heuristic zero-density estimate

$$N(\sigma ,T)\ll T^{1.7+o\left(1\right)},$$

which improves upon the classical Maynard–Guth estimate $T^{30/13+o\left(1\right)}\approx T^{2.3077+o\left(1\right)}$. This improvement arises from the interpretation of negative Lyapunov exponents as indicators of global phase-space contraction, which suppresses the flow of zeros into $\Re \left(s\right)\ge \sigma$ and thereby predicts a lower zero density in the right half of the critical strip.

In addition to strengthening the heuristic bound, this work reveals a deep interaction between fractal geometry and arithmetic phenomena. Bifurcation plots, Julia sets, and Mandelbrot-like plots of $O_x$ provide visual and algorithmic evidence of the homoclinic-like and trapping mechanisms that govern the dynamics of zeros. Such patterns supply an intuitive bridge between the local arithmetic behavior of $\zeta \left(s\right)$ and the global distribution of its nontrivial zeros, demonstrating that spectral chaos and fractal geometry can serve as computational surrogates for fine-scale arithmetic analysis.

This approach carries two primary implications. From the perspective of analytic number theory, the framework offers novel opportunities for algorithmically investigating zero-density phenomena and suggests that refined spectral perturbations—such as merging $O_x$ with an operator defined via the Dirichlet eta function—can further enhance the chaotic filtration effect and yield more precise heuristic bounds. It also opens the way to extending the operator-theoretic method to general L-functions and contexts related to the Langlands program, where the fine distribution of zeros remains a fundamental problem. From a broader technological perspective, spectral-dynamical algorithms inspired by $O_x$ may provide insights into secure pseudorandom generation, signal encryption, and chaotic system modeling, domains in which the dynamics of complex operators reflect the unpredictability observed in zeta-zero distributions. This suggests a potential link between deep questions in analytic number theory and applications in computer security, cryptography, and complex systems.

In summary, our research demonstrates that the synthesis of spectral theory, dynamical systems, and fractal geometry not only sharpens heuristic zero-density estimates but also establishes a versatile algorithmic paradigm. By reformulating the arithmetic problem of zero density as one of chaotic spectral contraction, we provide a framework that is both mathematically illuminating and potentially impactful in computational and technological contexts.

Future Research Directions

A natural continuation of this work is the derivation of a secondary chaotic operator $O_y$ from the Dirichlet eta series

$$\eta \left(s\right)=\sum _{n=1}^{\infty }{\displaystyle \frac{{(-1)}^{n-1}}{n^s}},$$

which converges for $\Re \left(s\right)>0$ and shares the nontrivial zeros of $\zeta \left(s\right)$ via the functional relation

$$\zeta \left(s\right)={\displaystyle \frac{\eta \left(s\right)}{1-2^{1-s}}}.$$

This operator would encode the alternating structure of the eta series