Zero-Temperature, Mean-Field Theory of Atomic Bose-Einstein Condensates

We review the application of zero-temperature, mean-field theory to current experimental atomic Bose-Einstein condensates. We assess the validity of the approximations made by comparing the mean-field results with a variety of experimental data.


Introduction
The recent reports of Bose-Einstein condensation (BEC) in weakly interacting trapped alkali gases [1][2][3] has confirmed a property of bosons first predicted in 1924 by Bose [4] for photons, and in 1925 by Einstein [5] for atoms. The production of such condensates has opened the possibility of a new generation of atomic physics experiments on meso-or macro-scopic assemblies of atoms in the same quantum state.
The pure BEC phenomenon is manifested in the quantum statistical mechanics of noninteracting Bose particles; when encountered in nature it is modified by the effects of particle interaction. Such effects are quite severe in the case of superfluid 4 He, which heretofore has been thought of as the canonical example of BEC: the strong interactions between 4 He atoms allow only about 10 % of the atoms to occupy the condensate state [6]. However, it has recently been shown [7] that on the periphery of 4 He droplets, the condensate fraction approaches unity; this is because the atomic density vanishes at the extremity of the droplet, so the atom-atom interaction energy also tends to zero. The recent alkali BECs are distinguished from liquid 4 He in that the alkali atomic density is subject to wide-ranging experimental control, and has been indeed taken into the regime in which a dilute gas approximation is applicable, and in which it is reasonable to think of all atoms as occupying the same quantum state, as in the original concept of BEC. However, as will be seen below, the condensate state is much different from that which would describe a noninteracting gas. It is modified by effects of interactions that are encapsulated in a mean-field description, which can be thought of as the Hartree approximation to the wavefunction of a system of Bose particles. The quantitative treatment of such systems is the subject of this paper.
Before proceeding, we wish to note three additional aspects of the alkali systems that contrast with the case of liquid 4 He, as these make the general approach to the problem somewhat different from those traditionally used to treat superfluid systems.
First, the alkalis are confined by an external potential (a magnetic field or combination of magnetic and light fields), so their density is inhomogeneous. Thus the alkali BECs cannot be described adequately by a spatially uniform condensate wavefunction such as that which is used to describe bulk liquid 4 He. Not only must quantitative modelling methods be modified to treat inhomogeneous vs homogeneous BECs, but there are qualitative differences as well: for negative scattering lengths, a small long-lived BEC can exist in the inhomogeneous case [8,9], but not in a homogeneous system [10].
Second, as discussed by Cornell [11], the alkali BECs are intrinsically metastable . The equilibrium state of a confined alkali system at sub-microKelvin temperatures is a solid. However, the time for recombination of the gas is very long in the dilute limit, and is at least of the order of seconds in the current experiments.
Third, as shown elsewhere in this Special Issue [12][13][14], the ultracold collision physics of alkali BECs is exceedingly complex. Although the effects of collisions can be encapsulated in a few parameters (scattering lengths), the quantitative determination of these parameters is quite difficult and remains an active area of research. The work described in this paper utilizes these parameters as basic input, and it should be kept in mind that their values are subject to significant uncertainties, in no case less than 10 %. This paper presents a partial review of work we have undertaken to date in the field of modeling the alkali BECs. This work has been based on a zero-temperature, mean-field formulation of the quantum mechanics of an externally confined system of weakly interacting Bose particles. Many of the results of this theory, such as condensate geometries, lifetimes, and excitation frequencies, can be directly compared with the data of current experiments. We shall use this comparison to assess the validity of applying mean-field theory (MFT) to the current crop of experimental Bose-Einstein condensates (BECs).
The zero-temperature MFT equations presented here were first derived by Bogoliubov [15] many years ago in order to study the superfluid 4 He. The system to which they apply is assumed to be a weakly interacting, dilute gas of identical bosons, which, as noted above, does not provide a good description of liquid helium. However, it seems to fit the conditions present in a system of magnetically trapped gas of neutral alkali atoms. We emphasize that the previous statement should not be taken to be true a priori , but rather must be subjected to stringent experimental tests. We present here an overview of the comparison of MFT predictions with experiment.
The plan of the paper is as follows. In Sec. 2 we present a derivation of the Gross-Piteavskii (GP) and Bogoliubov equations (which we have been calling here the ''MFT'' equations). As part of the discussion we shall attempt to provide a detailed description of all of the approximations made in arriving at the MFT equations. In Sec. 3 we present the results of solving these equations for cases where comparison with experiment is possible. Sections 4 and 5 describe the algorithms and numerical procedures we have used to obtain the results presented in this paper and in previous work cited therein. Actual solution of the MFT equations for cases of specific experimental interest is a subject that has developed quite recently, and we believe there is considerable scope for enhancement of computational efficiency over that attained in current practice. Such enhancements will certainly be needed to go beyond the zero-temperature MFT description of BEC. Thus the material in Secs. 4 and 5 is presented at a level of detail needed to document our approach for use by those who can improve upon it.

Mean-Field Theory: Approximations and Derivations
In this section we present a somewhat detailed derivation of the basic zero-temperature MFT equations. These equations consist of the Gross-Pitaevskii equation, which describes the properties of the condensed part of the trapped atomic cloud, and the Bogoliubov equations, which describe properties of the non-condensed part. We shall present two derivations of the MFT equations. The first derivation uses a Bogoliubov transformation to cast the grand-canonical hamiltonian for a collection of interacting bosons into the form of a collection of noninteracting quasi-particles with the condensate becoming the vacuum state. The second derivation uses linear-response theory [16] performed on the time-dependent Gross-Pitaevskii equation (which is itself derived from a variational principle) to obtain the basic MFT equations. Before presenting these derivations, we shall first discuss the fundamental approximations made in modeling a cloud of cold, trapped atoms.

Fundamental Approximations
In the current generation of BEC experiments [1][2][3], a cloud of alkali atoms is optically pre-cooled and then magnetically trapped and evaporatively cooled to very low temperatures. The first major approximation leading to the MFT description is that the internal states of the atoms are ignored. All of the atoms must, however, reside in a particular hyperfine atomic ground state in order to remain trapped. The direction of the magnetic moment associated with the atom's internal state has been polarized to lie along the direction of the trap magnetic field at the site of the atom. Since the atoms are very cold and thus slowly moving, we assume that the magnetic moment of the atom adiabatically follows the local magnetic field [17]. Thus the energy of the interaction of the atom's magnetic moment ( atom ) with the external magnetic field has the form V trap (r ) = atom B (r ). (1) Another feature of this assumption is that collisions between atoms in the cloud do not change the atom's internal state. That is, all collisions are assumed to be elastic . In fact, most inelastic (spin-flip) binary collisions will cause both atoms to be ejected from the trap. This, in turn, limits the lifetime of the condensate. Such lifetimes can be predicted within MFT in a reasonably accurate way for comparison with experiment. Such comparisons are presented below. The true interaction potential between atoms in the cloud is quite complex. See, in this regard, Refs. [12][13][14] in this Special Issue. Most of this complexity is evident, however, only when the atoms are in close proximity. At the low temperature and density conditions present in the trap, all scattering events occur at extremely low energy. Consequently, the atoms rarely come close enough to each other to sample the complex nature of the inter-atomic potential. The atom-atom interaction is therefore well characterized by the s -wave scattering length, and the interaction potential may be written in the form: (2) where U 0 = 4ប 2 a/M, a is the s -wave triplet scattering length, and M is the atomic mass.
In the next section we present the derivation of the MFT equations using the Bogoliubov prescription which begins with the assumption that the atomic cloud can be approximated by a restricted grand-canonical ensemble.

Bogoliubov Prescription
Consider the many-atom system whose temperature is well below the condensation point and which is composed of a condensate plus thermal atoms. The grand canonical, many-atom hamiltonian, K = Ĥ Ϫ N where Ĥ is the many-body hamiltonian and N is the number operator, is written in terms of the field operator as follows: Ϫ ͵dr † (r) (r) where H 0 is the bare-trap hamiltonian, is the chemical potential, and V trap (r ) is the trap potential.
The boson field operators † (r ) and (r), respectively create and destroy an atom at position r and satisfy the commutation relations.
Under the Bogoliubov approximation, the condensate is assumed to contain most of the atoms so that N Ϫ N 0 << N 0 , where N 0 denotes the macroscopic occupation of the condensate and N denotes the total number of condensate plus thermal atoms. In this case, the field operator can be written as the sum of a c -number condensate wave function, ⌿ (r ), plus a small correction, (r ), where ⌿ (r ) satisfies the normalization condition ͵dr |⌿ (r )| 2 = N 0 .
Inserting Eq. (6) into Eq. (3) and neglecting terms in (r ) higher than quadratic yields the following expression for K .
The first term in the above equation is a c -number and the second and third terms will vanish identically if ⌿ (r ) satisfies the GP equation [18] [H 0 + U 0 |⌿ (r )| 2 ]⌿ (r ) = ⌿ (r ).
The Bogoliubov-approximate grand canonical hamiltonian [15], K B , then takes the form where ' is a c -number. The Bogoliubov hamiltonian is a sum of a quadratic form and a c -number and can be cast into the form of a collection of noninteracting quasi-particles by the following Bogoliubov transformation [19] (r ) = (u (r)␤ + * (r)␤ † ) (10) and † (r ) = (u * (r)␤ † + (r)␤ ) (11) where the ␤ are quasi-particle creation and destruction operators and the implicit assumption is made that the condensate wave function is not included in the sum. The quasi-particle operators satisfy the usual commutation relations for boson creation and destruction operators The reduction of K B to a collection of noninteracting quasi-particles occurs if the u and satisfy the following equations (after setting ⌿ (r ) = N 0 and where and the u and are square-integrable functions. The final form (to within a c -number) of K B is This hamiltonian has the form of a collection of noninteracting quasi-particles for which the condensate is the vacuum. Complete details of the derivation of the final form of K B are given in Ref. [19]

Linear-Response Theory
The time-dependent Gross-Pitaevskii equation can be derived from an action principle if a ''boson coherent state'' is used as the trial wave function [20], The factor |⌽ n (t )͘ is the usual Hartree many-body trial wave function-an n -fold product of one single-particle orbital ⌿ (r ,t ). The time-dependent GP equation describes the evolution of this orbital, The basic MFT equations can be obtained by a standard linear-response analysis of this equation [16]. To this end, we consider the effect of adding a weak, sinusoidal perturbation to the trap potential. The time-dependent GP equation then has the form The f Ϯ (r ) are the (possibly spatially dependent) amplitudes of the sinusoidal perturbation and p is the probe frequency.
To find the linear response of the condensate to the driving field, we shall assume that ⌿ (r ,t ) takes the form of a sum of an undisturbed ground-state part and a response part that oscillates at frequencies Ϯ p : Here, is interpreted as the chemical potential of the undisturbed ground state, and the condensate wavefunction is represented by the (scaled) condensate orbital g (r ). The functions u (r ) and (r) are the components of the condensate's linear response to the external disturbance that oscillate at frequencies Ϯ p .
After inserting Eq. (20) into Eq. (19), retaining only terms up to first-order in u (r ), (r ), and f Ϯ and equating like powers of e Ϯi p t , there result three equations that must be simultaneously solved for g (r ), u (r ), (r ), and . These equations describe the linear response of the condensate to weak external perturbation and have the following form To make the final connection with the basic equations of MFT, we show that Eqs. (22) and (23) can be solved by writing their solution as an expansion in the condensate normal modes. The equations that determine these modes are identical to the Bogoliubov equations. We now discuss there solution.
To find the normal modes of the condensate, we first set f Ϯ (r ) to zero in Eqs. (22) and (23). It is clear that the resulting equations will support square-integrable solutions only for discrete values of p , (we shall label them as ). The normal-mode equations thus have the form and where represents a set of quantum numbers. These equations are identical to Eqs. (13) and (14) if E = ប . To complete the connection between the quasi-particle excitation spectrum and the condensate response we now show how these normal modes describe the condensate linear response.
We define a normal mode as the following two-component object: With this definition, Eqs. (24) and (25) can be cast in the form, Where [20] orthonormal set where the scalar product of two normal modes is defined by (29) and the † denotes the transposed, complex-conjugated matrix.
The linear response equations, Eqs. (22) and (23), can be written, using this notation, as where The solution of the linear response equations is found by expanding both (r ) and g(r) in the normal modes Where the g are given by the following overlap integral Substituting these expansions into Eq. (30) yields a system of completely uncoupled equations to be solved for the c . The final solution is written as Note that the linear response diverges when the condensate is driven exactly on resonance. This unphysical behavior results from our neglect of loss processes and nonlinear effects.

Mean-Field Theory: Comparison with Experiment
Mean-field theory provides predictions for a variety of measurable condensate properties. These properties include condensate geometries, densities, and excitation frequencies. In this section, we shall discuss the comparison of the predicted values of these quantities with experiment. Measurements of condensate properties have been, to date, primarily performed on condensates formed in the traps of Refs. [1] and [3]. Comparisons with properties of condensates formed in the trap of Ref. [2] are not possible at present because measurements of these properties still have substantial uncertainties.
The geometry, density, and excitation frequency predictions of the MFT equations do not contain any adjustable parameters. The numerical constants that are input into the theory are the atomic mass (M ), the radial and axial trap frequencies ( and z ), the number of condensate atoms (N 0 ), and the scattering length (a ). All of these numbers are determined experimentally. The scattering length, in particular, is determined principally by photoassociation spectroscopic measurements [14,21].
Condensate lifetimes, on the other hand, depend critically on two-and three-body scattering event rates. These rates are quite difficult to determine accurately [13]. Therefore we will not consider lifetimes in the comparison of MFT with experiment.

1. Geometries and Densities
3.1.1 Geometries and densities were among the first condensate properties to be measured (see e.g., Ref. [1]). It is important to note that both of these 87 Rb condensate properties are determined indirectly. That is, after condensate formation, the trap potential is dramatically lowered, allowing the condensate to undergo a ballistic expansion, and then an absorption picture is taken. This picture, which essentially exhibits the velocity distribution of the original condensate, can then be used to extrapolate back to the original spatial density distribution.
The most quantitative comparison of theoretical and experimental geometry and density performed to date is that of Holland and Cooper (Ref. [22]). In this work, the time-dependent GP equation was solved by direct numerical integration. The time variations of the trap potential that occurred in the experiment of Ref. [1] were modeled and a prediction for the absorption picture was obtained. The agreement between theory and experiment for the velocity distribution appears to be at the 5 % level.
To date, no measurements of condensate geometries and densities of the 87 Rb condensate have been obtained in situ . However, it is interesting to see what the MFT predictions are for this case. Figure 1 shows the density profile for the 87 Rb condensate for the trap parameters of Ref. [1] for the case of N 0 = 2012 atoms. A plot of peak density as a function of condensate population is shown in Fig. 2. The peak density for this case is 5. 2 ϫ 10 13 atoms/cm 3 . The extrapolated experimental result is 3 ϫ 10 13 atoms /cm 3 but this number is accompanied by substantial error bars within which the MFT number falls.   Journal of Research of the National Institute of Standards and Technology

23 Na
To date, condensates have been formed in two different traps at MIT. In the optical-plug trap, atoms were prevented from undergoing majorana spin flips through the use of a blue laser focussed at the center of the trap [3]. The "cloverleaf" trap is an Ioffe-Pritchard trap and thus there is no zero of the magnetic field [23]. The condensates formed in each of these traps are very different from the 87 Rb one. The major difference is size-the 23 Na condensates contain on the order of 10 6 atoms.
Mean-field theory solutions for condensates of this size can be obtained via the "Thomas-Fermi" approximation [24]. The approximation amounts to the neglect of the kinetic energy term in the Gross-Pitaevskii equation. This reduces Eq. (8) to where the kinetic energy term, -(ប 2 / 2M) ٌ 2 (r), has been neglected. The wavefunction thus has the form The relationship between N 0 and is found by the normalization condition ͵dr͉(r)͉ 2 = 1.

(37)
A comparison of the Thomas-Fermi approximate solution with the basis-set solution is shown in Fig. 3. The comparison of the geometries and densities of the 23 Na condensates with the results of MFT suffers from the uncertainty in the value of the sodium scattering length. In contrast to the 87 Rb scattering length, until recently the 23 Na scattering length was known only to within a factor of two [14]. The MFT solutions depend only on the parameter N 0 a /l sho where l sho is the harmonic oscillator length scale. Uncertainty in a will lead to a similar uncertainty in MFT predictions.

Excitations
The excitation data provide the opportunity for the most quantitative comparison with MFT to date. Excitation frequencies are measured by forming a condensate, driving it weakly by oscillating the trap potential, waiting a specified delay time, and then probing the condensate. This cycle was repeated for increasing delay times forming a time history of condensate oscillations. The width of the condensate was observed to oscillate and the frequency of this oscillation measured. This experiment was recently performed both at JILA [25] and at MIT [26]. The comparison [27] with the JILA results is shown in Fig. 4. The agreement varies between 2 % and 5 %. Similar agreement with the theory [28] was found in the MIT experiment.

Summary
This paper has been written in the early days of quantitative modelling of dilute atomic BECs, and there have yet been relatively few stringent tests of the validity of MFT. However, it has been found to have good predictive and interpretive value. As the accuracy of the experiments improves, and the uncertainties in the values of the microscopic parameters are reduced, we expect that dilute atomic BECs will provide a new testing ground for the Bogoliubov approximation and its variants, which are among the cornerstones of the quantum theory of many-particle systems.

Solving the Gross-Pitaevskii Equation
The time-independent GP equation is written We expand (r ) in the basis of trap-potential eigenfunctions, i.e., where the sum extends over the N basis basis-set functions, represents a set of quantum numbers, and (r ) is an eigenfunction of the trap Hamiltonian, Inserting Eq. (39) into Eq. (38) and taking the scalar product with each basis-set function converts the problem to one of finding the simultaneous root of N basis nonlinear algebraic functions having the following form The coefficients C ( 1 , 2 , 3 , 4 ) are the integrals of the product of four trap wavefunctions; several methods for evaluating these integrals are described in Appendix 5. These nonlinear equations can be solved by a standard Newton-Raphson [29] technique provided a good guess at the solution is known. Equations (42) are not sufficient to determine the set {a } since is also an unknown. If we choose a value of , and a corresponding value of N 0 , a solution may be found that will, in general, not be normalized. This solution may be normalized by the simultaneous rescal-ing of the basis-set coefficients {a }, and the number of condensate atoms, N 0 .
For example, find an arbitrary solution N A (r ) of the time independent GP equation then find the norm of the solution where A is some constant. The set {a A }, the basis-set coefficients for the solution N A (r ), and the condensate population N A may be rescaled as follows to give a normalized solution N B (r ) for the new GP equation So we find an arbitrary solution N A (r ), then find a normalized solution N B (r ) from N A (r ) by rescaling both the wavefunction and the number of atoms in the condensate. This method is also discussed in Refs. [30,31].
To ensure that the final normalized solution is actually the lowest energy state of the condensate, we simulate the adiabatic growth of the nonlinear term in Eq. (38) by first finding a solution for a value of close to the noninteracting ground state energy, which corresponds to a very small number of particles in the condensate, and so the condensate wavefunction is very close to that of the trap ground state. The initial values of the basis-set coefficients are those of the trap ground state, i.e., {a} ={1,0,0, . . .}. The Newton-Raphson algorithm is then used to solve Eqs. (42). The chemical potential is then incremented, corresponding to a increase in condensate size, and the basis-set coefficients found previously, for the smaller condensate, are used as the initial values for this new solution. This process, which is repeated for increasing values of , gives not only the lowest energy solution, but also generates solutions for a large range of , and corresponding N 0 , with little wasted computational time.

Solving the Bogoliubov Equations
In this section we describe how Eqs. (13) and (14) are numerically solved using a finite basis consisting of trap eigenfunctions. The solution is found by first expanding u (r ) and (r), in the eigenstates of H 0 to produce a set of algebraic equations. These equations are then cast into a generalized matrix eigenvalue equation having the form Ax = w Bx . In this equation, A and B are square matrices of order 2N basis (where N basis denotes the number of basis-set states used), x is an 2N basis -component column vector containing the coefficients in the basisset expansion of u (r ) and (r), and w is the generalized eigenvalue. We shall treat the case of a spherically symmetric harmonic trap potential for simplicity. It is important when solving the normal-mode equations using a finite basis that Eq. (21) for g (r ) also be solved within this set. The where = ␤ប trap and = ␣ trap . The full set of equations that must be solved to find the ␣ consists of those given just above for all values of l and m . The matrix element w ss' (lm) is given by where 2s+l, l, m (r ) = u sl (r )Y lm (r)/r defines the reduced radial wavefunction u sl (r ). The spherical symmetry of the condensate ground state imposes a particularly simple structure upon the equations for a slm () and b slm () . These equations divide into subsets of equations that are self-contained within the subspace of fixed angular momentum and magnetic quantum number, i.e., there is no coupling between a slm The quantum number n is defined as the number of nodes found in the radial wavefunctions, U nlm (r ) and V nlm (r ). Furthermore, n runs from -ϱ to ϱ because the normal frequencies occur in pairs of opposite sign. That this must be the case can be understood by examining Eqs. (24) and (25). If → -and if u and are exchanged, then Eqs. (24) and (25) remain unchanged. Hence, if there exists a solution of the normal-mode equations for positive , then a negative frequency solution also exists having the roles of u and reversed. Indeed, we find that normal frequencies occur in equalmagnitude pairs of opposite sign in our solutions. A negative value of n , therefore, labels a normal-mode solution obtained by performing the above transformation on the corresponding positive-n solution. The n = l = m = 0 normal-mode solution deserves special comment. This solution always has a frequency of zero (i.e., 000 = 0) with u 000 (r ) and 000 (r ) being proportional to g (r ). This condensate-mode solution has zero norm, see Eq. (29).
Moving the terms in Eqs. (50) and (51) containing ␣ to the right-hand-side gives where I 2 is the (2 ϫ 2) unit matrix and where s runs from 0 to N basis -1.
The final transition to the form of a generalized eigenvalue problem is accomplished by defining the following (N basis ϫ N basis ) matrices and (w (lm) ) ss' = w ss' Using these matrices, Eq. (58) may be written as a single matrix equation where c (nlm) is given by . c (nlm) = . .
and where the matrix A (lm) can be expressed as a tensor product of matrices: and B is given by where I N basis is the unit matrix of order N basis . The square matrix on the left of the ᮏ symbol is of order N basis and the matrix on its right is of order 2. The full matrix is of order 2N basis with elements that can be indexed as C (s,k)(s'k') = A ss' B kk' , where s runs from 0 to N basis -1 and k runs over 0, 1.

Appendix B. Integration Methods
The basis-set expansion method requires the evaluation of a set of integrals, the integrals of four basis-set functions over all space. In general the basis-set functions would be those for an anisotropic trapping potential, though in two of the cases of interest [1,2] they are eigenfunctions of the cylindrically symmetric harmonic oscillator. We will now consider how these integrals may be evaluated.

Cylindrical Symmetric Traps
The integral which must be evaluated is given by where the basis-set functions (r ), for a harmonic trapping potential are given by where = (x 2 + y 2 ) 1/2 and ␣ ,z = (M ,z /ប) 1/2 . These functions may be rewritten, with substitutions = (␣ ) 2 and = ␣ z z , to give the normalization constant remains unchanged, and the integral operator transforms, so the integral has been broken down into three onedimensional integrals, i.e., ϫU( 1 , 2 , 3 , 4 ) ⌶ (n z 1 , n z 2 , n z 3 , n z 4 )

Completely Anisotropic Traps
The integral to be evaluated for a completely anisotropic trap, is given by where the basis-set functions (r ) are now given by: where ␣ x, y, z = (M x, y, z /ប) 1/2 .These functions may be written, i (r ) = N i e -(␣ 2 x x 2 +␣ 2 y y 2 +␣ 2 z z 2 )/2 ϫ H n x i (␣ x x ) H n y i (␣ y y ) H n z i (␣ z z ).
With the substitutions ␣ x x = u , ␣ y y = and ␣ z z = w , this reduces Eq. (80) to 3 , n x 4 )⌶ (n y 1 , n y 2 , n y 3 , n y 4 ) ϫ ⌶ (n z 1 , n z 2 , n z 3 , n z 4 ) (84) we will now outline some of the methods we may use to evaluate these integrals.
Equation (78) may be evaluated by use of an addition formula for Hermite polynomials [32] to give n 3 !n 4 ! min (n3,n4) t = 0 2 t + s-1 ⌫ (s-n 1 ) ⌫ (s-n 2 ) t ! (n 3 -t)! (n 4 -t )! ϫ ⌫ (s-n 3 -n 4 +2t ), ⌶ (n 1 ,n 2 ,n 3 ,n 4 ) = if n 1 +n 2 +n 3 +n 4 is even; where s = (n 1 +n 2 +n 3 +n 4 -2t+1)/2. By expanding the Laguerre polynomials in Eq. (77) as finite series [33], an even more complicated expression can be derived for the cylindrically symmetric case. While these results are exact, they are of little practical use for computational calculations in standard extended precision arithmetic, since numerical accuracy is rapidly lost in the summations of terms which vary over several orders of magnitude. An alternative method for evaluating these integrals is that of gaussian quadratures which is described next.

Gaussian Quadrature Methods
For simplicity ⌶ (n 1 , n 2 , n 3 , n 4 ) = Using a gaussian quadrature integration technique [29] reduces this integral to a finite sum ⌶ (n 1 , n 2 , n 3 , n 4 ) = 1 ͙2 N () where g (h) i is the i th Gauss-Hermite quadrature point, and w (h) i is the corresponding weight. This sum is exactly equal to the integral provided N Quad () > 2 max(n z ) [29].
Applying the same technique to L (m 1 ) where g (l) i is the ith Gauss-Laguerre quadrature point, and w (l) i is the corresponding weight. This sum is exact provided N Quad () > 2 max(n )+max(m ) [29]. This means that the choice of basis-set functions fixes the number of quadrature points needed, e.g., for initial calculations of a cylindrically symmetric ground state condensate, we may use n = {0, . . . , 6}, m = 0, and n z = {0, . . . , 6}, fixing N Quad () = 13, and N Quad () = 13.