Project 1, deadline March 22, 2024

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

Jan 10, 2024

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_{ho}=\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_0 < 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_{ho} = 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} \end{array} \end{equation} $$

where (S) stands for spherical 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} \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 \eqref{trap_eqn} 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{align} \sum_{i < j}^{N} V_{ij} \equiv \sum_{i = 1}^{N}\sum_{j = i + 1}^{N} V_{ij}, \label{_auto2} \end{align} $$

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} \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} \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} \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} \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} \end{equation} $$

for the above trial wave function of Eq. (5) and defined by the terms in Eqs. (6) and (7).

Find first the local energy the case with only the harmonic oscillator potential, that is we set $ a=0 $ and discard totally the two-body potential.

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} \end{equation} $$

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

Our next step involves the calculation of 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

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

with the above trial wave function of Eq. (5). We rewrite

$$ \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], \end{equation*} $$

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

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

$$ \begin{equation*} f(r_{ij})= \exp{\left(u(r_{ij})\right)}, \end{equation*} $$

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

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

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_{kj})\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.

Note: In parts 1b, 1c, 1d, 1e and 1f you will develop all computational ingredients needed by studying only the non-interacting case. We add the repulsive interaction in the final two parts, 1g and 1h. The reason for doing so is that we can develop all programming ingredients and compare our results against exact analytical results.

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): Finding the best parameter(s)

When we performed the calculations in parts 1b) and 1c), we simply plotted the expectation value of the energy as a function of the parameter $ \alpha $. For large systems, this means that we end up with spending equally many Monte Carlo cycles for values of the energy away from the minimum. We can improve upon this by using various optimization algorithms. The aim of this part, still using only the non-interacting case, is to add to our code either a steepest descent algorithm or a stochastic gradient optmization algorithm in order to obtain the best possible parameter $ \alpha $ which minimized the expectation value of the energy.

Project 1 e): A better statistical analysis

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

A useful strategy here is to write your expectation values to file and then have a Python code which does the final statistical analysis. Alternatively, you can obviously write addition functions to be used by your main program and perform the final statistical analysis within the same code.

Project 1 f): Parallelizing your code

Before we add the two-body interaction, our final computational ingredient is to parallelize our code. With this last ingredient we have obtained a code framework which contains the essential elements used in a Variational Monte Carlo approach to a many-body problem. Dealing with a non-interacting case only till now allows us to continuously check our results against exact solutions.

You should parallelize your code using MPI or OpenMP.

Project 1 g): The repulsive interaction

We are now ready to include the repulsive two-body interaction.

We turn to the elliptic trap with a repulsive interaction. We fix, as in Refs. [1,2] below, $ a/a_{ho}=0.0043 $. We 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

$$ \begin{equation*} 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|). \end{equation*} $$

What is the expression for $ \gamma $? Choose the initial value for $ \beta=\gamma = 2.82843 $ and compute ground state energy using the trial wave function of Eq. (5) using only $ \alpha $ as variational parameter. Vary again 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. Benchmark your results with those of Refs. [1,2].

Project 1 h): 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 canvas to hand in your projects, log in at http://canvas.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.