Methods for a Nonuniform Bose Gas.

We review mathematical methods for the treatment of a system of Bose particles with nonuniform density. The use of the pseudopotential is explained, especially with respect to negative scattering lengths. It is emphasized that the delta-function potential produces no scattering in three dimensions, and should not be used in the Bogoliubov self-consistent field method, which is variational in nature. A common misuse of the Bogoliubov method at finite temperatures is pointed out. A Gaussian variational method is proposed.

1. Introduction

The observation of Bose-Einstein condensation in atomic traps [1]–[3] has renewed theoretical interest in studying a system of Bose particles with a nonuniform density. The experiments involve atoms of mass m condensed in a harmonic trap characterized by frequency ω. The relevant parameters, as discussed more fully in Refs. [4] and [5], have the following orders of magnitude: where a is the s-wave scattering length, λ the thermal wavelength, L the size of the system, and n is the average particle density. Thus, na3 << 1, a/L << 1, a/l << 1, and we can treat the system as a low-temperature weakly-interacting dilute gas using well-known methods [6].

In this paper, we review the following theoretical methods, and point out some common errors and misconceptions:

The pseudopotential method: One replaces the potential by a δ-function, and, with the help of a Bogoliubov transformation, obtains a weak-coupling expansion. This method seems adequate to describe current experiments.

The self-consistent field method: This is variational, and capable of treating strongly interacting cases in principle. It reduces to the pseudopotential method in appropriate limits; but one should not use a δ-function potential in a variational calculation, because such a potential in three spatial dimensions has no effect, if treated exactly.

The Gaussian variational method: One uses a Gaussian trial wave function, and obtains results similar to the self-consistent field method. It has the advantage of easily adaptable to the description of time-dependent phenomena.

All three methods are closely related to one another. Although we will not use Feynman graphs, it might be relevant to note that the Bogoliubov transformation corresponds to summing “one-loop” graphs. The self-consistent field method consists of readjusting parameters in the one-loop sum in a variational sense. The Gaussian method includes all one-loop contribution, plus parts of higher-loop contributions. Of these methods, the straightforward one-loop summation yields an expansion in na3. On the other hand, the accuracy of the other methods cannot be ascertained.

We do not attempt a comprehensive review, and the references cited are by no means complete. In a subject with such a long history, it is difficult to be aware of all relevant sources, and we apologize for any omissions [7].

2. Pseudopotential Method

In low-temperature calculations, one can replace the interatomic potential with a δ-function potential (4πaℏ2/m)δ(), where a > 0 is the s-wave scattering length. Such a potential should be treated in low-order perturbation theory only, because the δ-function potential in three spatial dimensions has no effect if treated exactly [8].

Consider the Schrödinger equation where E is 2m/ℏ2 times the energy. Taking the Fourier transform of both sides leads to where ϕ() is the Fourier transform of ψ (). This must hold for all , including k2 = E. Therefore ψ(0) = 0. This condition is automatically satisfied for all partial waves except s-waves. For s-waves, this seems to require that ψ() be discontinuous at = 0; but because of the singular nature of the potential at = 0, one should be careful. By treating the δ-function as the limit of a square well, one can verify that in three dimensions the s-waves are the same as those of a free particle except at = 0, where they drop to zero discontinuously. This behavior, however, does not lead to any scattering, nor level shift.

As a more careful analysis [9] shows, the δ-function potential must be regarded as the first term in a pseudo-potential, which correctly reproduces all scattering phase shifts, and depends on these phase shifts as input parameters. For the simple case of a hard-sphere potential, for which there is only one parameter a, the hard-sphere diameter, the pseudopotential has the form

The differential operator (∂/∂r)r removes spurious singularities in the wave function, which would otherwise make the ground-state energy diverge. In low-order perturbation theory, one can replace this operator by a simple subtraction procedure for the energy.

Using the first term in the pseudopotential, and making a Bogoliubov transformation, one can calculate exactly the first few terms of an expansion for the ground-state energy per particle of a uniform hard-sphere Bose gas [10, 11]:

This shows that the energy is not analytic at a = 0, and therefore cannot be continued to negative values of a.

Fetter [12] has given a systematic treatment of the nonuniform case based on the Bogoliubov transformation. For application to the atomic-trap experiments, one can also adapt the results of the uniform case using a Thomas-Fermi approximation [5].

What happens when the scattering length is negative? In this case, the two-body scattering is dominated by the attractive part of the potential. The attraction may or may not be strong enough to produce a two-particle bound state; but it will lead to an N-body bound state, for a sufficiently large N. The many-body problem is very different from that of the repulsive case, for it depends on the details of the potential. Consider two potentials with the same negative scattering length, one everywhere attractive, and the other having a repulsive core. The many-body system with the purely attractive potential will collapse spatially, whereas the one with hard core will yield extensive energies.

When the scattering length is negative, we must therefore use an actual potential, for example, one with a hard core plus an attractive tail. We can replace the hard-core part by a pseudopotential, and do the Bogoliubov transformation. In this manner, it has been shown [6] [13] that a dilute Bose system with such a potential exhibits a first-order phase transition separating a gas phase from a liquid phase. The transition line terminates in a critical point. The Bose-Einstein condensation occurs in the liquid phase, and divides it into liquid I and liquid II. In short, one reproduces qualitatively the phase diagram of 4He. Since the calculation is perturbative, one can follow the collapse of the phase diagram to that of the ideal Bose gas when the potential is turned off.

The results described above pertains to an infinite Bose system in the thermodynamic limit, and may be modified for a finite system. The authors of Ref. [2] observed Bose condensation in an atomic trap filled with in 7Li, which has a negative scattering length. In this case, the smallness of the system may have rendered the first-order transition indistinct. It is also possible that the observed phenomenon is metastable.

3. Variational Methods and the Energy Gap

A variational approach to the ground-state energy of a Bose gas was introduced by Girardeau and Arnowitt [14], who use (wisely) an actual potential instead of the δ-function. Their results are for a uniform gas, but otherwise very similar to what we obtain later in the self-consistent field method. The excitation spectrum in this approach, however, contains an energy gap, contrary to general expectations [15], and to the perturbation-theory results in the hard-sphere case [10].

A variational method is designed to yield the best ground-state energy, and may not yield the excitation spectrum accurately. In this instance, we know that the gap should have been filled with phonons, which can be described using general principles. We may therefore view the variational principle as a way to obtain an approximate ground-state wave function, and build the phonon states upon this ground state.

Feynman [16] argues that the low-lying excited states of a Bose system are purely phonon states, owing to the statistics, and suggests the following one-phonon wave function of momentum : where Ψ0 is the ground-state wave function. The excitation energy is where S is the liquid structure factor where with a the annihilation operator of a free particle of momentum ℏ, and Ω the volume. From general principles [17,18], we expect where c is the sound velocity as computed from the compressibility of the ground state. This leads to the phonon spectrum

Feynman’s argument can be sharpened by recasting it as the statement that the longitudinal sum-rules are saturated by the phonon states [17] [18]. All the equations above are verified in the dilute hard-sphere Bose gas [10] [18].

Other types of states that can be built upon Ψ0 are those describing superfluid flow [10] [18]: where α () is the superfluid phase, which is related to the superfluid velocity through

With this wave function one can describe quantized vortices.

4. Self-Consistent Field Method: Ground State

The self-consistent field method was first used by Bogoliubov in his treatment of superconductivity, and is therefore also known as the “Bogoliubov method.” We shall follow a concise formulation due to De Gennes [19]. A recent review of similar methods is given by Griffin [20].

We use a quantized-field representation with field operator Ψ (), and canonical conjugate iΨ†():

The Hamiltonian is where where Vext() is the external potential, and μ is the chemical potential. We displace the field operator by its ground-state expectation: where ψ () is the displaced field operator, and

The function ϕ () describes a nonuniform condensate, and ψ () annihilates particles not in the condensate.

Variational trial states are taken to be the eigenstates of an effective Hamiltonian that is quadratic in the field operators, chosen to be a mean-field version of the actual Hamiltonian. There are two variational functions ρ (,) and Δ (,), which will turn out to be ρ (,) = 〈ψ†() ψ ()〉, and Δ (,) = 〈ψ ()ψ ()〉. That is, ρ and Δ are, respectively, the direct and off-diagonal density correlation functions of the uncondensed particles.

The effective Hamiltonian is chosen to be where H(0) is independent of ψ, and H(1) and H(2) terms are respectively linear and quadratic in ψ: where

For arbitrary ρ and Δ, we can diagonalize the effective Hamiltonian through a generalized Bogoliubov transformation: where η and are annihilation and creation operators satisfying

To preserve the canonical commutators for ψ (), we require

We now require Heff to have the form or, equivalently,

These conditions determine the functions v, u and ϕ. The condition for ϕ comes from the requirement that in Heff the terms linear in ψ() vanish. The resulting equations are

To determine ρ and Δ, let |Φ〉 be the lowest eigenstate of Heff:

We require that the ground-state energy be at a minimum with respect to variations in ρ and Δ: where the brackets denote expectation value with respect to |Φ〉. In calculating the expectation value above, we use Wick’s theorem:

The calculation is facilitated by noting that δ 〈Heff〉 = 0, since Φ is an eigenstate. Thus we can use the equivalent condition

The results are as follows:

When these are substituted into Eqs. (29), (30), and (31), we have a system of coupled nonlinear integro-differential equations.

5. The Uniform Case

We put Vext = 0, assume V(,) = V( − ), and seek solutions with uniform density by putting where x, y and ϕ0 can be taken to be real. The condition [Eq. (26)] becomes , which can be satisfied by putting

This parametrization simplifies the calculations. For example, the liquid-structure factor of the ground state can be represented in the form where n0 is the condensate density given through

The Bogoliubov result is recovered by neglecting the sums, which is equivalent to putting n0 ≈ n: where ϵ = ℏ2k2/2m, and V is the Fourier transform of the interaction potential.

6. Reduction to δ-Function Potential

The equations in the last section are rather complex, and difficult to solve even numerically. To gain some insight, we reduce them to a simpler form by using a δ-function potential. We shall show that this potential has no effect if self-consistency is enforced, at least in the case of a uniform system. We put where

The functions ρ and Δ now depend only on :

Note that ρ is the depletion of the condensate due to interactions. Equations (29), (30), and (31) reduce to where we have suppressed the dependence of the functions. The dilute uniform hard-sphere gas is recovered by setting ρ = Δ = 0, for they give higher-order contributions.

For ρ = Δ = 0, Eq. (45) is the “nonlinear Schrödinger equation” [21] [22], which has been widely used in numerical studies of the condensate and its excitations [23]. One should be aware, however, that ρ and Δ may be important, especially for the excitations.

In the uniform case, with Vext = 0, we obtain, in the notation of the last section, where = ℏ2k2/2m, and the chemical potential μ has been eliminated in terms of the density n. The self-consistency conditions [Eq. (44)] then lead to where

The integrals are cut off at Λ, because otherwise Δ would be divergent. Solving for Δ from the second equation, we obtain in the limit Λ → ∞

The first equation then gives, in the limit Λ → ∞,

That is, there is no depletion of the condensate. It is then straightforward to show that the system is an ideal Bose gas. This result is reassuring; it shows that the method is able to verify that the δ-function potential has no effect.

7. Self-Consistent Field Method: Finite Temperatures

To extend the variational approach to finite temperatures, we can use the eigenstates of Heff to calculate the partition function, and then minimize the free energy where the brackets now denote thermal average with respect to Heff, and S is the entropy. De Gennes [19] argues that the quantity is stationary, and implements the variational principle by minimizing F − Feff. But the assertion is correct only if S were defined with respect to the eigenvalue spectrum of Heff. In fact, as we explain below, S is defined with respect to the energy levels 〈α|H|α〉, where |α〉 are the eigenstates of Heff. The finite-temperature results in Ref. [19] are therefore incorrect.

A proper variational principle for the free energy [24] is based on the inequality

This shows that, to calculate the trial partition function, we must use the energy spectrum

In contrast, the other scheme leads to a free Bose gas with single-particle energy ϵ, the eigenvalues from Eqs. (30) and (31). Because of the simplicity, it is widely used in the literature; but it does not follow from a variational principle.

As an example of the difference between the two schemes, the calculation of the partition of a hard-sphere Bose gas in Ref. [25] corresponds to using Eq. (55), and yields the virial coefficients in the gas phase. On the other hand, a calculation based on the other scheme would give a free Bose gas above the transition temperature.

8. Gaussian Variational Method

We go to a representation in which the field operator Ψ () and its conjugate have the forms where φ () is a c-number function. The more familiar representation, in which Ψ is diagonal, would be awkward to use in a non-relativistic case, in which Ψ and iΨ† are canoncial conjugates. (A similar situation happens in a relativistic fermion field obeying the Dirac equation.)

The state of the system is described by a wave functional Ψ [φ]. We use a Gaussian form for a trial wave functional: where C is a normalization constant and φ0 and G(,) are the variational parameters. Expectation values are given by functional integrals, for example where Φ is normalized to unity. It is easy to show that

The variational principle states that

The equation obtained for G and φ0 will not be given directly. To make contact with the self-consistent field method, we quote the results

For further analysis, it is convenient to introduce it easy to show that ρ (,) = ∫ d3zY(,)Y(,), and Δ (,) = − ∫ d3zX(,)Y(,). Now expand X and Y in terms of a basis {b()}: where X and Y are real. We then find

Comparison with Eq. (36) leads to the identification

This shows the equivalence between this method and the self-consistent field method.

One of the advantages of the Gaussian method lies in the possibility of generalization to the “time-dependent Hartree-Fock” method [26,27], which has been used successfully in nuclear physics. One of us (P.T.) plans to discuss this in a separate publication.