Non-equilibrium Atomic Condensates and Mixtures: Collective Modes, Condensate Growth and Thermalization

The non-equilibrium dynamics of trapped ultracold atomic gases, or mixtures thereof, is an extremely rich subject. Despite 20 years of studies, and remarkable progress mainly on the experimental front, numerous open question remain, related to the growth, relaxation and thermalisation of such systems, and there is still no universally-accepted theory for their theoretical description. In this paper we discuss one of the state-of-the-art kinetic approaches, which gives an intuitive picture of the physical processes happening at the microscopic scale, being broadly applicable both below and above the critical region (but not within the critical region itself). Specifically, the Zaremba-Nikuni-Griffin (ZNG) scheme provides a self-consistent description of the coupling between the condensate and the thermal atoms, including the collisions between these two subsystems. It has been successfully tested against experiments in various settings, including collective modes (e.g. monopole, dipole and quadrupole modes), topological excitations (solitons and vortices) and surface evaporative cooling. Here, we show that it can capture two important aspects of non- equilibrium dynamics for both single-component and two-component BECs: the Kohn mode (the undamped dipole oscillation independent of interactions and temperature) and (re)thermalization leading to condensate growth following sudden evaporation. Our simulations, performed in a spherically-symmetric trap reveal (i) an interesting two-stage dynamics and the emergence of a prominent monopole mode in the evaporative cooling of a single component Bose gas, and (ii) the long thermalization time associated with the sympathetic cooling of a realistic two-component mixture. Related open questions arise about the mechanisms and the nature of thermalization in such systems, where further controlled experiments are needed for benchmarking.

The ability to experimentally control the many-body interactions [5], the temperature of the system [6], the trapping potentials [7,8] as well as the coherent coupling between different hyperfine levels [9], makes BECs an ideal system to probe fundamental problems. In addition, the increasingly precise manipulation also brings about a growing interest in the applications of BECs to create quantum devices, such as an atomic SQUID [10,11] or a matter-wave interferometer [12][13][14].
Modelling the finite-temperature dynamics of a condensate is therefore of significant value. For a theoretical model to be useful, it should capture the essential physics with minimal input parameters. At the same time, it must also remain easy to use by solving either analytically or numerically with feasible amount of computational effort [61]. In this sense, the 'Zaremba-Nikuni-Griffin' (ZNG) model [62][63][64] is very successful in modelling existing experiments on single-component Bose gases for a broad temperature range (typically T /T c < 0.9, or T > T c ), with potential corrections at very low temperatues arising from the quasiparticle nature of the excitations. By construction, this model treats the interactions between the condensate and the thermal cloud fully self-consistently, so including all collisional processes between them, and the respective back-actions during their coupled dynamical evolution. In that respect, it provides an intuitive picture of the dynamics happening at the microscopic scale and allows us to explore dynamics ranging from the collisionless evolution to the hydrodynamic regime. In particular, it takes into account three key aspects, which are not typically simultaneously accounted for in full by other finite-temperature models (such as 'classical' or 'c-field' methods) -see e.g the reviews [61,[64][65][66][67][68]: (i) the dynamics of the thermal cloud (for example, this is only approximately included in so-called 'classical field methods' [61,65,67,69,70]); (ii) the spatially-dependent collisions and dissipation (the latter is typically ignored in stochastic treatments, but see e.g. the related stochastic soliton decay [38]); and (iii) the conservation of the total number of atoms (note that this is guaranteed by construction here; for alternative explicitly number-conserving approaches, see Ref. [71] and references therein).
From a physics perspective, there are at least two stringent tests that a 'good' finite temperature nonequilibrium theory should satisfy, namely the ability to (i) reproduce the Kohn mode, i.e. the undamped dipole oscillations for a shifted trap, an exact feature occuring irrespective of the interaction strength (this is essential for precision measurement of collective modes); (ii) predict the full dynamical thermalization between the condensate and thermal cloud, achievable only through the fully self-consistent coupling of the dynamics of those two subsystems, and thus model the condensate growth process.
The ZNG model satisfies both criteria, being valid both below and above the critical temperature (in the latter case it just reduces to the usual kinetic Boltzmann equation) [64]. Here, it should be noted that, while the model does not account for the physics of the actual region of critical fluctuations close to T c (and so could not, for example, model phenomena associated with spontaneous defect formation a la Kibble-Zurek [43][44][45][46][47][48]) because of its explicit symmetrybreaking ansatz, it can nonetheless accurately predict condensate growth once a small condensate seed is added to the system [54], with the ensuing results being independent of the (small) seed size. The interesting and related topic of condensate formation starting from an ultracold thermal bosonic system, has been thoroughly reviewed in [72]. The extension of this ZNG model beyond the single-component Bose gas is rather involved. While the theoretical framework of the ZNG model has been established for a spinor condensate [73,74] (with coherent coupling between hyperfine levels) and two-component condensates [75,76] (with incoherent coupling between the two components), dynamical simulations of the two-component ZNG model have only been reported recently [77] and those for a spinor condensate remain unreported. Much work is still needed to perform systematic studies of two-component Bose gases, and to compare the experimental findings with the model predictions.
The aim of this work is both to demonstrate that the ZNG model can capture the essential physics of non-equilibrium Bose gases at finite temperature, and also to use it to identify various deeper physical questions, where further experimental work is required to understand in detail. Such questions include the thermalization process during/after evaporative cooling, whether such a process involves just one, or more, timescales, the extent to which a binary mixture actually thermalizes on experimentally-relevant timescales, and more broadly the process of sympathetic cooling, for which we are not aware of any modelling describing the coupled dynamics once both components start exhibiting condensation.
After presenting a brief introduction of the ZNG model (section 2), including its recent extension to a two-component Bose-Einstein condensates [75,76], we briefly explain how the model is solved numerically in practice in section 3. We then demonstrate the application of the model to study (i) the Kohn mode (section 4) and (ii) thermalization (section 5), for both the single-component and the two-component cases. These examples serve as stringent tests of the validity of the ZNG model in capturing the essential physics.

The Zaremba-Nikuni-Griffin (ZNG) Model
The ZNG model was first obtained for a homogeneous Bose gas by Kirkpatrick and Dorfman [78,79], and subsequently derived for an inhomogeneous singlecomponent Bose gas [62][63][64], an inhomogenous spinor gas [73,74] and an inhomogeneous two-component Bose-Bose mixture [75,76]. We present here a brief description of the two-component Bose gas model (which is reduced to the single-component Bose gas model if the inter-component interaction is switched off) in a way that it can also be easily understood by the non-experts.
We consider an interacting bosonic binary system described by the Hamiltonian and the two-body interactions are given bŷ whereΨ j ≡Ψ j (r) Ψ † j ≡Ψ † j (r) is the bosonic annihilation (creation) operator for an atom of component-j with mass m j , which obeys the usual commutation relationships for bosons, The s-wave collisions between atoms in different components are encompassed by g kj = 2π 2 a kj /m kj , where a kj defines the scattering length between atoms in components j and k, and m −1 defines the reduced mass. The atoms are confined in a harmonic potential V j (r) = 1 2 m j [ω 2 r,j r 2 + ω 2 z,j z 2 ] with radial and axial angular frequencies, ω r,j and ω z,j , respectively. For simplicity, we also assume isotropic traps (ω r,j = ω z,j = ω j ) in our analysis. Also, for subsequent discussions that involve only a single component, we omit the subscript j to simplify the notation.
Central to the ZNG methodology is the symmetrybreaking ansatz, and the Beliaev decomposition of the Bose field operatorΨ j into its average, non-zero, value φ j = Ψ j denoting the condensate wavefunction, and a fluctuation operatorδ j , where angular brackets . . . denote the broken-symmetry ensemble average ‡. It is precisely because of this decomposition (and the assumption of a non-zero expectation value for the Bose field operator) that the ZNG model cannot account for condensate growth from a purely thermal system (i.e. without introducing a numericallyconvenient non-zero seed). The condensate density n c,j and the thermal cloud densityñ j for atoms of component j are then given separately by the wavefunction, and the fluctuation operator via the diagonal noncondensate densitỹ A kinetic model is developed in [75,76], in which we identify the condensate field φ j (r) = Ψ † j (r) and the thermal cloud density as the only slowly-varying relevant quantities. The triplet anomalous averages δ † jδ jδj and δ † kδ kδj , as well as the off-diagonal noncondensate density δ † kδ j for k = j, are treated perturbatively via adiabatic elimination [80,81]. In addition, the pair anomalous averages δ jδj and δ kδj are discarded [82] as they do not generate energyconserving contributions (to order g 2 ). The end result is that a condensate field φ j (r) obeys a dissipative Gross-Pitaevskii equation while the Wigner distribution function of the thermal atoms obeys a quantum Boltzmann equation Compared to the usual zero-temperature case, the condensate atoms now experience an effective potential U j c corrected by the presence of the thermal cloud, while the thermal atoms are modelled as classical particles moving in an effective potential U j n . These potentials include both the external potential V j and the mean-field contributions, which are related to the condensate density n c,j (r) = |φ j (r)| 2 , and the thermal cloud densityñ j (r) = dp/(2π ) 3 f j (p, r, t), as U j c = V j + g jj (n c,j + 2ñ j ) + g kj (n c,k +ñ k ), (9a) U j n = V j + 2g jj (n c,j +ñ j ) + g kj (n c,k +ñ k ).
With these effective potentials, locally a condensate atom of component j has energy where is the non-equilibrium chemical potential for component j and v c,j = m j n c,j Im(φ * j ∇φ j ) defines the superfluid velocity of component j with momentum p j c = m j v c,j . On the other hand, a thermal atom of component j with momentum p has a Hartree energy Equation (9a) and (9b) encapsulate the mean-field effects between the condensates and the thermal clouds, as well as the mean-field effects between the different components. As the condensates and the thermal clouds evolve in time, the changes in densities cause the clouds to exert a force on each other through the mean-field potentials and equations (6) and (8), leading to a damping effect, without explicit consideration of collisions. As a consequence, calculating equations (9a) and (9b) dynamically allows us to simulate the collisionless evolution and study the Landau damping [83].
The collisional integrals C ·· 22 , C ·· 12 and C kj 12 in equation (8) are important to establish the full thermal equilibrium of the system starting from a nonequilibrium state, and they are responsible for the collisional damping. In particular, the thermal-thermal collisional integral, 5 7 (1 + δ kj ) dp 2 dp 3 dp 4 (14) between thermal atoms of component k and component j (including both k = j and k = j) vanishes when the thermal atoms are in local thermodynamic equilibrium, with the inverse temperature β = (k B T ) −1 , the normal fluid velocity v n,j and the normal fluid chemical potentialμ j , all as a function of position r and time t.
Here, k B is the Boltzmann constant. On the other hand, the thermal-condensate collisional integral, 2 4 (1 + δ kj )n c,k dp 2 dp 3 dp 4 2 4 (1 + δ kj )n c,j dp 2 dp 3 dp 4 (including both k = j and k = j) leads to a change in the number of condensate atoms of component j through the source term such that the total number of atoms (condensate + thermal cloud) remains unchanged. The C kj 12 integral, and consequently the source term R kj , only vanishes when the condensate atoms and the thermal atoms reach a local diffusive equilibrium (i.e. µ j c =μ j ) [63,84].
In contrast to a spinor BEC, a two-component mixture admits a condensate-exchange collision that comes from a perturbative treatment of the normal pair average δ † kδ j . The corresponding collisional integral C kj 12 = 2πg 2 kj n c,k n c,j dp 1 dp 2 for k = j and the source term describe a collisional process whereby (in the forward process) one thermal atom of component k collides with a condensate atom of component j and promotes the condensate atom to the thermal cloud while itself being cooled and condenses into the condensate. When evaluated with respect to realistic mixtures at thermal equilibrium, this condensate-exchange collision has dominant collisional rates [75,76], but its impact in dynamical situations remains a subject of study.

Solving the ZNG model in practice
Despite the two equations (6) and (8) that summarize the model being simple-looking, solving the equations for an arbitrary non-equilibrium situation is a numerically-challenging task and, to the best of our knowledge, very few groups have achieved it. Its first numerical implementation is carried out and documented by Jackson and Zaremba [85]. This code is subsequently adapted for the various studies being carried out at Newcastle [33-35, 39, 40, 42]. Arahata and Nikuni reported its application to study first and second sound in a highly-elongated trap [86]. A parallel version sped up with OpenMPI has been developed by Märkle [87] to study surface evaporative cooling [60]. We have recently developed a parallel version of the two-component ZNG model, for which we speed up the computation with OpenMP. We have used our new code to study the collective modes and rethermalization dynamics of both a single-component Bose gas and a two-component mixture, some of which are reported in this work. In particular, our recent application of the ZNG model [77] to study the counterflow dipole oscillation of a two-component mixture is the first of its kind and helps to identify the use of the dipole oscillation to map out the miscibleimmiscible transition of a mixture. Currently, we are developing a new numerical code for use on an Nvidia graphics processing unit. Most recently, Straastma at JILA has exploited the spherical symmetry and implemented the ZNG approach for an isotropic trap to study the damping of monopole oscillation below the critical temperature [88].
The basic algorithm to solve the ZNG model is outlined in [85]. We first obtain the equilibrium distributions at a finite temperature T . In this case, both the source terms and the collision integrals vanish, hence we can set φ j (r, t) = φ j (r)e −iµ j t/ with chemical potential µ j and use a semi-classical Hartree-Fock approximation for the thermal cloud and the textbook result of Bose function, We then obtain the condensate and thermal density profiles by solving equations (6) (through imaginarytime propagation) and (20) self-consistently. The equilibrium condensate wavefunctions and the thermal cloud densities are then set as the initial conditions of equations (6) and (8) with appropriate modifications. In a typical dynamical simulation, the dissipative Gross-Pitaevskii equation (6) is solved with the highlyefficient Fourier split-step method while the quantum Boltzmann equation (8) is solved with the direct simulation Monte Carlo (DSMC) method [89,90]. A large number of test particles (typically of the order of millions) are generated according to the Bose- These test particles are then evolved in time according to Newton's equation of motion using the symplectic leapfrog method. This solves the left-hand side of equation (8).
In order to simulate the collisions and to calculate the source terms, i.e. solving the right-hand side of equation (8) and estimating equations (17) and (19), the test particles are binned into spatial cells, where the cells have adaptive volumes [91,92] to improve the speed and accuracy of the computations. Collisions are simulated by randomly selecting pairs of test particles belonging to the same cell and checking if they are going to collide using the acceptancerejection method. The exact forms of the probabilities used in the acceptance-rejection method have been given in [64,85]. It is important to note that the same probabilities are also used to compute the source terms (17) and (19), hence accurate estimates of the source terms are accompanied by a large number of checks on collision events. In other words, while it is possible to speed up the simulation of the C ·· 22 collisions through the scaling of probabilities [90,92], this is no longer feasible for the thermal-condensate collisions (C ·· 12 and C kj 12 processes). From our experience, the quality of a full ZNG simulation that includes all collisional processes is largely determined by how well we have simulated the collisions and calculated the source terms.
In the next two sections, we will study the dipole oscillations and rethermalization of both a singlecomponent Bose gas as well as a Bose-Bose mixture, by solving the full ZNG model, where all collision processes have been included, unless it is stated otherwise.

Example 1: Dipole oscillation / Kohn mode
A very stringent test of any finite-temperature model of an interacting Bose gas is the dipole mode, or more commonly known as the Kohn mode [93,94]. For a harmonically-trapped gas, the center-of-mass (COM) degree of freedom is decoupled from all other internal degrees of freedom, hence the COM of a Bose gas would exhibit an undamped dipole oscillation at the trap frequency, independent of the interactions among the particles or the temperature of the gas. In the presence of a thermal cloud at finite temperature, the Kohn mode manifests as the undamped in-phase oscillation of the condensate and the thermal cloud [17], while the out-of-phase oscillation displays damping and frequency shift [17,95,96]. The Kohn mode thus serves as a very accurate way to measure the trap frequency [97].

Kohn Mode for a Single Condensate
In order to correctly reproduce the undamped Kohn mode for a single-component condensate, it is crucial to include the dynamics of both the condensate and the thermal cloud, and to couple them self-consistently. To see this, we can write down the coupled equations of motion of the COM along the z-axis as where define the COMs of the condensate (23a) and the thermal cloud (23b) and N c = dr n c (r) and N t = drñ(r) (24) give the number of condensate and thermal atoms respectively. F c,t is the force acting on the condensate due to the thermal cloud, and vice versa for F t,c . The COM of the whole atomic cloud, including both the condensate and the thermal atoms, is stated in terms of and its equation of motion is obtained by simply adding up (22a) and (22b). Neglecting the thermal cloud dynamics [i.e. omitting (22b)] or not coupling the condensate and the thermal cloud self-consistently (e.g. F c,t = −F t,c ) would therefore not lead to the correct equation of motion for z tot , In figure 1, we show that the ZNG model can correctly capture the Kohn mode described by equation (26) [98]. We start with 10 5 87 Rb atoms at thermal equilibrium (temperature T = 250 nK, 60% condensate) in an isotropic trap with angular trap frequency ω = 2π × 200 Hz. The scattering length is chosen to be a Rb = 99a 0 . At the beginning of our dynamical simulation, we rigidly shift the condensate and the thermal cloud by one harmonic trap length (ℓ Rb = 0.76µm) but in opposite directions along the z-axis. As the clouds oscillate out of phase, they exert a force on each other, causing damping in their oscillations (figure 1a) [17,95,96]. However, the total COM, shown in figure 1b, displays an undamped oscillation that gives no hint of the internal dynamics. Moreover, the oscillation frequency corresponds precisely to the trap frequency.

Kohn Mode for a Binary Mixture
For a two-component mixture, the dipole oscillation is very useful in measuring the miscible-immiscible transition [77].
where F ij is the sum of all forces acting on the ith component due to the jth component. A self-consistent model necessarily imposes the restriction F ij = −F ji , hence the weighted COM, would oscillate at the trap frequency if ω 1 = ω 2 . We demonstrate in figure 2 that the ZNG model can capture the Kohn mode for a binary mixture. We have 87 Rb and 41 K with 10 5 atoms each in isotropic traps with frequency ω Rb = ω K = ω = 2π × 200 Hz. The scattering lengths are a Rb = 99a 0 , a K = 60a 0 and a Rb−K = 20a 0 [99,100]. The temperatures of both components are 250 nK, which lead to 57% and 63% condensate fractions for 87 Rb and 41 K respectively. Similar to the single-component study, we shift the condensates and the thermal clouds by one harmonic trap length of 87 Rb (ℓ Rb = 0.76µm) in opposite direction at the beginning of our dynamical simulation. The subsequent dynamics of the condensates (solid lines in figure 2a), the thermal clouds (dashed lines in figure 2b) as well as the individual components as a whole (figure 2b), are distinctly different from those of a simple harmonic oscillator. In constrast, the weighted COM (figure 2c) clearly behaves like a simple harmonic oscillator oscillating at the trap frequency, proving that the ZNG model can reproduce the Kohn mode for a two-component mixture.
To the best of our knowledge, this is the only numerically-viable approach which accurately models the Kohn mode at finite temperature. Classical field (e.g. 'PGPE') [101], or stochastic (e.g. SGPE / SPGPE) methods by construction violate this, due to the lack of a dynamical handling of the abovecut-off atoms. Other common methods for manybody quantum systems, e.g. exact diagonalization approach [102] or multi-configuration time-dependent Hartree (MCTDH) method [103,104], are to date limited to investigations at zero-temperature or studies of equilibrium properties, whereas positive-P method [105,106] often has stability issues at long simulation times. It remains to be seen how important the violation of the Kohn mode actually is in practice in terms of other observables.

Example 2: Thermalization / Condensate Growth
We next explore the interesting problem of thermalization, with direct relevance to the condensate growth from a partially-condensed initial state. For an out-of-equilibrium Bose gas, the interactions between the particles tend to redistribute energy among them. Except for the very few cases (e.g. the integrable one-dimensional systems [58,107] or the Boltzmann monopole mode [108,109]), the non-equilibrium

Single-component thermalization
We first consider the evaporative cooling of a singlecomponent Bose gas [51]. At the final stage of the cooling process, a radio-frequency sweep is applied to quickly remove atoms with an energy above a certain cutoff energy. A non-equilibrium situation is thus created, where the atoms will interact and rethermalize in a completely-isolated environment, with a growth in the condensate number solely due to the internal dynamics of the Bose gas. A ZNG-type simulation of the above scenario has been carried out by Bijlsma, Zaremba and Stoof [54] using the ergodic approximation for the thermal cloud and the Thomas-Fermi approximation for the condensate (see also the quantum kinetic treatment of Davis et al [53]). We report here results from full ZNG simulations without relying on these approximations. In particular, we show that (i) our numerical implementation can establish a thermal equilibrium state from a highly nonequilibrium situation, (ii) there is a smaller increase in the condensate number if the cutoff energy is too low and (iii) there is a two-stage dynamics that have not been discussed before in the past studies [51][52][53][54].  Our simulation begins with 10 6 87 Rb atoms in an isotropic harmonic trap with angular trap frequency ω = 2π × 200 Hz at a temperature of 700 nK (critical temperature T 0 c = 900 nK for a non-interacting gas). The scattering length is chosen to be a Rb = 99a 0 . A self-consistent calculation with the Hartree-Fock approximation yields a condensate fraction of approximately 23%. We first generate approximately 7.7×10 6 test particles, which are distributed according to the Bose-Einstein distribution, to simulate the thermal cloud. Subsequently, any test particle with energy ε(r, v) [= 1 2 m Rb v 2 + U Rb n (r), where v is the velocity of the test particle, r is the position of the test particle and m Rb is the mass of a 87 Rb atom] above a cutoff value E cut = k B T cut is removed from our simulation. This truncation process represents the rapid quench in experiments. It leaves approximately 6.7 × 10 6 test particles in the trap if T cut = 3200 nK, but only approximately 2.4×10 6 test particles if T cut = 1200 nK. Equation (6) and (8) Figure 5. Effective temperature T eff of the thermal cloud (left axis, green diamonds), equilibrated temperature Teq (left axis, blue circles) and BEC fractions (right axis, red squares) for different energy cutoff k B Tcut. Same parameters as figure 4. The effective temperature measures the width of the thermal cloud in velocity space. The equilibrated temperature is determined as the temperature that will give the same condensate number and total atom number from an equilibrium calculation as the dynamical simulation.

BEC fraction
thermal cloud distribution. An example of such an initial thermal cloud distribution as a function of the speed |v| for T cut = 2400 nK is shown as the thin blue line in figure 3. Instead of a sharp jump in the distribution (see e.g. figure 3 of [53]), our distribution increases gradually once the speed falls below a cutoff value. This apparent difference stems from the fact that we are showing the truncated distribution as a function of the speed (hence integrating out the spatial dependence of energy ε(r, v) on the effective potential U Rb n (r)), while [53] showed the distribution as a function of the energy. For this same reason, while T cut = 2400 nK corresponds to a maximum speed of 21.4 mm/s, the true maximum speed in the simulation is slightly below 20 mm/s. The removal of the high-energy particles decreases the total energy of the system. The thermal cloud then rethermalizes to produce a growth of the condensate number (figure 4) through the C 12 collisional process, and repopulates the high-energy modes (i.e. the exponentially-decreasing tail of the thick green line in figure 3) through the C 22 collisional process. It is interesting to note that, while the truncation process always leads to an increase in the condensate number, the magnitude of increase at the end of the simulation saturates at around 2000 nK as the cut becomes deeper. Eventually, this trend is reversed, and the condensate starts decreasing as the cut is deeper than 1800 nK (figure 4a), in qualitative agreement with experiments (see, e.g. [1]). Nevertheless, the condensate fraction, computed as the ratio of the condensate number over the total atom number at the end of simulation (ωt = 30), increases monotonically as we consider a deeper cut ( figure 4b and figure 5), similar to our earlier findings with surface evaporative cooling [60].
In contrast to previous studies on condensate growth, which typically describe the rethermalization dynamics by a single thermalization time scale τ (with possibly also an onset time τ onset if the initial condensate fraction is negligible), our simulation data displays two thermalization time scales. This is revealed by a double exponential fit to the condensate number N c (t) as a function of time, The first time scale, τ 1 , is associated with the rapid growth in the condensate number. For our current parameters, this is approximately given by ωτ 1 ≈ 3 (vertical dashes in figure 4b). The second time scale τ 2 is an order of magnitude larger and its origin is not yet fully resolved. Similar numerical condensate growth curves that could possibly be fitted with the double exponential form (29) are also found in earlier works that employ the classical field approach [57] or the ZNG approach [26,54], even though the double exponential fit is not explicitly mentioned in any of those studies. In figure 6, we show the excellent agreement between the double exponential fit (blue dashes) and the simulated data (magenta solid line) for T cut = 2400 nK, with fitted time scales ωτ 1 = 4 and ωτ 2 = 27. Similar agreement between the numerical data and the fitted curves are also obtained for the other values of T cut . In Appendix A, we further compare the residuals using different fitting functions to convincingly demonstrate that the double exponential fit (29) is indeed the best choice. We note that the two-stage dynamics reported here is quite different from the various two-stage dynamics that have been studied before. Most of the previous works deal with the early stage of condensate growth where the condensate is negligible.
For examples, Kagan et al [50,55] consider the formation of a quasi-condensate with phase-fluctuations that die out to produce the truly phase-coherent condensate; Bose-stimulated process is considered to arrive at a relaxation process with an onset time [51]; the two-stage dynamics observed in [110] is attributed to the slow ergodic mixing at the very early stage of condensate growth; or a two-stage condensation dynamics due to the different transverse and axial trapping energies [111,112].
In comparison, our simulations start off with an appreciable condensate fraction (23%) and is accompanied by the emergence of a monopole oscillation at the later stage.
This appearance of the monopole oscillation is reminiscent of the quadrupole mode observed in the growth of an elongated condensate [113][114][115]. The top panel of the inset of figure 6 shows the effective temperature T eff of the thermal cloud, as a function of the simulation time t, where T eff measures the average kinetic energy of the N tp test particles.
The bottom panel of the inset shows the widths of the condensate and the thermal cloud normalized to the initial values, defined as σ 2 (t) = dr r 2 n(r, t), (31) and n(r, t) = n c (r, t) orñ(r, t). The oscillatory structure with an angular frequency of approximately 2ω appearing after τ 1 hints at the presence of a monopole mode. Further evidence of the monopole oscillation is given in Appendix B, where we extract the monopole and quadrupole oscillations from a double exponential fit of the condensate widths as well as the thermal cloud widths, but the extracted data are not able to provide conclusive evidence with regards to the potential relationship between the second time scale τ 2 and the monopole oscillation. In view of the recent results on the Boltzmann monopole mode [88,108,109], whether this two-stage dynamics remains valid for an anisotropic trap would be an interesting subject for future studies. At this point, it is worth noting that the increase in the effective temperature T eff of the thermal cloud (inset of figure 6) might not be apparent at first sight. This is because we typically associate cooling with a decrease in the temperature. The key idea is that it only makes sense to use T eff , which measures the width of the thermal cloud in the velocity space, to characterize the degree of coldness when the system reaches equilibrium. In fact, this increase of T eff in time occurs naturally as the speed distribution of the thermal atoms relaxes from a truncated distribution (thin blue line in figure 3), which has a smaller width because of the lack of high-energy population, to an equilibrium distribution (thick green line in figure 3), which has an exponential tail that extends far into the high-energy region. A plot of T eff at the end of our simulation versus the cutoff energy (green diamonds in figure 5) shows that a deeper cut consistently yields a cooler equilibrium system.
On the other hand, there is also an interesting increase in the spatial width of the thermal cloud ( figure 6). This is because thermal atoms further away from the trap center are more likely to be removed by the truncation process, hence a pressure difference is created to push the thermal cloud outwards once the dynamical simulation commences. The outward expansion of the thermal cloud in turn can drive the monopole mode oscillation of the condensate through the mean field potential. An immediate consequence of this picture is that, the deeper is the cut (i.e. smaller T cut ), the larger is the pressure difference, hence the larger is the oscillation amplitude (see also Appendix B).
In spite of the long relaxation time τ 2 , we find that it is still possible to assign an equilibrium temperature T eq to the system at the end of the equilibration process (ωt = 30) §. This is done by searching for a temperature that yields the same condensate number and total atom number in our equilibrium calculations. The condensate density and the thermal cloud density obtained from our equilibrium calculation (dashes and dash-dots) are plotted against the density profiles at ωt = 30 of the dynamical simulation (solid lines) in figure 7 and show remarkable agreement. The decrease in T eq for deeper cut (figure 5) is consistent with our intuition and reaffirms the validity of the ZNG model.

Two-component thermalization
We now consider the sympathetic cooling of a twocomponent mixture [116]. In a typical experimental situation, the first component can be easily cooled, e.g. § We expect the system to reach its true equilibrium after about 10τ 2 . However, numerically simulating up to this time would require a computational time of one to two months. Since the condensate number is expected to increase by at most ten percent if we perform simulation from ωt = 30 to ωt = 300, we believe that our simulations running up to ωt = 30 is sufficient for the present purpose. via the application of a radio-frequency sweep, while the second component is cooled by being in thermal contact with the first component. Atoms of different components can therefore collide elastically and exchange energy as long as the two components overlap in space. The thermalization rate of such a scenario has been investigated both experimentally [117] and theoretically [117][118][119][120][121]. However, these theoretical estimates often rely on approximations to simplify the calculations, such as (i) the high-temperature approximation where all atoms follow a Maxwell-Boltzmann distribution, or (ii) the omission of mean-field contribution, where the condensate wavefunction is identified as the ground state of a harmonic oscillator.
A thorough investigation that takes into account the presence of both condensates self-consistently is still missing. We aim to eventually fill this gap by performing systematic studies of sympathetic cooling using our two-component ZNG model.
Due to the many different collisional processes involved (eight collisional 9000 10000 BEC number processes in a two-component mixture [75,76] versus two collisional processes in a single-component Bose gas), the dynamics is much more complicated and also much more interesting to study.
In figure 8, we show that our two-component ZNG model can establish an equilibrium solution in a sympathetic cooling setup. We start with a 87 Rb-41 K mixture, each with 20,000 atoms in an isotropic harmonic trap with identical trap frequency ω = 2π × 200 Hz. Initially, there is no thermal contact between the two components (i.e. inter-component scattering length a Rb−K = 0). The intra-component scattering lengths are a Rb = 99a 0 and a K = 60a 0 . The two components have different initial temperatures, T Rb = 160 nK and T K = 180 nK, and different condensate fractions (51% for 87 Rb and 43% for 41 K). The critical temperatures are T c,Rb = 226 nK and T c,K = 233 nK when taking into account of the finite-size corrections and mean-field corrections [97].
In order to initiate the cooling process, we linearly ramp up the inter-component scattering length a Rb−K from 0 to 20a 0 in 10 ms, and maintain the value afterwards. In this case, the mixture remains miscible [g 2 Rb−K /(g Rb g K ) = 0.3] such that there is a large spatial Rb, BEC K, BEC Rb, thermal K, thermal overlap of the two components (see figure 9 for the density profiles at t = 500/ω ≈ 400 ms). For realistic values of scattering lengths, the thermalization time scales are of the order of seconds, which translate into a very long computational time . In order to speed up the equilibration process within our simulations, we increases the inter-component thermalthermal collisional cross-section by 50 times, and the inter-component thermal-condensate collisional crosssection by 10 times when sampling the collision events.
In practice, this means that we make the substitutions in equation (6) and in equation (8) for k = j, κ 22 = 50 and κ 12 = 10.
In addition, we omit the condensate-exchange process (C kj 12 and R kj 12 ) in our simulation as its time scale is typically an order of magnitude shorter than the other collisional processes [75,76]. atoms, as 87 Rb atoms sympathetically cool 41 K atoms. Figure 8(b) shows the corresponding change in the effective temperatures of the thermal clouds (30). The fact that the condensate numbers saturate to finite values and the two effective temperatures converge shows that an equilibrium situation is established at the end of our simulation (ωt = 500).
Both condensate numbers in figure 8(a) are well fitted by an exponentially-decaying function, with the thermalization time scales ωτ Rb = 83 and ωτ K = 76. These time scales are to be compared with the prediction based on Maxwell-Boltzmann distribution [117], where σ Rb−K = 4πa 2 Rb−K is the cross-section, the factor κ 22 takes into account that we have scaled up our collisional probabilities to speed up the computation, . This yields the estimate ωτ ≈ 40 that is comparable to the simulated time scales. This agreement could stem from the fact that we have scaled up the inter-component thermal-thermal collisional probabilities, making all collisional processes to have comparable time scales within our simulations.
It is interesting to see if such an agreement remains when we restore the true collisional probabilities. In figure 10, we compare simulations with scaled probabilities (solid lines, κ 22 = 50, κ 12 = 10) and those with the true collisional probablities (κ 22 = 1, κ 12 = 1) that includes (dashed lines) or omits the condensate-exchange C 12 collisions (dotted lines). Our numerical results reveal that the condensate-exchange collision leads to faster thermalization at the initial stage, but it has less impact at the later stage, where the condensate numbers vary at comparable time scales in the presence or absence of the condensate-exchange collisions. In order to deduce the true thermalization time, the condensate numbers are plotted in terms of the rescaled timet, wheret = t for the simulation with scaled probabilities, butt = t/10 for the simulations with true probabilities. The similarity of the curves implies that the true thermalization time scale is likely to be of the order of ωτ = 800, shorter than the ωτ = 2000 predicted by equation (35).
As we use realistic numbers in our simulation, we expect the two-component thermalization time in experiments to be of the order 800/ω ≈ 0.6s, which is in fact comparable to typical waiting times in experiments. Taking into account that the overlap of two components is reduced due to the relative shift of trap minima, our simulation results therefore suggest that the two components may not necessarily have fully thermalized to a common temperature in experiments. This interesting question clearly merits further detailed investigation, both from a theoretical and an experimental perspective. An important factor to consider here is that on one hand one seeks a large inter-component interaction strength g 12 in order to speed up the C kj 12 (k = j) collisions that dominate the process of sympathetic cooling; however, a larger g 12 (> 0) leads to enhanced phase separation, thus minimizing the effective overlap over which thermalization / cooling takes place, resulting in a slower sympathetic cooling rate. A critical balance between those two competing mechanisms is here needed to optimize the crossthermalization and the sympathetic cooling process efficiently.

Outlook
We have demonstrated that our numerical implementation of the ZNG model can capture the essential physics of non-equilibrium Bose-Einstein condensates at finite temperature, namely the Kohn mode and the rethermalization dynamics associated with condensate growth for both single and binary atomic gases.
Revisiting the well-studied single-component problem, but for the specific case of a spherically-symmetric (isotropic) trap, we have observed two interesting features which deserve further attention, also from an experimental point of view. Firstly, we found that, following a rapid evaporative cooling quench, the system grows to a state with higher condensate fraction on two distinct timescales, a rapid one associated with rapid condensate number growth, and a much slower one. Analogous, yet physically distinct, two-stage condensate formation dynamics have been previously reported in various different contexts, e.g. associated with the excitation of a quadrupolar mode for highlyelongated systems. Parallel to this, we have found that the condensate growth process in an isotropic trap naturally excites the monopole mode, which in this geometry is long-lived and so could be experimentally observable (whereas in our case, excitation of the quadrupolar mode is suppressed). It would thus be very interesting to see an experimental study of controlled condensate growth, based on rapid evaporative cooling, as a function of trap aspect ratio. Extrapolating beyond our findings, one could envisage a situation whereby the condensation process leads either to a quadrupolar mode excitation for very elongated traps (for which monopole excitation is suppressed), or to a monopolar excitation for isotropic traps, with intermediate geometries having a variable amount of excitation and decay timescales for the different excitation modes. It would also be interesting to study whether the excited modes are in-phase or out-of-phase for the condensate and thermal cloud, a problem significantly more complicated than the previously conducted controlled excitation experiments [15,88,109], as the condensate fraction and number are constantly changing during the condensate growth. The possiblity of exciting different modes following shock-cooling based on system geometry thus appears to be well-worth analyzing further.
Interesting open problems can also be explored using the two-component ZNG model, with the model predictions benchmarked against experimental results. Two particularly attractive research directions concern the controlled studies of the collective modes of a mixture and the sympathetic cooling in the presence of a partially-condensed Bose gas.
The former has gathered increasing experimental interest [122][123][124][125] recently, but a major experimental challenge remains in reducing the relative shift of the trapping-potential minima. This relative shift can arise if the two components experience different spring constants of the harmonic potentials, or if the two components have different masses [126]. Any experimental achievement in minimizing/eliminating the trap sag would immediately open up the possibility to study the effect of inter-component interaction on collective modes, such as the number-dependent miscible-immiscible transition [77] or the bifurcation of single-component collective modes [127]. In addition, the presence of the thermal clouds, which could possibly reach the hydrodynamic regime, leads to a 'four-fluid' model of the binary mixture with more complicated dynamics, a glimpse of which can be found in [75,76].
On the other hand, while sympathetic cooling has become a routine procedure in cold-atom experiments, the theoretical investigation of the sympathetic-cooling rates below the critical temperature is still quite limited. Interesting questions that revolve around the thermalization time scales include the role of the condensates [119] and the effect of the spatial overlap between the two components (which is in turn affected by the trap sag or the miscibility). The eight collisional processes with disparate collisional time scales [75,76] also lead to the fascinating possibility of quasiequilibrium dynamics [63]. How the different collisional processes can enhance or suppress each other, and how this problem can be explored experimentally, are open questions to be answered.
Finally, we would like to highlight that our atomistic simulation of the two-component ZNG model are not limited to study a Bose-Bose mixture. With the appropriate modifications to the quantum Boltzmann equations (8) and the mean-field potentials (9a) and (9b), we can also study the non-equilibrium dynamics of a Bose-Fermi mixture, including the case where both the bosons and the fermions exist in the superfluid phase [122], or a partially-condensed Bose gas immersed in a spin-polarized Fermi sea [128,129].  The condensate growth curves of a singlecomponent Bose gas are fitted with three different functions:

Acknowledgments
(i) A two-stage formula due to Bose-stimulated process [51], allowing for an initiation time,  Since the residuals of the double exponential functions (red dash-dots) always remain small compared to the other two functions, we conclude that our simulated data is indeed well described by a double exponential function.

Appendix B. Extraction of monopole and quadrupole modes
We extract the monopole and quadrupole oscillation from the widths of the condensate as well as the thermal clouds. The background σ 2 bg (t) is first determined by fitting the radial width σ 2 (t)/3 (31) by the double exponential function (A.3), but with the time scales τ 1 and τ 2 as input parameters, where their values are determined by the fit of condensate numbers.
The oscillation amplitude is larger for the deeper cut because of the greater pressure that pushes the thermal atoms outwards at the start of the simulation (see the discussion near the end of section 5.1). We further determine the dominant oscillation frequencies by In general, the monopole mode displays oscillation with angular frequency between 2ω and 2.5ω, in reasonable agreement with the expected monopole frequency of a pure thermal cloud (2ω) [130] or a pure condensate ( √ 5ω) [131]. The oscillation amplitude decays in time, but it is difficult to extract a decay time scale due to the lack of regular decaying pattern, associated with the underlying condensate growth.
Because of numerical fluctuations generated in our simulations, we can also see small amplitude quadrupole oscillation that displays more irregularity than the monopole oscillation. The fast Fouriertransformed fluctuations display peak between ω and 1.5ω, in reasonable agreement with the quadrupole mode frequency √ 2ω of a pure condensate [131]. Due to the complicated nature of the excitation mechanism of such processes, occurring on top of a constantly cooling sample undergoing condensate growth, we are unable at the present time to extract more information related to the induced in-phase and out-of-phase oscillations, which have proven crucial in interpreting experimental findings under controlled excitation schemes [15,16,19,21,22,28,29]. A more detailed study of the coupled excitation amplitudes under different cooling protocols and geometries remains an interesting question for further studies.