# A numerical approach to the ground and excited states of a Bose-Einstein condensed gas confined in a completely anisotropic trap

## I INTRODUCTION

The experimental achievement of Bose-Einstein condensation (BEC) in dilute alkali-metal gases confined in magnetic traps has generated tremendous interest in the behavior of the inhomogeneous, weakly interacting, and dilute Bose gas. At low temperatures, the confined Bose gases have been shown to be well-described by mean-field linear-response theories (MFLRT) based on the Bogoliubov approximation, which assumes that the number of condensate atoms $`N_0`$ is a substantial fraction of the total number . The finite-temperature extensions of this theory, such as the Hartree-Fock-Bogoliubov and Popov approximations, also have been successfully applied to these systems, though the microscopic basis for these approaches and their agreement with experiment remain somewhat uncertain.

The usual approach taken within MFLRT is to first solve the (time-independent) non-linear Schrödinger equation for a given number of atoms in the condensate; the resulting wavefunction and chemical potential are then used in the linear(ized) equations for the quasiparticle excitations in order to obtain the eigenmodes of the system. At finite temperatures, the significant depletion of the condensate requires that this procedure be iterated to self-consistency. The magnitude of the non-linear, or self-energy, term appearing in the Schrödinger equation for the ground state is proportional to the condensate density. Thus, a fully three-dimensional solution of the non-separable equation for the condensate presents a number of numerical challenges, particularly for large values of $`N_0`$.

The magnetic trap used to confine these gases is well approximated by a harmonic potential, in which the frequencies ($`\omega _x,\omega _y,\omega _z`$) are generally incommensurate. All of the published experimental realizations to date have employed traps with either spherical or cylindrical symmetry, which are more conducive to numerical study. The number of experimental groups which have obtained BEC of confined alkali-metal gases is steadily rising, however. In order to consider a range of possible trap geometries, as well as to permit future investigations of any time-dependent properties, it is necessary to generalize the numerical calculations to allow for complete anisotropy. The purpose of the present paper is to demonstrate that the numerical obstacles may be overcome and to outline a robust and computationally efficient procedure within the MFLRT.

The approach taken by us is based on a discrete variable approximation (DVR) to the equations governing the statics of the condensate. The DVR has several advantages over the methods used by others for the current problem. Most standard techniques discretize the solution to the partial differential equations (PDE) for the condensate and excitations either in physical (grid) space or in function (basis set) space. Grid-based methods  have the advantage that all local interactions between particles have a local representation while the kinetic energy takes on a sparse and typically banded structure. This sparsity is a consequence of (typically) low-order finite difference approximations to the derivatives in the PDE, and is crucial for the implementation of iterative techniques  developed for large, sparse linear systems. The disadvantage of grid methods is that it is often difficult to approximate the derivatives that appear in the PDE’s to sufficient accuracy without resorting to high-order differences or very small step sizes. While an accurate representation of the kinetic energy is arguably of minor importance for the ground state when $`N_0`$ is large, it is crucial for the subsequent determination of excitations. Expansions in function space (often called spectral and/or pseudospectral methods ) are typically superior since it is possible to “analytically” differentiate the functions without further approximation. These methods, however, require the evaluation of potential energy matrix elements by quadrature, leading to a non-diagonal and dense matrix representation of local operators.

The DVR exploits the dual relationship between certain orthogonal polynomials (such as the classical orthogonal polynomials) and the points and weights of a Gauss quadrature. By using an appropriate set of polynomials as a basis for expanding the solution to the PDE, it is possible to maintain almost all of the advantages of a grid-based approach as well as the global convergence of a finite basis set. In the DVR, any local potential energy operator is diagonal and therefore easy to compute. The multidimensional kinetic energy also has a sparse representation because it is a commuting sum of one-dimensional operators. The sparseness is not as structured as is the case for more traditional approaches, but the expense of generating the kinetic energy matrix elements is mitigated by the iterative method used to solve the PDE: these terms need only be evaluated once. It is also worth noting that the expansion coefficients of the solutions to the PDE are trivially related to the values of the solution on the quadrature points which serve as the physical grid. Perhaps the most important point to stress, though, is that the solutions of the PDE scale very efficiently with basis set size. This will become very important for three-dimensional problems where the size of the various matrices could easily be as large as 100 000 by 100 000.

The usual starting point for the theoretical treatment of the inhomogeneous, weakly interacting, dilute Bose gas at zero temperature is the time-independent MFLRT in the Bogoliubov approximation. Since the derivation of the resulting equations may be found in many places, only the central results are presented here. In the absence of inhomogeneities other than the confining potential, all wavefunctions may be assumed to be real. In the weakly interacting and low-density limit relevant to the experiments of interest, the interatomic interactions may be approximated by a two-body delta-function contact potential with a coupling constant $`g`$. Minimizing the grand canonical potential for interacting bosons and then linearizing the resultant equations for the amplitude of the condensate $`\psi (𝐫)`$ and the excitations, $`u(𝐫)`$ and $`v(𝐫)`$, leads respectively to the Gross-Pitaevskii (GP) for the condensate and the Bogoliubov equations for the excitations:

$$\widehat{L}\psi (𝐫)=\mu \psi (𝐫);$$
(1)   
$`\left(\widehat{L}\mu +V_\mathrm{H}\right)u_n(𝐫)V_\mathrm{H}v_n(𝐫)`$ $`=`$ $`ϵ_nu_n(𝐫);`$ (2)   
$`\left(\widehat{L}\mu +V_\mathrm{H}\right)v_n(𝐫)V_\mathrm{H}u_n(𝐫)`$ $`=`$ $`ϵ_nv_n(𝐫),`$ (3)   
where the non-linear Schrödinger operator is written as

$$\widehat{L}=\frac{\mathrm{}^2}{2M}\stackrel{}{}^2+V_{\mathrm{trap}}(𝐫)+V_\mathrm{H}.$$
(4)   
The non-linearity of the GP equation is due to the mean-field, or Hartree, potential $`V_\mathrm{H}=g|\psi (𝐫)|^2`$; at long wavelengths, the coupling constant $`g4\pi \mathrm{}^2a/M`$, where $`a`$ and $`M`$ are respectively the $`s`$-wave scattering length and the mass of the atom. Only the case $`a>0`$ will be considered here. The chemical potential, $`\mu `$, fixes the number of atoms $`N_0`$ in the condensate (the contribution from the excitations is ignored in the Bogoliubov approximation), and $`ϵ_n`$ are the collective excitations of the system. Note that the condensate is the zero-energy solution of the Bogoliubov equations, $`\psi (𝐫)=u_0(𝐫)=v_0(𝐫)`$; the excitation frequencies are measured with respect to the ground state energy $`\mu `$.

The confinement due to the magnetic trap is well-described by a completely anisotropic harmonic oscillator potential:

$$V_{\mathrm{trap}}=\frac{M}{2}\left(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2\right).$$
(5)   
The analysis of the GP and Bogoliubov equations is considerably simplified by re-expressing the trap frequencies $`\omega `$ and coordinates in terms of the scaled variables $`(\omega _x,\omega _y,\omega _z)=\omega _0(1,\alpha ,\beta )`$ and $`(x,y,z)d_0(x,y,z)`$, where $`d_0=\sqrt{\mathrm{}/M\omega _0}`$ is the characteristic oscillator length in the $`\widehat{x}`$-direction. As a result, the coupling constant becomes $`g4\pi \eta _0`$, where we define $`\eta _0N_0a/d_0`$. The Schrödinger operator is then

$$\widehat{L}=\frac{1}{2}\stackrel{}{}^2+\frac{1}{2}\left(x^2+\alpha ^2y^2+\beta ^2z^2\right)+4\pi \eta _0\psi ^2(𝐫).$$
(6)   
All energies (including $`\mu ,ϵ_n`$) are now given in trap units $`\mathrm{}\omega _0`$. With these choices, the condensate and excited-state wavefunctions are normalized to unity in the rescaled (dimensionless) coordinates:

$`{\displaystyle 𝑑𝐫\psi ^2(𝐫)}`$ $`=`$ $`1;`$ (7)   
$`{\displaystyle 𝑑𝐫\left[u_n^2(𝐫)v_n^2(𝐫)\right]}`$ $`=`$ $`1.`$ (8)   
In the next section we discuss the numerical methods used to solve these equations.

## II NUMERICAL METHODS

The discrete variable representation is particularly useful for the problem of trapped Bose-condensed atoms. As discussed in the Introduction, all potential energy operators in the Hamiltonian are local and therefore have a diagonal matrix representation, while the kinetic energy still has a fairly simple and sparse structure. Since only the low-lying modes are relevant to the static properties of the weakly interacting Bose gas at low temperatures, the sparse Hamiltonian matrices are ideally suited to iterative techniques for the determination of eigenvalues. Since these methods are dominated by matrix-vector multiplies, a structured and sparse matrix offers considerable computational savings.

An important feature of the GP equation is that the interaction potential depends non-linearly on the solution. Thus, for a given $`N_0`$, the lowest eigenvalue $`\mu `$ and eigenvector $`\psi (𝐫)`$ must be found self-consistently. In practical terms, a robust iterative method is required. Indeed, as discussed below, no solution to the GP equation may be found simply by straight iteration when the universal scaling parameter $`\eta \alpha \beta \eta _0\genfrac{}{}{0pt}{}{>}{}10`$. Since $`a/d_0`$ is typically of order $`10^310^2`$, only a few thousand atoms in the trap (several orders of magnitude fewer than are experimentally relevant) would be sufficient to prevent the solution to the GP equation by iterative methods. We have developed a variant of a well-known technique, called direct inversion in the iterative subspace (DIIS) to render convergent the self-consistency process for the large $`N_0`$ of practical interest. DIIS complements alternative numerical approaches, such as the method of steepest descents (imaginary-time propagation), which have been successfully applied to the trapped Bose condensates .

### A DVR Techniques

The basic ideas of a discrete variable representation are quite old. The earliest applications  were designed to simplify the calculation of certain classes of matrix elements that appeared in finite basis set, variational calculations. Subsequent authors began to use the DVR more directly and, in some instances, to view it as the primary representation for the problem under consideration. The viewpoint of the current authors is that the DVR follows naturally from a particular choice of a finite basis set, one which is mathematically linked to Gauss quadratures.

Let us consider a basis of functions, $`\left\{\varphi _n(x),n=1,N\right\}`$, perhaps satisfying some set of boundary conditions over a finite or infinite interval, which is complete enough for expanding any unknown function over that interval to sufficient accuracy. For simplicity we assume that these functions form an orthonormal basis for the space. What we seek is a complementary set of “coordinate eigenfunctions”, $`\left\{u_i(x),i=1,N\right\}`$, and a generalized quadrature consisting of roots and weights $`\left\{x_k;w_k,k=1,N\right\}`$, such that

$$u_i(x_k)=\delta _{i,k}.$$
(9)   
We now expand the unknown $`u_i(x)`$ in the set of functions $`\varphi _n(x)`$,

$$u_i(x)=\underset{n=1}{\overset{N}{}}\varphi _n(x)\varphi _nu_i,$$
(10)   
and assume it is possible to evaluate the overlap integral, $`\varphi _nu_i`$, sufficiently accurately using the generalized quadrature rule so that

$`\varphi _nu_i=w_i\varphi _n(x_i);`$ (11)   
$`u_i(x)=w_i{\displaystyle \underset{n=1}{\overset{N}{}}}\varphi _n(x)\varphi _n(x_i).`$ (12)   
The result (11) follows directly from Eq. (9) and the integration rule for a Gauss quadrature:

$$fg_a^b𝑑xw(x)f(x)g(x)\underset{k}{}w_kf(x_k)g(x_k),$$
(13)   
where $`w(x)`$ is a non-negative weight function. Moreover, the delta-function property (9) of the coordinate eigenfunctions is exactly satisfied by Eq. (12) at the quadrature points $`x_k`$, since the $`\varphi _n`$ form a complete set:

$$\underset{n=1}{\overset{N}{}}\varphi _n(x_k)\varphi _n(x_i)=\frac{\delta _{i,k}}{w_i}.$$
(14)   
It should be underlined that the $`u_i(x)`$ defined by Eq. (12) are highly localized in the vicinity of the quadrature points but are not true delta-functions, since they are finite for $`xx_k`$.

Since the coordinate eigenfunctions $`u_i(x)`$ are defined as continuous functions of $`x`$, it is possible to differentiate and tabulate them using the known properties of the $`\varphi _n(x)`$. The crucial step, of course, is to be able to find the set of basis functions and the related generalized quadrature which justify Eqs. (11,12).

While no conditions have thus far been placed on the set $`\varphi _n(x)`$, other than orthonormality and completeness, in practice the requirement that the $`u_i(x)`$ satisfy Eq. (9) and Eqs. (11,12) limits the basis functions to the classical orthogonal polynomials. These are defined by a three term recursion relation of the form

$$\beta _j\varphi _j(x)=(x\alpha _j)\varphi _{j1}(x)\beta _{j1}\varphi _{j2}(x),$$
(15)   
with the properties

$$\varphi _i\varphi _j=\delta _{i,j},\varphi _0(x)=\mathrm{const}.$$
(16)   
where the inner product is defined in Eq. (13). Since Eq. (15) may be interpreted as the definition of the coordinate operator $`x`$ in the polynomial basis, the orthonormal eigenvectors of this tridiagonal matrix must diagonalize the coordinate operator. It is therefore not a coincidence that the associated eigenvalues are the generalized Gauss quadrature points associated with the original basis $`\varphi _n(x)`$. For this set of $`N`$ functions there is an associated generalized Gauss quadrature consisting of $`N`$ points and weights which assure us that Eqs. (11,12) are satisfied exactly.

There are as many ways to define the $`u_i(x)`$ as there are classical orthogonal polynomials. The most natural choice, however, is to use the Lagrange interpolating functions, which explicitly diagonalize Eq. (15). The Lagrange polynomials may be defined at the associated Gauss quadrature points by:

$$v_i(x)=\underset{k=1}{\overset{N}{}}^{}\frac{xx_k}{x_ix_k},$$
(17)   
where the prime denotes exclusion of the point $`x_i`$ in the product. The Gauss-Legendre quadrature points $`x_k`$ are defined on the interval \[-1:1\]; with the associated weights, one obtains the eigenfunctions:

$$u_i(x)=v_i(x)/\sqrt{w_i}.$$
(18)   
These polynomials have the “delta-function property” by construction, and are normalized so that they form an orthonormal set under the generalized $`N`$-point Gaussian quadrature:

$$u_iu_j=\delta _{i,j}.$$
(19)   
This definition ensures that there is an equality between exact integration and integration by quadrature for any integrand which may be represented as a polynomial of degree $`(2N1)`$ or smaller. In the present work, both Lagrange functions as well as a DVR based on a Hermite polynomial basis are considered.

Since the $`u_i(x)`$ diagonalize the coordinate operator, the matrix element of any operator $`𝒪(x)`$ which is a local function of $`x`$ satisfies

$$u_i𝒪(x)u_j=\delta _{i,j}𝒪(x_i).$$
(20)   
The DVR basis not only considerably simplifies the evaluation of many matrix elements of the Hamiltonian, but leads to a sparse representation as well. For many large matrix problems the only practical methods of diagonalization or matrix inversion require the operation of the Hamiltonian matrix on some known vector. Sparsity is a key ingredient to performing this operation efficiently.

While the result (20) is an identity within a particular $`N`$th-order Gauss quadrature, it is only exact when the product of $`u_i(x)`$, $`u_j(x)`$ and the local operator $`𝒪(x)`$ is a polynomial of degree $`2N1`$ or smaller . The present formalism, therefore, is particularly conducive to the solution of the GP equation in the limit of large $`N_0`$. In this case, the contribution of the kinetic term to the total energy is small. The square of the condensate wavefunction may then be written in the so-called Thomas-Fermi (TF) approximation:

$$\psi ^2(𝐫)\left[\mu _{\mathrm{TF}}\frac{1}{/}2(x^2+\alpha ^2y^2+\beta ^2z^2)\right]/4\pi \eta _0,$$
(21)   
where the normalization condition (7) yields the chemical potential in the TF limit:

$$\mu _{\mathrm{TF}}=\frac{1}{/}2\left(15\alpha \beta \eta _0\right)^{2/5}\frac{1}{/}2\left(15\eta \right)^{2/5},$$
(22)   
in units of $`\mathrm{}\omega _0`$. In the TF limit, both the interaction and confining potentials appearing in the Schrödinger operator $`\widehat{L}`$ (6) are evidently second-order polynomials. For a finite but large number of atoms, however, there is a small exponential tail near the boundary of the condensate cloud resulting from the finite kinetic energy. As shown in Sec. III, this exponential behavior may be effectively captured even when using a comparatively low-order basis.

For the present three-dimensional calculation, the condensate and excited-state wavefunctions are expanded using Cartesian coordinates in a product basis of coordinate functions for each dimension; in the condensate case, for example, one writes:

$$\psi (𝐫)=\underset{i,j,k}{}c_{ijk}u_i(x)u_j(y)u_k(z).$$
(23)   
The components of the three-dimensional kinetic energy $`T_{x,y,z}\frac{1}{/}2_{x,y,z}^2`$ has the following representation in the DVR product basis:

$`u_iu_ju_kTu_lu_mu_n`$ (24)   
$`=T_{i,l}\delta _{j,m}\delta _{k,n}+T_{j,m}\delta _{i,l}\delta _{k,n}+T_{k,n}\delta _{i,l}\delta _{j,m}.`$ (25)   
Thus, the Hamiltonian matrix separates into a sum of dense, one-dimensional kinetic energy matrices times Kronecker delta functions in the remaining variables, plus a purely diagonal matrix associated with the potential terms. While multidimensional expansions of this kind (23) exist in other coordinate systems, the choice of a product basis in Cartesian coordinates yields a separable kinetic energy operator which ensures a sparse matrix representation of the Hamiltonian in the multidimensional DVR space. The expansion coefficients $`c_{ijk}`$ found by diagonalization are proportional to the values of the wavefunction at the appropriate quadrature points.

### B Iteration and the DIIS Method

The central numerical difficulty in the solution of the GP equation is the non-linearity associated with the two-body interactions. Since the magnitude of the Hartree potential is proportional to the number of atoms in the condensate, the self-consistent solution when $`N_0`$ is large can become problematic. Many calculations reported to date have “inverted” the search; that is, the value of $`\mu `$ is fixed while the wavefunction and associated value of $`N_0`$ are obtained which satisfy the appropriate boundary and normalization conditions. Such a procedure is straightforward in one dimension (though somewhat less so in two), since a direct numerical integration of the GP equation can be implemented. In three dimensions, a root search procedure may be devised to accomplish the same thing, but the calculation would be extremely time-consuming. With the method of steepest descents, the solution of the three-dimensional GP equation could be obtained by imaginary-time propagation. Though this approach has been shown to yield accurate results for large numbers of atoms , it is not obvious that time-dependent techniques should be most suitable for the solution of an eigenproblem.

An elegant and (as shown below) inexpensive method for the solution of the GP equation is to fix $`N_0`$ and to self-consistently solve for $`\mu `$ and $`\psi (𝐫)`$. Unfortunately, a direct iteration of Eq. (1) does not usually converge to the global energy minimum, even for relatively small condensates ($`N_010^3`$, $`\eta 1`$). Initially, we attempted several simple schemes in order to improve convergence. One approach was to use the TF approximation (21) to begin the iteration sequence. This would seem to be an excellent starting point when the condensate is large and the TF solution is a reasonable approximation. In practice, we have found that the poor behavior of the TF wavefunction near the condensate boundary did not, in general, make this a viable procedure. Even correcting the TF result using a boundary-layer perturbation approach did not seem to greatly improve convergence.

Evidently, a robust numerical process is required which rapidly damps out the errors as the iteration sequence proceeds. One rather crude approach is to use a linear combination of the solutions from the $`(i1)^{\mathrm{th}}`$ and $`i^{\mathrm{th}}`$ steps to initiate the iteration process for the $`(i+1)^{\mathrm{th}}`$ step ; the proper linear combination may be found by numerical experimentation. The approach used by the present authors, known as Direct Inversion in the Iterative Subspace (DIIS), is a more systematic version of the above that uses the information from all of the previous iterations. The DIIS method is well known in the quantum chemistry community, where it is used to accelerate the convergence of self-consistent field calculations. As discussed below, for small to intermediate numbers of atoms ($`10^3N_010^6`$), DIIS is found to be at least as computationally efficient as the method of steepest descents.

DIIS begins by defining an error which characterizes the convergence of the iteration process. While the error $`\widehat{e}`$ may be defined in any number of ways, one that appears to be particularly stable numerically is the commutator:

$$\widehat{e}[\widehat{H},\widehat{\rho }],$$
(26)   
where $`\widehat{\rho }(𝐫,𝐫^{})`$ is the density matrix $`|\psi \psi |`$. Evidently, $`\widehat{e}`$ vanishes at convergence. Note that both the density matrix and the error are matrices having $`S^2`$ elements, where $`S`$ is the dimension of the Hilbert space. The procedure is to expand the current expression for the Hamiltonian and error as a linear combination of the $`m`$ previous values

$`\widehat{H}^{m+1}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{m}{}}}a_i\widehat{H}^i;`$ (27)   
$`\widehat{e}^{m+1}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{m}{}}}a_i\widehat{e}^i,`$ (28)   
subject to the constraint that

$$\underset{i=1}{\overset{m}{}}a_i=1.$$
(29)   
Minimizing the squared norm of $`\widehat{e}^{m+1}`$ with the constraint (29), one obtains $`m+1`$ linear equations:

$$\underset{j=1}{\overset{m}{}}B_{i,j}a_j\lambda =0,$$
(30)   
where

$$B_{i,j}\widehat{e}_i\widehat{e}_j,$$
(31)   
and $`j=1,2,\mathrm{},m`$. The Lagrange multiplier $`\lambda `$ enforcing the constraint yields the squared norm of the error.

The practical implementation of the DIIS algorithm, especially for large Hamiltonian matrices, deserves some additional comment. Although it appears to be necessary to store the Hamiltonian matrix for each iteration, this is not the case. The only part of the Hamiltonian matrix that varies from iteration to iteration is the nonlinear potential, and this matrix is diagonal in the DVR. Inserting Eq. (31) into (26), the matrix elements $`B_{i,j}`$ are simply given by:

$$B_{i,j}=2\phi _i\phi _j\psi _i\psi _j2\phi _i\psi _j\psi _i\phi _j,$$
(32)   
where $`\phi _i(𝐫)\widehat{H}^i\psi _i(𝐫)`$; it is important to note that $`\widehat{H}^i`$ is the Hamiltonian matrix constructed from $`\psi _{i1}(𝐫)`$, whose minimal eigenvector is $`\psi _i(𝐫)`$. On each iteration, therefore, it is only necessary to store the solution vector as well as the vector which represents the operation of the Hamiltonian matrix on the solution vector. Since these objects are both of dimension $`S`$ rather than $`S^2`$, they may easily be stored in central memory even on current workstations as long as the number of DIIS iterations does not become too large.

To summarize, only the vectors $`\psi (𝐫)`$ and $`\phi (𝐫)`$ need to be stored at each iteration step. From these, the scalar products (32) required to set up the set of linear equations (29) and (30) may be generated. Once the $`m`$ unknowns $`a_i`$ as well as $`\lambda `$ have been found, a new nonlinear potential may be constructed from the past guesses,

$$V_\mathrm{H}^{m+1}=\underset{i=1}{\overset{m}{}}a_iV_\mathrm{H}^i,$$
(33)   
and the next self-consistent field cycle initiated. Furthermore, by storing and not destroying all previous elements of $`B_{i,j}`$, only an additional row and column need be computed at each step of the iteration. Since the dimension of $`B`$ is assumed to be small, the additional storage should not be a problem on most workstations.

With the TF expression for the condensate wavefunction as the initial guess, the DIIS algorithm yields a convergent solution for significantly larger numbers of atoms ($`N_010^5`$) than was possible with a direct iterative procedure. As $`N_0`$ continues to grow, however, the number of DIIS iterations to achieve convergence increases until eventually the process fails. It appears that small errors in the intermediate solution of the GP equation, particularly in the ‘tail’ region close to the condensate surface, are amplified by the nonlinearity of the potential during the iteration process. Combining DIIS with a more realistic initial guess, we have been able to achieve convergence for the large number of atoms ($`N_010^6`$) relevant experimentally. Two schemes used to improve the starting approximation for the Hartree potential have proved particularly valuable. By slowly increasing the number of atoms, the converged solution for a slightly smaller condensate is an excellent choice. Alternatively, one may solve a modified GP equation with an exaggerated kinetic energy contribution, then slowly ramp down the degree of exaggeration (this technique is similar in spirit to simulated annealing). A simple modification of the original DIIS procedure was also found to be quite helpful. The DIIS procedure is broken into cycles with some maximum number of iterations allowed per cycle. At the end of each cycle, the DIIS procedure is restarted using the best available solution. By numerical experimentation it was found that once the root mean square error (or alternatively $`\sqrt{\lambda }`$) is reduced to $`10^210^3`$, it becomes possible to restart DIIS and to rapidly reduce errors to $`10^610^7`$ with very few additional iterations.

### C Interpolation of the Wavefunction

As discussed above, in the limit of very large $`N_0`$ where the TF theory is believed to be valid, both the interaction and confining potentials are well-approximated by low-order polynomials in $`x,y,z`$. Since an $`N`$-point Gauss quadrature is able to integrate a $`(2N1)`$th-order function exactly, it should be possible to capture the essential behavior of the condensate wavefunction (squared) using only a very small number of quadrature points. One is therefore left with a somewhat surprising conclusion: when the number of atoms becomes very large, the resulting self-consistency problem should simplify considerably. A very coarse DVR grid has two obvious advantages. The reduced dimension of the Hamiltonian matrix would accelerate the solution of the eigenproblem at each DIIS iteration. In addition, the fewer degrees of freedom for spatial variations should enable the self-consistent solution for the condensate to converge with fewer DIIS iterations. The crucial unknown is whether the exponential tail at the condensate boundary, due to the finite kinetic energy contribution, varies sufficiently rapidly to invalidate the low-order approximation applicable in the TF limit.

Suppose that the gross features of the condensate wavefunction could be found using a low-order Gauss quadrature, corresponding to only a few DVR points in each spatial direction. The self-consistent solution obtained using such a ‘coarse grid’ would then make an excellent initial guess for a more accurate ‘fine grid’ calculation, if the interpolation between grids could be implemented successfully. By increasing the number of DVR points in the mesh once or perhaps a few times, it should be possible to rapidly converge the solution of the GP equation for virtually arbitrary numbers of atoms. The interpolated wavefunction on the quadrature points $`(x_l,y_m,z_n)`$ associated with the fine grid is approximately

$$\psi (x_l,y_m,z_n)\underset{ijk}{}c_{ijk}u_i(x_l)u_j(y_m)u_k(z_n),$$
(34)   
where the sum is taken over coarse-grid expansion coefficients and coordinate eigenfunctions. Since the DVR points corresponding to the two different Gauss quadratures are not generally coincident, the values of the coarse-grid coordinate basis functions at the fine grid points must be obtained explicitly using either Eq. (12) or (17).

### D Extraction of Eigenvalues and Eigenvectors

The direct extraction of all the eigenvalues and eigenvectors of Eqs. (1-3) becomes computationally impractical for very large matrices. In three-dimensional systems, a sizeable number of basis functions (and therefore quadrature points) are usually required in order to adequately represent the wavefunctions over all space. This is a particularly important consideration for the excitations, whose rapid spatial variations can only be captured by high-order polynomials. As a result, the dimension of the matrices can become very large, even after block-diagonalizing according to parity. A full diagonalization of the eigenproblem, using standard and widely-available routines based on variants of the Givens-Householder method, places considerable demands on both storage and cpu-time. A full diagonalization is not necessary since only the lowest eigenvalue ($`\mu `$) of the GP equation is required, and it is sufficient to determine merely the lowest-lying excitations of the Bogoliubov equations at zero temperature.

When the matrices are too large to fit in memory, iterative methods must be used in order to extract the relevant low-lying states. These techniques require the frequent operation of the Hamiltonian on a vector, and are practical only if there are relatively few non-zero matrix elements. Indeed, a sparse representation of the Hamiltonian is a crucial feature of the DVR and the separability of the kinetic energy operator. In addition, standard iterative procedures return reasonable approximations only to the extremal eigenvalues, in the spirit of a Rayleigh-Ritz variational method. The GP equation (1) yields a positive-definite, real, symmetric, and sparse matrix whose lowest eigenvalue may be found using variants of the well-known Lanczos  or Davidson  algorithms.

The Davidson method, which is widely used by quantum chemists, develops a small orthonormal subspace of vectors which is adequate to describe the eigenpairs of physical interest. The Hamiltonian is projected into this $`p`$-dimensional subspace and the Rayleigh-Ritz variational principle ensures that an upper bound to the lowest $`n`$ eigenvalues is obtained from the process. Since the computational procedure in the application of the Davidson algorithm may be unfamiliar to some readers, it is summarized below:

1. Choose an initial set of m orthonormal trial vectors $`𝐛_i`$ , where $`mn`$.
2. Calculate $`𝐡_i`$ = $`\widehat{H}𝐛_i`$, the effect of the Hamiltonian matrix on these vectors.
3. Calculate the Hamiltonian matrix $`𝐛_i|𝐡_j`$ from the information generated in steps 1 and 2.
4. Solve the resultant small eigenvalue problem for the current values of the wavefunctions, $`𝚿_q^{\mathrm{cur}}`$ and eigenvalues, $`E_q^{\mathrm{cur}}`$ by a direct method.
5. Calculate the residuals,

$$𝐫_q=\left(E_q^{\mathrm{cur}}\widehat{H}\right)𝚿_q^{\mathrm{cur}}$$
(35)   
and test for convergence. If all roots are converged stop. If not select the unconverged residuals for further improvement.
6. Using the set of residuals solve the equation,

$$\left(\stackrel{~}{𝐇}E_q^{\mathrm{cur}}\right)𝚽_𝐪=𝐫_𝐪$$
(36)   
where $`\stackrel{~}{𝐇}`$ is some approximation to the exact $`\widehat{H}`$.
7. Schmidt orthonormalize the $`𝚽_q`$ to the $`𝐛_i`$ and then append them to the $`𝐛_i`$ to enlarge the vector space.
8. Calculate $`𝐡_i`$ = $`\widehat{H}𝐛_i`$, for the appended vectors and the additional matrix elements needed to border $`𝐛_i|𝐡_j`$.
9. Return to step 4.

In the original Davidson algorithm, which was designed for diagonally dominant matrices, the approximate Hamiltonian used in step 6 was the diagonal. Then the calculation of the solution to Eq. (36) is trivial. For diagonally dominant matrices, this Jacobi preconditioning method produces a quickly converging set of new vectors, and the entire calculation is dominated by step 2 in the sequence above. Although the DVR has the desirable feature that it produces a sparse representation, the matrix is unfortunately far from being diagonally dominant. Without an alternative preconditioning technique, the calculation can require several ($`S`$) full cycles which becomes computationally prohibitive.

Our approach to the preconditioning problem is based on using a separable approximation to $`\stackrel{~}{𝐇}`$, i.e. a Hamiltonian which can be written as a sum of commuting operators in the three-dimensional space. For simplicity, we have chosen to use the bare trap Hamiltonian, but other more complicated choices can be made, subject to the separability requirement. Eq. (36) then becomes:

$$\left(\stackrel{~}{𝐇}_0^1+\stackrel{~}{𝐇}_0^2+\stackrel{~}{𝐇}_0^3E_q^{\mathrm{cur}}\right)𝚽_q=𝐫_q.$$
(37)   
In order to solve this equation, we transform to the basis set which diagonalizes the separable Hamiltonian, recognizing that the transformation can be written as a product of three one-dimensional unitary transformations. Thus, the residual vector on the right hand side of Eq. (35) is transformed from the DVR to the diagonal representation, an operation which scales as the size of the three-dimensional basis set times the one-dimensional basis set. This linear scaling with basis set size is a consequence of the separable form of the approximate Hamiltonian, yet preserves the most desirable features of the DVR. Once the residual has been transformed to the diagonal representation, Eq. (36) is trivially solved and then the solution vector is transformed back to the DVR representation. This procedure is analogous to using fast Fourier transforms (FFT) to solve multidimensional PDE’s. The FFT is more efficient (but less general) due to special character of the Fourier transformation but the separability is a key factor to its properties. By using this separable preconditioning, we have been able to reduce the number of Davidson iterations to manageable size. In addition, the vectors from a previous calculation for a smaller condensate are found to be excellent trial vectors to initiate the calculation. The result is a robust procedure which is quite efficient for large matrices.

The matrix associated with the Bogoliubov equations (2,3), in contrast, is non-symmetric and has real eigenvalues occurring in positive-negative pairs. This matrix must be addressed using routines such as those based on an Arnoldi or a non-symmetric Davidson approach. Since the low-lying excitations of the Bogoliubov equations are those closest in magnitude to zero (relative to the chemical potential) and are therefore the least extremal, it is most convenient to decouple the original equations using the following linear combinations:

$`f_n(𝐫)`$ $``$ $`u_n(𝐫)+v_n(𝐫);`$ (38)   
$`g_n(𝐫)`$ $``$ $`u_n(𝐫)v_n(𝐫).`$ (39)   
The resulting eigenproblem then becomes

$$\left(\widehat{L}\mu +2V_\mathrm{H}\right)\left(\widehat{L}\mu \right)f_n(𝐫)=ϵ_n^2f_n(𝐫),$$
(40)   
where $`ϵ_0=0`$ corresponds to the ground-state energy relative to the chemical potential. Once the functions $`f_n`$ are found, the left eigenvectors $`g_n`$ may be obtained directly by utilizing

$$g_n(𝐫)=ϵ_n^1\left(\widehat{L}\mu \right)f_n(𝐫)$$
(41)   
subject to the normalization condition $`f|g=1`$, obtained from Eq. (8). While the resulting eigenproblem remains non-symmetric (since the condensate density does not commute with the operator $`\widehat{L}`$), the dimension of the matrix has been reduced by a factor of two. Furthermore, the eigenvalues $`ϵ_n^2`$ are positive-semidefinite and real . In practice, the composite operator on the left side of Eq. (40) is never determined explicitly; doing so at each DIIS iteration not only would be time-consuming, but also would lead to a dense matrix. Rather, at each Arnoldi/Davidson iteration the vector is multiplied in turn by the two distinct sparse operators.

## III results and discussion

The techniques discussed in the previous section are particularly useful under two circumstances: when the potential is completely anisotropic, and when the number of condensate atoms is very large. For illustrative purposes, results are presented for Bose-condensed sodium atoms confined in traps with three different geometries, all of which have been realized experimentally. The completely anisotropic trap considered here is a TOP trap with angular frequencies in the natural ratio  $`(\omega _x,\omega _y,\omega _z)=\omega _0^\mathrm{A}(1,\sqrt{2},2)`$ where $`\omega _0^\mathrm{A}=354\pi `$ rad/s; BEC in such a system has been recently observed by Phillips et al.. Other geometries considered are the cigar-shaped Ioffe-Pritchard trap of W. Ketterle with cylindrical symmetry and $`(\omega _x,\omega _y,\omega _z)\omega _0^\mathrm{C}(1,13.585,13.585)`$, where $`\omega _0^\mathrm{C}33.86\pi `$ rad/s, and the approximately spherical ‘4D’ potential of L. Hau with $`\omega _0^\mathrm{S}87`$ rad/s. All the calculations assume an $`s`$-wave scattering length for sodium of $`a=(52\pm 5)a_0`$, where $`a_0`$ is the Bohr radius.

### A Condensate

The solution of the GP equation (1) on the coarse and fine grids employed 60 and 180 harmonic oscillator basis functions (Hermite polynomials with a Gaussian prefactor) in each spatial direction, respectively. Coordinates were rescaled through $`\{x,y,z\}\{x,y/\alpha ,z/\beta \}`$; the non-linear coefficient becomes $`g=4\pi \eta =4\pi \alpha \beta \eta _0`$ in order to ensure the proper normalization of the condensate density (7). With this choice, the condensate cloud becomes almost spherical at large particle numbers, with an approximate radius given by the TF expression $`R_{\mathrm{TF}}=(15\eta )^{1/5}`$, in units of the characteristic trap length $`d_0=\sqrt{\mathrm{}/M\omega _0}`$. DVR points significantly beyond the TF radius were ignored. The ground state was assumed to have totally even parity, so the GP equation was solved in a single octant. The basis functions were eigenfunctions of a trap with frequencies reduced from their actual values by a factor of 20, in order to decrease the density of points. The resulting Hilbert spaces for the coarse and fine grids then have dimension $`\mathrm{3\hspace{0.17em}566}`$ and $`\mathrm{18\hspace{0.17em}685}`$ respectively, for all geometries and condensate numbers; both sizes are considerably reduced from the impractical values of $`\mathrm{216\hspace{0.17em}000}`$ and $`\mathrm{5\hspace{0.17em}832\hspace{0.17em}000}`$ one naively might have obtained.

In Table I, the chemical potentials obtained numerically using the coarse grid are given as a function of the number of atoms $`N_0=2^q`$ in the condensate for the spherical, cylindrical, and anisotropic geometries described above. The TF values, from Eq. (22), are included for comparison. The number of DIIS iterations required to yield a convergent solution to the GP equation increases with the number of atoms; in the cylindrical case, DIIS failed to converge for $`N_0=2^{19}=\mathrm{524\hspace{0.17em}288}`$ and greater. Choosing a finer grid always reduces the number of DIIS iterations; although increasing the number of points yields additional degrees of freedom for the Hartree potential, the variations of the condensate density (particularly in the surface region) are better captured. In the cylindrical case, which has a cigar shape, the coarse grid has too few DVR points in the two strongly-confining directions to adequately capture the behavior of the condensate at large $`N_0`$. For $`q>18`$, therefore, initial grids were chosen to contain a larger number of points not exceeding $`10^4`$. Since more points imply a larger eigenproblem at each iteration, practical calculations require as coarse a grid as possible.

The convergence of the chemical potentials with the number of condensate atoms is shown in Fig. 1. Both the cylindrical and spherical traps give rise to relatively slow convergence of the chemical potential to the TF value compared with the anisotropic case. This trend is expected for the spherical trap, where the confinement is extremely weak. In the cylindrical case the condensate is strongly confined, but only in the radial direction where the motion of the atoms is more or less frozen. Since the axial kinetic energy can be large, however, the cylindrical trap is effectively loose and one-dimensional. In contrast, the fully anisotropic trap considered here is relatively tight in all directions, and the chemical potential converges to the TF value more rapidly.

In the TF limit, most relevant quantities, such as the mean condensate radius $`R_{\mathrm{TF}}=(15\eta )^{1/5}d_0`$ and the chemical potential $`\mu _{\mathrm{TF}}=\frac{1}{/}2(15\eta )^{2/5}\mathrm{}\omega _0`$, are functions of the universal scaling parameter $`\eta =\alpha \beta (N_0a/d_0)`$. Similarly, the first-order correction to the TF chemical potential $`\mu _{\mathrm{TF}}`$, taking into account the average kinetic energy  and contributions from the potential energies , is proportional to $`R_{\mathrm{TF}}^2\mu _{\mathrm{TF}}^1`$. The fits to the numerical data of $`\mathrm{\Delta }\mu (1\mu _{\mathrm{TF}}/\mu )\gamma /\mu _{\mathrm{TF}}`$, shown in Fig. 1, are in reasonable agreement with this behavior (except for the cylindrical case for very small numbers of atoms $`N_010^4`$). One might naively expect, therefore, that the convergence of the chemical potential to the TF limit as a function of $`\eta `$ (which is proportional to the coefficient of the non-linear term in the GP equation) would be independent of trap geometry. As shown in Fig. 1, this is in fact not the case. The data for the cylindrical case converge far more slowly than those for the other geometries, though the magnitudes of the TF chemical potentials for a given value of $`\eta `$ are identical (note that the values of $`\mathrm{\Delta }\mu `$ shown in the figure for the cylindrical case are reduced from their actual values by an order of magnitude in order to facilitate comparison). The kinetic energy contribution to the chemical potential (and therefore to the total energy) is evidently strongly dependent on trap geometry. Thus, a large non-linear term in the GP equation does not necessarily mean that the TF approximation adequately represents the system.

In order to verify the accuracy of the solution obtained on the coarse grid, the condensate wavefunction was interpolated onto a $`\mathrm{18\hspace{0.17em}685}`$-point grid derived from a DVR basis with 180 polynomials in each direction (refer to Sec. II C for details on the interpolation technique). The solution of the GP equation was again converged on the fine grid. For small to intermediate values of the nonlinear coupling ($`\eta 3000`$), a single interpolation between the two meshes was sufficient to yield a solution on the fine grid with only a few ($`10`$) additional iterations. As $`\eta `$ increased, two or more successive interpolations and re-convergences were generally required before the solution on the finest grid could be obtained. In all cases, the values of the chemical potential at the two extremes were found to be identical to at least three decimal places.

The condensate wavefunctions obtained on the coarse and fine grids are compared in Fig. 2, for the case of $`N_0=2^{19}=\mathrm{524\hspace{0.17em}288}`$ atoms in the completely anisotropic trap. The coarse-grid condensate profile along a given axis appears rather crude, particularly in the surface region where the wavefunction varies rapidly; for a given basis, convergence is only obtained at the DVR points. Nevertheless, this wavefunction interpolated is virtually indistinguishable from the converged solution on the fine mesh. As $`\eta `$ increases, however, the length of the tail at the surface of the condensate shortens. Eventually, the coarse grid will have too few points in this crucial region to adequately capture the rapid variations. In this case, a large jump in mesh size results in an interpolated wavefunction that more poorly represents the self-consistent result.

### B Excitations

The excitations of a condensate in a fully anisotropic harmonic trap have been completely classified , and have been explicitly obtained in the low-density  and TF  limits. The states are polynomials of order $`N=l+m+n`$ and are labeled by the total parity $`P=(1)^{l+m+n}`$, where the quantum numbers $`(l,m,n)`$ represent the order of the polynomials along the $`(x,y,z)`$-directions in the non-interacting limit. In the strongly-interacting (or hydrodynamic) regime, there are four odd and four even parity low-lying modes with energies $`ϵ=\sqrt{l+\alpha ^2m+\beta ^2n}`$ in units of $`\mathrm{}\omega _0`$, where $`(l,m,n)`$ can be either $`0`$ or $`1`$. The ground state, relative to the chemical potential, has quantum numbers $`(0,0,0)`$. The only states with $`N=1`$ are the odd-parity dipole modes, where the center of mass oscillates with the three trap frequencies. For $`N=2`$, there are six quadrupole oscillations with even parity and stationary center of mass. Three of these have energies given by the expression above, with $`(l,m,n)=(1,1,0)`$ and cyclic permutations. The other three are the solutions of the secular equation:

$$\left|\begin{array}{ccc}3ϵ^2& 1& 1\\ 1& 3ϵ^2/\alpha ^2& 1\\ 1& 1& 3ϵ^2/\beta ^2\end{array}\right|=0.$$
(42)   
For $`(\alpha ,\beta )=(\sqrt{2},2)`$, the geometry considered here, one readily obtains $`ϵ=\sqrt{8\pm 4\sqrt{2}}`$ and $`\sqrt{5}`$. In the fully anisotropic case, all the excitation energies (except those for the odd-parity dipole modes) decrease with $`N_0`$ .

The low-lying excitation energies for a completely anisotropic trap have been computed numerically, and are shown in Fig. 3. All the calculations were obtained using the $`\mathrm{3\hspace{0.17em}566}`$-point grid. While small finite-size and coarse-graining effects are present, as evidenced by the small fluctuations in the ground-state energy $`ϵ=0`$, the results closely match the TF predictions described above. As expected, the odd-parity dipole modes are independent of the number of atoms in the trap. All the other frequencies depend strongly on $`N_0`$, decreasing from their non-interacting values in agreement with perturbative calculations .

While the frequencies of all the lowest-lying modes attain their large-number values by $`N_010^6`$ in the relatively strong anisotropic trap, the convergence to the hydrodynamic limit slows as the quantum numbers increase and the confinement is weakened. In Fig. 4, the low-lying excitations $`ϵ_{nl}`$ of the loose spherical trap are shown as a function of the number of atoms in the condensate, up to $`N_0=2^{22}4.2\times 10^6`$. In the hydrodynamic limit, the energies are given by $`ϵ_{nl}=\sqrt{l+2n(n+l+3/2)}`$ . The lowest number-dependent mode $`ϵ_{02}`$ agrees with its hydrodynamic value $`\sqrt{2}`$ to less than a percent by $`N_010^6`$. The numerical value of $`ϵ_{04}`$, in contrast, differs from its limiting value of $`2`$ by approximately $`2\%`$ even when the number of atoms is as large as $`4\times 10^6`$ atoms. Evidently, the magnitude of the excitation energy, relative to the ground state $`\mu `$, does not alone provide a sufficient indication of its convergence to the hydrodynamic limit. A similar number-dependence for the higher-lying excitations of a cylindrically-symmetric condensate has also been recently obtained .

## IV conclusions

A numerical procedure is introduced for the investigation of an interacting Bose gas at zero temperature confined in a completely anisotropic trap. The central feature of the technique is the use of the discrete variable representation (DVR) as the primary basis for the calculations. The DVR combines the best features of grid and basis-set techniques. All local operators are diagonal, so the evaluation of interaction matrix elements becomes trivial. While the kinetic energy has a more dense representation, it may be evaluated analytically to high accuracy by exploiting the underlying polynomial basis used to define the DVR. Furthermore, the kinetic energy needs to be evaluated only once.

In the present method the condensate density is determined self-consistently; for fully three-dimensional systems, this approach is considerably more efficient than conventional root-search algorithms. At each iteration, the ground-state wavefunction is obtained by iteratively diagonalizing the sparse GP Hamiltonian, using either a Lanczos or Davidson method. Convergence of the self-consistent solution to the GP equation is substantially hastened by employing direct inversion in the iterated subspace (DIIS). As the non-linear coefficient of the GP equation becomes very large, it often becomes necessary to employ more sophisticated techniques, including number or kinetic energy ramping and multigrid interpolation. In general, DIIS becomes more expensive as the number of atoms increases; it is conceivable that an alternative method such as imaginary-time propagation becomes more efficient than DIIS in the regime $`N_0\genfrac{}{}{0pt}{}{>}{}10^6`$.

The convergence of the chemical potential and the low-lying collective excitations is investigated as a function of trap geometry. The chemical potential is found to approach its Thomas-Fermi value more slowly as the confinement weakens or the degree of anisotropy becomes more pronounced; the convergence does not scale universally with the magnitude of the non-linear coefficient. The excitations of a completely anisotropic condensate have been calculated numerically, and for large numbers of atoms agree with the Thomas-Fermi predictions . For a very weak spherical trap, the collective frequencies converge to their hydrodynamic values  more slowly.

###### Acknowledgements.

The authors would like to acknowledge useful discussions with T. Bergeman, C. W. Clark, R. J. Dodd, M. Edwards, A. L. Fetter, E. Hagley, W. D. Phillips, and E. Tiesinga. This work was supported by the U.S. office of Naval Research.