\documentstyle[a4wide,12pt]{article}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\OP}[1]{{\bf\widehat{#1}}}

\begin{document}
\section*{Introduction to numerical projects}

Here follows a brief recipe and recommendation on how to write a report for each
project.
\begin{itemize}
\item Give a short description of the nature of the problem and the eventual 
numerical methods you have used.
\item Describe the algorithm you have used and/or developed. Here you may find it convenient
to use pseudocoding. In many cases you can describe the algorithm
in the program itself.

\item Include the source code of your program. Comment your program properly.
\item If possible, try to find analytic solutions, or known limits
in order to test your program when developing the code.
\item Include your results either in figure form or in a table. Remember to
       label your results. All tables and figures should have relevant captions
       and labels on the axes.
\item Try to evaluate the reliabilty and numerical stability/precision
of your results. If possible, include a qualitative and/or quantitative
discussion of the numerical stability, eventual loss of precision etc. 

\item Try to give an interpretation of you results in your answers to 
the problems.
\item Critique: if possible include your comments and reflections about the 
exercise, whether you felt you learnt something, ideas for improvements and 
other thoughts you've made when solving the exercise.
We wish to keep this course at the interactive level and your comments can help
us improve it.
\item Try to establish a practice where you log your work at the 
computerlab. You may find such a logbook very handy at later stages
in your work, especially when you don't properly remember 
what a previous test version 
of your program did. Here you could also record 
the time spent on solving the exercise, various algorithms you may have tested
or other topics which you feel worthy of mentioning.
\end{itemize}



\section*{Format for electronic delivery of report and programs}
%
The preferred format for the report is a PDF file. You can also
use DOC or postscript formats. 
As programming language we prefer that you choose between C/C++ and Fortran90/95.
You could also use Java or Python as programming languages. 
Matlab/Maple/Mathematica/IDL are not accepted, but you can use them to check your results where possible.
Finally, 
we do prefer that you work together. Optimal working groups consist of 
2-3 students, but more people can collaborate. You can then hand in a common report. 



\section*{Project 2, Variational and Diffusion Monte Carlo for BEC with $^{87}$Rb, deadline  june 6  12pm (midnight)}
The spectacular demonstration of Bose-Einstein condensation (BEC) in gases of
alkali atoms $^{87}$Rb, $^{23}$Na, $^7$Li confined in magnetic
traps \cite{anderson95,davis95,bradley95} has led to an explosion of interest in
confined Bose systems. Of interest is the fraction of condensed atoms, the
nature of the condensate, the excitations above the condensate, the atomic
density in the trap as a function of Temperature and the critical temperature of BEC,
$T_c$. The extensive progress made up to early 1999 is reviewed by Dalfovo et
al.\cite{dalfovo99}.

A key feature of the trapped alkali and atomic hydrogen systems is that they are
dilute. The characteristic dimensions of a typical trap for $^{87}$Rb is
$a_{h0}=\left( {\hbar}/{m\omega_\perp}\right)^\frac{1}{2}=1-2 \times 10^4$
\AA\ (Ref.~\cite{anderson95}). The interaction between $^{87}$Rb atoms can be well represented
by its s-wave scattering length, $a_{Rb}$. This scattering length lies in the
range $85 < a_{Rb} < 140 a_0$ where $a_0 = 0.5292$ \AA\ is the Bohr radius.
The definite value $a_{Rb} = 100 a_0$ is usually selected and
for calculations the definite ratio of atom size to trap size 
$a_{Rb}/a_{h0} = 4.33 \times 10^{-3}$ 
is usually chosen \cite{dalfovo99}. A typical $^{87}$Rb atom
density in the trap is $n \simeq 10^{12}- 10^{14}$ atoms/cm$^3$ giving an
inter-atom spacing $\ell \simeq 10^4$ \AA. Thus the effective atom size is small
compared to both the trap size and the inter-atom spacing, the condition
for diluteness (i.e., $na^3_{Rb} \simeq 10^{-6}$ where $n = N/V$ is the number
density). In this limit,
although the interaction is important, dilute gas approximations such as the
Bogoliubov theory \cite{bogoliubov47}, valid for small $na^3$ and large
condensate fraction $n_0 = N_0/N$, describe the system well. Also, since most
of the atoms are in the condensate (except near $T_c$), the Gross-Pitaevskii
equation \cite{gross61,pitaevskii61} for the condensate describes the whole gas
well. Effects of atoms excited above the condensate have been incorporated
within the Popov approximation \cite{hutchinson97}. 

The aim of this project is to use Variational Monte
Carlo (VMC) and Diffusion Monte Carlo (DMC) methods and evaluate 
the ground state energy of
a trapped, hard sphere Bose gas for different numbers of particles
with a specific
trial wave function. See Ref.~\cite{abinitio} for a discussion of VMC and DMC. 

This wave function is used 
to study the sensitivity of condensate and 
non-condensate properties to the hard sphere radius and the number 
of particles.
The trap we will use is  a spherical (S) or 
an elliptical (E) harmonic trap in three dimensions given by 
 \begin{equation}
V_{ext}({\bf r}) = 
\Bigg\{
\begin{array}{ll}
        \frac{1}{2}m\omega_{ho}^2r^2 & (S)\\
\strut
	\frac{1}{2}m[\omega_{ho}^2(x^2+y^2) + \omega_z^2z^2] & (E)
\label{trap_eqn}
\end{array}
\end{equation}
where (S) stands for symmetric and 
\begin{equation}
    H = \sum_i^N \left(
        \frac{-\hbar^2}{2m}
        { \bigtriangledown }_{i}^2 +
        V_{ext}({\bf{r}}_i)\right)  +
        \sum_{i<j}^{N} V_{int}({\bf{r}}_i,{\bf{r}}_j),
\end{equation}
as the two-body Hamiltonian of the system.
Here $\omega_{ho}^2$ defines the trap potential strength.  In the case of the
elliptical trap, $V_{ext}(x,y,z)$, $\omega_{ho}=\omega_{\perp}$ is the trap frequency
in the perpendicular or $xy$ plane and $\omega_z$ the frequency in the $z$
direction.
The mean square vibrational amplitude of a single boson at $T=0K$ in the 
trap (\ref{trap_eqn}) is $<x^2>=(\hbar/2m\omega_{ho})$ so that 
$a_{ho} \equiv (\hbar/m\omega_{ho})^{\frac{1}{2}}$ defines the 
characteristic length
of the trap.  The ratio of the frequencies is denoted 
$\lambda=\omega_z/\omega_{\perp}$ leading to a ratio of the
trap lengths
$(a_{\perp}/a_z)=(\omega_z/\omega_{\perp})^{\frac{1}{2}} = \sqrt{\lambda}$.

We represent the inter boson interaction by a pairwise, hard core potential
\begin{equation}
V_{int}(|{\bf r}_i-{\bf r}_j|) =  \Bigg\{
\begin{array}{ll}
        \infty & {|{\bf r}_i-{\bf r}_j|} \leq {a}\\
        0 & {|{\bf r}_i-{\bf r}_j|} > {a}
\end{array}
\end{equation}
where ${a}$ is the hard core diameter of the bosons.  Clearly, $V_{int}(|{\bf r}_i-{\bf r}_j|)$
is zero if the bosons are separated by a distance $|{\bf r}_i-{\bf r}_j|$ greater than $a$ but
infinite if they attempt to come within a distance $|{\bf r}_i-{\bf r}_j| \leq a$.

Our trial wave function for the ground state with $N$ atoms is given by
\be
\Psi_T({\bf R})=\Psi_T({\bf r}_1, {\bf r}_2, \dots {\bf r}_N,\alpha,\beta)=\prod_i g(\alpha,\beta,{\bf r}_i)\prod_{i<j}f(a,|{\bf r}_i-{\bf r}_j|),
\label{eq:trialwf}
\ee
where $\alpha$ and $\beta$ are variational parameters. The single-particle wave function is proportional
to the harmonic oscillator function for the ground state, i.e.,
\be
   g(\alpha,\beta,{\bf r}_i)= \exp{[-\alpha(x_i^2+y_i^2+\beta z_i^2)]}.
\ee
For spherical traps we have $\beta = 1$ and for non-interacting bosons ($a=0$) we have
$\alpha = 1/2a_{ho}^2$.
The correlation wave function is 
\be
   f(a,|{\bf r}_i-{\bf r}_j|)=\Bigg\{
\begin{array}{ll}
        0 & {|{\bf r}_i-{\bf r}_j|} \leq {a}\\
        (1-\frac{a}{|{\bf r}_i-{\bf r}_j|}) & {|{\bf r}_i-{\bf r}_j|} > {a}.
\end{array}
\ee  
\section*{Exercise 1: Variational Monte Carlo study of the BEC ground state}

\begin{enumerate}
\item[1a)] Find analytic expressions for the local energy 
\be
   E_L({\bf R})=\frac{1}{\Psi_T({\bf R})}H\Psi_T({\bf R}),
   \label{eq:locale}
\ee
for the above 
trial wave function of Eq.~(\ref{eq:trialwf}).
Compute also the analytic expression for the drift force to be used in the diffusion part given by
\be
  F = \frac{2\nabla \Psi_T}{\Psi_T}.
\ee
\item[1b)] Write a Variational Monte Carlo program which uses standard Metropolis sampling 
and compute the ground state energy 
of a spherical harmonic oscillator ($\beta = 1$) with no interaction.     
Use natural units and make an analysis of your calculations using both the analytic expression for the 
local energy and a numerical calculation of the kinetic energy using numerical derivatives.
The only variational parameter is $\alpha$. Perform these calculations for $N=10, 50 $ and
$100$ atoms. Compare your results with the exact answer. 
\item[1c)] We turn now to the elliptic trap with a hard core interaction. 
We fix, as in Ref.~\cite{dubois2001} $a/a_{ho}=0.0043$. Introduce lengths in units 
of $a_{ho}$, $r\rightarrow r/a_{ho}$ and energy in units of $\hbar\omega_{ho}$.
Show then that the original Hamiltonian can be rewritten as 
\be 
   H=\sum_{i=1}^N\frac{1}{2}\left(-\nabla^2_i+x_i^2+y_i^2+\gamma^2z_i^2\right)+\sum_{i<j}V_{int}(|{\bf r}_i-{\bf r}_j|).
\ee
What is the expression for $\gamma$?
Choose the initial value for $\beta=\gamma = 2.82843$ and set up a VMC program
which computes the ground state energy using the trial wave function of Eq.~(\ref{eq:trialwf}). 
using only $\alpha$ as variational parameter.
Use standard Metropolis sampling and vary the parameter $\alpha$ in order to find a 
minimum. Perform the calculations for $N=10,50$ and $N=100$ and compare your results to those from the ideal case in the previous exercise. 
\item[1d)] We repeat exercise 1c), but now we replace the brute force Metropolis algorithm with 
importance sampling based on the Fokker-Planck and the Langevin equations. You will need this part for the next exercise as well.
If the parallel cluster setup works (hopefully!), you should also parallelize your code in much the same
fashion as you did with project 1.
\end{enumerate}
\section*{Exercise 2: Diffusion Monte Carlo study of the BEC ground state}
The DMC method is based on rewriting the 
Schr\"odinger equation in imaginary time, by defining
$\tau=it$. The imaginary time Schr\"odinger equation is then
\be
   \frac{\partial \psi}{\partial \tau}=-\OP{H}\psi,
\ee
where we have omitted the dependence on $\tau$ and the spatial variables
in $\psi$.
The wave function $\psi$  is again expanded in eigenstates of the Hamiltonian 
\be
    \psi = \sum_i^{\infty}c_i\phi_i,
    \label{eq:wexpansion_dmc}
\ee 
where 
\be 
   \OP{H}\phi_i=\epsilon_i\phi_i, 
\ee 
$\epsilon_i$ being an eigenstate of $\OP{H}$. 
A formal solution of the imaginary time Schr\"odinger equation is  
\be 
   \psi(\tau_1+\delta\tau)=e^{-\OP{H}\delta\tau}\psi(\tau_1) 
\ee 
where the state  $\psi(\tau_1)$ 
evolves from an imaginary time $\tau_1$ to a later time $\tau_1+\delta$.  
If the initial state $\psi(\tau_1)$ is 
expanded in energy ordered eigenstates,
following Eq.~(\ref{eq:wexpansion_dmc}), then we obtain
\be
    \psi(\delta\tau)=\sum_i^{\infty}c_ie^{-\epsilon_i\delta\tau}\phi_i.
\ee

Hence any initial state, $\psi$, 
that is not orthogonal to the ground state $\phi_0$ 
will evolve to the ground state in the long time limit, that is
\be
  \lim_{\tau\rightarrow\infty}\psi(\delta\tau)=c_0e^{-\epsilon_0\tau}\phi_0.
\ee
This derivation shares many formal similarities with that given 
for the variational principle discussed in the previous
sections. However in the DMC method the imaginary
time evolution results in excited states decaying exponentially fast, 
whereas in the VMC method any excited state contributions 
remain and contribute
to the VMC energy. 

The DMC method is a realisation of the above derivation in position space.
Including the spatial variables as well, the above equation reads
\be
  \lim_{\tau\rightarrow\infty}\psi({\bf R}, \delta\tau)=
   c_0e^{-\epsilon_0\tau}\phi_0({\bf R}).
   \label{eq:wexpansion_dmc2}
\ee

By introducing a constant offset to the energy, 
$E_T=\epsilon_0$, 
the long-time limit of Eq.~(\ref{eq:wexpansion_dmc2}) 
can be kept finite. If the Hamiltonian is separated into the
kinetic energy and potential terms, the imaginary time Schr\"odinger equation,
takes on a form similar to a diffusion equation, namely
\be
   -\frac{\partial \psi({\bf R}, \tau)}{\partial \tau}=
    \left[\sum_i^{N}-\frac{1}{2}\nabla^2_i\psi({\bf R}, \tau)\right]
    +(V({\bf R})-E_T)\psi({\bf R}, \tau).
    \label{eq:dmcequation1}
\ee
This equation is a diffusion equation where  the wave function $\psi$
may be interpreted as the density of diffusing particles (or ``walkers''), 
and the term  $V({\bf R})-E_T$ 
is a rate term describing a potential-dependent increase 
or decrease in the particle density. 
The above equation may be transformed into a form suitable for Monte Carlo methods, but this leads to a very inefficient algorithm. The
potential $V({\bf R})$  
is unbounded in coulombic systems and hence the rate term
$V({\bf R})-E_T$  can diverge. 
Large fluctuations in the particle density then
result and give impractically large statistical errors. 

These fluctuations may be substantially reduced by the incorporation 
of importance sampling in the
algorithm. 
Importance sampling is essential for DMC methods, 
if the simulation is to be efficient. 
A trial or guiding wave function $\psi_T({\bf R})$, which closely
approximates the ground state wave function is introduced.
This is where typically the VMC result would enter, see also
discussion below
A new distribution is defined as 
\be
   f({\bf R}, \tau)=\psi_T({\bf R})\psi({\bf R}, \tau),
\ee
which is also a solution of the Schr\"odinger equation when  
$\psi({\bf R}, \tau)$ 
is a solution. 
Eq.~(\ref{eq:dmcequation1}) consequently modified to
\be
   \frac{\partial f({\bf R}, \tau)}{\partial \tau}=
    \frac{1}{2}\nabla\left[\nabla -F({\bf R})\right]f({\bf R}, \tau)
    +(E_L({\bf R})-E_T)f({\bf R}, \tau).
    \label{eq:dmcequation2}
\ee
In this equation we have introduced the so-called force-term $F$,
given by
\be
   F({\bf R})=\frac{2\nabla \psi_T({\bf R})}{ \psi_T({\bf R})},
\ee
and is commonly referred to as the ``quantum force''. 
The local energy $E_L$ is defined as previously
\be
    E_L{\bf R})=-\frac{1}{\psi_T({\bf R})}
                \frac{\nabla^2 \psi_T({\bf R})}{2}+V({\bf R})\psi_T({\bf R}),
\ee
and is computed, as in the VMC method, with respect to the trial wave function.

We can give the following interpretation to Eq.~(\ref{eq:dmcequation2}).
The right hand side of the importance sampled DMC equation
consists, from left to right, of diffusion, drift and rate terms. The
problematic potential dependent rate term of the non-importance 
sampled method is replaced by a term dependent on the difference 
between the local
energy of the guiding wave function and the trial energy. 
The trial energy is initially chosen to be the VMC energy of 
the trial  wave function, and is
updated as the simulation progresses. Use of an optimised 
trial function minimises the difference between the local 
and trial energies, and hence
minimises fluctuations in the distribution $f$ . 
A wave function optimised using VMC is ideal for this purpose, 
and in practice VMC provides the best
method for obtaining wave functions that accurately 
approximate ground state wave functions locally. 
The trial wave function may be also
constructed to minimise the number of divergences in , 
unlike the non-importance sampled method where divergences in the coulomb
interactions are always present. 

\begin{enumerate}

\item[2a)] The last part of this project deals with the implementation of the DMC algorithm 
using the short time approximation and importance sampling, see section 3.2 of Ref.~\cite{abinitio}.

Here you need as input your results from the VMC calculation for $\Psi_T$ and $E_T$. 
Perform these calculations for the same number of atoms as in Exercise 1.
Since we are only interested in the ground state, we can restrict the number of walkers to at most 
ten. Study the satbility of the results as a function of the number of walkers.

Compare also your results with the VMC results from the previous exercise. 
    


\end{enumerate}

\section*{Literature}
\begin{thebibliography}{999}%\parskip=0pt\itemsep 0pt%
\footnotesize
\bibitem {anderson95}%1
M.H. Anderson, J.R. Ensher, M.R. Matthews, C.E. Wieman, and E.A. Cornell, {\em Science} {\bf
269}, 198 (1995).

\bibitem{davis95}%2
K.B. Davis, M.-O. Mewes, M.R. Andrews, N.J. van Druten, D.S. Durfee, D.M. Kurn, and W. Ketterle,
{\em Phys. Rev. Lett.} {\bf 75}, 3969 (1995).

\bibitem{bradley95}%3
C.C. Bradley, C.A. Sackett, J.J. Tolett, and R.G. Hulet, {\em Phys. Rev. Lett.} {\bf 75}, 1687
(1995); C.C. Bradley, C.A. Sackett, and R.G. Hulet, {\em ibid.} {\bf 78}, 985 (1997).

\bibitem{dalfovo99}%4
F. Dalfovo, S. Giorgini, L. Pitaevskii, and S. Stringari, {\em Rev. Mod. Phys.} {\bf 71}, 463
(1999).

\bibitem{bogoliubov47}%5
N.N. Bogoliubov, {\em J. Phys. (Moscow)} {\bf 11}, 23 (1947).

\bibitem{gross61}%6
E.P. Gross, {\em Nuovo Cimento} {\bf 20}, 454 (1961).

\bibitem{pitaevskii61}%7
L.P. Pitaevskii, {\em Ah. Eksp. Teor. Fiz.} {\bf 40}, 646 (1961) [{\em Sov. Phys. JETP} {\bf 13},
451 (1961)].

\bibitem{hutchinson97}%8
D.A.W. Hutchinson, E. Zaremba, and A. Griffin, {\em Phys. Rev. Lett.} {\bf 78}, 1842 (1997).



\bibitem{abinitio} B.~L.~Hammond, W.~A.~Lester and P.~J.~Reynolds, Monte Carlo methods
in Ab Inition Quantum Chemistry, (World Scientific, Singapore, 1994), chapters 2 and 3.

\bibitem{dubois2001} J.~L.~DuBois and H.~R.~Glyde, Phys.~Rev.~{\bf A63}, 023602 (2001). 
\end{thebibliography}



\end{document}