Project 1, deadline March 22

Computational Physics I FYS4411/FYS9411, Department of Physics, University of Oslo, Norway

Date: Jan 25, 2019

Introduction

The spectacular demonstration of Bose-Einstein condensation (BEC) in gases of alkali atoms $^{87}$Rb, $^{23}$Na, $^7$Li confined in magnetic traps 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$.

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\ . 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. A typical $^{87}$Rb atom density in the trap is $n \simeq 10^{12}- 10^{14}$ atoms per cubic cm, 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 ($na^3_{Rb} \simeq 10^{-6}$ where $n = N/V$ is the number density).

Many theoretical studies of Bose-Einstein condensates (BEC) in gases of alkali atoms confined in magnetic or optical traps have been conducted in the framework of the Gross-Pitaevskii (GP) equation. The key point for the validity of this description is the dilute condition of these systems, that is, the average distance between the atoms is much larger than the range of the inter-atomic interaction. In this situation the physics is dominated by two-body collisions, well described in terms of the $s$-wave scattering length $a$. The crucial parameter defining the condition for diluteness is the gas parameter $x(\mathbf{r})= n(\mathbf{r}) a^3$, where $n(\mathbf{r})$ is the local density of the system. For low values of the average gas parameter $x_{av}\le 10^{-3}$, the mean field Gross-Pitaevskii equation does an excellent job. However, in recent experiments, the local gas parameter may well exceed this value due to the possibility of tuning the scattering length in the presence of a so-called Feshbach resonance.

Thus, improved many-body methods like Monte Carlo calculations may be needed. 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.

This trial 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 one, two and finally three dimensions, with the latter given by

$$ \begin{equation} V_{ext}(\mathbf{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} \tag{1} \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}({\mathbf{r}}_i)\right) + \sum_{i<j}^{N} V_{int}({\mathbf{r}}_i,{\mathbf{r}}_j), \label{_auto1} \tag{2} \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 (1) is $\langle x^2\rangle=(\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}$. Note that we use the shorthand notation

$$ \begin{equation} \sum_{i < j}^{N} V_{ij} \equiv \sum_{i = 1}^{N}\sum_{j = i + 1}^{N} V_{ij}, \label{_auto2} \tag{3} \end{equation} $$

that is, the notation $i < j$ under the summation sign signifies a double sum running over all pairwise interactions once.

We will represent the inter-boson interaction by a pairwise, repulsive potential

$$ \begin{equation} V_{int}(|\mathbf{r}_i-\mathbf{r}_j|) = \Bigg\{ \begin{array}{ll} \infty & {|\mathbf{r}_i-\mathbf{r}_j|} \leq {a}\\ 0 & {|\mathbf{r}_i-\mathbf{r}_j|} > {a} \end{array} \label{_auto3} \tag{4} \end{equation} $$

where $a$ is the so-called hard-core diameter of the bosons. Clearly, $V_{int}(|\mathbf{r}_i-\mathbf{r}_j|)$ is zero if the bosons are separated by a distance $|\mathbf{r}_i-\mathbf{r}_j|$ greater than $a$ but infinite if they attempt to come within a distance $|\mathbf{r}_i-\mathbf{r}_j| \leq a$.

Our trial wave function for the ground state with $N$ atoms is given by

$$ \begin{equation} \Psi_T(\mathbf{r})=\Psi_T(\mathbf{r}_1, \mathbf{r}_2, \dots \mathbf{r}_N,\alpha,\beta) =\left[ \prod_i g(\alpha,\beta,\mathbf{r}_i) \right] \left[ \prod_{j<k}f(a,|\mathbf{r}_j-\mathbf{r}_k|) \right], \label{eq:trialwf} \tag{5} \end{equation} $$

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.,

$$ \begin{equation} g(\alpha,\beta,\mathbf{r}_i)= \exp{[-\alpha(x_i^2+y_i^2+\beta z_i^2)]}. \label{_auto4} \tag{6} \end{equation} $$

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

$$ \begin{equation} f(a,|\mathbf{r}_i-\mathbf{r}_j|)=\Bigg\{ \begin{array}{ll} 0 & {|\mathbf{r}_i-\mathbf{r}_j|} \leq {a}\\ (1-\frac{a}{|\mathbf{r}_i-\mathbf{r}_j|}) & {|\mathbf{r}_i-\mathbf{r}_j|} > {a}. \end{array} \label{_auto5} \tag{7} \end{equation} $$

Project 1 a): Local energy

Find the analytic expressions for the local energy

$$ \begin{equation} E_L(\mathbf{r})=\frac{1}{\Psi_T(\mathbf{r})}H\Psi_T(\mathbf{r}), \label{eq:locale} \tag{8} \end{equation} $$

for the above trial wave function of Eq. (5). Find first the local energy the case with only the harmonic oscillator potential, that is we set $a=0$. Use first that $\beta =1$ and find the relevant local energies in one, two and three dimensions for one and $N$ particles with the same mass.

Compute also the analytic expression for the drift force to be used in importance sampling

$$ \begin{equation} F = \frac{2\nabla \Psi_T}{\Psi_T}. \label{_auto6} \tag{9} \end{equation} $$

Find first the equivalent expressions for the just the harmonic oscillator part in one, two and three dimensions with $\beta=1$.

Next, we will find the local energy for the full problem in three dimensions. The tricky part is to find an analytic expressions for the derivative of the trial wave function

$$ \frac{1}{\Psi_T(\mathbf{r})}\sum_i^{N}\nabla_i^2\Psi_T(\mathbf{r}), $$

with the above trial wave function of Eq. (5). We rewrite (and we can use the same general expressions for projects 2 and 3)

$$ \Psi_T(\mathbf{r})=\Psi_T(\mathbf{r}_1, \mathbf{r}_2, \dots \mathbf{r}_N,\alpha,\beta) =\left[ \prod_i g(\alpha,\beta,\mathbf{r}_i) \right] \left[ \prod_{j<k}f(a,|\mathbf{r}_j-\mathbf{r}_k|) \right], $$

$$ \Psi_T(\mathbf{r})=\left[ \prod_i g(\alpha,\beta,\mathbf{r}_i) \right] \exp{\left(\sum_{j<k}u(r_{jk})\right)} $$

where we have defined $r_{ij}=|\mathbf{r}_i-\mathbf{r}_j|$ and

$$ f(r_{ij})= \exp{\left(\sum_{i<j}u(r_{ij})\right)}, $$

with $u(r_{ij})=\ln{f(r_{ij})}$. We have also

$$ g(\alpha,\beta,\mathbf{r}_i) = \exp{\left[-\alpha(x_i^2+y_i^2+\beta z_i^2)\right]}= \phi(\mathbf{r}_i). $$

Show that the first derivative for particle $k$ is

$$ \begin{align*} \nabla_k\Psi_T(\mathbf{r}) &= \nabla_k\phi(\mathbf{r}_k)\left[\prod_{i\ne k}\phi(\mathbf{r}_i)\right]\exp{\left(\sum_{j<m}u(r_{jm})\right)} \\ &\qquad + \left[\prod_i\phi(\mathbf{r}_i)\right] \exp{\left(\sum_{j<m}u(r_{jm})\right)}\sum_{l\ne k}\nabla_k u(r_{kl}), \end{align*} $$

and find the final expression for our specific trial function. The expression for the second derivative is (show this)

$$ \begin{align*} \frac{1}{\Psi_T(\mathbf{r})}\nabla_k^2\Psi_T(\mathbf{r}) &= \frac{\nabla_k^2\phi(\mathbf{r}_k)}{\phi(\mathbf{r}_k)} + 2\frac{\nabla_k\phi(\mathbf{r}_k)}{\phi(\mathbf{r}_k)} \left(\sum_{j\ne k}\frac{(\mathbf{r}_k-\mathbf{r}_j)}{r_{kj}}u'(r_{ij})\right) \\ &\qquad + \sum_{i\ne k}\sum_{j \ne k}\frac{(\mathbf{r}_k-\mathbf{r}_i)(\mathbf{r}_k-\mathbf{r}_j)}{r_{ki}r_{kj}}u'(r_{ki})u'(r_{kj}) \\ &\qquad + \sum_{j\ne k}\left( u''(r_{kj})+\frac{2}{r_{kj}}u'(r_{kj})\right). \end{align*} $$

Use this expression to find the final second derivative entering the definition of the local energy. You need to get the analytic expression for this expression using the harmonic oscillator wave functions and the correlation term defined in the project.

Project 1 b): Developing the code

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 and one dimension. 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. The only variational parameter is $\alpha$. Perform these calculations for $N=1$, $N=10$, $100$ and $500$ atoms. Compare your results with the exact answer. Extend then your results to two and three dimensions and compare with the analytical results.

Project 1 c): Adding importance sampling

We repeat part b), 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. Run the calculations for the one, two and three-dimensional systems only and without the repulsive potential. Study the dependence of the results as a function of the time step $\delta t$.
Compare the results with those obtained under b) and comment eventual differences.

Project 1 d): A better statistical analysis

In performing the Monte Carlo analysis we will use the blocking and bootstrap techniques to make the statistical analysis of the numerical data. Present your results with a proper evaluation of the statistical errors. Repeat the calculations from exercise c) and include a proper error analysis. Limit yourself to the three-dimensional case only.

Project 1 e): The repulsive interaction

[e)] We turn now to the elliptic trap with a repulsive

interaction. We fix, as in Refs. [1,2] below, $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

$$ 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}(|\mathbf{r}_i-\mathbf{r}_j|). $$

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. (5) 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 exercises. In actual calculations employing the Metropolis algorithm, all moves are recast into the chosen simulation cell with periodic boundary conditions. To carry out consistently the Metropolis moves, it has to be assumed that the correlation function has a range shorter than $L/2$. Then, to decide if a move of a single particle is accepted or not, only the set of particles contained in a sphere of radius $L/2$ centered at the referred particle have to be considered. Benchmark your results with those of Refs. [1,2].

Project 1 f): Finding the best parameter

Repeat the previous calculations by varying the energy using the conjugate gradient method or similar methods to obtain the best possible parameter $\alpha$.

Project 1 g): Onebody densities

With the optimal parameters for the ground state wave function, compute again the onebody density with and without the Jastrow factor. How important are the correlations induced by the Jastrow factor?

Literature

J. L. DuBois and H. R. Glyde, H. R., Bose-Einstein condensation in trapped bosons: A variational Monte Carlo analysis, Phys. Rev. A 63, 023602 (2001).
J. K. Nilsen, J. Mur-Petit, M. Guilleumas, M. Hjorth-Jensen, and A. Polls, Vortices in atomic Bose-Einstein condensates in the large-gas-parameter region, Phys. Rev. A 71, 053610 (2005).

Introduction to numerical projects

Here follows a brief recipe and recommendation on how to write a report for each project.

Give a short description of the nature of the problem and the eventual numerical methods you have used.
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.
Include the source code of your program. Comment your program properly.
If possible, try to find analytic solutions, or known limits in order to test your program when developing the code.
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.
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.
Try to give an interpretation of you results in your answers to the problems.
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.
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.

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 or as an ipython notebook file. As programming language we prefer that you choose between C/C++, Fortran2008 or Python. The following prescription should be followed when preparing the report:

Use Devilry to hand in your projects, log in at http://devilry.ifi.uio.no with your normal UiO username and password.
Upload only the report file! For the source code file(s) you have developed please provide us with your link to your github domain. The report file should include all of your discussions and a list of the codes you have developed. The full version of the codes should be in your github repository.
In your github repository, please include a folder which contains selected results. These can be in the form of output from your code for a selected set of runs and input parameters.
Still in your github make a folder where you place your codes.
In this and all later projects, you should include tests (for example unit tests) of your code(s).
Comments from us on your projects, approval or not, corrections to be made etc can be found under your Devilry domain and are only visible to you and the teachers of the course.

Finally, we encourage you to work two and two together. Optimal working groups consist of 2-3 students. You can then hand in a common report.