Experimental Observations on the Uncomputability of the Riemann Hypothesis
Mathematics Department, University of Auckland
PDF (with full size equations)
Abstract: This paper seeks to explore whether the Riemann hypothesis falls into a class of putatively unprovable mathematical conjectures, which arise as a result of unpredictable irregularity. It also seeks to provide an experimental basis to discover some of the mathematical enigmas surrounding these conjectures, by providing Matlab and C programs which the reader can use to explore and better understand these systems (see appendix 6).
Fig 1: The Riemann functions and : absolute value in red, angle in green. The pole at z = 1 and the non-trivial zeros on x = ½ showing in as a peak and dimples. The trivial zeros are at the angle shifts at even integers on the negative real axis. The corresponding zeros of show in the central foci of angle shift with the absolute value and angle reflecting the function’s symmetry between z and 1 - z. If there is an analytic reason why the zeros are on x = ½ one would expect it to be a manifest property of the reflective symmetry of .
The Riemann hypothesis,  remains the most challenging unsolved problem in mathematics at the beginning of the third millennium. The other two problems of similar status, Fermat’s last theorem  and the Poincare conjecture  have both succumbed to solutions by Andrew Wiles and Grigori Perelman in tours de force using a swathe of advanced techniques from diverse mathematical areas. Fermat’s last theorem states that no three integers a, b, c can satisfy for n > 2. The Poincare conjecture states that any 3-manifold (a space locally like n-dimensional space, such as a sphere, torus or Klein bottle) on which any loop can be continuously shrunk to a point is a 3-sphere. Both of these, despite being difficult problems, have a justifiable case that a solution ought to exist, and that they are not undecidable propositions.
The anticipation in the mathematical community thus remains largely focused on the notion that the Riemann hypothesis is in-principle a provable proposition, which is also eminently plausible, and indeed confirmed by all numerical instances so far discovered. Nevertheless the problem has resisted all attempts to close in on it from areas at least as diverse as the other two, so one might not be entirely foolish to suggest that there may be mathematical barriers to proving the Riemann hypothesis that are fundamentally different.
"Mathematicians have tried in vain to discover some order in the sequence of prime numbers
but we have every reason to believe that there are some mysteries
which the human mind will never penetrate." Euler 1751
The Riemann hypothesis states that the Zeta function   has all its non-trivial zeros on the vertical line x = ½ in the complex plane. This like the previous two problems is an intuitively obvious result, which has great plausibility, since billions of its zeros do lie on the line, but is it, if it is logically equivalent to potentially unprovable, or even undecidable propositions about large primes?
The relationship between the zeta sum formula and the product of primes was discovered by Leonhard Euler and is equivalent to prime sieving.
Zeta also has relationships with the primes at real integer values of s, where the Euler product formula can be used to show that the probability that s integers are relatively prime is 1/ .
While both sides of  are convergent only for real(s) > 1, using the above construction, we can see that we can extend convergence to real(s) > 0 using Dirichlet’s Eta  :
In 1739 Euler gave a bizarre proof which would be frowned on by teachers of modern mathematics, by playing with seemingly impossible values of the real zeta function in mischevious ways: , so differentiating .
Substituting x = -1, he then produced . But this is by definition , so expressing zeta in terms of eta, as above, he got , which caused Abel to declare "The divergent series are the invention of the devil, and it is a shame to base on them any demonstration whatsoever." However over a hundred years later in 1859 Riemann would show that these are precisely the values of eta and zeta realized by analytic continuation.
Riemann  analytically extended the zeta function to the entire complex plane, except the simple pole at z = 1, by considering the integral definition of the gamma function. By extending the resulting path integral to the entire complex plane, he then established the functional equation:
thus enabling to be extended to using reflection about the line x = 1/2 (see appendix 5 for Lanczos approximation to gamma used in generating zeta in the complex plane).
We will now give a complete derivation of the analytic continuation of  to make the process as clear as possible. We start with the following integral expression using the gamma function in the form of a Mellin transform:
Now let , where C is a path from infinity round a small circle not enclosing any of the poles at and back to infinity, as the integral round the circular region tends to zero for Re(s)>1, as the radius tends to zero. So , and this gives a uniformly convergent integral function of s providing an analytic continuation of over the entire plane.
If we now take another path Cn enclosing a larger square with corners as shown right we have:
since the outer integral tends to zero and the residue pairs at Ī2inπ separating the two contour paths above are
We can now derive the other symmetrical forms of the functional equation (check):
This reflectivity relation then became the basis for Riemann's function, which is symmetrical about x = 1/2 and whose zeros are identical to those of .
thus enabling to be extended to using reflection about the line x = 1/2.
The trivial zeros then arise from the zeros of the sine on the negative real axis and the non-trivial from itself in the critical strip.
By careful contour integration Riemann then established that the irregularly-spaced zeros between 0 and t grow according to . More precisely, the number of zeros between 0 and t has been found based on contour integration round the zeros and poles using functional equation and the Cauchy argument principle to be:
The argument of the product is the sum of the arguments, but the last term , measuring fluctuations about the average, which jumps by 1 at each zero, and declines with derivative , grows extremely slowly in the average as . Numerical calculations confirm that S grows very slowly: |S(T)| < 1 for T < 280, |S(T)| < 2 for T < 6800000, and the largest value of |S(T)| found so far  is aroud 3.2. It is thus hard to predict the eventual behavior of the zeros, for extremely large numbers are required before this factor becomes significant and hence the fact that very large computed zeros are on the critical line doesn’t necessarily establish the likelihood that all non-trivial zeros are on the line.
By turning to the logs of and Riemann establishes a formula for the primes:
where F(x) is the number of primes less than x growing by ½ at primes.
Using a Fourier transform, he establishes that the primes can be counted by forming an infinite series of integrals over the zeros of , remarking in the process his hypothesis that the zeros of were real, or equivalently, those of were on x = ½:
The first term corresponds to the pole, the second to the non-trivial zeros, the third to the trivial zeros and the last is . By the Möbius inversion formula we arrive at: , where is the prime counting function, the Möbius function distinct primes and 0 otherwise (appendix 4).
The simplest version of this formula and the one we shall use  is: , where the von Mangoldt function and 0 otherwise (appendix 4), are the zeros of , and the summation is performed over zeros of increasing .
To derive this  , we start with Perron's formula, again based on contour integration:
Now if we move the contour away to the left to Re(s) = -N for a large odd integer N for each zero and pole of we have order 1 so the values are at s = 0, x at s = 1, and for a zero of order m. Hence we have
Now the integrand tends to 0 as N tends to infinity, and the zero contributes , giving a summed weighting of , so we get
Applying this formula numerically is effectively assuming the Riemann hypothesis, because, in performing a calculation, we are working on the basis that all the zeros we are using are on x = ½, and since it is the real part of which determines the magnitude of the contribution of each of these terms to the fluctuations in the prime distribution, the prime distribution that emerges is a practical application of the Riemann hypothesis.
The effect of the process is to create an integral transform of the distribution of primes, with the terms involving the non-trivial zeros forming a series of superimposed wave functions very similar to a Fourier series as can be seen from fig 2. This integral transform converges rapidly for the first few primes but more slowly as the primes grow larger.
Fig 2: The prime counting function iterated (above) for 10, 40 and 100 zero pairs of in the range of 1 to 20, and (below) for 100, 400 and 1000 zero pairs in the range 100 to 120 shows evidence of wave superposition of contributions from each of the successive zeros, and of the increasing number of zeros required to resolve higher primes corresponding to the logarithmically closer spacing of higher zeros.
Looking from the other end of the transform, we can examine how the product formula converges to the zeta zeros as the number of primes in the product increases. There are two ways of doing this. In fig 2b is shown the Fourier sin transform of the prime counting function of fig 2, showing coincidence between the major fluctuations of the transform and the zeta zeros.
Fig 2b (Top left): Fourier transform of the prime counting function, showing coincidence with zeta zeros (blue above) confirms the zeta zeros are in a sense 'holograms' of the primes . (Bottom left) Function showing the fluctuations of the Riemann estimate about compared with the same function substituting Li(x) for (red). RH has been shown to be equivalent to the statement (for all x >= 2657) because the explicit formula shows the magnitude of the oscillations of primes around their expected position is controlled by the real parts of the zeros of the zeta function. Littlewood proved that there exists a value where without determining what these values would be. Skewes then produced an unimaginably huge upper bound, of which Hardy said "the truth has defeated not only all the evidence of the facts and of common sense, but even a mathematical imagination as powerful as that of Gauss"  . The bound came down more recently thanks to te Riele and later Hudson and Bayes to . (Top right) A similar conjecture by Mertens that (see Mobius function) which would have proved the Riemann hypothesis was also found false at a value of around by Odyzko's colleague Herman te Riele. (Bottom right) Lindelof hypothesis that grows more slowly that any fixed power of t remains unproven despite being 'easier' than RH.
Let us now examine how zeros with real part larger or smaller than 1/2 would affect the explicit formula for the prime distribution.
Fig 2c: The explicit formula gives a strong indication of why all the zeros may have to be on the critical line. When they are, (centre), the individual contributions of pairs result in quadratically increasing amplitudes of the oscillating wave function (below, summed lower right black). When these are off the line (above left and right) the prime distribution is no longer constant between primes because the summed pairs have the wrong trends in amplitude. This effect remains pronounced (lower left), even if only one zero is moved off the critical line to form a pair Ī0.1 off, consistent with the symmetry about the critical line of (see fig3h for expanded view). (Lower right) Shifting the real parts alters the amplitude trend of successive prime fluctuations in the summed terms (0.2 blue, 0.5 black, 0.8 cyan).
When we examine the formula for the summed wave functions:
we can see that , the increased amplitude of the summed terms at primes as shown in fig 2d depends only in the imaginary values of the zeros, as is confirmed from the above formula for the summed conjugate pairs .
By contrast, in the terms in the explicit formula:
the leading term xp, which is inverse quadratic for x=1/2 is pivotal in the sum of the wave functions remaining constant on intervals in the zeta prime distribution. Explaining qualitatively why this is so, would pertain, not just to the first few zeros, but to all zeros, solving RH.
Fig 2d: (Right) When the terms are summed, they have bursts of amplitude at the primes, forming a duality to the Fourier transform in fig 2b, even when the real parts are off x=1/2 (0.2 green, 0.5 black, 0.8 yellow). Shifting the real part off x=1/2 affects their relative amplitudes at each x but not their peaks at the primes and prime powers. (Left) However when even a single zero is moved off the critical line to form a symmetric pair consistent with the symmetry of , the convergence to a constant function on non-primes is disrupted by a rogue amplitude. The effect remains pronounced even when the zero is the 400th one.
While the Fourier sin transform of fig2b differs from the process forming the zeta zeros, so can't easily be used to form a proof about the locations of the zeros, and the product formula equals the zeta function only for real(z)>1 invoking the same problems, the explicit formula is a precise result of the zeta-prime relationship which was proven by Riemann himself to be a step square wave function on primes and prime powers. Therefore showing the constraints necessary to produce a step square wave superposition in a similar manner to determining the coefficients of a square wave Fourier series would require the real parts to be on the critical line would establish RH. On the other hand, this might require the consequence of RH on the prime distribution to be true, to prove it.
A second view of the other end of the zeta transform can be seen from examining the Euler product itself. In fig 3 are the truncated products for primes less up to 2, 3, and powers of 10 up to a million terms. Each of the terms in the product formula is finite in the critical strip except for the pole at 1. Since , each term is also periodic with period with minima on the critical line x = ½ at .
Fig 2e: (Top left) Fourier sin transform of integer step function (black) overlaid on that of the prime counting function (green) and actual zeta zeros (blue). (Top right) Step function minima distribution blue closely coincides with zeta zero distribution (green). (Lower right) Difference between the two. (Lower left) Applying the explicit formula to the step minima at real part 1/2 results in a function with integer local maxima with neighbourhood peaks at primes. (Inset) The step function and summed Mangoldt prime distribution compared.
Part of the intrinsic difficulty of RH is that we are trying to compare an irregular transform consisting of zeta zeros with the irregular prime distribution. Thus to a certain extent the zeros and primes are mutually encoded, so that it is difficult to establish RH without knowing the prime distribution and vice versa. To look for a comparative test function with more regular properties, we now examine the Fourier sin transform of the integer step function as shown in fig 2e, which I discovered accidentally has intriguing properties. Attempting to construct an eta-like associated Dirichlet function for the zeta analogue of this distribution, which would be convergent on the critical strip, and thus reveal its exact zeros, would be very difficult, because the k-th coefficient, instead of being 1, would be the total number of possible factorizations of k, which may grow faster than the fixed powers of k in the denominator. The Fourier sin transform gives a smoother transformed function than the Mangoldt psi prime counting function, and its minima can be easily found. Bizarrely, the minima are irregular and correspond closely with zeta zeros, with an overall distribution closely following the actual zeros up into the thousands. When the explicit formula generating the prime distribution from the zeta zeros is applied to the transformed step minima, we gain an ascending function with integer local maxima and neighbourhood peaks at the primes, subtly combining both features. The fact that we gain a meaningful distribution at all suggests there should be an analytic continuation like the zeta function underlying the duality.
Fig 3: Convergence of the sum formula in the critical strip (from left) is compared with that of the product formula (from right) in the range t =10-40 for the absolute value of . The pattern of the zeros is clearly manifest in the multiplicative superposition of the terms in the product formula in which the zeta zeros form the combined minima, although generally not coinciding with the minima of any one of the prime terms.
The narrowing of gaps between the zeros for increasing t thus results entirely from the phase relationships between successive prime periods centered on the x axis. For any finite product of the positions of the zeta zeros on the critical line neither represent minima nor zeros of the finite product, which declines to zero along the line where is a zero of , but grows to a peak for intermediate values before declining to infinitesimal values for negative x.
While the sum and product formulas for large finite values closely approximate one another for x > ½ the behavior of the finite product formula even for p in the millions becomes increasingly different from , for x < ½, with a set of escalating peaks and troughs increasing slowly in number with , so that even for primes up to a million there are still only a small number of these peaks and troughs between any two zeros of, causing a significant deviation from , even for exponentially large p.
The lack of a proof of the Riemann hypothesis doesn't just mean we don't know all the zeros are on the line x = 1/2 , it means that despite all the zeros we know of lying neatly and precisely smack bang on the line x = 1/2 , no one knows why any of them do, for if we had a definitive reason why the first zero 1/2 + 14.13472514 i has real value precisely 1/2 we would have a reason to know why they all do. Neither do we know why the imaginary parts have the values they do.
So let's have a look at the dynamics of the product formula in the critical strip, noting that the equivalence between this and zeta strictly holds only for x > 1. We can see immediately that on the critical line, the product formula has become unstable and doesn't converge consistently to 0. If we look at the first zero, we can see some very odd things are happening.
Fig 3b: The orbit of the first zero for increasing terms in the product formula initially has erratic values for the first few primes (a) driven by prime fluctuations, but as the primes increase enters into an unstable orbit which grows to a peak (b), with successive higher peaks which return to values close to 0 (c). Looking at the absolute value (d) of individual terms (green) and the cumulative product (blue), we can see that oscillating trends in the prime wave functions are now causing exponential bursts in the cumulative product. Viewing this in terms of the log of the product shows an approximate quadratic growth in the amplitude as the number of terms increases.
As shown in fig 3b, the initial erratic prime-driven values of the product, begin to settle into a more regular trend, which leads eventually to regular oscillations in the terms resulting in exponentiating peaks tending to unbounded values while intervening values tend to 0. We will thus get differing answers for the limit depending on how the product terms are grouped.
As we move up the zeros, this process becomes more complicated with a vastly longer sequence of erratic steps forming a fractal orbit resembling a stochastic process.
Fig 3c: Two higher zeros showing extended erratic transients in the cumulative product. Fractal expansions of the orbit of the zero with y=121412.
By the time we reach values as high as 121412, the erratic transients have become so long that it is unclear that for realizable values of primes, the orbit pattern will have settled into quasi-regular oscillations.
An indication of how high the values of primes would have to be to see any resolution of the orbits of larger zeros can be seen from fig 3e, where even a zero as small as 523 takes until primes of the order of 3 million to begin to enter the oscillatory phase.
Fig 3d: When we evaluate the cumulative product up to the 1,642,052th prime 26299991, we find the first zero y~14 (top) has grown to a peak of around 10 million, while the zero y~523 (middle) has only begun to enter its first oscillatory burst around the 200,000th prime of around 3 million and y~121412 is as yet showing no signs of having fully explored its fractal dynamics.
The trend to exponential bursting is portrayed in fig 3e, where a complex plot of the values of the product and zeta are compared for differing numbers of cumulative terms.
Fig 3e: Cumulative product and zeta compared at the first zero for 100, 2513 and 84270 prime terms. The value of the product at 84270 is not 0 but 6.25. Exponentially larger values occur at successive peaks, as shown in figs 3b, 3d.
Now we turn to looking at the additive representation, where each of the added functions f(z)=1/(kz) is a complex exponential varying exponentially in modulus with x and sinusoidally in angle with y. The effect of additive superposition is akin to a wave superposition of each of these functions. Since the eta representation of zeta is valid in the critical strip, we can easily investigate how rapidly the zeta function converges to its zeros, as the number of iterations is increased. When we do this we find that for a given zero, there are a series of erratic transients followed by a steady winding toward zero past a critical iteration number. However as we go up to higher zeros, these transients become longer and longer and involve winding in and out of a series of non-zero values resulting from the combined wave interference of the successive iterates. In fig 3f we can see just how complicated this process has become for the 2000th zero.
Fig 3f: Convergence of zeta to 0 at the four zeros starting at the 10000th zero shows erratic transients involving winding in and out of around 20 separate false zeros before eventual convergence. There is no consistent trend between successive zeros.
Video of the "Dance of the Zeros" - the sequences of the 401 zeros from the 99,600th to the 100,000th done for Eta(z), as it gives more consistent convergence.
Video of a smooth transition from the 99996th zero up the line x=1/2 to the 100000th.
Although it is tantalizingly obvious that all the zeros lie enticingly on the line x=1/2, it would appear extremely unlikely that a general argument can explain why ALL of them do because the spiraling behavior involves the real and imaginary parts equally and the imaginary parts of the zeros form an uncomputable sequence of values. Similarly although we know every Collatz sequence for positive integers ends in a 3 cycle, proving all do is unsolved. If we accept we can‚Äôt find a formula for the roots of a degree 5 polynomial, thanks to Galois, expecting to find those of zeta, despite its apparent symmetry may be a delusion.
We can make an estimate of how rapidly zeta converges to zero for successive zeros, we do this we arrive at a highly erratic relationship for individual zeros. This proves to be a result of the way zeta is expressed in terms of eta, as revealed in the smooth trend for eta with a power law of ~13x0.68.
Fig 3g: (Above) The number of iteration steps in the eta-derived zeta series required to get 5 steps with 0.005 of 0 varies erratically from one zero to the next, but this is a disguised effect of the presence of the 1/(1-21-z) term so becomes a smooth curve for eta (below).
We can also explore the source of the false convergences to non-zero values forming the spirals in the above sequences. As can be seen from the next image plotting the real and imaginary components and the angle and modulus of the individual terms, major shifts of convergence coincide with interference effects, when a number of terms in sequence have a similar angle due to constructive interference in the waves of angle in the imaginary direction for each term in the series, thus contributing a systematic shift in values, while steady trends in angle tend to cancel. Close to these values, a change from n-z to (n+1)-z causes a near perfect frame shift of the angle to a multiple of pi later permitting a cumulative position shift. This also causes a problem for RH because of the lack of an obvious relationship between a sum of square roots n1/2 defined by x=1/2 and the logarithmic variation of the imaginary waves defined by Cis(yln(n)). The mode-locking phases in fig 3h can be calculated directly by finding where the waves match phase:
This corresponds closely to the series of the mode shifts (yellow) in real and imaginary parts (blue and red).
Fig 3h: Trends in real and imaginary parts of the sequence of terms converging to the 20000th zero are compared with thrends in the angle and modulus of the individual terms. Showing convergence shifts correspond to angular mode-locking.
Although the Riemann hypothesis (RH) has not been proved, all of the 1013 of zeros so far found in the range up to 1020 lie on the critical line. It has been established that an infinite number of zeros lie on the critical line, that over 40% of the zeros do, and that all but an infinitely small proportion lie within of x = 1/2.
However these results do not guarantee RH is true. Littlewoood who had been given RH as a PhD thesis topic by Hardy and proved an infinite number of zeros are on the critical line said  “There is no evidence whatever for it (unless one counts that it is always nice when any function has only real roots). One should not believe things for which there is no evidence.”
Fig 4: Above: t values for the zeta zeros. Below: the distribution of primes. A complementary relationship results from the superposition of the product formula terms, whose periodicities are . On the average, the zeros are distributed as and the primes, from the prime number theorem as .
We have seen that S(t) grows extremely slowly with t so that major fluctuations in the zeros might not emerge with the large numbers so far computed. Other properties of the zeta function, such as changes in the topology of real( )=0 between 121414 and 121416, shown in fig 5, emerge only with moderately large numbers. RH is equivalent to the conjecture that the prime counting function , where . In another classic result closely related to the zeta zeros, it has been proven that changes sign infinitely often, although the difference is negative for all calculated primes. Skewes’ bounds  for a change of sign of assuming RH (see fig 2b), and not assuming it, show such changes can occur far beyond numbers so far computed. Although lower computer bounds of 1.39822 × 10316 where there are at least 10153 consecutive such integers near this value without assuming RH have been established, these are still astronomical by comparison with the known zeta zeros, so further anomalies in zeta zeros could appear.
Fig 5: Left: Curves of (blue) Li(x) (red) and x/log(x) (green) showing the fit of the estimates.
Right: Region near t=121415 may be the first place where the curves wrap around each other.
The zeros have slight errors of position, off x = ½, due to the limit on the number of terms in the sum formula used.
On the other hand the Riemann hypothesis might be unprovable yet conditionally true, like theories of physics, such as relativity and quantum theory, which, if all acid tests to refute the theory’s predictions fail, remains valid until a counterexample is found under new physical conditions. If RH were found to be formally undecidable, demonstrating the inability to prove it false would be grounds to declare it true, as it would have been proven that no counterexample, i.e. offline zero, could exist.
Finally it might turn out to be both unprovable and uncertain, because it can only be resolved by a computation that suffers Turing machine halting undecidability. Whether many cellular automata will terminate, and the Collatz conjecture, have similar unprovability problems associated with unpredictable intermittencies associated with growth and decay. If RH proved unprovable on analytical arguments alone, such as the obvious internal symmetry of the functional equation manifest so obviously in , RH might simply prove to be logically equivalent to all its complementary formulations in terms of primes and number theory, so that it could only be proven if and only if these arithmetic results could be proved independently. For example, the prime number theorem, which was first proved based on the zeta function now has elementary proofs, showing the belief that such number theory results could be proved only by analytic techniques was unfounded. Conrey9 admits as much inclosing his review in Notices of the AMS:
A major difficulty in trying to construct a proof of RH through analysis is that the zeros of L-functions behave so much differently from zeros of many of the special functions we are used to seeing in mathematics and mathematical physics. For example, it is known that the zeta-function does not satisfy any differential equation. ... It is my belief that RH is a genuinely arithmetic question that likely will not succumb to methods of analysis. There is a growing body of evidence indicating that one needs to consider families of L-functions in order to make progress on this difficult question.
Nevertheless, although Conrey has placed his faith in the various L-functions which generalize zeta, associated with structures such as the elliptic curves pivotal in solving Fermat’s last theorem, and the families associated with finite fields, for which Andre Weil had proved RH, which can be associated with orthogonal, unitary and symplectic symmetry types, he notes an ultimate impasse:
“There is a growing body of evidence that there is a conspiracy among L-functions – a conspiracy that is preventing us from solving RH”.
This arises from the inability to eliminate a spurious zero near 1, the Landau-Siegel zero, which stymies predictions.
“The Riemann hypothesis is a very elusive thing. You may remember in Peer Gynt there is a mystical character, the Boyg, which bars Peer Gynt’s way wherever he goes. The Riemann hypothesis resembles the Boyg!”
“Did you know John Nash, the protagonist of the film ‘A Beautiful Mind’?” “I did, and I had enormous respect for him. He solved three very difficult mathematical problems and then he turned to the Riemann hypothesis, which is deep mystery. By comparison, Fermat's is nothing. With Fermat's - once they found a connection to another problem - they could do it. But the Riemann hypothesis, there are many connections, and still it cannot be done. Nash tried to tackle it and that's when he broke down.”
Another interesting indication of this impasse, which also highlights irregular features in zeta, akin to randomness, arises if we examine the L-function , where , is the Möbius L-function (appendix 4). An equivalent of RH is that , which would guarantee the Möbius function above would converge for x > ½, and show there were no poles (and hence no zeta zeros), however Mertens’ conjecture that was proved false and is in serious doubt. The above equivalence of RH to M(x) means is as likely to have a 1 as a -1, thus behaving essentially like a random function, which Chaitin  noted in considering whether RH might be unprovable. Denjoy’s probabilistic argument for RH is precisely this - comparing the spacing of the zeros with a random walk - that if is a random sequence of 1’s and -1’s then the simple random walk, with probability 1. Hence this says the parity of the number of prime factors of a number varies randomly. The problem as Terrence Tao has pointed out  is that the primes show both pseudo random and ordered behavior, making a proof difficult, because the process is then not able to be captured in a finite symbolic description.
Finally one fundamentally important universality property: somewhere in the critical strip the zeta function fits arbitrarily closely any smooth non-zero function in a small neighbourhood, . This is done by first approximating the function by a finite product of primes in the product formula and then showing the total product, i.e. zeta, comes arbitrarily close to this for a suitable neighbourhood of 3/4+iT.
A major breakthrough was thought to have happened when Montgomery made contact with the physicists studying nuclear energy levels at Princeton and found that the pair correlations in the gaps between the zeta zeros followed the same Gaussian unitary ensemble statistics as chaotic quantum systems and energy levels of large nuclei:
Montgomery was also taken to discuss the twin prime conjecture (appendix 3) with Gödel, who likewise worried that this might be undecidable. The GUE statistic and its time-reversible real variant, the GOE, or Gaussian orthogonal ensemble, appear in many forms of quantum system whose classical analogue is chaotic. The corresponding form for fermions, rather than bosons is the symplectic GSE. These include both the many body problem of nuclear energetics, highly excited atoms in a magnetic field and the quantum stadium problem. One of the greatest moments in this interaction of fields9,16 was the Keating’s development of a formula for the zeta moments, using characteristic polynomials of unitary matrices (see appendix 1).
Initially this correspondence between fields caused great excitation in both the mathematics and physics communities, and a number of eminent researchers tried, so far in vain, to prove RH by discovering a system of random Hermitian matrices whose eigenvalues would be real and might correspond to the zeros of thus showing they had to be real and hence those of would be on x = ½. However this program has so far not borne fruit, despite many concerted attempts8.
Fig 6: Left: The spacing between zeta zeros 100-10,000 is compared with a Gaussian Unitary Ensemble distribution (red), showing coincidence of the two statistics, and a Gaussian Orthogonal Ensemble (green) characteristic of the Wigner distribution of atomic nuclear energies. Inset: Pair correlations for the first 10^5 zeros compared with the theoretical GUE distribution. Right: Quantum dot stadium eigenvalues display both (a) GOE and (b) GUE forms of random matrix distribution under increasing magnetization of the electron’s orbit.
In attempting to create a convergence between Hermitian operators and the zeta function, researchers have constructed a variety of candidates, some very complex and others deceptively simple. Berry has presented one of the most straightforward of these, the semi-classical operator H=xp, and attempting to modify it to establish an operator having correspondence with zeta, demonstrating several putative connections between this and the zeta zeros. However the space on which this acts is not elucidated and the complex plane would need to be ‘sewn up’ in Berry’s own words into a region which makes the dynamics quantally bound. Secondly there is no elucidated relationship between the primes and the periodic orbits of the Riemann dynamics.
Connes has constructed a Hermitian operator whose eigenvalues are the non-trivial Riemann zeros. His operator is the transfer (Perron-Frobenius) operator of a classical transformation. Berry comments that such operators formally resemble quantum hamiltonians, but these usually have very complex non-discrete spectra with singular eigenfunctions. Connes gets a discrete spectrum by making the operator act on an abstract space here the primes acting on the Euler product are built in using a space of p-adic numbers and their units. The proof of the Riemann hypothesis is then transformed into establishing the proof of a certain classical trace formula.
Selberg has constructed a zeta function related to hyperbolic motion on constant curvature surfaces generated by discrete groups. The product formula is not over primes, but over all primitive periodic orbits for the motion of the surface considered.
where are the lengths of the orbits, and s is complex. This function like the Euler product is defined only for real(Z(s)) > 1, but can be analytically continued to the entire complex plane:
As a result, Z(s) has both trivial zeros at 1, 0, -1, -2 etc. and a set of non-trivial zeros putatively on the critical line x = ½. Z(s) has a similar trace formula to the Weil explicit formula for sums over the zeros of zeta. The correspondence between primes and periodic orbits, provides a correspondence between zeros and eigen-momenta in which ln(p) corresponds to the orbital period , resulting in an equivalent expression of the prime/periodic orbit number theorem:
Fundamental to the problem are two issues. The first is that the duality already seen in the relationship between zeta and the prime products is already the duality transform one is seeking. That is the system that decodes the zeta zeros is the distribution of numerical primes itself, so seeking an analogue from other mathematical areas cannot necessarily simplify the problem.
Secondly these GUE systems may show similarities in their statistics to zeta’s zeros, because they share overall features combining structured constraints and pseudo-randomness in common with the primes and zeta zeros, without necessarily being isomorphic to them. In a sense there is a regress occurring, in which attempting to model GUE systems to zeta’s zeros results in more elaborate mathematical constructions which share zeta’s characteristics but neither provide a breakthrough in proving RH, nor result in a real valued quantum operator.
The quantum stadium is a direct analogue of the classical chaotic stadium billiard which displays the classical butterfly effect of chaos - sensitive dependence on initial conditions - and for almost all orbits produces a dense trajectory filling the stadium as shown in fig 7 (a). Within this classical system is a dense set of repelling periodicities, any arbitrarily small deviation from which results in a dense orbit, or a differing periodicity.
The quantum versions of this system behave in a fundamentally different manner. While the initial stages of a trajectory follow the classical picture, after a limited period of time, called the quantum break time, they have a cumulatively increasing probability of entering one of the eigenvalues of the system. These eigenvalues turn out to correspond to the closed orbits of the classical system, which have now become probability maxima of the quantum system because wave spreading has effectively compensated for sensitive instability of the orbit, resulting in wave-periodicity and so-called scarring of the quantum wave function by probability maxima along these closed orbits, which also extend to fractal eigenstates of open chaotic systems.
Moreover, unlike the eigenvalues of ordered quantum systems such as the Lyman, Balmer and Paschen series of orbital electrons, whose energy separation converges to zero at infinity, the chaotic eigenvalues display energy separation statistically distributed as a GOE or GUE suppressing small energy transitions between eigenvalues. Semi-classical simulations of such systems, using a small classical wave packet, generally give similar results, showing the suppression of chaos and the separation of eigenvalues is directly caused by wave spreading.
In systems like the quantum stadium, the closed orbits and their eigenvalues are playing a role similar to the primes in that they are orthogonal or uncoupled to one another, are determined by constraints which result in a discrete spectrum and form an irrationally related subset of the phase space. Primes among the numbers behave similarly in that they have no common factors, form a discrete spectrum having no consistent rational formulation and act as a set of discrete generators of all the other integers. Thus the correspondence may be analogical but not homologous.
Fig 7: Quantum chaos: The classical stadium billiards is chaotic. A given trajectory has sensitive dependence on initial conditions. As well as space-filling chaotic orbits (a)  , the stadium is densely filled with repelling periodic orbits, three of which are shown in black in (d). Because they are repelling, neighbouring orbits are thrown further away, rather than being attracted into a stable periodic orbit, so arbitrary small deviations lead to a chaotic orbit, causing almost all orbits to be chaotic. The quantum solution of the stadium potential well (b)  and (d)  shows ‘scarring’ of the wave function along these repelling orbits, thus repressing the classical chaos, through probabilities clumping on the repelling orbits. A semi-classical simulation (c) shows why this is so. A small wavelet bounces back and forth, forming a periodic wave pattern, because even when slightly off the repelling orbit the wave still overlaps itself and can form standing wave constructive interference when its energy and frequency corresponds to one of the eigenvalues of a periodic orbit, even though the orbit is classically repelling. The quantum solution is scarred on precisely these orbits (d). This causes resonances such as absorption peaks of a highly magnetically excited atom (e) to coincide with the eigenfunctions of the repelling periodic orbits, just as the orbital waves of an atom constructively interfere with themselves, in completing an orbit to form a standing wave, like that of a plucked string. The result is that, over time, in the quantum system, although the behaviour may be transiently chaotic, it eventually settles into a periodic solution. Experimental realizations such as the scanning tunneling view of an electron on a copper sheet bounded by a stadium of carefully-placed iron atoms (f)  , confirm the general picture, although, in this experiment, tunneling leaked the wave function outside too much to demonstrate proper scarring. The semi-classical approach matches closely to the full quantum calculation (g).
The end result is that for a variety of closed quantum systems, wave spreading eventually represses classical chaos by scarring, causing the periodic eigenfunctions to become eventual solutions of any time-dependent problem, although the initial trajectory behaves erratically, just as does an orbit in the classical situation. For example, a periodically kicked quantum rotator  ,  will stochastically gain energy, just as in the classical situation, until a quantum break time  , after which it will become trapped in one of the quantum solutions. A highly excited atom in a magnetic field will have its absorbance peaks at the periodic solutions, and quantum tunneling will likewise use scarred eigenvalues as its principle modes of tunneling  ,  .
Evidence supporting differences between these two types of system comes from studies of the fractal dimension of the graph of zeta zero gaps for large zeros, which shows a Hurst exponent of 0.095 corresponding to a fractal dimension of ~1.9, with anti-persistence, indicating large gaps are followed by smaller ones, self-similarity over a wide range of values and significant differences from corresponding GUE systems. When corresponding block sizes of zeros and random matrices are used, Hurst exponents for the zeros and matrices are 0.34 and 0.65, suggesting fundamental differences in fractal structure.
Fig 8: Left: Pattern of differences between successive zeros of zeta in the range from 1021 for 104 successive zeros determined by Odlyzko  shows a Hurst exponent corresponding to a fractal dimension of 1.9 with pronounced negative persistence, differing significantly from corresponding statistics of GUE random matrices. Right: Julia set of zeta by ‘inverse iteration’ of the zeros  highlighting their superimposed first 20 successive inverse images, shows all zeros are in the basin of attraction of the fixed attractor . The positions of the 1st and 2nd non-trivial zeros can be seen as small annuli just to the right of the centre dark cleft in the first two Julia bulbs shown in inset.
Quicktime movie of bifurcations in the RZ Julia set
Quicktime movie of wave function method on RZ Julia highlighting the zeros in frame 1.
Quicktime movie of Gaussian wave function highlighting the non-trivial zeros in frame 3.
Matlab code for Riemann Zeta and Gamma functions as in Fig 36
Matlab code for Riemann Zeta and Gamma function Julia sets as in Fig 36
Matlab code for Riemann Zeta colour coded Julia sets as in Fig 37c,d
Matlab code for Riemann Zeta and Gamma function parameter planes as in Fig 37c,d
Fig 8b: (a) The Riemann Zeta function , showing its pole at 1 the trivial zeros at even negative real values and the non-trivial zeros on the line x=1/2. (b) The Julia set of , highlighting eventually fixed points in the internal basin mapping to , to which the zeros are also mapped. Inset in right overlap of eventually fixed point and non-trivial zero showing their proximity, with the eventually fixed point to lying on the curve where .
As illustrated in fig 8b, the Julia set of (Woon) forms the boundary between the basin of attraction of and the attracting fixed point . The first six non-trivial zeros of the function, from on, famous for the Riemann hypothesis - that they are on all on the line x=1/2 - lie in the basin of attraction of this fixed point. In fact all the zeros of the function do, including the trivial ones at , as all are mapped to 0, which iterates to .
Although the Julia set of c = 0 appears to be connected, this is not the case, because all the trivial zeros on the negative real axis also iterate to the attracting fixed point and they are accompanied by an infinite fractal series of complete island copies of the original Julia set caught in the V's in the negative real half plane (see right).
To test the question of parameter planes as a measure of the bifurcations of a family of Julia sets, we now examine the parameter planes of the function . The most outstanding critical value of is the value 1, to which all z with positive real values tend as . This however is a case of self-organized criticality, as it projects the entire positive half plane directly on to the neighbourhood of the pole at 1. The Julia set is thus extremely sensitive to small changes in c, undergoing explosions of the positive half-plane for arbitrarily small changes in c.
Fig 8c: (a) Parameter plane of iterating from the critical value 1 [-40,6] x [-40,40]. (b) Enlarged region [-20.5,-15.5] showing complex bifurcations. (c,d) A tiny fractal 'lake' barely visible around c = 0, shown expanded [-0.1,1] x [-1.1,1.1] (c) and larger (d) . This corresponds to c for which the positive half-plane diverges due to the pole at 1. Inset in (c) at half the scale is the corresponding 'lake' iterating from 0 corresponding to c for which the zeros of diverge. (e) Sample Julia sets on the real line, with colour chosen to distinguish divergence, convergence to the critical strip and convergence into positive half plane - shades of blue for points tending to , red for points asymptotic to a periodicity in the critical range, green asymptotic to a periodicity with real parts entirely positive, and grey to black for non-periodic points (largely absent). The range from -15 to -0.1 is similar to 0.6 with no apparent bifurcations. (f) An enlargement of the region shown # in (a). The sensitive changes in the Julia set from 0 to 0.0125 correspond to a transition off right of the island in (d).
In fig 8c we show the parameter plane for the critical value 1, in which the black regions indicate the critical value remains finite and colours indicate divergent iteration to , along with corresponding illustrations of Julia sets, which highlight both complex bifurcations on the real line in the interval [-21,-15], and explosions of the positive half plane on either side of 0, at ~-0.005 and ~0.001. Both these features correspond closely to the parameter plane, which shows both complex fractal structure in the former range and a tiny fractal 'lake' around 0 connected by a dendritic thread of further 'lakes' winding in the imaginary direction to the divergent region.
Fig 8d: Shading the bulbs demonstrates their periodicity, as confirmed by the Julia set portraits, which display the corresponding rotational periodicities in their spirals , confirming an upward set of odd periodicities 3, 5, 7, 9 ... and a downward set of integer periodicities 3, 4, 5, 6 ..., each with mediant fractality viz (3,4)=7, (4,5)=9, (3,5)=8, (5,7)=12.
In fig 8d we demonstrate that the same mediant winding sequences appear on the bulbs of the black region as in the standard Mandelbrot set of polynomial and transcendentals functions. The right hand plane coloured black in fig 37c(a) is iterating to the fixed point 1 + z because - = 1, so + z = 1 + z. But this is true for all z with large positive real part, so the iteration is fixed on 1 + z. We can thus colour the black region according to how many steps it takes to reach within epsilon of a fixed point or periodicity and colour by the number of steps in blue and add redness for the period. This immediately shows up the periodicities of the bulbs neighbouring the boundary which can be confirmed to correspond to Julia sets with rotational periodicity the same number, confirming the sequences of periods of the bulbs and the mediant relationship in the fractal progression.
Intriguingly the Farey tree of mediants appear in one variant of RH. Farey sequences consist of all fractions with denominators up to n in order of magnitude ‚Äď viz . Notice that each fraction is the mediant of its neighbours (i.e. ).
Two versions of RH state Farey sequences are as regular as possible  :
Fig 8e: Examination of the region [-16,-15]x[-1,-2], marked * in fig 8c(b), shows the role of the parameter planes on the 'critical' values 0 (taken by iterating from the zero at -2) and 1 (taken by iterating from 999) classify the Julia sets of . For comparison the critical saddle at ~ -15.34 with critical value 0.5206 is also plotted in this domain, displaying a classic period multiplying bulb. In the centre are shown the local parameter planes for these three values. Top left: The two parameter planes overlapped to show unions and intersection with locations of Julia set parameters 1 - 10 (A is a computational artifact). Only 0 and 1 show classifying attributes, with c values in displaying connected central regions (red) containing the attracting fixed point for the zeros. By contrast, c values in display bounded periodicities in the positive half-plane (in this domain red). Values of c in both sets display both features, those in one, one feature, and those in neither, neither feature, confirming the classification.
A second parameter plane iterating from the value 0, which, although it is not strictly a critical value is the value of all the zeros of f. This provides a similar overall profile and a fractal 'lake' centered on 1, corresponding to the bifurcations as the zeros cross the pole at 1 and escape the internal basins of the Julia set, which has already become disconnected into an infinite set of connected islands. We investigate the relationship between these further in fig 8e.
Fig 8f: Examination of the region [-16,-15]x[-1,-2], marked * in fig 37c(b), The Mandelbrot set of ¬†from the critical saddle at ~ -15.34 with critical value 0.5206 displays a classic bulb. Selecting a point in the period 3 subregion generates a beautiful example of period 3 Douady rabbits in the corresponding Julia set, confirming the relationship between the two and the capacity of the zeta function to locally model any function.
A second parameter plane iterating from the critical saddle at ~ -15.34 shows how the Julia peroperties of the zeta function can follow the same rules as our original quadratic and finally demonstrate the universality of the dark heart and its bulb structures as an atlas of Julia dynamics.
Fig 8g: Left Mandelbrot set of using the critical point z = 1.3828 + 42.2909i provides a period 3 Julia set at the location -0.375467 + 42.981618i starred in the period 3 sub bulb.
Finally the jewel in the crown! When we go to z = -17.3739, the first critical point on the x-axis with negative critical value -3.7436 we find, sitting in the inner bay of the central area, a perfect quadratic black heart, almost exactly as the quadratic heart of the Mandelbrot set of f(z)=z2+c does. These hearts are perfectly replicated in the fractal bays of all the dendrites up and down the critical strip.
When we go further negative to the first critical point with negative critical value -3.7436 at z = -17.3739 we find, sitting in the inner bay of the central area, a perfect quadratic black heart. The companion paper “Fractal Geography of the Riemann Zeta Function” shows the dark hearts of every type of critical point of zeta, both real and unreal for both additive and multiplicative parameter planes.
The real critical points vary from the miniscules through the transitional critical point at ~ -15 to the exponentiating vast criticals, each with their own Mandelbrot set and fractal dark heart satellites and the unreals in varying positions to the right of the critical line (see zeros of dzeta fig 15) each have their own set of dark hearts which contribute to the Julia dynamics.
Fig 8h: The quadratic heart marked by the black star for the critical point at -17.3739, with a
'Douady rabbit' period 3 Julia kernel chosen by selecting a point in the period 3 bulb.
The parameter plane of provides further source for investigation, as it preserves the zeros of zeta. As can be seen in the following figure, the zeros correspond to preimages of 0 in the fractal copy of the central basin accompanying each major dendrite in the multiplicative parameter plane. This is because the iteration at the zeta zero c maps and on to subsequent c iterates of 0. The zeta zeros and their pre-images are the only points eventually mapped to the origin, so they are the points forming the foci of the fractal pre-images of the main cleft. However their individual dynamics has no regular pattern because the origin is then iterated on in a manner unique to each c.
Fig 8(f): Multiplicative parameter plane showing relationship between zeros of zeta and boundary bays on the parameter plane forming pre-images of the central bay about 0. As we move up the zeta zeros to 134, their relationships become complex (near right), with some having dynamics enclosed by that of neighbouring ones. As the the imaginary values become larger at 410 (far right) several zeros may be enclosed.
Another dynamical view of the zeta function is provided by taking the function as the values of a vector field, the so-called 'holomorphic flow' Broughan and Barnett (2003), in which successive non-trivial zeros have a complex encoding in terms of sources, sinks and rotational flow, without any decoded pattern to date.
Fig 8g: Holomorphic flow (left) in the region between zeros at 282.4651148, which is a weak sink, and 283.2111857, which is a source, is reflected in the multiplicative parameter plane (centre) with the zeros superimposed, showing the lower one is enclosed, while that of the upper zero is part of a main dendrite spanning several zeros. The neighbourhood of the zero faithfully contains a pre-image of the main cleft surrounding the origin, further enlarged lower right, which demonstrates the same dynamics as the full-size multiplicative Julia set of the zero.
Fig 8h: (Above) Holomorphic sinks in the first 500 zeros all show a similar enclosed basin profile in the multiplicative parameter plane
contrasts with irregular quantitative variation in source/sink strength in the range up to the first sink (below).
Significantly the multiplicative parameter plane differs fundamentally from the additive one in that the major dendrites, rather than following the irregular and increasingly close-knit pattern of the zeros, have a regular frequency of around 9.04.
Fig 8i: While the dendrites of the additive parameter plane (a) [around 300i] follow the pattern of the zeros (b) and become more closely spaced with increasing imaginary values (see fig 4), in the multiplicative parameter plane (c) the zeros are absorbed into a regular pattern of large dendrites with a periodicity of around 9.04, which do not become more closely spaced at larger values (e) [121412i]. The reason can be seen by examining the function , which shares the zeta zeros, and also has the larger periodicity. We can see the relationship by examining how the parameter plane maps the 'critical' point z= 1000, viz , which coincides with for z=c. If this has a negative value, the point will escape into the unbounded region while if it is positive it will remain finite.leading to the periodic dendrties of the multiplicative parameter plane. The close correspondence between the function and the parameter plane of is emphasized in the overlap (f) of the two [around 300i].
The periodic behaviour of in this vertical region is confirmed in the next figure.
Fig 8j: (Left) A complex plot of the image under of the line passing from 0.5+250i to 0.5+350i shows irregular cycles through the zeta zeros. (Centre), the corresponding plot for 4+250i to 4+350i shows quasiperiodicity, with the real part regularly oscillating (Lower Right), rather than the irregular oscillations of the critical line (Upper right).
Sound (.wav) file of the regular and exotic oscillations, moving up from 3+250i to 0+75250i
A further means of exploring the zeta zeros is through the Newton's function of or :
Fig 8k: (Left: Newton Julia of zeta with trivial zeros red and non-trivial green. (Top centre) Basins of attraction of the non-trivial zeros as expressed by the Newton function of , with the basin of the third zero above and below highlighted in red. As can be seen in the enlargements from the centre (lower left) and one zero above and below (right) each fractal bud contains a complete fractal representation of all the non-trivial zeta zeros as is the case with Newton's cubic roots of unity. The blue regions in the centre of the features correspond to points jumping outside the escape boundary. They would reduce or disappear if the escape boundary were set higher. These points are not escaping and do still iterate to a zeta zero. In effect a fractal basin of every zeta zero is contained within each fractal feature.
We can thus portray the basins of attraction of the zeros and investigate the consequences for off critical line zeros in terms of the resulting fractal basin structure. As can be seen in the above figure looking at the basins of the Newton function for , the basins of each of the zeros are fractally reflected an infinite number of times in the boundaries of each. Given the reflective symmetry of , this would require a fractal connection along the critical line between the symmetrical pair of zeros.
A fundamental reason for the difficulty of proving RH, shared with a number of other unresolved conjectures arises from the irregular nature of both the prime and zero distributions. Analytic solutions to integration and differential equations problems frequently do not exist because the properties of the flow curves contain irregularities that prevent their properties being captured in a finite analytic formula. Many such systems can enter chaos and form complex systems at the edge of chaos through repeated bifurcation. Many problems in number theory and related areas also result in unproven conjectures, because the system is too irregular to permit a finite formula which can be proved to hold asymptotically. In such a situation the conjecture may become effectively undecidable because its proof would require computing all of an infinite set of cases – a computation that would terminate only if a counterexample to the conjecture appeared in the computation. Such a conjecture then becomes a kind of Turing halting problem.
The Collatz, Ulam, Syracuse, or 3n+1 conjecture  is a much simpler conjecture than the Riemann hypothesis, which nevertheless has resisted proof, despite the passing of the Poincare conjecture and Fermat’s last theorem. It seeks to prove that the iteration n a positive integer, always consists of a finite sequence converging to 1. Notably of this problem, Paul Erdös said “Mathematics is not yet ready for such problems”, so why mathematicians should expect RH to be solved outright, as a sine qua non, is a little puzzling.
The difficulty again is that, while every integer appears to eventually end up in the period 3 cycle some numbers take a very tortuous route to get there, involving alternating climbing and falling, with no established regularities to the pattern.
To understand why a simple e.g. inductive formula as not been found, which can extrapolate a solution for all positive integers, we can embed the integer iteration into a wider real or complex number iteration:
Fig 9: Above the sequence for 27 has 113 steps rising as high as 9232, before descending to the period 3 cycle.
Below: Cobweb plot of the iteration embedded in the graph of an ascending cosine function.
The real version of this iteration on integers is shown in fig 9 for n = 27 demonstrating it has the same effect as the 3n+1 iteration. We will examine the complex dynamics further below. The reduced sequence using (3n+1)/2, which shortens the 3n+1 sequence one unit at every odd integer, as 3n+1 is always even, will give a similar process.
Fig 11: Above: The 3n+1 problem appears to always end in the period 3 solution 4>2>1 for all positive integers but has three different periods of 2, 5, and 18 for negative integers. Below: Scatter plot of iteration numbers shows the Moire-pattern like structured irregularity that has given it the name of the ‘hailstone’ sequence.
Although the 3n+1 sequence appears to consistently converge to 4 > 2 > 1 for all positive integers and does have a general trend of largest iterates m(n) growing on average with n2, greater and greater numbers occur, forming path records which significantly exceed this estimate as shown in fig 10. The current highest known, the 88th path record is 1,980976,057694,848447 which reached 64,024667,322193,133530,165877,294264,738020 before eventually entering the 4 > 2 > 1 cycle. The successive path records 2(2), 3(16), 7(52), 15(160), 27(9232) etc. vary erratically along:
Although all positive integers appear to have eventual period 3, negative numbers are distributed erratically, with three different periods 2, 5 and 18. For positive integers, it has been proved that there are no non-trivial cycles with length less than 275,000, but the problem remains unsolved.
Conway has proved that Collatz type problems can be formally undecidable. The Turing halting issue is particularly troublesome, because a serial computation over all integers will continue to compute if either all integers have period 3, or the process strikes a divergent, or non-periodic sequence unless the process is arbitrarily cut off, in which case it may miss a new path record, and will terminate for certain only if it discovers another periodicity.
A hint of why this problem may be unprovable comes when we embed it in the complex plane and consider the analytic discrete dynamics. In this situation, it becomes immediately apparent that, while all the positive integers lie within basins of attraction of the attracting period 3 cycle to which the positive integers are eventually periodic, the negative periodicities have the modulus of their derivative greater than 1 and are thus repelling. The negative integers are thus in a state of so-called self-organized criticality, in which they are eventually periodic to an unstable repelling cycle on the boundary of the Julia set for which any small analytic perturbation would cause a catastrophic bifurcation of the orbit as shown in fig 12.
While the integer 3n+1 process can be neither said to be stable nor unstable, as it is logically confined to discrete integers, the complex analytic discrete process is clearly unstable. Thus while an elementary or inductive proof that all positive integers converge might be plausible, this is logically connected to the more highly disordered negative sequence, which looks a more highly resistant entity to proof.
Fig 12: Top the Julia set of the ascending cosine function modeling the 3n+1 iteration, showing positive integers lie within basins of attraction of the attracting period 3 cycle. Even integers occur in basins of attraction of the attracting period 3 cycle, coinciding with periodicities of the cosine function, while odd integers occur in smaller offset basins of attraction because the odd periodicity of the cosine is repelling. By contrast, negative integers occur in eventually periodic orbits leading to an unstable repelling periodic cycles on the boundary of the Julia set, for which any -small perturbation would cause the iteration to either go to infinity or to another attracting basin, or to a non-periodic orbit in the Julia set. Second row show the positions for -17 (period 18) -1 (period 2) and 27 (113 steps to period 3). Centre shows the offset of the basin for 27 and the enlarged location (left) of -1 on the boundary of the Julia set. At right is the corresponding location for the associate ascending cosine if (3n+1)/2 is used instead of 3n+1, and at bottom is the corresponding Julia set for this function, which shows equivalent behavior. The large central bulbs and their fractal replicates are the basin of attraction of the fixed attractor z = 0.
To see how severe this unpredictability can become for arithmetic functions, we only need to turn to the aliquot sequence,  , in which a number is iteratively replaced by the sum of its proper divisors (see appendix 4). Although it is plausible that all numbers eventually converge to periodic solutions, it has not been established that none have infinite aperiodic (e.g. divergent) sequences. Sequences of individual numbers such as 138 vary erratically, growing for a time, then apparently shrinking to much smaller values, before expanding beyond the limits of computation.
Most numbers end in a prime followed by 1 and 0, but others end in other periods: 1 (perfect number e.g. 6>6), 2 (amicable numbers e.g. 220<>284, as Pythagoras said true friendship was comparable to these), or higher (sociable e.g. 1264460 period 4). Eventually-periodic numbers are called aspiring. Sequences from different numbers may unite. The smallest perfect numbers are 6, 28, 496, 8128, 33550336.
Fig 13: Aliquot sequences for increasing integers. (Red): peak numbers exceeding Matlab’s integer bound of . (Blue) iteration lengths before periodicity. (Green) eventual period. Period 0 in the green graph also indicates overflow of the computational bound.
All perfect numbers so far discovered have the form , as noted by Euclid and proved by Euler. The 44th is which has 9808358 digits. Cycles with 4, 5, 6, 8, 9 and 28 members are known. Five numbers below 1000: 276, 552, 564, 660, and 966 have not yet been established to converge because they have exceeded current computational bounds. There are 81 below 10,000 and 906 below 100,000. About 1% of all integers are estimated to be beginning numbers (key numbers) of such a hypothetical open-end sequence.
Fig 14: Log size vs iterations graphs of 2856 and 980460 showing eventual periods 2 and 28.
2. Zeta Bernoulli Numbers
The first computer program using Babbage’s calculation engine was to produce the Bernoulli numbers, which also define real integer values of zeta13.
The twin prime conjecture - that there are always greater prime pairs - remains unresolved although Euler proved by elementary argument that there are an infinite number of primes (consider the product of all primes + 1 if this is not the case). The predicted distribution is:
Fig 14: Divisor function , totient function , Möbius function and Mangoldt function all have Dirichlet series sharing the zeta zeros as zeros or poles.
4. Other Dirichlet series for zeta
Initially I used the Lanczos approximation to the gamma function to extend the domain to :
with and calculated independently by Paul Godfrey as constants from the relation:
where are coefficients of the Chebyshev polynomial matrix. See Matlab toolbox for details.
A satisfactory alternative I have used subsequently is truncated to the same number of product terms as used in the sum formula for zeta.
7. A Formula for Depicting the Derivative of Zeta
Fig 15: The derivative of the Zeta function iterated to 1002 terms
One can easily use a derivative approximation for any of the functions investigated, however, here is a formal derivative of the key functions:
This can be extended to the derivative of Xi:
First Zeta Critical Points (dzeta approximation)
z = -15.3364, z(z) = 0.5206,
z = -17.3739, z(z) = -3.7436,
z = -19.4031, z(z) = 33.8083
z = -21.42790, z(z) = -374.418694 (uncomputable Inf overflow from here)
z = -23.44919, z(z) = 4988.005569
z = -25.46771, z(z) = -78673.950161
z = -27.48401, z(z) = 1449688.849563
z = 2.463354349136353 + 23.29765892028809i, z(z) = 0.9289 + 0.0308i
z = 1.286577939987183 + 31.7081241607666i, z(z) = 0.7073 + 0.0110i
z = 2.306453704833984 + 38.48979568481445i, z(z) = 0.9919 - 0.0931i
z = 1.382872581481934 + 42.29097747802734i, z(z) = 0.7930 + 0.1273i
z =-1.209025621414185+65.92226409912109i, f(z) = 0.7776 - 0.2061i
A short list of some zeta zeros mentioned in the dark heart research, which can be easily copied and pasted into the text fields are as follows:
1-12: 14.13472514, 21.02203964, 25.01085758, 30.42487613, 32.93506159, 37.58617816, 40.91871901, 43.32707328, 48.00515088, 49.77383248, 52.97032148, 56.4462477
125-130: 278.2507435, 279.2292509, 282.4651148, 283.2111857, 284.835964, 286.6674454
287-293: 523.9605309, 525.0773857, 527.9036416, 528.4062139, 529.8062263, 530.8669179, 532.688183
171382-171390: 121412.139210209, 121412.990421458, 121414.488895067, 121414.739043607, 121415.047364581, 121415.640550747, 121416.302522095, 121416.823543637, 121417.618749154
282.4651148, 391.4560836, 446.8606227, 527.9036416 H,
637.3971932, 653.6495716, 681.8949915, 762.7000333
For further zeros, see:
 Conrey J B (2003) The Riemann Hypothesis http://www.ams.org/notices/200303/fea-conrey-web.pdf
 Notices of the AMS Feb 2006 p223
 Marcus du Sautoy, The Music of the Primes, Harper Collins, 2003.
 Liu B, Zhang G, Dai J, Zhang H (1998) Eigenvalues and Eigenfunctions of a Stadium-Shaped Quantum Dot Subjected to a Perpendicular Magnetic Field Chin. Phys. Lett. 15/9 628-30.
 Berry M, Keating H (1999) H=xp and the Riemann Zeros http://www.phy.bris.ac.uk/people/berry_mv/the_papers/Berry306.pdf
 Connes, Alain (1996) Trace Formula in Noncommutative Geometry and the Zeros of the Riemann Zeta Function C.R. Acad. Sci. Paris 323, 1231-6.
 Bogomolny E (2007) Riemann zeta function and quantum chaos arXiv:0708.4223v1
 Casati G, Maspero G, Shepelyanski D (1999) Quantum Fractal Eigenstates Physica D 131 311-6.
 Alisa Bokulich (2008) Can Classical Structures Explain Quantum Phenomena? Brit. J. Phil. Sci. 59 217–235
 Gutzwiller, M.C. (1992). Quantum Chaos. Scientific American 266 78 - 84.
 Heller E, Tomsovic S (1993) Postmodern Quantum Mechanics Physics Today July 38-46.
 Moore F, Robinson J, Bharucha C, Sundaram B, Raizen M, (1995) Atom optics realization of the quantum δ-kicked rotor Physical Review Letters 75/25 4598-4601.
 Raizen M, Moore F, Robinson J, Bharucha C and Sundaram B (1996) An experimental realization of the quantum δ-kicked rotor Quantum Semiclass. Opt. 8 687–692.
 Berry, M (1989) Quantum physics on the edge of chaos New Scientist 29 Oct http://www.fortunecity.com/emachines/e11/86/edgechaos.html
 Wilkinson P, Fromhold T, Eaves L, Sheard F, Miura N Takamasu T (1996) Observations of 'scarred' wavefunctions in a quantum well with chaotic electron dynamics Nature 380 608-610.
 Monteiro T, Delande D, Connerade J (1997) Have quantum scars been observed? [+reply] Nature 387 863-864.
 King C (2009) Exploding the Dark Heart of Chaos http://www.math.auckland.ac.nz/~king/DarkHeart/DarkHeart.htm
 J. B. Conrey, D.W. Farmer, J.P. Keating, M.O. Rubinstein & N.C. Snaith (2005), Integral moments of L-functions, Proc. London. Math. Soc., 91, 33-104.
 J.B. Conrey, D. Farmer, J.P. Keating, M.O. Rubinstein & N.C. Snaith (2008), Lower order terms in the full moment conjecture for the Riemann zeta function, J. Number Theory 128, 1516-1554.
43. Odlyzko, Andrew (2002) (PDF), Zeros of the Riemann zeta function: Conjectures and computations