Approximate Solutions of the Nonlinear Schrödinger Equation for Ground and Excited States of Bose-Einstein Condensates

I present simple analytical methods for computing the properties of ground and excited states of Bose-Einstein condensates, and compare their results to extensive numerical simulations. I consider the effect of vortices in the condensate for both positive and negative scattering lengths, a, and find an analytical expression for the large-N0 limit of the vortex critical frequency for a > 0, by approximate solution of the time-independent nonlinear Schrödinger equation.


Introduction
The recent production of trapped atomic vapors of 87 Rb [1], 7 Li [2], and 23 Na [3] at the phase space densi-ties necessary to generate Bose-Einstein condensation (BEC) in the ground state of magnetic traps has generated significant interest in this area of physics. Previously, BEC has only been encountered in the He II phase of liquid 4 He, and there has also been accumulating evidence for BEC of an exciton gas in cuprous oxide [4]. However, it now seems that the fraction of atoms that participate in BEC in He II is of the order of 10 %, because of the presence of strong atomic interactions [5,6]. In the dilute alkali gas BECs, such interactions are much weaker, and experimental techniques have been used to produce samples in which the condensate fraction is very close to unity. This makes the alkali systems most attractive for the study of the essential phenomena of BEC, especially since the condensate density and temperature, and perhaps even the intrinsic atomic interactions [7], are subject to external experimental control.
However, the alkali condensates are always produced in particle traps, and this endows them with an essentially inhomogeneous character quite distinct from the homogensous nature of 4 He II. Detailed understanding of the properties of these systems requires [8][9][10][11][12][13][14][15][16] that this inhomogeneity be taken into account. The purpose of this paper is to present simple analytical methods which have been found to give results that are in good qualitative agreement with those of extensive numerical calculations based on mean-field theory, i.e., the basis-set methods of Refs. [12,16]. These methods are thus expected to be useful to experimentalists who wish to explore alternative trap designs, or to modellers who wish to analyze large-scale numerical results.

Overview of Dilute Gas Bose-Einstein Condensates
For 87 Rb and 23 Na the scattering lengths, a , are known to be positive [3,17,18]. In standard mean-field theory, this leads to a stable condensate wave function which is larger than the noninteracting ground state of the trap, and which can be modeled well by a mean field equation [8][9][10][11][12][13][14][15][16]. For 7 Li, a is known to be negative [2], which produces a net effective attractive interaction between the particles. A homogeneous condensate with a < 0 is predicted to be unstable [19]. However, a spatially confined condensate has been shown to be metastable [9]; and, when the average number of atoms in the condensate, N 0 , is sufficiently small, the condensate lifetime can be quite long. The potential energy of the confined system acts as a barrier to the collapse of the condensate [9,11]. Other discussions of the negative scattering length case have addressed the general issues of energetic stability and the possibility of a transition to a denser phase [20][21][22]. Theoretical predictions of the sizes and lifetimes of Bose condensed gases can now be compared directly with experiments, as will be discussed below. The effect of vortex formation [11][12][13][14] may also be examined in the present framework.
In this paper, I examine the general cases of the two approximate methods which have been most used to study Bose condensed systems, i.e., the Thomas-Fermi model [8,11], and the variational method [11]. While the ground state properties of such Bose condensed gases have been studied in some detail [8][9][10][11][12][13][14][15][16], most of the methods used have been computationally intensive. I show here that approximate methods give good quantitative results for the ground state properties of these gases, and also for vortex solutions of the equations of motion.

Basic Results of Mean-Field Theory
For the low atomic energies and densities achieved in the experiments of Refs. [1][2][3] the basic structure of a BEC should be well approximated by mean field theory, i.e., the Ginzburg-Pitaevskii-Gross (GPG) energy functional [23], described in Eq. (1). In such cases the atom-atom interaction is dominated by the effects of s -wave collisions. The interaction potential may thus be encapsulated in the form V (r , r ') = U 0 ␦ (r -r '), where U 0 = 4ប 2 a /M . Solutions of the resulting many-particle mean-field Schrödinger equation provide information about the size, density, chemical potential, and lifetime of the condensate as a function of N 0 . The evolution of a condensate, after release from the weak trap of Ref. [1] for time of flight studies, has been described using the time-dependent solutions of the Gross-Pitaevskii (GP) equation in Ref. [10]. Here I consider the condensates formed in the tight trap [1], which should be adequately represented by a time-independent formalism. The Gross-Pitaevskii-Gross (GPG) energy functional [23] for this system may be written where (r ) is the common wavefunction for each atom in the condensate, and [ ] is the energy per particle in the system. Minimizing [ ] with respect to variations in the wavefunction (r ), while enforcing the normalization condition, ͐dr | (r )| 2 = 1, gives the timeindependent GPG equation [23].
where V trap (r ) is usually an anisotropic harmonic oscillator potential [1,2]. Baym and Pethick [11] have identified the dimensionless parameter for a cylindrically symmetric harmonic system as where ␣ trap = (M /ប) 1/2 , and V trap (r ) = M( 2 2 + z 2 z 2 )/2. I use Eq. (1) to model a condensate using a variational trial wave function, and in the limit of sufficiently large we will also examine the Thomas-Fermi model, using Eq. (2), where the kinetic energy of the condensate atoms is neglected.
An important property, of both current and future experimental interest, is the lifetime of the condensate with respect to population loss from the trapped state. Such losses derive from spin-flip collisions that force the spin-flipped atoms out of the trap. These collisions are considered here to be extrinsic to the kinetics of BEC formation, and are treated by a simple rateequation approach. The two-body population loss rate for a condensate is given by Eq. (20) in Ref. [24]: where n 0 (r ) and n' (r) are the respective condensate and noncondensate number densities, and ␣ is the two-body loss rate coefficient associated with spin-flip collisions. Since, by the use of forced evaporative cooling, the noncondensate atoms may be "cut" away [1], forcing n '(r )~0, Eq. [4] reduces to Finding a solution of the GP equation, or minimum of the GPG energy functional, together with the loss-rate expression, enables one to estimate the condensate lifetime in units of ␣ [16,25]. While the ground-state properties of BECs are of primary interest, the possibility of vortex formation has also been considered [11][12][13][14]. A steady-state vortex wave function may be written as [26] (r ) = 1/2 (r ) e iS(r) , where (r ) is the density of the system. Here, S (r ), the velocity potential, changes by 2m along any closed loop around the vortex, with m being the vortex winding number. For the cylindrically symmetric harmonic traps of Refs. [1,2], a wavefunction of the form in Eq. (6), with the z -axis as the vortex line, is related to states of definite z -component of the total angular momentum. The corresponding eigenfunction may be written as where {n , n z , m } are the radial, axial, and angular momentum (z -component) quantum numbers. Equation (7) is equivalent to a winding number m vortex state. However, for an interacting system, the function f n n z m (r ) will also be a function of N 0 . As pointed out in Refs. [11][12][13][14], a rotating system should cause the formation of a vortex line or ring with critical rotation frequency, crit . In a homogeneous rotating system, the critical frequency is given by [26] where R is the radius of the system, and the radius of the vortex core. The critical frequency may be calculated directly from the equation ⌬E -L и = 0, where ⌬E is the energy gap between ground and vortex states, L is the vortex angular momentum, and is the rotational angular velocity of the system [26]. In the case of a cylindrically symmetric harmonic system, L = N 0 mប , and ⌬E = N 0 ( m -0 ), where m is the energy per particle in the m th vortex state, and it is a function of N 0 . This gives a critical frequency of Solving the GPG equation, and finding the energy per particle of the condensate, one may calculate the critical frequency for formation of such vortices in the system.

Simplified Variational Model
For a noninteracting cylindrically symmetric harmonic system, the eigenfunctions of the Schrödinger where ␣ ,z = (M⍀ , z /ប) 1/2 , H n (x ) is the n-th order Hermite polynomial, and L n (m) (x ) is the m-th associated Laguerre polynomial of order n [27]. In the noninteracting case (i.e., U 0 = 0) ⍀ = and ⍀ z = z , however, using Eqs. (10) as variational trial wavefunctions for the interacting system, with ⍀ and ⍀ z treated as arbitrary variational parameters, we find the GPG energy functional, where ␤ n n z m = 2ͱ 2 is a dimensionless constant depending on the state being considered, normalized so ␤ 0 00 = 1, and is important for analysis of vortex states. While Eq. (11) is simple to analyze numerically, as I show later, it may also be used to give approximate solutions in the limit of >> 1. While it is easy to examine the perturbative limit, | l << 1, it is of little experimental interest, since | | Ն 1 in all experiments to date [1].
I minimize Eq. (11) with respect to the variational parameters, ⍀ and ⍀ z , to find, given ⍀ = /⌬ , found from Ѩ /Ѩ⍀ = 0, where one has When N 0 is sufficiently large, ⍀ << and ⍀ z << z , so one needs only to examine terms of highest-order in and z . This is directly equivalent to neglecting the kinetic-energy terms in Eq. (11). This approximation then gives These expressions may be substituted back into Eq. (11) in order to evaluate the energy per particle of the system, and one may then calculate the critical frequency for vortex formation. For a > 0, the lowest vortex critical frequency occurs for the ground state to {n = 0, n z = 0, m = 1} transition, giving where 1/e is the radius of the 1/e population density of the ground-state condensate in the z = 0 plane. Unlike the homogeneous case of 4 He II, where crit~l n(R / )/R 2 from Eq. (8), the trapped condensate has crit~1 /R 2 in the large N 0 limit, where R = 1/e is a characteristic size of the condensate in the plane of rotation. For a < 0, the kinetic energy is never negligible, and becomes more important as the condensate population approaches the point of collapse. Since the kinetic energy may not be neglected in comparison with the selfinteraction energy, and the potential energy is of roughly the same order as the kinetic, the shrinkage of the condensate spatial distribution is not severe enough to make the potential energy much smaller than the kinetic energy. Hence, Eq. (14) must be solved numerically. Figure 1 shows plots, for the parameters appropriate to the experiment of Ref. [2], of the vortex critical rotation frequency calculated from numerical solution of Eq. (14) for transitions, m = 0 to m = 1 (long-dashed line), m = 0 to m = 2 (medium-dashed line), m = 0 to m = 3 (short-dashed line). Also shown is the frequency for an m = 0 to m = 1 transition (solid line), calculated by the basis-set method. As can be seen, the lowest critical frequency, and therefore the first transition to occur when the system is rotated, should be the m = 0 to m = ϱ transition. The maximum stable population for a condensate with negative scattering length is given by ⌬ = 0 from Eq. (15). Assuming ⍀ z ≈ z , this would give = ͭ 2n + m + 1 ␤ n n z m ͮ N max 000 .   Equation (21) gives N max 000~3 000 for parameters appropriate to the experiment of Ref. [2]. However, one finds N max 000~1 500 by numerical solution of Eq. (14), close to the value found in Refs. [9,[12][13][14]. Figure 2 shows the maximum population of a vortex state in the trap of Ref. [2], comparing the results of Eq. (22) (solid line), numerical solution of Eq. (14) (squares), and results from Ref. [13] (circles). As N 0 approaches the maximum stable level, many of the properties exhibit extreme behavior, i.e., the two-body loss rate has been shown to increase rapidly as the condensate population approaches the point of collapse [16].
The two-body loss rate is given, subtituting Eq. (10) into Eq. (5) by This shows, because R (2) n n z m (N 0 ) < R (2) n n z (m + 1) (N 0 ), that vortex states have longer lifetimes than those states without, or with less, angular momentum. Figure 3 shows a plot of the loss rate from a 7 Li condensate in the experimental configuration of Ref. [2], the total loss rate [including three-body recombination (which was found to be negligible except very close to the point of collapse)] from Ref. [16] (solid line), and the two-body loss rate calculated from Eq. (23) (broken line). While the variational approach clearly does not show as rapid a growth in the loss rate as the critical population, it does demonstrate qualitatively the same behavior.

Thomas-Fermi Model
For a large condensate, >> 1, with a > 0 we find meaning the kinetic energy term in Eq. (2) may be neglected [8,11]  Journal of Research of the National Institute of Standards and Technology condensate with a < 0, since the kinetic energy is the stabilizing factor in preventing the collapse of the condensate to a more dense phase [9,11]. For a > 0, however, this condition may be satified for a fairly small condensate population [8]. Previously, this method has been used to study ground state condensates in the large N 0 limit [8,11], but the possibility of studying vortex states has only been considered peripherally [12]. While Eq. (24) may still be true in the case of a vortex, we may not neglect the effect of the angular momentum on the spatial distribution of the condensate, since the angular momentum imposes a "centrifugal barrier." This reduces Eq. (2) to where the kinetic energy not related to the angular momentum has been neglected, and the angular dependence, i.e., the phase factor (which does not effect the density distribution), of the wavefunction has been ignored. For a cylindrical system one has ␥ 2 = m 2 , and for a spherical system one has ␥ 2 = l (l+1), where l is the total angular momentum quantum number. This gives a Thomas-Fermi wave function of the form where V eff (r ) = V trap (r ) + ប 2 ␥ 2 / 2M 2 is the effective trap potential. The relationship between N 0 and is found by the conventional normalization condition, ͐dr | | 2 = 1. Figure 4 shows a plot of the number density of a 2000 atom 87 Rb condensate for parameters appropriate to Ref. [1], the wave function calculated by basis-set expansion as in Ref. [12,16] (solid line), Thomas-Fermi (brokenline), and variational (dotted line) methods. In the simplest case, a harmonic trap with ␥ = 0, many of the properties of the condensate may be calculated explicitly [8,11], i.e., the two-body loss rate is given by where trap 3 = x y z = 2 z . Figure 5 shows the twobody loss rate from the ground state of the trap of Journal of Research of the National Institute of Standards and Technology to 200 s reported in Ref. [12]. The longer maximum lifetime of 280 s compared to 200 s occurs because I have neglected three-body recombination, while it was considered in Ref. [12]. The condensate has a lifetime, calculated by the basis-set method, because of spin-flip collisions in the range 39 s to 340 s. As we would expect the Thomas-Fermi model gives a shorter lifetime than the basis-set method, since the condensate is more tightly confined (within a finite region of space), c.f. Fig. 4. I have also calculated the vortex critical frequencies for parameters appropriate to the trap of Ref. [1] for comparison to the results obtained earlier for the variational model, and those presented in Refs. [12][13][14].  (14,15), to minimize the GPG energy functional), and 14 Hz (Eq. (20), the large N 0 approximation), compared to the 26 Hz reported in Ref. [13,14]. While one would not expect the Thomas-Fermi approximation to give accurate results for relatively small N 0 , as can be seen from Figs. (4 and 5), there is reasonable agreement with the basis-set calculations for as few as 2000 atoms.

Conclusions
I have shown that the variational method, using the wave functions of Eq. (10), gives very good agreement with the results found by much more computationally intensive methods [8][9][10][12][13][14]16]. For a > 0, it is clear that the variational and Thomas-Fermi methods may be used to examine the full range of condensate populations, at least when ground states and vortex states are being considered, to provide good estimates of the condensate properties, i.e., lower bounds on lifetimes and peak densities. The Thomas-Fermi approximation also has the advantage of being analytically solvable for the most common cases of experimental interest, i.e., harmonic traps as in Refs. [1,2]. It may also be solved numerically with little difficultly for traps of more arbitrary geometries [3]. For a < 0, the variational methods do not give as good agreement to the numerical calculations [13,14,16], especially as the critical population levels are approached. While qualitative behavior is retained, i.e., maximum population levels are in good agreement, other properties of the condensate, e.g., twobody collisional loss rates, do not show the same extreme behavior near the critical population levels.
The critical frequency for vortex formation in a large condensate of repulsive atoms was found to scale as crit ϳ 1/R 2 in the limit of large N 0 , where R is the characteristic radius of the condensate in the plane of rotation, compared to a homogeneous system in which crit ϳ ln(R / )/R 2 . The formation of vortices, and scaling of the vortex critical frequency, should be examined experimentally in the near future. While the approximate methods used here to calculate the vortex critical frequencies give reasonable results in comparison to the methods of Refs. [12][13][14], the scaling behaviour at large N 0 appears to be very similar for all the methods. In the case of negative scattering length, the formation of vortices has been considered as a mechanism for forming condensates containing greater numbers of atoms [13,14]: an m = 1 vortex may contain ϳ 4000 atoms, compared to ϳ 1300 atoms for the ground state in the trap of Ref. [2]. I have shown, however, that the formation frequency for such vortices is lowest for the m = ϱ state, so the actual formation of such structures in experiments must be studied more fully. Indeed a time-dependent calculation may be required to examine the condensate response to rotation near the critical vortex formation frequency. It is possible that vortices may be formed during the evaporative cooling cycle, as high energy atoms are "cut" away, and low energy atoms are ejected because of the rapid increase in loss processes [16,[20][21][22] for lower energy states.
Acknowledgments I thank K. Burnett, M. Edwards, and C. W. Clark for helpful comments.