# Second-order, number-conserving description of non-equilibrium dynamics in finite-temperature Bose-Einstein condensates

###### Abstract

While the Gross–Pitaevskii equation is well-established as the canonical dynamical description of atomic Bose-Einstein condensates (BECs) at zero-temperature, describing the dynamics of BECs at finite temperatures remains a difficult theoretical problem, particularly when considering low-temperature, non-equilibrium systems in which depletion of the condensate occurs dynamically as a result of external driving. In this paper, we describe a fully time-dependent numerical implementation of a second-order, number-conserving description of finite-temperature BEC dynamics. This description consists of equations of motion describing the coupled dynamics of the condensate and non-condensate fractions in a self-consistent manner, and is ideally suited for the study of low-temperature, non-equilibrium, driven systems. The -kicked-rotor BEC provides a prototypical example of such a system, and we demonstrate the efficacy of our numerical implementation by investigating its dynamics at finite temperature. We demonstrate that the qualitative features of the system dynamics at zero temperature are generally preserved at finite temperatures, and predict a quantitative finite-temperature shift of resonance frequencies which would be relevant for, and could be verified by, future experiments.

###### pacs:

03.75.Kk 67.85.De 05.30.Jp## I Introduction

Even at zero temperature, typical atomic Bose-Einstein condensates (BECs) contain a small non-condensate fraction due to inter-atomic interactions. In real experiments, which necessarily take place at finite temperatures and may involve dynamical depletion of the condensate Gardiner (2002) — for instance, due to non-equilibrium dynamics induced by changes in applied external fields, a non-negligible non-condensate fraction often arises. A full understanding of the dynamics of such systems requires a fully-dynamical, finite-temperature theoretical description which goes beyond the mean-field, zero-temperature Gross–Pitaevskii equation (GPE). Due to the complexity of such descriptions, the non-equilibrium dynamics of atomic BECs in the presence of a significant non-condensate fraction (whether thermal or dynamical in origin) remains a largely open problem Proukakis and Jackson (2008); Gasenzer and Pawlowski (2008).

In systems where the non-condensate fraction is primarily thermal in origin, a variety of theoretical descriptions have been developed and successfully applied to non-equilibrium atomic BEC dynamics. These include symmetry-breaking descriptions, which are based on a perturbative expansion about a mean field; for example, the Hartree–Fock–Bogoliubov–Popov (HFBP) description Hutchinson et al. (1997); Dodd et al. (1998); Reidl et al. (1999); Giorgini (1998), and the Zaremba–Nikuni–Griffin (ZNG) description Zaremba et al. (1999); Griffin et al. (2009); Jackson and Zaremba (2001, 2002); Williams et al. (2002); Jackson et al. (2009, 2007). Other successful descriptions have been obtained in the context of -field methods, reviewed in Blakie et al. (2008), which describe the highly occupied modes of the system as a classical field. These descriptions include the projected Gross–Pitaevskii equation (PGPE) Davis et al. (2001); Davis and Blakie (2006); Simula and Blakie (2006), the stochastic projected GPE (SPGPE) Gardiner et al. (2002); Bradley et al. (2008) and stochastic GPE (SGPE) Stoof (1997); Proukakis and Jackson (2008); Cockburn and Proukakis (2009); Cockburn et al. (2010); Weiler et al. (2008); Neely et al. (2010); Fialko et al. (2012), and the truncated Wigner PGPE Lobo et al. (2004); Norrie et al. (2006); Scott et al. (2008); Martin and Ruostekoski (2012).

However, the non-equilibrium dynamics of systems in which a low-temperature BEC is driven by an applied external field leading to dynamical depletion of the condensate have not been as widely investigated. The prime example of such systems are atomic BEC analogs of generic quantum chaotic systems; e.g., the kicked accelerator Ma et al. (2004); Schlunk et al. (2003a), kicked harmonic oscillator Billam and Gardiner (2009); Gardiner et al. (1997, 2000), and kicked rotor Zhang et al. (2004); Liu et al. (2006); Shepelyansky (1993); Reslen et al. (2008); Monteiro et al. (2009); Billam and Gardiner (2012). These systems offer an excellent test bed for exploring generic issues of quantum chaos Schlunk et al. (2003b); Shepelyansky (1993), quantum superposition Weiss and Teichmann (2008), quantum resonances Billam and Gardiner (2009); Reslen et al. (2008); Monteiro et al. (2009), dynamical instability and dynamical depletion Gardiner et al. (2000); Monteiro et al. (2009); Zhang et al. (2004); Reslen et al. (2008), and even entropy, thermalization and integrability Polkovnikov et al. (2011); Yukalov (2011); Yurovsky and Olshanii (2011). A particular issue in the description of these systems is the occurence of nonlinear quantum resonances Monteiro et al. (2009): values of the system parameters at which the system resonantly absorbs energy from the driving. Such resonances are an extremely generic feature of these systems, and are sensitive to the value of the nonlinearity (product of -wave scattering length and condensate occupation ). Consequently, one expects the dynamical depletion of the condensate caused by the driving to rapidly and fundamentally alter the subsequent dynamics of the system close to these resonances, compared to the predictions of the Gross–Pitaevskii equation Gardiner (2002).

The rapid and severe back-action of condensate depletion on the non-equilibrium dynamics of these systems presents considerable challenges for the above-mentioned finite-temperature descriptions. In particular, the PGPE and SPGPE are fundamentally restricted to a high-temperature regime, and lack a quantum treatment of the pair-excitation processes which drive condensate depletion. In principle the truncated-Wigner PGPE may be used to model non-equilibrium dynamics of driven systems at low (and zero) temperatures Martin and Ruostekoski (ress). However, in order to go beyond a close-to-equilibrium approximation in which the condensate remains in its initial state, one must obtain the dynamics of the condensate within such a treatment from an ensemble average of individual GPE trajectories. While such a treatment is entirely possible, and can for example be conducted by determining the single-particle density matrix in a number-conserving fashion Wright et al. (2011), the issues of spurious thermalization and inaccurate long-time dynamics in the truncated Wigner method (see, e.g., Refs. Blakie et al. (2008); Sinatra et al. (2002)) would necessitate a very careful approach when applying such a treatment to the non-equilibrium dynamics of driven systems.

In order to comprehensively describe such systems one must thus self-consistently capture the back action of the increasing population of, and the subsequent dynamics in, low-lying non-condensate excitations on the condensate. To achieve this, a description which consistently divides the system into well-defined condensate and non-condensate fractions at all times would appear to be necessary (as opposed to a -field description in which the condensate fraction must be extracted by subsequent ensemble averaging). In that they comprise a set of coupled equations for the condensate and the non-condensate components, symmetry-breaking descriptions would thus appear to offer a more appropriate treatment of low-temperature driven systems. However, all such dynamical descriptions which have been applied to BEC dynamics to date have contained an approximate treatment of pair-excitation processes that neglects the anomalous average; while suitable at higher temperatures, this leads to descriptions which incompletely model the phonon character of elementary excitations at low temperatures Allen et al. (ress). Furthermore, the assumption of a symmetry-broken condensate state, and hence overall number-non-conservation, in any symmetry-breaking method presents an additional potential problem; since the prime contribution to the departure from GPE dynamics in low-temperature driven systems results from, and can be extremely sensitive to, changes in the number of condensate particles, conservation of the total atom number appears to be necessary in order to obtain a fully self-consistent description of such systems as realized in atomic BEC experiments (in which the total atom number is finite and fixed) Gardiner and Morgan (2007). In a formal sense, coupled equations of motion for the condensate and non-condensate within a symmetry-breaking treatment do not contain any terms to explicitly preserve the orthogonality of the condensate and non-condensate modes (see Section II.1); in the context of low-temperature driven systems, this means an ambiguity in interpreting the solutions of those equations for anything other than short times.

In view of the above considerations number-conserving descriptions Castin and Dum (1998); Gardiner (1997); Girardeau and Arnowitt (1959); Girardeau (1998); Morgan (2000), in which the system is partitioned into manifestly orthogonal condensate and non-condensate parts using the Penrose-Onsager criterion Penrose and Onsager (1956), are, in principle, ideally suited to this regime. However, until recently only the first-order number-conserving description of Gardiner Gardiner (1997) and Castin and Dum Castin and Dum (1998) offered a description which was numerically tractable in fully time-dependent form. Application of this first-order number-conserving description to the -kicked-rotor BEC Zhang et al. (2004); Liu et al. (2006); Reslen et al. (2008) and -kicked-harmonic-oscillator BEC Gardiner et al. (2000); Rebuzzini et al. (2007) revealed a general tendency for rapid, unbounded growth. Unfortunately, as the equation describing the condensate in this first-order treatment is simply the GPE, this rapid, unbounded growth can be viewed as a consequence of linearized instabilities present in the GPE, which would be removed by a higher-order description containing a self-consistent back-action of the non-condensate on the condensate Gardiner (2002). Thus, similar to the methods discussed above, the first-order description is also of limited use in describing the long-time dynamics of low-temperature, driven BEC systems.

However, such a self-consistent back-action is present within a second-order number-conserving description, such as presented by Gardiner and Morgan in Ref. Gardiner and Morgan (2007) and previously used highly successfully by Morgan Morgan et al. (2003); Morgan (2004, 2005) to calculate, using a linear response treatment, the excitation frequencies of a BEC at finite temperature as measured in experiments at JILA Jin et al. (1996, 1997) and MIT Mewes et al. (1996); Stamper-Kurn et al. (1998). Recently Billam and Gardiner (2012), the first fully time-dependent implementation of the second-order number-conserving equations of motion was applied to an initially zero-temperature -kicked-rotor BEC. This application revealed that the second-order description’s self-consistent back-action does indeed damp out the unbounded growth in non-condensate fraction seen in the first-order description. Consequently, the time-dependent form of the second-order number-conserving description provides an excellent tool for studying the non-equilibrium dynamics of driven BECs at low temperatures.

The purpose of this paper is primarily to present, in detail, analytic and numerical techniques which allow one to evolve the second-order number-conserving equations of motion at finite temperatures, enabling a systematic exploration of non-equilibrium BEC dynamics at finite temperatures. A restricted subset of these techniques, limited to zero-temperature initial conditions, were used to obtain the results in Ref. Billam and Gardiner (2012), but were not described in detail in that work. We also use these techniques to conduct a finite-temperature study of the non-equilibrium dynamics of the -kicked-rotor BEC, extending the results of Ref. Billam and Gardiner (2012) and previous works. Remarkably, we observe that many of the sharp nonlinear resonance resonance features observed at zero temperature in Refs. Billam and Gardiner (2012); Monteiro et al. (2009) are qualitatively preserved at temperatures sufficient to cause significant initial condensate depletion. In particular, our exploration of the effects of initial temperature predicts a finite-temperature shift in the system’s resonance frequencies which could be experimentally measured and verified.

The remainder of the manuscript is structured as follows: We begin by introducing, in Section II, the second-order number-conserving description of the dynamics and present the second-order equations of motion in their most general form. We give a detailed discussion of the self-consistent properties and regime of validity of the resulting description at both zero and finite temperatures. In particular, we discuss potential problems finding equilibrium initial conditions at higher temperatures. In Section III we develop a numerical method to evolve initial states according to the second-order equations of motion. This method explicitly includes the nonlinear, non-local terms coupling the condensate and the quasiparticle modes. It is these terms which maintain orthogonality between the condensate and non-condensate, and which are particularly difficult to deal with in the context of general numerical integration schemes. In Section IV we apply the method to systematically explore the effects of finite initial temperatures in the -kicked-rotor BEC. We also extend the analysis of this system given in Ref. Billam and Gardiner (2012) by further exploring the contrast between first- and second-order descriptions, and analyzing the evolution of the single-particle von Neumann entropy. We predict a shift of resonance frequencies at finite temperature which could feasibly be detected in experiments. Section V comprises the conclusions. A technical appendix follows, in which we present details of our numerical technique.

## Ii Second-order, number-conserving theoretical description

### ii.1 Motivation

We consider a system of bosonic atoms of mass , confined by an external potential , and interacting via pair-wise -wave contact interactions. The Hamiltonian for such a system is given by

(1) |

where and are second-quantized field operators obeying standard equal-time bosonic commutation relations. Here, the single-particle Hamiltonian is

(2) |

where

(3) |

for -wave scattering length .

Computing the full many-body dynamics of such a system directly from the
Hamiltonian Eq. (1) is, in
general, an intractable problem. However, in the case where is large and
the majority of the system is Bose-Einstein condensed, one can obtain an
approximation to the full many-body dynamics via a perturbative description, in
which the non-condensate fraction constitutes the small parameter. To develop
such a description, we adopt the definition of Bose-Einstein condensation —
applicable to a finite-size and interacting system as described by
Eq. (1) — given by Penrose and
Onsager Penrose and Onsager (1956). In this definition the condensate mode,
is identified as a single macroscopically-occupied eigenstate
of the single-particle density matrix^{1}^{1}1Note that in
Eq. (4) the brackets denote an
expectation value with respect to the full many-particle density matrix; at
finite-temperature this involves a thermal, as well as a quantum, average.

(4) |

Consequently, the condensate mode satisfies

(5) |

where the condensate mode occupation, , is taken to be much greater than all other single-particle mode occupations [given by the other eigenvalues of ] at any time.

Using this definition, one can explicitly partition the field operator into a condensate part, , and a non-condensate part, :

(6) |

Here, the condensate annihilation operator annihilates an atom in the condensate mode . To maintain hermiticity of the single-particle density matrix one must ensure that the part of the field operator describing the condensate mode, , is explicitly orthogonal to the part describing the non-condensate, , at all times. Formally, one defines

(7) |

where

(8) |

and hence one sees that and are “orthogonal” in the sense that

(9) |

This partition is number-conserving, in the sense that the system remains in a state of fixed total atom number.

This number-conserving partition can be contrasted the symmetry-breaking partition used in other finite temperature descriptions, which is based on a partition of the field operator into a finite expectation value (or order parameter), , and an operator-valued fluctuation ,

(10) |

Two related consequences of the symmetry-breaking partition [Eq. (10)] are a breaking of the overall phase-symmetry of the Hamiltonian, and a general lack of orthogonality between the condensate and non-condensate parts of the system. Specifically, having constructed a perturbative description in terms of and one generically finds that, in the case of an inhomogeneous order parameter, and are not orthogonal in the sense of Eq. (9). As a result, one cannot uniquely reassemble a Hermitian density matrix for the system from and , as one can from and . This lack of an explicit, time-independent orthogonality between condensate and non-condensate [in the sense of Eq. (9)] also means that one cannot unambiguously determine the condensate population, , and the non-condensate population, , associated with the symmetry-breaking partition. In the context of driven low-temperature BECs, where the dynamics are extremely strongly influenced by changes in , this ambiguity in the condensate population constitutes a genuine disadvantage of the symmetry-breaking approach.

A number-conserving description of the system’s dynamics is obtained by treating the non-condensate part perturbatively, by means of an expansion in terms of a suitable fluctuation operator Gardiner (1997); Castin and Dum (1998); Gardiner and Morgan (2007). Neglecting the non-condensate completely yields a zeroth-order description in the form of the Gross–Pitaevskii equation (GPE); this takes the form [we explicitly denote the condensate mode in this zeroth-order description by for clarity in later sections]

(11) |

where , and is given by

(12) |

As in Refs. Billam and Gardiner (2012); Gardiner and Morgan (2007); Gardiner (1997), is defined in such a way as to be generally time-dependent. However, as is explicitly real, the dynamics resulting from its time evolution consist only of an irrelevant global phase. When considering, as we do in this paper, the evolution of a stationary, equilibrium initial condition subject to a time-dependent perturbation, it is most convenient to work with a time-independent GPE eigenvalue given by

(13) |

where represents the stationary, equilibrium state of the GPE. Since several subtly different eigenvalues appear in the subsequent development, here we have introduced the convention that the subscript index to denotes the order (with respect to the non-condensate fluctuation operators) of the functional defining , while the bracketed superscript index to denotes the order of approximation of the condensate wavefunction appearing inside the functional definition. Furthermore, if any eigenvalue is time dependent we always denote this explicitly with a argument. In the case where a time argument is absent, the corresponding eigenvalue should be understood as having been evaluated for the appropriate equilibrium initial condition.

To work with the time-independent eigenvalue in the GPE description, one simply makes the replacement in Eq. (11). Note that while arises naturally as a nonlinear eigenvalue during the development of a zeroth-order number-conserving equation of motion, at this level of approximation it is equivalent to the chemical potential, , which would be introduced as a Lagrange multiplier to determine the average particle number in a symmetry-breaking approach.

In contrast to the GPE, the second-order number-conserving description we use in this paper provides a description of both the condensate and the non-condensate, and consists of mutually-coupled equations for both, which we outline in the following section.

### ii.2 Equations of motion

Conducting a self-consistent, second-order, number-conserving expansion — as detailed in Ref. Gardiner and Morgan (2007)— leads to a number-conserving generalized GPE (GGPE) for the dynamics of the condensate mode :

(14) |

Here we have introduced the convention, adopted generally hereafter, of omitting explicit time-dependences wherever this aids clarity. The dynamics of the non-condensate enter the GGPE [Eq. (14)] through the normal [] and anomalous [] averages. These are respectively defined by

(15) |

and

(16) |

where

(17) |

is the number-conserving fluctuation operator in which a perturbative expansion has been conducted. Note that the corresponding small parameter is , where is the occupation of the non-condensate. We have also introduced the complex, time-dependent “GGPE eigenvalue” , which is given by

(18) |

It is essential to emphasize that this “eigenvalue” is in general time-dependent, and that Eq. (18) should be taken to be evaluated using the values of , , etc. at a particular time. As we demonstrate in the subsequent section, this time dependence of is crucial in order to obtain self-consistent number dynamics. However, as in the case of the GPE [Eqs. (11) and (13)] it is considerably more convenient to work with a real, time-independent eigenvalue, and in Section III we will reformulate the second-order equations of motion in order to consistently (up to second order) replace with such a quantity.

The GGPE must be coupled to a consistent set of equations of motion for the non-condensate; at second-order these consist of the number-conserving modified Bogoliubov–de Gennes equations (MBdGE):

(19) |

where

(20) |

and

(21) |

The modification in these equations — with respect to the ordinary BdG equations obtained in a symmetry-breaking description — is the appearance of the projector which explicitly enforces the orthogonality of the condensate and non-condensate, and of the time-dependent “GPE eigenvalue” . This “GPE eigenvalue” is obtained by substituting the GGPE condensate mode into the GPE [Eq. (11)] in place of the GPE condensate mode , and allowing the eigenvalue to be generically time-dependent. That is,

(22) |

where we have resurrected time-arguments for clarity. As in the zeroth-order, GPE, case — where it was possible, and convenient for the case of a stationary, equilibrium initial state, to replace the time-dependent GPE eigenvalue with a time-independent GPE eigenvalue — is an explicitly real quantity. It will thus later be convenient, for the case of a stationary, equilibrium initial state, to replace the time-dependent “GPE eigenvalue” with the time-independent “GPE eigenvalue” given by

(23) |

Note, however, that the replacement of with must be accomplished consistently with the replacement of with a real, time-independent quantity as discussed previously. The details of this consistent replacement are outlined in Section III.

The generalized Gross–Pitaevskii equation [Eq. (14)] and the modified Bogoliubov–de Gennes equations [Eq. (19)] complete the fully dynamical, second-order, number-conserving description obtained by Gardiner and Morgan in Ref. Gardiner and Morgan (2007) and previously used within a linear response treatment by Morgan Morgan et al. (2003); Morgan (2004, 2005). In this paper, we develop these equations into a form where their fully dynamical time evolution can be realized numerically, as in Ref. Billam and Gardiner (2012).

### ii.3 Discussion

Before we outline our method for the simultaneous numerical solution of Eqs. (14) and (19), a few comments regarding the properties of the second-order equations of motion and their regime of validity are in order. Firstly we note that, in general, the diagonal part of the anomalous average, , is ultraviolet-divergent and must be appropriately renormalized; this procedure is described in detail in Refs. Gardiner and Morgan (2007) and Morgan (2000). One exception is in the case of quasi-one-dimensional (quasi-1D) systems, such as the one we consider in Section IV, where renormalization is not necessary.

Secondly, the number-conserving equations of motion used here are derived by expanding the total Hamiltonian in powers of the number-conserving fluctuation operator [Eq. (17)]. This operator is advantageous for three primary reasons: (1) it scales proportionally to the number of non-condensate atoms, which we wish to treat as a small parameter; (2) it avoids the need to expand inverse-square-root number-operators when expanding the Hamiltonian, which is particularly convenient; and (3) it is a well-defined fluctuation operator in the sense that . It should be noted that the commutation relations of

(24) |

give rise to quasiparticles which are, in principle, only approximately bosonic. However, when restricted to the (quadratic) order of approximation we consider, the corresponding quasiparticles are indeed exactly bosonic. Hence, the favorable properties of mentioned above make it a preferable choice for developing the second-order equations of motion Morgan (2000); Morgan et al. (2003); Morgan (2004, 2005); Gardiner and Morgan (2007); Gardiner and Billam (ress).

Note also that the appearance of higher-than-second-order fluctuation terms in the equations of motion has been prevented by working within a consistent Gaussian approximation Gardiner and Morgan (2007); that is, all quadratic products of operators are required to take the form of pair averages. This constitutes a Gaussian approximation in that all higher-order moments of the fluctuation distribution are assumed to be describable in terms of and Köhler and Burnett (2002). As was explicitly demonstrated by Morgan Morgan (2000) this is consistent within the second-order number-conserving description (see also Refs. Gardiner and Morgan (2007); Proukakis and Jackson (2008)).

A key property of the second-order number-conserving equations of motion is their number-self-consistency: in contrast to and — which can both be considered as low-order approximations to the chemical potential — is a complex eigenvalue. The meaning of the imaginary part of can be understood by considering the (implicit) time dependence of , which is given (to quadratic order) by

(25) |

Thus the time-dependent, imaginary part of acts to keep the condensate mode normalized to unity despite the growth or decay of the condensate population. This illustrates the presence of number-self-consistency in the dynamical coupling between the GGPE and MBdGE. In contrast, the first-order description Gardiner (1997); Castin and Dum (1998) — which consists of the same MBdGE [Eq. (19)] coupled to the ordinary GPE [Eq. (11)] — is not in general number-self-consistent, since the condensate population is fixed. The first-order description can only be considered to be number-self-consistent when viewed as the limit of the second-order description as Gardiner and Morgan (2007). The zeroth-order description, consisting of the GPE [Eq. (11)] alone, is trivially number-self-consistent, as it ignores the growth and decay of the condensate altogether. The number-self-consistency of the zero-, first-, and second-order number-conserving descriptions are illustrated schematically in Fig. 1.

A final feature of the second-order, number-conserving equations of motion is that the non-condensate is described by a MBdGE [Eq. (19)] which does not contain pair averages of the non-condensate operators; the terms are in no way altered from the first-order description and the terms appearing in the MBdGE consist only of a GPE Hamiltonian and the “GPE eigenvalue” of Eq. (23). Importantly, however, the GPE Hamiltonian and the “GPE eigenvalue” are, in the second-order description, evaluated in terms of the second-order GGPE wavefunction . This leads to an apparent inconsistency in that the spinors and are no longer exact, zero-energy solutions of the MBdGE, as they would be for the GPE wavefunction . This has the unfortunate side-effect of making the identification of self-consistent initial conditions difficult at high temperatures, particularly in inhomogeneous systems. However, this problem can be viewed purely as a consequence of applying the theory outside of its regime of validity, as argued by Morgan in his study of condensate excitations at finite-temperature Morgan et al. (2003); Morgan (2004, 2005); in this work he demonstrated that the second-order description remains self-consistent at high temperatures (approaching ) provided one restricts oneself to a linear response treatment, in which a self-consistent equilibrium solution is not necessary. A fully dynamical treatment, which requires a self-consistent equilibrium initial condition, is thus restricted to lower temperatures (such that ) than a linear response treatment.

## Iii Fully time-dependent numerical implementation

### iii.1 Reformulation of equations of motion

#### iii.1.1 Elimination of complex, time-dependent “eigenvalue”

In this Section we develop a numerical method for evolving the combined GGPE and MBdGE system of Eqs. (14) and (19). In order to do so, it is of great convenience to eliminate the imaginary part of the GGPE “eigenvalue” [Eq. (18)]; this can be done by describing the condensate using a mode function normalized to the condensate population, . Hence, we define

(26) |

in terms of which the GGPE can be re-expressed [using Eq. (25) for the number evolution] as

(27) |

where the adjusted GGPE eigenvalue, , is given by

(28) |

This is explicitly real, and should be evaluated explicitly in terms of the stationary, equilibrium initial condition of the GGPE at .

Note that in order to cast the GGPE in the above form, with time-independent eigenvalue , we have formally moved into a frame which cancels the time-dependence of the real part of . However, transforming the MBdGE into the same frame leads to additional terms of the form (where refers to the equilibrium initial condition) appearing in the operator . Specifically, such terms appear with

(29) |

where all quantities except on the right hand side are time-dependent. Since these terms are small, and formally of the same (cubic) order with respect to the fluctuation operator as terms already omitted from the description of the non-condensate at second order Gardiner and Morgan (2007), they should be neglected to maintain a consistent description. We thus cast the MBdGE as

(30) |

a form which emphasizes the significant structural analogies between the GGPE and the MBdGE.

#### iii.1.2 Quasiparticle decomposition

In numerical studies of non-equilibrium dynamics at finite temperature one
generally wishes to begin from a finite-temperature equilibrium condition and
explore the dynamics resulting from, e.g., driving or an abrupt quench event in
a fully-time dependent way^{2}^{2}2One could in principle start from a
non-equilibrium initial condition, but we consider only the equilibrium case
here.. In the second-order number-conserving description such a thermal and
dynamical equilibrium state corresponds to a self-consistent, stationary
solution of the equations of motion in which the elementary quasiparticle
excitations of the system are populated according to the appropriate thermal
Bose distribution.

The appropriate quasiparticle basis in the second-order number-conserving description is the basis of Bogoliubov quasiparticle modes which diagonalize the stationary MBdGE Gardiner and Morgan (2007):

(31) |

In terms of the quasiparticle mode functions and , the fluctuation operator can be expanded as

(32) |

where and are quasiparticle creation and annihilation operators which, at this order, are assumed to have bosonic commutation relations (see Section II.3). Assuming all time-dependence to reside in the quasiparticle mode functions and , with the quasiparticle creation and annihilation operators and time-independent, the MBdGE [Eq. (30)] take the form

(33) |

and

(34) |

At initial thermal equilibrium, the quasiparticle creation and annihilation operators have the following pair averages:

(35) | |||

(36) |

where is the Bose-Einstein factor Gardiner and Morgan (2007); Morgan (2004, 2000)

(37) |

The term represents a finite-size correction, given by Morgan (2000, 2004); Mandl (1988)

(38) |

This leads to quasiparticle expressions for the non-condensate normal and anomalous averages [Eqs. (15) and (16)]

(39) | ||||

(40) |

We reiterate that, by choosing a scheme where all time dependence resides in the quasiparticle mode functions and , the quasiparticle populations appearing in Eqs. (39) and (40) remain fixed in time.

#### iii.1.3 Introduction of the spinor

The primary difficultly in simulating the coupled GGPE-MBdGE system numerically is the problem of orthogonalization: both equations contain terms which function to maintain orthogonality between the condensate and non-condensate. This is in contrast to the case of the first-order GPE-MBdGE system, where the GPE evolves in isolation from the MBdGE; this de-coupling in the first-order system means the evolution of the MBdGE can be computed by ignoring the projector terms in the MBdGE throughout the evolution, and simply projecting the final state orthogonally to the condensate Gardiner et al. (2000). However, if one were to similarly ignore the projectors during the evolution of the second-order GGPE-MBdGE system, one would then have to re-orthogonalize both the condensate and quasiparticle modes with respect to an unknown basis at the end of the evolution. Consequently, the projection terms, and the non-local integrals they involve, must be explicitly included in any numerical method. In this Section we develop such a method, by re-casting the GGPE and MBdGE in a form which exploits their apparent symmetries, and allows us to include the projection terms using a split-step technique.

After substituting the quasiparticle expressions for the normal and anomalous average [Eqs. (39) and (40)] into the GGPE [Eq. (27)], the GGPE can be recast in the form

(41) |

and the MBdGE [Eqs. (33) and (34)] can be recast in the form

(42) | ||||

(43) |

where

(44) | |||

(45) | |||

(46) |

and

(47) |

where . This reformulation of the problem allows one to write the coupled evolution of the condensate wavefunction and the first quasiparticle modes as a nonlinear matrix equation in a ()-dimensional spinor space:

(48) |

Here the vector is defined by

(49) |

and the operator is defined by

(50) |

In any actual calculation, this spinor space is rendered finite-dimensional by the need for a finite quasiparticle momentum cut-off . However, may, in principle be arbitrarily large.

### iii.2 Operator-splitting scheme for time-evolution

#### iii.2.1 Separation of position and momentum terms

As we have already accounted for all creation and annihilation operators
through the quasiparticle decomposition, each entry in the matrix defining
can be thought of as an operator in the first-quantized
sense^{3}^{3}3That is, an operator which could appear on the right side of a
single-particle Schrödinger equation.. From an analytic
perspective, this notation seems to achieve little more than ‘tidying’ —
abstracting away much of the detail. Importantly, however, all the operators
which are off-diagonal in the spinor space [the operators ] are
diagonal in the position representation, since they ultimately consist of a
multiplication by a spatially-varying function. In contrast, all the operators
which have off-diagonal components in the position representation [that is, the
kinetic energy operator implicitly contained in and hence
in ] appear only on the diagonal in the spinor space.

From a numerical perspective this arrangement is extremely useful, as it makes the evolution amenable to a split-step approximation. This is achieved by splitting into the sum of a term representing the linear single-particle evolution, , and a term representing nonlinear parts of the evolution, . The linear, single-particle term is defined by

(51) |

where represents the identity matrix, and we have also defined

(52) |

The nonlinear term is then defined simply by

(53) |

Note that this matrix contains diagonal entries of the form and where

(54) |

Written in this form, contains all kinetic energy terms; while these terms are not diagonal in the position representation, they do all lie on the diagonal in the spinor space. Consequently, the evolution due to the kinetic energy terms can be computed for each mode function separately. In contrast, while does contain off-diagonal elements in the spinor space, each element of is diagonal in the position representation. Consequently, the evolution due to these terms can be computed straightforwardly over a discrete spatial grid.

In the spinor space, the split-step approximation for the evolution of the system over short times is given by

(55) |

where each of the three evolution operators on the right is assumed to act instantaneously, and the symmetrization reduces the error to the order of Javanainen and Ruostekoski (2006). Since is diagonal in the spinor space, we have

(56) |

#### iii.2.2 Analytic form for position-space evolution

The spatial operator , is more problematic because it is not diagonal in the spinor space, and because it contains the nonlocal, nonlinear terms necessary to preserve orthogonality between the condensate and non-condensate. However, each of its entries is diagonal in the position representation, meaning we only need exponentiate in the spinor space in order to obtain an operator we can evaluate and use for short-time propagation at a given position: such an operator — in distinct contrast to which consists of an application of to each mode individually — couples the mode functions by acting on all modes at once. The numerical advantage of casting the evolution in terms of such an operator is that it almost entirely overcomes the problem of retaining orthogonality between the condensate and the quasiparticle modes, which poses severe difficulties for other schemes. In practice we do find that a correction for numerical round-off error is still required in order to preserve complete orthogonality over long times; this can be achieved by explicitly orthogonalizing the quasiparticle modes with respect to the condensate at the end of each time-step. However, our experience is that applying such a correction after every time-step does not cause decay of the total atom number, indicating that the need for such a correction does indeed stem from the numerical round-off error inherent in finite-precision arithmetic. Indeed, the problem of loss of orthogonality between vectors due to finite-precision arithmetic is well-known in the context of, e.g., the Gram-Schmidt process Golub and Van Loan (1996).

In principle, the evolution due to the term can be obtained, for a given , by exponentiating the matrix numerically at each spatial grid point. However, with the necessary number of floating point operations necessary to diagonalize the matrix scaling at least as Golub and Van Loan (1996), this leads to an algorithm which is too slow to be practicable. To overcome this limitation, we have computed a general analytic expression for for arbitrary : this reduces the scaling of the computational effort to (i.e., equivalent to a matrix-vector multiplication). This expression is explicitly derived in Appendix A, where we also show that, at any specific point (and time ) on the computational grid, the time evolution due to can be approximated locally, after numerical evaluation of the non-local integrals appearing in the functions over the entire position-space grid, by the following coupled equations:

(57) |

(58) |

(59) |

where

(60) |

and we have defined

(61) |