Exploring Quantum , Classical and Semiclassical Chaos in the Stadium Billiard

This paper explores quantum and classical chaos in the stadium billiard using Matlab simulations to investigate the behavior of wave functions in the stadium and the corresponding classical orbits believed to underlie wave function scarring. The simulations use three complementary methods. The quantum wave functions are modeled using a cellular automaton (CA) simulating a Hamiltonian wave function with discrete (square pixel) boundary conditions approaching the stadium in the classical limit. The classical orbits are computed by solving the reflection equations at the classical boundary thus giving direct insights into the wave functions and eigenstates of the quantum stadium. Finally a simplified semi-classical algorithm is developed to show the comparison between this and the quantum wave function method. Fig 1: Evidence of scarring: Superposition of the first 3300 iterations of a wave function initiated by two different generating functions (a and c in fig 4) with wavelength chosen to be a fraction of the associated classical repelling orbit. Upper two figures show superposition of phase amplitudes and lower two the corresponding superimposed probability distributions of the squared amplitudes, showing evident scarring. Generally the probability distributions give better resolution of the scarring. Introduction: Estimation of the quantum wave functions in chaotic systems, from the nucleus of atoms from helium to uranium to the eigenfunctions of a wave-particle in the Bunimovich stadium (the classical Coliseum stadium shape consisting of a rectangle or square capped off with semi-circular discs) proved to be initially intractable and even in the post-modern era of revived semi-classical methods computationally complex to the point of being counter-intuitive. Part of the purpose of this paper is to simplify understanding these quantum solutions so their dynamics can be more directly appreciated and understood. A full quantum solution to a quantum chaotic wave function, given its boundary constraints proved illusive for the founders of quantum theory in the case of the helium nucleus. In 1970 ii Gutzwiller iii developed a semi-classical method using a trace formula, which integrates the Green s function of wave propagation over all coordinates so that it can be expressed in terms of the repelling orbits hidden within the classically chaotic version of the quantum system. The trace formula for a particle of mass m and momentum p = 2mE inside a box with arbitrarily shaped walls is g(E) = 1 E En = g(E) + 1 2 ih n ml


Introduction
Estimation of the quantum wave functions in chaotic systems, from the nucleus of atoms from helium to uranium to the eigenfunctions of a wave-particle in the Bunimovich stadium [1,2] (the classical Coliseum stadium shape consisting of a rectangle or square capped off with semi-circular discs) proved to be initially intractable and even in the post-modern era of revived semi-classical methods computationally complex to the point of being counter-intuitive.Part of the purpose of this paper is to simplify understanding these quantum solutions so their dynamics can be more directly appreciated and understood.A full quantum solution to a quantum chaotic wave function, given its boundary constraints proved illusive for the founders of quantum theory in the case of the helium nucleus.In 1970 Gutzwiller developed a semiclassical method using a trace formula, which integrates the Green's function of wave propagation over all coordinates so that it can be expressed in terms of the repelling orbits hidden within the classically chaotic version of the quantum system [3,4].
The trace formula for a particle of mass m and momentum p = √ 2mE inside a box with arbitrarily shaped walls is where ḡ(e) is a smooth function giving the mean density of states, summed over all orbits γ of length l γ with v γ the phase shift counting the focal points and twice the number of reflections off the walls and over all retracings k of these orbits [5].The stability matrix M γ records the sen- Quantum chaos: (a) the stadium is densely filled with repelling periodic orbits [9], three of which are shown in black in (d).The quantum solution of the stadium potential well (b) [4] and (d) [5] 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.The quantum solution is scarred on precisely these orbits (d).This causes it 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.Over time, although the behavior may be transiently chaotic, the quantum system eventually settles into a periodic solution or a linear combination of these.
Experimental realizations, such as a magnetized electron in a quantum dot (e) [10] [Gaussian orthogonal ensemble (GOE) and Gaussian unitary ensemble (GUE) statistics], or the scanning tunneling view of an electron on a copper sheet bounded by a stadium of carefully-placed iron atoms (f) [11], confirm the general picture, although quantum tunneling in the latter case leaked the wave function outside too much to demonstrate proper scarring.The semi-classical approach matches closely the quantum calculation (g).
sitivity of the orbit to changes in initial conditions.From the left hand side of the equation it can be seen that if the sum converges (this is generally possible only with some difficulty and depends on reordering the sum terms), the eigenvalues will appear as singularities in the sum.The general similarity of the trace formula to the Riemann zeta function and this function's Gaussian unitary ensemble (GUE) type statistics shared by quantum chaotic systems (see panel (e) in Figure 1) has led to extensive attempts to solve quantum chaotic systems through the paradigm of the zeta zeros (for a detailed discussion on the quantum chaos connection with the Riemann zeta function see section 7).Both quantum and semi-classical approaches have been used to model systems such as the stadium billiard.In the quantum approach a wave packet is propagated throughout the potential well using a time dependent Schrödinger equation via Fourier transform techniques.This wave propagation approach is discretely modeled in the cellular automaton described in section 2. Alternatively the semiclassical approach uses the least repelling periodic orbits in the classical stadium to generate amplitudes based on the path lengths of all classical trajectories connecting a given location to another in the region.
Heller and Tomsovic [5][6][7] among others have used the semi-classical orbit-based approach to develop the theory of scarring of the quantum wave [8].The authors comment: Accurate excited-state eigenvalues have been computed from knowledge of relatively few periodic orbits.However, deep questions about the convergence of the trace formula and its modifications remain unanswered.[5, p.40]A particularly complete expose of the methods involved is provided in their Physical Review E paper [7].

The Quantum Cellular Automaton
To get a more direct picture of the process of wave spreading and wave function scarring [see panels (a-d) in Figure 2], this simulation depends on a simple wave propagation formula for a cellular automaton which preserves the essential nature of Hamiltonian dynamics and the Green's function governing transmission of wave amplitudes from point to point.Somewhat paradoxically, although the local rules are defined on a rectangular grind with only N, S, E and W neighbors, the process well models wave transmission at a variety of angles and under situations where the curved surfaces of the potential well boundaries refocus the wave fronts continuously in various ways.Again this is possible although the local rules at the boundary provide for exclusively horizontal and vertical reflection.
In effect, the wave period and implicit frequency of the generating wave gives the energy and momentum of the transmitted wave-particle and the individual pixels on the boundary provide a simulation of the smaller quantum features at the atomic level occurring in any real physical implementation of the quantum stadium, as exemplified in even more rudimentary discrete form, in the atomic sized stadium shown in panel (f) of Figure 1.
The cellular automation iterates the wave function according to the rule: (2) As in panel (e) of Figure 2, black positions whose four neighbors are white are all iterated in one step, followed by the white positions iterated in terms of their newly updated black neighbors.This provides a discrete iteration modeling a Hamiltonian process resembling action-angle variables with an angle shift of π 2 .
In keeping with the quantized nature of the system, the boundary of the region consists of a discrete series of vertical and horizontal edges corresponding to a discrete region of square cells approaching the curved stadium shape in the infinite limit.Vertical and horizontal edges have a reflecting rule according to which the relevant internal boundary cells are doubled and the external out of region boundary cells are omitted.The rule is illustrated by the following edge and corner situations: The process is also capable of handling applied forces by integrating their effect: The simulation uses a variety of generating functions to set up excitations whose period corresponds to a fraction of the orbit length of a classical repelling orbit in the stadium.Broadly they are of the type: combining one or more plane cosine waves of appropriate period lying along an orbit path clipped by being multiplied by a circular Gaussian, although circular and unidirectional traveling wave packets can also be established by clipping a sector of an emerging circular wave.Six of the scars resulting from superimposing successive probabilities measured as the square of the amplitude are shown in panels (b-g) of Figure 3.The predominant repelling orbits clearly show as probability peaks coinciding with the corresponding repelling orbits in just the manner of the semi-classical simulations in Figure 3 (see also [5]).In panels (a*-h*) are shown the corresponding generating functions for each of the images (a-h) in Figure 3.
The images in Figure 3 were generated using the first 3300 iterations of the generating function, however when the scars (b) and (f) were iterated on from 3300 to 6600 iterations the scar remained stable, as shown in Figure 4, verifying that it is not just a transient feature of the expanding generating functions.
These patterns of scarring can also be compared with other potential wells such as a rectangular one which leads to ordered wave functions and the lemon shaped potential well, which provides an example of so-called soft chaos in which ordered quasi-periodic solutions coexist with chaotic ones.In Figure 5 the generating functions of scars (a) and (b) in Figure 3 are applied again for 3300 iterations in the rectangular well resulting in ordered wave functions (a, b), and a vertical generating function in the lemon well produces an ordered wave function corresponding to the quasi-periodic stable orbit superimposed below (c, d).
Significantly the cellular automaton is able to well handle oblique generating functions with periods not reducible directly to an integer number of iterations.In Figure 2 are shown both the probability superposition scars and the amplitude superposition scars, confirming that even when an oblique wave with an irrational period is modeled, the cellular automaton is able to successfully superimpose amplitudes in phase over thousands of iterations.automaton can also well-model other wave experiments in physics such as the double-slit experiment illustrated in panel (a) of Figure 7.
It is also possible to compare the effects of different boundaries directly by choosing a circular wave generator of period a fraction of the overall dimensions (both dimensions for stadium and rectangle and the vertical dimension for the lemon).In panels (c-e) of Figure 7 the three boundaries are compared and the scarring characteristic of the stadium shows up clearly and quite elegantly suggestive of several superimposed orbit solutions.
The cellular automaton can also be easily interrogated to produce an autocorrelation function by multiplying the states of the generating wave function ψ 0 (x, y) by the propagated wave function ψ t (x, y) and integrating over all cells in the array.In panels (f, g) of Figure 7 are shown the autocorrelation function for the bow tie generator up to 3300 iterations consisting of about 7 sweeps across the length of the stadium and the fast Fourier transform, giving back the energy/frequency/period of the generated wave function.

Orbits in the Classical Stadium
To complement the wave cellular automaton and draw superimposed classical orbits, a ray drawing algorithm was developed for the same matrix of elements used in the wave cellular automaton algorithm.This also enables portrayal of several features of the associated classical chaotic system.In panels (h-k) of Figure 7 are shown how sensitive dependence on initial conditions (in this case a small angle difference in the initial orbit leads on multiple reflection in the curved sections of the stadium to exponentially increasing angular divergences between the orbits, even-Figure 6: A more detailed long-term evolution of scarring of a wave packet with twice the frequency and half the Gaussian variance for the same sequence of orbits as in Figure 3, superimposing stages: 0-3300, 3300-6600, 6600-9900, 9900-13200 for (f) and (g).The long-term evolution suggests that the systems are converging to a linear combination of the principal orbits due to the initial waveform exciting more than one orbital solution.Hence (b) and (e), being two of the strongest orbits, retain their excitation almost unchanged through 9900 iterations, while for example (g) already shows features of (b) in its initial excitations.Classical orbits in the stadium display the classic features of chaotic billiards.(h) Sensitive dependence on initial conditions illustrated for two closely related initial conditions (red and green) show exponentiating divergence with each reflection.Chaos also causes the space to be densely permeated with repelling periodic orbits, two of which are illustrated in (i) and (k).Because they are repelling, neighboring orbits are thrown further away, rather than being attracted into a stable periodic orbit as in an ordered system.Consequently almost all orbits are recurrent and densely fill phase space (j).
tually ending in essentially a random relationship with one another.Letting an orbit run for a large number of reflections also shows that most orbits are not periodic, but rather densely permeate the phase space of the stadium.Although the repelling periodic points are rare, they are dense which means there are periodic orbits arbitrarily close to any trajectory.Most orbits are recurrent so that they eventually come arbitrarily close to a given trajectory, as shown in panel (a) of Figure 9.
To more closely model the classical chaotic stadium, a second version of the ray-drawing algorithm was developed which rather than operating on a matrix of pixels, solves the reflection equations at the boundary on a classical continuous stadium.In turn this algorithm was refined to search for and display the closed orbits through a given point for reflection periods up to 20.Because these are testing numerically for a closed orbit by close approx-imation, they are really detecting close recurrence in a finite number of reflections, rather than an exact periodic solution.
As a dramatic increase is made in the number of reflections, such an orbit will eventually break stability, as the periodic solution is repelling, and it may over time recur closely to other repelling orbits.In panels (a-c) of Figure 9, the computed solution to the bow tie 4-period in panel (b) of Figure 3 is iterated for 200 steps in each direction, showing noticeable recurrences to another period 4 skewed bow tie and a rotationally symmetric period 8 orbit.
In panels (a-h) of Figure 8 sample orbits displaying reflective symmetry in the axes, rotational symmetry by half a revolution and asymmetric orbits are shown illustrating the variety of closed orbits and their symmetries.(c,  d), symmetric by vertical reflection (e, f), and asymmetric (g, h).Supplements 3-5: Quicktime movies of symmetric period 16 orbits passing through (0,0), (0,1) or (2,0).(i) Distribution of lengths of 6444 close orbit recurrences of up to 20 reflections generated from positions (0,0), (2,0) and (0,1), having rotational and reflective symmetry shaded to give relative densities.(j) Corresponding distribution of orbit length differences of the ranked length array.The right hand peak is an artifact caused by missing orbits with higher reflection numbers.(k) Distribution of length differences by frequency shows a descending law in which larger gaps are rarer.This pattern is repeated in the inset when the outliers are removed by choosing only the values from 2000 to 5000.This distribution differs markedly from the distribution of primes (see panel (e) of Figure 11) and is unlikely to be an artifact of the orbit selection because it is clear there are sequences of orbits with orbit lengths converging in the limit, such as those exemplified in (a), which tend asymptotically to the circumference of the stadium.According to the trace formula given by Equation 1 the phase shifts and sensitivity of the orbits also have to be included before a GUE distribution is arrived at, however, as shown in (l) a roughly GUE distribution, with prominent outliers corresponding to the shortest groups of orbits, does arise when the inverses of the orbit lengths are taken as a measure of energy.
back to the starting trajectory for 2 up to 20 reflections.The algorithm initially uses 80000 steps of angle through a sample right angle of directions from a given point and tests for local minima of a metric combining the distance from the initial vertex and divergence of angular direction.The local minima are then scanned down a further four orders of magnitude to find as close a recurrence to a closed orbit as possible given the numerical resolution of the floating point process.These solutions are then plotted as closed orbits and their lengths and sensitivity to initial conditions are determined as an initial step towards calculations based on the trace formula given by Equation 1.
The distribution of the lengths was also investigated as a comparison with the prime numbers whose separations increase and Riemann zeta zero distributions whose separations decrease with increasing imaginary value and whose step differences have the statistics of a GUE.Panels (i-l) of Figure 8 highlight the differences between the distribution of orbits and these two.
The orbit length profile shows some regions that are sparsely filled, and others which are more densely filled.
Certain small finite values of length, such as the perimeter distance of the stadium and its multiples have an infinite number of orbits of length tending to these values.The statistic of gap distances shows a descending power law trend, which is repeated internally, for example when the orbits from 3000 to 5000 are chosen, which, from panels (i, j) of Figure 8, is a relatively consistent central section of the orbit lengths not suffering from anomalies of small value divergences or missing large orbits due to cutting off the search at 20 reflections.The perimeter-angle coordinate system described in panels (d-g) of Figure 9 makes it possible to investigate the classical billiard as a discrete dynamical system.In this system orbits over a finite number of reflections become homologous in a symbolic dynamical representation in terms only of which curved and straight sections the orbits reflect off (e.g.rtblltrrb...).This makes it possible to reduce summations over all orbits, in forming classical or semiclassical correlation functions and wave plots, to summations over a single member from each family.The dynamic can thus be described in terms of the Markov partition of phase space determined by the orbits which strike or emanate from the four junctions at (±1, ±1), labelled in red, magenta, green and blue dots in panels (d-g) of Figure 9.

The Semi-Classical Approach
The semi-classical approach retrieves very comparable results to a full quantum representation, panel (g) in Figure 1, by tracing a key member of each family of orbits, in the partition of the classical billiard, whose momentum divides the length by the time p = l/t and applying an os-cillating wave propagator to the member of the form [7]: ×e ıS j (q,q ;t)/ −ıπv j /2 (6 where d is the number of degrees of freedom, S is the action, with j weighted according to the product of the orthogonal stability matrices: along the orbit, giving a correlation function for a single orbital family member as follows: e ı[S γ (q,q ;t)+...]/ −ıπv/2−1/(2A 0 )(...) (9) where . . .denote further terms, and To explore a simplified semi-classical process, illustrating some of the principles of the above approach, we developed a Matlab program superimposing successive reflections of a 360 degree fan of 2 12 − 2 18 rays from the centre of the stadium with an oscillating propagator having period a fraction of the linear dimensions of the stadium.Cells are summed by superimposing the amplitude of each propagator at the point its trajectory crosses a given cell, weighted by the stability factor m 11 as above.
In theory this should have a similar effect to the selection of the orbit with the fractional momentum arriving at a given cell at a given time, since, as we are adding the cumulative effects over successive times, orbits with an integral number of wavelengths over a given path will interfere constructively and those not maintaining phase will effectively cancel one another out.
A more refined method would be to consider the heteroclinic orbits passing between two Gaussian distributions around the initial and final states, by selecting from the fractal intersections of the forward orbits from the initial state and the backward orbits from the final state those with a fractional momentum/period.One can then integrate this process over all the final states to give a time evolving distribution.
In panels (h, i) of Figure 9 the semi-classical method for a propagator having 1/60th of the width of a 600 pixel stadium is compared with the cellular automaton wave function method with a circular wave packet of wavelength a 1/54th fraction of the 594 cell width of the stadium.Although the profiles of the two simplified methods and their scarring patterns are not identical, both display similar features of scarring by wave coherence, illustrating how the two methods can highlight the quantum suppression of chaos.The semi-classical version also shows prominent effects of caustics caused by the curved ends as verified in panels (j-o) of Figure 9.
As shown in panels (j-o) of Figure 9, with the stability factor m 11 inverted to highlight unstable orbit sections, because it is developed from a single fan of rays, the corresponding plots for successive reflection paths using superimposed absolute (or squared) amplitudes in panels (j-n) do not manifest closed orbit scarring, but rather display periodic versions of the caustics found when classical rays are superimposed in panel (o).We have thus made a complete transition from classical chaos to caustic order.

Matlab Files
Source Matlab files for the project can be downloaded as Supplement 6.

Summary
The central method introduced in the paper, provides a straightforward cellular automata model for a semiclassical approach to quantum chaos, more natural than the semi-classical methods of [5][6][7][8], demonstrating quantum scarring of the wave function around repelling classical orbits in fine detail.It thus both replicates the quantum and semi-classical approaches of Tomsovic and Heller and provides a transparent intuitively clear way of imaging and envisaging how quantum chaos actually comes about.The classical orbit simulations further test some of the assumptions of the semi-classical theory using rays, confirming some of the statistical results while providing a significantly different profile from other previous results.

Appendix: The Quantum Chaos Connection with the Riemann Zeta Function
A variety of lines of evidence [12] link phenomena of physics, and in particular quantum chaos, to the Riemann zeta function expressible both as a harmonic Dirichlet sum and an Euler product over primes [13].This has in turn led to hopes that quantum chaos might help solve the Riemann hypothesis that the non-trivial zeros of zeta, those not lying regularly on the negative real axis, all lie on the critical line x = 1 2 .A major breakthrough was thought to have been made when it was discovered that pair correlations in the gaps between the zeta zeros followed the same GUE statistics as chaotic quantum systems and energy levels of large nuclei The GUE statistic, and its time-reversible real variant, the Gaussian orthogonal ensemble (GOE), appear in many forms of quantum systems whose classical analogue is chaotic.These include the many body problem of nuclear energetics, highly excited atoms in a magnetic field and the quantum stadium problem.One of the defining moments in this interaction of fields was Keating's development [14], along with Berry [15] of a formula for the zeta moments, using characteristic polynomials of unitary matrices, based on the operator H BK = −ı x d dx + 1 2 , which however has a continuous spectrum and cannot yield the hypothetical operator with eigenvalues the nontrivial zeros of the Riemann zeta function [16].Initially this correspondence between fields caused great excitation in both the mathematics and physics communities, and a number of eminent researchers tried to prove Riemann hypothesis by discovering a system of random Hermitian matrices whose eigenvalues would be real and might correspond to the zeros of the related function thus showing they had to be real and hence those of ζ(s) would be on x = 1 2 .However this program has so far not borne fruit.
In attempting to create a convergence between Hermitian operators and the zeta function, researchers have constructed a variety of candidates.Berry [17] has presented one of the most straightforward of these, the semiclassical 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 would make the dynamics quantally bound.Secondly, there is no elucidated relationship between the primes and the periodic orbits of the Riemann dynamics.In an ingenious application of the approach, Berry [18,19], on the basis the zeta zeros correspond to eigenvalues of a hypothetical Hamiltonian whose classical trajectories are chaotic and without time-reversal symmetry successfully predicted the longer range deviations of the zeta zeros from the GUE, which received confirmation in Odlyzko's high accuracy super-computer calculation of higher zeta zeros [20].
However, an actual candidate Hamiltonian has remained elusive and despite several further candidates such as H = x + 1 x p + 1 p , Berry and Keating as of 2011, although showing this candidate has the same asymptotic mean spectral density as the Riemann zeros, state: We are not claiming that our hamiltonian has an immediate connection with the Riemann zeta function.This is ruled out not only by the fact that the mean eigenvalue density dif-fers from the density of Riemann zeros after the first terms, but by a more fundamental difference in the periodic orbits.For [the above formula] there is a single primitive periodic orbit for each energy E; and for the conjectured dynamics underlying the zeta function, there is a family of primitive orbits for each 'energy' t, labelled by primes p, with periods log p.This absence of connection with the primes is shared by all variants of xp, and our analysis gives no hint of a resolution of this difficulty.[21, p.13]In a more experimental vein Berry has in 2012 developed far field radiative profiles for two zeta variants [22].
As an indication of the nature of the problem, Schumayer and colleagues have constructed potentials with energy eigenvalues equal to the prime numbers and to the zeros of the zeta function and show that these turn out to be multifractals [23].
Connes has constructed a Hermitian operator whose eigenvalues are the non-trivial Riemann zeros [24].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 where 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 [25].The product formula is not over primes, but over all primitive periodic orbits (ppo) for the motion of the surface considered 1 − e −l p (s+m) (14) where l p 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, . . .and a set of non-trivial zeros putatively on the critical line x = 1 2 .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 T p , 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 com-bining 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 Riemann hypothesis, 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 panel (j) of Figure 7. Within this classical system is a dense set of repelling periodicities, any arbitrarily small deviation from which results in a dense orbit, or transition to 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 [26].
In systems like the quantum stadium, the closed orbits 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 generic, but not homologous.
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.The search for quantum systems reflecting the primes and/or zeta zeros may thus suffer from the same problem facing researchers into L-functions that generalize zeta, associated with structures such as the elliptic curves and modular forms pivotal in solving Fermat's last theorem, the primes of other number fields, differential Maass and automorphic forms, referred to by Brian Conrey [27]

in his AMS review:
There is a growing body of evidence that there is a conspiracy among L-functions -a conspiracy that is preventing us from solving Riemann hypothesis! [27, p.351]The program of seeking to prove the Riemann hypothesis using common properties of L-functions seems to suffer from the same regress we have seen above [28].Really the disparate L-functions, rather than shedding new light on the problem, are in a deep sense more elaborate analytic complexifications of the root prime-zeta relationship.The L-functions of elliptic curves and other number fields, for example, are more complicated representations, whose product formulae still ultimately depend on the prime distribution of the integers.Gaussian primes for example, although they come in three types, are ultimately definable in terms of the integer primes.Likewise elliptic curve L-functions are defined in terms of quadratic prime products.Even the varied Maass forms, modular differential functions satisfying the hyperbolic Laplace wave (17) So ultimately we come back to the long-proven fact that the non-trivial zeros lie on the critical line x = 1 2 if and only if the primes are distributed with fluctuations limited by order x 1 2 , since from Riemann himself [13] and Ingham [29], the supremum of real parts of the zeros is the infimum of numbers β such that the error in the prime number theorem π(x) ∼ x 2 dt ln t is O(x β ).Given the fact that the zeros also have to be symmetric around the critical line due to the analytic extension the fixing of the zeros to the critical line simply reflects the fact that the primes are as evenly distributed as they can be given the fact that they cannot.Thus the tantalizing and elegant simply of Riemann hypothesis belies the fact that it is the real integer primes' root asymptotic behavior determining the entire picture.The orbits of chaotic quantum systems provide an intriguing parallel to the primes as 'orthogonal generators',  but to form a full analogue of Riemann hypothesis we need an equivalent of the prime sieving that generates the Euler product relation defining zeta as a complex harmonic series, because the product formula does not hold in the critical strip and the analytic continuation is possible only for the harmonic sum form.This underlying property at the root of what it is to be prime remains distinct from the notion of a quantum candidate having real eigenvalues, and may explain the difficulty of finding a suitable Hamiltonian, which would, in effect, have to sieve the dynamics in terms of the eigenvalue orbits.It thus remains unclear whether quantum chaotic systems can generate a fully-fledged L-function for which a generalized Riemann hypothesis would be applicable.This brings us full circle.Quantum theory poses a unique challenge to our notions of reality.Far from the quasi-theory of knowledge the Copenhagen interpretation asserts it to be, aspects of quantum mechanics, from bit to it, to fuzzy quantum logic and the concept of entanglement appear to be more fundamental even than our classical notions of points and sets.However primes and integers also go as close to the foundations of existence Quanta | DOI: 10.12743/quanta.v3i1.23 as we can envisage.From the nothingness of 0 to the unity of 1 we derive addition and the integers by simple concatenation 1 + 1 = 2 and multiplication and hence the primes, and the prime distribution determining Riemann hypothesis, as a recursive form of concatenation 2 × 3 = 2 + 2 + 2. Whether quantum theory is the foundation of numbers or vice versa remains an enigma to be resolved.

Figure 1 :
Figure1: Quantum chaos: (a) the stadium is densely filled with repelling periodic orbits[9], three of which are shown in black in (d).The quantum solution of the stadium potential well (b)[4] and (d)[5] 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.The quantum solution is scarred on precisely these orbits (d).This causes it 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.Over time, although the behavior may be transiently chaotic, the quantum system eventually settles into a periodic solution or a linear combination of these.Experimental realizations, such as a magnetized electron in a quantum dot (e)[10] [Gaussian orthogonal ensemble (GOE) and Gaussian unitary ensemble (GUE) statistics], or the scanning tunneling view of an electron on a copper sheet bounded by a stadium of carefully-placed iron atoms (f)[11], confirm the general picture, although quantum tunneling in the latter case leaked the wave function outside too much to demonstrate proper scarring.The semi-classical approach matches closely the quantum calculation (g).

Figure 2 :
Figure 2: Evidence of scarring: Superposition of the first 3300 iterations of a wave function initiated by two different generating functions (a and c in Figure 3) with wavelength chosen to be a fraction of the associated classical repelling orbit.Panels (a) and (b) show superposition of phase amplitudes, whereas (c) and (d) the corresponding superimposed probability distributions of the squared amplitudes, showing evident scarring.Generally the probability distributions give better resolution of the scarring.(e) The wave function cellular automaton iterates in two stages.First in the transition 1 → 2 white positions move using black neighbour values.Then in the transition 2 → 3 black positions move.Finally in the transition 3 → 4 white positions move again.This enables a discrete simulation of a Hamiltonian process.Supplement 1: Quicktime movie of the wave function cellular automaton operating.

Figure 3 :
Figure 3: Scars of the wave stadium function illustrated corresponding to 3,300 iterations of 6 distinct repelling periodic orbits in the classical system.Each of the wave functions in panels (a-h) is iterated from the corresponding generating function in panels (a*-h*), whose wavelength and direction is chosen to correspond to a fraction of the orbit length.The scars are formed by superimposing the squared wave functions at each complete revolution of the phase angle (corresponding to 11 steps of iteration because of the choice of stadium size and initial generator).The rectangular orbit scar in panel (g) also shows resonance with the bow tie in panel (b).Its orbit-length-dividing period of 1.2 relative to panel (a) is close to the bow tie's 1.29 and it closely shares the vertical section of its orbit with panel (b).(h) A vertical wave generator corresponding to the neutral vertical orbit, produces an ordered wave function.(a*-h*) The generating functions for the scars in panels (a-h) consisted of radial Gaussian wave packets of plane cosine waves aligned along the classical repelling orbits shown with wave periods a fixed fraction of the orbit length.

Figure 4 :
Figure 4: The long-term stability of the scars in the cellular automaton simulation is confirmed when the cases (b) and (f) in Figure 3 are run for a further 3300 iterations and the initial transient states of the generator are thus omitted entirely.(a, b) Iterations 1-3300 superimposed.(c, d) Iterations 3300-6600 superimposed.In these figures the absolute wave amplitudes have been added rather than the squared amplitudes to illustrate various methods of scar highlighting.

Figure 5 :
Figure 5: (a, b) The corresponding non-chaotic wave functions for generators a*, b* in Figure 3 in a rectangular potential well do not show scarring, but form regular ordered solutions.(c, d) Non-chaotic wave function corresponding to a quasi-periodic orbit in a lemon-shaped potential well with the quasi-periodic orbit superimposed.(e) The cellular automaton also well models the semi-classical view of the spreading of a traveling wave packet.A small traveling wavelet, generated by clipping sector of a circular ripple with a Gaussian envelope, bounces back and forth, ultimately forming a periodic wave pattern, due to internal resonances generated by its characteristic frequency and wave spreading.Wave spreading compensates for the exponential spreading of the classical repelling orbit and the wave function can form standing wave constructive interference when its energy and consequently its frequency and period corresponds to a fraction of the length of the periodic orbit.Supplement 2: Avi movie of the spreading of the traveling wave packet.

Figure 7 :
Figure 7: (a) A modification of the rectangular potential well consists of a discrete model of a double slit interference experiment, showing the cellular automaton well models circular wave fronts emerging from a double slit and forming constructive and destructive interference fringes.(b) A vertical wave generator produces an ordered wave function on the stadium.(c-e) Probability distributions for a circular wave packet of period a fraction of the overall dimensions for 3300 iterations for the stadium (c), lemon (d) and rectangle (e) display the differences between chaotic and ordered wave forms.(f) Autocorrelation function of the generated bow tie scar in panel (b) of Figure 3. (g) Fast Fourier transform giving frequency domain.(h-k)Classical orbits in the stadium display the classic features of chaotic billiards.(h) Sensitive dependence on initial conditions illustrated for two closely related initial conditions (red and green) show exponentiating divergence with each reflection.Chaos also causes the space to be densely permeated with repelling periodic orbits, two of which are illustrated in (i) and (k).Because they are repelling, neighboring orbits are thrown further away, rather than being attracted into a stable periodic orbit as in an ordered system.Consequently almost all orbits are recurrent and densely fill phase space (j).

Figure 8 :
Figure 8: A collection of classical 16-step close recurrences approximating periodic repelling orbits, some with reverse orbit symmetrization to 32 steps.Repelling orbits symmetric by horizontal reflection (a, b), symmetric by a rotation by 180 degrees (c, d), symmetric by vertical reflection (e, f),and asymmetric (g, h).Supplements 3-5: Quicktime movies of symmetric period 16 orbits passing through (0,0), (0,1) or (2,0).(i) Distribution of lengths of 6444 close orbit recurrences of up to 20 reflections generated from positions (0,0), (2,0) and (0,1), having rotational and reflective symmetry shaded to give relative densities.(j) Corresponding distribution of orbit length differences of the ranked length array.The right hand peak is an artifact caused by missing orbits with higher reflection numbers.(k) Distribution of length differences by frequency shows a descending law in which larger gaps are rarer.This pattern is repeated in the inset when the outliers are removed by choosing only the values from 2000 to 5000.This distribution differs markedly from the distribution of primes (see panel (e) of Figure11) and is unlikely to be an artifact of the orbit selection because it is clear there are sequences of orbits with orbit lengths converging in the limit, such as those exemplified in (a), which tend asymptotically to the circumference of the stadium.According to the trace formula given by Equation 1 the phase shifts and sensitivity of the orbits also have to be included before a GUE distribution is arrived at, however, as shown in (l) a roughly GUE distribution, with prominent outliers corresponding to the shortest groups of orbits, does arise when the inverses of the orbit lengths are taken as a measure of energy.

Figure 9 :
Figure9: (a) Extended orbit of 200 reflections in both directions from a close homoclinic recurrence generating in (b) the repelling bow tie orbit representation (black) later shows heteroclinic recurrence to other orbits, a related period 4 skewed bow tie (blue) and a period 8 orbit (red).Because the dynamics is conservative and area in phase space is conserved, the chaotic process consists of horseshoe saddles and repulsion in one eigenspace is compensated by attraction on the other.(c) Single directed orbit of 500 reflections from (0, 1) with direction (1, −3) can be transformed by the coordinate transformation shown into a momentum-position representation in which P is the cosine of the angle of the ray to the tangent and Q is the perimeter distance e.g. from the mid-point of the right circular arc.This 2D coordinate system results in the discrete dynamic below for 50000 iterations, in which phase space is ergodically (pseudo-randomly) filled except for a diminishing pair of central regions corresponding to the neighborhood of vertical stable orbits.(d) The one-step partition of phase space induced by colouring all the orbits fanning off (±1, ±1), colour-coded as in the discussion.(e) the corresponding partition for four steps.(f) The region mapped from the red ellipse in two reflection steps becomes fractally recurrent after six steps (g).(h) Semiclassical superimposed periodic ray trace of the stadium from a fan of rays at the centre with propagator having period a fraction of the stadium dimensions for 4 reflections.(i) The wave function cellular automaton, with a circular wave packet generator, having a similar fractional period, relative to the stadium width, to the semiclassical example.(j-o) Semiclassical superimposed periodic ray trace of the stadium from a circular wave generator at the centre with period a fraction of the stadium dimensions for (j) 2, (k) 3, (l) 4, (m) 5, and (n) 8 reflections.Rather than highlighting the scars around existing periodic orbits, the method highlights periodic versions of the successive caustics caused by repeated reflection also evident in (o), the classical ray trace for 4 reflections, shown for comparison.

Figure 10 :
Figure 10: (a)  The spacing between neighbouring zeta zeros 10 2 − 10 4 is compared with a GUE distribution 32π −2 s 2 e −4s 2 /π (red), showing coincidence of the two statistics, and a GOE 1 2 πse −πs 2 (green) characteristic of the Wigner distribution of atomic nuclear energies.The relationship is short range only and deviates for larger spacings[20] such as 100 apart (c).(b) Pair correlations for the first 10 5 neighbouring zeros compared with the theoretical GUE distribution.(c) Berry's predictions using semi-classical analysis (dashed line) of longer-range deviations of the zeta zeros from the GUE were confirmed numerically (black circles) using Odlyzko's numerical calculations.

Figure 11 :
Figure 11: Correlation coefficients of primes (a, c) and zeta zeros (b, d).While prime pairs have little correlation in their gap sizes, zeta zeros show consistent correlation, varying both as one moves up the zeros and as one examines longer range effects.(a, b) Graph of correlation coefficients using the first 20000 primes and zeros correlating differences of adjacent values with those j units further on.(c, d) Colour map where (i, j) is the correlation coefficient of the 500 pairs (X, Y) where X(k) = (v(k), v(k + 1)) and Y(k) = (v(k + j), v(k + j + 1)), k = i, . . ., i + 500.(e, f) Distributions of the first 600 primes (e) and zeta zeros (f).
function ∆ = −y 2 ∂ 2 ∂x 2 + ∂ 2 ∂y 2 ,have Euler products which are again generated in terms of the integer primes p, as illustrated by the cubic Mass form L(z, ϕ × χ) =

Figure 12 :
Figure 12: (a) Pattern of differences between successive zeros of zeta in the range from 10 21 for 10 4 successive zeros determined by Andrew 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.(b) The semi-classical potential (dashed line), and the fractal potential (solid line) supporting the first two hundred zeros of ζ(s) as energy eigenvalues [23].Inset is the difference.

Figure 13 :
Figure 13: A menagerie of L-functions: Riemann zeta, a Dirichlet L-function modulo 5, Dedekind zeta and a Hecke L-function of the complex primes of the Gaussian integers as an extension field, the L-function of the elliptic curve y 2 + xy = x 3 + x 2 −210x−441, third and fourth degree Maass forms and the modular discriminant all appear to have their non-trivial zeros on their weighted critical line, but so far have not shed new light on Riemann hypothesis [28].