\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 Monte Carlo for BEC with $^{87}$Rb, deadline  april 16  12am (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 the Variational Monte
Carlo (VMC) method 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.

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 importance sampling
\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 derivation.
Compare the CPU time difference.  You should also  parallelize your code.
The only variational parameter is $\alpha$. Perform these calculations for $N=10$, 
$100$ and $500$ 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 Refs.~\cite{dubois2001,jon2005} $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. 
Discuss your results and comment on eventual differences between importance sampling and brute force sampling.

In the statistical analysis of your results you should use for example the blocking technique from
Project 1.
Your code should reproduce the results of Refs.~\cite{dubois2001,jon2005}.

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

\bibitem{jon2005}   J.~K.~Nilsen, J.~Mur-Petit, M.~Guilleumas, M.~Hjorth-Jensen, 
  and A.~Polls, Phys.~Rev.~{\bf A71}, 053610 (2005).

\end{thebibliography}

\end{document}





