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

TL;DR

This paper develops exact stochastic Ehrenfest relations for hot Bose-Einstein condensates, accounting for damping mechanisms, and validates the theory through simulations of center of mass fluctuations in a trapped system.

Contribution

It introduces a novel analytical framework for describing number and energy damping in hot Bose-Einstein condensates using stochastic Ehrenfest relations.

Findings

Close agreement between simulations and analytical results

Foundation for exploring hot Bose-Einstein condensates analytically

Enhanced understanding of damping mechanisms in BECs

Abstract

Describing partially-condensed 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. We test the theory by applying it to the centre of mass fluctuations of a harmonically trapped prolate system, finding close agreement between c-field simulations and analytical results. The formalism lays the foundation to analytically explore experimentally accessible hot Bose-Einstein condensates.

Equations

$$H=\int d^3\mathbf{r}{\psi }^{\ast }(\mathbf{r},t)\left(-\frac{{\hslash }^2{\nabla }^2}{2m}+V(\mathbf{r},t)+\frac{\text{varg}}{2}|\psi (\mathbf{r},t){|}^2\right)\psi (\mathbf{r},t),$$

$$i\hslash \frac{\partial \psi (\mathbf{r},t)}{\partial t}=\frac{\delta H}{\delta {\psi }^{\ast }(\mathbf{r},t)}=L\psi (\mathbf{r},t),$$

$$L\psi (\mathbf{r},t)\equiv \left(-\frac{{\hslash }^2{\nabla }^2}{2m}+V(\mathbf{r},t)+\text{varg}|\psi (\mathbf{r},t){|}^2\right)\psi (\mathbf{r},t).$$

$$N=\int d^3\mathbf{r}{\psi }^{\ast }(\mathbf{r},t)\psi (\mathbf{r},t),$$

$$A(t)\equiv ⟨\hat{A}⟩=\int d^3\mathbf{r}⟨\psi |\hat{A}|\mathbf{r}⟩\psi (\mathbf{r},t)=\int d^3\mathbf{r}{\psi }^{\ast }(\mathbf{r},t)⟨\mathbf{r}|\hat{A}|\psi ⟩,$$

$$\frac{d\mathbf{R}(t)}{dt}$$

$$\frac{d\mathbf{P}(t)}{dt}$$

$$\frac{d\mathbf{L}(t)}{dt}$$

$$\frac{dH(t)}{dt}$$

$$\frac{dN(t)}{dt}$$

$$H_{\mathrm{sp}}\equiv -\frac{{\hslash }^2{\nabla }^2}{2m}+V(\mathbf{r}).$$

$$\hat{\mathcal{P}}$$

$$\mathcal{P}\left\{F(\mathbf{r})\right\}\equiv \sum _{n\in C}{\varphi }_n(\mathbf{r})\int d^3{\mathbf{r}}^{\prime }{\varphi }_n^{\ast }({\mathbf{r}}^{\prime })F({\mathbf{r}}^{\prime })=\int d^3{\mathbf{r}}^{\prime }\delta (\mathbf{r},{\mathbf{r}}^{\prime })F({\mathbf{r}}^{\prime }),$$

$$\delta (\mathbf{r},{\mathbf{r}}^{\prime })\equiv \sum _{n\in C}{\varphi }_n(\mathbf{r}){\varphi }_n^{\ast }({\mathbf{r}}^{\prime })=⟨\mathbf{r}|\hat{\mathcal{P}}|{\mathbf{r}}^{\prime }⟩.$$

$$\mathcal{Q}\left\{F(\mathbf{r})\right\}$$

$$\psi (\mathbf{r},t)=\sum _{n\in C}{\alpha }_n(t){\varphi }_n(\mathbf{r}),$$

$$\frac{\bar{\delta }}{\bar{\delta }\psi (\mathbf{r},t)}=\sum _{n\in C}{\varphi }_n^{\ast }(\mathbf{r})\frac{\partial }{\partial {\alpha }_n},\frac{\bar{\delta }}{\bar{\delta }{\psi }^{\ast }(\mathbf{r},t)}=\sum _{n\in C}{\varphi }_n(\mathbf{r})\frac{\partial }{\partial {\alpha }_n^{\ast }}.$$

$$\frac{\bar{\delta }F[\psi ,{\psi }^{\ast }]}{\bar{\delta }\psi (\mathbf{r},t)}={\mathcal{P}}^{\ast }\left\{\frac{\delta F[\psi ,{\psi }^{\ast }]}{\delta \psi (\mathbf{r},t)}\right\},\frac{\bar{\delta }F[\psi ,{\psi }^{\ast }]}{\bar{\delta }{\psi }^{\ast }(\mathbf{r},t)}=\mathcal{P}\left\{\frac{\delta F[\psi ,{\psi }^{\ast }]}{\delta {\psi }^{\ast }(\mathbf{r},t)}\right\}.$$

$$i\hslash \frac{\partial \psi (\mathbf{r},t)}{\partial t}=\frac{\bar{\delta }H}{\bar{\delta }{\psi }^{\ast }(\mathbf{r},t)}=\mathcal{P}\left\{L\psi (\mathbf{r},t)\right\},$$

$$\{F,G\}$$

$${\partial }_t\psi (\mathbf{r})$$

$$\{N,H\}$$

$$\frac{\partial W[\psi ,{\psi }^{\ast }]}{\partial t}=$$

$$+\int d^3\mathbf{r}\left[\frac{\gamma k_{\mathrm{B}}T}{\hslash }\frac{{\bar{\delta }}^{(2)}}{\bar{\delta }\psi (\mathbf{r})\bar{\delta }{\psi }^{\ast }(\mathbf{r})}+\text{h.c.}\right]W[\psi ,{\psi }^{\ast }]$$

$$+\int d^3\mathbf{r}\int d^3{\mathbf{r}}^{\prime }\frac{k_{\mathrm{B}}T}{\hslash }\epsilon (\mathbf{r}-{\mathbf{r}}^{\prime })\left[\frac{\bar{\delta }}{\bar{\delta }\psi (\mathbf{r})}\psi (\mathbf{r})\frac{\bar{\delta }}{\bar{\delta }{\psi }^{\ast }({\mathbf{r}}^{\prime })}{\psi }^{\ast }({\mathbf{r}}^{\prime })+\text{h.c.}\right]W[\psi ,{\psi }^{\ast }]$$

$$-\int d^3\mathbf{r}\int d^3{\mathbf{r}}^{\prime }\frac{k_{\mathrm{B}}T}{\hslash }\epsilon (\mathbf{r}-{\mathbf{r}}^{\prime })\left[\frac{\bar{\delta }}{\bar{\delta }\psi (\mathbf{r})}\psi (\mathbf{r})\frac{\bar{\delta }}{\bar{\delta }\psi ({\mathbf{r}}^{\prime })}\psi ({\mathbf{r}}^{\prime })+\text{h.c.}\right]W[\psi ,{\psi }^{\ast }].$$

$$\gamma$$

$$V_{\epsilon }(\mathbf{r})$$

$$\epsilon (\mathbf{r})\equiv \frac{\mathcal{M}}{(2\pi {)}^3}\int d^3\mathbf{k}S(k)e^{i\mathbf{k}\cdot \mathbf{r}},$$

$$\mathcal{M}$$

$$\text{SS}d\psi (\mathbf{r})$$

$$⟨dW(\mathbf{r},t)d{W}^{\ast }({\mathbf{r}}^{\prime },t)⟩=\frac{2\gamma k_{\mathrm{B}}T}{\hslash }\delta (\mathbf{r},{\mathbf{r}}^{\prime })dt,$$

$$⟨dU(\mathbf{r},t)dU({\mathbf{r}}^{\prime },t)⟩=\frac{2k_{\mathrm{B}}T}{\hslash }\epsilon (\mathbf{r}-{\mathbf{r}}^{\prime })dt.$$

$$\frac{\partial W[\psi ,{\psi }^{\ast }]}{\partial t}=$$

Full text

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

Rob G. McDonald1,2, Peter Barnett1,2, Fradom Atayee1,2, Ashton S. Bradley1,2*,

1 Department of Physics, University of Otago, Dunedin, New Zealand

2 The Dodd-Walls Centre for Photonic and Quantum Technologies

• ashton.bradley@otago.ac.nz

Abstract

**Describing partially-condensed 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. We test the theory by applying it to the centre of mass fluctuations of a harmonically trapped prolate system, finding close agreement between c-field simulations and analytical results. The formalism lays the foundation to analytically explore experimentally accessible hot Bose-Einstein condensates. **

Contents

1. 1 Introduction

2. 2 Background

3. 2.1 Gross-Pitaevskii equation

4. 2.2 Stochastic projected Gross-Pitaevskii equation

5. 2.2.1 Projector

6. 2.2.2 Equations of motion

7. 2.2.3 Energy cutoff and the single-particle basis

8. 3 Stochastic Ehrenfest relations

9. 3.1 Ito form of the SPGPE

10. 3.2 Functional Change of Variables

11. 3.3 One-body operators

12. 3.4 Noise correlations for multiple moments

13. 3.5 Finite-temperature stochastic Ehrenfest relations

14. 4 Application to centre of mass fluctuations

15. 4.1 Stochastic equations for centre of mass

16. 4.2 Analytic and numerical solutions

17. 5 Discussion and Conclusions

18. 5.1 Discussion

19. 5.2 Conclusions

20. A Projector terms

21. A.1 General operators

22. A.2 One-body operators

23. B Centre of mass equations

24. B.1 Position and momentum

25. B.2 Dimensionless variable $z(t)$

26. B.3 Numerical evaluation of projector terms

1 Introduction

A system of Bose atoms with temperature $T$ undergoes a dramatic change in behavior at the critical temperature for the formation of a Bose-Einstein condensate (BEC), $T_{c}$. Far below the BEC transition $T\ll T_{c}$, a nearly pure BEC forms, consisting of a highly occupied many-body quantum state; in this regime dilute gas BEC are renowned both for their high experimental control, and precise theoretical description [1]. At temperatures $T\gg T_{c}$ thermal energy dominates, the quantum statistics are unimportant, and a Boltzman description captures the physical properties of the atoms. When $T\sim T_{c}$ the quantum statistics of the atoms are decisive, despite appreciable thermal energy. In a hot BEC $T\lesssim T_{c}$, competition between thermal and interaction effects leads to fragmentation of the condensate, and formation of vortices, solitons, and phononic excitations. A cooling quench across the transition can inject interesting excitations into the BEC that form as remnants of the broken $U(1)$ symmetry [2, 3], and rich turbulent dynamics develop from the competition between thermal, quantum, and interaction effects, posing a challenge for theory.

At low temperatures that contain significant condensate and non-condensate fraction, mean field theory provides a useful description, upon which the GPE and its generalizations are based. The Zaremba-Nikuni-Griffin $U(1)$ symmetry breaking approach has proven well suited for practical calculations in this low-temperature regime [4], having the virtue that the interactions between condensate and thermal cloud, and their respective dynamics, are all included in the dynamical description. However, it’s strength for low temperatures presents a limitation at high temperatures: despite notable successes, e.g. for collective modes [5], the symmetry-breaking ansatz has limited scope for describing strongly fluctuating systems containing large non-condensate fraction. Fortunately, the scope of the Gross-Pitaevskii equation (GPE) upon which ZNG is based goes much further than mean-field theory: GPE-like field equations appear naturally in phase-space representations of Bose gases [6], suggesting a generalization of quantum optical open systems theory [7] for describing hot Bose gases. Indeed, various generalized Gross-Pitaevskii theories have been developed for the high temperature regime, describing many modes that are weakly occupied, interacting, and partially coherent [8, 9]. Devoid of the symmetry-breaking ansatz, these approaches have advantages for high-temperature work.

One approach capable of describing experiments across the phase transition is known as the stochastic projected Gross-Pitaevskii equation (SPGPE) [10, 11, 12, 13, 14]. The SPGPE was developed as a synthesis of quantum kinetic theory [15] and the projected Gross-Pitaevskii equation [16], and provides a 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 [17]. The complete SPGPE includes a number-damping reservoir interaction (of Ginzburg-Landau type) described in previous works [3, 17], 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 [9, 13, 18, 19], yet due to technical challenges, its physical effects have thus far been largely unexplored.

Ehrenfest’s theorem [20] 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 [21]. 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 stochastic Ehrenfest relations (SERs) 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 [21] by explicitly 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 treatments of a broad class of low-energy excitations in BEC including solitons, vortices, and collective modes, for both number- and energy-damping decay channels. As a test we apply the Ehrenfest relations to the centre of mass fluctuations of a harmonically trapped system tightly confined along two spatial dimensions. We find that the analtyic solution of the SER for the centre of mass is in close agreement with SPGPE simulations. The SPGPE theory of spinor and multicomponent systems [14] 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 SERs for the SPGPE. We then derive a set of simple fluctuation-dissipation relations describing equilibrium ensembles. In Sec. 4 we test the formalism on the centre of mass motion of a harmonically trapped quasi-1D system. In Sec. 5 we discuss links between our work and relevant linterature, and give our conclusions.

2 Background

2.1 Gross-Pitaevskii equation

The Gross-Pitaevskii equation is the equation of motion for a scalar complex field evolving according to the Gross-Pitaevskii Hamiltonian

[TABLE]

where $V(\mathbf{r},t)$ is an external potential and the interaction strength $\varg=4\pi\hbar^{2}a_{s}/m$ is the two-body interaction strength in the cold-collision regime [22] via the $s$-wave scattering length $a_{s}$ and atomic mass $m$. The Gross-Pitaevskii equation may be generated by taking the functional derivative of the Gross-Pitaevskii Hamiltonian

[TABLE]

where the Gross-Pitaevskii operator is

[TABLE]

The total particle number

[TABLE]

is conserved under the Hamiltonian evolution. Observables $A(t)$ of the system that are given by the expectation value of an operator $\hat{A}$ are defined by

[TABLE]

for example $\mathbf{R}(t)=\int d^{3}\mathbf{r}\;\psi^{*}(\mathbf{r})\mathbf{r}\psi(\mathbf{r})$ is the mean position of the atoms. The internal cancellation of s-wave interaction forces renders the Ehrenfest relations for the GPE identical to those of the Schrodinger equation [23]:

[TABLE]

for the position $\mathbf{R}(t)$, momentum $\mathbf{P}(t)$, angular momentum $\mathbf{L}(t)$, energy $H(t)$, and number $N(t)$ respectively. The focus of the present article is to develop a detailed theory of modifications to these collective equations that occur due to reservoir interactions in a Bose-Einstein condensate.

As is often convenient when changing representations, here we are using $|\psi\rangle$ to represent the abstract state ket associated with the pure-state wavefunction $\psi(\mathbf{r})=\langle\mathbf{r}|\psi\rangle$. In the context of phase-space methods introduced below, it should be noted that we do not use this notation for the many-body quantum state, but just as a convenient shorthand for the wavefunction appearing in individual trajectories. Many-body expectation values are calculated as ensemble averages over trajectories of the stochastic field theory [9]. When generalizing the above definitions [(6)—(10)] to the many-body theory as described below, the average $\langle\cdot\rangle$ refers only to the wavefunction expectation over position coordinates. The combination of wavefunction and ensemble average will be denoted by $\langle\langle\cdot\rangle\rangle$.

2.2 Stochastic projected Gross-Pitaevskii equation

The derivation of the complete SPGPE for scalar BEC appeared in [10]. A later generalization allows for the possibility of multiple components and spins [14]. The derivation involves finding a master equation for the coherent region density operator. A series of approximations valid at high temperartures [21] then gives the high temperature regime master equation, which is mapped to an equivalent equation of motion for the multimode Wigner distribution function $W[\psi,\psi^{*}]$ using quantum to classical operator correspondences [7, 24]. The result of this is a generalized Fokker-Planck equation (FPE) 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 [25]. To establish notation and formulate the problem, we now summarize the key results.

2.2.1 Projector

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(\mathbf{r},t)\equiv V(\mathbf{r})+\delta V(\mathbf{r},t)$. The time independent part defines the single-particle Hamiltonian

[TABLE]

The basis of representation is chosen to be the eigenstates of this Hamiltonian, satisfying $H_{{\scriptstyle{\rm sp}}}\phi_{n}(\mathbf{r})=\epsilon_{n}\phi_{n}(\mathbf{r})$, where $n$ denotes the set of quantum numbers required to completely describe the single-particle basis states, with corresponding energy eigenvalues $\epsilon_{n}$.

The projector can be written in operator form as

[TABLE]

where the coherent region $C$ is defined as the set of basis states beneath the cutoff: $C\equiv\{n:\epsilon_{n}\leq\epsilon_{\rm c}\}$. In position space, this becomes an integral operator projecting an arbitrary function $F(\mathbf{r})$ into $C$, with action on the function $F(\mathbf{r})$ written as

[TABLE]

where we have made use of the projected Dirac-delta distribution

[TABLE]

The orthogonal projector $\hat{{\cal Q}}=1-\hat{{\cal P}}$, has the action on the function $F(\mathbf{r})$

[TABLE]

where the incoherent region $I$ is defined as $I\equiv\{n:\epsilon_{n}>\epsilon_{\rm c}\}$. Since $\hat{\cal P}$ is Hermitian, for any two states $|F\rangle$, $|G\rangle$, we have $\langle F|\hat{{\cal P}}|G\rangle=\langle F|\hat{{\cal P}}G\rangle=\langle\hat{\cal P}F|G\rangle$.

The classical field is represented as a sum over the basis states $\phi_{n}(\mathbf{r})$ with weight $\alpha_{n}(t)$

[TABLE]

where the field is now restricted to exist entirely within the coherent region. Without reservoir interactions, the projected Gross-Pitaevskii equation (PGPE) [16] is obtained by taking the projected functional derivative of $H$, where these are defined by

[TABLE]

The projected functional derivative is related to the regular functional derivative by

[TABLE]

In particular, the equation of motion for the c-field Eq. (16) is given by the PGPE:

[TABLE]

describing Hamiltonian evolution for a field that is formally restricted to the $C$ region. Defining the caonical momentum $\Pi(\mathbf{r},t)\equiv i\hbar\psi(\mathbf{r},t)^{*}$, and projected functional Poisson bracket of any two functionals $F[\psi,\Pi]$, $G[\psi,\Pi]$, as

[TABLE]

the PGPE follows the expected Hamiltonian structure

[TABLE]

Provided $\delta V\equiv 0$, $N$ and $H$ are both formally conserved by the PGPE: $\partial_{t}N=\{N,H\}=0$, $\partial_{t}H=\{H,H\}=0$. A short account of the former is instructive. Evaluating the projected derivatives, we have

[TABLE]

Since ${\cal P}$ is hermitian, and ${\cal P}\psi\equiv\psi$ for projected fields, we then have $\{N,H\}\equiv 0$. As usual, additional conserved quantities may exist when the confining potential possesses additional symmetries, and numerical conservation requires that the basis also respect the symmetry. In practice this requires that the basis modes are eigenstates of the the symmetry generator.

2.2.2 Equations of motion

The PGPE gives Hamiltonian evolution of the coherent region field in isolation. 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. Applying the truncated Wigner approximation to the Born-Markov master equation for the Bose field leads to the Fokker-Planck equation [10]

[TABLE]

The number-damping interction is parametrised by the number-damping rate $\gamma$ describing two-body collisions that transfer particles between $C$ and $I$. For a thermal Bose reservoir the dimensionless rate $\gamma$ is [11]

[TABLE]

where $\lambda_{\rm dB}=\sqrt{2\pi\hbar^{2}/mk_{\rm B}T}$ is the thermal de Broglie wavelength, $\beta=1/k_{\rm B}T$, $\Phi[z,x,a]=\sum_{k=0}^{\infty}z^{k}/(a+k)^{x}$ is the Lerch transcendent.111A consequence of the high-energy cutoff ($\epsilon_{\rm c}\sim 3\mu$) used in the SPGPE is that $\gamma$ is approximately independent of position, simplifying the SPGPE. Typical experimentally relevant values are of order $10^{-5}-10^{-4}$ [11]. Energy damping is described by the potential

[TABLE]

where $\mathbf{j}(\mathbf{r})=i\hbar\left[\psi(\mathbf{r})\nabla\psi^{*}(\mathbf{r})-\psi^{*}(\mathbf{r})\nabla\psi(\mathbf{r})\right]/2m$ is the current, and the epsilon function

[TABLE]

where $S(k)\equiv|\mathbf{k}|^{-1}$ is the scattering kernel. For identical bosons the rate constant $\mathcal{M}$ is given by [10, 13]

[TABLE]

while for distinct reservoir species, it is reduced by $1/2$ [14]. The scattering potential Eq. (25) describes the reservoir interaction involving transfer of energy and momentum without transfer of particles.

Direct solution of the Fokker-Planck equation of motion (2.2.2) is impractical numerically. Working with this FPE involving only derivatives up to second order, the Wigner function evolution equation can be mapped to an equivalent SDE for $\psi(\mathbf{r})$ [26], provided the diffusion matrix is positive semi-definite. For the FPE Eq. (2.2.2) this condition is always satisfied, and the stochastic projected Gross-Pitaevskii equation (SPGPE) is [10, 11, 14]

[TABLE]

where the (complex) number-damping noise has the non-zero correlations

[TABLE]

and the (real) energy-damping noise has non-zero correlations

[TABLE]

The notation $\SS$ indicates that the SDE is to be interpreted as a stochastic integral in Stratonovich form [26], and thus at a given time $t$ the noises are not independent of the fields.

2.2.3 Energy cutoff and the single-particle basis

Before presenting our results, we remark that there exists a variety of interpretations in the literature of a central feature of $C$-field theory: what constitutes a cutoff, and how should it be imposed in practice? In some cases any convenient numerical basis of representation is truncated. However, truncation in an arbitrary basis will not usually generate a consistent energy cutoff, and some care should be taken when choosing a basis of representation. In particular, if a highly inappropriate basis of representation is used, then thermal noise driving the $C$-field modes will generate thermalization artifacts [21]. In our approach to the SPGPE we impose the cutoff in the single-particle basis generated by the time-independent confining potential. There is a clear physical motivation for this choice (clearly articulated by the pioneers of the PGPE): at sufficiently high energy — usually of order 2-3$\mu$ — the interacting many-body problem is approximately diagonal in the single-particle basis. Thus a well-chosen cutoff in the single-particle basis is also an approximate energy cutoff for the many-body system, precisely what is required for its separation into coherent and incoherent subspaces. These considerations will feature below in our discussion of equilibrium properties of the SPGPE.

3 Stochastic Ehrenfest relations

Ehrenfest relations for the stochastic field can be derived by casting the SPGPE in Ito form. The noises and fields are then uncorrelated, allowing change of variables using the rules of Ito calculus.

3.1 Ito form of the SPGPE

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

[TABLE]

which maps to the Ito SPGPE

[TABLE]

where the first term in the second line is called the Stratonovich correction. In contrast to the Stratonovich SPGPE, in the Ito SPGPE [denoted $(\textbf{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.

3.2 Functional Change of Variables

Consider any functional of the projected fields $A\equiv A[\psi,\psi^{*},t]$. Using the rules of Ito calculus, we can find the an SDE for $A$ in the form

[TABLE]

where we consistently include all terms up to order $dt$, including terms quadratic in the noises that generate second derivatives. We require the property of the epsilon function that it is diagonal in $k$-space:

[TABLE]

where the Fourier transform is defined as

[TABLE]

a simplification that we will make use of below. Using

[TABLE]

we rewrite (3.2) as a stochastic Ehrenfest equation with additional terms due to the projector

[TABLE]

where the noise correlations are

[TABLE]

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 trace222We refer to $dA^{\varepsilon}$ loosely as the energy-damping trace term, despite the fact it is not strictly a trace; it is the energy-damping counterpart to the number-damping trace term, which is a true trace. 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 explore the consequences of this formulation, providing explicit stochastic equations for particular observables of interest.

3.3 One-body operators

If we restrict our attention to one-body operators, the projected functional calculus assumes a particularly simple form. Consider a quantity that is the pure-state trace of a one-body operator

[TABLE]

where the operator $\hat{A}$ is independent of $\psi,\psi^{*}$. We emphasize that this quantity is spatial average for the pure state wavefunction associated with each trajectory, and does not represent the many-body expectation value found by tracing over the density matrix (the ensemble average). The non-zero projected functional derivatives are

[TABLE]

The projected functional derivatives corresponding to the operator $\hat{A}$ are equivalent to the regular functional derivatives of the totally projected operator $\hat{A}_{\cal P}\equiv\hat{{\cal P}}\hat{A}\hat{{\cal 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

[TABLE]

where the noise correlations are

[TABLE]

and the projector terms $dA^{\varepsilon}-Q_{A}^{\varepsilon}$ are given in Appendix A. We have expressed Eq. (43) in terms of the commutator $[\hat{A},\hat{H}_{{\scriptstyle{\rm sp}}}]$ to show the connection to the standard Ehrenfest relations.

3.4 Noise correlations for multiple moments

When considering more than one moment, one must take care when considering the noise terms that arise. Let $B=B[\psi,\psi^{*},t]$ be another moment of the SPGPE with evolution described by Eq. (37). The correlations between the noises corresponding to $A$ and the noises corresponding to $B$ are

[TABLE]

for number-damping and energy-damping respectively.

The special case of one-body operators allows significant simplification: for moments that are the expectation of a one-body operator the correlations reduce to

[TABLE]

for number-damping and energy-damping respectively.

3.5 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

[TABLE]

where the noise correlators are given by Eqs. (46),(47). 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:

1. i)

*Ensemble averages.—*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 [21]. Those Ehrenfest relations described the evolution of ensemble avarages, whereas in the present formulation we retain all noises in the collective equations. 2. ii)

*Multiplicative noise.—*Superficially, the multiplicative noise in the SPGPE Eq. (28) may appear to have been transformed into additive noise. However, we emphasize that multiplicative noises remain present in terms of the form $\sqrt{\langle\hat{A}\rangle}dW_{j}$, since $\langle\hat{A}\rangle$ is not a noise average. 3. iii)

*Additive noise.—*Reduction to additive noise can be achieved in special cases where the system can be well-described 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 [19]. 4. iv)

*Projector terms.—*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. 5. v)

*Thermal equilibrium.—*Formally it is known that trajectories of the SPGPE will approach a grand canonical ensemble [10] describing the thermal equilibrium properties of the system. The c-field ensemble is then drawn from the equilibrium Wigner functional

[TABLE]

In the SPGPE, the particular ensemble depends upon the ratio ${\cal M}/\gamma$, with canonical equilibrium reached in the limit ${\cal M}\gg\gamma$, and grand-canonical equilibrium reached in the opposite limit. The equilibrium ensemble satisfies a set of fluctuation-dissipation relations, as may be derived by applying the steady-state condition to Eq. (50). In each case there is a balance between thermal noise and the damping. For example, for pure number damping it is easily shown that $\langle\langle L-\mu\rangle\rangle=k_{\rm B}T{\cal N}$ for ${\cal N}\equiv\textrm{Tr}\;\hat{{\cal P}}$ modes in the $C$-region. The correct balance between damping and noise is manifest through the absence of the dissipation parameters from such expressions, consistent with their absence from the equilibrium ensemble Eq. (51). We note that ensemble (in)equivalence [27] and the choice of an optimal cutoff [28] are both active areas of research.

4 Application to centre of mass fluctuations

As a test of the formalism we consider centre of mass oscillations in a quasi-1D harmonically trapped system with frequency for transverse trapping $\omega_{\perp}$ such that $\omega_{\perp}\gg\omega$, the frequency for the axial trap. Our primary aim is to test the validity of our analytical approach against c-field simulations. The system we choose is particularly simple and thus amenable to analytic treatment of the SERs. However, we must take some care in interpreting the reservoir theory in this case, as the Bose gas in a harmonic trap obeys Kohn’s theorem.

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. While the second-quantized field theory satisfies Kohn’s theorem, the best available c-field reservoir theory describes the $I$ region as time independent, and hence violates Kohn’s theorem. C-field theory is thus currently best suited to systems for which a time-independent high-energy thermal reservoir is a reasonable approximation. We emphasize that the low-energy c-field is a dynamical treatment of both condensate and non-condensate; above a high-energy cutoff the dynamics are typically not included due to technical limitations.

While the time-independent reservoir approximation (TIRA) is not strictly applicable for scalar BEC in the purely harmonic trap, there are a number of physical systems where it is applicable: Kohn’s theorem does not apply to a scalar Bose gas held in a harmonic trap if the trap becomes non harmonic at high energy. The theorem is also inapplicable for a harmonically trapped system if the reservoir consists of second atomic species confined by a different trapping potential, as may occur during sympathetic cooling. Furthermore, any system that is not harmonically trapped will not obey Kohn’s theorem, and is thus potentially amenable to the TIRA. Examples systems in non-harmonic traps for which the theory is applicable include vortex decay in hard-wall confinement [46], soliton decay in a 1D toroidal trap (where the present approach was first used) [19], and persistent current formation in a 3D toroidal trap, where SPGPE simulations [17] compare well with experiment [44].

In this work our approach is simply to test the formalism on a simple model system within the TIRA by integrating out the spatial degrees of freedom to find effective stochastic equations of motion for the centre of mass. We stress that the TIRA approached used here is physically valid for (at least) two scenarios if immediate interest: non-harmonic trapping at high energies for a scalar BEC, and sympathetic cooling involving two BEC components in different harmonic traps [14]333The theory introduced here can be developed into a first-principles treatment for a particular atom number and temperature by self-consistent calculation of the reservoir interaction parameters [3, 17, 13].. An effective one-dimensional theory can be found in the prolate trapping regime where the reservoir remains three dimensional444Without this restriction the complexity of the reservoir interactions increases significantly and we do not pursue this regime further here., and the low-dimensional subspace is well-described by projecting onto the transverse ground state of the confining potential.

The one dimensional SPGPE assumes an identical functional form to the three dimensional SPGPE, but with dimensionally reduced damping and noise coefficients [18]. The number-damping term is only changed by the replacement $\varg\to\varg_{\rm 1D}$, as are the Hamiltonian terms. The scattering kernel in the 1D reduction of Eq. (26) becomes

[TABLE]

where $\textrm{erfcx}(q)\equiv e^{q^{2}}\textrm{erfc}(q)$ is the scaled complementary error function, and $a_{\perp}\equiv\sqrt{\hbar/m\omega_{\perp}}$ is the transverse harmonic oscillator length, much smaller than $a_{\omega}=\sqrt{\hbar/m\omega}$. The system is assumed to be sufficiently condensed that the centre of mass motion can be approximated using a Thomas-Fermi wavefunction ansatz. We do not require that this is a good approximation — indeed, the high-temperature regime requires that it must not be. We merely require that it should approximate the centre of mass motion. The Thomas-Fermi wave function allowing for arbitrary variations in the centre of mass position $x(t)$ and momentum $p(t)$ is

[TABLE]

Using this ansatz, the integrals for Eqs. (50) can be evaluated analytically.555We note that in the basis of harmonic oscillator modes used to represent the $C$-region the trace terms in the Ehrenfest relations for position and momentum vanish.

4.1 Stochastic equations for centre of mass

In the interest of brevity, we outline the derivation here, providing details in Appendix B. We first use the Thomas-Fermi ansatz Eq. (53) and evaluate the spatial integrals in (50), (50) to derive stochastic equations for $x$ and $p$ that are exact within the Thomas-Fermi approximation. The projector terms, given in Appendix B.1, severely constrain exact analytic progress beyond this point. Provided these terms may be safely neglected, we can arrive at an approximate set of equations that can be solved analytically. For consistent c-field simulations, neglecting the projector terms in the SERs 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. It is convenient to rewrite the equations of motion as a single equation of motion for the dimensionless complex variable

[TABLE]

and neglect terms of higher order than linear in $z(t)$, $z^{*}(t)$, to find the simpler stochastic equations given in Appendix B, Eq. (B). Assuming it is physically consistent to neglect projector terms in (B), (77a), we arrive at the equation of motion

[TABLE]

with noise correlations

[TABLE]

The equation of motion describes simple harmonic motion with two sources of damping and noise.

4.2 Analytic and numerical solutions

In this section we present an analytical treatment of the stochastic equations discussed in Section 4.1, and compare the results with numerical simulations of the 1D SPGPE. The analytic approach offers a simple intuitive picture of the dynamics of the centre of mass, for both number and energy damping mechanisms, and for both dissipation and noise.

Neglecting the projector terms, we can write the coupled differential equation as a vector SDE representing an Ornstein-Uhlenbeck process666The defining property of an Ornstein-Uhlenbeck process is that the drift and diffusion matrices are independent of the stochastic field. Since the noise is additive, the distinction between Ito and Stratonovich is irrelevant for Ornstein-Uhlenbeck processes.

[TABLE]

where $\mathbf{u}(t)=[x(t),p(t)]^{\intercal}$, and

[TABLE]

are the drift and diffusion matrices respectively, and $d\mathbf{W}(t)=[dW_{1}(t),dW_{2}(t)]^{\intercal}$ is a vector of independent real Wiener processes with correlations $\langle dW_{n}(t)dW_{m}(t)\rangle=\delta_{mn}dt$. The SDE has the formal solution

[TABLE]

where have assumed the initial state ${\mathbf{u}}(0)$ to be deterministic. The mean undergoes exponential decay $\langle{\mathbf{u}}(t)\rangle=\exp\left[-\boldsymbol{\Lambda}t\right]{\mathbf{u}}(0)$.

For a system with $x(0)=x_{0}$ and $p(0)=0$, the centre of mass position over time is given by

[TABLE]

where we have defined the frequency $\omega_{\gamma\varepsilon}=\sqrt{\omega^{2}-(\Lambda_{\gamma}-\Lambda_{\varepsilon})^{2}}$. This simple result has a clear interpretation in terms of the decaying oscillatory motion of the centre of mass with effective damping rates given by the SPGPE. Note that since we have arrived at a description with additive noise, the noise does not alter the average coordinate.

To study this regime of energy-damped dynamics numerically, we consider the steady-state correlations given by

[TABLE]

Fourier transforming with respect to $\tau$ in the steady-state, the fluctuation spectra are given by

[TABLE]

The steady-state correlations for position-position, momentum-momentum, and momentum-position are

[TABLE]

respectively. The steady-state spectra for position-position, momentum-momentum, and momentum-position are

[TABLE]

respectively.

We now test our analytical description against SPGPE simulations. While the SPGPE can be used for quantitative modelling [29, 17, 30] here our aim is simply to test our analytic solutions of the SPGPE using the SERs. We thus choose parameters that are representative of BEC experiments, and leave a first principles treatment of reservoir interactions [18] for future work.

We perform simulations of the 1D SPGPE [18] with the initial condition given by Eq. (53) with $x(0)=p(0)=0$. We use a chemical potential of $\mu=100\hbar\omega$, an energy cutoff of $\epsilon_{\rm c}=250\hbar\omega$, a temperature of $T=500\hbar\omega/k_{\rm B}$, an interaction strength of $\varg_{\rm 1D}=0.01\hbar\omega a_{\omega}$, an energy-damping rate of $\mathcal{M}=0.0005a_{\omega}^{2}$, and a number-damping rate of $\gamma=0.001$. We chose values for the damping rates such that $\Lambda_{\varepsilon}\approx\Lambda_{\gamma}$ and thus neither damping process is dominant over the other. Timescales are considered in units of the harmonic oscillator time period $t_{\omega}\equiv 2\pi/\omega$.

We compare our analytic solutions with the numerical data from equilibrated SPGPE simulations. In Fig. 2 we show the steady-state correlation functions for an ensemble of 5000 trajectories777When calculating the correlations from the numeric data, we have assumed that the system has reached equilibrium after five trap cycles $t=5t_{\omega}$, and used ergodic averaging over the remaining time interval $t=5t_{\omega}$. With respect to the dissipation timescale $(\Lambda_{\gamma}+\Lambda_{\varepsilon})^{-1}$, the timescale of averaging is equivalent to $t=2.78(\Lambda_{\gamma}+\Lambda_{\varepsilon})$. We see that the analytic and numeric results show excellent agreement for short times with differences becoming more pronounced for larger $\tau$. Similarly, the steady-state spectra for position-position, momentum-momentum, and momentum-position are shown in Fig. 3, where the numeric spectra are obtained using the Wiener-Khinchin theorem applied to the numeric steady-state correlations. The bulk oscillation seen in Fig. 2 is the Kohn mode, as seen from the peak at $\Omega\equiv\omega$ in Fig. 3. The spectral linewidth represents the decay time of the two-time correlation function, determined by $\Lambda_{\varepsilon}$ and $\Lambda_{\gamma}$, 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.

5 Discussion and Conclusions

5.1 Discussion

The SERs 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 [31], where a number of useful functional relations were derived for Wigner phase-space methods [6] 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 [32], has been used in several studies of finite-temperature BEC dynamics [33, 34, 35, 36, 37]. 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. [38] 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 energy-damping mechanism is absent from the Stoof SGPE, an interesting future direction is to identify further experimentally accessible regimes capable of distinguishing between the two approaches [30]. A related question receiving recent attention is that of ensemble equivalence [27]; in general the equilibrium ensemble generated by the SPGPE will be affected by the relative strengths of number and energy damping 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. The relative occupation at the cutoff thus determines the presence of spurious dynamics due to the projector, and this population can typically be chosen to be of order unity. Ultimately, large projector corrections generate spurious dynamics and for a consistent treatment the cutoff should be chosen such that the projector corrections are small.

5.2 Conclusions

We have developed a set of exact stochastic Ehrenfest relations (SERs) for the complete stochastic projected Gross-Pitaevskii equation [10], an equation of motion which has significant application in the study of finite-temperature Bose-Einstein condensates [3, 11, 9]. In addition to the number-damping process, the SPGPE contains a number-conserving dissipative mechanism that transfers energy between the system and reservoir. Our main result, the SERs, retain the stochastic nature of the SPGPE and explicitly contain terms that result from both number- and energy-damping processes.

We tested our stochastic Ehrenfest equations in two ways. Considering the centre of mass 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. Our chosen test system has the weakness that the centre of mass motion in a purely harmonic trap is not strictly amenable to the reservoir theory due to a violation of Kohn’s theorem, However, our treatment is physically relevant for non-harmonic trapping, multicomponent systems, and other systems of interest that physically violate Kohn’s theorem, provided a low-energy fraction is harmonically trapped. Indeed, since the thermal equilibrium properties involve small excursions from equilibrium, the confining potential is only required to be locally harmonic near the trap minimum and any number of non-harmonic effects may intrude at larger distances. The centre of mass motion thus provides an excellent formal and numerical test of the SERs, being one of the simplest states of motion to handle analytically.

We have shown that SERs can be used to obtain analytic equations that agree with numerical solutions of the full SPGPE and offer some physical insight into the open system dynamics. Future work will explore systems involving analytically tractable excitations such as vortex decay in hard-wall confinement [46], soliton [39] and phase-slip dynamics [40] in toroidal confinement, sympathetic cooling [41, 42], spinor BECs [43], and quantum turbulence in non-harmonic traps [44, 45, 46, 47].

Acknowledgements

We thank Matthew Reeves and Sam Rooney for valuable discussions.

Funding information

ASB was supported by the Marsden Fund (Contract UOO1726), and the Dodd-Walls Centre for Photonic and Quantum Technologies.

Appendix A Projector terms

A.1 General operators

For general operator $\hat{A}$, the drift projector terms are

[TABLE]

the number-damping, and energy-damping diffusion projector terms are

[TABLE]

and the epsilon term is

[TABLE]

A.2 One-body operators

For the special case where $\hat{A}$ is a one-body operator, the projector corrections take a simpler form. In this case we find the drift projector terms

[TABLE]

diffusion projector terms

[TABLE]

and the epsilon term

[TABLE]

Appendix B Centre of mass equations

Integrating over the Thomas-Fermi ansatz, the SERs for position and momentum take the form of a pair of coupled stochastic differential equations for the centre of mass position and momentum

[TABLE]

Here the noise correlations are

[TABLE]

and the drift and diffusion rates are

[TABLE]

derived from the effective one dimensional SPGPE [18] by integrating over the Thomas-Fermi ansatz. In terms of the dimensionless phase-space amplitude Eq. (54) the equation of motion becomes

[TABLE]

where the noises have non-zero correlations

[TABLE]

and the projector terms are given below. The projector terms are relatively simple for this SPGPE describing a BEC in an oblate parabolic trap with only one effective $C$-region dimension [18]. We first consider the position and momentum terms separately, and then give the terms for the dimensionless complex variable Eq. (54).

B.1 Position and momentum

Noting that

[TABLE]

the projector and epsilon terms for $x$ and $p$ are

[TABLE]

and

[TABLE]

B.2 Dimensionless variable $z(t)$

The projector corrections for the dimensionless variable phase-space amplitude $z(t)$ are

[TABLE]

B.3 Numerical evaluation of projector terms

With the exception of the epsilon correction Eq. (82g), all the cutoff terms are a result of mode mixing between the highest energy coherent mode and the lowest energy incoherent mode. In the limit of imposing a very high energy cutoff, the populations of the modes above the cutoff approach zero and hence all the cutoff terms go to zero. For a well-chosen but finite cutoff the integrals involving the overlap of the lowest energy incoherent mode $\phi_{n_{\rm c}+1}(x)$ and the coherent field wave function $\psi(x)$ should be small, as $\phi_{n_{\rm c}+1}(x)$ is highly oscillatory and the mode population is also small by definition. We claim that for a well-chosen cutoff the cutoff terms are small enough such that they may be neglected, and we justify this in two ways. Firstly, we consider the magnitude of a selection of the cutoff terms by calculating them numerically. Second, we consider the analytic solutions that can be found by neglecting the cutoff terms and show that these agree well with simulations of the 1D SPGPE.

When considering the cutoff terms involving mode mixing at the cutoff (i.e. all except the epsilon correction), we note that only the Hamiltonian cutoff term $q_{z}^{H}$ Eq. (82a) does not have one of the damping rates as a multiplying factor. As the damping rates have a typical value several orders of magnitude less than unity, it is reasonable to expect that of all these terms, the Hamiltonian cutoff term will be the largest. If we show that $q_{z}^{H}$ is small enough to be neglected then we can reason that $q_{z}^{\gamma}$ Eq. (82b), $q_{z}^{\varepsilon}$ Eq. (82c), $d_{z}^{\gamma}$ Eq. (82d), $d_{z}^{\varepsilon}|_{(1)}$ Eq. (82e), and $d_{z}^{\varepsilon}|_{(2)}$ Eq. (82f) are smaller still and so can certainly be neglected also.

The epsilon term $dz^{\varepsilon}$ Eq. (82g) is distinct from the other terms, as it is not a result of mode mixing at the cutoff. We expect the two terms in the epsilon term to almost cancel, as it is clear this is the case in the limit that the projector becomes the identity. This is an important result of the earlier step where we found that writing the SPGPE in Ito form resulted in an extra term; without this extra term the epsilon term would in general be non-negligible, and in fact diverges in the infinite cutoff limit. We thus monitor the magnitudes of $q_{z}^{H}$ and $dz^{\varepsilon}$ over the course of an ensemble of numerical trajectories.

We define the relative cutoff term magnitudes by

[TABLE]

noting that $|\dot{z}(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. 4, where the initial state is Eq. (53) 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 $\sim 0.1$ early in the dynamics, we see that once the system has equilibrated they approach a steady-state of order $\sim 0.01$.