Beyond Gross-Pitaevskii:local density vs. correlated basis approach for trapped bosons

We study the ground state of a system of Bose hard-spheres trapped in an isotropic harmonic potential to investigate the effect of the interatomic correlations and the accuracy of the Gross-Pitaevskii equation. We compare a local density approximation, based on the energy functional derived from the low density expansion of the energy of the uniform hard sphere gas, and a correlated wave function approach which explicitly introduces the correlations induced by the potential. Both higher order terms in the low density expansion, beyond Gross-Pitaevskii, and explicit dynamical correlations have effects of the order of percent when the number of trapped particles becomes similar to that attained in recent experiments.

The recent discovery of Bose-Einstein condensation (BEC) of magnetically trapped alkali atoms has generated a huge amount of theoretical investigations. A review of the present situation can be found in refs. 1,2 . Present experimental conditions are such that the atomic gas is very dilute, i.e., the average distance among the atoms is much larger than the range of the interaction. As a consequence, the physics should be dominated by two body collisions, generally well described in terms of the s-wave scattering length. The case of a positive scattering length is equivalent to consider a very dilute system of hard spheres (HS), whose diameter coincides with the scattering length itself. So, the Gross-Pitaevskii (GP) theory for weakly interacting bosons seems the logical tool to study these systems and most of the present days theoretical work is founded on (or has its starting point in) this theory 3 . However, in very recent experiments the number of trapped atoms has spectacularly increased reaching N values of the order 10 6 and 10 7 atoms 4 . Therefore, it seems logical to ask for a deeper study of the effect of the interatomic correlations and of the accuracy of the GP scheme in this new scenario.
In order to address these questions, we study the ground state of a system of Bose hard spheres trapped by a harmonic potential. More precisely, we consider hard spheres with a diameter of 52.9Å, corresponding to the s-wave triplet-spin scattering length of 87 Rb, in an isotropic harmonic trap characterized by an angular frequency ω/2π = 77.78Hz. We also examine the large N atomic sodium case of ref. 4 . We use and compare two methods: (i) a local density approximation (LDA) based on an energy functional derived by the lowdensity expansion of the energy of an uniform hard sphere gas; (ii) a correlated basis function (CBF) approach which explicitely takes into account the dynamical correlations induced by the potential and which is not, in principle, limited to purely repulsive interactions.
LDA theory. In the case of an uniform hard sphere Bose gas, a perturbation theory in the expansion parameter x = na 3 , where n is the density of the system, leads to the following low density (LD) expansion for the energy density 5 : The energy functional associated with the GP theory is simply obtained in LDA by keeping only the first term in the expansion (1): where the wavef unction ψ is normalized to the total number of atoms. By performing a functional variation of E GP [ψ] one finds the Euler-Lagrange equation, known as the GP equation where µ is the chemical potential. This equation has the form of a nonlinear stationary Schrödinger equation, and it has been solved for several types of traps using different numerical methods [6][7][8][9] .
A logical step further, in the spirit of LDA, is to include into the energy functional the next terms of the correlation energy (1). It is convenient to simplify the notation by expressing lengths and energies in harmonic oscillator (HO) units. The HO length is a HO = (h/mω) 1/2 , and the spatial coordinates, the energy and the wave function are rescaled as r = a HOr , E =hωE 1 and ψ(r) = (N/a 3 HO ) 1/2 ψ 1 (r), where ψ 1 (r) is normalized to unity. By taking into account the next terms of the correlation energy, LDA provides a modified GP (MGP) energy functional and a modified Gross-Pitaevskii equation whereā = a/a HO and µ 1 is the chemical potential in HO units. The MGP equation contains extra nonlinar terms in ψ 1 . CBF approach. CBF theory is a powerful tool to study strongly interacting many-body systems (for a review, see ref. 10 ). In particular, it was used to study the HS homogenous Fermi gas in order to ascertain the accuracy of low density type expansions 11 .
For N interacting bosons, at T=0 temperature, described by the hamiltonian where V ext (r i ) is the confining potential and V (r ij ) is the interatomic potential, the CBF ground-state wave function is where F (1, .., N) is a many-body correlation operator applied to the mean field w.f. ψ M F . The advantage of using a CBF basis lies in the fact that non perturbative effects, as the short range repulsion for the hard spheres, may be directly inserted into the correlation. The simplest correlation operator has the Jastrow form 12 where the Jastrow correlation function, f J (r), depends on the interparticle distance only. f J (r) is variationally determined by minimizing the expectation value of E CBF = ψ CBF |H|ψ CBF / ψ CBF |ψ CBF . E CBF may be evaluated either by MonteCarlo techniques or by cluster expansions of the needed one-and two-body densities, ρ 1 (r 1 ) and ρ 2 (r 1 , r 2 ). The energy per particle can be written as and the correlation energy, V corr , is where ρ 1 is normalized to N. ρ 2 may be calculated by cluster expansion in powers of the dynamical correlation, h(r) = f 2 J (r)−1, and the integral Hypernetted Chain (HNC) equations 13 may be used to sum infinite classes of cluster diagrams.
Given the low-density of the trapped bosons system, it looks plausible to start by a lowest order (LO) cluster expansion. In this approximation ρ CBF with respect to ρ 1 leads to the LO-correlated Hartree with s = |r − r 1 |. The optimal choice for the Jastrow factor would be the one satisfying the Euler equation δE CBF /δf J (r) = 0. Otherwise, parametrized functional forms may be chosen whose parameters are found through the minimization process. We adopt here the correlation function minimizing the lowest order energy of a homogenous Bose gas with a healing condition at a distance d (taken as variational parameter). For the HS case, f J (r < a) = 0 and f J (r > a) = u(r)/r where u(r) is solution of the Schrödinger-like equation −u" = K 2 u. f J (r) has the form 14 where the healing conditions, f J (r ≥ d) = 1 and f ′ J (r = d) = 0, are fixed by the relation: Results. We begin by briefly discussing the boson HS homogeneous case. It was shown in ref. 14 that the lowest order cluster energy of an infinite system of bosonic hard spheres, correlated by the Jastrow factor (14), asymptotically tends to the first term of the LD expansion (LD0) in (1) when the healing distance goes to infinity. We have numerically checked this behavior. The situation changes if all the many-body cluster terms are included via HNC summation, since the computed energy shows a clear minimum in d. Some results are given in Fig.1 for various x-values. The energies have been multiplied by 10 3(2,1) at x = 10 −6(−5,−4) , respectively. The Figure shows the energy estimates computed by retaining different LD expansion terms (LD1 and LD2 correspond to adding the second and third terms in (1)), the HNC energies and, at the highest x-value, the exact energy, evaluated by a diffusion MonteCarlo method 15 . The quality of the HNC results is excellent as they practically coincide with the exact ones at all the densities (the low density value at low x and the DMC one at large x). The convergence of the LD expansion at large x appears to be rather poor, pointing to the relevance of successive contributions.
The GP, MGP and CH/LO equations for ψ 1 have been solved by the steepest descent method 16 for the isotropic harmonic oscillator trap previously described. An initial trial state is projected onto the minimum of the functional by propagating it in imaginary time. In practice, one chooses an arbitrary time step ∆t and iterates the equation by normalizing ψ 1 to 1 at each iteration. When the number of particles becomes very large the time step, which governs the rate of convergence, should be taken accordingly small. The number of iterations substantially increases and it is convenient to start the convergence process from a reasonable wave function (for the GP and MGP equations we start from the ground state harmonic oscillator function which minimizes the GP energy functional, while, in the CH/LO case, we start from the GP solution). We find that the LO approximation in the finite system shows a behavior similar to the infinite case. In particular, the solution of the CHE/LO equation provides an energy (in reduced unities) very close to the GP one when the healing distance becomes very large. In fact, for N = 10 5 andā = 4.33 × 10 −3 , corresponding to the 87 Rb scattering length, we obtain E (LO) 1 /N = 12.57, 12.28 and 12.11 atd = 10, 12 and 15, respectively, while E (GP ) 1 /N = 12.10. This indicates that many-body effects can be recovered only by a HNC treatment also for a finite number of atoms. However, the solution of the fully correlated HNC Hartree equation, even if feasible in principle, is very cumbersome, computationally time consuming and numerically delicate. For this reason we have decided to estimate these effects by adopting a LDA-type approach to V corr . We approximate V corr ∼ V LDA corr where E hom HN C (ρ 1 ) is the HNC homogeneous gas energy per particle at density ρ 1 . The minimization of the energy gives the HNC correlated Hartree equation (CH/HNC) where we have again introduced the scaled unities and x = ρ 1 a 3 = Nā 3 |ψ 1 | 2 . The calculations have been performed for the 87 Rb scattering length. The scaled energies per particle and the root mean square radii are reported in Table I for particle numbers from 10 3 to 10 7 . The Table also shows the results obtained by neglecting the kinetic energy term in the GP equation. This approach, lousely called Thomas Fermi (TF) approximation, has been discussed in the literature and allows for deriving simple analytical expressions 6 . The differences between this TF approach and a rigorous one have been recently discussed 17,18 for spatially inhomogeneous Bose condensates. LDA has been used 1 to estimate corrections to the GP energy within the TF approximation and retaining only the first correction in (1). The second correction is negative and partially cancels the first one. For instance, the cancellations go from ∼ 15% for N = 10 4 to ∼ 40% at N = 10 6 if we just take the central densities, whereas the final energy is reduced by ∼ 15% at N = 10 6 and it is practically unaffected by the second correction at lower N-values.
As expected, the TF results are close to the GP ones when N becomes large. The differences between GP and MGP increase with the number of particles and are of the order of 4% for the chemical potential and 2.5% for the energy at N = 10 7 . The higher order terms in the low density expansion always have a repulsive effect. The same behavior is shown by the HNC results, which, however, are less repulsive than MGP at the large N values.
We notice that if one uses the GP solution to perturbatively estimate the MGP energy, then the correction is negative (at N = 10 7 , ∆E 1 = −4.54). The non linear character of the MGP equation is responsible for this discrepancy.
The density profile (normalized to unity) for N = 10 7 particles is given in Fig.2. For this large number of particles the TF and GP densities are close, whereas the more repulsive MGP and HNC solutions lower the central density, expanding the density distribution and providing a larger radius, as shown in Table I .
We have also considered a system of N = 1.5 × 10 7 Na atoms (a = 27.5Å) in a spherical trap having a frequency of 230 Hz. These conditions roughly correspond to those of the experiment described in ref. 4 . The results are shown in the last row of the Table and in Fig.2. The effects of the correlations are similar to those found in the large N Rb cases. The energy increases by ∼ 1% and the r.m.s. radius by ∼ 0.7% respect to GP. The HNC central density is slightly reduced.
In conclusion, we find that both higher order terms in the low density expansion (beyond the GP approach and evaluated in LDA) and explicit dynamical correlations (induced by the strong repulsion) have a not negligible effect when the number of trapped particles becomes large. So, they must be properly considered.
In this respect, CBF theory may play a preminent role. In addition, it allows for a fully microscopic investigation for any type of potential. This is particularly interesting in view of its application to atoms having negative scattering lengths, as 7 Li. In this case the potential exhibits an actractive part and the GP equation has a metastable solution only if N does not exceed some critical value N c . A recent study 19 , which makes use of an effective interaction, has shown that a new branch of Bose condensate may appear at higher densities. CBF may provide further insights into this problem having access to the full structure of the potential itself.
The introduction of a HNC energy functional derived by the homogeneous HS system provides a quick and probably reliable way of embodying correlation effects into the treatment of the trapped atoms. We see as a particular appeal of this approach the clear possibilty of extending it to non spherical traps, as in real experiments. Moreover, as already stated, potentials other than the simple hard sphere one may be easily handled.