Vortices in atomic Bose-Einstein condensates in the large-gas-parameter region

J. K. Nilsen, 1 J. Mur-Petit, 2 M. Guilleumas, 2 M. Hjorth-Jensen, and A. Polls Department of Physics, University of Oslo, N-0316 Oslo, Norway Departament d’Estructura i Constituents de la Matèria, Universitat de Barcelona, E-08028 Barcelona Center of Mathematics for Applications, University of Oslo, N-0316 Oslo, Norway PH Division, CERN, CH-1211 Geneve 23, Switzerland Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, USA sReceived 10 December 2004; published 19 May 2005 d


I. INTRODUCTION
Most 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 ͓1͔. The key point for the validity of this description is the dilute condition of these systems, i.e., the average distance between the atoms is much larger than the range of the interatomic 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͑r͒ = n͑r͒a 3 , where n͑r͒ is the local density of the system. For low values of the average gas parameter x av ഛ 10 −3 , the mean-field Gross-Pitaevskii equation does an excellent job ͑see, for example, Ref. ͓2͔ for a review͒. 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 Feshbach resonance ͓3,4͔.
Under such circumstances it is unavoidable to test the accuracy of the GP equation by performing microscopic calculations. If we consider cases where the gas parameter has been driven to a region were one can still have a universal regime, i.e., that the specific shape of the potential is unimportant, we may attempt to describe the system as dilute hard spheres whose diameter coincides with the scattering length. However, the value of x is such that the calculation of the energy of the uniform hard-sphere Bose gas would require to take into account the second term in the low-density expansion ͓5͔ of the energy density where m is the mass of the atoms treated as hard spheres. For the case of uniform systems, the validity of this expansion has been carefully studied using diffusion Monte Carlo ͓6͔ and hypernetted-chain techniques ͓7͔.
The energy functional associated with the GP theory is obtained within the framework of the local-density approxi-mation ͑LDA͒ by keeping only the first term in the lowdensity expansion of Eq. ͑1͒ is the confining potential defined by the two angular frequencies Ќ and z . The condensate wave function ⌿ is normalized to the total number of particles. By performing a functional variation of E GP ͓⌿͔ with respect to ⌿ * one finds the corresponding Euler-Lagrange equation, known as the Gross-Pitaevskii ͑GP͒ equation where is the chemical potential, which accounts for the conservation of the number of particles. Within the LDA framework, the next step is to include into the energy functional of Eq. ͑2͒ the next term of the low-density expansion of Eq. ͑1͒. The functional variation gives then rise to the so-called modified GP equation ͑MGP͒ ͓8͔ Rb atoms to tune their scattering length ͓3͔. Fully microscopic calculations using a hard-spheres interaction have also been performed in the framework of variational and diffusion Monte Carlo methods ͓10-13͔.
In this work we compare the results of the GP and MGP equations discussed above, Eqs. ͑4͒ and ͑5͒, with variational Monte Carlo ͑VMC͒ calculations for axially symmetric traps in a region ͑x Ͼ 10 −3 ͒, where the validity of the GP equation is not clear. We examine both the ground state and excited states having a vortex line along the z axis.
In the next section we present our numerical approaches together with a discussion of ground-state properties. In Sec. III we proceed to study several trial wave functions to describe the excited state with one vortex. A comparison between VMC and the GP and MGP equations is done. We summarize our results in Sec. IV.

II. NUMERICAL APPROACHES AND GROUND-STATE PROPERTIES
The starting point of the Monte Carlo calculations is the Hamiltonian for N trapped interacting atoms given by The two-body interaction V int ͉͑r i − r j ͉͒ between the atoms is described by a hard-core potential of radius a, where a is the scattering length. The atoms are thus treated as hard spheres. The next step is to define a trial wave function where F͑1,… , N͒ is a many-body correlation operator applied to the mean-field wave function ⌿ MF . The advantage of using a correlated trial wave function lies in the fact that nonperturbative effects, as the short-range repulsion between atoms may be directly incorporated into the trial wave function. The simplest correlation operator has the Jastrow form ͓14͔, F͑1, ... ,N͒ = ͟ iϽj f͑r ij ͒. ͑8͒ In our variational calculations we use a two-body correlation function, which is the solution of the Schrödinger equation for a pair of atoms at very low energy interacting via a hardcore potential of diameter a. The ansatz for the correlation function f͑r͒ reads This type of correlation, besides being physically motivated, has been successfully used in Refs. ͓10,11͔ to study both spherically symmetric and deformed traps. These authors have also explored the quality of this correlation function by comparing variational Monte Carlo ͑VMC͒ and diffusion Monte Carlo ͑DMC͒ calculations for the case of spherically symmetric traps ͓12͔, with a good agreement between the VMC and DMC results. The deformation of the trap is incorporated in the meanfield wave function ⌿ MF , which is taken as the product of N single-particle wave functions ͑r͒ = A͑␣͒ 1/4 exp͓− 1 2 ␣͑x 2 + y 2 + z 2 ͔͒, ͑10͒ where ␣ is taken as the variational parameter of the calculation, and A͑␣͒ = ͑␣ / ͒ 3/4 is the normalization constant. The parameter = z / Ќ is kept fixed and set equal to the asymmetry of the trap. In this way the mean-field wave function ⌿ MF has all the particles in the condensate, the latter being described by the wave function . The evaluation of the expectation value of the Hamiltonian with this correlated trial wave function provides an upper bound to the ground-state energy of the system This expectation value has been evaluated by the Metropolis Monte Carlo method of integration ͓15,16͔.
The energy obtained with the Hamiltonian of Eq. ͑6͒ can be directly compared to the output of the GP and MGP equations, see Eqs. ͑4͒ and ͑5͒. The Gross-Pitaevskii equations represent a mean-field description, with all the atoms in the condensate. In fact, the additional correlations, which are taken into account in the second-order term of the lowdensity expansion of the energy ͓see Eq. ͑1͔͒, are incorporated in the density functional and, therefore, in the solution of the MGP equation. In contrast, the Monte Carlo calculation explicitly incorporates the interatomic correlations, and therefore one could, in principle, find the natural orbits and extract the occupation of the condensate ͓10͔.
The GP and MGP equations have been solved by the steepest descent method ͓17͔ for the deformed harmonic oscillator trap previously described in Eq. ͑3͒. An initial deformed trial state is projected onto the minimum of the functional by propagating it in imaginary time. In practice, one chooses a small time step ⌬t and iterates the equation by normalizing ⌿ at each iteration. When the gas parameter becomes large, the time step, which governs the rate of convergence, should be taken accordingly small. Convergence is reached when the chemical potential becomes a constant independent of the position, see Eqs. ͑4͒ and ͑5͒.
For the comparison of the results obtained with the different GP-type equations and the variational Monte Carlo calculations, we consider a disk-shaped trap with = z / Ќ = ͱ 8, see Ref. ͓18͔. We have fixed the scattering length to a =35a Rb , with a Rb = 100a 0 , a 0 being the Bohr radius. We set the number of confined atoms to N = 500 in order to keep the amount of computing time acceptable when using the Monte Carlo method. All the numerical results are given in units of the harmonic oscillator length a Ќ = ͑ប / m Ќ ͒ 1/2 and the harmonic oscillator energy ប Ќ .
First we analyze the GP and MGP results reported in Table I. For a scattering length a =35a Rb , the corrections of the MGP approach to the chemical potential are of the order of 20%. The energy corrections are also relevant, and it is interesting to study the different contributions to the energy. The kinetic energy is given by while the harmonic oscillator energy due to the trapping potential reads and the interaction energies E 1 and E 2 are given by The virial theorem is used to establish a relation between the different contributions to the energy, viz., which serves as a proof of the numerical accuracy of the solution of the GP equations. The results in Table I show that this test is well satisfied by all calculations. Note that the kinetic energy associated with the meanfield descriptions is not negligible, indicating that the regime where the Thomas-Fermi approximation to the GP equation is valid has not been reached. In this limit, the chemical potential is where ā = a / a Ќ is the dimensionless scattering length, and the energy per particle E TF / N =5 TF / 7. In this approach we have E TF / N = 9.03ប Ќ and TF = 12.64ប Ќ . Both these values differ from the values reported in Table I. However, this approximation can still be used to estimate the peak value of the gas parameter, namely, which yields x TF pk = 0.023. At this rather large value of the diluteness parameter, the corrections brought by the MGP equation to the GP results are expected to be relevant ͓6,7,9͔.
However, x is low enough to allow for a mean-field approach ͑as it is the case of the MGP equation͒. For such density regimes, a mean-field approach provides a rather good description when compared to a microscopic calculation ͓8͔.
The variational Monte Carlo results are also given in Table I and show a close agreement with the results provided by the MGP equation. Note that in this approach, and using the Hamiltonian of Eq. ͑6͒, the potential energy is zero since the wave function is strictly zero inside the core. The total energy in this case is distributed between E HO and the true kinetic energy. Actually the only energies that can be directly compared to the GP results are the total and the harmonic oscillator energies.
The Monte Carlo results obtained with the Metropolis algorithm take into account the energy of 27 000 configurations, grouped in 90 blocks of 300 movements. At each Monte Carlo step we move all the particles and the acceptance is around 58%. A thermalization process is incorporated at the beginning of the Monte Carlo process and before each block. In the Monte Carlo calculation we have used the Pandharipande-Bethe prescription for the kinetic energy ͓16͔, which produces a smaller variance. To get a feeling for the numerical accuracy of our VMC results, we list here GP, MGP, and VMC results in the dilute limit. We employ N = 500 particles and a scattering length for 87 Rb considered by Dalfovo and Stringari ͓19͔, which in units of the oscillator parameter perpendicular to the z axis is 4.33ϫ 10 −3 . We obtain energies in units of the oscillator energy of 3.303 2151, 3.308 0392, and 3.324 1881 ͑10͒ for GP, MGP, and VMC calculations, respectively. The VMC results are for an optimum variational parameter ␣ = 0.475. Taking into account that the two-body correlation has been kept fixed, and that the only variational parameter is ␣, these results indicate that our ansatz for the variational wave function is a viable one. Actually, as the reader will note from the discussion below, this discrepancy of roughly 0.5% is of the same relative order as for the higher density cases reported here.
In the minimization process we keep fixed the parameter in the single-particle wave function of Eq. ͑10͒, i.e., we assume that the deformation of the trap is transferred to the wave function, and vary only ␣. At the minimum, ␣ = 0.7687. One can also explore the effects of the correlations in the density profiles. These profiles, which represent a column density defined according to n c ͑r Ќ ͒ = ͵ dz n͑r Ќ ,z͒ ͑ 20͒ and normalized such that 2 ͐ dr Ќ r Ќ n c ͑r Ќ ͒ = 1, are shown in Fig. 1 for the various approximations used in this work. The repulsive character of the correction term of the MGP equation translates into a decrease of the value of the column density at the origin and an increase of the size of the condensate ͓8,9͔. This gives a slightly more extended profile for the MGP approach compared to both the GP and the VMC results. As one can see from Fig. 1, there is a much better agreement between the Monte Carlo and MGP profiles than with the corresponding profile from the GP calculation, particularly at small values of the radial distance where the density is larger. The good agreement between VMC and MGP does not guarantee that these methods give a good description of the system. However, as it was shown in Ref. ͓11͔ for the case of spherical traps, the improvements introduced in the trial wave function by a diffusion Monte Carlo calculation, which, in principle, allows for an exact solution of the manybody problem, are rather small and the variational wave function of Eq. ͑10͒ provides a very good description of the system. Therefore we assume that the same will be true for deformed traps. Furthermore, for these values of the diluteness parameter, the MGP equation is very useful to calculate the energy, chemical potential, and density profiles of the ground state of the system for condensates with larger number of particles, which would be computationally prohibitive for a Monte Carlo calculation.

III. VORTEX STATES
The existence of these excited condensate states is crucial to studies of the superfluid behavior of trapped atomic condensates. In this section we study the effects of correlations in vortex states. We consider a singly quantized vortex line along the z axis. This means that all the atoms rotate around the z axis with a quantized circulation. The GP equation can easily be generalized to describe this kind of vortex states ͓2͔ by using the following ansatz for the condensate wave function where is the angle around the z axis and is an integer. This vortex state has a tangential velocity v = ប mr Ќ , ͑22͒ where r Ќ = ͱ x 2 + y 2 is the distance to the symmetry axis of the vortex. The number represents the quantum of circulation, and the total angular momentum along the z axis is given by Nប. Introducing the wave function of Eq. ͑21͒, in the GP energy functional of Eq. ͑2͒, one gets the corresponding GP energy functional for the vortex state which incorporates a centrifugal term in the density functional, arising from the quantized flow of atoms around the vortex core. This term defines a rotational energy The corresponding nonlinear Schrödinger equation obtained by functional variation is ͑25͒ Adding E 2 to the density functional and after performing a functional variation one gets the corresponding MGP equation for the vortex state.
Based on the virial theorem, one can again derive a relation between the different contributions to the energy 2E kin − 2E HO + 3E 1 + 9 2 E 2 + 2E rot = 0. ͑26͒ The thermodynamic critical angular frequency ⍀ c required to produce a vortex of vorticity is obtained by comparing the energy of the system in the rotating frame with and without the vortex ͓20͔ A main feature of a vortex state is the hole ͑core of the vortex͒ that appears in the center of the density profile along the rotation axis. From Eq. ͑25͒, it is clear that the solution of this equation has to vanish on the z axis because of the presence of the centrifugal term. The size of the core is characterized by the healing length.
For the microscopic description of the vortex state we use an Onsager-Feynman-type trial wave function ͓21͔ where ⌿ 0 ͑1,… , N͒ is the ground-state wave function. The phase factor ͚ j j depends on the angular variables of the particles and is the equivalent to the phase factor introduced in the mean-field description of Eq. ͑21͒. The function f͑r Ќ ͒ modulates the density as a function of the radial coordinate r Ќ . We examine three types of f͑r Ќ ͒. In the first ansatz we use the simple option In the second case we consider, where d is a variational parameter. Note that for d = 1, the behavior of f 2 ͑r Ќ ͒ for small r Ќ coincides with the behavior of f 1 ͑r Ќ ͒. Finally, the third function is that of Ref. ͓22͔, which has been used in the context of quantum liquids, where d is again a variational parameter. These three trial wave functions describe a singly quantized vortex state ͑ =1͒, whose axis lies in the z direction and with a tangential velocity field v = ប / mr Ќ . The evaluation of the expectation value of the Hamiltonian ͓Eq. ͑6͔͒ with these wave functions is equivalent to calculate the mean value of the Hamiltonian with ⌿͑1,… , N͒ = ͟ j f͑r Ќ,j ͒⌿ 0 ͑1,… , N͒. In this way the rotational contribution to the energy has been directly incorporated in the Hamiltonian. Minimizing this new problem provides the best energy and wave functions inside this subspace of wave functions. In the context of liquid 4 He there have been attempts to perform a full minimization allowing for a more general phase function. The analysis indicates that the present procedure provides very accurate results ͓23͔.
We start by discussing the GP and MGP results ͑obtained by the steepest descent method ͓17͔ as done for the ground state as well͒ with an initial condensate wave function It is worth mentioning, as a check of the numerical procedure, that starting with f 2 ͑r Ќ ͒ or f 3 ͑r Ќ ͒ to modulate the condensate wave function we converge to the same results as with f 1 ͑r Ќ ͒.
As expected, the presence of the vortex increases the chemical potential. Also E HO has a small increase, related to the enlargement of the profile because of the presence of the vortex hole. These results are listed in Table II. Although the MGP corrections to the energy are sizable and of the same order as those in the ground state, the critical frequency, ⍀ GP = 0.29 Ќ , is barely affected as both energies, the energy of the vortex state and the ground-state energy, are shifted by similar amounts, yielding ⍀ MGP = 0.24 Ќ .
The GP and MGP profiles for the vortex state are shown in Fig. 2. As a consequence of the repulsive character of the MGP corrections, the central density of the GP ground-state density profile is higher than the MGP one and, therefore, the depth of the hole around the z axis is larger in the GP approach. However, the healing length is almost the same.
As can be seen from Table III, the Monte Carlo results for the energies are in good agreement with the MGP ones for all the trial wave functions considered. This table shows two types of calculations. In the first three rows we list the energies obtained by keeping ⌿ 0 equal to the ground-state wave function and performing the minimization with respect to the parameter d in the modulating function, except in the case of f 1 , which has no variational parameters. In the second set of results, we perform a minimization allowing to vary also the harmonic oscillator parameter ␣ of the wave function ⌿ 0 . The changes in ␣ and d do not yield significant changes in the computed energy.
The density profiles seem to be more sensitive to the modulating function, as one can see from Fig. 2. These profiles correspond to the case where the ground-state wave function ⌿ 0 is kept fixed when we minimize the energy of the vortex state. For f 3 we obtain a radial structure, which is not present in the mean-field approach ͓22͔. The MGP profile shows a broader surface region than the VMC profiles. In the core of the vortex, the MGP profile looks very similar to the VMC results with the modulating function f 2 of Eq. ͑30͒. These two results exhibit a smaller healing length than the VMC calculation which employs f 1 .
From the variational point of view, the best description of the vortex should correspond to the wave function that provides the minimum energy. According to this criterion, this corresponds to the trial wave function built with the modulating function f 1 of Eq. ͑29͒.
Finally, in Fig. 3 we plot the density profiles for all VMC calculations, with and without vortices. We note that they all provide a similar healing length and that the asymptotic behavior is almost equal for both the ground state and the vortex states.

IV. CONCLUSIONS
We have compared the results of the Gross-Pitaevskii ͑GP͒ and the modified Gross-Pitaevskii ͑MGP͒ equations to ab initio variational Monte Carlo calculations for Bose-Einstein condensates of atoms in deformed traps. We have studied both the ground state and excited states having a vortex line along the z axis. The interatomic potential has been characterized by a hard-sphere potential with a radius that coincides with the scattering length used in the GP and MGP equations.
We have performed the calculations for 500 particles. The parameters characterizing the trap and the scattering length have been chosen to reach values of the gas parameter where the MGP calculations provide corrections of the order of 20% compared to the GP results. It is indeed very interesting that even at such values of the gas parameter one can still describe the system in terms of mean-field approaches. We find, for example, an excellent agreement between the MGP and VMC results, especially for the energies of the ground state and the vortex states. The MGP and VMC density profiles for the ground state are also in good agreement. The situation is different for the vortex state. Three different trial wave functions produce similar energies but slightly different profiles. In the core of the vortex, the MGP profile is close to the profiles obtained with the ansatzes f 1 and f 2 of Eqs. ͑29͒ and ͑30͒, respectively. These functions yield also the lowest energies. Whether a diffusion Monte Carlo ͑DMC͒ calculation will show a similar trend remains to see. We are planning DMC studies of the systems discussed here. Our preliminary DMC calculations for the energy of the ground state show little change with respect to the VMC results and, hence, a very good agreement with the MGP results.
In summary, we would like to point out that the good agreement between the VMC and MGP is rather encouraging and allows for further MGP explorations of vortex states in condensates with both a larger number of interacting atoms and large scattering lengths.  3. Column density n c ͑r Ќ ͒ as a function of the distance to the z axis, comparing the VMC profiles for the vortex state corresponding to the different Onsager-Feynman ansatzes ͑lines with symbols as in Fig. 2͒ and the ground state ͑dashed line with full circles͒. The trap parameters and the scattering length are the same as in the two preceding figures. The radial distance is given in units of a Ќ = ͑ប / m Ќ ͒ 1/2 . The column density is dimensionless. See text for further details.