Ultracold atoms at unitarity within quantum Monte Carlo

In VMC the energy is calculated as the expectation value of the Hamiltonian with an approximate many-body trial wave function. In the more accurate DMC method the ground state energy is obtained by evolving the wave function in imaginary time so that it decays towards the ground state. Projector methods such as DMC suffer from a fermion sign problem, which is evaded by making the fixed-node approximation, and importance sampling is introduced to reduce the statistical noise. The importance-sampled fixed-node fermion DMC algorithm was first used by Ceperley and Alder to study the electron gas.Ceperley:PRL:1980

The key quantity in VMC and DMC calculations is the trial wave function, which controls both the accuracy that can be obtained and the statistical efficiency of the computation. In VMC the accuracy of the energy estimate depends on the whole of the trial wave function, while in DMC it depends only on the form of its nodal surface, as the DMC algorithm (in principle) gives the lowest energy compatible with the fixed nodal surface. In practice the DMC energy estimate also shows some dependence on the timestep used and a very weak dependence on the population of particle configurations. Improving the trial wave function tends to reduce these biases and improve the statistical efficiency of the calculations.

The VMC algorithm generates particle configurations distributed according to

where $\Psi_{\rm T}$ is a real trial wave function and ${\bf{R}}$ is the 3$N$-dimensional vector of the coordinates of the $N$ particles. DMC generates configurations distributed according to

where $\Phi({\bf R})$ is the DMC wave function. The total energy in both VMC and DMC is calculated from

where $E_{\rm L}({\bf R})=\Psi_{\rm T}^{-1}\hat{H}\Psi_{\rm T}$ is the local energy, and $p=p_{\rm V}$ in VMC and $p=p_{\rm D}$ in DMC.

DMC expectation values of operators which do not commute with the Hamiltonian depend on the entire trial wave function, not just its nodal surface. To reduce this bias, at the expense of increasing the noise, one can use the extrapolation approximation,Ceperley:book:1986

where $\left<\hat{A}\right>_{\rm DMC}$ and $\left<\hat{A}\right>_{\rm VMC}$ are the DMC and VMC expectation values of operator $\hat{A}$, respectively. The quantity $\left<\hat{A}\right>_{\rm DMC}-\left<\hat{A}\right>_{\rm VMC}$ gives a measure of the accuracy of the trial wave function.

The construction of accurate trial wave functions for Fermi atomic gases is not straightforward. The variation of the wave function must be described very accurately at small inter-particle separations where the interaction potential varies very rapidly. The binding energy of an isolated molecule is vanishingly small at resonance, but the exact value of $\xi$ for the gas is certainly smaller than the BCS mean-field value of 0.59, and therefore the interactions between molecules are very important at unitarity. The exact wave function for an isolated pair of opposite spin fermionic atoms at resonance decays as the inverse distance between the particles, and we must describe the deviations from this behavior in the gas phase. We conclude that it is necessary to provide a good description of both the long- and short-ranged behavior of the pairing function to obtain accurate results for the system at unitarity.

The simulations are performed with a finite number of particles, and we wish to obtain results which accurately reflect those that would be obtained with an infinite number. A great deal of experience has been gained in performing QMC calculations for electron gases, and the finite size effects are greatly reduced by averaging the energies obtained at different wave vectorsLin:PRE:2001 (“twist averaging”) and extrapolating the averaged energies to the infinite system limit. The QMC studies of ultracold atoms reported so far have not employed very large numbers of particles and have not used twist averaging, and it is not clear whether the results are converged with respect to system size. Although we have not used twist averaging in the calculations reported here, we have used systems with 128 particles, which is approximately twice the largest number used in previous DMC calculations.Carlson:PRL:2003 ; Astrakharchik:PRL:2004 ; Chang:PRA:2004 ; Astrakharchik:PRL:2005 ; Chang:PRL:2005 ; Carlson:PRL:2005 ; Gezerlis:PRL:2009

Another problem is that we really want the solution for a delta-function potential, but for computational reasons we use a well of finite width. The ground state of the many-particle system with a delta-function potential is a molecular gas because bound states with more than two particles cannot exist in this case but, with a finite well-width, clusters of particles can form at high densities. In practice this instability of the gas phase occurs only for densities where $\mu r_{s}\ll 1$, and we work at much lower densities. Nevertheless, it is clear that the results of QMC calculations will depend on both the well width and the particle density.

We used a singlet-pairing BCS-like wave function consisting of a determinant of identical pairing orbitals, $\varphi\left(r_{ij}\right)$, each of which is a function of the separation of an up- and a down-spin particle. The pairing orbital is represented by a sum of polynomial terms,

We set $\varphi$ to zero for $r_{ij}$ greater than the cutoff length $L_{\rm P}$. The third power in Eq. (8) was chosen to ensure that the local energy is continuous at $L_{\rm P}$. The $\gamma_{n}$ are optimizable parameters, but the value of $\gamma_{1}$ is determined by the condition that $\varphi$ is cuspless at the origin. We tested various values of $N_{\rm P}$ and chose $N_{\rm P}=4$ for the results presented here. The optimized value of the cutoff length of $L_{\rm P}=3.5$ a.u. is 3.5 times the average distance between particles. We also tested orbitals represented by linear combinations of Gaussian orbitals, linear combinations of plane waves, and linear combinations of Gaussian orbitals, plane waves, and polynomials. Gaussian orbitals also appeared to be a useful choice, but we chose the polynomial basis as the orbitals and their derivatives can be evaluated much more rapidly.

The determinant of pairing orbitals is multiplied by a Jastrow factor of the form $e^{J}$, with

where $N_{\rm s}$ is the number of stars of symmetry-related reciprocal lattice vectors ${\bf{G}}$, and $N_{\rm J}$ is the order of the polynomial. The $\lambda_{s}$ and $\theta_{n}$ are optimizable parameters, although $\theta_{1}$ is determined by the condition that $J$ is cuspless at the origin, and we set the polynomial part of $J$ to zero for $r_{ij}>L_{\rm J}$.Drummond:PRB:2004 After some testing we chose $N_{\rm s}=4$ an $N_{\rm J}=8$, and we used an optimized value of $L_{\rm J}=0.86$ a.u.

We also applied backflow a transformation to the determinant of pairing orbitals.Feynman:PR:1954 ; Feynman:PR:1956 The particle coordinates ${\bf{r}}_{i}$ are replaced by collective coordinates ${\bf{x}}_{i}({\bf{R}})={\bf{r}}_{i}+{\bm{\zeta}}_{i}({\bf{R}})$, where ${\bm{\zeta}}_{i}({\bf{R}})$ is the backflow displacement of particle $i$, which depends on all the particle positions. The backflow displacement is given by

We have used the form

where $L_{\rm B}$ is a cutoff length and $N_{\rm B}$ is the order of the polynomial in the backflow expansion and the $\rho_{n}$ are optimizable parameters, with $\rho_{1}$ determined by the condition that $\eta$ is cuspless at the origin.Lopez-Rios:PRE:2006 We chose $N_{\rm B}=5$ and used an optimized value of $L_{\rm B}=1.04$ a.u.

The wave functions were optimized within a VMC procedure by minimizing the mean absolute deviation of the local energies from their median value. We found this approach to be superior to variance minimization schemes.Kent:PRB:1999 ; Drummond:PRB:2005 The polynomial term in the Jastrow factor is much more important than the plane-wave part. The Jastrow and backflow functions can, in principle, operate between both parallel and anti-parallel spin particles, although the correlation effects between the parallel-spin particles beyond the exchange interaction already included in the determinant are small. We did, however, find a small lowering of the VMC energy when we allowed the plane wave parameters in the Jastrow factor to be non-zero for parallel-spin particles. We did not include parallel-spin terms for the polynomials. The wave function contains a total of 28 variable parameters.

We used the CASINO codecasino for all of our QMC calculations. We performed some test calculations with 32 and 64 particles, but all of our results reported in this paper were obtained with 128 particles. We used a time step of 0.015 a.u. for all the DMC results presented in this paper. Test calculations using a timestep of 0.03 a.u. did not change the total energy within the statistical error bars achieved. We used a mean population of 3200 configurations, and tests indicated that the population control bias with this number of configurations is negligible.

When evaluating the $\xi$ parameter it is not immediately obvious whether to use the non-interacting energy $E_{\rm NI}$ of the finite system studied, the infinite system, or some other value. Energies of non-interacting systems for various particle numbers and fcc and simple cubic (sc) simulation cells are shown in Fig. 1. $E_{\rm NI}$ oscillates in an irregular manner about the infinite system value as the particle number $N$ is increased and converges rather slowly with $N$. Earlier DMC calculations of $\xi$ used sc cells and particle numbers from $N=14$ to 66.Carlson:PRL:2003 ; Astrakharchik:PRL:2004 ; Chang:PRL:2005 ; Carlson:PRL:2005 The error in $E_{\rm NI}$ for the 66-particle system is quite small (0.5%), although it increases to nearly 5% for the sc cell with 80 particles. The rate of convergence for the fcc cell is similar to that for the sc cell.

The convergence of the interacting energy with system size has not been well-studied, but using more particles is expected to improve the results. The momentum distribution of the non-interacting system has a discontinuity at the Fermi momentum, while that of the interacting system is smooth, see Sec. IV.3, so that the kinetic energy of the non-interacting system varies rapidly in k-space, leading to rapid fluctuations in the kinetic energy with system size. It therefore seems likely that the interacting energy will converge faster with system size than the non-interacting energy.

Our values of $\xi$, and those from other calculations and some experiments, are given in Table 1. The DMC energy is bounded from above by the VMC energy and from below by the exact energy and, as expected, the VMC values of $\xi$ are a little larger than the DMC ones. The results reported in Table 1 were obtained using the value of $E_{\rm NI}$ for the finite system studied of 1.08307 a.u. per particle, compared with the infinite-system result of 1.10495 a.u. per particle. As discussed above, one might argue that it is more appropriate to calculate $\xi$ using the infinite-system value of $E_{\rm NI}$, as the interacting energy is expected to converge faster with system size than the non-interacting energy. Using the infinite-system non-interacting energy gives values of $\xi=0.4428(1)$ in VMC, and $\xi=0.4253(1)$ in DMC, which are about 2% lower than those reported in Table 1. Our DMC value of $\xi$ is very similar to previous ones.Carlson:PRL:2003 ; Astrakharchik:PRL:2004 ; Chang:PRL:2005 ; Carlson:PRL:2005

We evaluated the spatially and rotationally averaged pair correlation functions (PCFs) for the parallel and anti-parallel spin pairs, which are shown in Figs. 2 and 3. The difference in our VMC and DMC results was negligible on the scales of Figs. 2 and 3, and consequently the effect of using the extrapolation of Eq. (7) is negligible. Figures 2 and 3 show the DMC data itself, not fits to the data. The noise in the data is very small, but can just be resolved in Fig. 2 at small $r/r_{s}$, where the number of counts is small. The parallel-spin PCF shown in Fig. 2 shows a hole largely confined to the region $r/r_{s}<2$, which is essentially an exchange hole, and the PCF is very similar to the non-interacting (Hartree-Fock) result. This is consistent with the fact that we found only a very weak parallel-spin Jastrow factor. The anti-parallel-spin PCF (Fig. 3) shows a very strong enhancement for $r/r_{s}<1$ arising from the pairing. The behavior at small $r/r_{s}$ is not shown as it depends strongly on the well width. The anti-parallel PCF dips below unity in the region $1<r/r_{s}<2$. The PCFs are similar to those reported in Fig. 3 of Astrakharchik et al.Astrakharchik:PRL:2004 and Fig. 1 of Chang and Pandharipande.Chang:PRL:2005

The one-body density-matrix (OBDM) may be written as

and the two-body density-matrix (TBDM) as

where $p({\bf R})$ is the VMC or DMC probability distribution of Eqs. (4) or (5), ${\bf r}_{1}$ and ${\bf r}_{1}^{\prime}$ are $\alpha$-spin particle coordinates, ${\bf r}_{2}$ and ${\bf r}_{2}^{\prime}$ are $\beta$-spin particle coordinates, and $N_{\alpha}=N_{\beta}=N/2$ is the number of particles of each spin type.

We have evaluated the translationally and rotationally averaged density matrices, which we denote by $\rho_{\alpha}^{(1)}(r)$ and $\rho_{\alpha}^{(2)}(r)$, respectively. The momentum density $n(k)$ is the Fourier transform of $\rho_{\alpha}^{(1)}(r)$, but we evaluate it directly in Fourier space, which is a somewhat better numerical approach. Our data for $n(k)$ at unitarity are broadly similar to those presented in Fig. 2 of Ref. Carlson:PRL:2003, and Fig. 2 of Ref. Astrakharchik:PRL:2005, , although in detail there are some differences. Our data show a monotonic decrease of $n(k)$ with increasing momentum, in common with the results of Ref. Astrakharchik:PRL:2005, , but in conflict with Ref. Carlson:PRL:2003, , which show a peak below the Fermi momentum. Carlson et al.Carlson:PRL:2003 show VMC and DMC data for 14 and 38 particles. The VMC and DMC data are in good agreement above $k_{\rm F}$, but below $k_{\rm F}$ there are substantial differences. The differences between our VMC and DMC data are very small, and would be barely visible on the scale of Fig. 4. This suggests that our trial wave functions are superior to those of Ref. Carlson:PRL:2003, . Astrakharchik et al.Astrakharchik:PRL:2005 only report DMC data in the form of a fit, rather than giving the calculated values. Our momentum density is a little closer to the BCS form than the curve shown in Ref. Astrakharchik:PRL:2005, .

The OBDM is shown in Fig. 5. Again, the VMC and DMC data are virtually indistinguishable so that extrapolation is unnecessary. Our calculated OBDM is very similar to that shown in Fig. 1 of Ref. Astrakharchik:PRL:2005, .

The condensate fraction $c$ is related to the translationally and rotationally averaged TBDM by

where $\Omega$ is the volume of the simulation cell and $N/2$ is the number of pairs of particles in the system. VMC and DMC data for the TBDM are shown in Fig. 6. In this case we find some differences between the VMC and DMC results. We obtained a condensate fraction of $c=0.57$ in VMC and $c=0.54$ in DMC, see Fig. 6, so that the value obtained using the extrapolation of Eq. (7) is 0.51. Our DMC values are lower than those of $0.61(2)$ for 38 particles and $0.57(2)$ for 66 particles reported by Astrakharchik et al.Astrakharchik:PRL:2005

As mentioned in Sec. II, we require $\mu r_{s}\gg 1$ for dilute conditions. This can be satisfied by, for example, fixing $r_{s}$ and choosing the effective width of the potential well $2/\mu$ to be sufficiently small, or by fixing $2/\mu$ and choosing $r_{s}$ to be sufficiently large. We have tried both of these approaches, but did not obtain a smooth variation of the energy when reducing the range of the interaction, at least partly because the wave function optimization becomes more difficult. We obtained smoother results when increasing the value of $r_{s}$ while keeping $\mu=12$, as shown in Table 2, and the optimizations worked well in these calculations. The value of $\xi$ slowly increases with $\mu r_{s}$, but this behavior could arise from the fixed-node error inherent in the trial wave functions. Note that the differences between the VMC and DMC energies increase with $\mu r_{s}$, indicating that the trial wave functions are becoming less accurate. Increasing the size of the simulation cell extends the range of the trial wave function and makes it more difficult to represent.