% -*- LaTeX -*-
\chapter{Decisions about computation}

\section{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
%$\RP{b;\x_0,\k_0}(t)$,
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$:
\begin{equation}
|\k_0| = {1\over\hbar} \sqrt{2\,m \left(\E - V(\x_0)\right)} .
\end{equation}
(The symbol $\E$ serves as a reminder that we are not talking about
the expectation value of the energy of the wave~function,
$\bra<\psi|\H\ket|\psi>$, 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 $\alpha$ such that $k_x=|\k_0|\cos\alpha$ and
$k_y=|\k_0|\sin\alpha$.

This leaves us with four independent parameters,\footnote{Of 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.} $x_0$, $y_0$, $\alpha$ and $t$.  While this is still
relatively many, the values of those parameters can be explored
reasonably well within the constraints described above.

\section{Choosing the potential}

The discussion of the previous section shows the fact that working
with a two-dimensional 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,
\begin{equation}
V = \mu {x^2\over2} + \left( y - {x^2\over2} \right)^2 ,
\end{equation}
where the parameter $\mu$ is set to $0.1$.
\begin{figure}[p]
\begin{center}
\leavevmode\epsfysize=\imagesize\epsfbox{plots/nelson.eps}
\end{center}
\caption{Contour plot of the Nelson potential}
\label{fig:nelson}
\end{figure}
\begin{figure}[p]
\begin{center}
\leavevmode\epsfysize=\imagesize\epsfbox{plots/delfi.eps}
\end{center}
\caption{Contour plot of the Delfi potential}
\label{fig:delfi}
\end{figure}
As \fig{nelson} shows, the contours of Nelson potential are
horn-shaped and relatively narrow, growing longer and increasingly
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.

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,
\begin{equation}
V = {1\over4}\left(x^2 + y^2\right)^2
	+ \left(2\over3\right)^{3\over2} x^3 (y-1)
	- y \left(x^2 + y^2\right)
	+ {1\over2}\left(x^2 + 3 y^2\right) ,
\end{equation}
depicted in \fig{delfi}.  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~\cite{delfi}.

\section{Choosing parameters}
\label{sec: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.

\subsection{Energy}

\begin{figure}[p]
\begin{center}
%\leavevmode\epsfbox{plots/???}
\leavevmode\vbox to 8in{\vss}
\end{center}
\caption{The \Et\ plot for the Delfi potential,
reproduced from \protect\cite{delfi}}
\begin{center}\size{5}{0pt}\selectfont
  \copyright\ 1986, Massachusetts Institute of Technology.
  Reproduced with permission.
\end{center}
\label{fig:Et}
\end{figure}

In order to obtain interesting results, the value for energy must be
chosen carefully.  An invaluable tool in this endeavor is the \Et\
plot of the periodic orbits of the Delfi potential (\fig{Et}),
reproduced from \cite{delfi}.  The \Et\ 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 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 \Et\ 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\'e 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,\footnote{Divergent in the sense
that a small perturbation grows rapidly with time, i.e. the Lyapunov
exponent in the neighborhood of the trajectory is large.} 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 \approx 0.3$ and $\E \approx 0.9$.  Unfortunately, available time
did not permit me to generate return-probability data for both of
these regions, and I ultimately decided to work with $\E=0.3$.

\begin{figure}[p]
\begin{center}
\leavevmode\epsfysize=\imagesize\epsfbox{plots/poin/poincare30a.eps}
\end{center}
\caption{The ($x=0$, $p_x>0$) Poincar\'e section of the Delfi potential
         for $E=0.3$}
\label{fig:poincare-first}
\label{fig:poincare,x}
\end{figure}
\begin{figure}[p]
\begin{center}
\leavevmode\epsfysize=\imagesize\epsfbox{plots/poin/poincare30a_l.eps}
\end{center}
\caption{Locations of periodic orbits in \protect\fig{poincare,x},
         identified by family}
\label{fig:poincare,x;legend}
\end{figure}
\begin{figure}[p]
\begin{center}
\leavevmode\epsfysize=\imagesize\epsfbox{plots/poin/poincare30b.eps}
\end{center}
\caption{The ($y=0$, $p_y>0$) Poincar\'e section of the Delfi potential
         for $E=0.3$}
\label{fig:poincare,y}
\end{figure}
\begin{figure}[p]
\begin{center}
\leavevmode\epsfysize=\imagesize\epsfbox{plots/poin/poincare30b_l.eps}
\end{center}
\caption{Locations of periodic orbits in \protect\fig{poincare,y},
         identified by family}
\label{fig:poincare,y;legend}
\label{fig:poincare-last}
\end{figure}

Two Poincar\'e sections for $E=0.3$ are shown in \fig{poincare,x}
and~\ref{fig:poincare,y}.  Since the points corresponding to unstable
periodic orbits are hard to locate in a Poincar\'e plot, all periodic
orbits have been identified in \fig{poincare,x;legend}
and~\ref{fig:poincare,y;legend}.  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 \Et\ plot) is visible in the
Poincar\'e plots; however, its period is too large for it to be
relevant in our analysis, so it wasn't labelled.

\subsection{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.)

\subsection{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
\begin{equation}
\delta_i > \pi { (x_i)_{\text{max}} - (x_i)_{\text{min}} \over N_n }.
\tag{\ref{eq:delta,x}}
\end{equation}
Within this constraint, a smaller value of $\delta_i$ is preferable.

For the sake of simplicity, I decided to use the same value of
$\delta$ in both directions.  Since the Delfi potential is wider than
it is tall at the range of energies under consideration, $\delta$ is
more constrained by $x$ then by $y$.  For $E=0.3$, the extents are
\begin{xalignat}{2}
x_{\text{min}} &= -0.606526 , &
x_{\text{max}} &=  1.250443 .
\end{xalignat}
Inserting these values into equation~\eq{delta,x} determines that the
the constraint is narrowly satisfied when $\delta = 0.023$.

\subsection{Planck constant}

After we have fixed the value of $\delta$, we can rewrite
equation~\eq{delta,p} to express a requirement for $\hbar$ alone,
\begin{equation}
\hbar > \delta \; p_{\text{max}} = \delta \; \sqrt{2mE} .
\end{equation}
Again, a lower $\hbar$ is preferable as long as this constraint is
satisfied: since the size of the discretization steps in momentum
representation is proportional to $\hbar$, 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 $\hbar$ is obtained by inserting numbers into the
inequality and setting $\hbar$ to barely satisfy the requirement, which
results in $\hbar = 0.0179$.

\subsection{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 $\delta$ and
$2\pi\over N\delta$, respectively.  According to
equations~\eq{gaussian} and~\eq{gaussian,k}, the width of the
position-representation wave~function for a Gaussian wave packet is
$b\over\sqrt{2}$, and its width in the momentum representation is
$1\over b\sqrt{2}$.  Putting all of this together, we get
\begin{equation}
{\,\,{b\over\sqrt2}\,\, \over \delta} =
	{\,\,{1\over b\sqrt2}\,\, \over {2\pi\over N\delta}}
\end{equation}
and, after solving for $b$,
we obtain $b = \delta\sqrt{N\over2\pi} = 0.104\sqrt2 = 0.147$.

\subsection{Time step}

Finally, one must select a value for the time step $\Delta t$.  In
order for the approximation of equation~\eq{separated} to be accurate,
one would suspect that $\Delta t$ needs to be small compared to
$\hbar/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 $\Delta t_?$ was acceptable, I calculated the time evolution for
a few trajectories, using the parameters chosen above, for $\Delta
t=\Delta t_?$ and $\Delta t={\Delta t_?\over2}$.  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 $\Delta t=0.01$ and $\Delta 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\cdot10^{-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, $\Delta t = 0.01$ is small enough to be used
in the final calculations.

\section{A sample evolution}

The evolution of a sample Gaussian wave packet, using the parameters
selected here, is plotted in \fig{evol-first}---\ref{fig:evol-last}.
The packet parameters used were $x=0$, $y=0$, $\alpha =
\tan^{-1}(p_y/p_x) = 0$.

\begin{figure}[p]
\vbox{\hbox{%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_0.eps}%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_50.eps}%
}\hbox{%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_100.eps}%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_150.eps}%
}\hbox{%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_200.eps}%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_250.eps}%
}}%
\caption{Time-evolution from a sample initial state}
\label{fig:evol-first}
\end{figure}
\begin{figure}[p]
\vbox{\hbox{%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_300.eps}%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_350.eps}%
}\hbox{%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_400.eps}%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_450.eps}%
}\hbox{%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_500.eps}%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_550.eps}%
}}%
\caption{Time-evolution from a sample initial state (continued)}
\end{figure}
\begin{figure}[p]
\vbox{\hbox{%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_600.eps}%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_650.eps}%
}\hbox{%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_700.eps}%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_750.eps}%
}\hbox{%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_800.eps}%
\epsfxsize=.5\textwidth \epsfbox{plots/generic/prop_850.eps}%
}}%
\caption{Time-evolution from a sample initial state (continued)}
\label{fig:evol-last}
\end{figure}
