Temporal non-equilibrium dynamics of a Bose Josephson junction in presence of incoherent excitations

The time-dependent non-equilibrium dynamics of a Bose-Einstein condensate (BEC) typically generates incoherent excitations out of the condensate due to the finite frequencies present in the time evolution. We present a detailed derivation of a general non-equilibrium Green's function technique which describes the coupled time evolution of an interacting BEC and its single-particle excitations in a trap, based on an expansion in terms of the exact eigenstates of the trap potential. We analyze the dynamics of a Bose system in a small double-well potential with initially all particles in the condensate. When the trap frequency is larger than the Josephson frequency, $\Delta>\omega_J$, the dynamics changes at a characteristic time $\tau_c$ abruptly from slow Josephson oscillations of the BEC to fast Rabi oscillations driven by quasiparticle excitations in the trap. For times $t<\tau_c$ the Josephson oscillations are undamped, in agreement with experiments. We analyze the physical origin of the finite scale $\tau_c$ as well as its dependence on the trap parameter $\Delta$.


I. INTRODUCTION
Time-dependent phenomena are at the heart of the physics of ultracold atomic gases and, in particular, Bose-Einstein condensates (BEC). The creation of BECs in optical traps [1,2] has opened a box of Pandora for the study of temporal dynamics. This is because the relevant time scales in these systems allow for direct access to the time evolution, like quench dynamics [3,4] or DC and AC Josephson oscillations [5,6], but also because most transport experiments are done by applying time-dependent external fields, like ratchet potentials [7], or by expansion experiments [1,2], since stationary transport cannot be driven easily in the neutral, finite-size systems.
There are several powerful numerical ways which target dynamical behaviour of superfluids. For instance, in one-dimensional systems, quench and relaxation dynamics as well as temporal spread of correlations have been treated by time-dependent density matrix renormalization group and nonequilibrium dynamical mean-field theory for bosons [8,9]. Also the exact quantum dynamics of a 1D bosonic Josephson junction (BJJ) has been investigated by solving many-body Schrödinger equation of motion by using the multiconfigurational time-dependent Hartree for bosons (MCTDHB) method [10]. This work demonstrated that methods based on the standard semiclassical Gross-Pitaevskii description are insufficient for capturing the complicated dynamics of the BJJ, and a many-body approach is necessary [10]. This is because a time-evolving system will always create incoherent excitations above the BEC (which we will call quasiparticle (qp) excitations), even if it was initially prepared with all particles in the BEC, provided the frequencies characterizing its time dependence are larger than the lowest qp excitation energy, characterized by the trap frequency. If the BEC is initially prepared in an excited state, for instance, by a spatial modulation of the BEC amplitude or by an occupation imbalance of a BEC in a double-well potential, qp excitations out of the condensate are necessary to thermalize the system, i.e., to reach a thermal distribution of all states and to accommodate the entropy of the thermal state. Hence, one may expect that even in a non-integrable system one can prevent thermalization by making the trap small enough, so that the trap frequency exceeds the characteristic frequency of the time evolution. Multiple undamped (non-thermalizing) [5,11] as well as strongly damped (thermalizing) [12] Josephson oscillations of BECs in double-well potentials have been observed experimentally in traps with different sizes and particle numbers. The above discussion exemplifies the potential necessity of considering incoherent excitations in any temporal evolution of a BEC.
A first, general description of Bose many-body systems beyond the Gross-Pitaevskii means field dynamics was developed in the seminal works by Beliaev [13], Kadanoff and Baym [14], Kane and Kadanoff [15], and Hohenberg and Martin [16]. Kinetic theory which includes collisions between then condensate and non-condensate part of the system was developed by Griffin, Nikuni and Zaremba [17]. Further approaches to a quantum kinetic description can be found in Refs. [18][19][20] and references therein. However, the wide applicability of these methods is somewhat hampered by the large number of variables inherent to the spatial representation, and/or (as well as in numerical methods [10]) suitability for only moderate number of atoms.
In the present paper we lay out a general, tractable method for the quantum dynamics of Bose-Einstein condensed systems in the presence of non-condensate excitations. It is based on a systematical expansion of the bosonic field operator in the eigenstates of a noninteracting Bose system in an arbitrary trap potential [21], where the operators for the lowest-lying states are replaced by their respective, time-dependent condensate amplitudes and the full quantum dynamics of the excited states is retained. We derive the coupled equations of motion for the classical condensate propagator and the quantum propagator of the excitations. By using this eigenstate representation, the spatial dependence is cast into the eigenstate wave functions and the number of degrees of freedom can be limited to the number of trap states to be considered, thus keeping the quantum problem tractable even for complicated trap potentials, depending only on a few system parameters. We can treat any number of particles.
While our formalism is completely general, in this work we apply it to a system of an interacting Bose gas confined in a double-well potential with an initial population imbalance of the BECs in the left and in the right potential well. This system is established to exhibit Josephson oscillations [5,6,22]. The Josephson effect, well known from superconductors [23,24] and superfluids [25,26] has its perculiarities in BEC systems [22,[27][28][29][30][31]. Apart from already mentioned experiments it has been also observed in 1D Bose Josephson junctions on a chip [32] and arrays of Bose-Einstein condensates [33]. The BJJ is usually described in terms of the oscillating condensate population imbalance where N 1 (N 2 ) is the number of particles in the left (right) well. The dynamics of z(t) is very sensitive to interactions. One distinguishes two regimes: a delocalized regime with the time average z(t) = 0, and a self-trapped (ST) regime with z(t) = 0 [5,22]. The transition from the delocalized to the self-trapped regime occurs when the initial population imbalance z(0) exceeds a critical value z c (0), or equivalently, the interaction strength between bosons becomes greater than a certain critical value, which depends on z(0).
Once incoherent single-particle excitations are taken into account, the semiclassical dynamics of the BJJ drastically changes [10,21,34,35]. First of all, the self-trapped state is destroyed by qp excitations [10,21]. Secondly, we have shown that the nonequilibrium qp dynamics does not set in instantaneously as the Josephson coupling is switched on, but with some delay, described by a time scale τ c [21]. For t > τ c quasi-particles are excited in an avalanche manner leading to fast Rabi-like oscillations of the qp occupation numbers and of the condensates between the wells. We provide a detailed analysis of the origin of τ c and its dependence on systems parameters.
The paper is organized as follows. In Section II we present the model Hamiltonian, and main parameters of our problem are derived. In Section III we unveil the general formalism in terms of Keldysh Dyson equations for condensed and noncondensed particles.
In Section IV we apply our formalism to a double-well potential and consider the effect of quasi-particles on the BJJ within Bogoliubov-Hartree-Fock approximation. We discuss our results for particle imbalances, and inverse τ c versus various system parameters in Section V. Main conclusions and perspectives are given in Section VI.

II. THE MODEL HAMILTONIAN
A system of weakly interacting bosons, trapped in an external double-well potential ( Fig. 1), is most generally described by the Hamiltonian whereΨ(r, t) is a bosonic field operator, and we assumed a contact interaction between the bosons with g = 4πa s /m (a s is the s-wave scattering length). V ext is the external double-well trapping potential.
We now expand the field operatorΨ(r, t) in terms of a complete basis {ϕ − , ϕ + , ϕ 1 , ϕ 2 ...ϕ M } of the exact single-particle eigenstates of the double well potential V ext (x, t > 0) after the coupling between the wells is turned on. Here ϕ − is the ground state wave-function of the double well, ϕ + is the first excited state wave function, ϕ 1 is the second excited state and so on. The field operator hence readŝ where the first two terms correspond to the usual two mode approximation for a condensate in a double well [22,29,31] with φ 1 (r) = (ϕ − (r) + ϕ + (r))/ √ 2 and φ 2 (r) = (ϕ − (r) − ϕ + (r))/ √ 2. Since the ϕ + (r), ϕ − (r) are (approximately) parity symmetric or antisymmetric about the center of the double-well, respectively, the φ 1 (r) and φ 2 (r) are localized wave functions in the left or right well. The destruction operators in these left and right levels are approximated by their complex amplitudes, α = 1, 2, where N α and θ α are the number of particles and the phase of the condensate in the left (right) well of the potential. This semiclassical treatment of the BECs neglects phase fluctuations. It is justified when the BEC particle number is sufficiently large, N α ≫ 1, e.g., for the experiments [5]. The applicability of the semiclassical approximation has been discussed in detail in Refs. [34][35][36] and has been tested experimentally in Ref. [37]. The bosonic destruction operatorsb n refer to bosons out of the condensate, their quasiparticle dynamics will be treated fully quantum mechanically. For simplicity we assume that the wave functions ϕ n extend over both wells and are real. In the numerical evaluations we limit the number of levels which can be occupied by the quasi-particles to M = 5.
At time t = 0 the barrier between the wells is suddenly lowered, and a Josephson link is established between the two condensate reservoirs. The Hamiltonian for our system for Finally, mixing between the BECs and the quasiparticle system is described by Here It is convenient to introduce two interaction constants so that processes proportional to J ′ describe quasiparticle assisted Josephson tunneling.
H mix is then Hereafter all our energy scales are counted from E 0 . The main parameters in our model are, hence, the interlevel spacing ∆, the quasiparticle assisted Josephson coupling J ′ and the interaction K (Fig. 1).
In the next chapters we proceed to the derivation of equations of motion for the condensates and the quasiparticles described by the quantum operatorsb,b † .

A. Quantum kinetic theory: Dyson equations
In this Section we formulate the quantum kinetic theory of the Bose Josephson junction in terms of nonequilibrium Keldysh Green's functions [39] defined along the Keldysh time contour. To shorten the notations we denote the classical part of our field operator (3) as Ψ 0 (r, t) and the quantum part asΨ(r, t). We define correspondingly two single particle Green's functions, one for the quasiparticles G and one for the condensed particles C.
, and T C denotes time ordering along the Keldysh contour, i.e. each of the bosonic Green's functions G and F is a 2 × 2 matrix in Keldysh space. The condensate propagator C is classical with trivial time ordering, For the derivation of the equation of motion for G and C, we proceed in the standard way (see, for instance Refs. [17,40]). We first write down the associated Dyson equations for both Green's functions on the Keldysh contour, which read as follows where C d2 ≡ dr 2 C dt 2 and the subscript C denotes the time integration along the Keldysh contour. For convenience we decompose the self-energies into the first order Hartree-Fock contribution S HF (Σ HF ) and the higher-order part S (Σ) which accounts for collisions.
Collisional self-energies, which are responsible for dissipation and equilibration can be derived within self-consistent second-order approximation.
The bare 2 × 2 propagator G 0 is defined as where and 1 is the 2 × 2 identity matrix.
The Hartree-Fock self-energies correspond to the one-particle irreducible diagrammatic contributions of the first order in the interaction strength g Parametrizing the Keldysh contour in terms of real time variable we obtain from Eqs. (20) and (21) Here we introduced the single particle excitations spectral function and analogously and for the corresponding self-energies not containing the Hartree-Fock contributions. The "lesser" and "greater" Green's functions are defined as usual [40] It is also convenient to introduce the so-called statistical propagator (which is just a Keldysh component of the quasiparticle Green's function G K (1, 1 ′ ) divided by two), and rewrite our equations in terms of the spectral function and the statistical propagator which contain the information about the spectrum and the occupation number, respectively.
We hence get and where Π is defined as Eqs. (34) and (35) were obtained by taking the difference and the sum of Eq. (27) for the greater and lesser propagator, and then inserting the definition of the spectral function (28) and the statistical function (33), respectively. The Hartree-Fock self-energies in terms of the spectral and statistical functions read where we used It is worth mentioning that these self-energies only contain the symmetrized two-point propagators instead of the time-ordered ones, the non-condensate and anomalous density are then expressed in terms of the propagators G < and F < evaluated at equal points in space and time correspondingly.
We now derive equation of motion for the condensate and the single-particle excitations from our general Dyson equations.

B. Equation of motion for the condensates
Using the bosonic field operator (3), we straightforwardly obtain for the propagators (18) and (19), the spectral function (28) and the statistical function (33) the following expressions where (41) and (42) in the Dyson equation for the condensate Green's function (26) we obtain the Dyson equation for C αβ (t, t ′ ) where the space dependence is absorbed in the parameters E 0 , U, U ′ , K, J and J ′ , In this equation we adopt the sum convention for all greek indices, i.e., greek indeces appearing in a term twice are summed over from 1 to 2, and γ αβ = S > αβ − S < αβ accounts for collisions. The bare propagator is given by The Hartree-Fock self-energies read with α = 1, 2. In Fig. 2 where we used, that S HF αγ (t, t ′ ) = S HF αγ (t)δ(t − t ′ ). The equation of motion for the time dependent condensate amplitude a α (t) can be obtained by taking the upper left component of Eq. (51) and then dividing by a * β (t ′ ), where S HF αβ and W HF αβ are the upper left and the upper right components of S HF αβ defined in Eq. (50), i.e., Similarly, γ S αβ and γ W αβ are the upper left and the upper right components of γ αβ , respectively.

C. The equation of motion for the single particle excitations
We follow here the procedure analogous to the one described in the previous section for the evolution Eqs. (34) and (35) for non-condensate particle spectral and statistical functions. Inserting Eqs. (41) and (42) into Eqs. (34) and (35) where Π nm = (Σ > nm + Σ < nm )/2 and Γ nm = Σ > nm − Σ < nm and the sum convention for all latin indices is implied, where the sum runs over all positive integers. The bare single-particle propagator and the Hartree-Fock self-energy are given by where ǫ n , U ′ and K αβnm are defined in Eqs. (11), (12) and (14). The first-order contributions to the quasiparticle self-energy in Eq. (57) are shown in Fig. 3.
where we used that Σ HF nm (t, t ′ ) = Σ HF nm (t)δ(t − t ′ ). Computing the evolution of F nm for equal time arguments requires special attention. In order to do it properly, we take the sum and the difference of Eq. (59) and its hermitian conjugate, and evaluate it at equal times (see details in the next section). Note, that the components of A nm (t, t) are fixed for all times due to the bosonic commutation relations. One must also obey particle and energy conservation by applying conserving approximations -for details see Appendix.
Eqs. (51), (58) and (59) constitute the general equations of motion for the condensate and the non-condensate (spectral and statistical) propagators, respectively. They are coupled via the self-energies which are functionals of these propagators and must be evaluated self-consistently in order to obtain a conserving approximation (see appendix). The timedependent Hartee-Fock self-energies describe the dynamical shift of the condensate and the single-particle levels due to the dynamical change of their occupation numbers and their interactions. The higher order interaction terms on the right-hand side of the equations of motion describe inelastic quasiparticle collisions. They are, in general, responsible for quasiparticle damping, damping of the condensate oscillations and for thermalization. In the present paper, we are interested in how the non-equilibrium condensate oscillations induce single-particle excitations and how the latter act back on the condensate dynamics. Hence, we do not consider quasiparticle collision effects here. The results with collisions will be published elsewhere [38].

IV. FINAL EQUATIONS IN TERMS OF BOGOLIUBOV-HARTREE-FOCK AP-PROXIMATION
We now neglect collisional terms and treat the nonequilibrium Bose Josephson junction within Bogoluibov-Hartree-Fock approximation. From Eqs. (52) and (59) we get It is not necessary to consider the equations for the spectral function A nm , since the only non-condensate related quantities appearing in the self-energies S HF αβ and Σ HF nm are the components of the statistical function F nm . This implies that eq. (60) and (61) decouple completely from the equations for A nm (this will not be the case when we take collisions into account).
Taking the difference of Eq. (61) with its hermitian conjugate we get The sum of Eq. (61) with its hermitian conjugate gives for its upper right component. Here we assumed the notations The self-energies Σ HF nm (t) and Ω HF nm (t) are given by where α = 1, 2 and N α = a * α a α . Evaluation of eq. (62) and (63) at equal times gives Note, that in our previous work we considered a special case of the functions F G nm and F F nm , or their equivalents, being diagonal in the energy space, i.e. being proportional to δ nm . In this work we do not adhere to this approximation and solve Eqs. (67) and (68) for all n and m.
Plugging the expressions for S HF αβ and W HF αβ into Eq. (60) we obtain The equation for a 2 (t) is obtained from Eq.(69) by replacing a 1 by a 2 and visa versa.
We now solve the system of equations (67), (68) and (69) numerically. We assume that initially a α (0) = N α (0)e iθα(0) , with z(0) = 0.6 and initial phase difference θ 1 (0) − θ 2 (0) = 0, and the total particle number equal to 500000 (we have also solved for different particle numbers, for instance, for N = 5000, but did not notice any crucial effect produced by the different particle number -see Fig. 10). It is also convenient to renormalised the other parameters [21]: u = NU/J, For an initially delocalised junction we assume u = u ′ = 5 and for an initially self-trapped we take u = u ′ = 25. We then analyse the behaviour of our system depending on k, j ′ and ∆ (∆ is expressed in units of J).

V. RESULTS
We now solve the system of equations numerically. We take the initial conditions with z(0) = 0.6 and initial phase difference θ 1 (0) − θ 2 (0) = 0. The total particle number is taken equal to N = N 1 (0)+N 2 (0) = 500000, with initially all particles in the condensate. It is convenient to define the dimensionless parameters [21] u = NU/J, k = NK/J, j ′ = NJ ′ /J, For an initially delocalised junction we assume u = u ′ = 5 and for an initially self-trapped junction we take u = u ′ = 25. We then analyze the behaviour of our system depending on k, j ′ and ∆ (∆ is expressed in units of J in the following).
In the absence of coupling to the qp excitations (J ′ = K = 0) or as long as the qp states are not occupied, the equations of motion (51), (58), (59) reduce to Eq. (52) and, hence, for a BEC in a double-well potential, to the well-known two-mode approximation for the BEC dynamics [22], exhibiting undamped Josephson oscillations [21]. Our main result is, however, the appearance of a new characteristic time scale τ c , which marks the onset of nonequilibrium qp dynamics. It may be surprising at first sight that for a discrete level spectrum τ c can take a non-zero value [21], i.e., quasiparticles are not excited immediately. This is because quasiparticles cannot be excited in low-order time-dependent perturbation theory, if the Josephson frequency is less than the (interaction-renormalized) level spacing, ω J < ∆ ef f , even though the initial state with population imbalance z(0) > 0 is a highly excited state of the system [21]. τ c depends on system parameters, in particular, the energy splitting ∆ and the quasi-particle assisted Josephson coupling j ′ . We now analyze in detail this new physics.
In Fig. 4(a) the particle imbalance of the initially delocalized junction is shown. One can clearly distinguish two regimes of slow and fast oscillations, respectively. Comparing with Fig. 4(c) reveals that the onset of fast oscillations is marked by an abrupt, avalanche-like population of the single-particle excited states, i.e., it coincides with τ c . For times t < τ c there is only a small, virtual qp population. In this regime, our formalism reduces to the two-mode approximation [22] for the condensates, and the dynamics resembles Josephson oscillations. One may, therefore, conjecture that the fast oscillation regime for t > τ c is dominated by Rabi oscillations of the single-particle occupation numbers. This will be confirmed below by identifying the oscillation frequencies as the Rabi frequencies for t > τ c .
The drastic change is also seen in the phase portrait of the junction in Fig. 4(b), where the point on the phase trajectory corresponding to τ c is marked by a thick dot. Before τ c the time average of the phase difference is equal to zero, as expected in the semiclassical case.
However, after τ c the time average of the phase difference is equal to π. We refer to this kind of behaviour as to "0 − π" crossover. It may be understood heuristically as follows. One can think of the two-mode BEC system as a driven oscillator, where the qp system provides the driving force. For t < τ c this "oscillator" has no driving force and oscillates at its resonance frequency, i.e., the Josephson frequency ω J [21,22], with a phase shift θ = 0. However, once in the regime for t > τ c , the oscillator is periodically driven at the Rabi frequency, far above its resonance frequency, resulting in a phase shift of θ = π.
We observe similar behaviour for the initially self-trapped Bose Josephson junction, shown Eventually, it will become impossible to excite quasiparticles and τ c will tend to infinity. We note in passing that the uppermost (5th) qp level is pushed upward by a level repulsion effect for t < τ c , but is pulled down by interactions once it gets populated.
The afore-mentioned fast oscillations can now be identified as Rabi oscillations by comparing their frequency with the renormalized level spacing ∆ ef f . To that end we plot in t > τ c , respectively, and also the Fourier transform of the BEC level E c1 (t) for t > τ c .
Similar results are obtained for the other levels. It is seen that in the qp dominated regime, t > τ c , both the condensate and the qp levels oscillate predominantly with frequencies given by the renormalized level spacing, ∆ ef f , and with the energy spacing between condensate and the uppermost level, 5∆ ef f . Note that the uppermost level has a large oscillation amplitude, because its oscillation amplitude is not bounded from above by other levels. This unambiguously shows that the fast oscillations are Rabi oscillations.
We now study the detailed dependence of the characteristic time scale τ c on other parameters. The dependence of the inverse τ c on the interelevel spacing ∆ and quasiparticle-assisted Josephson coupling j ′ for the initially delocalised junction is shown in Fig. 8. One can see resonance behaviour of the τ c : it diverges for certain values of ∆. The first divergence happens at ∆ = ∆ (1) res , which depends on j ′ . It turns out that this ∆ E c1 and E c2 oscillations. This period T cond = T jos /2 ≈ 1, where T jos is the oscillation period of particle imbalance z(t) before quasiparticles get excited (this is because E c1,2 depends on the absolute value of |z(t)| ). Since E c1 and E c2 oscillate in antiphase, τ c at ∆ (1) res corresponds to the minimal difference between the condensate energies. Next resonance happens when τ c ≈ T cond /2 + T cond . In some curves (for instance, curves 3,4,5 and 8) the third resonance can be also seen (it would correspond for τ c ≈ 5T cond /2).
For the initially self-trapped junction, only one resonance can be found. The value of τ c at ∆ res is equal to approximately 0.3, which is a half period of the self-trapped particle imbalance oscillations T ST = 0.6. This is the maximum value of τ c which we can achieve for these initial conditions. In this case we can also plot a phase diagram, i.e. ∆ res versus j ′ ( Fig. 10(a) ). Then below the curve τ c will take finite values, while above the curve τ c will be infinite and the junction will be dominated by the semiclassical dynamics. Note, that this phase diagram almost does not depend on the number of particles. Another interesting fact, is that ∆ res (j ′ ) − ∆ res (j ′ − 10) = 3ω ST , where ω ST = 2π/T ST ≈ 10.5. For completeness we also plot the j ′ -dependent values of ∆ (1) res of the initially delocalised junction ( Fig.10(b)). In Fig. 11 we present the dependence of the inverse τ c versus quasiparticle assisted coupling j ′ for various values of interaction k for the initially delocalised junction. One can see, when k is finite, it is more difficult to excite quasiparticles, and increasing k results in larger values of τ c . Interestingly, the curves for τ −1 c collapse onto a single curve by shifting j ′ by a k-dependent constant (see Fig. 12).

VI. CONCLUSIONS AND DISCUSSION
We studied in detail the nonequilibrium dynamics of a double well Bose Josephson junction in two cases: the junction is initially delocalised, and the junction is initially self-trapped in a semiclassical meaning described in Ref. [22]. The nonequilibrium is coursed by the abrupt switching on of the Josephson coupling between the well, which corresponds to the experimental situation [5]. We treat our system within Bogoliubov-Hartree-Fock first order approach and develop within it the nonequilibrium equations of motion for the condensate and quasi-particle parts of our Hamiltonian.
Our main finding is that the nonequilibrium dynamics does not set in immediately at time when the Josephson coupling is turned on, but after a certain period which we call τ c . Moreover, we find when it does set in at τ c , it happens in an avalanche way due to fast quasiparticle oscillations between the discrete energy levels of the system. The τ c depends on many of the system parameters and corresponds to a moment in time when the highest condensate eigenenergy is coincident with the eigenenergy of the first excited level. Since we find regimes with infinite τ c , this explains experimental finding that a junction can maintain undamped Josephson oscillations, although the barrier between the walls was ramped up in a sudden way [5].
We find also, that when nonequilibrium physics establishes, the junction inevitably switches to a π junction, a junction, whose time-dependent phase difference is equal to π when averaged over time. This is a necessary condition to sustain strong Rabi-like oscil-lations between the wells. In the future we aim to study how the nonequilibrium Josephson junction will equilibrate and thermalise (if at all), in which case we will consider collisions between the quasi-particles in the full second order approximation.
Appendix: Particle and energy conservation For any closed system the particle number and energy are conserved quantities. The particle number conservation arises as a consequence of the invariance of the Hamiltonian under global phase change. The mean particle number can be expressed in terms of the condensates wave functions and the statistical function, as follows and it can be proven to be constant by making use of the equations of motion for a α and F G nn . It is well known, that the exact solution of an isolated system exclude dissipation. One should therefore take care, that self-energies are derived within a conserving approximation (otherwise the non-physical dissipation can take place). If (like in our case) the self-energies S and Σ are derived by means of a conserving approximation, for instance from a Luttinger-Ward functional [41,42], the mean energy can be written as and can be proven to be conserved. The energy of the condensate fractions is H cond (t) = i 2 Tr E αβ 1 + 1 2 S HF αβ (t) C βα (t, t) + and the energy of the single particle excitations is H exc (t) = i 2 Tr ǫ n δ nm 1 + 1 2 Σ HF nm (t) F mn (t, t) Here again we sum over all indices appearing twice in every term.