The Condensate Wave Function of a Trapped Atomic Gas

We discuss various properties of the ground state of a Bose-condensed dilute gas confined by an external potential. We devote particular attention to the role played by the interaction in determining the kinetic energy of the system and the aspect ratio of the velocity distribution. The structure of the wave function near the classical turning point is discussed and the drawback of the Thomas-Fermi approximation is explicitly pointed out. We consider also states with quantized vorticity and calculate the critical angular velocity for the production of vortices. The presence of vortex states is found to increases the stability of the condensate in the case of attractive interactions.


Introduction
In this paper we discuss some relevant properties of the ground state of a Bose-condensed atomic gas confined by an external potential. Our starting point is the Gross-Pitaevskii equation which gives the proper Schrödinger equation for the order parameter of an inhomogeneous dilute Bose-condensed gas at zero temperature. Using this equation it is possible to discuss various ground state properties of the system: the form of the atomic cloud, the role of the interatomic potential, the velocity distribution, and so on. All these quantities are essential for the interpretation of the recent experiments on Bose-Einstein condensation in ultracold alkali atom gases [1][2][3]. An important feature is the profound difference between systems interacting with repulsive and attractive forces. In the latter case, in particular, the stationary solution given by the Gross-Pitaevskii equa-tion is of metastable type. If the number of atoms is too large such a solution becomes unstable and the system collapses. In this paper we will also discuss some rotational properties of the system, in particular the structure of vortices and the critical angular frequency needed to generate a rotational instability.
The Gross-Pitaevskii equation for the order parameter (r ) = ͗(r)͘ has the well known form [4]: where V ext is the external confining potential, which is usually chosen in the form of an anisotropic harmonic well. The role of interactions is accounted for by the non-linear term and is parametrized by the s -wave scattering length a . The quantity is the chemical potential and is fixed by imposing the proper normalization, N =͐ dr , to the density of the system = ͉ ͉ 2 . The Gross-Pitaevskii equation ignores interaction effects due to the atoms outside the condensate. This is an excellent approximation for a dilute Bose gas at low temperature where the depletion of the condensate is negligible. An important question to discuss concerning the ground state of a trapped Bose gas is the role of the interatomic potential. At first sight one would in fact expect that the role of interactions be negligible in a dilute system, where the usual expansion parameter a 3 is extremely small. Actually it turns out that the interaction can have a deep influence on the solution of Eq. (1), where its effect turns out to be fixed by the adimensional parameter Na /a HO , where a HO = ͙ប /(m HO ) is the harmonic oscillator length. This parameter can be indeed rather large despite the smallness of a 3 . The final result is that the system is still fully Bose condensed, but the structure of its wave function can be strongly affected by the interatomic forces.
For the above reasons it is useful to discuss the solution of the Gross-Pitaevskii equation in two relevant limits: the noninteracting model and the strongly repulsive limit, which corresponds to the Thomas-Fermi approximation. We will make a comparison with the exact numerical solution of the Gross-Pitaevskii equation in order to point out the role of the interaction. We will devote special attention to the structure of the condensate wave function near the boundary, close to the classical turning point, and we will finally study the case of quantized vortices.

The Noninteracting Model
When the scattering length a vanishes, the problem reduces to the solution of a one-body Schrödinger equa-tion. Let us put a = 0 in Eq. (1) and take the external potential as an anisotropic harmonic oscillator: The ground state wave function becomes where = z / Ќ . The Gaussian has different transverse and vertical widths. In particular one has ͗x 2 ͘ = ͗y 2 ͘ = (1/2)a 2 Ќ and ͗z 2 ͘ = (1/2) Ϫ1 a 2 Ќ . The chemical potential is (1 + /2)ប Ќ and coincides with the energy per particle, while the kinetic energy per particle has the simple form An interesting quantity to discuss is the ratio ͙͗p 2 z ͘/͗p 2 x ͘ which provides a measure of the aspect ratio , characterizing the anisotropy of the velocity distribution. Using the wave function Eq. (3) one finds Values of the aspect ratio different from 1 reflect a peculiar and unique feature of Bose-Einstein condensation.

The Strongly Repulsive Limit (Thomas-Fermi Approximation)
The opposite limit is obtained when the interaction is so strong, or the number of particles so large, that the kinetic energy term ᭞ 2 can be neglected in the Gross-Pitaevskii equation, Eq. (1). It corresponds to very large values of the dimensionless parameter Also in this case the solution of the Gross-Pitaevskii equation is trivial and the wave function has the form: if the right hand side is positive, and = 0 elsewhere. The chemical potential is easily calculated by imposing the normalization condition ͐ dr = N . One finds = 1 Due to the different scaling properties of the wave function with respect to the variable z (compare Eqs. (3) and (7)), the aspect ratio ͙͗p 2 z ͘/͗p 2 x ͘ in this case is equal to differently from the noninteracting case Eq. (5).
The wave function Eq. (7) is expected to approximate well the exact solution of the Gross-Pitaevskii equation Eq. (1) for large N , apart from the structure of the surface region where the exact wave function has to vanish smoothly. Some relevant observables, as the kinetic energy, can be significantly affected by this surface structure, as we will see in Sec. 5.

Solution of the Gross-Pitaevskii Equation
The Gross-Pitaevskii equation, Eq.
(1) can be solved numerically [5,6]. We transform the differential equation in a functional minimization and use a steepest descent method to solve the minimization problem on a grid of points. As an example of atoms with repulsive interaction we choose 87 Rb, as in the experiment of Ref.
[1]. For the s -wave triplet-spin scattering length we use a = 100a 0 , where a 0 is the Bohr radius. The asymmetry parameter is taken = z / Ќ = ͙8 and the axial frequency z /2 = 220 Hz. The corresponding characteristic length is a Ќ = 1.222 ϫ 10 Ϫ4 cm.
Results for the chemical potential and the energy per particle are shown in Table 1. Both quantities are expressed in units of ប Ќ , or of the equivalent temperature ប Ќ /k B = 3.73 nK. The partial contributions to the energy per particle coming from the kinetic energy (kin), the harmonic oscillator potential (HO) and the internal potential energy (pot) are also given. The N = 1 case coincides with the noninteracting anisotropic harmonic oscillator: in this case the total energy per particle coincides with the analytic value (1 + /2) = 2.414. When N increases the repulsion among atoms tends to lower the central density, expanding the cloud of atoms towards regions where the trapping potential is higher. A typical profile of the condensate wave function is plotted along the x -axis for N = 5000 in Fig. 1. The exact minimization of the Gross-Pitaevskii functional (solid line) is compared with the noninteracting case (dashed line) and the Thomas-Fermi limit (dot-dashed). When N is large we observe an increase of both interaction and harmonic oscillator potential energy per particle (the latter effect follows from the expansion of the cloud). Conversely, the kinetic energy per particle decreases because the density distribution is flattened. In the strongly repulsive limit, N → ϱ, one should find that the internal potential energy is much greater than the kinetic energy. Indeed the convergence towards this limit turns out to be rather slow as we will show in the next section.
Another interesting quantity which can be easily calculated from the ground state wave function is the aspect ratio of the velocity distribution, that is the ratio ͙͗p 2 z ͘/ ͗p 2 x ͘. This quantity is equal to ͙ in the noninteracting case and should approach in the strongly repulsive limit. The numerical results, as a function of N , are shown in Fig. 2. The two limiting cases are shown as dashed lines. One clearly sees that the convergence to the value 2.828 = is very slow; the aspect ratio remains well below the asymptotic value even for N = 20000. The aspect ratio measured in Ref.
[1] is estimated to be about 50 % larger than the noninteracting value, while the number of particles is of the order of 5000. The agreement with our results is good, even if one has to consider that the experimental estimate implicitly assumes a ballistic expansion of the atoms after switching off the external trap. The effects of the interaction on the expansion of the gas should be explicitly taken into account in order to draw more definitive conclusions.
As an example of atoms with attractive interaction we choose 7 Li, as in the experiment of Ref. [2]. For the s -wave triplet-spin scattering length we use a = Ϫ 27a 0 . The axial frequency reported in Ref. [2] is z /2 = 117 Hz and the corresponding characteristic length is a Ќ = 2.972 ϫ 10 Ϫ4 cm. The transverse frequency is z / 2 = 163 Hz, so that the asymmetry parameter is = z / Ќ = 0.72. The first important point to stress is that Gross-Pitaevskii functional has no global minimum for negative scattering length. This reflects the tendency of the system to collapse. For spatially inhomogeneous systems, however, the zero-point energy can exceed the attractive potential, producing local minima of the functional when the density of atoms is not too high.
The most striking difference with respect to the repulsive case is that here the central density of the cloud increases rapidly with N , as shown in Fig. 3. This is the effect of adding more and more attractive potential energy. When the central density reaches a certain critical limit the system collapses and the solution of the Gross-Pitaevskii equation does not converge anymore. In 7 Li, with the input parameters given above, the critical number N turns out to be about 1400. In [6] we have found that a stationary solution of the Gross-Pitaevskii equation with larger values of N can be obtained if a vortex line is present in the system. The possible occurrence of vortices in these trapped Bose gases will be discussed in the Sec. 6.

Wave Function at the Boundary
As one clearly sees in Fig. 1, the Thomas-Fermi approximation fails to reproduce the structure of the order parameter at the surface of the atomic cloud in the case of positive scattering length (repulsive interaction). Several measurable quantities can be significantly affected by the behavior of the wave function in this region. In order to provide a good model for these quantities one has to go beyond the Thomas-Fermi approximation. One of these relevant observables is the kinetic energy In fact the Thomas-Fermi approximation (7) for the wave function is not appropriate to evaluate E kin ; it produces a logarithmic divergence in the integrand of Eq.
(9), occurring at the classical turning point, where V ext = . This reveals that the evaluation of E kin requires higher accuracy in the description of the boundary region. In order to provide the proper description of the condensate wave function near the boundary we have recently proposed [7] a suitable expansion of the Gross-Pitaevskii equation near the classical turning point. The resulting analysis allows one to obtain the proper expansion for the kinetic energy. We briefly sketch here the case of isotropic traps (a Ќ = a z = a HO ). Let R be the boundary of system spherical system, determined by the equation = V ext (R ). Near this point one can carry out the expansion where F = m 2 HO R is the modulus of the attractive external force evaluated at r = R . Close to the boundary, where |r Ϫ R | << R, the Gross-Pitaevskii equation takes the form Let us now introduce the dimensionless variable is a typical thickness of the boundary giving, as we will see later, the distance from the classical radius R where the Thomas-Fermi approximation starts failing. Then we introduce the dimensionless function defined by where the radius R is related to N by the equation

Vortices
The structure of vortices in a trapped Bose gas can be naturally investigated in the present formalism. Let us consider states having a vortex line along the z -axis and all the atoms flowing around it with quantized circulation. One can write the axially symmetric condensate wave function in the form where (r ) = ͙(r) is the modulus, while the phase S acts as a velocity potential: = (ប/m)᭞S. By choosing S = , where is the angle around the z -axis and is an integer, one has vortex states with tangential velocity with r 2 Ќ = x 2 + y 2 . The number is the quantum of circulation, and the angular momentum along z is Nប .
If the complex wave function Eq. (18) is used in the derivation of the Gross-Pitaevkii equation, one gets which differs from Eq. (1) only for the addition of a centrifugal potential. This new term forces the solution to vanish on the z -axis for 0. For noninteracting particles one falls again in the case of the stationary Schrödinger equation for the anisotropic harmonic potential. For instance the = 1 solution has the form (r ) ϰ r Ќ exp ͫϪ 1 To get the energy per particle for the 0 states one has simply to sum ប Ќ to the energy per particle of the ground state without vortices.
In the interacting case the kinetic energy can not be neglected even for large N , since it determines the structure of the vortex core. In particular, the balance between the kinetic energy and the interaction energy fixes a typical distance over which the condensate wave function can heal. For a dilute Bose gas the healing length is given by where is the density of the system. In the case of a vortex it corresponds to the distance over which the wave function increases from zero, on the vortex axis, to the bulk density. For the trapped atoms in the N → ϱ limit one finds R Thus the healing length is small compared with the size of the cloud if R is much bigger than a Ќ . The critical angular velocity required to produce vortex states is easily calculated once the energies of the states with and without vortices is known. One has to compare the energy of a vortex state in frame rotating with angular frequency ⍀ , that is (E Ϫ ⍀L z ), with the energy of the ground state with no vortices. Since the angular momentum per particle is ប , the critical angular velocity is In the noninteracting case the difference of energy per particle is simply ប Ќ , so that ⍀ c = Ќ . We have solved numerically the Gross-Pitaevskii equation (20) both for rubidium and lithium. In Fig. 6 we show the wave function of a cloud of 5000 87 Rb atoms; the = 1 wave function (Fig. 6b), which corresponds to atoms flowing around the z -axis with angular momentum Nប , is compared with the = 0 ground state (Fig. 6a). The atoms are pushed away from the axis forming a toroidal cloud. From the energy of the vortex states we calculate the critical angular velocity, through Eq. (24). The results for = 1 are shown in Fig. 7. The critical angular velocity decreases rapidly with N . For N > 5000 it is less than 40 % of the noninteracting value, given by the transverse angular frequency Ќ of the trap. The healing length is the distance over which the wave function grows from zero to the bulk value. In the limit of large systems it can be approximated by Eq. (22) with equal to the density in the central part of the toroidal distribution. Both the estimate of and ⍀ c Coming back to the question of the stability for negative scattering length, we notice that, when the local minimum associated with wave functions of the form shown in Fig. 3 disappears, nothing prevents a priori the existence of other local minima associated with different configurations. Such configurations should have local density lower than the critical one. A natural way to obtain a favourable situation is to move the atoms away from the z -axis, conserving the total number of particles. This happens in the presence of a vortex. In Fig. 8 we show the wave function for 1000 7 Li atoms with no vortices (Fig. 8a) and with an axial vortex of unit circulation (Fig. 8b).We use the same units in both cases, so one can see that the maximum value of the wave function inside the toroidal distribution of the vortex is approximately a factor two lower than the central value in the state with no vorticity (the density is four times smaller). The critical angular frequency for the formation of the vortex state in Fig. 8 is 1.12 times the transverse angular frequency of the trap. In systems with attractive interaction the critical angular velocity is larger than for noninteracting particles, while the opposite is true for repulsive interaction. This is because it costs internal potential energy to lower the average density, as the vortex does, for attractive interactions. However, once a vortex is created, the corresponding state is more stable than in the absence of vorticity: one can put more atoms inside the rotating cloud before reaching the critical density for the final collapse. Indeed we find local minima of the Gross-Pitaevskii functional for N much larger than 1400 if > 0. For = 1 we find a critical value of N Ӎ 4000; for = 2 and 3 we find critical values of 6500 and 8300, respectively. It is worth mentioning that the number of particles in the condensate reported in the experimental work of Ref. [2] is an order of magnitude higher than the critical value for the stability of the Gross-Pitaevskii solution without vorticity (N Ӎ 1400). The presence of vortices might explain the large size of the observed Bose-condensed gas. Further experimental data are needed to draw more definitive conclusions.