Exploring Quantum, Classical and Semiclassical Chaos in the Stadium Billiard
Chris King
Mathematics Department – University of Auckland – 9 Jul 2009
PDF 2013 Quanta 3/1 16-31 DOI:10.12743/quanta.v3i1.23

 

You can now run the wave function CAs in my open source application CA2D for Mac!

 

Abstract: 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 [i]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  inside a box with arbitrarily shaped walls is

 [1]

where  is a smooth function giving the mean density of states, summed over all orbits  of length  with  the phase shift counting the focal points and twice the number of reflections off the walls and over all retracings k of these orbits. The stability matrix records the sensitivity 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 [iv] and this function’s GUE type statistics shared by quantum chaotic systems (fig 2e) has led to extensive attempts to solve quantum chaotic systems through the paradigm of the zeta zeros.

 

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 CA below. Alternatively the semi-classical 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 [v] [vi] among others have used the semi-classical orbit-based approach to develop the theory of scarring of the quantum wave [vii]. 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. A particularly complete expose of the methods involved is provided in their Physical Review paper [viii].

 

Fig 2: Quantum chaos: (a) [ix], the stadium is densely filled with repelling periodic orbits, three of which are shown in black in (d). The quantum solution of the stadium potential well (b) [x] and (d) [xi] 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 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) [xii] ((a) and (b) give GOE and GUE statistics respectively), or the scanning tunneling view of an electron on a copper sheet bounded by a stadium of carefully-placed iron atoms (f) [xiii], 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).

 

The Quantum Cellular Automaton

To get a more direct picture of the process of wave spreading and wave function scarring, 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 grid 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 of fig 2(f).

 

Fig 3: The wave function CA iterates in two stages. First in a > b white positions move using black neighbour values. Then in b>c black positions move. Finally in c>d white moves again. This enables a discrete simulation of a Hamiltonian process [xiv].

 

The cellular automation iterates the wave function according to the rule:

 [2]

As in figure 3, 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 .

 

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 illustrated by the following two edge and corner situations:

  [3]

 

The process is also capable of handling applied forces by integrating their effect:

 

   [4]

 

The simulation uses a variety of generating function 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:

 [5]

 

combining one or more ‘plane cosine waves of appropriate period lying along an orbit path clipped being multiplied by a circular Gaussian, although circular and uni-directional 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 fig 4. The predominant repelling orbits clearly show as probability peaks coinciding with the corresponding repelling orbits in just the matter of the semi-classical simulations already noted in the references and discussion. In fig 5 are shown the corresponding generating functions for each of the images in fig 4.

 

Fig 4: Scars of the wave stadium function illustrated corresponding to 3,300 iterations of 6 distinct repelling periodic orbits in the classical system.  The wave function is iterated from the corresponding generating function in fig 6, 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 (g) scar also shows resonance with the bow tie (b). Its orbit-length-dividing period of 1.2 relative to (a) is close to the bow tie's 1.29 and it closely shares the vertical section of its orbit with (b).(Bottom) (h) A vertical wave generator corresponding to the neutral vertical orbit, produces an ordered wave function.

 

The images in fig 4 were generated using the first 3300 iterations of the generating function, however when the scars (b) and (f) in fig 4 were iterated on from 3300 to 6600 iterations the scar remained stable , as shown in fig 8, verifying it is not just a transient feature of the expanding generating functions.

 

Fig 5: The generating functions for the scars in fig 4 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.

 

Fig 6: (Left):The corresponding non-chaotic wave functions for generators a, b in fig 4 in a rectangular potential well do not show scarring, but form regular ordered solutions.  (Right): Non-chaotic wave function corresponding to a quasi-periodic orbit in a lemon-shaped potential well with the quasi-periodic orbit superimposed below.

 

These patterns of scarring can also be compared with other potentials 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 fig 6 the generating functions of scars (a) and (b) in fig 4 are applied again for 3300 iterations in the rectangular well resulting in ordered wave functions, and a vertical generating function in the lemon well produces an ordered wave function corresponding to the quasi-periodic stable orbit superimposed below.

 

Fig 7: 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.

Significantly the CA is able to well handle oblique generating functions with periods not reducible directly to an integer number of iterations.  In fig 1 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 CA is able to successfully superimpose amplitudes in phase over thousands of iterations.

 

Fig 8: The long-term stability of the scars in the CA simulation is confirmed when the cases (b) and (f) in fig 4 are run for a further 3,300 iterations and the initial transient states of the generator are thus omitted entirely. Top: Iterations 1-3300 superimposed. Bottom: 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.

 

The CA can also be used to model the simulations of traveling wave packets illustrated from other research. In fig 7 is shown a series of images of iterations of a right-moving wave packet giving a similar evolution profile to the one generated by Heller and Tomsovic in fig 2(c).

 

Fig 8b: 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 fig 4, 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 ititial excitations.

 

The CA can also well-model other wave experiments in physics such as the double-slit experiment illustrated in fig 9.

 

Fig 9: A modification of the rectangular potential well consists of a discrete model of a double slit interference experiment, showing the CA well models circular wave fronts emerging from a double slit and forming constructive and destructive interference fringes.

 

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 fig 10 the three boundaries are compared and the scarring characteristic of the stadium shows up clearly and quite elegantly suggestive of several superimposed orbit solutions.

 

 

Fig 10: Probability distributions for a circular wave packet of period a fraction of the overall dimensions for 3300 iterations for the stadium, rectangle and lemon display the differences between chaotic and ordered wave forms.

 

The CA can also be easily interrogated to produce an autocorrelation function  [6] by multiplying the states of the generating wave function  by the propagated wave function  and integrating over all cells in the array. In fig 11 is 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.

 

Fig 11: (Left) Autocorrelation function of the generated ‘bow tie’ scar in fig 4(b). (Right) FFT giving frequency domain.

 

Orbits in the Classical Stadium

To complement the wave CA and draw superimposed classical orbits, a ray drawing algorithm was developed for the same matrix of elements used in the wave CA algorithm. This also enables portrayal of several features of the associated classical chaotic system.

 

Fig 12: Classical orbits in the stadium display the classic features of chaotic billiards. Top left: 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 below.  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 (top right).

 

In fig 12 are shown how sensitive dependence on initial conditions (in this case a mall angle difference in the initial orbit leads on multiple reflection in the curved sections of the stadium to exponentially incr3asing angular divergences between the orbits, eventually 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 give trajectory, as shown in fig 15

 

Fig 13: A collection of classical 16-step close recurrences approximating periodic repelling orbits, some with reverse orbit symmetrization to 32 steps. Top row symmetric by horizontal reflection, second row by a rotation by 180^o, third by vertical reflection and the fourth row asymmetric.

 

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 approximation, they are really detecting close recurrence in a finite number of reflections, rather than an exact periodic solution.

 

Thus is a dramatic increase is made in the number of reflections, such an orbit will eventually break stability as the periodic solution is repelling, and may over time recur closely to other repelling orbits. In fig 15, the computed solution to the ‘bow tie’ 4-period in fig 4(b) 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 fig 13 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.

 

The algorithm was set to compute all close recurrences 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 [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 Gaussian unitary ensemble (GUE).

 

Fig 14: (Above) 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. (Below) 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.

 

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 fig 14, 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.

Nevertheless when the inverses of orbit lengths are taken to give a measure of energy as inversely proportional to wavelength, a GUE-like statistical distribution does arise, suggesting the GUE statistics seen in experiments on the quantum stadium do result from a collection of eigenvalues each corresponding to a closed orbit.

 

Fig 15: (left) 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 f primes 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 the one’s exemplified top left in fig 11, which tend asymptotically to the circumference of the stadium. According to the trace formula [1] the phase shifts and sensitivity of the orbits also have to be included before a GUE distribution is arrived at, however, as shown (right) 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.

 

The perimeter-angle coordinate system described in fig 16 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.

 

Fig 16: (Top left) Extended orbit of 200 reflections in both directions from a close homoclinic recurrence generating the repelling bow tie orbit representation (bottom left 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. (Top right) 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.

 

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 between curved and straight sections at (±1, ±1), labelled in red, magenta, green and blue dots in fig 16.

 

Fig 17: (Top left) The one-step partition of phase space induced by colouring all the orbits fanning off (±1, ±1), colour-coded as in the discussion and fig 16. (Top right) the corresponding partition for four steps. (Bottom left) The region mapped from the red ellipse in two reflection steps becomes fractally recurrent after six steps (bottom right).

 

The Semi-Classical Approach

The semi-classical approach retrieves very comparable results to a full quantum representation, fig 2(g), 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 and applying an oscillating wave propagator to the member of the form ^[viii]:

 [7]

where S is the action, with j weighted according to the product of the orthogonal stability matrices:

 [8]

along the orbit, giving a correlation function for a single orbital family member as follows:

[9]

 

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^o fan of 2¹² – 2¹⁸ 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  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 fig 18 the semi-classical method for a propagator having wavelength 1/60^th of the width of a 600 pixel stadium is compared with the CA 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 fig 19.

 

Fig 18: (Top) 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. (Bottom) The wave function CA, with a circular wave packet generator, having a similar fractional period, relative to the stadium width, to the semiclassical example above.

 

As shown in fig 19, because it is developed from a single fan of rays, the corresponding plots for successive reflection paths using superimposed absolute (or squared) amplitudes (a-e) do not manifest closed orbit scarring, but rather display periodic versions of the caustics found when classical rays are superimposed (f). We have thus made a complete transition from classical chaos to caustic order.

 

Fig 19: 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 (a) 2, (b) 3, (c) 4, (d) 5, and (e) 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 (f), the classical ray trace for 4 reflections, shown for comparison.

 

Matlab Files

 

Source matlab files for the project can be downloaded from:

http://www.dhushara.com/DarkHeart/QStad/QStad.zip

Further matlab files for the zeta function and discrete chaos are available at:

http://www.dhushara.com/DarkHeart/RH2/RH.htm

http://www.dhushara.com/DarkHeart/DarkHeart.htm

 

References