Generalized Mean Fields for Trapped Atomic Bose-Einstein Condensates.

We describe generalized time-dependent mean-field equations for partially condensed samples of trapped and evaporatively cooled atoms. These equations give a way of investigating the various order parameters that may be present as well as the existence of a mean value of the field due to condensed atoms. Our approach provides us with a closed system of self-consistent equations for the order parameters present. The equations we derive are shown to reduce to other treatments in the literature in various limits. We also show how the equation of motion method allows us to construct a formalism that can handle the evolution of these mean fields due to two-loop kinetics.

1. Introduction to Generalized Mean Fields

A great deal of interest has been generated in Bose-Einstein Condensation (BEC) by the recent experimental observations of this phenomenon in evaporatively cooled alkali gases [1], [2], [3]. These experiments have, in particular, given great impetus to the study of inhomogeneous Bose gases. There are, of course, a fair number of treatments of the Bose gas in the literature, which make use of the various tools of the many body trade and a thorough account of the special features of the physics of homogeneous Bose-condensed liquids is given in [4]. Most treatments deal with Bose-condensed systems that are static and either in, or close to, thermal equilibrium. The majority of studies have also, until recently, been restricted to homogeneous gases, with the notable exceptions of the work of Pitaevskii [5] and Fetter [6]. However, due to the fact that BEC has been experimentally achieved in a trap, there are several more recent studies that have addressed the inhomogeneous case. In this paper we shall present an equation of motion approach [7] that is well-suited to handling both time-dependent and non-equilibrium cases for a trapped (i.e., inhomogeneous) gas. We shall also show how this approach reduces to the conventional T ≠ 0 theory in the appropriate limits.

At non-zero temperatures the atomic assembly does, of course, consist of excitations above the condensate which have to be taken properly into account. Perhaps the best known way to do this is based on the finite-temperature Hartree-Fock-Bogoliubov (HFB) equations [8]–[11]. These equations include the effects of the mean value of the atomic field, the effects of thermally excited atoms, as well as the so-called anomalous averages of atomic operators [12]. These anomalous averages represent strong pair correlations between the particles, that can arise when the effective interactions are attractive. The HFB equations can be derived in a variety of ways, e.g., via Green functions, variational principles [8], [11] and through an equation of motion approach. An important discussion of the HFB theory has been recently given by Griffin [13], who describes how the behavior of an inhomogeneous Bose-condensed gas can be obtained for the full temperature range within the context of the Popov approximation [14] (explained below).

In this paper we present what could be considered a generalization of his approach to include the possibility of yet further order parameters of the type discussed by Lalöe [15]. In doing this, we also lay the ground for a fully non-equilibrium technique. We have chosen to work with the equation of motion method, as opposed to other approaches like the non-equilibrium Green function technique (Keldysh Method) [16], [17], as it is somewhat easier to produce kinetic equations and also to give a link to simpler equilibrium theories. Furthermore, the equation of motion approach makes it easier to look at “one-loop” effects away from equilibrium with arbitrary initial conditions. We emphasize that our techniques are aimed at the dilute gases that are being produced by the evaporative cooling method [18]. For those, a kinetic equation approach should be a very accurate and convenient technique. It is also possible to set up efficient simulation techniques based on Monte-Carlo wavefunction methods for a generalized density matrix approach.

We still need, however, to construct a proper self-consistent starting point for the kinetics and employ a consistent decoupling approximation. We believe this to be equivalent to the approach advocated in [19] for the calculation of the time variation of the superconducting order parameter in a Fermi liquid. Central to such an approach lies the idea that the main effects of the excitations can be described through the mean fields (“one-loop” diagrams), whereas all higher order correlation effects can be summed up into kinetic terms terms (i.e., “two-loop” effects). We shall not be able to give a full account of the kinetic terms here, but we will show clearly how they enter our formulation and how they can be computed. The full details of the kinetic terms will be given elsewhere [20].

We start our analysis by considering a group of atoms trapped by an external field. The spatial and temporal behavior can both be fully described, in second quantization, in terms of the boson field operator . We will first outline our approach in terms of , before giving our results for the occupation number representation in the subsequent sections, the latter treatment being better suited to a computational approach to the problem. For the dilute gases under consideration, we assume that the dominant interatomic interactions occur pairwise. The hamiltonian of the system can be thus written in the form Here is the interaction between the particles, which, for the dilute bose gas, is conventionally taken to be Here U0 is related to the scattering length a via [21]. The potential problem of applying this in a self-consistent approach in a 3-dimensional gas is described by Huang et al. [22]. In fact our approach does not rely on this approximation, as will become clear in Sec. 2. For the point of discussion, let us however assume, for the moment, that it is valid; this will enable us to give our initial discussion in terms that link easily with previous accounts.

The operators satisfy the equal-time boson commutation relations The Heisenberg equation of motion for then becomes This exact operator equation is the starting point for most quantum-mechanical treatments of a Bose system. We are, in the main, interested in investigating the temporal variation of the coherent (mean-field) relative to the incoherent parts of the boson field operator. We thus follow the usual route of decomposing as [6] Here expresses the non-vanishing macroscopic value of the Bose field associated with broken symmetry due to the presence of BEC. We note that there are some subtle extra issues that arise in this procedure in the case of small trapped assemblies of atoms. For the point of view of most of this article we shall assume that this is a well-defined quantity. The reader should however bear in mind that in small systems (, t) may evolve in amplitude and phase [23], [24]. In this paper we shall develop a method for determining this evolution.

Using this decomposition we note that we can re-express the interaction term of Eq. (5) in the form

The first term in the above expression corresponds to the mean field term which, if it were the only term retained, would convert the mean value of Eq. (5) into the usual Nonlinear Schrödinger Equation (NLSE) or Gross-Pitaevskii equation [25] This is the simplest mean field equation describing the atomic assembly and is only strictly valid when all the atoms are condensed, e.g., at temperatures close to zero. The NLSE has already been studied in some detail along with the link to the T = 0 elementary excitations (phonons) in isotropic [26] – [29] and anisotropic [30], [31] traps.

The next four terms appearing on the right hand side of Eq. (7) represent the interaction between condensate and excited states. Let us now discuss their meaning. Consider substituting those terms into the operator equation of motion Eq. (5) and taking the average of that expression. The first two terms in the resulting mean-field equation correspond to the direct and exchange Hartree-Fock (HF) terms. The next two terms represent the Bogoliubov correction. In the Popov approximation [13], [14] of the HFB equations, the anomalous averages of the form are assumed to be negligible, so that the second of the Bogoliubov corrections is dropped. This approximation can be shown to give accurate results for a homogeneous gas with purely repulsive interactions.

Let us now substitute Eq. (7) in full (i.e., including the term) into the operator Eq. (5) and consider its mean value. This gives the exact equation of motion Conventionally, the terms are handled in the self-consistent approximation [13] This implies that such terms vanish trivially in the finite temperature generalization of the NLSE, given by Eq. (9). However, the terms contain the effects of more complex correlations between the particles, which are neglected in the HFB approximation. We have chosen to analyze the effect of these correlations, thus extending the pure mean-field treatments (e.g., Ref. [13]). Such correlations of three particle operators are transiently produced during collisions and thus generate the collisional feed term for the condensate, as discussed in [32], [33]. For the dilute gas, these kinetic effects can be dealt with, as discussed below and in [20]. We note once again that we have chosen to do the calculation in this way, because such treatment gives a full equation of motion that is valid outside thermal equilibrium and can handle essentially arbitrary initial states. The Keldysh non-equilibrium approach [16], [17] should clearly produce results entirely equivalent to ours in the dilute gas limit. One should however bear in mind that the results from two such different formalisms may be difficult to reconcile at first sight.

Let us consider the behavior of triple product averages over time scales that are long compared to the duration of a typical binary atomic collision. A self-consistent evolution may then produce non-vanishing mean values of the triples, which would mean they have to be considered on an equal footing with the other order parameters of the system (i.e., BEC or pairing). Such a possibility has already been discussed in [15], where a cluster expansion approach is used, instead of the usual perturbative expansion in the interatomic potential, for the description of dilute degenerate gases with short-range repulsive interactions. In this paper we derive a self-consistent set of equations that allows for the inclusion of possible non-vanishing triple product mean values, which we shall henceforth refer to as “triplets” in the spirit of Lalöe [15]. Whether such an order parameter can be present in a real system and how it might affect the behavior of the system will be the subject for future study.

We shall therefore procede in a generalized mean-field approach for which triplets are retained as a possibility. We are particularly interested in laying the ground for the study of such self-consistent solutions and their role in the “macroscopic” behavior of a gas. When looking at the structure of Eq. (9) it is clear that in order to obtain a closed set of equations for the various possible order parameters, we shall need the evolution of the averages of products of two and three fluctuation operators δ(†)(, t). In our self-consistent coupled equations, all operators being averaged over have the same time label which we shall not, for the sake of convenience, give explicitly in the equations. In expressing the equations for such averages we have used suffixes 1 and 2 to indicate whether the boson field operators are acting on spatial coordinates or ′ respectively. For the binary products, corresponding to what we shall term the generalized density matrix, we find the following equations of motion here denotes a Kronecker delta.

These equations are again exact, although fairly useless without any approximations. We note that, apart from averages of products of three fluctuation operators, terms with products of four such operators appear on the right-hand side of the above equations. These represent more complex correlations between particles undergoing collisions and we will make appropriate approximations for them as we proceed.

In principle, one could work out the equations for the rate of change of both these quantities which would in turn depend on higher order correlations. It is, however, essential to terminate this coupling to higher and higher order parameters by making some suitable truncation, which is a very natural procedure in a weakly-coupled gas [7]. In this paper, we have chosen to include the possibility of mean values of triple products in our model and reduce higher particle correlations to lower ones (via an approximate version of Wick’s theorem). By including all correlations up to the triplets, we can describe the effects of the kinetics on the evolution of the condensate in a fully self-consistent mean-field approach. In order to obtain the desired closed set of coupled equations within this framework, we need of course to consider the equations of motion for averages of products of three fluctuation operators, which is done in Sec. 4.

The method described in this paper can, in principle, be used for calculating the relative importance of these various order parameters. The possible appearance of a pairing or other order parameters in the case of attractive interactions has been anticipated for some time [34]. In [35] Stoof has discussed the possibility of the anomalous average becoming the important order parameter in the case of a homogeneous gas with attractive interactions. We would like to examine the closely related issue for a trapped inhomogeneous gas.

In the next section we repeat the above analysis in the occupation number formalism and discuss the various approximations made in more detail. We then show how the equations can be cast in a generalized density matrix form (Sec. 3) and derive equations of motion for the triplets (Sec. 4). We further illustrate the consistency of the derived equations with simpler ones found in the literature (Sec. 5) and argue how our analysis can be used to set the basis for a full non-equilibrium theory.

2. Mean Field Equations for a Trapped Atomic System

In the occupation number representation the hamiltonian for pairwise interactions has the form In this equation the operator contains the kinetic energy as well as the trap potential and represents the interaction potential between a pair of particles. This is defined by Here the ψ are the spatial parts of the Fock space decomposition of given by Analogous to Eq. (6) we now decompose the time-dependent single-particle operators according to This decomposition is extremely useful as it allows the application of Wick’s theorem for decomposing higher order correlations into lower order products.

Initially we consider the temporal dependence of the mean field which yields the exact equation1 [analogous to Eq. (9)]

Here we have defined the one-body density matrix (t) and the two-body or pair density matrix (t) by their respective matrix elements2 We note that (t) and (t) are both n × n matrices, where n is the number of single-particle states.

The equation of motion for the single-particle density matrix (t) is given by with We now proceed to decompose operators a(†) (t) according to Eq. (16). Hence (upon suppressing temporal dependence) the exact form of Eq. (20) contains terms with the following structure (with the appropriate in-dices): and their hermitian conjugate quantities. At this point we will simplify our expressions by assuming that quantities of the form 〈c†c†cc〉 can be approximated by their mean-field values according to the conventional form of Wick’s theorem In this approximation, we are essentially assuming that the effects of higher order correlations (such as 〈c†c†c†cc〉) on the mean values of 〈c†c†cc〉 can be neglected. Within our mean-field approximation defined by Eq. (22), we obtain the following equations of motion Here we have defined the time-dependent Hartree-Fock hamiltonian and the pairing field by their respective matrix elements [11] as follows:

We note that in an alternative (variational) approach [11] (see Appendix A), these may be defined as where is the energy functional of the system.

The equations of motion for the pair matrix can be derived along the same lines as above. Consideration of the term 〈[Ĥ, a (t)a (t)]〉 gives rise to the correlation (where r, m, n represent indices that are being summed over). By using the decomposition equation [Eq. (16)] for obtaining correlations of fluctuation operators and imposing the boson commutation rules for the c-operators, we obtain the term , which is handled in the mean-field approximation in an analogous fashion to Eq. (22).

Hence, the equation analogous to Eq. (12) reads The parameters z, , and 〈c(†)cc〉 can be treated (self-consistently) as suitable order parameters for the system, provided they are well-defined. This could be true either in a completely static or a slowly-varying situation. It is however also possible that the kinetic correlation terms 〈c(†)cc〉 are established on time-scales much shorter than the mean evolution times of the parameters z, , and . This is equivalent to saying that the collision duration relevant to the evolution of 〈c(†)cc〉 must be much smaller than the mean time between collisions. This should hold for dilute weakly-interacting Bose gases for which nta3 ≪ 1 where nt is the number density in the trap.

In the next section we write these equations in a more compact generalized matrix form. We then derive equations of motion for the triplets (Sec. 4) and show that if the latter are neglected (Sec. 5), the equations derived in this paper are consistent with simplified versions appearing in the literature.

3. The Equations of Motion in a Generalized Density Matrix Form

Equations (23) and (27) can be cast in a generalized matrix form where the matrices M and N have been defined in terms of their elements and represents the transpose of matrix N.

It is straightforward to show that these can be written in matrix notation as where all matrices are of the form 2n×2n. Here we have defined the generalized boson density matrix R by and the matrix η by where 11 is the unity matrix of the n-dimensional Fock space.

For convenience, we have also defined the boson quasiparticle hamiltonian We note that the matrices R and H are not the same as those defined by Blaizot and Ripka [11] as discussed in Appendix A.

We have further defined the matrix which gives the effect of the triplets on the single-particle and pair excitations.

Thus far, we have derived equations of motions for and which involve triplets. In order to get a closed set of mean values, we need the equations of motion for the triplets which we derive below.

4. Equations of Motion for the Triplets

Before deriving the equations of motion for triplets of the form 〈c(†)cc〉 we would like to point out that they can be used in two ways. Firstly, they allow us to construct a self-consistent generalized mean-field theory. Furthermore, they offer the possibility of including kinetic effects of condensate formation.

We define the time-dependent n×n×n rank 3 triplet tensors and by their respective elements and

The equations of motion for the triplets can be obtained by following the procedure of Sec. 2. We note that this leads to very complex expressions, so the account will be kept very brief here.

Consider initially the equation of motion fot the triplet γ. Upon imposing the decomposition equation [Eq. (16)], we note that the right hand side of the resulting equation contains—apart from the simpler correlations—correlations of four and five single-particle fluctuation operators. We assume that the effects of correlations of six or more operators (present on the right hand side of higher equations of motion) are negligible when looking at the behavior of the system during a typical collision. This means that we set all terms 〈cccccc〉 ≃ (〈cc〉)3, (〈ccc〉)2 (appearing in the equations of motion) to zero, irrespective of their indices or whether they contain products of normal or anomalous averages. Correlations of four fluctuation operators are handled in the mean-field approximation Eq. (22). In treating the equation of motion for we also use Wick’s theorem to decompose correlations of five fluctuation operators in the form

It is essential that the above decomposition is in accord with the mean field approximation defined by Eq. (22). Strictly speaking Eq. (22) should contain an extra term on the right hand side due to the effects of five particle operator correlations, and these would have to be included in a treatment that went beyond the mean-field approach.

Writing the resulting equation of motion in more compact notation we obtain where we have defined the following quantities

Similarly we have

There are a range of ways in which we can make use of these equations. Firstly, we could solve Eqs. (17), (23), (27), (40), and (46) self-consistently as a generalization of the HFB equations that includes the extra anomalous averages. Exploring the phase diagram for the various order parameters will, however, in itself be an enormous task. It will, nonetheless, be an essential one in the study of these new evaporatively cooled assemblies. In the first instance an examination of the HFB-Popov behaviour will be an important advance. The next step should be the study of how 〈cc〉 and 〈c(†)cc〉 arise in the presence of attractive interactions. We believe that this will be an important theme in the next few years. We can also use these equations to investigate the damping of time dependent solutions that come from one-loop processes, e.g., the Landau damping discussed by Payne and Griffin [36] for the homogenous case. This means we can explore the range of validity of the simpler NLSE approaches to condensate time evolution.

These equations can also be used to treat condensation kinetics. It should however be pointed out that a completely general treatment that also includes kinetics of excited states, further requires us to look at the equations of motion of terms 〈c(†)c (†)cc〉, where the algebraic complexity becomes quite formidable.

5. Reduction to Simpler Equations

We now neglect the triplets and obtain simpler equations in appropriate limits, some of which can be found in the literature.

5.1 The Nonlinear Schrödinger Equation

If we further ignore single-particle and pair excitations in the time-dependent equation for z (t) [Eq. (17)] and restore explicit time-dependence we end up with the mean-field equation where the time-dependent mean-field Hartree-Fock hamiltonian (0)(t) is defined by One should note that the interaction potential appearing in this equation is the original one. The replacement of this by the two-body T-matrix is straightforward but relies on a decoupling—ladder approximation—in the equations of motion. We shall give the explicit account of this in a future publication. With this replacement Eq. (51) becomes the analogue of the NLSE [Eq. (8)]. in the occupation number representation we note that Eq. (15) can be equivalently written in the form This confirms that we can reconstruct the NLSE by multiplying Eq. (51) by and summing over all indices i. The factor has been included so as to re-scale the time-dependent coefficients z(t) which obey the normalization condition At T = 0 the whole system is condensed in the lowest energy eigenstate and since according to Eq. (54) all the time-dependence has been included in the coefficients z (t), we can re-express these as where the z are now simply constant complex numbers and µ is the ground state chemical potential. Clearly Eq. (51) then reduces to the time-independent NLSE [11] Here we see that we can trivially transform from the time-dependent to the time-independent picture, merely by setting the time-dependent term on the left-hand side to zero and simultaneously replacing the time-dependent mean-field Hartree-Fock hamiltonian by the time-independent one given by We note that this remains true even at non-zero temperatures, when the second term also contains correlations due to the excited particles.

5.2 Static Equations

In this section we will show that Eqs. (23) and (27) are also consistent with the T > 0 static HFB equations that we have additionally derived from a variational principle approach in Appendix A. Neglecting all time dependence in Eqs. (28) and (29), we thus obtain the static matrix equations However, to solve these equations self-consistently, we also need an equation normalizing the vector . We use the substitution in Eq. (17) to obtain the static equation where the chemical potential µ appears as discussed earlier.

This can be re-expressed as where and

More generally, we obtain the matrix equations [11] which give us the values of all z for fixed values of h and ∆. Equations (60), (61) and (66) constitute the full set of static Hartree-Fock-Bogoliubov equations which can be solved self-consistently. If we set ∆ = = 0, we thus obtain the static self-consistent HFB description of a condensate with elementary excitations in the Popov approximation [14], which further reduces to the time-independent NLSE for = 0.

6. Discussion

In this paper we have developed a time-dependent self-consistent generalized mean-field method for describing a trapped and partially condensed atomic system. Central to our approach lies the operator decomposition equation [Eq. (6)] (or equivalently Eq. (16)) into a mean field and a fluctuating part. This decomposition is reasonable, as we do not make any assumptions about the initial value of the wavefunction Φ (, t), or its evolution. Consequently, our treatment should remain valid over a wide range of conditions, e.g., for cases close to the transition temperature Tc. Our equations allow for the possibility of other order parameters and give us a way of investigating their possible role in an evaporatively cooled gas.

In our treatment, we have considered an inhomogeneous gas, which implies that the excited states of the system are strongly influenced by the finite size of the trapped assembly. This natural length scale will, of course, lead to the presence of a gap in the excitation spectrum, as opposed to the gapless spectrum of the homogeneous gas.

Let us initially assume that the triplets are small in comparison with the usual HFB order parameters. In the case of moderately-sized condensates, the solutions of the HFB equations will hence be a good basis, as long as the kinetic processes produce collisional widths that are less than the separation of the levels. However, as the number of atoms in the trap is increased, we will approach the so-called quasi-homogeneous limit for the condensate. This happens when the healing length for the gas is much smaller than the size of the condensate. In that regime, it becomes essential to describe the mean fields by means of local densities [i.e., z(r), (r), and (r)]. We thus obtain locally well-defined quasiparticles (in the sense that their collisional widths are smaller than the characteristic frequencies) in regions specified by the healing length of the atomic assembly under consideration. This approach has been employed by Dorfman and Kirkpatrick to derive transport coefficients for the dilute Bose gas [32].

The discussion above relies on assumption that the HFB states for the trap as a whole constitute the most sensible state. This, in turn, implies that we should, in order to be consistent, assume that any kinetic rates are less than the frequency separation of the HFB excitations, which seems reasonable for the case of condensates with modest numbers of particles. Our equations are, however, not limited to the pure HFB regime. For example, in the case of repulsive interactions, both pairing and triplets may be treated purturbatively and the analysis extended beyond HFB.