Dynamics of Hot Bose-Einstein Condensates: stochastic Ehrenfest relations for number and energy damping

Describing high-temperature Bose gases poses a long-standing theoretical challenge. We present exact stochastic Ehrenfest relations for the stochastic projected Gross-Pitaevskii equation, including both number and energy damping mechanisms, and all projector terms that arise from the energy cutoff separating system from reservoir. Analytic solutions for the center of mass position, momentum, and their two-time correlators are in close agreement with simulations of a harmonically trapped prolate system. The formalism lays the foundation to analytically explore experimentally accessible hot Bose-Einstein condensates.

tractable approach for numerical simulations of hot matter-wave dynamics that includes all significant reservoir interaction processes. The SPGPE describes the evolution of a high-temperature partially condensed system within a classical field approximation, is valid on either side of the critical point, and has been used to quantitatively model the BEC phase transition [3], and BEC dynamics observed in high-temperature experiments [16]. The complete SPGPE includes a number-damping reservoir interaction (of Ginzburg-Landau type) described in previous works [3,16], and an additional interaction involving exchange of energy with the reservoir [2]; the latter is a number-conserving interaction that can in principle have a significant influence on dissipative evolution [8,12,17,18], yet due to technical challenges, its physical effects have thus far been largely unexplored.
Ehrenfest's theorem [19] relates, for example, the time derivative of the expectation (with respect to the wavefunction) of momentum and position of a particle with the potential experienced by the particle. The Ehrenfest relations can be extended to the Gross-Pitaevskii fluid -essentially without modification -due to the net cancellation of two-body interaction forces. Previously, development of the number-damping SPGPE theory was aided by Ehrenfest relations [20]. Such relations for ensemble averages provide an essential validity check for numerical work, and physical consistency tests for the reservoir theory, including ensemble averages expected in equilibrium.
In this work we derive stochastic collective equations for the SPGPE including the effects of both number-damping and energy-damping. These equations form the extension of Ehrenfest relations to finite-temperature stochastic field theory of Bose-Einstein condensates. They extend beyond the scope of previously derived relations of this type [20] by explictly retaining all noises and cutoff terms. Importantly, for many one-body operators the multiplicative noise in the energy-damped SPGPE is transformed into additive noise in the collective equations. This simplification opens the way for analytical treaments of a broad class of low-energy excitations in BEC including solitons, vortices, and collective modes, for both number-and energy-damping decay channels. The SPGPE theory of spinor and multicomponent systems [13] enables these techniques to be applied to systems where dissipation can only proceed via energy damping.
The paper is structured as follows. The GPE, PGPE, and SPGPE are introduced in Sec. 2, and the Ehrenfest relations for the GPE and PGPE briefly reviewed. In Sec. 3 we apply Ito's formula for change of variables to derive stochastic Ehrenfest relations for the SPGPE. As a first application,m in Sec. 4 we apply the formalism to the center of mass motion of a harmonically trapped quasi-1D system. We find that neglecting corrections arising due to the projector results in an analytically tractable equation of motion for the system position and momentum, taking the form of a vector Ornstein-Uhlenbeck equation. By comparing analytic solutions of this to full simulations of the SPGPE, we find validation of both the relevant stochastic Ehrenfest relation and the neglect of the projector corrections. In Sec. 5 we discuss links between our work and relevant linterature, and give our conclusions.

Gross-Pitaevskii equation
The Gross-Pitaevskii equation is the equation of motion for a scalar complex field evolving according to the Gross-Pitaevskii Hamiltonian where V(r, t) is an external potential and the interaction strength = 4π 2 a s /m is the two-body interaction strength in the cold-collision regime [21] via the s-wave scattering length a s and atomic mass m. The Gross-Pitaevskii equation is obtained by taking the functional derivative of the Gross-Pitaevskii Hamiltonian where the Gross-Pitaevskii operator is The total particle number in a system governed by the Gross-Pitaevskii equation is Observables A(t) of the system that are given by the expectation value of an operatorÂ are defined by The internal cancellation of s-wave interaction forces renders the Ehrenfest relations for the GPE identical to those of the Schrodinger equation [22]: for the position R(t), momentum P(t), angular momentum L(t), energy H(t), and number N(t) respectively.

Projection Operators
The SPGPE requires a clear formulation of a projection operator that formally and numerically projects the nonlinear GPE dynamics into a low-energy subspace. To define the projector, the external potential is split into a time-invariant part and a time-dependent part V(r, t) ≡ V(r) + δV(r, t). The time independent part defines the single-particle Hamiltonian The basis of representation is chosen to be the eigenstates of this Hamiltonian, satisfying H sp φ n (r) = n φ n (r), where n denotes the set of quantum numbers required to completely describe the basis, with corresponding energy eigenvalues n .
The projector can be written in operator form aŝ where the coherent region C is defined as the set of basis states beneath the cutoff: C ≡ {n : n ≤ c }.
In position space, this becomes an integral operator projecting an arbitrary function F(r) into C, with action on the function F(r) written as where we have made use of the projected Dirac-delta distribution The orthogonal projectorQ = 1 −P, has the action on the function F(r) where the incoherent region I is defined as I ≡ {n : n > c }. SinceP is Hermitian, for any two states |F , |G , we have F|P|G = F|PG = P F|G . In position space, this gives a useful property under integration where the complex conjugate projector is Analagous relations hold for the orthogonal projectorQ.

Projected Gross-Pitaevskii equation
The classical field is represented as a sum over the basis states φ n (r) with weight α n (t) where the field is now restricted to exist entirely within the coherent region. The projected Gross-Pitaevskii equation (PGPE) is obtained by taking the projected functional derivative of H, where these are defined byδ The projected functional derivative of a functional is related to the regular functional derivative bȳ describing Hamiltonian evolution for a field that is formally restricted to the C region. Defining the caonical momentum Π(r, t) ≡ i ψ(r, t) * , and projected functional Poisson bracket of any two functionals F[ψ, Π], G[ψ, Π], as the PGPE follows the usual Hamiltonian structure Provided δV ≡ 0, N and H are both conserved by the PGPE: ∂ t N = {N, H} = 0, ∂ t H = {H, H} = 0. As usual, additional conserved quantities may exists when the confining potential possesses additional symmetries. We emphasize that the conservation laws are a strict consequence of the projector in Eq. (21), and thus any numerical implementation must impose the projector to high precision in a suitable basis to accurately solve the PGPE.

Number-Damped Projected Gross-Pitaevskii equation
The PGPE gives evolution of the coherent region field in isolation, that is, there is no interaction between the coherent region and incoherent region. In SPGPE c-field theory, the incoherent region is usually assumed to be thermalised and is thus treated as a reservoir that acts as a damping mechanism for the coherent region. In the number-damped PGPE the noise in SPGPE is neglected, as are the energy damping terms, and the damping mechanism is parametrised entirely by the number-damping rate γ describing two-body collisions that transfer particles between C and I. Formulated in the convenient rotating frame defined by the reservoir chemical potential µ, the damped PGPE takes the form Lacking any thermal noise, the damped PGPE will evolve an arbitrary (non-trivial) initial condition to a ground state with chemical potential µ. For a thermal Bose reservoir the dimensionless parameter γ is [10] where It is not a priori evident that γ should assume such a simple form. However, an advantage of choosing a high energy cutoff c ∼ 3µ, as required within the SPGPE formalism, is that γ is independent of position to a very good approximation [10]. Typical experimentally relevant values are of order 10 −5 − 10 −4 .

Stochastic projected Gross-Pitaevskii equation
The derivation of the complete SPGPE for scalar BEC appeared in [9]. A later generalization allows for the possibility of multiple components and spins [13]. The derivation involves finding a master equation for the coherent region density operator. A series of approximations valid at high temperartures [20] then gives the high temperature regime master equation, which is mapped to an equivalent equation of motion for the multimode Wigner distribution function W[ψ, ψ * ] using quantum to classical operator correspondences [6,23]. The result of this is a generalized Fokker-Planck equation of motion for the Wigner function that includes third-order functional derivates. An equation of motion for a quasi-probability distribution can only be mapped to an SDE if it takes the form of a Fokker-Planck equation, that is, only contains functional derivatives up to second order and has a positive semi-definite diffusion matrix. Futher progress can be made by neglecting the third-order derivatives, an approximation known as the truncated Wigner approximation. Such an approximation is physically well justified at high temperatures where thermal and classical noise dominate over quantum noise [24]. Applying the truncated Wigner approximation leads to the Fokker-Planck equation where γ is defined in Eq. (25), and we have introduced the energy-damping potential with the coherent region current and the epsilon function and S (k) ≡ |k| −1 is the scattering kernel. The rate constant M is given by [9,12] The epsilon function has the useful property under integration that we will use below. Having arrived at a genuine Fokker-Planck equation involving only derivatives up to second order, we are now able to map the Wigner function evolution equation to an equivalent SDE for ψ(r) [25], provided the diffusion matrix is positive semi-definite. For the FPE Eq. (26) this condition is always satisfied, and the stochastic projected Gross-Pitaevskii equation (SPGPE) is [9, 10, 13] where the complex number-damping noise has the non-zero correlations and the real energy-damping noise has non-zero correlations The notation (S) indicates that the SDE is to be interpreted as a stochastic integral in Stratonovich form [25], and thus at a given time t the noises are not independent of the fields.

Ito form of the SPGPE
We can most easily change variables in a stochastic differential equation using results from Ito calculus. Recasting the SPGPE in Ito form involves reordering the projected function derivatives in the final term of the FPE and mapping the equation to a new SDE.
The FPE is equivalent to which maps to the Ito SPGPE where the first term in the second line is called the Stratonovich correction. In contrast to the Stratonovich SPGPE, in the Ito SPGPE [denoted (I)] the noises are always independent of the fields at a given time t. This formulation has distinct advantages for formal manipulation that we exploit below.

Functional Change of Variables
Consider any functional of the projected fields A ≡ A[ψ, ψ * , t]. Using the rules of Ito calculus, we can find the an SDE for A in the form where we consistently include all terms up to order dt, including terms quadratic in the noises that generate second derivatives. Usinḡ we rewrite (37) as a stochastic Ehrenfest equation with additional terms due to the projector where the dW i (t) are standard real Wiener processes with unit correlations and the Hamiltonian, number-damping, and energy-damping drift projector terms are given in Appendix A.1. The number-damping and energy-damping each have a corresponding drift, diffusion, and trace 1 term. It is worth remarking that that the general expression obtained above allows for any observable. Furthermore, it is well-suited for significant simplification via an ansatz for the wavefunction, provided the integrals may be evaluated. In the remainder of this paper we provide explicit stochastic equations for particular observables of interest, and also carry out the ansatz approach for describing Kohn-mode oscillations of a highly oblate system in a parabolic trap.

One-body operators
If we restrict our attention to one-body operators, the projected functional calculus assumes a particularly simple form. Consider an observable that is the expectation of a one-body operator where the operatorÂ is independent of ψ, ψ * . The non-zero projected functional derivatives arē The projected functional derivatives corresponding to the operatorÂ are equivalent to the regular functional derivatives of the totally projected operatorÂ P ≡PÂP; the change of variables is then simple to construct. Again writing in the form of a wave function average with additional projector terms, the SDE of the observable A is where the projector terms are given in Appendix A.

Finite-temperature stochastic Ehrenfest relations
We consider the Ehrenfest relations for position, momentum, angular momentum, grand canonical energy, and coherent region particle number. The complete set of stochastic Ehrenfest relations for the SPGPE are This set of equations are our main result. They take the form of generalized Ehrenfest relations with additional damping and noise terms arising from the reservoir coupling processes. They provide a starting point for finding analytic descriptions of hot BEC dynamics, and also provide tests for numerical consistency of SPGPE simulations. We make the following remarks: i) Neglecting energy-damping and taking the average over the noise (in Ito form the fields and noises are uncorrelated), we immediately recover the Ehrenfest relations for the number-damped SPGPE as found in [20]. Those Ehrenfest relations described the evolution of ensemble avarages, whereas in the present formulation we retain all noises in the collective equations.
ii) Superficially, the multiplicative noise in the SPGPE Eq. (32) may appear to have been transformed into additive noise. However, we emphasize that multiplicative noises remain present in terms of the form Â dW j , since Â is not a noise average.
iii) Reduction to additive noise can be achieved in special cases where the system can be welldescribed by a suitable physically motivated ansatz wavefunction. If such an ansatz is available, reduction to an additive noise SDE provides a significant simplification that can enable analytical progress [18].
iv) The projector terms are all consistently accounted for, and in general contribute additional noises. However, provided the basis of projection is properly chosen, their effect is typically only a small correction. Testing whether such terms are negligible provides a useful consistency test for a well-chosen cutoff.

Approach to equilibrium
The first point deserves further discussion. Collective equations of the Ehrenfest type may be obtained for the number-damped SPGPE, including the projector term, by both neglecting all energy damping terms, and averaging over the remaining number-damping noises [20]. For example, the evolution of the grand canonical Hamiltonian K reduces to where denotes both wavefunction expectation values, and averages over the noise, and Tr is the traced noise average. We can use this to gain an intuitive picture of the grand-canonical ensemble furnished by the number-damped SPGPE. For a time-independent potential, the grand-canonical Hamiltonian K ≡ H − µN has two notable regimes. First, in the absence of noise, obtained by formally setting T = 0 in Eq. (46), the evolution will obey dK/dt ≤ 0, eventually taking the field to the ground state ψ 0 satisfying Lψ 0 ≡ µψ 0 . Second, when T > 0, the steady state satisfies a consequence of the fluctuation-dissipation relation. The fact that Eq. (47) depends only on T in equilibrium is a fundamental property of the ensemble, and also serves as a very useful numerical validation of any simulations: all equilibrium properties must be independent of γ.

Kohn mode
As a first application of the formalism we consider Kohn oscillations in a quasi-1D harmonically trapped system with frequency for transverse trapping ω ⊥ such that ω ⊥ ω, the frequency for the axial trap. Kohn's theorem states that in a harmonically trapped system the centre of mass undergoes bulk oscillations about the trap centre at the trapping frequency. Since the SPGPE reservoir theory describes the I region as time independent, the system it describes cannot satisfy the Kohn theorem. In any case, we refer to center of mass motion as the Kohn mode. We emphasize that the timeindependent reservoir approximation is compatible with systems of physical iterest. In particular, for a scalar BEC, the harmonic trapping potential may have post-harmonic corrections at high energy. Furthermore, a scenario that occurs in sympathetic cooling involves one component harmonically trapped, and another thermal component held in a different trap so that the system does not satisfy the conditions of the Kohn theorem. It should also be clear that the Thomas-Fermi ansatz cannot describe any of the non-condensate contained in the c-field. However, it contains the essential feature of the Kohn mode in the absence of dissipation, namely that the center of mass moves without changing the spatial density profile. Our aim is to integrate out the spatial degrees of freedom and find effective stochastic equations that include the effects of damping and thermal noise on the Kohn mode.
A simple low-dimensional theory can be found in the regime where the reservoir remains three dimensional 2 , and the low dimensional subspace is well-described by projecting onto the transverse ground state of the confining potential. The one dimensional SPGPE takes an identical functional form to the three dimensional SPGPE, but with dimensionally reduced damping and noise coefficients [17]. The number-damping term is only changed by the replacement → 1D , as are the Hamiltonian terms. The scattering kernel in the 1D reduction of Eq. (29) becomes where erfcx(q) ≡ e q 2 erfc(q) is the scaled complementary error function, and a ⊥ ≡ √ /mω ⊥ is the transverse harmonic oscillator length, much smaller than a ω = √ /mω. The system is assumed to be sufficiently condensed that the center of mass motion can be approximated using a Thomas-Fermi wavefunction ansatz. The Thomas-Fermi wave function allowing for arbitrary variations in the centre of mass position x(t) and momentum p(t) is Using this ansatz, the integrals for Eqs. (45) can be evaluated analytically within the Thomas-Fermi approximation. 3

Stochastic equations
Integrating over the Thomas-Fermi ansatz, the stochastic Ehrenfest relations for position and momentum take the form of a pair of coupled stochastic differential equations where we have derived the number-damping and energy-damping drift and diffusion constants from the effective one dimensional SPGPE [17] and the projector terms are given in Appendix B.1. The projector terms Eqs. (76), (77) severely constrain exact analytic progress beyond this point. However, provided these terms may be safely neglected, we arrive at an approximate set of equations that can be solved analytically. For consistent c-field simulations, neglecting the projector terms in the stochastic Ehrenfest relations will always be a good approximation, since for a well-chosen cutoff the mode population near the cutoff will be relatively small. Since we are considering equilibrium states, we can assume that the values of x(t) and p(t) will also be small relative to the characteristic harmonic oscillator length and momentum scales; we test our assumptions by solving the SPGPE numerically. For this purpose we can rewrite the equations of motion as a single equation of motion for the dimensionless complex variable and neglect terms of higher order than linear in z(t), z * (t), to find with projector terms given explicitly in Appendix B.2. The first term describes simple harmonic motion, and the second and third terms damp the Kohn mode amplitude. The third term vanishes when energy and number damping occur at the same rate; typically this is not the case, and the term distorts the spiral decay of the phase-space amplitude described by the second term, as shown in Figure  2.

Analytic and Numerical solutions
In this section we present an anlytical treatment of the stochastic equations derived in Section 4.1, and compare the results with numerical simulations of the 1D SPGPE.
Neglecting the projector terms, we can write the coupled differential equation as a vector SDE representing an Ornstein-Uhlenbeck process 4 where u(t) = [x(t), p(t)] , and are the drift and diffusion matrices respectively, and dW(t) = [dW 1 (t), dW 2 (t)] is a vector of independent real Wiener processes with correlations dW n (t)dW m (t) = δ mn dt. The SDE has the formal solution where have assumed the initial state u(0) to be deterministic. The mean undergoes exponential decay u(t) = exp [−Λt] u(0). For a system with x(0) = x 0 and p(0) = 0, the centre of mass position over time is given by where we have defined the frequency ω γε = ω 2 − (Λ γ − Λ ε ) 2 . We now test our analytical description against SPGPE simulations. While the SPGPE can be used for quantitative modelling [16,26,27] here our aim is to test our analytic solutions of the SPGPE using the stochastic Ehrenfest relations. We choose parameters that are representative of BEC experiments, and leave a first principles treatment of reservoir interactions [17] for later work.
We perform simulations of the 1D SPGPE [17] with the initial condition given by Eq. (49) with x(0) = p(0) = 0. We use a chemical potential of µ = 100 ω, an energy cutoff of c = 250 ω, a temperature of T = 500 ω/k B , an interaction strength of 1D = 0.01 ωa ω , an energy-damping rate of M = 0.0005a 2 ω , and a number-damping rate of γ = 0.001. We chose values for the damping rates such that Λ ε ≈ Λ γ and thus neither damping process is dominant over the other. Timescales are considered in units of the harmonic oscillator time period t ω ≡ 2π/ω. We perform simulations for an ensemble of 500 trajectories with an initial displacement of x(0) = 0.2R to test the exponential decay rate. The result of this is shown in Fig. 3, where we see that the analytically predicted decay rate agrees well with simulations of the complete SPGPE.
To explore this regime of energy-damped dynamics futher further, we also consider the steadystate correlations given by Fourier transforming with respect to τ in the steady-state, the fluctuation spectra are given by The steady-state correlations for position-position, momentum-momentum, and momentum-position are respectively. The steady-state spectra for position-position, momentum-momentum, and momentum- position are respectively. We compare these analytic solutions with the numerical data from equilibrated SPGPE simulations. In Fig. 4 we show the steady-state correlation functions for an ensemble of 5000 trajectories 5 . We see that the analytic and numeric results show excellent agreement for short times with differences becoming more pronounced for larger τ. Similarly, the steady-state spectra for position-position, momentum-momentum, and momentum-position are shown in Fig. 5, where the numeric spectra are obtained using the Wiener-Khinchin theorem applied to the numeric steady-state correlations. The bulk oscillation seen in Fig. 4 is the Kohn mode, as seen from the peak at Ω = ω in Fig. 5. The spectral linewidth represents the decay time of the two-time correlation function, determined by Λ ε and Λ γ , the energy and number damping rates. Again we see that the analytic and numeric results show good agreement. To assess the validity of neglecting projector terms in our analytical treatment, we evaluate their contribution numerically in Appendix B.3, finding that indeed the projector correction is negligible.

Discussion
The stochastic Ehrenfest relations derived here using projected functional calculus contain many projector terms, posing a technical challenge. As an application of projected functional calculus, our approach bears comparison with that of Opanchuk et al [28]. In that work a number of useful functional relations were derived for Wigner phase space methods [5] involving projected fields. The action of the projector was reduced to the action of matrix operations, as is always possible in a finite basis. While not without its own technical challenges, our approach has the advantage that the equations of motion have an obvious link with the continuum limit recovered for a high cutoff.
An alternate stochastic reservoir theory, namely the Stoof SGPE [29], has been used in several studies of finite-temperature BEC dynamics [30][31][32][33][34]. While the SGPE lacks a projection operator of the form used in the present work, our aims are similar in spirit to the work of Duine et al. [35] applying path-integral techniques to the Stoof SGPE to find effective stochastic equations for a reduced set of variables. However, the explicit high-energy cutoff in the SPGPE necessitates a different technical approach, with significant differences appearing in the resulting stochastic equations. As the energydamping mechanism is absent from the Stoof SGPE theory, an interesting future direction is to identify further experimentally accessible regimes capable of distinguishing between the two approaches [27]. A related question recieving recent attention is that of ensemble equivalence [36]. In the context of the present work, the relative importance of number and energy damping in a given system will influence the ensemble generated by the SPGPE.

Conclusions
We have developed a set of exact stochastic Ehrenfest relations for the complete stochastic projected Gross-Pitaevskii equation [9], an equation of motion which has significant application in the study of finite-temperature Bose-Einstein condensates [3,8,10]. In addition to the number-damping process, the SPGPE contains a number-conserving dissipative mechanism that only involves energy transfer between the system and reservoir. Our main result, the stochastic Ehrenfest relations, retain the stochastic nature of the SPGPE and explicitly contain terms that result from both dissipative processes.
In applications of the SPGPE, the energy cutoff must be chosen such that the bulk of the coherent region modes are significantly occupied compared to the modes in the incoherent region and at the cutoff energy. The relative occupation at the cutoff is what determines the presence of spurious dynamics due to the projector, and this population is chosen to be of order unity. Ultimately, the cutoff should be chosen such that the projector corrections are small. We verified this in two ways. Considering the motion of a finite-temperature quasi-1D condensate near equilibrium, we tracked the size of the largest projector corrections and saw they are indeed small. We also compared the steady-state correlations of position and momentum to analytic solutions derived by neglecting the projector corrections, finding excellent agreement.
We have shown that stochastic Ehrenfest relations can be used to obtain analytic equations that agree with numerical solutions of the full SGPE. Future work will use this formalism to analytically explore experimentally accessible systems. The stochastic Ehrenfest relations provide a useful starting point for describing a range of dissipative dynamics in hot BEC including soliton evolution [37], phase-slip dynamics [38], sympathetic cooling [39,40], spinor BECs [41], and quantum turbulence [42][43][44][45].

A Projector terms A.1 General operators
For general operatorÂ, the drift projector terms are the number-damping, and energy-damping diffusion projector terms are and the epsilon term is

B Kohn mode projector terms
We present the projector terms for the SPGPE describing a BEC in an oblate parabolic trap with only one effective C-region dimension [17]. We first consider the position and momentum terms separately, and then give the terms for the dimensionless complex variable Eq. (53).

B.1 Position and momentum
Noting that the projector and epsilon terms for x and p are and q H p = 1 N 2mω(n c + 1) Re α * n c (t) dxφ * n c +1 (x)ψ(x)n(x) , (77a) q γ p = 1 γ N 2mω(n c + 1) Im α * n c (t) dxφ * n c +1 (x)ψ(x)n(x) , (77b) epsilon term would in general be non-negligible, and in fact diverges in the infinite cutoff limit. We thus monitor the magnitudes of q H z and dz ε over the course of an ensemble of numerical trajectories. We define the relative cutoff term magnitudes by noting that |ż(t)| is strictly non-zero for harmonic motion. If these values remain significantly less than unity, then we may conclude that the effects of the cutoff terms are negligible. For the simulations of section 4.2, example relative cutoff term magnitudes are shown in Fig. 6, where the initial state is Eq. (49) with x(0) = p(0) = 0 and we have taken an ensemble average over 1000 trajectories. While the relative magnitudes can reach as high as ∼ 0.1 early in the dynamics, we see that once the system has equilibrated they approach a steady-state of order ∼ 0.01.