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} $$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*} $$as
$$ \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.
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.
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.
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.
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.
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.
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].
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?
Here follows a brief recipe and recommendation on how to write a report for each project.
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:
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.