Quantum-mechanical return probability of a wave packet and its relationship to classical periodic orbits

by

Albert Dvornik

Submitted to the Department of Physics in partial fulfillment of the requirements for the degree of

Bachelor of Science in Physics

at the

MASSACHUSETTS INSTITUTE OF TECHNOLOGY

February 1998

\begin{matrix}c \\ ◯ \end{matrix} Albert Dvornik, MCMXCVIII. All rights reserved.

The author hereby grants to MIT permission to reproduce and distribute publicly paper and electronic copies of this thesis document

in whole or in part, and to grant others the right to do so.

Author . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Department of Physics

January 26, 1998

Certified by. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Michel Baranger Professor of Physics Thesis Supervisor

Accepted by . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . June L. Matthews Professor of Physics

Chairman, Undergraduate Thesis Committee Quantum-mechanical return probability of a wave packet

and its relationship to classical periodic orbits

by

Albert Dvornik

Submitted to the Department of Physics

on January 26, 1998, in partial fulfillment of the

requirements for the degree of

Bachelor of Science in Physics

Abstract

Quantum-mechanical return probabilities are suggested as a way to analyze the recurrent behavior of quantum-mechanical dynamical systems. An analogy is drawn between the maxima of return probabilities for a system and the periodic orbits of the analogous classical system.

A numerical model for modeling the time evolution of a quantum-mechanical state, by iterated application of short-term time-evolution operators to a discretely sampled wave function, is introduced and evaluated. The requirements of such a computational model are analyzed.

This computational model is then applied to the study of return probabilities for a small problem targeted to run on workstation-class computers. Return probabilities are calculated for Gaussian wave packets whose coordinates and momenta correspond to the classical energy E = 0.3. The resulting data demonstrates a high degree of correspondence between the maxima of return probabilities and the corresponding classical periodic orbits, but also demonstrates some of the ways in which discrepancies may arise.

Thesis Supervisor: Michel Baranger

Title: Professor of Physics

Quantum-mechanical return probability of a wave packet

and its relationship to classical periodic orbits

by

Albert Dvornik

Submitted to the Department of Physics

on January 26, 1998, in partial fulfillment of the

requirements for the degree of

Bachelor of Science in Physics

Abstract

Quantum-mechanical return probabilities are suggested as a way to analyze the recurrent behavior of quantum-mechanical dynamical systems. An analogy is drawn between the maxima of return probabilities for a system and the periodic orbits of the analogous classical system.

A numerical model for modeling the time evolution of a quantum-mechanical state, by iterated application of short-term time-evolution operators to a discretely sampled wave function, is introduced and evaluated. The requirements of such a computational model are analyzed.

This computational model is then applied to the study of return probabilities for a small problem targeted to run on workstation-class computers. Return probabilities are calculated for Gaussian wave packets whose coordinates and momenta correspond to the classical energy E = 0.3. The resulting data demonstrates a high degree of correspondence between the maxima of return probabilities and the corresponding classical periodic orbits, but also demonstrates some of the ways in which discrepancies may arise.

Thesis Supervisor: Michel Baranger

Title: Professor of Physics

Acknowledgments

I am deeply indebted to a large number of people who helped, in one way or another, with the completion of this work. I hope that I have not omitted any of them, and that they will not be offended if I have.

My thesis advisor, Michel Baranger, was both very helpful in the development of the material, and far more understanding than I had any right to deserve. Robert Cruikshank and Abbe Cohen read an advanced draft and helped improve the clarity of my writing. The calculations were speeded up by the FFTW software package, written by Steven Johnson and Matteo Frigo, which was brought to my attention by Dr. Power.

Much gratitude goes to my friends and relatives, who have all offered much emotional support, at the times when I needed it the most. My parents, in particular, supported me against all odds, which is yet another thing that I won’t ever be able to repay them for. I must also thank the past and present staff of the Undergraduate Physics Office for the large amount of help I have received there over the years.

Finally, Craig Fields came up with the flipbook idea. Since I am the one writing this paragraph, I disavow any credit for it, and reject all responsibility as well.

Contents

1 Introduction 8 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2 A few notes on notation . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Modelling the dynamics of a quantum-mechanical system 11 2.1 The time evolution operator . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Means of computing the time evolution . . . . . . . . . . . . . . . . . 12 2.3 Discretization of space . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 Return probabilities 17 3.1 Periodic orbits and return probabilities . . . . . . . . . . . . . . . . . 17 3.2 Gaussian wave packets . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3 Momentum inversion and time reversal . . . . . . . . . . . . . . . . . 19

4 Limitations of a discrete model 21 4.1 A formal problem: aliasing . . . . . . . . . . . . . . . . . . . . . . . . 21 4.2 A physical problem: good sampling . . . . . . . . . . . . . . . . . . . 23

5 Decisions about computation 26 5.1 What to calculate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 5.2 Choosing the potential . . . . . . . . . . . . . . . . . . . . . . . . . . 27 5.3 Choosing parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.4 A sample evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 6 Analysis of results 40

6.1 Return probability plots . . . . . . . . . . . . . . . . . . . . . . . . . 40 6.2 Return probabilities of the Delfi potential . . . . . . . . . . . . . . . . 41 6.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

A Return probability plots for the Delfi potential with E=0.3 54

B Software used for wave-function propagation 124 B.1 Choosing the programming language . . . . . . . . . . . . . . . . . . 124 B.2 Complex number conventions . . . . . . . . . . . . . . . . . . . . . . 125

Bibliography 127

List of Figures

5-1 Contour plot of the Nelson potential . . . . . . . . . . . . . . . . . . 28 5-2 Contour plot of the Delfi potential . . . . . . . . . . . . . . . . . . . . 28 5-3 The E–τ plot for the Delfi potential, reproduced from [Kagarlis, 1986] 30 5-4 The (x = 0, p_{x} > 0) Poincaré section of the Delfi potential for E = 0.3 32 5-5 Locations of periodic orbits in fig. 5-4, identified by family . . . . . . 32 5-6 The (y = 0, p_{y} > 0) Poincaré section of the Delfi potential for E = 0.3 33 5-7 Locations of periodic orbits in fig. 5-6, identified by family . . . . . . 33 5-8 Time-evolution from a sample initial state . . . . . . . . . . . . . . . 37 5-9 Time-evolution from a sample initial state (continued) . . . . . . . . 38 5-10 Time-evolution from a sample initial state (continued) . . . . . . . . 39

6-1 Return probabilities for t = 3.8 . . . . . . . . . . . . . . . . . . . . . 41 6-2 Return probabilities for t = 7.4 . . . . . . . . . . . . . . . . . . . . . 43 6-3 A plot of the classical Family 2 trajectory at E = 0.3 superimposed on

the return probabilities for t = 7.4 . . . . . . . . . . . . . . . . . . . . 43 6-4 A plot of the classical Family 9 trajectory at E = 0.3 superimposed on

the return probabilities for t = 7.4 (along with Family 2). . . . . . . . 44 6-5 Time-evolution starting from the “narrow” bifurcation of Family 9 . . 45 6-6 Time-evolution starting from the “wide” bifurcation of Family 9 … 46 6-7 Return probabilities for t = 7.7 . . . . . . . . . . . . . . . . . . . . . 48 6-8 Return probabilities for t = 8.0 . . . . . . . . . . . . . . . . . . . . . 49 6-9 A plot of the classical Family 4 trajectory at E = 0.3 superimposed on

the return probabilities for t = 8.0 . . . . . . . . . . . . . . . . . . . . 49 6-10 Return probabilities for t = 11.6 . . . . . . . . . . . . . . . . . . . . . 50 6-11 Return probabilities for t = 7.0 . . . . . . . . . . . . . . . . . . . . . 52 6-12 A plot of the classical Family 4 trajectory at E = 0.3943 superimposed

on the return probabilities for t = 7.0 . . . . . . . . . . . . . . . . . . 52

Chapter 1

Introduction

1.1 Motivation

One of the particularly useful ways of studying classical chaotic systems is the examination of their periodic orbits. It seems sensible to assume that an examination of a similar concept for a quantum-mechanical chaotic system might yield equally illuminating results. The work that follows examines a possible way of achieving this, and applies it to a particular example sized to be computationally feasible.

First, I will discuss the time-evolution operator in its usual form, and then demonstrate several approximations well-suited to be used in a numerical model that uses a discretized wave function. This provides the background needed to define the quantum-mechanical return probabilities of a wave packet. I will suggest a possible relationship between those return probabilities and the periodic orbits of an analogous classical system.

I will also discuss some of the limitations of the numerical model I used, and mention ways of avoiding those limitations. This numerical model will then be applied to the study of a specific problem: the behavior of a quantum-mechanical system based on the Delfi potential [Kagarlis, 1986] for the energy E = 0.3. Finally, the results of this computation will be discussed.

1.2 A few notes on notation

I have attempted to express the ideas behind this text as clearly as I could. However, in order to avoid loading the text with repetitive explanations of what the reader may already be familiar with, I adopted some notational conventions that are in more or less common use in the physics community.

Vector quantities are typeset in boldface, so x denotes a vector while x_{i} is a scalar. An integral over a vector quantity is assumed to be an integral over the entire set of values of that quantity; thus, ∫ d^{D}x is just shorthand for ∫_{−∞}^{∞}··· ∫_{−∞}^{∞} dx_{1} … dx_{D}.

Since some of the covered topics invite confusion between state variables or functions and the corresponding operators, the operators are identified by a circumflex (“hat”), for example: Â |α⟩ = A |α⟩.

I will primarily use the Dirac bra-ket notation because it seems to clarify certain questions, such as whether a state is being considered time-dependent or fixed. Unless otherwise noted, all kets under consideration will be viewed in the Schrödinger “picture”: the basis vectors are held time-independent, and the ket corresponding to an evolving state changes with time.

Particular ket symbols will be used for certain commonly used eigenkets. The energy eigenfunctions of a system are written as |n⟩, with Ĥ |n⟩ = E_{n} |n⟩. The ket |x⟩ is an eigenket of position with the eigenvalue x. (More verbosely, |x⟩ is a simultaneous eigenket of all position operators x̂_{i}, and the respective eigenvalues of each of those operators are the components of the vector x.) Likewise, |p⟩ is an eigenket of momentum with the eigenvalue p.

I will occasionally make use of the convention, originated by Einstein, of implied summation over a repeated index variable that occurs on only one side of the equation. This helps with the readability of long expressions with many sums, where writing out the summations explicitly would be cumbersome. However, since I prefer writing out sums where it is possible to do so, the use of this convention could go unnoticed and cause some confusion. To prevent that, I will make a note of its use in the text surrounding the equations.

Finally, matrix mechanics sometimes blurs the distinction between infinite sums and integrals over a basis. For instance, the completeness requirement for a basis is usually written as ∑_{a} |a⟩ ⟨a| = 1, but when applied to the |x⟩ basis, it really means to say ∫ |x⟩ ⟨x| d^{D}x = 1. Those two forms are clearly discrete and continuous formulations of the same idea, but the sum and the integral are in fact not equivalent in the strict mathematical sense. However, the relationship between the two is trivial, and it is usually clear whether summation or integration should be applied for a particular set of kets. Consequently, I will often write a summation over a set of basis vectors where I really mean “summation or integration, as appropriate,” and I hope I will be forgiven for it.

This notation will also be used to extend the Einstein convention to an implied integral over a set of basis vectors when expressions with matching indices (like |x_{i}⟩ ⟨x_{i}|) occur on one side of an equation. As with the regular implied summation, any use of this will be noted in the text.

Modelling the dynamics of a

quantum-mechanical system

2.1 The time evolution operator

The dynamical behavior of the state of an isolated quantum-mechanical system is described by the Schrödinger equation,

iℏ\frac{∂}{∂t} |ψ(t)⟩ = Ĥ |ψ(t)⟩ .  (2.1)

Clearly, fully specifying the state of the system at any given time uniquely determines the state of the system for all other times.

The equation (2.1) combined with the definition of derivative yields an expression relating the ket at time t to the ket for an infinitesimally later time t + dt:

|ψ(t + dt)⟩ = |ψ(t)⟩ − \frac{i}{ℏ}Ĥdt |ψ(t)⟩ = (1 − \frac{i}{ℏ}Ĥdt) |ψ(t)⟩ .  (2.2)

The expression in parentheses is evidently an operator which propagates the ket in time by an infinitesimally small amount dt. If the Hamiltonian is time-independent, then successive infinitesimal propagations can be compounded to form an operator that effects an arbitrary finite time difference:

Û(t) = lim_{N→∞} (1 − \frac{i}{ℏ}Ĥ \frac{t}{N} )^{N} = exp (−\frac{iĤt}{ℏ}) .  (2.3)

The operator Û is commonly referred to as the time-evolution operator^{1} since, when applied to an initial ket, it traces out the dynamical evolution of that state in time:

|ψ(t)⟩ = Û(t − t_{0}) |ψ(t_{0})⟩ .  (2.4)

Since Ĥ is hermitian, Û is unitary, i.e. Û^{†}Û = ÛÛ^{†} = 1. This is important because it means time evolution conserves probability. Another property, which is easy to verify by inspection, is that Û(−t) = Û^{†}(t).

2.2 Means of computing the time evolution

In principle, the time-evolution operator enables us to determine the exact wave function at any point in time from a given set of initial conditions. However, the form expressed in equation (2.3) is not immediately useful for numerical modeling. In order to be stored in a computer, the state of the system must be approximated by a finite set of numeric values. This can be accomplished either by decomposing it into components with respect to some orthonormal basis, or by simply sampling the wave function (which is, of course, a decomposition into the set of simultaneous eigenfunctions of {x̂_{i}}). Since the number of components in the decomposition will be finite, the basis in question will not be complete, but even an incomplete basis should yield acceptable results for some class of “well-behaved” functions.

The eigenfunctions of Ĥ form a basis that is particularly easy to work with. Let |n⟩ be an eigenfunction of the Hamiltonian with the eigenvalue E_{n}. Then |n⟩ will also be an eigenfunction of Û, propagating as

1Sometimes the term propagator is also used, although the word is more commonly used in the literature for the kernel K(x′, t′; x, t) = ⟨x′| Û(t′ − t) |x⟩ of the x representation of the evolution operator written in integral form.

|n; t⟩ = Û(t) |n⟩ = e^{−iE_{n}t/ℏ} |n⟩ .  (2.5)

The eigenfunctions of Ĥ are shown once again to be stationary states, changing only by an overall phase under time evolution. The decomposition of the wave function for any time t in terms of this basis is thus trivially related to the decomposition at t = 0. (Finding the eigenfunctions of Ĥ for a general Hamiltonian is, of course, much harder, but can be accomplished numerically.) This approach has already been used successfully, for example in [Baranger, 1993]. I attempted to determine if satisfactory results could be obtained using a different method.

For a basis other than the eigenfunctions of Ĥ, treating the exponential in equation (2.3) directly as a series sum is generally impractical for numeric calculations, due to the computational costs of applying Ĥ^{n}. For example, for decomposition into |x⟩, applying Ĥ^{n} involves taking 2n-th derivatives, which becomes increasingly hard to do well for a sampled function as n grows large.

Let us assume that the Hamiltonian can be separated into the form Ĥ = K̂ + V̂, where the kinetic energy operator is the usual K̂ = p̂^{2}/2m, while the operator corresponding to potential energy is given by V̂ = V (x̂_{1}, …, x̂_{n}), where V (x_{1}, …, x_{n}) is an arbitrary hermitian (and hence real) function of coordinates alone. This will clearly be true whenever the system corresponds to a single particle in a potential well.

Because of their structure, K̂ is particularly simple to apply to the simultaneous eigenfunctions of {p̂_{i}}, while V̂ is easy on the eigenfunctions of {x̂_{i}}. More specifically, K̂ is diagonalized in the momentum representation, while V̂ is diagonalized in the position representation:

K̂ = ∑_{p} |p⟩ \frac{p^{2}}{2m} ⟨p|  V̂ = ∑_{x} |x⟩ V (x) ⟨x|  (2.6)

Switching from one of these representations to the other is covered in most quantum mechanics textbooks [Sakurai, 1985; Gasiorowicz, 1974], and involves a pair of Fourier transforms,

⟨p|ψ⟩ = ∫ ⟨p|x′⟩ d^{D}x′ ⟨x′|ψ⟩ = (2πℏ)^{−D/2} ∫ e^{−ip·x′/ℏ} ⟨x′|ψ⟩ d^{D}x′,
⟨x|ψ⟩ = ∫ ⟨x|p′⟩ d^{D}p′ ⟨p′|ψ⟩ = (2πℏ)^{−D/2} ∫ e^{+ip′·x/ℏ} ⟨p′|ψ⟩ d^{D}p′.

(2.7)

This observation suggests the idea of separating the evolution operator into “kinetic” and “potential” parts. K̂ and V̂ do not commute, so e^{−i(K̂+V̂)t/ℏ} is not exactly separable into a product of exponentials each containing only V̂ and K̂. However, if t is made small enough compared to ℏ/E, the effect of that non-commutativity decreases with t^{2} and becomes small compared to the first-order terms in the evolution operator. For a small Δt, we can thus write

Û(Δt) ≈ e^{−iV̂ Δt/ℏ}e^{−iK̂Δt/ℏ} ≈ e^{−iK̂Δt/ℏ}e^{−iV̂ Δt/ℏ}.  (2.8)

The two forms above can be combined into a single symmetric form, which is beneficial for two reasons. First, a symmetric form will preserve Û^{†}(t) = Û(−t). Secondly, it effectively halves Δt compared to the original expressions. The result^{2} is

Û(Δt) ≈ Û_{1}(Δt) ≡ e^{−iV̂ Δt/2ℏ}e^{−iK̂Δt/ℏ}e^{−iV̂ Δt/2ℏ}.  (2.9)

Arbitrarily long evolution can be accomplished by iterating application of Û_{1} for acceptably short fragments of the required time interval. Also diagonalizing the parts of Û_{1} analogously to equation (2.6) yields

|ψ(NΔt)⟩ ≈ |x_{i_{1}}⟩e^{−iV (x_{i_{1}})Δt/2ℏ}⟨x_{i_{1}}|p_{j_{1}}⟩e^{−i^{\frac{p{2 \atop j_{1}}}{2m}} Δt/ℏ}⟨p_{j_{1}}|x_{i_{2}}⟩e^{−iV (x_{i_{2}})Δt/ℏ}⟨x_{i_{2}}|p_{j_{2}}⟩ …

… ⟨p_{j_{l−1}}|x_{i_{l}}⟩e^{−iV (x_{i_{l}})Δt/ℏ}⟨x_{i_{l}}|p_{j_{l}}⟩e^{−i^{\frac{p{2 \atop j_{l}}}{2m}} Δt/ℏ}⟨p_{j_{l}}|x_{j_{l+1}}⟩ …

… ⟨x_{j_{N}} |p_{j_{N}} ⟩e^{−i^{\frac{p{2 \atop j_{N}}}{2m}} Δt/ℏ}⟨p_{j_{N}} |x_{i_{N+1}}⟩e^{−iV (x_{i_{N+1}})Δt/2ℏ}⟨x_{i_{N+1}}|ψ(t = 0)⟩ (2.10) where summation (more correctly, integration) over each pair of x_{i} and p_{j} is implied by notational convention. Note that only the two exponentials at the ends have 2 in the denominator; pairs of successive instances of e^{−iV̂ Δt/2ℏ} have been joined everywhere else.

2Of course, we could also symmetrize to e^{−iK̂Δt/2ℏ}e^{−iV̂ Δt/ℏ}e^{−iK̂Δt/2ℏ}, but the other form is more convenient if we want the original and final state to be expressed in position representation.

The occurrences of ⟨p_{j}|x_{i}⟩ and ⟨x_{i}|p_{j}⟩ in equation (2.10), under implied integration, perform a change of basis, and can be identified with the Fourier transforms in equation (2.7). The rest of equation (2.10) consists of multiplications of the wave function by a complex function. Since both of these operations are relatively simple to perform numerically, equation (2.10) should provide a suitable origin for a computational model of wave function propagation.

2.3 Discretization of space

Before we continue, it is convenient to rewrite equations (2.7), making the wave functions into explicit functions of the positions or momenta. It turns out that it’s more useful to write the momentum-space wave function as a function of the wave vector, k = p/ℏ, rather than the momentum:

Ψ(x) = ⟨x|ψ⟩  Φ(k) = ℏ^{D/2} ⟨p = ℏk|ψ⟩  (2.11)

(The factor of ℏ^{D/2} is required for Φ(k) to be properly normalized.^{3}) Using these substitutions, equations (2.7) become

Φ(k) = (2π)^{−D/2} ∫ e^{−ik·x′}Ψ(x′)d^{D}x′,
Ψ(x) = (2π)^{−D/2} ∫ e^{+ik′·x}Φ(k′)d^{D}k′.

(2.12)

We can now consider which representation should be used for the wave function in equation (2.10). The simplest choice is to take the values of the positionrepresentation wave function at a discrete rectangular grid, i.e. to consider the set of values

3If ⟨ψ|ψ⟩ = 1, then ∫ ⟨ψ|p⟩ ⟨p|ψ⟩ d^{D}(p/ℏ) = 1/ℏ^{D}.

Ψ_{l_{1}...l_{D}} = Ψ(x_{1} = δ_{1}l_{1}, …, x_{D} = δ_{D}l_{D}),  (2.13)

when each index l_{i} takes its values from some range of integers. A normalization constant was intentionally omitted from this expression, since it is convenient to ignore it when comparing discrete and continuous representations of the same function. However, it should be noted that this definition sets the normalization of discretized wave functions to (δ_{1} … δ_{D})^{−1/2} instead of 1.

Let N_{i} be the number of distinct values taken by l_{i}. In that case, the momentumrepresentation wave function can be approximated by the discrete Fourier transform,

Φ_{m_{1}...m_{D}} = \frac{1}{\sqrt{N_{1} … N_{D}}} ∑_{l_{1}} ⋯ ∑_{l_{D}} e^{−2πi(\frac{l_{1}m_{1}}{N_{1}} +···+ \frac{l_{D}m_{D}}{N_{D}} )}Ψ_{l_{1}...l_{D}},  (2.14a)
Ψ_{l_{1}...l_{D}} = \frac{1}{\sqrt{N_{1} … N_{D}}} ∑_{m_{1}} ⋯ ∑_{m_{D}} e^{+2πi(\frac{l_{1}m_{1}}{N_{1}} +···+ \frac{l_{D}m_{D}}{N_{D}} )}Φ_{m_{1}...m_{D}}, (2.14b)

where the discrete wave function is related to the continuous one by

Φ_{m_{1}...m_{D}} ≈ Φ(k_{1} = \frac{2π}{N_{1}δ_{1}} m_{1}, …, k_{D} = \frac{2π}{N_{D}δ_{D}} m_{D}).  (2.15)

The number of meaningful distinct values taken by m_{i} will also be N_{i}, since the indices are cyclic (that is, Φ_{...m_{i}...} = Φ_{...(m_{i}±N_{i})...}). Not very surprisingly, Φ_{m_{1}...m_{D}} is also normalized to (δ_{1} … δ_{D})^{−1/2}.

An algorithm known as the Fast Fourier Transform (FFT) provides an efficient way to calculate the expressions (2.14), particularly in the case where all sizes N_{n} are powers of 2. The details of FFT theory are outside the scope of this thesis, but extensive literature exists on the subject [Brigham, 1974; Nussbaumer, 1982]. Very recently, a blazingly fast FFT implementation that adaptively exploits the strong points of different processor architectures was developed by Matteo Frigo and Steven G. Johnson [Frigo and Johnson, 1997]. I started using it after it outperformed my own hand-optimized code, and I highly recommend it.

Return probabilities

3.1 Periodic orbits and return probabilities

One of the most powerful tools in analyzing a chaotic classical conservative system is a study of its periodic trajectories. Such a study, which may involve constructing an E–τ plot and locating the periodic trajectories in the Poincaré sections, is capable of revealing all significant facts about the dynamics of a system.

In the quantum-mechanical case, there are no trajectories as such. Furthermore, the evolution of a localized state cannot be periodic in most potentials.^{1} However, there is still a way to obtain a measure of the periodicity of the evolution of a wave function.

The return probability of a state is defined as

R(t) = |⟨ψ| Û(t) |ψ⟩|^{2} .  (3.1)

The return probability is the probability that the state ψ, after evolving for some amount of time, can again be identified with the original state. If the initial wave function is localized, then a high return probability (for t ≉ 0) indicates that the wave packet has moved around in the potential, but stayed localized and returned to a state similar to the initial one. This suggests that when the ψ is made to resemble the classical particle with some definite position and momentum, the classical periodic orbits ought to correspond to maxima of R. Previous numerical experiments have shown this conjecture to be true [Baranger, 1993].

1Specifically, unless the energy eigenvalues are degenerate or the differences between them are commensurate, any state that is periodic (up to a phase) can be composed of at most two energy eigenstates. Since any state that is localized to a small area of the potential needs to consist of many energy eigenstates, such states cannot be periodic in any system in which the energy steps are incommensurate.

3.2 Gaussian wave packets

Heisenberg uncertainty relations state that a state cannot be localized arbitrarily well in both position and momentum. A state that aspires to resemble a classical particle therefore needs to be both somewhat localized in x and somewhat localized in p. A convenient example of such a state is the Gaussian wave packet |b; x_{0}, k_{0}⟩, defined by

⟨x|b; x_{0}, k_{0}⟩ = A e^{ik_{0}·x} e^{− \frac{(x−x_{0})^{2}}{2b^{2}}} ,  (3.2)

where A = π^{− \frac{D}{4}} b^{− \frac{D}{2}} is a normalization constant. The probability density of observing a Gaussian wave packet at a specific position forms a Gaussian distribution of width^{2} \frac{b}{\sqrt{2}}, centered around x_{0}. Since the momentum representation is

⟨k|b; x_{0}, k_{0}⟩ = A′ e^{ik·x_{0}} e^{\frac{−(k−k_{0})^{2}b^{2}}{2}} ,  (3.3)

the probability density of observing a specific momentum also forms a Gaussian distribution, whose width is \frac{1}{b\sqrt{2}}. As a result, decreasing b makes the packet better localized in position at the expense of making it less localized in momentum, and vice versa.

2The the standard deviation of each of x_{i} is σ_{i} = b/\sqrt{2}, while the standard deviation of |x| is given by σ_{|x|} = b\sqrt{D/2}.

3.3 Momentum inversion and time reversal

An interesting (and useful) symmetry is related to the fact that changing the sign of k_{0} is reverses the direction of time-evolution. More formally, we assert that two states with opposite momenta propagate as

⟨x| Û(t) |b; x_{0}, −k_{0}⟩ = ⟨x| Û(−t) |b; x_{0}, k_{0}⟩^{∗} .  (3.4)

In order to demonstrate this, let us assume that the stationary states have been defined such that their wave functions are real,^{3} so that ⟨n|x⟩ = ⟨x|n⟩^{∗} = ⟨x|n⟩. Then we can write

⟨x_{1}| e^{− \frac{iĤt}{ℏ}} |x_{2}⟩ = ∑_{n} ⟨x_{1}|n⟩ e^{− \frac{iEnt}{ℏ}} ⟨n|x_{2}⟩

= ∑_{n} ⟨n|x_{1}⟩ e^{− \frac{iEnt}{ℏ}} ⟨x_{2}|n⟩

= ∑_{n} ⟨x_{2}|n⟩ e^{− \frac{iEnt}{ℏ}} ⟨n|x_{1}⟩

= ⟨x_{2}| e^{− \frac{iĤt}{ℏ}} |x_{1}⟩

= ⟨x_{1}| e^{+ \frac{iĤt}{ℏ}} |x_{2}⟩^{∗}  (3.5)

The states with opposite momenta satisfy ⟨x|b; x_{0}, −k_{0}⟩ = ⟨x|b; x_{0}, k_{0}⟩^{∗}, so

⟨x| Û(t) |b; x_{0}, −k_{0}⟩^{∗} = (∫ ⟨x| Û(t) |x′⟩ ⟨x′|b; x_{0}, −k_{0}⟩ d^{D}x′)^{∗}
= ∫ ⟨x| Û(−t) |x′⟩ ⟨x′|b; x_{0}, k_{0}⟩ d^{D}x′
= ⟨x| Û(−t) |b; x_{0}, k_{0}⟩ ,  (3.6)

which is equivalent to equation (3.4).

One of the easily deducible implications of this idea is that the return probabilities are symmetric under momentum reversal, since

3This will be possible in any system that is invariant under time-reversal. If the time-reversal invariance does not hold, the imaginary unit enters into the wave functions of stationary states by more than just a constant phase factor, and therefore those wave functions cannot usually be made real.

⟨b; x_{0}, −k_{0}| Û(t) |b; x_{0}, −k_{0}⟩ = ∫ ⟨b; x_{0}, −k_{0}|x′⟩ ⟨x′| Û(t) |b; x_{0}, −k_{0}⟩ d^{D}x′

= ∫ ⟨b; x_{0}, k_{0}|x′⟩^{∗} ⟨x′| Û(−t) |b; x_{0}, k_{0}⟩^{∗} d^{D}x′ = ⟨b; x_{0}, k_{0}| Û^{†}(t) |b; x_{0}, k_{0}⟩^{∗}

= ⟨b; x_{0}, k_{0}| Û(t) |b; x_{0}, k_{0}⟩ .  (3.7)

Of course, this is a very convenient symmetry, since it allows us to reduce the amount of necessary computation by a factor of two.

Limitations of a discrete model

4.1 A formal problem: aliasing

As I mentioned in section 2.3, the indices of Φ_{m_{1}...m_{D}} are cyclic, so in the discrete representation the space of momenta has the topology of a torus rather than that of a Euclidean space. What causes this behavior and what are its effects?

On the purely formal level, the cyclic nature of the momenta stems from the fact that the exponentials of equation (2.14a) are periodic. However, there is a more fundamental reason having to do with the fact that Ψ_{l_{1}...l_{D}} is discretely sampled. Let us consider a wave function of the form Ψ(x) = e^{ik_{0}x_{n}} (for some fixed index n). This is a momentum eigenfunction with the momentum pointed in the direction of the nth coordinate axis. Its discrete counterpart, along the lines of equation (2.13), is

Ψ_{...l_{n}...} = e^{i k_{0} l_{n}δ_{n}},  (4.1)

where all indices other than l_{n} are immaterial. If the wave number k_{0} is increased by \frac{2π}{δ_{n}}, the discrete values become

Ψ′_{...l_{n}...} = e^{i(k_{0}+ \frac{2π}{δn})l_{n}δ_{n}} = e^{i(k_{0} l_{n}δ_{n}+2πl_{n})} = e^{i k_{0} l_{n}δ_{n}}.  (4.2)

There is no change in the sampled values of Ψ! If the phase of Ψ changes by some amount between two successive points, the phase of Ψ′ will change by that amount plus one full period, and since we have no way of examining the data between the successive points, the extra period will not change the value of the wave function. This makes it clear that discretization inherently removes the distinction between momenta whose components differ by \frac{2π}{δ_{i}}ℏ. All of such momenta will contribute to the same value of Ψ_{l_{1}...l_{D}} and become indistinguishable from each other, which is known as aliasing.

This is a serious issue, since those differing momenta contribute very differently to the wave function propagation. There are two basic ways of resolving the conflict: one can ensure that only a single one of the set of indistinguishable momenta can be available to the system, or one can abandon the conventional kinetic-energy operator K̂ = p̂^{2}/2m and adopt one that is consistent with the periodic boundary conditions mentioned above.

Choosing a different kinetic-energy operator is an intriguing proposition, and it would perhaps be particularly interesting to consider an “infinitely mirrored” system for which both coordinates and momenta have periodic boundary conditions. However, the drawback is that this kind of periodic boundary conditions isn’t commonly studied in examinations of chaotic systems in classical mechanics. This means that any serious interpretation of such a quantum-mechanical system would require a fullfledged examination of the behavior of the corresponding periodic classical system, which would be at least as involving as the quantum-mechanical analysis itself. Since my work had been progressing in a different direction, diverting my attention to a study of periodic Hamiltonians would have required more time than I had available.

For a non-periodic K̂, each value that can be taken by an index of Φ_{m_{1}...m_{D}} must be made to represent only one of the possible corresponding values of the momentum. The simplest way to achieve this is to require Φ(k) to be approximately zero everywhere outside some specific range of momenta. If the ratio δ_{i}/ℏ is then adjusted to satisfy

\frac{δ_{i}}{ℏ} < \frac{2π}{(p_{i})_{max} − (p_{i})_{min}} ,  (4.3)

then any two accessible values of p_{i} can only differ by less than \frac{2π}{δ_{i}}ℏ, which means that no aliasing will occur.

The astute reader has probably noticed that, just as the transform (2.14a) makes the indices of Φ cyclic, the opposite transform (2.14b) makes the indices of the resulting Ψ cyclic. This implies that Ψ needs to be rescued from an analogous aliasing problem. Since the distance between the two adjacent points in the ith dimension is δ_{i}, the indistinguishable points will be separated by multiples of N_{i}δ_{i}. (The same result can be obtained by a more rigorous analysis parallel to the one performed for Φ.) Thus, if the positions are restricted to a limited area, the condition

δ_{i} > \frac{(x_{i})_{max} − (x_{i})_{min}}{N_{n}}  (4.4)

ensures that any value of each of the indices corresponds to exactly one possible position.

For any bound state in a conservative potential, the positions and momenta are already restricted by the fact that the total energy of the system is conserved. The wave function may tunnel into the classically forbidden regions of either positions or momenta, but it will decay exponentially and eventually become insignificantly small. Thus, given a conservative system whose energy states are all bound, and given some fixed value of its total energy, we can find a set of ranges of positions and momenta outside of which the position-representation and momentum-representation wave functions practically vanish. The equations (4.3) and (4.4) then provide the conditions required to avoid aliasing problems.

4.2 A physical problem: good sampling

The discussion in the previous section is based solely on the mathematical properties of discrete Fourier transforms, and it merely insures the self-consistency of the computational model. It turns out that there are additional physical requirements which must be considered, and which give rise to constraints that are stronger than the conditions of the previous section.

Given the physical nature of our model, positive momenta are equally likely as the negative ones. It is therefore most reasonable to assume that the available range of momenta is centered around p = 0, i.e. that (p_{i})_{min} = −(p_{i})_{max}. Without any loss of generality we can also assume that the available range of coordinates is also centered around x = 0, since we are free to choose the origin for the potential.

If a wave function is localized in momentum, then its phase in coordinate representation precesses with the wavelength approximately equal to the de Broglie wavelength

λ = \frac{2πℏ}{p}.  (4.5)

If we compare this with condition (4.3), we see that the condition is equivalent to

δ_{i} < \frac{πℏ}{p_{max}} = \frac{λ_{min}}{2} ,  (4.6)

i.e. the distance between the grid points is allowed to be as large as half of the shortest possible wavelength. Since the real and the imaginary part of the wave function are periodic with a period that is approximately equal to the wavelength, this requirement is clearly not strong enough, since taking only two points for each period of the function is almost guaranteed to lose most of the meaningful data. A commonly used “rule of thumb” for this kind of a numerical model of a wave function is that the grid spacing must be no larger than – λ = λ/2π, which results in a stronger condition

\frac{δ_{i}}{ℏ} < \frac{1}{p_{max}} .  (4.7)

Of course, since the state under consideration is also localized in position, the same problem will occur for the wave function in momentum representation, which will precess with the “wavelength”

Λ = \frac{2πℏ}{x}  (4.8)

where Λ has units of momentum (not wave number). An analogous chain of reasoning leads to the condition

δ_{i} > π \frac{(x_{i})_{max} − (x_{i})_{min}}{N_{n}} .  (4.9)

The same result can be achieved by other means, for example by considering the quantization of the area in phase space bounded by any classical trajectory. This illustrates that it is a fundamental requirement linked to the fact that the wave function that resembles a classical particle is localized in both position and momentum.

The two new conditions, (4.7) and (4.9), are both much stronger than the conditions that come from the mathematical properties of discrete Fourier transforms. This means that their restrictions apply within the restrictions of the discrete grid; in either representation, large areas of the grid will be inaccessible to the wave function in order to insure adequate sampling in the other representation. In each dimension, less than 1/π of the range covered by the grid is actually available. This may at first seem like a waste of good grid points, but it is impossible to avoid: it is a fundamental limitation of discretization, in position and momentum representations, of a wave function which is localized in momentum and position.

Decisions about computation

5.1 What to calculate

The return probabilities for a set of Gaussian wave packets can obviously be examined in a large number of possible ways. However, the available computational hardware consisted of mid-range workstations, which limited my options considerably.

The return probability of a Gaussian wave packet depends on the values of t, b, x_{0} and k_{0}. This means that the parameter space to be explored has 2D +2 degrees of freedom, where D is the dimension of the model. Since one-dimensional conservative problems are exactly integrable, the smallest interesting dimensionality for study of chaos is D = 2. However, six independent parameters is still too many to be practical given the computational constraints, and some meaningful subset of the parameter space needs to be selected for analysis.

The ultimate purpose of our model is to compare the behavior of return probability values with the data about classical periodic orbits. Of the six parameters, b clearly has the least interesting influence on the return probability. Therefore, we can just set it to some sensible fixed value that ensures reasonable behavior of the computational model.

Furthermore, we note that energy is conserved across any classical periodic orbit of a conservative system. This means that the return probabilities should be considered together with other return probabilities for the same energy. Consequently, if we need to examine only a portion of the full parameter space, it makes sense to consider only one fixed value of the classical energy.

Fixing the value of energy means that the absolute value of k_{0} is uniquely determined by the value of x_{0}:

|k_{0}| = \frac{1}{ℏ}\sqrt{2 m (E_{classical} − V (x_{0}))}.  (5.1)

(The symbol E_{classical} serves as a reminder that we are not talking about the expectation value of the energy of the wave function, ⟨ψ| Ĥ |ψ⟩, but about the classical energy of the associated point-like classical particle, p^{2}/2m+V (x).) The wave vector is then a function of direction alone, which can in the two-dimensional case be represented as an angle α such that k_{x} = |k_{0}| cos α and k_{y} = |k_{0}| sin α.

This leaves us with four independent parameters,^{1} x_{0}, y_{0}, α and t. While this is still relatively many, the values of those parameters can be explored reasonably well within the constraints described above.

5.2 Choosing the potential

The discussion of the previous section shows the fact that working with a twodimensional system was my only reasonable choice. However, there are many interesting two-dimensional potential functions that give rise to chaotic behavior, and not all are well-suited to the method of propagation used here.

This point is nicely illustrated by the Nelson potential,

V = μ\frac{x^{2}}{2} + (y − \frac{x^{2}}{2} )^{2} ,  (5.2)

where the parameter μ is set to 0.1. As fig. 5-1 shows, the contours of Nelson potential are horn-shaped and relatively narrow, growing longer and increasingly Figure 5-1: Contour plot of the Nelson potential Figure 5-2: Contour plot of the Delfi potential curved as the energy increases. This, however, means that the ranges of available x and y coordinates quickly become large, compared to the relatively small area of the classically allowed region. Consequently, the width of the “belt” of the allowed region at x = 0 becomes only a few times larger than the size of the grid that can be supported by the available computer hardware, and the quality of the simulation suffers as a result.

1Of course, the four coordinates of the parameter space could be chosen in a number of other ways. This one will prove convenient, since it allows plotting the return probabilities in a particularly illuminating manner.

This shows that the ideal potential for this model of discrete computation will stay in a compact, roughly rectangular area throughout the “interesting” range of energies. With this in mind, I chose to examine the Delfi potential,

V = \frac{1}{4} (x^{2} + y^{2})^{2} + (\frac{2}{3})^{\frac{3}{2}} x^{3}(y − 1) − y (x^{2} + y^{2}) + \frac{1}{2} (x^{2} + 3y^{2}) , (5.3)

depicted in fig. 5-2. Clearly the geometry of this potential fits a rectangular region much better. The classical behavior of the Delfi potential was first investigated by Marios Kagarlis in his bachelor’s thesis [Kagarlis, 1986].

5.3 Choosing parameters

The final task remaining before the start of actual computation is to decide upon the values of various computation parameters. Since the reasonable values for many of the parameters depend in part on the values of other parameters, the most fundamental ones need to be considered first.

5.3.1 Energy

In order to obtain interesting results, the value for energy must be chosen carefully. An invaluable tool in this endeavor is the E–τ plot of the periodic orbits of the Delfi potential (fig. 5-3), reproduced from [Kagarlis, 1986]. The E–τ plot identifies the classical families of periodic orbits whose influence is likely to exhibit itself in the return probabilities for some value of the energy. It makes it possible to avoid the energies where the period of one of the families is close to an integral multiple of the Figure 5-3: The E–τ plot for the Delfi potential, reproduced from [Kagarlis, 1986]

\begin{matrix}c \\ ◯ \end{matrix} 1986, Massachusetts Institute of Technology. Reproduced with permission.

period of another, which makes the features of either hard to distinguish. The plot also distinguishes between the stable and unstable trajectories, allowing us to choose an energy for which both stable and unstable trajectories are to be found and their effects on return probabilities compared.

However, an analysis of the E–τ plot is not sufficient. Even if the periods of two trajectories are extremely close, one may be able to distinguish their return probabilities if the shapes of the trajectories are sufficiently different. I therefore also studied a number of Poincaré plots and periodic trajectory plots for various energies.

Finally, for some energies the quantum-mechanical propagation has undesirable features that are hard to observe in an examination of the classical behavior. For example, if a periodic orbit is not only unstable but very strongly divergent,^{2} a wave packet that starts out resembling a particle on this trajectory may spread out too rapidly to generate an identifiable maximum. This suggests that it would also be worthwhile to examine numerically the time evolution of a wave packet for a few of the more interesting sets of initial conditions, such as those corresponding to periodic orbits.

The combination of the methods outlined above suggests that the two most interesting regions for a quantum-mechanical comparison lie at E_{classical} ≈ 0.3 and E_{classical} ≈ 0.9. Unfortunately, available time did not permit me to generate returnprobability data for both of these regions, and I ultimately decided to work with E_{classical} = 0.3.

Two Poincaré sections for E = 0.3 are shown in fig. 5-4 and 5-6. Since the points corresponding to unstable periodic orbits are hard to locate in a Poincaré plot, all periodic orbits have been identified in fig. 5-5 and 5-7. In the x = 0 section, the vertically oscillating trajectory which belongs to Family 1 is spread all around the energy boundary of the plot. Finally, a stable period quadrupling of Family 4 (which is not present in the E–τ plot) is visible in the Poincaré plots; however, its period is too large for it to be relevant in our analysis, so it wasn’t labelled.

2Divergent in the sense that a small perturbation grows rapidly with time, i.e. the Lyapunov exponent in the neighborhood of the trajectory is large.

Figure 5-4: The (x = 0, p_{x} > 0) Poincaré section of the Delfi potential for E = 0.3 Figure 5-5: Locations of periodic orbits in fig. 5-4, identified by family Figure 5-6: The (y = 0, p_{y} > 0) Poincaré section of the Delfi potential for E = 0.3

Figure 5-7: Locations of periodic orbits in fig. 5-6, identified by family

5.3.2 Grid size

The size of the grid is limited exclusively by the power of the available computing resources. A larger grid will clearly provide better accuracy, but at some point the calculation time becomes prohibitively slow. I have found that a 256x256 grid is the largest reasonable size for a piecemeal calculation distributed across a number of mid-range workstations. (Of course, given the current rate of speed improvements, this will probably be out of date by the time you read it.)

5.3.3 Grid spacing

The choice of energy determines the range of classically allowed positions for which V (x, y) < E. Given the end points of that range, the grid spacing must satisfy the previously mentioned condition

δ_{i} > π \frac{(x_{i})_{max} − (x_{i})_{min}}{N_{n}} .  (4.9)

Within this constraint, a smaller value of δ_{i} is preferable.

For the sake of simplicity, I decided to use the same value of δ in both directions. Since the Delfi potential is wider than it is tall at the range of energies under consideration, δ is more constrained by x then by y. For E = 0.3, the extents are

x_{min} = −0.606526,  x_{max} = 1.250443.  (5.4)

Inserting these values into equation (4.9) determines that the the constraint is narrowly satisfied when δ = 0.023.

5.3.4 Planck constant

After we have fixed the value of δ, we can rewrite equation (4.7) to express a requirement for ℏ alone,

ℏ > δ p_{max} = δ \sqrt{2mE}.  (5.5)

Again, a lower ℏ is preferable as long as this constraint is satisfied: since the size of the discretization steps in momentum representation is proportional to ℏ, smaller values mean that a larger number of grid points will cover the region where the wave function is not negligibly small.

The customary convention of setting m = 1 will be used here. A numeric value for ℏ is obtained by inserting numbers into the inequality and setting ℏ to barely satisfy the requirement, which results in ℏ = 0.0179.

5.3.5 Gaussian packet width

The discrete representations of a Gaussian wave packet will be the most accurate when its width in the position representation takes up about as many discrete grid points as its width in the momentum representation. The spacings of the two grids are δ and \frac{2π}{Nδ}, respectively. According to equations (3.2) and (3.3), the width of the position-representation wave function for a Gaussian wave packet is \frac{b}{\sqrt{2}}, and its width in the momentum representation is \frac{1}{b\sqrt{2}}. Putting all of this together, we get

\frac{\frac{b}{\sqrt{2}}}{δ} = \frac{\frac{1}{b\sqrt{2}}}{\frac{2π}{Nδ}}  (5.6)

and, after solving for b, we obtain b = δ\sqrt{\frac{N}{2π}} = 0.104\sqrt{2} = 0.147.

5.3.6 Time step

Finally, one must select a value for the time step Δt. In order for the approximation of equation (2.9) to be accurate, one would suspect that Δt needs to be small compared to ℏ/E. But how small is “small”? There is no particular rule to be applied here, as the answer likely depends on the nature of the chosen potential. Instead, in order to determine if a particular value Δt_{?} was acceptable, I calculated the time evolution for a few trajectories, using the parameters chosen above, for Δt = Δt_{?} and Δt = \frac{Δt_{?}}{2} . If the difference between these two calculations decreases slowly, then the error in each of the calculations can be assumed to be, at worst, not much larger than the difference itself.

For Δt = 0.01 and Δt = 0.05, I found that the relative differences between the amplitudes were always smaller than 10^{−3}, after a time-evolution of duration T = 10. This is not great, but it isn’t particularly bad, either: if we make the conservative assumption that the relative error in amplitude can be as large as 5 · 10^{−3}, the errors in return probabilities due to the size of the time step will still be smaller than 1%. This is very likely smaller than the error due to the fact that the wave function is discretized. Therefore, Δt = 0.01 is small enough to be used in the final calculations.

5.4 A sample evolution

The evolution of a sample Gaussian wave packet, using the parameters selected here, is plotted in fig. 5-8—5-10. The packet parameters used were x = 0, y = 0, α = tan^{−1}(p_{y}/p_{x}) = 0.

Figure 5-8: Time-evolution from a sample initial state 37 Figure 5-9: Time-evolution from a sample initial state (continued)

38 Figure 5-10: Time-evolution from a sample initial state (continued)

39

Analysis of results

6.1 Return probability plots

A convenient way to plot the return probabilities is to gather together the return probabilities that share the same t and x_{0}, but whose momentum is oriented in different directions. If we made a polar plot of those probabilities, the directions with high probabilities would stand out from the rest, forming “spikes” pointed in the directions in which the return is most likely. Instead of making just one plot, we can draw around each point x_{0} a polar plot of return probabilities for wave packets originating from that point, which will simultaneously display all of the return probability data for a specific E_{classical} and t.

The primary advantage of doing this is that the maxima of the polar plots give the entire figure a sense of directional continuity similar to the effect of the arrows of a quiver plot. This visual grouping corresponds to a group of wave packets whose return probabilities for the specified propagation duration are high, and whose centers approximately lie on a curved line. (For an example, see fig. 6-8.) Our primary interest, motivated by discussion of chapter 3, is to examine how closely these groupings correspond to the classical periodic orbits for the values of t when the the return probabilities reach their local maxima.

6.2 Return probabilities of the Delfi potential

After writing software to implement the computational model described in the preceding chapters, I computed a set of values of return probabilities for the Delfi potential at the energy E_{classical} = 0.3, and the values of the other parameters chosen as described in section 5.3. All of the results are plotted in Appendix A. Let us examine the features of those plots and compare them with the data obtained for the classical problem.

6.2.1 Family 1

At t = 0, the values of the return probabilities are actually all equal to 1, since the packet hasn’t yet evolved at all. As the packet is evolving, they rapidly become very small. For times below t = 3.0, packets haven’t been able to propagate far enough to get back near their original position, so the return probabilities stay close to zero. After t = 3.0, however, some of the probabilities begin to increase, reaching

Figure 6-1: Return probabilities for t = 3.8 a local maximum for t = 3.8. The data displayed in fig. 6-1 is fairly easy to interpret. All of the maxima of return probabilities are exactly vertical, and their values are highest at x = 0. This suggests that any wave packet with E = 0.3 and x ≈ 0 that is initially moving in the vertical direction will continue to oscillate in the vertical direction almost unchanged, spreading out relatively slowly, until the spreading makes it indistinguishable.

An examination of the classical behavior of the Delfi potential indicates that one of the families of the periodic orbits, identified as Family 1 in [Kagarlis, 1986], does in fact consist of vertical oscillations with x = 0. (The trajectories of Family 1 do stay vertical for any energy, since \frac{∂V (x,y)}{∂x} |_{x=0} = 0 for all y.) The E–τ plot shows that the family is unstable for E = 0.3, and confirms that the period is close to the return time of the quantum model. The trajectory for E = 0.3 is very easy to find, since it is known to be vertical, and its period is found to be τ = 3.865, which is fairly close to the return time of the associated maximum of return probability.

6.2.2 Family 2

As the return time is increased above t = 3.8, the return probabilities decrease again, practically vanishing at t = 4.5. The plots start showing visible data again slightly below t = 6. The return probabilities increase more slowly than previously, and the shape exhibited by the return probability plots changes noticeably, which should be remembered for later analysis.

The first clear maximum of the return probabilities occurs for t = 7.4. According to the E–τ plot, this corresponds to Family 2, which starts as small horizontal oscillations at low energies. Since its trajectories aren’t trivial like those of Family 1, finding its trajectory for E = 0.3 is somewhat harder. One possible solution involves the monodromy matrix method, which is how the E–τ plot was generated [Kagarlis, 1986]. A somewhat more brute-force approach can be made using Poincaré plots. Since Family 2 is unstable, it cannot be easily located on a Poincaré plot for E = 0.3, but it is easy to find it at a lower energy where it is stable and follow its movement as the energy changes. The trajectory for E = 0.3, whose period is

Figure 6-2: Return probabilities for t = 7.4

Figure 6-3: A plot of the classical Family 2 trajectory at E = 0.3 superimposed on the return probabilities for t = 7.4 τ = 7.451, is displayed in fig. 6-3 along with the return probabilities.

In examining fig. 6-3, one finds that the return probabilities match the periodic orbit reasonably well in the left part of the potential: the maxima lie close to the trajectory, and the direction of the probabilities is almost tangential to the curve. On the right, however, the curve and the maxima end up visibly apart.

A part of the reason for this is optical: at the points for which a small perpendicular element crosses the main one, the main element is darker and therefore looks longer. However, the region where the return probabilities are large certainly does extend much further to the right than the classical trajectory would suggest. Can we explain what causes this?

A look at the E–τ plot provides a possible answer: in the region of interest, Family 2 has a period doubling, Family 9. At E = 0.3, the period of Family 9 is still close to twice the period of Family 2. If the the two doubled “halves” of Family 9 lie close enough to each other, a packet which is initially on one of them could overlap its

Figure 6-4: A plot of the classical Family 9 trajectory at E = 0.3 superimposed on the return probabilities for t = 7.4 (along with Family 2).

Figure 6-5: Time-evolution starting from the “narrow” bifurcation of Family 9 The points between t = 0 and t = 4.2, when the packet moves to the left, turns and returns on the “wide” branch, are not shown.

Figure 6-6: Time-evolution starting from the “wide” bifurcation of Family 9 The points between t = 0 and t = 4.2, when the packet moves to the left, turns and returns on the “narrow” branch, are not shown. propagated counterpart which is, half a period later, following the other “half”.

In order to verify if this guess is correct, we must actually find Family 9, which can be done in an analogous way as with Family 2. The period of the orbit, which is also unstable, is τ = 14.877, and the trajectory itself is displayed in fig. 6-4.

The geometry of Family 9 does explain the spread of the maxima. The explanation is not as complete as one might wish, since the narrower bifurcated part of the family seems to influence return probabilities more than the wider one. In order to explain this, we need to examine the evolution of the wave functions that correspond to classical Family 9.

Packets that start in the “narrow” branch of the bifurcation propagate as would be expected, so that at t = 7.4 they are localized at approximately the corresponding point in the “wide” branch (fig. 6-5). However, the packets starting in the “narrow” branch initially spread much faster. At t ≈ 6.6, when the packet is being reflected from the side of the potential, it is almost uniformly spread over a long segment. As a result, after the turning point, interference can be observed in the wave function. In this case, the values of the wave function close to its original location increase dramatically, so that the return probability for t = 7.4 is found much higher than it would otherwise be.

This demonstrates that, while the correspondence between return probabilities and the classical periodic orbits is usually very good, there are cases where quantal phenomena—in this case, spreading and interference—affect the return probabilities in a very visible way, creating marked differences between cases which would be expected to be similar if one were to rely on solely intuition about classical mechanics.

6.2.3 Family 1 revisited

As the return probabilities associated with Family 2 decrease for t > 7.4, the evidence of two other families begins to appear. The first of them is a recurrence of Family 1 (after two full orbits), which peaks at t = 7.7. This is, as expected, very close to twice the period found classically for Family 1, 2τ = 7.731. The vertical elements of the plot (fig. 6-7) once again follow the shape of the classical trajectory; the curved

Figure 6-7: Return probabilities for t = 7.7

ones are a part of an unrelated family.

6.2.4 Family 4

The curved flow visible in the background of fig. 6-7 reaches its own maximum at t = 8.0. This maximum can be identified with Family 4, which is a period doubling of Family 1, but which at this energy has a very different shape.

Unlike the other families, Family 4 is stable for E = 0.3, and is therefore easy to locate on a Poincaré plot. The trajectory has a period of τ = 8.030, and matches the return probability data extremely well.

Figure 6-8: Return probabilities for t = 8.0

Figure 6-9: A plot of the classical Family 4 trajectory at E = 0.3 superimposed on the return probabilities for t = 8.0

6.2.5 Family 1, once again

After t = 8, the return probabilities once again decline to near-zero values. The only remaining feature that emerges before t = 12, where my calculation ended, was the third appearance of Family 1. Once again, the maxima of the return probabilities are aligned vertically and grouped around x = 0. The value of the return probabilities has decreased since the second appearance of this family, just as it has decreased between the first and second ones; since the wave packet is spreading, the overlap decreases with each orbit. However, the third maximum is still easily noticeable, particularly so since no other features exhibit themselves in this time interval.

Figure 6-10: Return probabilities for t = 11.6

6.2.6 A shadow of Family 2

Let us now return our attention to the time interval between t = 6.5 and t = 7.3. As the probabilities increase, the maxima arrange themselves into a shallow horizontal shape, but at t ≈ 7.2, the shape suddenly transforms to become much more similar to Family 4. A careful examination shows that the probabilities of the original shape do in fact exhibit a local maximum for t = 7.0. But where does this maximum come from? There are no other periodic trajectories to be found with the appropriate period or the shape similar to the one displayed in fig. 6-11.

Once again, the E–τ plot provides a clue. While it displays no trajectories with period close to 7.0 for E = 0.3, it shows that Family 2 reaches a minimum in energy for τ ≈ 7.0, with the energy of E ≈ 0.4. Clearly, this energy would be inaccessible in the classical model, but wave packets are composed of many energy eigenstates and therefore could be found at a range of energies, just as they could be found at a range of positions or momenta.

I tried to look for a resemblance between the return probabilities and the classical trajectory corresponding to the local energy minimum of Family 4, which occurs for E ≈ 0.3943 and τ ≈ 6.948. The result, shown in fig. 6-12, shows a very strong similarity between the two.

We can verify that the spread in energy for a Gaussian wave packet with E_{classical} = 0.3 is large enough to cover the energy of the classical trajectory, E = 0.3943. While we can’t easily recover the probability distribution in energy (which is a decomposition of the packet into energy eigenfunctions), it is easy to calculate both the average of the energy, E_{0} = ⟨ψ| Ĥ |ψ⟩, and its variance, σ_{E}^{2} = ⟨ψ| (Ĥ − E_{0})^{2} |ψ⟩, by numerically applying Ĥ = K̂ + V̂. For one of the test packets used in calculating the return probabilities,^{1} the values obtained are E_{0} = 0.32, σ_{E}^{2} = 5.8 · 10^{−3}, and σ_{E} = 0.076. Thus, the minimum of Family 4 is only one standard deviation away from the mean energy, which means that our explanation makes sense. It must be concluded that the features observed in return probabilities can be caused by classical behavior that is classically inaccessible to the wave packet under consideration!

1Specifically, the parameters of the packet that was used were E = 0.3, x = 0, y = 0.18, α = tan^{−1}(p_{y}/p_{x}) = 30^{∘}.

Figure 6-11: Return probabilities for t = 7.0

Figure 6-12: A plot of the classical Family 4 trajectory at E = 0.3943 superimposed on the return probabilities for t = 7.0

6.3 Conclusions

The material in this chapter demonstrates that, at least in the case under consideration, there exists a strong correspondence between the local maxima of the return probabilities of a single-particle quantum-mechanical system and the periodic orbits of the analogous classical system. However, since wave packets do not possess a single, definite value of the energy, the return probabilities can in a smaller degree reflect the behavior of the classical system at a different energy. The quantum-mechanical system also exhibits certain kinds of behavior which are not present in its classical counterpart, such as packet spreading, interference and tunneling.

Another point of note is that there is relatively little observed difference between the influence of stable and unstable orbits. The areas with relatively large return probabilities are found to be slightly narrower for the stable orbits (Family 4) than for unstable ones, but the difference is slight. The packets corresponding to strongly unstable periodic orbits may have their probability shift to other periodic orbits (Families 9 and 2), but this effect—if present at all—wasn’t found to be strong. In general, the quantum-mechanical system seems to be much less affected by stability. (Unfortunately, this is unlikely to make quantum-mechanical models useful for study of unstable families of classical orbits, because good results can be obtained with much less computationally intensive classical models, such as the monodromy method.)

In general, I found that using explicit propagation of discretized wave functions works reasonably well as a means of calculating return probabilities for Gaussian wave packets. While relatively computationally expensive, it would probably be particularly useful when computing the stationary states is particularly hard, or in the cases where a fast hardware-based way of computing discrete Fourier transforms is available.

Appendix A

Return probability plots for the

Delfi potential with E=0.3

The following values of parameters were used to generate the plots included in this appendix:

E_{classical} = 0.3
N_{1} = N_{2} = 256
δ_{1} = δ_{2} = 0.023  [grid spacing]
ℏ = 0.0179
b = 0.104\sqrt{2} [wave packet width parameter]
Δt = 0.01
m = 1

The points for which the return probabilities are calculated are spaced by 0.06 in both x and y directions, and the directions of momentum are spaced apart by 10^{∘}. Return probability plots were generated for times separated by 0.1 = 10 Δt.

All return probability data was scaled up to bring out the features that would otherwise be too small to see. The data was multiplied by a scaling factor, f, times one half of the distance between plot points (in this case, 0.03). However, the areas where the data is small do not become more interesting when f is increased, so all of the following plots use the same value, f = 1.37. Most of the blank plots have been omitted.

55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123

Appendix B

Software used for wave-function

propagation

B.1 Choosing the programming language

MIT’s public computing system, named Athena, consists of about five hundred workstations in a dozen computer clusters that are made freely available to the MIT community. The software used for calculation of return probabilities was intended to be able to be run on them. Unfortunately, FORTRAN compilers were available only for the oldest, slowest of these machines at the time when I started writing my simulation software. This made using FORTRAN impractical.

One notable disadvantage C in numerical modeling of quantum mechanics is the lack of availability of complex numbers. Of course, it isn’t hard to explicitly write out the complex-number operations in terms of their real and imaginary parts, but it makes longer expressions nearly unreadable and typing errors almost impossible to find. How can this be avoided?

Some C compilers provide complex numbers as an extension to the standard C language. The one example available to me is gcc, the GNU C compiler. Gcc is widely used on UNIX systems, comparatively easy to port, and free. Most importantly to me, it has been made available on all current Athena platforms.

Thus, gcc seemed the compiler of choice for working on Athena. However, relying on it would cause problems for using my code in a different environment. A possible alternative would be to choose C++, where I could implement complex numbers as an object class and define arithmetic operations for them myself. (Some C++ compilers already provide a complex class, but creating my own would help both performance^{1} and portability.) C++ isn’t as widely used as C, but this approach will work with any C++ compiler.

The two approaches described above are not mutually exclusive. C++ is closely enough related to C that it’s easy to make most ANSI C code compile under a standard C++ compiler. It is easy to design a few simple syntactical conventions that can be used as props to access the underlying implementation of complex numbers. These props will be implemented completely differently for C++ than for extended C, but the differences can be restricted to a header file. The rest of the software can then simply use a unified set of conventions that will work under both gcc and C++.

B.2 Complex number conventions

The following are the conventions I decided to implement and use in my code, in order to hide the language-dependent details of the implementation of complex numbers:

• The complex numbers should be declared as type complex. Addition, subtraction, multiplication, division and negation (unary minus) applied to numbers of this type use complex-number arithmetic.

• To simplify changes required to adapt the code to more powerful numeric types (such as long double), type real should be used as the real number type corresponding to complex. (Currently, real is always an alias for double.)

• The symbol I represents the imaginary unit. (Lower-case i is not used to avoid conflicts with a commonly used index variable.)

1In at least one case, the compiler’s complex class failed to specify rules for complex-times-real multiplication, which means that the compiler must generate code for real-to-complex conversion and a complex multiplication. In theory, the inefficiency could get fixed by a compiler’s optimization pass– but in this case it doesn’t.

• Real and imaginary parts of an expression can be extracted with Re(expr) and Im(expr), respectively.

• The complex conjugate of an expression can be written as conj(expr).

In addition, several handy functions for common complex-number operations are defined as follows:

Cexp(z) returns e^{z}.
Cabs(z) returns |z| = \sqrt{(Re z)^{2} + (Im z)^{2}}.
Csqabs(z) returns |z|^{2} = (Re z)^{2} + (Im z)^{2}.
Csqrt(z) returns \sqrt{z}.
Carg(z) returns arg z = tan^{−1} (\frac{Im z}{Re z} ).

Bibliography

[Baranger, 1993] Baranger, M. (1993). A calculation of return probabilities by decomposing the wavefunction into eigenfunctions of energy. Currently unpublished notes.

[Berge et al., 1984] Berge, P., Pomeau, Y., and Vidal, C. (1984). Order within chaos: towards a deterministic approach to turbulence. ?

[Brigham, 1974] Brigham, E. O. (1974). The Fast Fourier Transform. Prentice-Hall.

[Frigo and Johnson, 1997] Frigo, M. and Johnson, S. G. (1997). The Fastest Fourier Transform in the West. Technical Report MIT-LCS-TR-728, MIT Laboratory for Computer Science. (Available at http://theory.lcs.mit.edu/~fftw/, along with additional FFTW information and the source code.)

[Gasiorowicz, 1974] Gasiorowicz, S. (1974). Quantum Physics. John Wiley and Sons.

[Kagarlis, 1986] Kagarlis, M. A. (1986). Classification and topology of the families of periodic trajectories in a generic two-dimensional potential. B.S. thesis, MIT.

[Messiah, 1966] Messiah, A. (1966). Quantum Mechanics. John Wiley and Sons.

[Nussbaumer, 1982] Nussbaumer, H. J. (1982). Fast Fourier Transform and Convolution Algorithms. Springer-Verlag, 2nd edition.

[Press et al., 1994] Press, W. H. et al. (1994). Numerical Recipes in C. Cambridge University Press, 2nd edition.

[Sakurai, 1985] Sakurai, J. J. (1985). Modern Quantum Mechanics. Benjamin/ Cummings Publishing.

[Taylor, 1976] Taylor, T. T. (1976). Mechanics: Classical and Quantum. Pergamon Press.