Universal dynamics on the way to thermalization

Suppose we start with a cold, dilute, incoherent homogeneous gas in three dimensions, at a temperature above the BEC phase transition. Let us apply a cooling quench to this gas by removing particles from the high-momentum tail of the distribution and suppose that, appropriate collisions provided, the system re-equilibrates to thermal equilibrium. The final temperature is then given by the total energy after the quench, which we assume to be below the BEC critical temperature. We note that the question of whether a system can thermalize at all [24] can be addressed by exloring which types of stable attractors exist.

Two different scenarios are possible: If a sufficiently small amount of energy is removed, the subsequent scattering of particles into the low-energy modes builds up a thermal Rayleigh−Jeans distribution, $n(k)\sim {{k}_{{\rm B}}}T/{{k}^{2}}$, in a quasi-adiabatic way. The chemical potential increases, and a fraction of particles is deposited in the lowest mode, forming a BEC. During this process, tangles of defect lines can be found in the Bose field by filtering out short-wavelength fluctuations [17].

In the second scenario, given a sufficiently strong cooling quench cutting away (e.g., all particles above a cutoff, k_c), the remaining overpopulation of the modes just below k_c results in a vigorous transport towards lower energies that has the form of a strong-wave-turbulence inverse cascade [25]. This cascade induces a long-lived, power-law single-particle spectrum, $n(k)\sim {{k}^{-\zeta }}$, with an exponent, $\zeta =5$, that is distinctly larger than the exponent $\zeta =2$ of the thermal Rayleigh–Jeans distribution. The emergence of the strong cascade can be explained by the dominance of vortical superfluid flow around line defects over compressible, longitudinal sound excitations and density fluctuations in the respective regime of wavelengths. The power law, ${{k}^{-5}}$, signals that the system evolves near a distinct non-thermal fixed point (NTFP) where the evolution critically slows down. This power law has been predicted by using a nonperturbative Greenʼs function as well as renormalization-group techniques [2, 26, 27]. This particular power can also be traced back to the flow pattern around the vortex lines [25], which now become visible without filtering out short-wavelength fluctuations. The two possible paths to BE condensation are shown schematically in figure 1. Whether the system, during the condensation dynamics, approaches the NTFP or moves in a direct way to thermal equilibrium depends, for a closed system, on the initial conditions (i.e., on the strength of the cooling quench). We refer to the condensation process which takes the 'detour' via the NTFP as hydrodynamic BE condensation because of the dynamical scale separation of incompressible and compressible components of the particle current.

To reveal these dynamics, we study a dilute Bose gas in the classical-wave limit, statistically sampling the classical equation of motion that has the form of the Gross–Pitaevskii equation (GPE) for the classical field, $\phi ({\bf x},t)$ (see appendix A.1 for details). We compute the evolution in a cubic box with periodic boundary conditions. The initial field in momentum space, $\phi ({\bf k},0)$, is sampled from a distribution that allows its phase to be random, while $|\phi ({\bf k},0){{|}^{2}}$ has Gaussian fluctuations around a momentum distribution, $n({\bf k},0)$, which is flat between $k\;=\;|{\bf k}|=0$ and k ≃ k₀, and falls off as ${{k}^{-\alpha }}$ for k ⪆ k₀ (top left panel of figure 3).

Zoom In Zoom Out Reset image size

We consider an initial power-law falloff of n(k) that is close to or steeper than that expected by a self-similar solution of the wave Boltzmann equation, $\alpha \simeq 2.4$ [10], corresponding to an inverse particle cascade in weak wave turbulence theory [6, 7]. $\alpha =2$ would correspond to a stable initial thermal distribution at finite chemical potential, and thus any $\alpha \gt 2$ is required to set off a re-equilibration to a lower temperature. In figure 2 we show the ensuing time evolution of the condensate occupation number, $n(0,t)$, for the different α. In each case, the evolution leads to a BEC characterized by a nonvanishing $n(0,t)$. As we chose the same $n(0,0)$ and $k_{0}^{\alpha }$ for each α, the final condensate fraction, $n(0,t)/N$, of the total particle number N grows with α. This is because the larger α cuts off more high-momentum particles and leaves less energy to be thermally redistributed. Power-law growth $n(0,t)\sim t$, predicted in [6], and later ∼t³ [9], is seen for near-thermal $\alpha =2.5$, while for large α the late-time growth reduces to ∼t². As a result, strong cooling quenches qualitatively modify the way the BEC grows, and a slowing down of this process is seen.

Zoom In Zoom Out Reset image size

Beyond the zero mode, the occupation spectrum of the nonzero momentum modes allows us to follow the dynamical process of strongly wave-turbulent particle transport to lower energies, and to identify a further smoking gun for the approach of the non thermal fixed point. In figure 3 we show the time evolution of the single-particle distribution, $n(k,t)$, over the absolute momenta, $k=|{\bf k}|$, at four different times and for different α. During the initial evolution, ($t\lesssim {{10}^{2}}\tau$), many particles gradually move to lower wave numbers, while at the same time relatively few particles deposit the surplus energy in the high-momentum modes, refilling them again. According to [10], a weak wave-turbulence inverse particle cascade with $n(k)\sim {{k}^{-2.4}}$ is expected within the range of the validity of kinetic theory. After a while, however, the spectra developing from the different initial α differ strongly. For $\alpha \gtrsim 3$, the distribution develops a bimodal structure, with $n(k)\sim {{k}^{-5}}$ in the infrared (IR) and $n(k)\sim {{k}^{-2}}$ in the (UV). At very long times, this bimodal structure decays towards a global $n(k)\sim {{k}^{-2}}$ (not shown). For $\alpha \lesssim 3$, the distribution directly reaches a thermal Rayleigh–Jeans scaling, $n(k)\sim T/{{k}^{2}}$.

To interpret our results in the context of NTFPs, we decompose kinetic-energy spectra as in [28] (see the appendix for details) into contributions from incompressible and compressible flow and quantum pressure fluctuations. In figure 4(a), we show the evolution of the different components, ${{n}_{\delta }}(k)$, $\delta \in \;\{{\rm i},{\rm C},{\rm q}\}$, for the weak initial quench, $\alpha =2.5$: The incompressible component, n_i, which accounts for vortical flow with a velocity vector changing transversally to its direction is roughly equal in size to the compressible component, n_c, of the longitudinal density fluctuations (i.e., sound excitations). At the same time, the quantum pressure part, n_q, is insignificant on all scales. For $t\lesssim {{10}^{3}}\tau$, due to the absence of phase coherence [25], the resulting spectra do not add up to the single-particle spectrum, $n(k)\ne {{n}_{{\rm i}}}(k)+{{n}_{{\rm C}}}(k)+{{n}_{{\rm q}}}(k)$. n(k) grows in the regime of low momenta, meaning that phase coherence is established and growing, and a condensate fraction appears (figure 2). For the case of the strongly nonthermal initial distribution, $\alpha =10$, the evolution is shown in figure 4(b). In contrast to figure 4(a), two macroscopic flows can be observed: one to the UV, and one to the IR. Conservation of particle and energy immediately imply that, when sent out from the regime of intermediate frequencies, energy is deposited in the UV, while particles are predominantly transferred to the IR. This leads to an inverse particle cascade with approximately k-independent radial particle flux, $Q(k)\equiv Q$, and a corresponding direct energy cascade to the UV [25]. The inverse particle cascade reflects strong wave turbulence, characterized by $n(k)\sim {{k}^{-5}}$ [27]. The decomposition in figure 4(b) makes clear that this power law is caused by incompressible excitations alone [25, 29], establishing a dominantly ideal hydrodynamic (superfluid) BE condensation process. In the UV, the excitations follow a thermal $n(k)\sim {{k}^{-2}}$ and are dominated by n_c and n_q.

Zoom In Zoom Out Reset image size

The above results show that during the hydrodynamic condensation process, incompressible flow temporarily dominates in the IR regime at the expense of compressible excitations. The opposite occurs for the compressible excitations in the UV. This dynamical scale separation and bimodal power-law distribution is a signature of the system approaching the NTFP.

In figure 5 we show, for intermediate times, the three-dimensional distribution of points where the density falls below $0.2\%$ of the average density $\bar{n}$, for the systems quenched with $\alpha =2.5$ and $\alpha =10$. We filtered out modes with wavenumbers larger than $k\xi =0.45$, but in the hydrodynamic case the high-momentum fluctuations barely distort the figure. The vortex tangles corroborate the findings of [17] for both cases of α. However, a remarkable difference exists in the distribution of the phase angle, $\varphi ({\bf x})$, of the Bose field, as can be inferred from figure 4. While after weak quenches, the circular flow has strong longitudinal (compressible) fluctuations, in the evolution passing the NTFP, macroscopic quantized vortical flow and distinct Vinen tangles [30] appear in the superfluid.

Zoom In Zoom Out Reset image size

The evolution of the vortex distribution, together with the phase coherence building up in the gas, which eventually leads to a fully coherent BEC on a low-temperature background, allows us to get a complementary understanding of the approach to and departure from the NTFP. In figure 6, we show the evolution of the system, again starting from the different initial quenches labelled by α, in a reduced phase space defined by two different characteristic length scales, in units of the healing length, ξ: The coherence length, ${{l}_{{\rm C}}}$, is a measure for the mean spatial falloff of phase coherence in space. We define it as the integral over the angle-averaged first-order coherence function, ${{g}^{(1)}}(r)$ (see appendix A.3). The correlation length, ${{l}_{{\rm d}}}$, measures the mean decay of vorticity in the proximity of the vortex line defects in the system. For an isolated circular vortex ring, ${{l}_{{\rm d}}}$ is proportional to the diameter of the ring and thus has a similar relevance as the mean distance between vortices and antivortices in a two-dimensional superfluid. The trajectories of the quench-cooled Bose gas in the $({{l}_{{\rm d}}},{{l}_{{\rm C}}})$-plane clearly exhibit the NTFP. Lying in the lower left corner, it corresponds to a strongly coherent gas with a maximum mutual separation and minimum bending of the vortex filaments. The NTFP corresponds ideally to a single vortex ring of maximum extent within the volume, which is known to be a metastable configuration decaying only slowly via bending and sound generation [22, 31]. The marked time steps show the critical slowing down of the system approaching the NTFP, while colliding vortex rings are seen to reconnect and form larger rings. After a long period in the vicinity of the fixed point, the system eventually departs towards final thermal equilibrium. This process corresponds to the shrinking of the last remaining vortex ring by transferring energy to the incoherent sound excitations and particles to the condensate mode. The NTFP sits at the crossroads of attractive and repulsive directions in our reduced phase space.

Zoom In Zoom Out Reset image size

Our findings suggest that a transition occurs, from a direct evolution towards thermal equilibrium to an evolution that first approaches the nonthermal fixed point, at a value of α between 3 and 4. Each of the trajectories eventually reaches a thermal configuration in the phase with a condensate present. Hence, the transition is expected to occur smoothly. As seen in figure 3, for $\alpha =3$ the occupation number spectrum already indicates the formation of vortex excitations, while the relatively high occupation of the higher-energy modes prevents a clear separation of incompressible and compressible fractions. While a prediction of the transition value of α is beyond the scope of the present work, this value depends on the system parameters, and the total particle number and energy after the quench, which determine the final temperature and condensate fraction. Since nonlinear interactions are involved in the redistribution process, the transition value is also expected to depend on the particular nonequilibrium distribution of energy and particle numbers across the different momenta.

Since the dynamics we are interested in exclusively affects the low-momentum, strongly populated field modes, we employ the so-called classical field method, which yields, within numerical accuracy and the classical-wave approximation, exact results for the time-evolving observables [39, 40]. For this, initial field configurations, $\phi (k,{{t}_{0}})$, are sampled from Gaussian probability distributions and then propagated according to the classical equations of motion. At the end of the time evolution, correlation functions are obtained from ensemble averages over the set of sampled paths. We study the dynamics of a dilute Bose gas by statistically sampling the classical equation of motion $(\hbar =1)$,

$$\begin{eqnarray}&&{\rm i}{{\partial }_{t}}\phi ({\bf x},t)=\left[ -\frac{{{\nabla }^{2}}}{2\;m}+g|\phi ({\bf x},t){{|}^{2}} \right]\phi ({\bf x},t).\end{eqnarray} \tag{ A.1 }$$

We consider a gas of N atoms in a box of size L³, with periodic boundary conditions and mean density $\bar{n}=N/{{L}^{3}}$. Lengths are measured in units of the healing length $\xi ={{(2\;mg\bar{n})}^{-1/2}}$ and time in units of $\tau =m{{\xi }^{2}}$. Simulations were done on a cubic grid with 256³ points. We have explored the dependence on grid size in our previous work [25, 41]. There is no visible dependence of the power laws found in the IR if the chosen grid is sufficiently large, as it is here.

The initial field in momentum space, $\phi ({\bf k},0)=\sqrt{n({\bf k},0)}{\rm exp} \{i\varphi ({\bf k},0)\}$, is parametrized in terms of a randomly chosen phase, $\varphi ({\bf k},0)\;\in \;[0,2\pi )$, and a density, $n({\bf k},0)=f(k){{\nu }_{{\bf k}}}$, with ${{\nu }_{{\bf k}}}\geqslant 0$ drawn from an exponential distribution, $P({{\nu }_{{\bf k}}})={\rm exp} (-{{\nu }_{{\bf k}}})$, for each ${\bf k}$. We choose the structure function,

$$\begin{eqnarray}&&f(k)=\frac{{{f}_{\alpha }}}{k_{0}^{\alpha }+{{k}^{\alpha }}},\end{eqnarray} \tag{ A.2 }$$

for different values of α, with cutoff ${{({{k}_{0}}\xi )}^{\alpha }}=0.2/0.{{44}^{\alpha }}$ and normalization ${{f}_{\alpha }}=400/0.{{44}^{\alpha }}$. We compare results for a range of different cooling quenches defined by the power laws $\alpha =2.5,...,10.0$, varying the total number between $N={{10}^{9}}$ ($\alpha =2.5$) and $N=4.3\times {{10}^{8}}$ ($\alpha =10$).

To interpret our results in the context of superfluid turbulence, we analyze kinetic-energy spectra as in [42]. Using the polar coordinate representation of the complex field $\phi =\sqrt{n}\;{\rm exp} \{i\varphi \}$ and the definition ${\bf v}=\nabla \varphi /m$ of the velocity field, the kinetic energy, ${{E}_{{\rm kin}}}=\int {\rm d}{\bf x}\;\langle |\nabla \phi ({\bf x},t){{|}^{2}}\rangle /2$, is split: ${{E}_{{\rm kin}}}={{E}_{{\rm v}}}+{{E}_{{\rm q}}}$, into the 'classical' component ${{E}_{{\rm v}}}=\int {\rm d}{\bf x}\;\langle |\sqrt{n}{\bf v}{{|}^{2}}\rangle /2$, with a 'quantum-pressure' component, ${{E}_{{\rm q}}}=\int {\rm d}{\bf x}\;\langle |\nabla \sqrt{n}{{|}^{2}}\rangle /2$. Radial particle spectra

$$\begin{eqnarray}&&{{n}_{\delta }}(k)=\frac{1}{2}\int {\rm d}\Omega \;\langle |{{{\bf w}}_{\delta }}({\bf k}){{|}^{2}}\rangle ,\quad \delta ={\rm v},{\rm q}\end{eqnarray} \tag{ A.3 }$$

in terms of the generalized velocities, ${{{\bf w}}_{{\rm v}}}=\sqrt{n}{\bf v}$ and ${{{\bf w}}_{{\rm q}}}=\nabla \sqrt{n}$, are furthermore decomposed by ${{{\bf w}}_{{\rm v}}}={{{\bf w}}_{{\rm i}}}+{{{\bf w}}_{{\rm C}}}$ into 'incompressible' ($\nabla \cdot {{{\bf w}}_{{\rm i}}}=0$) and 'compressible' ($\nabla \times {{{\bf w}}_{{\rm C}}}=0$) parts to distinguish the vortical superfluid and sound excitations of the gas, respectively. Note that the factor $\sqrt{n}$ in ${{{\bf w}}_{{\rm v}}}$ regularizes the velocity field in the UV [42]. The incompressible component reflects transverse motion (i.e., the flow which changes perpendicularly to its direction). It arises mainly from superfluid vortical flow. The compressible component accounts for the longitudinal flow corresponding to sound wave excitations.

The evolution of the integrated fraction, $\Delta (t)={{N}_{c}}({{k}_{\lambda }},t)/{{N}_{i}}({{k}_{\lambda }},t)$, with ${{N}_{\delta }}({{k}_{\lambda }},t)={{\int }_{|{\bf k}|\lt {{k}_{\lambda }}}}{\rm d}{\bf k}\;{{n}_{\delta }}({\bf k},t)$, of compressible to incompressible occupations below a momentum scale, ${{k}_{\lambda }}\xi =0.35$, is shown in figure A1 . Initially, $\Delta \sim 0.5$. Starting from $\alpha \leqslant 3$, Δ stays approximately constant, while for $\alpha \gt 3$, Δ decreases for $t\gtrsim {{10}^{2}}\tau$ towards zero before increasing again when thermal equilibrium is approached. The inset shows $\Delta (t=606)$ as a function of α in which one identifies a transition to a separation of the components, depending on the strength of the initial cooling quench.

Zoom In Zoom Out Reset image size

The coherence length, ${{l}_{{\rm C}}}$, is defined as the integral of the angle-averaged first-order coherence function

$$\begin{eqnarray}&&{{g}^{(1)}}(r)=\int {\rm d}\Omega \frac{\langle {{\phi }^{*}}({\bf x})\phi ({\bf x}+{\bf r})\rangle }{\sqrt{\langle n({\bf x})\rangle \langle n({\bf x}+{\bf r})\rangle }},\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&{{l}_{{\rm C}}}=\int {\rm d}r\;{{g}^{(1)}}(r),\end{eqnarray} \tag{ A.5 }$$

which becomes independent of ${\bf x}$ in the ensemble average. In contrast to ${{r}_{{\rm coh}}}=\int {\rm d}r\;{{r}^{3}}\;{{g}^{(1)}}(r)/\int {\rm d}r\;{{r}^{2}}\;{{g}^{(1)}}(r)$, our choice of the coherence length, ${{l}_{{\rm C}}}$, does not enlarge insignificant contributions at large r.

The vortex-correlation length is defined in terms of the vorticity. The curl of the velocity field, ${\bf u}({\bf x})={\bf rot}\;{\bf v}({\bf x})=\nabla \times \nabla \varphi ({\bf x})/m$, vanishes except at the positions of topological defects. At these phase defects, it yields the quantization and direction of the vortex line. Based on the vorticity, ${\bf u}({\bf x})$, we define angle-averaged correlations of the vorticity as

$$\begin{eqnarray}&&C({\bf x},r)=\int {\rm d}{{\Omega }_{{\bf r}}}\;{\bf u}({\bf x}){\bf u}({\bf x}+{\bf r}).\end{eqnarray} \tag{ A.6 }$$

For small longitudinal distances, r, along the vortex filament, the vorticity, ${\bf u}$, points in approximately the same direction as at ${\bf x}$, rendering $C({\bf x},r)$ positive. In general, however, the vorticity shrinks with growing r, and the first transverse zero, ${{r}_{0}}({\bf x})$, marks the smallest distance, at which the vorticity points predominantly in the opposite direction due to the curving of the same vortex line or the proximity of other vortex lines. Therefore, ${{r}_{0}}({\bf x})$ has a meaning similar to the pair distance of anticirculating vortices in two-dimensional gases, and works equally well for vortex tangles and vortex rings. ${{l}_{{\rm d}}}$ is defined as the average of the transitions through zero,

$$\begin{eqnarray}&&{{l}_{{\rm d}}}={{\mathcal{N}}^{-1}}{{\int }_{\mathcal{C}}}{\rm d}l\;{{r}_{0}}(l),\end{eqnarray} \tag{ A.7 }$$

taking the normalization integral over all vortex lines, and with $\mathcal{N}={{\int }_{\mathcal{C}}}{\rm d}l$. Accounting for the transition through zero for each ${\bf x}$ has the advantage that, for small and large vortex rings, r₀ is in effect calculated separately and averaged with a weight corresponding to the length of each ring. If one took the transition through zero of the averaged correlation function, $c(r)=[\int {{{\rm d}}^{3}}x\;C({\bf x},r)]/\int {{{\rm d}}^{3}}x\;|{\bf u}({\bf x}){{|}^{2}}$, the positive contribution of a large ring at some distance, r, could be of a similar order as the negative contributions of smaller rings (see figure A2 ). If both contributions almost cancel each other, small fluctuations around zero could lead to discontinuous jumps of r₀ with time.

Zoom In Zoom Out Reset image size