Universal dynamics on the way to thermalisation

It is demonstrated how a many-body system far from thermal equilibrium can exhibit universal dynamics in passing a non-thermal fixed point. As an example, the process of Bose-Einstein (BE) condensation of a dilute cold gas is considered. If the particle flux into the low-energy modes, induced, e.g., by a cooling quench, is sufficiently strong, the Bose gas develops a characteristic power-law single-particle spectrum $n(k)\sim k^{-5}$, and critical slowing down in time occurs. The fixed point is shown to be marked by the creation and dilution of tangled vortex lines. Alternatively, for a weak cooling quench and particle flux, the condensation process runs quasi adiabatically, passing by the fixed point in far distance, and signatures of critical scaling remain absent.


Introduction
Systems like a ferromagnet or a bosonic quantum gas can undergo a second-order phase transition in the Ehrenfest classification, at which their relevant physical properties become independent of microscopic details. This independence leads to the concept of universality, which has become extremely successful in classifying and characterizing matter in thermal equilibrium. Many different phenomena can be characterized in terms of a few classes governed by the same critical properties. Generalizing this concept to universal dynamics far away from thermal equilibrium requires us to develop an understanding of the structure of the space of possible nonthermal states. In addition to the well-known dynamics near thermal fixed-points [1], critical regions such as fixed points, critical lines, and surfaces can affect the evolution far from equilibrium. When a closed system approaches such critical configurations, memory about its particular initial state would be partially lost. Critical slowing down in the actual time evolution would be observed since the largest scales, not only in space but also in time, dominate the system dynamics. Thermal states are attractive, stable fixed points in this framework. Different systems would be related by means of the universality class that their dynamical evolution falls into. Predictions for the behavior of very different physical systems could be obtained on the basis of comparatively few exemplary measurements. When looking at fundamental science applications, this could link the dynamics observed in very different contexts.
Here we show that such universal time evolution is possible, using the example of a strongly shock-cooled three-dimensional normal-fluid Bose gas, which approaches a nonthermal fixed point [2,3] before it eventually proceeds to a thermal distribution below the BE critical temperature. The cooling quench puts the system far out of equilibrium so that it can undergo an evolution passing the vicinity of a critical point, where it gets stuck for a long time due to the critically enhanced role of ultra-low-energy modes. Looking at the critical configuration in real space, we find it to be dominated by a dilute ensemble of quasi-topological vortex line defects, which only very slowly decay via interactions with sound.
In this article, we demonstrate how the superfluid turbulence period is set into the context of universal dynamics. Depending on whether the condensation evolution passes close to or further away from a nonthermal fixed point, the process appears in qualitatively different forms (figure 1). Depending on the strength, α, of an initial cooling quench, the gas can thermalize directly to a Bose-Einstein condensate, or it can first approach and critically slow down near a nonthermal fixed point (NTFP). Therefore, it is characterized by a self-similar particle spectrum, ∼ − n k k ( ) 5 , and a diluting ensemble of vortex rings that are growing in size.

Adiabatic versus turbulent condensation dynamics
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 k T k ( ) B 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, , with an exponent, ζ = 5, that is distinctly larger than the exponent ζ = 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.

Semiclassical simulations
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, ϕ t x ( , ) (see appendix A.1 for details). We compute the evolution in a cubic box with periodic boundary conditions. The initial field in momentum space, ϕ k ( , 0), is sampled from a distribution that allows its phase to be random, while ϕ k | ( , 0)| 2 has Gaussian fluctuations around a momentum distribution, n k ( , 0), which is flat between = = k k | | 0 and k ≃ k 0 , and falls off as α − k for k ≳ k 0 (top left panel of figure 3).

Results: evolution of the zero mode
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, α ≃ 2.4 [10], corresponding to an inverse particle cascade in weak wave turbulence theory [6,7]. α = 2 would correspond to a stable initial thermal distribution at finite chemical potential, and thus any α > 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 t (0, ), for the different α. In each case, the evolution leads to a BEC characterized by a nonvanishing n t (0, ). As we chose the same n (0, 0) and α k 0 for each α, the final condensate fraction, n t N (0, ) , 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 t t (0, ) , predicted in [6], and later ∼t 3 [9], is seen for near-thermal α = 2.5, while for large α the late-time growth reduces to ∼t 2 . As a result, strong cooling quenches qualitatively modify the way the BEC grows, and a slowing down of this process is seen.

Evolution of the single-particle momentum spectrum
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 k | |, at four different times and for different α. During the initial evolution, ( τ ≲ t 10 2 ), 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 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 α ≳ 3, the distribution develops a bimodal structure, with ∼ − n k k ( ) 5 in the infrared (IR) and ∼ − n k k ( ) 2 in the (UV). At very long times, this bimodal structure decays towards a global ∼ − n k k ( ) 2 (not shown). For α ≲ 3, the distribution directly reaches a thermal Rayleigh-Jeans scaling, ∼ n k T k ( ) 2 .

Hydrodynamic condensation
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 k ( ), δ ∈ {i, c, q}, for the weak initial quench, α = 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 10 3 , due to the absence of phase coherence [25], the resulting spectra do not add up to the single-particle spectrum, . 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, α = 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, , and a corresponding direct energy cascade to the UV [25]. The inverse particle cascade reflects strong wave turbulence, characterized by ∼ − n k 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 k ( ) 2 and are dominated by n c and n q . 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. 3. A flow picture of the approach of the nonthermal fixed point

Defect formation and dilution
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 n, for the systems quenched with α = 2.5 and α = 10. We filtered out modes with wavenumbers larger than ξ = k 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, φ 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.

Representation in a reduced phase space
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 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 r ( ) (1) (see appendix A.3). The correlation length, l 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 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 l ( , ) D 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. 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 α = 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.

Conclusions
The process of Bose-Einstein condensation in a quench-cooled, dilute cold gas can show features of universal dynamics. Provided a sufficiently strong cooling quench, the condensing system passes by a partially attractive nonthermal fixed point (NTFP) where it is critically slowed down. The approach of the fixed point is marked by the appearance of incompressible flow around tangled vortex lines. In this regime, particles cannot be deposited quickly enough into the zero mode, so they form an excess population with a characteristic power-law falloff within the low-energy modes. In contrast, slow, near-adiabatic condensation can exhibit the appearance of vortical motion, which is, however, distorted by strong compressible sound excitations. The critical slowing down of the phase coherence length and vortex distance provide smoking guns for the detection of the universal dynamics in experiments. A complete characterization of NTFP, in terms of a full set of critical exponents and thus universality classes, including anomalous dimensions, is most desirable, because it would expand the theory of critical phenomena far away from equilibrium. Understanding the possible different paths to a BEC is of fundamental interest beyond the realm of ultracold gases, from the phenomenology of the solid state up to the highest energies (e.g., in heavy-ion collisions [32][33][34][35][36] or earlyuniverse evolution [2,26,32,37,38]).

A.1. Semiclassical simulations
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, ϕ 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 = We consider a gas of N atoms in a box of size L 3 , with periodic boundary conditions and mean density = n N L 3 . Lengths are measured in units of the healing length ξ = − mgn (2 ) 1 2 and time in units of τ ξ = m 2 . Simulations were done on a cubic grid with 256 3 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.

A.2. Hydrodynamic decomposition
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 ϕ φ = n i exp { } and the definition φ =  m v of the velocity field, the kinetic energy, , is split: , with a 'quantumpressure' component, ) parts to distinguish the vortical superfluid and sound excitations of the gas, respectively. Note that the factor n in w 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