\documentstyle[quiet,twoside,fullpage,12pt]{report}
\newcommand{\PSbox}[3]{\mbox{\special{psfile=/tmp/liju/#1}\hspace{#2}\rule{0pt}{#3}}}
\newcommand{\psA}[1]{\PSbox{#1 hoffset=-40 voffset=-105 hscale=45 vscale=50}{7.5cm}{7.2cm}}
\newcommand{\psB}[1]{\PSbox{#1 hoffset=-98 voffset=-200 hscale=100 vscale=100}{15cm}{14.4cm}}
\newcommand{\psC}[1]{\PSbox{#1 hoffset=-115 voffset=-330 hscale=105 vscale=160}{15cm}{21.4cm}}
\newcommand{\ee}{\begin{equation}}
\newcommand{\ed}{\end{equation}}
\newcommand{\mm}{\begin{eqnarray}}
\newcommand{\md}{\end{eqnarray}}
\newcommand{\nn}{\nonumber\\}
\newcommand{\av}[1]{< #1 >}
\newcommand{\pdf}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}}
\newcommand{\f}[2]{\frac{#1}{#2}}
\setlength{\parskip}{\baselineskip}
\setlength{\parindent}{0em}
\title{Quantum Scattering Problem on Computers\\
--- a term project for 22.113}
\date{\today}
\author{Li Ju}
\begin{document}
\maketitle
\tableofcontents
\newpage
\chapter{Introduction}
As a term project for MIT course 22.113 {\em Nuclear and Atomic
Collision Phenomena}, I investigated the computational approaches to
various kinds of quantum scattering problems. My motivation is first
to understand the analytic formulations on which these computational
methods are made possible, and then to implement them on the computer
and get correct results. Two characteristic problems were chosen:
1) The quantum scattering problem for arbitary central pair potential.
2) Time-dependent 1D and 2D Schrodinger Equation.
These problems, in general, are not very complicated. But we can learn
something from these nice ``little'' ones, and solve more complex
problems in the future.\\
Reference Books:
1) {\em Computational Physics}\\ by Steven E.Koonin, 1986 Addison-Wesley.
2) {\em Computational Methods in Physics \& Engineering}\\ by Samual S.M.Wong,
1992 Prentice-Hall.
3) {\em Numerical Recipes in C --- The Arts of Scientific Computing}.\\
2nd Edition, 1992 Cambridge University Press.
\chapter{Central Pair Scattering}
\section{Theoretical Formulation}
Denote the pairwise central potential as $V(r)$, the problem for
two-body quantum scattering can be routinely transformed into the one
body Schrodinger equation:
$$(\frac{-\hbar^{2}}{2\mu}\nabla^{2}+V(r))\psi(\vec{r}) = E\psi(\vec{r}) $$
where $\mu$ is the effective mass
\ee \mu = \frac{mM}{m+M} \ed
Let $k = \sqrt{2\mu E}/\hbar$. The asymptotic behaviour of the
scattered plane wave is
\ee \psi(\vec{r}) \longrightarrow e^{ikz}+f(\theta)\frac{e^{ikr}}{r} \ed
and the corresponding differential cross section is
\ee \sigma(\theta) = \frac{dN/d\Omega}{J_{in}} = |f(\theta)|^{2} \ed
Because $V(r)$ is pairwise, $[{\cal H ,J}] = 0$. The angular momentum
and the Hamiltonian can be diagonalized concurrently, so
\ee \psi(r,\theta) = \sum_{l=0}^{\infty}\frac{u_{l}(r)}{r}P_{l}(\cos \theta)
\ed
with $u_{l}(r)$ satisfies
\ee \frac{d^{2}u_{l}(r)}{dr^{2}}+[k^{2}-\frac{l(l+1)}{r^{2}}-\frac{2\mu V(r)}
{\hbar^{2}}]u_{l}(r) = 0 \ed
The general solution to this equation is of the form
\ee u_{l}(r) = B_{l}rj_{l}(kr)+C_{l}rn_{l}(kr) \ed
where $j_{l}(r)$ and $n_{l}(r)$ are the regular and irregular (Neumann
function) spherical Bessel functions. The asymptotic behaviour of
$u_{l}(r)$ at $kr\gg 1$ is $$ u_{l}(r) \longrightarrow
\frac{B_{l}}{k}\sin(
kr-\frac{l\pi}{2})-\frac{C_{l}}{k}\cos(kr-\frac{l\pi}{2}) $$
\ee = \frac{a_{l}}{k}\sin(kr-\frac{l\pi}{2}+\delta_{l}) \ed
so
\ee \psi(r,\theta) \longrightarrow \sum_{l=0}^{\infty}i^{l}(2l+1) a_{l} P_{l}(
\cos \theta) \frac{\sin(kr-l\pi/2+\delta_{l})}{kr} \ed
There's also an frequently used expansion of plane wave into spherical
harmonics: $$ \exp(ikr\cos \theta) = \sum_{l=0}^{\infty}
i^{l}(2l+1)j_{l}(kr)P_{l}(
\cos \theta) $$
\ee \longrightarrow \sum_{l=0}^{\infty} i^{l}(2l+1) P_{l}(
\cos \theta) \frac{\sin(kr-l\pi/2)}{kr} \ed
Compare with (2.2), we see that in order to achieve the predicted
scattering behaviour, there should be
\ee a_{l}\sin(kr-l\pi/2+\delta_{l})-\sin(kr-l\pi/2) = f_{l}e^{ikr} \ed
and $f_{l}$ should not be a function of $r$, so
$$a_{l}\frac{e^{ikr-il\pi /2+i\delta_{l}}-e^{-ikr+il\pi /2-i\delta_{l}}}
{2i}-\frac{e^{ikr-il\pi /2}-e^{-ikr+il\pi /2}}{2i} = f_{l}e^{ikr} $$
The $e^{-ikr}$ terms on the right hand side must vanish, so
\ee a_{l} = e^{i\delta_{l}} \ed
and we get
\ee f_{l} = \exp(\frac{-il\pi}{2})\frac{e^{2i\delta_{l}}-1}{2i}
= (-i)^{l}e^{i\delta_{l}}\sin\delta_{l} \ed
So the scattered wave is \ee f(\theta)\frac{e^{ikr}}{r} =
\sum_{l=0}^{\infty}(2l+1) e^{i\delta_{l}}\sin\delta_{l} P_{l}(
\cos \theta) \frac{\exp(ikr)}{kr} \ed
and the differential cross section is just
\ee \sigma(\theta) =
\frac{1}{k^{2}}|\sum_{l=0}^{\infty}(2l+1)
e^{i\delta_{l}}\sin\delta_{l} P_{l}( \cos \theta)|^{2}
\ed
and the total cross section is
\ee \sigma =
\frac{4\pi}{k^{2}}\sum_{l=0}^{\infty}(2l+1)
\sin^{2}\delta_{l}
\ed
Although the sums above extends over all $l$, we can select a
reasonable cutoff in the sense that when $$
\frac{l_{max}(l_{max}+1)\hbar^{2}}{2\mu r_{max}^{2}} \gg E $$ the
$l_{max}$ partial wave decrease so quickly even at the surface of the
potential field that it became negligible, which leads to the
approximation of
\ee l_{max} \sim kr_{max} \ed
This estimate is usually slightly low.
\newpage
\section{Computational Method}
To find the phase shift for each partial wave, we solve the
equation
\ee \frac{d^{2}u_{l}(r)}{dr^{2}}
+[k^{2}-\frac{l(l+1)}{r^{2}}-\frac{2\mu V(r)}
{\hbar^{2}}]u_{l}(r) = 0 \ed numerically using the {\bf Numerov Method},
which is a particular simple and efficient method in intergrating
2nd-order differential equations.
{\bf Numerov Method}: For equations with the general form
\ee \frac{d^{2}y}{dx^{2}}+k(x)y = S(x) \ed
the method gives
$$ (1+\frac{h^{2}}{12}k_{n+1})y_{n+1}-
2(1-\frac{5h^{2}}{12}k_{n})y_{n} +
(1+\frac{h^{2}}{12}k_{n-1})y_{n-1} $$
\ee =
\frac{h^{2}}{12}(S_{n+1}+10S_{n}+S_{n-1})+O(h^{6}) \ed
Notice the local error of $O(h^{6})$, which is one order more accurate
than the fourth-order Runge-Kutta method.
Since for our particular problem, $S(x)=0$, it's a linear differential
equation. And we know that $u_{l}(r=0) = 0$, so the choice of
$u_{l}(r=h)$ became totally arbitary, because we can always normalize
the entire partial wave later if neccesary. (However, $u_{l}(r=h)$ is
often chosen to be a small number to avoid overflow at greater $r$).
We intergrate the equation to a certain point $r_{1}$ where $V(r_{1})$
is small enough so that the wave function can be regarded as
approaching asymptotic behaviour. So
\ee u_{l}(r_{1}) \approx Akr_{1}
[\cos(\delta_{l})j_{l}(kr_{1})-\sin(\delta_{l})n_{l}
(kr_{1})] \ed
We continue intergrating to a larger radius $r_{2}>r_{1}$, it also has
\ee u_{l}(r_{2}) \approx Akr_{2}
[\cos(\delta_{l})j_{l}(kr_{2})-\sin(\delta_{l})n_{l} (kr_{2})] \ed
We can interpolate $\delta_{l}$ from these two points, after some
manipulations we get
$$ G =
\frac{r_{1}u_{l}(r_{2})}{r_{2}u_{l}(r_{1})}$$
$$\delta_{l} =
\tan^{-1}(\frac{Gj_{l}(r_{1})-j_{l}(r_{2})}
{Gn_{l}(r_{1})-n_{l}(r_{2})})$$
\newpage
\section{Results}
\psA{a2.ps} \psA{a1.ps}
\psA{a3.ps} \psA{a4.ps}
In this section we'll mainly show the computation results for hard
sphere potential. (Although arbitary central potential is easily
applicable, there don't exist much to compare with.) In classical
mechanics we've shown that the total cross section for hard sphere
collision is $\pi r^{2}$, with $r$ as the radius of the sphere. In
quantum mechanics, however, the result is quite different: For $kr\ll
1$, $\sigma_{t} \rightarrow 4\pi r^2$, and there will only be
significant contribution from $l=0$ so the scattering is isotropic,
i,e, the differntial cross section is a constant. This is clearly
shown in Fig.1.
When $kr\gg 1$, $\sigma_{t} \rightarrow 2\pi r^2$, and the
differential cross section is highly condensed in the front. The
computation result is shown in Fig.2.
The intermediate state of $kr = 1$ is shown in Fig.3.
Fig.4 is the $\sigma_{t}/\pi r^{2}$ versus $kr$ curve, it start at
$kr=0$ with the value of 4, and then come down to $2$ as $kr
\rightarrow \infty$. This curve is not directly available from
theory, but should be quite important. We can call it our first
score.
\begin{center}
\psA{a5.ps}
\end{center}
Fig.5 is the scattering by a Guassian potential:
\ee V(r) = \exp(-(kr)^{2}/\alpha^2) \ed
Here $\alpha$ is chosen to be 1, corresponding to situation of hard
sphere collison with $kr = 1$. One should note the unusual tail
inflextion at $\theta \simeq 1$.
\chapter{Time Dependent Schrodinger Equation}
The Time Dependent Schrodinger Equation (TDSE) is
\ee i\hbar\frac{\partial \psi}{\partial t}
= {\cal H} \psi = (\frac{-\hbar^{2}}{2\mu}\nabla^{2}+V(\vec{r}))\psi
\ed
Using reduced length and time unit we transform the equation into
\ee i\frac{\partial \psi}{\partial t}
= {\cal H} \psi = (-\nabla^{2}+V(\vec{r}))\psi
\ed
Such partial differential equation with 1st-order time derivative is
classified as {\em parabolic}, with respect to {\em elliptic}
equations (0th order) and {\em hyperbolic} equations (2nd order).
Another kind of parabolic differential equations is the diffusion
equation.
Under discretization, the 2nd order derivative becomes
\ee
\pdd{\phi(x)}{x} \longrightarrow \f{1}{h^2}\{\phi(x+h)-2\phi(x)+\phi(x-h) \}
+ {\cal O}(h^2) \ed
and the general solution for (3.2) can be written as
\ee
\psi(x,t) = \exp(-it{\cal H}) \psi(x,0)
\ed
The propagation operator can be approximated in the small time limit
as
\ee \exp(-i\tau{\cal H})\psi(x,0) = (1-i\tau{\cal H})\psi(x,0)
+{\cal O}(\tau^{2}) \ed The right hand side of the equation can be
directly evaluated using (3.3). It is called the {\em explicit} scheme
for the solution of parabolic equations. However, this scheme has
severe limitations, that is, when the time constant $\tau$ get large,
the scheme become unstable. This can be easily understood if we
consider the eigenfunctions of ${{\cal H}}$ under current boundary
conditions:
\ee {\cal H} \psi_{n} = \epsilon_{n} \psi_{n} \ed
Since ${\cal H}$ is discretized, the number of eigenmodes for (3.6)
isn't infinity as in continuum space, but equals to the $N$, the
degree of discretization in space domain. Let's consider the highest
eigenvalue (energy level), under rough approximation
\ee \epsilon_{N} \approx (\f{\pi}{h})^2 \ed
because $\pi/h$ is the edge of the first Brillouin zone, and we treat
it as a free particle. We can see that when $\tau\epsilon_{N}$ is in
the order of unity, the expansion (3.5) become invalid and those high
eigenvalued modes (actually noise) would be endlessly amplified until
the system breaks down. So the criterea of stabilty for the explicit
scheme is
\ee \tau(\f{\pi}{h})^2 \ll 1\ed
Since $h$ is usually chosen to be small in order to minimize the error
in (3.3), we had to pick very small $\tau$ to stablize the scheme,
thus making the observation of a time evolution process very
difficult.
There is a much better scheme called the {\bf Crank-Nicholson method}:
we can approximate the propagation operator of (3.4) in the {\bf Cayley
form} as
\ee
\exp(-i\tau{\cal H}) = \f{1-\f{1}{2}i\tau{\cal H}}{1+\f{1}{2}i\tau\cal
H} + {\cal O}(\tau^3) \ed We observe that this form is one order more
accurate than the explicit scheme. What's more, it's an unitary
operator and therefor the amplitude of each eigenvalue component will
always remain the same as the initial amplitude, because
\ee \left | \f{1-\f{1}{2}i\tau\epsilon_{n}}{1+\f{1}{2}i\tau\epsilon_{n}}
\right | = 1 \ed
This property is called {\em norm conserving} and it is the
fundamental attribute of the Schrodinger wave equation. One direct
conclusion is that the scheme is always {\bf stable}, even when
$\tau\epsilon_{N} \sim 1$, provided that the initial amplitude of the
high $k$ components are small, they will still cause little error
later. What count are those components in the low $k$ regime with
large amplitudes, where $\tau\epsilon$ is small and the expansion is
valid.
However, this is an {\em implicit} scheme because the Cayley operator
can't be directly evaluated. We can rewrite (3.9) as
\ee
\{ 1+\f{1}{2}i\tau{\cal H} \} \psi(x_{j},t+\tau)
= \{ 1-\f{1}{2}i\tau{\cal H} \} \psi(x_{j},t) \ed
For one dimensional case
\ee {\cal H} = -\frac{d^2}{dx^2} + V(x) \ed
so
\ee \{ 1+\f{i\tau}{2}V(x_{j})-\f{i\tau}{2}\frac{d^2}{dx^2} \}
\psi(x_{j},t+\tau)
= \{ 1-\f{i\tau}{2}V(x_{j})+\f{i\tau}{2}\frac{d^2}{dx^2} \}
\psi(x_{j},t) \ed
After some rearrangements we get the following equation
\ee
A_{i}^{-}\psi_{i-1}^{\prime}+A_{i}^{0}\psi_{i}^{\prime}
+A_{i}^{+}\psi_{i+1}^{\prime} = b_{i}
\ed
with
\ee A_{i}^{-} = A_{i}^{+} = \f{-i\tau}{2h^{2}} \ed
\ee A_{i}^{0} = 1 + \f{i\tau}{2}V_{i}+\f{i\tau}{h^{2}} \ed
\ee b_{i} = (2-A_{i}^{0})\psi_{i}-A_{i}^{+}\psi_{i+1}-A_{i}^{-}\psi_{i-1} \ed
This forms a set of $N\times N$ linear equations, which usually
requires $N^{2}$ operations even in applying the inversion matrix each
time. Fortunately, the following algorithm (Gaussian elimination and
back-substitution) provides a very efficient solution (of the order
$N$ operations) for the tri-diagonal system as (3.14).
Using {\bf C} language convention, let's represent the space in an
$N+1$ matrix: the Dirichlet boundary condition is imposed on
$\psi^{\prime}_{0}$ and $\psi^{\prime}_{N}$, which might be time
dependent. We assume that solution of (3.14) satisfies relations of
the form
\ee \psi^{\prime}_{i+1} = \alpha_{i}\psi^{\prime}_{i}+\beta_{i} \ed
and substitute it into (3.14)
\ee
A_{i}^{-}\psi^{\prime}_{i-1}+A_{i}^{0}\psi^{\prime}_{i}
+A_{i}^{+}(\alpha_{i}\psi^{\prime}_{i}+\beta_{i}) = b_{i} \ed This is
of the same form as (3.18) except one regresion of index, so comparing
with the assumed relation for $\psi_{i-1}^{\prime}$, we identify the
backward recursion relationship for $\alpha_{i-1}$ and $\beta_{i-1}$
to be
\mm \gamma_{i} && = -1/(A_{i}^{0}+A_{i}^{+}\alpha_{i}) \nn
\alpha_{i-1} && = \gamma_{i}A_{i}^{-} \nn
\beta_{i-1} && = \gamma_{i}(A_{i}^{+}\beta_{i}-b_{i})
\md
In order to satisfy the boundary condition at $\psi^{\prime}_{N}$, we
let
\mm \alpha_{N-1} =0& \beta_{N-1} = \psi^{\prime}_{N} \md
and make a backsweep from here to get $\alpha_{i},\beta_{i}
(i=N-2,0)$. After that we do a forward sweep from $i=0$ to $N-1$ to
get the new wave function.
\section{One Dimensional Scattering of a Wave Packet}
There are several possibilities in selecting the initial form of the
wave function. The recommendation of Goldberg, Schey, and Schwartz is
to use a Guassian wave packet:
\ee \Psi(x,t=0) = \exp(ik_{0}x-(x-x_{c})^{2}/2\sigma_{0}^{2}) \ed
The wave packet has a group velocity of $k_{0}$, but it will also
spread out by itself(dispersion) as time evolves. We let it collide
with an square barrier potential to see what happens.
The total length of the scenario was chosen to be $1$, and was divided
into 1000 equal parts, so $h=0.001$.$\sigma_{0} = 40$ and $k_{0}$ was
chosen to be $50\pi$, and $\tau$ was chosen to be $2\times 10^{-6}$,
so the expansion parameter in the spectrum center is about
$k_{0}^{2}\tau
\sim 0.05$, which is marginally satisfactory since the Cayley form
perserves 2nd order accuracy. The height of the potential barrier was
in the order $10^4$ to match the kinetic energy, and the width to be
of 0.05, thus about 50 mesh sites.
Shown in Fig.1 is the case of $V=50000$, which is about twice the
energy of the incoming wave --- so it's mostly bounced back when
hitting the barrier.
Fig.2 is the case where the kinetic energy is one time higher than the
potential barrier, so it passed through quite easily.
In Fig.3 $V=50000$, which is the same as the first case. But instead
of Gaussian form, the initial wave pack is a square. Since there are
usually more high-frequency components contained in a discontinous
function, there are more wave components passing through than the
first case.
\newpage \psC{pic1.ps} \newpage \psC{pic2.ps} \newpage \psC{pic3.ps}
\section{Two dimensional TDSE}
There are complications in the two dimensional TDSE: we can no longer
solve the implicit equations using the efficient back-substitution
method for tri-diagonal systems as in the one dimensional case,
because the discretization
\ee \nabla^{2}\psi_{i,j} = \frac{1}{h^2}(\psi_{i+1,j}+\psi_{i-1,j}
+\psi_{i,j+1}+\psi_{i,j-1}-4\psi_{i,j}) + {\cal O}(h^2) \ed
involves more than 3 variables. However, we can always seperate ${\cal
H}$ into ${\cal H}_{x}+{\cal H}_{y}$, and apply each Cayley operator
onto the system seperately. We'll analyze the accuracy of such
factorization down below:
Since
\ee {\cal H} = -\frac{d^2}{dx^2}-\frac{d^2}{dy^2}+V(x,y) \ed
Let's define \mm {\cal H}_{x} &&= -\frac{d^2}{dx^2} + V_{1}(x,y) \nn
{\cal H}_{y} &&= -\frac{d^2}{dy^2} + V_{2}(x,y) \nn
V(x,y) &&= V_{1}(x,y)+V_{2}(x,y) \md
so \ee {\cal H} = {\cal H}_{x}+{\cal H}_{y} \ed
and
\mm
\exp(-i\tau{\cal H}) &&= \exp(-i\tau{\cal H}_{x}-i\tau{\cal H}_{y}) \nn
&&=1-i\tau{\cal H}_{x}-i\tau{\cal H}_{y}+\f{1}{2}
(i\tau{\cal H}_{x}+i\tau{\cal H}_{y})^2+{\cal O}(\tau^3) \nn
\md
Because
\mm
\f{1-\f{1}{2}i\tau{\cal H}_{x}}{1+\f{1}{2}i\tau{\cal
H}_{x}} &&=
(1-\f{1}{2}i\tau{\cal H}_{x})[1-\f{1}{2}i\tau{\cal
H}_{x}+(\f{1}{2}i\tau{\cal
H}_{x})^2]+{\cal O}(\tau^3) \nn
&&= 1-i\tau{\cal H}_{x}+\f{1}{2}(i\tau{\cal H}_{x})^2+{\cal O}(\tau^3)
\md
so
\mm
(\f{1-\f{1}{2}i\tau{\cal H}_{x}}{1+\f{1}{2}i\tau{\cal
H}_{x}})
(\f{1-\f{1}{2}i\tau{\cal H}_{y}}{1+\f{1}{2}i\tau{\cal
H}_{y}}) =
[1-i\tau{\cal H}_{x}+\f{1}{2}(i\tau{\cal H}_{x})^2]
[1-i\tau{\cal H}_{y}+\f{1}{2}(i\tau{\cal H}_{y})^2]
+{\cal O}(\tau^3) \nn
=1-i\tau{\cal H}_{x}-i\tau{\cal H}_{y}
+\f{1}{2}(i\tau{\cal H}_{x})^2
+\f{1}{2}(i\tau{\cal H}_{y})^2
-\tau^2 {\cal H}_{x}{\cal H}_{y}
+{\cal O}(\tau^3)
\md
Comparing with (3.27)
\ee \exp(-i\tau{\cal H}) =
(\f{1-\f{1}{2}i\tau{\cal H}_{x}}{1+\f{1}{2}i\tau{\cal
H}_{x}})
(\f{1-\f{1}{2}i\tau{\cal H}_{y}}{1+\f{1}{2}i\tau{\cal
H}_{y}}) + \f{\tau^2}{2}[{\cal H}_{x},{\cal H}_{y}] +
{\cal O}(\tau^{3}) \ed
We can see that the error caused by this two-step process is in the
order of $\tau^2$, thus destroying one order of accuracy.
Nevertheless, we'll adopt this method because otherwise the cost would
become unbearable.
One may intuitively suggest picking $V_{1}(x,y)$ and $V_{2}(x,y)$ such
that $[{\cal H}_{x},{\cal H}_{y}] = 0$. Indeed, that's the purpose of
the author at the first time to split $V(x,y)$ into two parts.
However, after some manipulations it's found to be generally
impossible, because the commutation of $d^2/dx^2$ with a scaler
function would generate free $d/dx$ operators, which can't be offset
by simple manipulations. Another idea is to symmetrize the right hand
side of (3.30) and use the fact that
\ee \exp(-i\tau{\cal H}) =
(\f{1-\f{1}{2}i\tau{\cal H}_{x}}{1+\f{1}{2}i\tau{\cal H}_{x}})
(\f{1-\f{1}{2}i\tau{\cal H}_{y}}{1+\f{1}{2}i\tau{\cal H}_{y}})
+(\f{1-\f{1}{2}i\tau{\cal H}_{y}}{1+\f{1}{2}i\tau{\cal H}_{y}})
(\f{1-\f{1}{2}i\tau{\cal H}_{x}}{1+\f{1}{2}i\tau{\cal H}_{x}})+ {\cal
O}(\tau^{3}) \ed
However, that destroy the {\em norm conserving}
property of the simple scheme (3.30), because the sum of two unitary operators
isn't necessarily an unitary operator.
Next we'll analyze the influence of finite discretization in both
space and time to our computation. One may immediately remember the
dispersion relation derived by Huang in our first solid state physics
course, and how the finite discretization of space (lattice sites)
influence the behaviour of a should-be-continous wave. In general,
suppose a wave component
\ee \psi = \exp(ikx-i\omega t) \ed
in free propagation (without potential). Under finite discretization
of space:
\mm {\cal H}\psi &&= \f{-1}{h^2}[\exp(ik(x+h)-i\omega t)+\exp(ik(x-h)-i\omega t)-2\exp(ikx-i\omega t)] \nn
&&= \f{-1}{h^2}(\exp(ikh)+\exp(-ikh)-2) \psi \nn
&&= \f{2-2\cos(kh)}{h^2} \psi \nn
&&= \epsilon\psi
\md
and under further discretization of time by putting it into the Cayley
form (recursion relation), we get the following dispersion relation:
\ee
\exp(-i\omega\tau) = \f{1-\f{1}{2}i\tau\epsilon}{1+\f{1}{2}i\tau\epsilon}
\ed
or
\mm \omega(k) &&= \f{2}{\tau} \tan^{-1}(\f{\tau-\tau\cos(kh)}{h^2}) \nn
&&= k^2 - \f{1}{12}h^2k^4+(\f{1}{360}h^4-\f{1}{12}\tau^2)k^6+{\cal O}(k^8)
\md
Obviously in the long wave length limit we have our
\ee w(k) \simeq k^2 \ed
behaviour back. The question is how long is long, and that is answered
by the other terms in (3.34).
As an example of 2D TDSE implementation, the following scenario was
devised: A $1.0\times 1.0$ square was devided into a $500\times 500$
mesh array. At $t=0$ the initial condition was a steady plane wave flow
in the $X$ direction
\ee \psi(x,y) = \exp(ikx-i\omega t) \ed
$k = 10\pi$ so the wavelength $\lambda = 1/5 = 100$ mesh points, which
was considered to be enough for a sinusoidal function. $\tau$ was
chosen to be $2\times 10^{-6}$ and $\omega$ can be calculated using
dispersion relation (3.34). The boundary condition of (3.36) was
imposed on 4 sides of the square, to simulate the enviroment of an
infinite ``river'' of steady flow. (However, this would cause the
biggest trouble later on.) Such situation would continue for $300\tau$
to test the code, if it's working, then the plane wave would continue
to be a plane wave and $|\psi(x,y,t<300\tau)| = 1$. At $t=300\tau$
something happens: a cylinder barrier of radius = 50 mesh points,
height = 300000 suddenly jump out of nothing in the center, and stirs
the flow.
\psA{surf11.ps} \psA{contour11.ps}
\newpage
\begin{center}
\psA{surf16.ps} \psA{contour16.ps} \\
.\\
.\\
.\\
\psA{surf21.ps} \psA{contour21.ps} \\
\end{center}
\newpage
\begin{center}
\psA{surf26.ps} \psA{contour26.ps} \\
.\\
.\\
.\\
\psA{surf31.ps} \psA{contour31.ps} \\
\end{center}
\newpage
\begin{center}
\psA{surf36.ps} \psA{contour36.ps} \\
.\\
.\\
.\\
\psA{surf41.ps} \psA{contour41.ps} \\
\end{center}
\end{document}