Time-dependent condensation of bosonic potassium

We calculate the time-dependent formation of Bose--Einstein condensates (BECs) in potassium vapours based on a previously derived exactly solvable nonlinear boson diffusion equation (NBDE). Thermalization following a sudden energy quench from an initial temperature $T_\mathrm{i}$ to a final temperature $T_\mathrm{f}$ below the critical value and BEC formation are accounted for using closed-form analytical solutions of the NBDE. The time-dependent condensate fraction is compared with available $^{39}$K data for various scattering lengths.


Introduction
The time-dependent mathematical modelling of nonequilibrium phenomena is an outstanding problem in physics: How does an isolated quantum many-body system thermalize that is originally far from equilibrium? This problem is particularly relevant and interesting for ultracold Bose gases, where in the course of thermalization a transition to a Bose-Einstein condensate (BEC) can occur below a critical temperature T c which it is not yet fully understood [1].
Ever since and even before the discovery of BEC in alkali atoms, numerical nonequilibrium theories have been developed to account for the time-dependence of thermalization and condensate formation. Predictions such as [2] for sodium and rubidium had been made. For 23 Na, some of these theoretical results for the timedependent condensate growth [3] have been compared with experimental data [4].
For 87 Rb, data have been compared to a numerical model based on quantum kinetic theory [5]. Here, cooling into the quantum-degenerate regime was achieved by continuous evaporation for a duration of several seconds, whereas in the sodium experiment a short radiofrequency pulse was applied to remove high-energy atoms.
Only recently, new experimental results for potassium have become available with a sufficiently detailed time resolution: Thermalization and condensate formation is investigated in a homogeneous 3D Bose gas of 39 K with tunable interactions and near-perfect isolation in a cylindrical optical box [1].
As a transparent and analytically solvable model for the time-dependence of thermalization and condensate Institut für Theoretische Physik der Universität Heidelberg, Philosophenweg 16, D-69120 Heidelberg, Germany, EU formation of bosonic systems, a nonlinear boson diffusion equation (NBDE) has been derived in Ref. [6], applied to ultracold atoms [7,8], and solved exactly with the necessary boundary conditions at the singularity [9]. We have used the model for the calculation of the time-dependent condensate fraction in 23 Na [10], where we have tested it successfully against the historical MIT data [4]. In the present work, we perform a related calculation for 39 K to compare with the newly available Cambridge data for tunable scattering lengths [1].

The model
First we briefly review our previously developed [6,7,8,9,10] nonlinear model that we use in this work to compute the nonequilibrium evolution following a sudden energy quench, and time-dependent condensate fractions in equilibrating Bose gases of ultracold atoms with a focus on a comparison with recent data [1] for 39 K. The cold-atom vapour in a trap is described as a timedependent mean field with a collision term and we consider only s-wave scattering with the corresponding scattering lengths. The N -body density operatorρ N (t) is composed of N single-particle wave functions of the atoms which are solutions of the time-dependent Hartree-Fock equations plus a time-irreversible collision term K N (t) that causes the system to thermalize through random two-body collisions ( = c = 1) with the self-consistent Hartree-Fock mean fieldĤ HF (t) of the atoms in a trap that provides an external potential.
The full many-body problem is reduced to the onebody level in an approximate version for the ensembleaveraged single-particle density operatorρ 1 (t). Its diagonal elements represent the probability for a particle to be in a state |α with energy ǫ α ρ 1 (t) α,α = n(ǫ α , t) ≡ n α (ǫ, t) . ( The total number of particles is N = α n α , and we neglect here the off-diagonal terms of the density matrix. As discussed in Refs. [6,8,10], the occupation-number distribution n α (ǫ, t) in a finite Bose system obeys a Boltzmann-like collision term, where the distribution function depends only on energy and time. A condition for such a reduction to 1 + 1 dimensions is spatial and momentum isotropy. The corresponding assumption of sufficient ergodicity has been widely discussed in the literature [11,12,13,14,15,16]. For the thermal cloud of cold atoms surrounding a Bose-Einstein condensate (BEC), it is expected to be essentially fulfilled, even though the condensate in a trap is spatially anisotropic. Different spatial dimensions enter our present formulation through the density of states, which depends also on the type of confinement. The model calculations in this work are for a 3D system, and in line with the 39 K experiment [1], we investigate results for the density of states of bosonic atoms confined in a cylindrical optical box. The collision term for the single-particle occupationnumber distribution of the energy eigenstates ǫ α , n ≡ n α ≡ n(ǫ α , t) , can be transformed [6] into a nonlinear partial differential equation In this nonlinear boson diffusion equation (NBDE), the drift term v(ǫ, t) accounts for dissipative effects, whereas D(ǫ, t) mediates diffusion of particles in the energy space: It accounts for the broadening of the distribution functions, and in particular, for the softening of the sharp cut at ǫ = ǫ i that signifies the quench, as well as for the diffusion of particles into the condensed state. The many-body physics is contained in these transport coefficients, which depend on energy, time, and the second moment of the interaction. It has been shown [10] that the stationary solution n ∞ (ǫ) of the NBDE for t → ∞ equals the Bose-Einstein equilibrium distribution n eq (ǫ) provided the ratio v/D has no energy dependence, requiring lim t→∞ [−v(ǫ, t)/D(ǫ, t)] ≡ 1/T . The chemical potential is µ ≤ 0 in a finite Bose system.
In the limit of energy-independent transport coefficients, the nonlinear boson diffusion equation for the occupation-number distribution n(ǫ, t) becomes Again, the thermal equilibrium distribution n eq is a stationary solution with µ ≤ 0 and T = −D/v. The NBDE with constant transport coefficients preserves the essential features of Bose-Einstein statistics that are contained in the bosonic Boltzmann equation. It is one of the few nonlinear partial differential equations with a clear physical meaning that can be solved exactly through a nonlinear transformation, as outlined in Refs. [6,7,8]. The resulting solution is where the time-dependent partition function Z(ǫ, t) obeys a linear diffusion equation which has a Gaussian Green's function G(ǫ, x, t). The partition function can then be written as an integral over the Green's function and an exponential function F (x) F (x) depends on the initial occupation-number distribution n i according to The definite integral over the initial conditions taken at the lower limit in Eq. (8) drops out in the calculation of n(ǫ, t) and can be replaced [9] by the indefinite integral with the integration constant set to zero without affecting the accuracy of the calculation. These modifications allow us to compute the partition function and the overall solution for the occupation-number distribution function Eq. (6) analytically.
With the free Green's function G ≡ G free of the linear diffusion equation, the physically correct solution with the Bose-Einstein equilibrium limit is attained in the UV region, but not in the IR [6]. To solve this problem, one has to consider the boundary conditions at the singularity ǫ = µ ≤ 0 [8]. They can be expressed as lim ǫ↓µ n(ǫ, t) = ∞ ∀ t. One obtains a vanishing partition function at the boundary Z(µ, t) = 0, and the energy range is restricted to ǫ ≥ µ. This requires a new Green's function [8] that equals zero at ǫ = µ ∀ t. It can be written as (9) and the partition function with this boundary condition becomes which is equivalent to The function F remains unaltered with respect to Eq. (8), but its argument is shifted by the chemical potential.
For fixed chemical potential µ and an initial temperature T i that differs from the final temperature T f = −D/v as in an energy quench, we have solved the combined initial-and boundary value problem exactly in Ref. [9] using the above nonlinear transformation [8] from Eq. (6). Expressions of the form [1−exp (−z/T i )] Ti/T f that appear in F (x) and thus, in the partition function Eq. (10), can be reformulated using the generalized binomial theorem which results in an infinite series expansion for the timedependent partition function (and also for its derivative), The analytical expressions for f Ti,T f k (ǫ, t) are combinations of exponentials and error functions, they are given explicitly in Ref. [9]. We can now compute the distribution function directly from Eq. (6). The convergence of the solutions has been studied in Ref. [10], indicating that -depending on the system and the specific observable -a result with k max = 10 − 40 expansion coefficients is indistinguishable from the exact solution with k max = ∞. It has also been shown that the series terminates in case T i /T f is an integer, and we shall explicitly use that property in this work, where the exact solution is already obtained for k max = 4. A brief summary of the nonlinear diffusion model for both, boson and fermion systems at low and high energies and temperatures is presented in Ref. [17].

Energy quench and thermalization
In the 39 K experiment of Ref. [1] with tunable interactions, a far-from-equilibrium atomic cloud is created by removing 77% of the atoms and ≃ 97.5% of the total energy through an energy quench below the critical temperature, thus setting the stage for subsequent condensate formation. This quench corresponds to cutting off the atoms with energy ǫ > ǫ i according to  Fig. 1, but using µ = 0 and a double-log-scale. Thermalization in the UV tail is significantly slower than below the evaporative cut and in particular, in the IR. Timesteps are t = 8, 30, 100, 300, and 4000 ms. The dot-dashed line is the Rayleigh-Jeans power law in the initial distribution.
Hence, the number of atoms just after the quench becomes with the 3D density of states In agreement with experiment [1], it is kept constant in the subsequent thermalization and condensate formation process. The equilibrium value of the ensuing time-dependent condensate fraction for t → ∞ can then be calculated from the ratio of the difference in particle content of the initial nonequilibrium distribution at finite µ i < 0 minus the final equilibrium Bose-Einstein distribution with temperature T f and vanishing chemical potential, divided by particle content of the initial distribution. The result is where the integration over the initial thermal distribution has an upper limit of ǫ = ǫ i , according to the cut due to the quench, θ (1 − ǫ/ǫ i ). The measured equilibrium condensate fraction [1] is N eq c /N i = 40(5) % for all scattering lengths that have been investigated, a = (100 − 800) a ∞ with the Bohr radius a ∞ .
The above Eqs. (14), (15) provide two conditions to calculate the initial value of the chemical potential µ i , and the cut ǫ = ǫ i that defines the quench. Inserting Eq. (14) into Eq. (15), an implicit equation for the initial chemical potential µ = µ i that does not depend on ǫ i is obtained, and with this value, the cut ǫ = ǫ i can directly be computed from Eq. (14).
In the 39 K-experiment, the initial temperature before the quench is T i = 130 nK, and the final temperature T f = 32(2) nK in accordance with energy-and particle-number conservation [1]. In view of the fact that these numbers as well as the experimental equilibrium condensate fraction of 40(5) % are approximate, we use a value of T f ≡ 32.5 nK in our model calculation. This considerably simplifies the analytical calculations, because the value of T i /T f is an integer, such that the series expansion of the exact analytical solution of the NBDE terminates already at k max = 4. (A calculation with T f = 32.0 nK and k max = 40 yields results which are hardly distinguishable).
With these values for the initial and final temperatures, together with a removal rate of 77% of the atoms in the quench and an equilibrium condensate fraction of 40%, we obtain the resulting values for the initial chemical potential and the cooling cut 1 as At and below the cut, the average initial occupation is n i (ǫ ≤ ǫ i ) ≥ 7.52. The nonequilibrium system is highly overoccupied and will move quickly into the condensed state to remove the surplus particles from the cloud. In the simplified model with constant transport coefficients, the values of v and D that are required for the nonequilibrium calculation are related to the equilibrium temperature T f = −D/v and the local thermalization time τ eq . The first relation is a consequence 1 removing 97.5% of the energy [1] requires a slightly larger value of the cut of the equality of the stationary solution of the nonlinear diffusion equation with the Bose-Einstein distribution as discussed above, whereas the thermalization time has been derived from an asymptotic expansion of the error functions in the analytical NBDE solutions [6] for θ-function initial conditions at the cut ǫ = ǫ i as τ eq ≡ τ ∞ (ǫ = ǫ i ) = 4D/(9v 2 ). For a general initial distribution with a cut at ǫ = ǫ i , it can be written as with f = 4 for fermions [6] and f = 4/9 for bosons in case of initial theta-function distributions. The analytical result for a truncated Bose-Einstein initial distribution has not been derived yet. It will be considerably shorter than the one for a theta-function distribution, because the shape of the initial distribution is much closer to the final BE-result than a box distribution.
Here we take f = 0.045 to compute the transport coefficients in order to be consistent with the thermalization time of τ exp eq ≃ 600 ms estimated from the data [1] for a = 140 a ∞ .
With a final temperature T f = −D/v = 32.5 nK and a condensate-formation time τ eq = f D/v 2 ≃ 600 ms for a scattering length of a = 140 a ∞ , the transport coefficients that are required for the analytical solution of the NBDE are thus obtained as v = −f T f /τ eq ≃ −0.00246 nK/ms .
The nonequilibrium calculation can now be performed with the above values for the drift and diffusion coefficients and fixed chemical potential. First we compute the distribution functions for quenching the distribution from T i to T f according to Eq. (6) with the timedependent partition function Z from Eq. (12). The results for thermalization in 39 K are shown in Figs. 1, 2. Analogous analytical solutions -albeit for different parameters -were shown in Fig. 2 of Ref. [10] to agree precisely with numerical Matlab-solutions of the NBDE using appropriate boundary conditions and hence, we can trust that the results are correct. Fig. 1 shows the time-dependent thermalization from the initial nonequilibrium distribution that is produced via quenching at ǫ i = 15.56 nK to the final Bose-Einstein equilibrium distribution with T f = 32.5 nK and fixed chemical potential. Using the above transport coefficients, results are shown for seven timesteps at t = 0.2, 4, 20, 60, 100, 200, and 400 ms with increasing dash lengths. Thermal equilibrium in the cloud is achieved quickly at t < 100 ms in the infrared, somewhat more time is needed at intermediate energies, and much more time in the ultraviolet. The equilibration in the UV thermal tail is more clearly visible in Fig. 2, which is a double-logarithmic plot with an additional timestep at 4000 ms. Thermalization in the cloud below the cut is seen to be essentially completed at the equilibration time τ eq , but it takes much longer in the far Maxwell-Boltzmann tail. The Rayleigh-Jeans power law in the initial distribution is clearly visible, and also power-law behaviour of the subsequent nonequilibrium evolution in limited energy intervals.

Time-dependent condensate formation
With the analytical NBDE solutions for a given chemical potential µ, we can now proceed to calculate the time-dependent transition to the condensate with the additional condition of particle-number conservation at each timestep. This requires a time-dependent chemical potential that eventually approaches zero, µ(t) → 0; from then on, particles start moving into the condensed state. The particle number in the thermal cloud at any time t is given by withñ(ǫ, t) the distribution function calculated using the value of µ ≡ µ(t) ≤ 0 that is computed at each timestep such that the particle number is conserved, Once µ(t) has reached the value zero, condensation starts, and and any further difference between the initial and the actual particle number is attributed to the condensate, The definition of a time-dependent chemical potential is debatable, as the distribution function for µ i ≤ µ(t) < 0 does not have the equilibrium form. However, for µ(t) < 0 we do not use the calculation to compute physical observables. It is only when equilibrium at µ(t) = 0 is reached and the phase transition to the condensate occurs that we start calculating the timedependent condensate fraction from the NBDE solution plus the condition of particle-number conservation. Although the transport equation cannot explicitly treat the buildup of coherence in the condensate because it does not consider phases, the nonequilibrium-statistical approach together with the condition of particle-number conservation thus allows to properly account for the time-dependent number of particles that move into the condensed state.
To evaluate the integral in Eq. (19) for the particle number in the thermal cloud once µ = 0 has been reached at t = τ ini , no analytical solution is presently known for finite time t < ∞, but it can be evaluated numerically. The resulting time evolution of condensate formation in 39 K is calculated first with the transport parameters D, v of table 1 for a scattering length The initial chemical potential is µ i = −0.67 nK and it increases with time towards zero, when condensation begins. We have tested [10] that no condensate forms if T f at the critical number density n c would remain above the critical temperature T c for bosonic atoms of mass m, This expression is a consequence of equilibrium-statistical mechanics -at the critical temperature, the chemical potential becomes zero. It implies that sufficient coherence for the phase transition to the condensate is achieved when the interparticle distance becomes equal to the thermal de Broglie wavelength. Since the Bose-Einstein equilibrium distribution is equal to the stationary limit of the NBDE for t → ∞, the critical temperature is an integral part of our nonequilibrium model, just as it is part of equilibrium statistics. Accordingly, the equilibration time τ eq is equivalent to the formation time of the thermalized condensate.
We can now compare the NBDE-solutions and the ensuing time-dependent (quasi-)condensate fractions with the data [1] for various scattering lengths, and here we focus on a/a ∞ = 140, 280, 400 and 800. With both transport coefficients D, v scaled linearly with the scattering lengths (table 1), the results are displayed in Fig. 3. The values of the transport coefficients have not been further optimzed with respect to the Cambridge data at the individual scattering lengths. A related log-logplot is shown in Ref. [17], where the nonlinear diffusion equation is discussed in a more general context of thermalization at low and high energies as well as temperatures for both, bosonic and fermionic systems.
Condensation starts at the initiation time τ ini , and rises according to the condensate-formation time τ eq (table 1) towards the equilibrium value N eq c /N i . The initiation time required from the data is larger than the calculated time to reach µ = 0 due to the assumption of constant transport coefficients, but considerably smaller than the condensate-formation time, τ ini ≪ τ eq . These calculations are performed with our C++ code that is based on the exact analytical NBDEsolutions for n(ǫ, t) from Eq. (6) with k max = 4 in the series expansion, and uses numerical integration methods to compute µ(t) and N c (t)/N i .
The resulting thermalization timescales are seen to be longest for the shortest scattering length a = 140 a ∞ , and decrease towards larger values of a/a ∞ = 280, 400, 800. As is evident in particular for a = 140 a ∞ , small deviations in the short-time behaviour may be due to our basic assumption of constant transport coefficients, which causes condensate formation to start with full strength and not gradually. This feature appears to be absent in the 87 Rb data [5], but it had also been observed when comparing with 23 Na data [10]. As an improvement, one would have to solve the full NBDE Eq. (3) numerically with time-and energy-dependent transport coefficients.
With constant transport coefficients, we confirm within the nonlinear diffusion model what has already been shown in the experimental work [1] based solely on their data: The time-dependent condensate fractions for different scattering lengths fall on a single universal curve when plotted as functions of t ′ = ta/(300 a ∞ ), before reaching the equilibrium limit for t ′ → ∞.
The 1/a-dependence of the timescales -including that of the initiation time τ ini , and the condensateformation time τ eq -implies that these are set by the inverse interaction energies [12,19], and not by the inverse cross sections [20], which would yield a 1/a 2 -scaling: A larger interaction energy causes faster initiation of the (quasi-)condensate formation, and also a more rapid thermalization towards the equilibrium value of the condensate fraction. The 1/a-rather than 1/a 2 -dependence is due to the emerging coherence between the highly occupied IR modes [21].
The nonequilibrium system shows self-similar scaling [1], with a net flow of particles towards the infrared (bottom-up) and into the condensed state. This is accompanied by a corresponding energy flow towards the ultraviolet carried by a relatively small number of atoms that build up the MB tail, see Figs. 1 and 2.

Conclusions
While statistical equilibrium for bosons and fermions is uniquely determined by temperature T and chemical potential µ in the Bose-Einstein and Fermi-Dirac distributions, we can account for the time-dependent approach to equilibrium via a nonlinear diffusion equation.
In the approximate case of constant diffusion and drift coefficients, we have solved the equation exactly through a suitable nonlinear transformation. The nonequilibrium system behaviour is determined through diffusion coefficient D and drift coefficient v, the initial and final temperatures, and the initial chemical potential. Particle-number conservation is ensured during the thermalization process via elastic two-body collisions. This is essential for the proper calculation of the timedependent condensate fraction.
We have applied this model to the thermalization of ultracold 39 K atoms with tunable interaction strength, and compared to recent Cambridge data. As in our earlier comparison with 23 Na MIT data for the timedependent condensate fraction, overall agreement with the experimental results is found. Minor deviations in the short-time behaviour at small scattering lengths might indicate the limits of a description that relies on constant transport coefficients in order to obtain an analytical solution of the problem.
The transport equation does not consider phases and -as other Boltzmann-type numerical approaches -does not explicitly treat the buildup of coherence in the condensate. Nevertheless, the nonequilibrium-statistical method together with the condition of particlenumber conservation allows to properly account for the time-dependent number of particles that move into the condensed state, because it includes the condition that condensate formation starts at the point in time where the chemical potential vanishes.
In agreement with the Cambridge data, the condensate-initiation time τ ini and the condensate-formation time τ eq are inversely proportional to the tunable scattering length, inducing shorter time constants for larger scattering lengths and a universal dependence of the (quasi-)condensate fraction on a time parameter that scales linearly with the scattering length. The condensateinitiation time is found to be always considerably shorter than the condensate-formation time.
It would be valuable to have precise experimental information about the time-dependent condensate formation in other cold-atom systems such as 87 Rb and 7 Li using a single-quench technique to remove the highenergy atoms under similarly controlled conditions in order to test the available numerical and analytical nonequilibrium-statistical approaches in detailed comparisons with data.