% -*- LaTeX -*-
\chapter{Modelling the dynamics of a quantum-mechanical system}

\section{The time evolution operator}

The dynamical behavior of the state of an isolated quantum-mechanical
system is described by the Schr\"{o}dinger equation,
\begin{equation}
i\hbar{\partial\over\partial t} \kpsi(t)= \H \kpsi(t) .
\label{eq:tdse}
\end{equation}
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~\eq{tdse} 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$:
\begin{equation}
\kpsi(t+dt)
	= \kpsi(t) - {i\over\hbar}\H dt \kpsi(t)
	= \left(1 - {i\over\hbar}\H dt\right) \kpsi(t) .
\end{equation}
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:
\begin{equation}
\U(t)
= \lim_{N\rightarrow\infty}
      \left( 1-{i\over\hbar}\H {t\over N} \right)^N
= \exp\left( -{i\H t\over\hbar} \right) .
\label{eq:propagator}
\end{equation}
The operator $\U$ is commonly referred to as the time-evolution
operator\footnote{Sometimes the term {\em propagator\/} is also used,
although the word is more commonly used in the literature for the
kernel $K(\x',t' ; \x,t) = \bra<\x'| \U(t'-t) \ket|\x>$ of the $\x$
representation of the evolution operator written in integral form.} since,
when applied to an initial ket, it traces out the dynamical evolution of
that state in time:
\begin{equation}
\kpsi(t) = \U(t-t_0) \kpsi(t_0).
\end{equation}

Since $\H$ is hermitian, $\U$ is unitary, i.e. $\U^{\dag}\U =
\U\U^{\dag} = 1$.  This is important because it means time evolution
conserves probability.  Another property, which is easy to verify by
inspection, is that $\U(-t) = \U^{\dag}(t)$.

\section{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~\eq{propagator} 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 $\{\hat 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 $\H$ form a basis that is particularly easy to
work with.  Let $\ket|n>$ be an eigenfunction of the Hamiltonian with
the eigenvalue $E_n$.  Then $\ket|n>$ will also be an eigenfunction of
$\U$, propagating as
\begin{equation}
\ket|n;t> = \U(t) \ket|n> = e^{-i E_n t / \hbar} \ket|n>.
\end{equation}
The eigenfunctions of $\H$ 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 $\H$ for a general Hamiltonian is, of
course, much harder, but can be accomplished numerically.)  This
approach has already been used successfully, for example in
\cite{baranger93}.  I attempted to determine if satisfactory results could
be obtained using a different method.

For a basis other than the eigenfunctions of $\H$, treating the exponential
in equation~\eq{propagator} directly as a series sum is generally
impractical for numeric calculations, due to the computational costs of
applying $\H^n$.  For example, for decomposition into $\ket|\x>$, applying
$\H^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 $\H
= \K + \V$, where the kinetic energy operator is the usual $\K =
\hat\p^2 / 2m$, while the operator corresponding to potential energy
is given by $\V = V(\hat x_1,\ldots,\hat x_n)$, where $V(x_1,\ldots,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 $\{\hat p_i\}$, while $\V$ is easy on the
eigenfunctions of $\{\hat x_i\}$.  More specifically, $\K$ is diagonalized
in the momentum representation, while $\V$ is diagonalized in the position
representation:
\begin{xalignat}{2}
\K &= \sum_\p \ket|\p> {\p^2 \over 2m} \bra<\p|
&
\V &= \sum_\x \ket|\x> V(\x) \bra<\x|
\label{eq:diagonal}\end{xalignat}
Switching from one of these representations to the other
is covered in most quantum mechanics textbooks \cite{sakurai,gasiorowicz},
and involves a pair of Fourier transforms,
\begin{equation}\begin{aligned}
\braket<\p|\psi> &= \int \braket<\p|\x'> d^D\x' \braket<\x'|\psi>
   = (2\pi\hbar)^{-{D/2}}
	\int e^{-{i \p\cdot\x' / \hbar}} \braket<\x'|\psi> d^D\x' ,
\\
\braket<\x|\psi> &= \int \braket<\x|\p'> d^D\p' \braket<\p'|\psi>
   = (2\pi\hbar)^{-{D/2}}
	\int e^{+{i \p'\cdot\x / \hbar}} \braket<\p'|\psi> d^D\p' .
\end{aligned}\label{eq:fourier}\end{equation}

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 / \hbar}$ is not exactly separable into a product of
exponentials each containing only $\V$ and $\K$.  However, if $t$ is made
small enough compared to $\hbar/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 $\Delta t$, we can thus write
\begin{equation}
\U(\Delta t) \approx e^{-i\V\Delta t / \hbar} e^{-i\K\Delta t / \hbar}
             \approx e^{-i\K\Delta t / \hbar} e^{-i\V\Delta t / \hbar} .
\end{equation}
The two forms above can be combined into a single symmetric form, which is
beneficial for two reasons.  First, a symmetric form will preserve
$\U^{\dag}(t)=\U(-t)$.
Secondly, it effectively halves $\Delta t$ compared to the original
expressions.  The result\footnote{Of course, we could also symmetrize to
$e^{-i\K\Delta t / 2\hbar} e^{-i\V\Delta t / \hbar} e^{-i\K\Delta
t / 2\hbar}$, but the other form is more convenient if we want the
original and final state to be expressed in position representation.} is
\begin{equation}
\U(\Delta t) \approx \UU(\Delta t) \equiv
    e^{-i \V \Delta t / 2\hbar}
    e^{-i \K \Delta t / \hbar}
    e^{-i \V \Delta t / 2\hbar} .
\label{eq:separated}
\end{equation}

Arbitrarily long evolution can be accomplished by iterating application of
$\UU$ for acceptably short fragments of the required time interval.  Also
diagonalizing the parts of $\UU$ analogously to equation~\eq{diagonal} yields
\begin{multline}
\ket|\psi(N\Delta t)> \approx
  |\x_{i_1}\>     e^{-iV(\x_{i_1})             \Delta t/2\hbar} \<\x_{i_1}|
   \p_{j_1}\>     e^{-i{\p_{j_1}^2\over 2m}     \Delta t/\hbar} \<\p_{j_1}|
   \x_{i_2}\>     e^{-iV(\x_{i_2})              \Delta t/\hbar} \<\x_{i_2}|
	\p_{j_2}\> \ldots \\ \ldots \<\p_{j_{l-1}}|
   \x_{i_l}\>     e^{-iV(\x_{i_l})              \Delta t/\hbar} \<\x_{i_l}|
   \p_{j_l}\>     e^{-i{\p_{j_l}^2\over 2m}     \Delta t/\hbar} \<\p_{j_l}|
	\x_{j_{l+1}}\> \ldots \\ \ldots \<\x_{j_N}|
   \p_{j_N}\>     e^{-i{\p_{j_N}^2\over 2m}     \Delta t/\hbar} \<\p_{j_N}|
   \x_{i_{N+1}}\> e^{-iV(\x_{i_{N+1}})         \Delta t/2\hbar} \<\x_{i_{N+1}}|
\psi(t=0)\>
\label{eq:prop_mess}
\end{multline}
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^{-i \V \Delta t /2\hbar}$ have been joined everywhere else.

The occurrences of $\braket<\p_j|\x_i>$ and $\braket<\x_i|\p_j>$ in
equation~\eq{prop_mess}, under implied integration, perform a change of
basis, and can be identified with the Fourier transforms in
equation~\eq{fourier}.  The rest of equation~\eq{prop_mess} consists of
multiplications of the wave~function by a complex function.  Since both of
these operations are relatively simple to perform numerically,
equation~\eq{prop_mess} should provide a suitable origin for a
computational model of wave~function propagation.

\section{Discretization of space}

Before we continue, it is convenient to rewrite
equations~\eq{fourier}, 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/\hbar$, rather than the momentum:
\begin{xalignat}{2}
\Psi(\x) &= \braket<\x|\psi> &
\Phi(\k) &= \hbar^{D/2} \braket<\p=\hbar\k|\psi>
\end{xalignat}
(The factor of $\hbar^{D/2}$ is required for $\Phi(\k)$ to be
properly normalized.\footnote{If $\braket<\psi|\psi> = 1$, then
$\int\braket<\psi|\p>\braket<\p|\psi> d^D(\p/\hbar) = 1/\hbar^D$.})
Using these substitutions, equations~\eq{fourier} become
\begin{equation}\begin{aligned}
\Phi(\k) &= (2\pi)^{-{D/2}} \int e^{-{i \k\cdot\x'}} \Psi(\x') d^D\x' , \\
\Psi(\x) &= (2\pi)^{-{D/2}} \int e^{+{i \k'\cdot\x}} \Phi(\k') d^D\k' .
\end{aligned}\label{eq:fourfun}\end{equation}

We can now consider
which representation should be used for the wave~function in
equation~\eq{prop_mess}.  The simplest choice is to take the values of
the position-representation wave~function at a discrete rectangular
grid, i.e. to consider the set of values
\begin{equation}
\Psi_{l_1\ldots l_D} =
	\Psi(x_1=\delta_1l_1,\ldots,x_D=\delta_Dl_D) ,
\label{eq:sample,x}
\end{equation}
%taken over some range for each of the indices $l_1,\ldots,l_D$.
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
$\left(\delta_1\ldots\delta_D\right)^{-1/2}$ instead of $1$.

Let
$N_i$ be the number of distinct values taken by $l_i$.  In that case,
the momentum-representation wave~function can be approximated by the
discrete Fourier transform,
\begin{alignat}{1}
% kludge: just {align} moves the second equation number, thinking that
% the line is too long.  But {alignat}{1} is happy.
\label{eq:dft}
\tag{\theequation a}\label{eq:dft,p}
\Phi_{m_1\ldots m_D} &=
    {1\over\sqrt{N_1\ldots N_D}} \sum_{l_1}\cdots\sum_{l_D}
	e^{-2\pi i\left({l_1m_1\over N_1}+\cdots+{l_D m_D\over N_D}\right)}
	\Psi_{l_1\ldots l_D} , \\
\tag{\theequation b}\label{eq:dft,x}
\Psi_{l_1\ldots l_D} &=
    {1\over\sqrt{N_1\ldots N_D}} \sum_{m_1}\cdots\sum_{m_D}
	e^{+2\pi i\left({l_1m_1\over N_1}+\cdots+{l_D m_D\over N_D}\right)}
	\Phi_{m_1\ldots m_D} ,
\end{alignat}
\stepcounter{equation}
where the discrete wave~function is related to the continuous one by
\begin{equation}
\Phi_{m_1\ldots m_D} \approx
    \Phi(k_1={2\pi\over N_1\delta_1}m_1,\ldots,k_D={2\pi\over N_D\delta_D}m_D).
\label{eq:sample,p}
\end{equation}
The number of meaningful distinct values taken by $m_i$ will also be
$N_i$, since the indices are cyclic (that is, $\Phi_{\ldots m_i\ldots}
= \Phi_{\ldots(m_i\pm N_i)\ldots}$).\label{txt:cyclic}
Not very surprisingly, $\Phi_{m_1\ldots m_D}$ is also normalized to
$\left(\delta_1\ldots\delta_D\right)^{-1/2}$.


An algorithm known as the Fast Fourier Transform (FFT) provides an
efficient way to calculate the expressions~\eq{dft}, 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 \cite{oran,nussbaumer}.  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 \cite{fftw}.  I started using it after it
outperformed my own hand-optimized code, and I highly recommend it.
