Universal dynamics and non-thermal fixed points in quantum fluids far from equilibrium

Closed quantum systems far from thermal equilibrium can show universal dynamics near attractor solutions, known as non-thermal fixed points, generically in the form of scaling behavior in space and time. A systematic classification and comprehensive understanding of such scaling solutions are tasks of future developments in non-equilibrium quantum many-body theory. In this tutorial review, we outline several analytical approaches to non-thermal fixed points and summarize corresponding numerical and experimental results. The analytic methods include a non-perturbative kinetic theory derived within the two-particle irreducible effective-action formalism, as well as a low-energy effective field theory framework. As one of the driving forces of this research field are numerical simulations, we summarize the main results of exemplary cases of universal dynamics in ultracold Bose gases. This encompasses quantum vortex ensembles in turbulent superfluids as well as recently observed real-time instanton solutions in one-dimensional spinor condensates.

In this brief tutorial review, we focus on universal dynamics of dilute Bose gases close to a non-thermal fixed point.Universality here means that the evolution after some time becomes to a certain extent independent of the initial condition as well as of microscopic details.The universal intermediate state, that develops, is determined only by symmetry properties and possibly a limited set of relevant quantities and/or functions pre-determined by the initial configuration.Generically, this allows categorizing systems into universality classes based on their symmetry properties and the family of far-from-equilibrium states the initial condition belongs to.
Our article is organized as follows.In Sect.II we introduce the main concepts of non-thermal fixed points.Sect.III contains a summary of the main theoretical approaches to describing non-thermal fixed points in ultracold quantum gases.In Sect.IV, we compare the analytical predictions with numerical simulations and discuss the role of non-linear (topological) excitations.Sect.V summarizes experimental results on non-thermal fixed points.We close our tutorial review with an outlook to future research in the field, see Sect.VI.

II. NON-THERMAL FIXED POINTS
The concept of non-thermal fixed points is motivated by the ideas of (near-)equilibrium renormalization group (RG) theory.Generalizing fixed points of RG flow equations, which characterize, e.g., critical phenomena in (thermal) equilibrium, non-thermal fixed points appear in time evolution flows out of equilibrium.This includes, in particular, universal, selfsimilar evolution and the transient appearance of largely scalefree spatial patterns.Associated with relaxation of closed systems, they are typically subject to conservation laws.In this chapter, we summarize the main concepts of non-thermal fixed points.

A. Universal scaling
In the RG framework, one studies a physical system in a way which resembles looking at it through a microscope at different resolutions.Close to a critical point, one typically observes that the system looks self similar, i.e., it does not change its appearance when varying the resolution.
As a simple example, consider a two-point correlation function C(x; s) of some locally measurable observable, which, if the system is homogeneous and isotropic, depends only on the distance x = |r 1 − r 2 | between two positions r i in space.The second argument, s, is a number that defines the resolution in units of a fixed length scale and represents the flow parameter of the RG.Changing the value of s, the correlation function C(x; s) should change accordingly.Self-similarity means that C(x; s) rescales as C(x; s) = s ζ f (x/s).This implies that the correlations are solely characterized by a universal exponent ζ and a scaling function f .
A fixed point of the RG flow equation corresponds to the case when the system becomes fully s-independent, which happens when f (x) ∼ x ζ .Typically, however, for a realistic physical system, the fixed point is partially repulsive.In this case, the scaling function f retains some information about characteristic scales, such as a correlation length ξ, and therefore does not assume a pure power-law form.The system's RG flow only approaches the fixed point but generically does not reach it before being driven away again.Consider, for example, a continuous phase transition in equilibrium, at which the correlation length diverges.The (fine-tuned) system can be precisely at the RG fixed point only in the thermodynamic limit, which allows having a diverging correlation length and thus a pure scaling form describing its correlations at any finite scale.
Taking the evolution time t as the scale parameter, the renormalization-group idea can be extended to the time evolution of non-equilibrium systems.The corresponding fixed point of the RG flow is called a non-thermal fixed point.In the scaling regime near a non-thermal fixed point, the evolution of the time-dependent version of the correlation function introduced above is determined by C(x; t) = t α f (t −β x), with now two universal exponents α and β that assume, in general, nonzero values.The associated correlation length of the system changes as a power of time, ξ(t) ∼ t β .Note that the time evolution taking power-law characteristics is equivalent to critical slowing down, here in real time.We remark that, depending on the sign of β, increasing the time t can correspond to either a reduction or an increase of the microscope resolution.
In general, the scaling exponents α and β, together with the scaling function f , allow us to determine the universality class associated with the fixed point [94], but they may not form a sufficient criterion for that.It is, in particular, expected that the evolution of very different physical systems far from equilibrium can be categorized by means of their possible kinds of spatio-temporal scaling behavior.A full classification of such universality remains an open problem.However, similar to the case of equilibrium critical phenomena, underlying symmetries of the system are expected to play a crucial role.

Direct thermalization
Self-similar dynamics Far from equilibrium Close to equilibrium FIG. 1. Schematics of different scenarios of thermalization.Within a subclass of far-from-equilibrium conditions all the states undergo the same self-similar evolution regime before reaching equilibrium.In contrast, a generic close-to-equilibrium initial state thermalizes directly without any universal scaling dynamics in between.Figure adapted from [48].
Although the evolving system, close to a non-thermal fixed point, forgets about many details of where it comes from, in analogy to equilibrium RG flows, the initial conditions of the flow are not entirely irrelevant.Whether a physical system will approach a non-thermal fixed point and show universal scaling dynamics, or which fixed point it will be able to reach, in general depends on the particular initial state.Going back to the RG analogy, one can imagine a space of all possible states.The evolution of one state to another can be represented as a trajectory in this space.A set of all the trajectories forms a flow in the state space, similar to a flow of coupling constants in the RG theory or to a phase portrait of some dynamical system.
While the asymptotic state is typically expected to correspond to one of the system's possible equilibrium configurations, there can be attractors near which the evolution is critically slowed down.These attractors are exactly the aforementioned non-thermal fixed points.Therefore, in general, the whole space can be divided into regions that are attracted to different non-thermal fixed points.At the same time, some initial conditions may not lead to a non-thermal fixed point at all but instead to direct thermalization, see Fig. 1.It is commonly accepted, however, that the key precondition for the system to reach universal self-similar scaling dynamics is an extreme out-of-equilibrium initial configuration characterized by either strong statistical fluctuations or a strong (inhomogeneous) mean field.

B. Self-similar transport
As a relevant example, consider the time evolution of a single-component dilute gas of bosonic atoms in three spatial dimensions, described by the classical Gross-Pitaevskii (GP) field equation of motion, particles removed by strong cooling quench initial distribution after quench e n e r g y p a r t i c l e s log k 2. Self-similar scaling in time and space close to a non-thermal fixed point.The sketch shows, on a double-logarithmic scale, the time evolution of the single-particle momentum distribution n(t, k) of a Bose gas for two different times t (solid and short-dashed lines).Starting from an extreme initial distribution marked by the red long-dashed line, being the result of a strong cooling quench, a bi-directional redistribution of particles in momentum space occurs as indicated by the arrows.Particle transport towards low momenta as well as energy transport to larger momenta are characterized by self-similar scaling evolutions in space and time according to , with universal scaling exponents α and β, different for both directions.Here, t ref is an arbitrary reference time within the temporal scaling regime.The infrared transport (green arrow) conserves the particle number, which is concentrated at small momenta.In contrast, the energy, being concentrated at high momenta, is conserved in the redistribution of short-wavelength fluctuations (blue arrow).See main text for details. Figure adapted from [96].
where M is the atom mass, and g = 4πa/M, with s-wave scattering length a, is a coupling constant.This coupling, multiplying the local density ρ(x, t) = |ψ(x, t)| 2 , quantifies the interaction 'potential'.Here and in the following, we choose natural units in which = 1.
The system can approach a non-thermal fixed point as the result of a strong initial cooling quench [96], see Fig. 2 as well as [94,[108][109][110].An extreme version of such a quench can be achieved, e.g., by first cooling the system adiabatically such that its chemical potential is 0 < −µ k B T , where the temperature T T c is just above the critical temperature T c separating the normal and the Bose condensed phases of the gas, and then removing all particles with energy higher than ∼ |µ|.This leads to a distribution that drops abruptly above a momentum scale Q, with zero-mode occupation n 0 and Heaviside function Θ, see the red dashed line in Fig. 2. If the corresponding energy is on the order of the ground-state energy of the post-quench fully condensed gas with uniform density ρ, Q 2 /2M |µ| gρ, then the majority of the energy of the gas after the quench is concentrated at the scale Q k ξ , the healing-length momentum scale k ξ = 8πaρ.
Most importantly, such a strong cooling quench leads to an extreme initial condition for the subsequent dynamics.
The post-quench distribution is strongly over-occupied at momenta k < Q, as compared to the final equilibrium distribution.This initial overpopulation of modes with energies ∼ Q 2 /2M induces inverse particle transport from intermediate to lower momenta, while energy is transported to higher wave numbers [94,108,109], as indicated by the arrows in Fig. 2. The overall transport, which subsequently develops is thus characterized by a bi-directional, in general non-local redistribution of particles and energy.This transport requires interactions, i.e., collisions between the particles in the gas, which give rise to energy and momentum exchange, allowing certain particles to loose momentum and energy while others speed up in their motion.This is illustrated in Fig. 2. For example, particles making up the over-occupation at intermediate momenta close to the scale Q, which are being transferred to increase the occupancy of modes of lower momenta, loose a considerable part of their kinetic energy (green arrow, note the logarithmic scales).Hence, in order for the total energy conservation to be satisfied, other particles need to be scattered to higher-momentum modes within the tail (blue arrow).
The evolution eventually becomes universal in the sense that it is then approximately independent of the precise initial conditions set by the cooling quench as well as of the particular values of the physical parameters characterizing the system.In the vicinity of a non-thermal fixed point, the momentum distribution of the Bose gas rescales self-similarly, within a certain range of momenta, according to with some reference time t ref .
The distribution shifts to lower momenta for β > 0, while transport to larger momenta occurs in the case of β < 0. A bi-directional scaling evolution is, in general, characterized by two different sets of scaling exponents.One set describes the inverse particle transport towards low momenta whereas the second set quantifies the transport of energy towards large momenta.

C. Scaling function
While the spatio-temporal scaling provides the 'smoking gun' for the approach of a non-thermal fixed point, in all cases examined so far, also power-law scaling of the momentum distribution, n(k) ∼ k −ζ , has been observed and reflects the character of the underlying transport, see Fig. 2. In both, the infrared (IR) regime of inverse transport to lower momenta (β > 0) and the ultraviolet (UV) range, in which a direct transport to higher momenta prevails (β < 0), the distribution function typically assumes a (potentially) different power-law form.At any finite time after the quench, both, the IR and the UV distributions are cutoff at some scale k Λ , below which n(t, k) flattens out, and k λ , above which it more steeply, e.g., exponentially falls to zero.Both, k Λ and k λ , in general vary in time as a result of the transport, as indicated in Fig. 2. Evaluated at a fixed reference time t ref , the fixed-point solution (3) further defines the universal scaling function f s (k) = n(t ref , k).Within a limited range of momenta, it satisfies the scaling hypothesis f s (k) = s ζ f s (sk), with an additional, in general independent scaling exponent ζ.A frequently used simple ansatz for the scaling function f s (Q) in the IR region is given by It interpolates between the universal power-law behavior f s (k) ∼ k −ζ for k > k Λ and the plateau region f s (k) ∼ const.below the running scale k Λ , see the inset of Fig. 2. Combining the spatio-temporal scaling form (3) with the scaling function (4) gives that the momentum scale evolves as k Λ (t) ∼ t −β , corresponding to a characteristic length scale growing as Λ (t) ∼ t β .

D. Conservation laws
Global conservation laws -applying within a certain, extended regime of momenta -strongly constrain the redistribution underlying the self-similar dynamics in the vicinity of the non-thermal fixed point.Hence, they play a crucial role for the possible scaling evolution as they impose scaling relations between the scaling exponents.For example, the conservation of the total particle number, d d k n(t, k) = N(t) ≡ N, with n(t, k) evolving according to Eq. ( 3), in d spatial dimensions, requires that α = dβ.
In a closed system, both, the total energy and particle number, need to be conserved by the transport.For the bidirectional transport sketched in Fig. 2, the inverse flow is dominated by particle-number conservation, while the highmomentum modes accumulate the major part of the kinetic energy.For this to be the case, the power-law exponents ζ of n(k) ∼ k −ζ can be within a certain range of values only [96,111].For example, in the simpler case that ζ is the same everywhere between the IR and UV cutoff scales, k Λ k k λ , one needs to have d < ζ < d + 2 in d spatial dimensions, for the particle, ∼ n(t, k), and energy distributions, ∼ k 2 n(t, k), to be dominated by IR and UV scales, k k Λ and k k λ , respectively.Note that, only if this condition is fulfilled, the bi-directional transport can separate particle number and energy, which is one of the preconditions for self-similar universal scaling dynamics to occur.In the opposite case, for values of ζ, which let both, particles and energy to be concentrated at either side of the spectral range, scaling evolution will come out differently.The ensuing shock-wave-type redistributions in momentum space have been discussed in detail in [96,111], in the context of the build-up and decay of weak wave turbulence in classical systems.

E. Coarsening and phase ordering
The self-similar transport in momentum space can emerge from rather different underlying physical configurations and processes.For instance, the dynamics can be driven not only by the conserved redistribution of quasiparticle excitations such as in weak wave turbulence [94,96] but also by the reconfiguration of spatial patterns like magnetization domains [101,102] or by the annihilation of (topological) defects populating the system [103,108].The latter dynamics can be considered as the buildup of an inverse superfluid turbulent cascade [98,99,103].In contrast, if defects are subdominant or absent at all, which is the case, e.g., for U(N) symmetric models in the large-N limit [112], the strongly occupied modes exhibiting scaling near the fixed point [94,96] typically reflect strong phase fluctuations not subject to an incompressibility constraint.These can be described, e.g., by the re-summed kinetic theory discussed in Sect.III A or a low-energy effective theory, see Sect.III B and Ref. [113].The associated scaling exponents are generically different for both types of dynamics, with and without patterns or defects [94,100,103].
The concept of non-thermal fixed points thus includes scaling dynamics which exhibits coarsening and phase-ordering kinetics [62,63] following the creation of defects and nonlinear patterns after a quench, e.g., across an ordering phase transition.In most cases so far, such coarsening phenomena have been discussed within an open-system framework, considering the system to be coupled to a particle or heat bath.It is understood that the coupling to an external bath, which is usually described by means of a driven-diffusive model, can in general be realized also within a closed system, where part of the system, e.g., the high-energetic modes assume the role of the bath.From this point of view, the theory of non-thermal fixed points includes that of coarsening and opens an approach for capturing the entire scaling dynamics within a closed system from first principles, see, e.g.[106,107].

III. ANALYTICAL APPROACHES TO NON-THERMAL FIXED POINTS
After introducing the basic concept of non-thermal fixed points in the previous section we are set to discuss various methods employed to describe the universal scaling dynamics, focusing on analytical approaches in the present section.More detailed presentations of the formalism can be found, e.g., in Refs.[96,113].

A. Re-summed kinetic theory
A non-thermal fixed point is characterized by algebraic scaling in space and time towards smaller wave numbers, i.e., greater lengths, as formalized by the scaling form (3) for the single-particle momentum distribution, with the typical scaling function (4) defining the shape of the distribution.This implies the characteristic length scale to scale as k Λ ∼ t −β .
Consider a field theory such as the GP model (1) of a single-component dilute superfluid.In quantized form, the bosonic field operators obey the standard commutation relations [ψ(t, x), ψ(t, y) † ] = δ(x − y), [ψ(t, x), ψ(t, y)] = 0.For simplicity we restrict ourselves to a homogeneous system, e.g., a gas in a box with periodic boundary conditions, which one may describe in terms of the energy eigenmodes of some leading-order quasiparticle Hamiltonian.In the periodic box, these are plane waves with wave number k, e.g., free particle excitations with energy, i.e., frequency ω(k) = k 2 /2M or collective (sound) modes with ω(k) = c s |k|, with speed of sound c s = (gρ 0 /M) 1/2 , for a flat mean density ρ 0 .
In the following, we will restrict ourselves to the universal scaling dynamics of two-point functions.A simple example is the momentum distribution n(t, k), cf.Eq. ( 2).In quantum field theory, the exact time evolution of (in general unequaltime) two-point correlators is governed by the Kadanoff-Baym equations, cf., e.g., [71,95,114].These are derived within the (Baym-Kadanoff-)Schwinger-(Mahanthappa-Bakshi-)Keldysh formalism [115], typically in a path-integral setting, involving a closed time path C from some initial time t 0 to infinity and back to t 0 , along which the above time ordering T C of the field operators is defined.
In writing down the equations for G, one hides the generic dependence on all the arbitrary high correlations developing in the dynamical evolution of the interacting system in expressing the equations in terms of G (and the one-point function ψ(x) ) only.This comes at the cost that the equations are, in general, represented by an infinite series of Feynman diagrams made up of G and bare vertices.While in principle exact, a solution of these integro-differential equations is quite involved in practice, which makes them cumbersome for a theoretical analysis.For both, analytical insight and numerical evaluations, one usually needs to truncate the diagrammatic series and then still approximate the equations further to exhibit the mechanisms relevant at a non-thermal fixed point.
In the latter step, a crucial observation is that the scaling dynamics is reached at late times and low momenta, suggesting a slow dependence of the function G(x, y) on the central-time direction t ∼ x 0 + y 0 .This suggests an approximate description, known as the gradient expansion, that takes into account only low orders of both temporal and spatial central-coordinate, x + y, derivatives.One decomposes the time-ordered Green's function G(x, y) ≡ F(x, y) − (i/2) sgn C (x 0 − y 0 )ρ(x, y), with the sign function evaluating to ±1 for x 0 later/earlier than y 0 on the path C, into its symmetric 'statistical' F and anti-symmetric 'spectral' ρ components [71,95].This helps separating the information about the occupation number of the quasiparticle eigenmodes of the system ρ carries information about the spectral character of the quasiparticles, in particular their energy ω(k) and stability, i.e., spectral widths.These are approximately independent of the central time and space, x +y, and Fourier transformed with respect to the relative coordinate x − y, the resulting function ρ(ω, k), to a first approximation, looks like a delta distribution δ(ω − ω k ), i.e., a spectral distribution evaluating the frequency ω to the eigenfrequency ω k = ω(k) of momentum mode k.Hence, all frequencies k 0 = ω can easily be integrated out, such that the dynamic equations are left to involve F and thus n, depending on the central time t and the momenta only.
The statistical function F also contains information about the (quasi)particle distribution n(t, k) and therefore about the statistical occupancy of mode k, which is obtained by frequency integration over the statistical function F(t; ω, k).This corresponds to its equal-time entries Sending the initial time t 0 → −∞ one derives, at leading order in the gradient expansion, a quantum Boltzmann equation (QBE), for the time evolution of the the occupation number distribution n(t, k) = ψ † (t, k)ψ(t, k) .Here, I[n](t, k) is a scattering integral.Restricting ourselves to the case of elastic 2 ↔ 2 scatterings, the latter takes the form with T kpqr being the scattering T -matrix, for which we will later present specific expressions, and the (d + 1)-dimensional delta distributions imply energy and momentum conservation, with k 0 = ω k .The collision kernel under the integral (6) describes the redistribution of the occupations n k = n(t, k) of momentum modes k with eigenfrequency ω k due to elastic 2 ↔ 2 collisions from modes q and r into k and p and vice versa.But note that also collective scattering effects beyond 2 ↔ 2 processes can be captured in the T -matrix using, e.g., the re-summation techniques discussed in the following.
In presence of a Bose condensate, the occupation numbers describe quasiparticle excitations.Their properties enter the scattering matrix and the mode eigenfrequencies.Here, we consider transport entirely within the range of a fixed scaling of the dispersion ω k ∼ k z , with dynamical scaling exponent z, such that processes leading to a change in particle number are suppressed.
Two classical limits of the QBE scattering integral I[n](t, k) exist.The usual, Boltzmann integral for classical particles is obtained in the limit of n(t, k) 1.In the opposite case of large occupation numbers, n(t, k) 1, termed the classicalwave limit, the scattering integral reads Here, the QBE reduces to the so-called wave-Boltzmann equation (WBE), which is the subject of the following discussion.It best suits our interests, viz., in the universal dynamics of a near-degenerate Bose gas obeying n(t, k) 1 within the relevant, infrared momentum region.

Scaling of the scattering integral and the T -matrix
In the kinetic approximation, scaling features of the system at a non-thermal fixed point are directly encoded in the properties of the scattering integral.For a general treatment that governs the cases of presence and absence of a condensate density, we focus on the scaling of the distribution of quasiparticles, in the following denoted by n Q (t, k), instead of the single-particle momentum distribution n(t, k).Note that, in the case of free particles, with dispersion ω(k) = k 2 /2M ∼ k z , i.e., of a dynamical exponent z = 2, they are identical, n Q ≡ n.For Bogoliubov sound with dispersion ω(k) = c s k and thus z = 1, the scaling of n Q differs from the scaling of n due to the k-dependent Bogoliubov mode functions characterizing the transformation between the particle and quasiparticle basis, n(t, k) Using a positive real scaling factor s, the self-similar evolution of the quasiparticle distribution at a nonthermal fixed point reads We remark that, by choosing the scaling parameter s = (t/t ref ) β , one obtains the scaling form stated in the example in (3).
As the scattering integral, in the classical-wave limit, is a homogeneous function of momentum and time, it obeys scaling, provided the scaling (8) of the quasiparticle distribution, according to with scaling exponent µ = 2(d + m) − z − 3α/β.Here, m is the scaling dimension of the modulus of the T -matrix, Generally, this scaling hypothesis for the T -matrix does not hold over the whole range of momenta.In fact, scaling, with different exponents, is found within separate limited scaling regions, which we discuss in the next section.
Besides the spatio-temporal scaling, we would also like to derive the spatial scaling form, in particular the exponent ζ defined in (4).Consider, for this, the simple example of a universal quasiparticle distribution at a fixed time t 0 , which, at least in a limited regime of momenta, takes the pure powerlaw form, with fixed-time momentum scaling exponent κ.This requires also the T -matrix to show spatial momentum scaling at a fixed instance in time, with m κ being, in general, different from m.Note n Q and thus Eq. ( 11) in realistic cases is regularized by an IR cutoff k Λ , recall the function (4), and, analogously, by a UV cutoff k λ , to ensure that the scattering integral stays finite in the limits k k Λ and k k λ .

Perturbative region: two-body scattering
For the non-condensed, weakly interacting Bose gas away from unitarity, the T -matrix is well approximated by As the matrix elements are momentum independent we obtain m κ = m = 0.It can be shown that Eq. ( 13) represents the leading perturbative approximation of the full momentumdependent many-body coupling function.
In presence of a condensate density ρ 0 ≤ ρ, sound wave excitations become relevant below the healing-length momentum scale k ξ = 2gρ 0 M. Within leading-order perturbative approximation, the elastic scattering of these sound waves is described by the T -matrix [96] |T kpqr | 2 = (2π) 4 (Mc s ) 4  kpqr Hence, for the Bogoliubov sound we obtain the scaling exponents m κ = m = −2.

Collective scattering: non-perturbative many-body T -matrix
The above perturbative results are in general applicable to the UV range of momenta.However, scaling behavior in the far IR regime, where the momentum occupation numbers grow large, requires an approach beyond the leading-order perturbative approximation as contributions to the scattering integral of order higher than g 2 (i.e., collective phenomena) are no longer negligible.
In order to correctly describe the infrared physics, one therefore has to take into account scattering collective effects.The latter can be achieved by performing a non-perturbative schannel loop re-summation, which is typically derived within the two-particle irreducible (2PI) effective action formalism [116].The re-summation procedure is schematically depicted in Fig. 3.For an N-component field subject to a U(N)symmetric interaction term ∼ gρ 2 /2 in the Lagrangian, depending on the total density ρ = N a=1 ψ † a ψ a , it is equivalent to a large-N approximation at next-to-leading order.As we will demonstrate in Sect.III B, it reflects that also the non-linear , with non-condensed particle density ρ nc .Note that p Ξ sets the scale separating the perturbative region at large p from the non-perturbative collective-scattering region within which the coupling assumes the form (16). Figure taken from [96].
term in the corresponding field equation, cf.(1) for N = 1, depends only on the total density and thus suppresses density fluctuations while the single-component densities ρ a are free to fluctuate.Irrespective of the actual value of N we can use this re-summation scheme to calculate an effective momentumdependent coupling constant g eff (k) that replaces the bare coupling g. (Hence, we neglect, for the first, the conditions for the appropriateness of the chosen approximation.)This effective coupling depends on the distributions n Q (t, k) and thus on momentum, and therefore changes the scaling exponent m of the T -matrix within the IR regime of momenta.In particular, g eff (k) becomes suppressed in the IR to below its bare value g.This ultimately leads to different temporal and spatial scaling of the (quasi)particle spectrum.
For free particles (z = 2) in d = 3 dimensions one obtains [96] where ω k − ω r and k − r are the energy (ω k = k 2 /2M) and momentum transferred in a 2 ↔ 2 scattering process, respectively, and the function on the right is symmetrized, by the permutation terms, in the momenta.At large momenta, the effective coupling is constant and agrees with the perturbative result, i.e., one finds g eff = g.However, below the characteristic momentum scale k Ξ = 2gρ nc M, the effective coupling deviates from the bare coupling g.Here, ρ nc = ρ tot − ρ 0 denotes the non-condensed particle density.Within a momentum range of k Λ k k Ξ , the effective coupling is found to assume the universal scaling form independent of both, the microscopic interaction constant g, and the particular value of the scaling exponent κ of n Q .Below the IR cutoff, i.e., for momenta k < k Λ , the effective coupling becomes constant again, see Fig. 4. Making use of the scaling properties of the effective coupling, we obtain γ κ = 0 in the perturbative regime and γ κ = 2 in the collective-scattering regime for free particles with z = 2.
Together with (15) this yields the corresponding scaling exponent of the T -matrix to be m κ = 2.The same analysis of the effective coupling can be performed for the Bogoliubov dispersion with z = 1.In contrast to free particles the scaling exponent of the T -matrix reads m κ = 0, see [96] for details.

Scaling analysis of the kinetic equation
To quantify the momentum exponent κ, cf.(11), leading to a bi-directional scaling evolution, we study the scaling of the quasiparticle distribution at a fixed evolution time.As the density of quasiparticles, and the energy density, are physical observables, they must be finite.Let us assume that the momentum distribution is isotropic, n Q (k) ≡ n Q (k), and obeys bare power-law scaling n Q ∼ k −κ .The exponent κ then determines whether the IR or the UV regime dominates quasiparticle and energy densities.For a bi-directional selfsimilar evolution the quasiparticle density has to dominate the IR and the energy density the UV.As briefly discussed in the introduction, this is possible within a window of exponents and κ is either the same or different in the IR and UV regions, the latter case being depicted in Fig. 2. Note that, also as introduced before, for ρ Q and Q to be finite, the quasiparticle distribution requires regularizations in the IR and the UV limits, in terms of k Λ and k λ , respectively.According to the scaling hypothesis the time evolution of the quasiparticle distribution is captured by (8), with universal scaling exponents α and β.Global conservation laws strongly constrain the form of the correlations in the system and the ensuing dynamics and thus play a crucial role for the possible scaling phenomena as they imply scaling relations between the exponents α and β.Conservation of the total quasiparticle density (18) Analogously, if the dynamics conserves the energy density (19), the relation has to be fulfilled.The scaling relations ( 21) and ( 22) cannot both be satisfied at the same time for nonzero α and β if z 0. This leaves us with two possibilities: Either α = β = 0, or the scaling hypothesis (8) has to be extended to allow for different rescalings of the IR and the UV parts of the scaling function.In the following, we denote IR exponents with α, β and UV exponents with α , β , respectively.Making use of the global conservation laws as well as of the power-law scaling of the quasiparticle distribution, n Q ∼ k −κ , one finds the scaling relations This implies ββ ≤ 0, i.e., the IR and UV scales, k Λ and k λ , rescale in opposite directions.We remark that these relations hold in the limit of a large scaling spectral region, i.e., for k Λ k λ .Note that energy conservation only affects the UV shift with exponent β , (23b), while particle conservation gives the relation (23a) for the exponent β in the IR.
With this at hand we are finally able to derive analytical expressions for the scaling exponents based on the kinetic theory approach.Performing the s-channel loop-re-summation, the non-perturbative effective coupling g eff can be derived, which depends on the occupation numbers itself, see Fig. 4 and [96] for details of the calculation.The aforementioned anomalous dimension η appears as a scaling dimension of the spectral function and takes into account the possibility to have more involved spectral distributions ρ(ω, k) than the mentioned deltafunction type of free quasiparticles.
As a result, one finds the general scaling relations for the (quasi)particle distributions, where ζ = κ − η + 2 − z.To show possible differences in the scaling behavior of the particle and quasiparticle distributions we added the relations for the particle distribution which scales as n(k) ∼ k z−2+η n Q (k) relative to the quasiparticle number, see the beginning of Sect.III A 1. Note that the momentum scaling of n(k) is characterized by the scaling exponent ζ according to (24d).
From a scaling analysis of the quantum Boltzmann equation, Eqs. ( 8) and ( 9), one obtains the scaling relation Employing the scaling properties of the T -matrix within the different momentum regimes, together with the global conservation laws of the system, one finds the scaling exponents by means of simple power counting to be We remark that the exponents (26b) are usually not observed as the UV region is dominated by a near-thermalized tail.
Cf. also Table II in [96] for a more general account of exponents in the cases of strong and weak wave turbulence.
The above analytic predictions are backed by various numerical results obtained previously and thereafter.The IR scaling exponent β = 1/z has been proposed based on numerical simulations in [117], and was assumed in [84,103].For a single-component Bose gas in d = 3 dimensions, the exponents governing the IR spatio-temporal scaling have been numerically determined to be α = 1.66 (12), β = 0.55(3), in agreement with the analytically predicted values [94].
During the early-time evolution after a strong cooling quench, an exponent ζ d + 1 was seen in semi-classical simulations for d = 3 in [94] and [113], for d = 2 in [108], and for d = 1 in [104,118].In numerical simulations, one has often observed the exponent to be close to d + 2 rather than d + 1, cf., e.g., [99,108].This, however, is due to vortex defects dominating the scaling, which, for low N, is the case in d = 2 and 3 spatial dimensions, as was discussed in [108].
In these studies, a power-law fall-off of the number distribution with ζ = d + 1 was observed in the compressible component only, viz., as soon as the incompressible component had become subdominant following the self-annihilation of the last vortex pair or ring.Also numerical implementations of the full kinetic equation, in d = 3 dimensions, resulted in ζ 4, see [119].
For the Bose gas, the exponents stated above are expected to be valid in d = 3 dimensions as well as in d = 2.The one-dimensional case is rather different due to kinematic constraints on elastic 2 ↔ 2 scattering from energy and particlenumber conservation, which require a more careful analysis, but do not necessarily exclude the predictions to apply.
In the numerical section below, we will demonstrate scaling near non-thermal fixed points in various settings, which go beyond the above analytical approach, as there, the dynamics will be strongly influenced by the appearance of nonlinear and topological excitations.Such excitations have not been taken into account in the basic analytic approach presented above.It is expected, though, that effective field theories can be formulated and analysed along similar lines, that have the potential to describe the scaling under the influence of such excitations.A first example has recently been proposed for the sine-Gordon model [106,107].Using a nonperturbative field-theoretic approach similar to the one summarized above, scaling exponents were predicted for different non-thermal fixed points of the sine-Gordon model [106].This comprises anomalous scaling with β = 1/(d + 2), α = dβ, and κ = 1/(2d + 2), values, which have been corroborated, within varying agreement in d = 2 and d = 3 dimensions, by simulations of a non-linear Schrödinger equation with Besselfunction non-linearity, obtained as the non-relativistic limit of the sine-Gordon equation of motion [107].

B. Low-energy effective field theory
While, in the previous section, collective phenomena that modify the properties of the scattering matrix were taken into account by means of a coupling re-summation scheme, alternative approaches are also available.For example, one can first reformulate the theory in terms of the relevant degrees of freedom, such that the resulting description becomes more easy to treat in the region relevant for the universal dynamics.Given that this mainly affects the low-momentum scales, it is suggestive to employ a low-energy effective field theory approach [120,121].Generally, this requires the key degrees of freedom to be identified, that describe the physics under consideration.In the following, we will briefly outline how this idea can be implemented to describe non-thermal fixed points in a quenched U(N)-symmetric multicomponent Bose gas with quartic interactions.For more details, see Ref. [113].
The crucial observation is that, similar to the singlecomponent Gross-Pitaevskii (GP) model ( 1), its Ncomponent generalisation defines a separation of energy scales between the collective modes and the free particle excitations.Upon adopting a density-phase representation of the field, Φ a = √ ρ a exp{iθ a }, the classical equation of motion reveals that, at low momenta, density fluctuations δρ a = ρ a −ρ (0) a around a mean density ρ (0)  a are suppressed by a factor of ∼ |k|/k Ξ compared to phase fluctuations θ a (around a constant background phase).Here, k Ξ = [2Mρ (0) g] 1/2 is the healing-length momentum scale associated with the total density ρ (0) = a ρ (0) a .The density fluctuations can therefore be integrated out yielding a low-energy effective action S eff [θ] of the model, which depends on the phase degrees of freedom only.This approximation represents a non-linear generalization of the (Tomonaga-)Luttinger Bose liquid [12,13,122].Furthermore, the system provides two types of low-energy modes: N − 1 Goldstone excitations with a quadratic freeparticle-like dispersion ω 1 (k) = ... = ω N−1 (k) = k 2 /2M, corresponding to relative phases between different components, and a single Bogoliubov quasiparticle mode with ω N (k) = [k 2 /2M(k 2 /2M + 2gρ (0) )] 1/2 related to the total phase.This suggests that the physics below the scale k Ξ is well described by the dynamics of gap-less quasiparticles, albeit of two different types.Whereas the single sound mode varies the total density, the N − 1 free quasiparticles represent relative density fluctuations, which locally redistribute the particles in the different components, while keeping the total density constant.We will re-encounter similar excitations in our discussion of a spin-1, three-component system in Sect.IV B.
The resulting low-energy effective action capturing these quasiparticles turns out to contain interaction terms with momentum-dependent couplings, which is in contrast to the coupling constant g in the underlying GP model (1).This in-dicates that the resulting model is non-local in nature, as is commonly expected for an effective theory [113].
Moreover, taking the large-N limit, this action becomes diagonal in component space up to O(1/N) corrections and thus breaks up into N independent replicas.This means that the phases θ a of the different components decouple in the limit of large N. Taking the limit N → ∞, the Bogoliubov mode is no longer present, which suggests that the relative phases are dominating the dynamics of the system, governing the spatial redistribution of relative particle densities between the components, which is not energetically suppressed by the interactions.The N → ∞ effective action in momentum space reads [113] S Here, C denotes again the Schwinger-Keldysh contour, D ab is a free inverse propagator, and we have introduced a short-hand notation, for the interaction terms.These matrix elements contain the momentum-depending coupling g 1/N (k) = gk 2 /2k 2 Ξ ≡ g G (k)/N, which can be compared with the effective coupling obtained by means of the s-channel re-summation for the GP model, cf.(16).The index G of the coupling refers to the relevant Goldstone excitations in the large-N limit.

Spatio-temporal scaling
To analyze the scaling behavior at a non-thermal fixed point we proceed as in Sect.III A by evaluating the WBE in Eq. (7).Instead of the quasiparticle distribution n Q we consider the distribution of phase-excitation quasiparticles, f a (t, k) = θ a (t, k)θ a (t, −k) , dropping in the following the indices to ease the notation.The scattering integral has two contributions, which arise from 3-and 4-wave interaction terms in the effective action action (27), The form of the 3-and 4-point scattering integrals can be inferred from the effective action to be where the corresponding T -matrices are defined by with interaction couplings Here, 'perm s ' denote permutations of the sets of momentum arguments.The scattering integrals scale, analogously to Eq. ( 9), with exponents where γ = 2(z − 1) is the scaling exponent of the effective coupling g eff (k) = s −γ g eff (sk).We remark that the subscript of the coupling is chosen as a general notation covering both cases, z = 2 as well as z = 1.
Using the scaling relation in (25) one can, in principle, derive a closed system of equations, from which the scaling exponents α and β can be inferred.However, since, for different values of the dimensionality d and the momentum scale of interest, one term in the scattering integral can dominate over the other one, it is more reasonable to analyze them independently.
To close the system of equations, an additional relation is required, which is provided either by quasiparticle number conservation, (21), or energy conservation, (22), within the scaling regime.Taking these constraints into account we obtain In the large-N limit (z = 2, γ = 2), the resulting scaling exponents read for both, 3-and 4-point vertices, and for the 4-point vertex, while, at the same time, for the 3-point vertex, no valid solution exists [113].We point out that the above exponents are equivalent to the respective exponents derived in the large-N re-summed kinetic theory for the fundamental Bose fields, for the case of a dynamical exponent z = 2, and a vanishing anomalous dimension η = 0, cf.Sect.III A 4.
One can ask whether both 3-and 4-wave interactions are equally relevant.To answer this question, a comparison of the spatio-temporal scaling properties of the scattering integrals, for a given fixed-point solution f (t, k), is required.Focusing on the conserved IR transport of quasiparticles, for which α = dβ, we obtain In the large-N limit, for which z = 2, one finds µ 3 = µ 4 .
Hence, the relative importance of the scattering integrals I 3 and I 4 should remain throughout the evolution of the system.

Scaling solution
In the remainder of this section we briefly discuss the purely spatial momentum scaling.The scaling of the QBE at a fixed evolution time t = t 0 implies κ = −µ κ,l , where µ κ,l is the spatial scaling exponent of the corresponding scattering integral, I l (t 0 , k) = s −µ κ,l I l (t 0 , sk).Power-counting of the scattering integrals, together with the above stated scaling relation, gives For a given κ l , and assuming the large-N limit (z = 2 and γ = 2), one finds that Hence, the 4-wave scattering integral is expected to dominate at small momenta, k → 0. This implies that, at the nonthermal fixed point, the quasiparticle distribution f (t, k) ∼ k −κ is characterized by the momentum scaling exponent κ = κ 4 = d + 1.This result appears to contradict the previous analysis of the spatio-temporal scaling, which, in the large-N limit, showed equal importance of I 3 and I 4 .We emphasize, however, that the scaling exponents α and β corresponding to the spatio-temporal scaling properties are obtained from relations, which are independent of the precise form of f (t, k) but only require the scaling relation Hence, the questions which vertex is responsible for the shape of the scaling function and which of the vertices dominates the transport can be answered independently of each other.See Ref. [113] for a detailed discussion of this point.8) and g = 0.56 (8).The inset shows the unscaled spectra at five di↵erent times.(b) Rescaled occupation number spectrum for an initial condition with a checker-board of 16 ⇥ 16 vortices with alternating winding numbers w = ±6.The slower evolution of the spectra is reflected in the distinctly smaller scaling exponents ↵ a = 0.40( 5) and a = 0.19 (5).In both panels, the solid black line indicates the scaling function (47), with, in (a), ⇣ g = 4.0(1), and (b), ⇣ a = 5.7 (3).
and ↵ a = 0.40 (5), again reflecting particle conservation, and a steeper fall-o↵ with exponent ⇣ a = 5.7(3).For both fixed points, the scaling of the spectra is a manifestation of the time evolution of the mean separation scale of defects (vortices) in the system.This is confirmed by investigating the mean defect distance `d(t), as seen in Fig. 6.The blue triangles are obtained by averaging the separations of defects, which grows larger in time as the ensemble dilutes.One clearly sees the t 1/2 power law reflected by the separation as well as by the spectra.Interestingly, one can also observe flows of the system from the anomalous fixed point to the Gaussian fixed point (green squares in Fig. 6).Clustering of vortices leads to the initial slow evolution of the system with ' 1/5, yet as the clusters decompose, they e↵ectively behave as a randomly distributed vortex ensembles of ⇡ 1500 elementary vortices, which eventually coarsen with ' 1/2.
We finally emphasise that vortex-anti-vortex annihilation was studied experimentally in a quasi-two-dimensional trapping potential, following the excitation of the system by means of a laser comb pulled through the disc-shaped Bose condensate of 87 Rb atoms [51].Vortices and anti-vortices were tracked separately, and the evolution of their mean distance corroborated the predicted scalings with both, ' 1/2 and ' 1/5.

B. Real-time instantons
A di↵erent system exhibiting self-similar scaling far from equilibrium is the spin-1 Bose gas in d = 1 spatial dimension [49,105], modeled by the Hamiltonian where = ( 1 , 0 , 1 ) T is the three-component bosonic spinor field representing the magnetic sub-levels m F = 0, ±1 of the F = 1 hyperfine manifold and M is the atom mass.q denotes the quadratic Zeeman field strength, which shifts the energies of the m F = ±1 components relative to the m F = 0 component.The term c 0 ⇢ 2 encompasses densitydensity interactions, where ⇢ = † • is the total density.Spin changing collisions are described by the term c 1 |F | 2 , with F = † f and f = ( f x , f y , f z ) being the generators of the so(3) Lie algebra in the three-dimensional fundamental representation.The Hamiltonian of the system is SO(3)⇥ U(1) or, for q , 0, SO(2) f z ⇥ U(1) symmetric.The mean-field phase diagram of the spinor gas spanned in the c 1 -q plane admits various distinct ground states.To prepare the system far from equilibrium, we quench q, such that the system crosses the second-order quantum phase transition line from the polar phase (c 1 < 0, q > 2⇢|c 1 |), showing no magnetization, to the easy-plane phase (c 1 < 0, 0 < q < 2⇢|c 1 |), in which the full SO(2) f z ⇥ U(1) symmetry is broken and which, in the ground state, exhibits magnetization in the F x -F y plane.
Following the quench, the system attempts adjusting to a new ground state, and instabilities form in the Bogoliubov spin eigenmodes of the complex fields m = | m | exp(i' m ), which in this case excite the transverse spin degree of freedom F ? ⌘ F x + iF y = |F ?| exp(i' L ), giving rise to struc-FIG.5. Time evolution of the occupation number spectrum n(t, k) = ψ * (k)ψ(k) .Lengths are measured units of healing length ξ h = (2Mρg) −1/2 with homogeneous density ρ ≡ |ψ(x)| 2 , and time in units of the corresponding interaction time t h = 2Mξ 2 h = √ h /c s , with speed of sound c s = (ρg/M) 1/2 .(a) Rescaled occupation number spectrum for an initial condition with N d = 2400 randomly distributed elementary vortices with winding numbers w = ±1.The rescaling is done with α g = 1.10(8) and β g = 0.56 (8).The inset shows the unscaled spectra at five different times.(b) Rescaled occupation number spectrum for an initial condition with a checker-board of 16 × 16 vortices with alternating winding numbers w = ±6.The slower evolution of the spectra is reflected in the distinctly smaller scaling exponents α a = 0.40 (5) and β a = 0.19 (5).In both panels, the solid black line indicates the scaling function (47), with, in (a), ζ g = 4.0(1), and (b), ζ a = 5.7(3).Figures taken from [103].

IV. NUMERICAL ANALYSIS OF NON-THERMAL FIXED POINTS
In this section, we present numerical simulations of dilute Bose gases prepared in far-from equilibrium initial states, and discuss the ensuing dynamics leading to non-thermal fixed points.
The theory of phase ordering kinetics deals with the relaxation of systems out of equilibrium into an ordered phase.Universal scaling of the system in time and space is associated with non-linear and topological excitations, which introduce time varying length scales into the system, growing as Λ (t) ∼ k Λ (t) −1 ∼ t β , with the universal exponent β.We simulate the dynamics of these gases using the semi-classical truncated Wigner approximation (TWA) [123], that is valid since the systems are in a regime of highly occupied modes.To this end, we consider the classical equations of motions of the respective system, given, e.g., for a one-component dilute Bose gas, by the Gross-Pitaevskii equation (1).To recover beyond-mean-field dynamics, we introduce noise in the Bogoliubov modes to the initial condition and propagate (1) across many noise realizations and average over them.Technically, the propagation is done by means of a pseudo-spectral split-step Fourier method, which ensures the conservation of crucial quantities such as particle number and energy.
In the following, we will present the examples of two systems exhibiting three distinct non-thermal fixed points.First, we will discuss the one-component Bose gas with unstable topological vortices written into the initial condition.We find two different non-thermal fixed points, which are set apart by distinct preparations and decay dynamics of the vortex ensemble [103].Subsequently, we illustrate the phenomenology of a non-thermal fixed point arising after a parameter quench in a spin-1 Bose gas, which excites topological defects in spin space, governing a characteristic length scale Λ of the system to grow algebraically in time.

A. Gaussian and anomalous fixed points in vortex gases
When studying the effects of non-thermal fixed points in the far-from-equilibrium dynamics of cold Bose gases, one typically simulates the time evolution following a strong cooling quench, i.e., from an initial condition, in which the momentum modes are equally occupied up to a maximum cut-off scale Q, recall Eq. (2).The corresponding complex phases of the field modes are chosen at random in each mode.As a result of the a strong quench, the short-time dynamics is characterized by the scattering of macroscopically occupied modes.At later times, strong phase and density fluctuations grow due to non-linear interactions, leading to shock waves, which are giving way to phase gradients forming vortices and anti-vortices.As a result, a length scale is introduced into the system via the mean separation of topological defects, which grows larger in time as vortices and anti-vortices annihilate in a pairwise manner and the defect ensemble dilutes.
Although a strong cooling quench is generically found to lead to a non-thermal fixed point, we wish to further our understanding of the effect of the vortex ensemble on the selfsimilar scaling of the system.Hence, we initialize the system with vortices in it, allowing us to maximize our control over the parameters such as: number of vortices, winding numbers and geometric distribution.The system is thus prepared as a homogeneous, fully phase-coherent state with quantum fluc- tuations included in the empty modes.Interestingly, one finds that, depending on the manner of preparation of the vortex initial condition, two distinct non-thermal fixed points can be observed.
For the first exemplary initial condition, leading to a socalled (near) Gaussian fixed point, the phase of the gas is imprinted with N d = 2400 elementary vortices with winding numbers w = ±1 in a spatially random manner.The ensuing dynamics show the dilution of defects, as the mean separation scale of the system grows with the annihilation of vortex-antivortex pairs [124].This is reflected in the momentum-space field-correlation function of the system, i.e., the occupation number spectrum where k = |k|, t ref is a reference time, f is a universal scaling function, and α and β are the universal scaling exponents, which for reasons of particle number conservation (U(1)-symmetry) are related by α = dβ in d spatial dimensions.Fig. 5a illustrates this scaling evolution.The scaling function depends on the scalar momentum modulus only and takes the form where the constants A and k Λ ∼ −1 Λ are evaluated, in line with Eq. ( 46), at t = t ref .The analytical predictions for a U(1) model with vortices [99,100] are corroborated by the Cf. the experimental results reported in [46].Figure from [103].
The second initial condition, leading to a so-called anomalous non-thermal fixed point, the existence thereof going beyond the analytical predictions, is obtained by imprinting an initial checker-board lattice of vortices with alternating winding numbers w = ±6, as seen in Fig. 6a.
As vortices with winding numbers |w| > 1 are unstable, they quickly decompose into elementary vortices.During the subsequent turbulent evolution, they are observed to form clusters of vortices of either circulation, such that they tend to screen each other.They thus combine to larger eddies and give rise to a quasi-classical turbulent flow [124].As a result, the dipole-pair formation and mutual annihilation of vortices and anti-vortices becomes strongly suppressed.It was shown in [103] that this slowed evolution can be modelled by assuming the vortices to decay predominantly via three-body collisions.The vortex dynamics results in a considerably slowed spatio-temporal rescaling of the correlations as compared to the above near the Gaussian fixed point.As seen in Fig. 5b, the spectra proceed to scale self-similarly with the same kind of universal scaling function, yet with exponents β a = 0.19 (5) and α a = 0.40 (5), again reflecting particle conservation, and showing a steeper fall-off with exponent ζ a = 5.7 (3).
For both fixed points, the scaling of the spectra is a manifestation of the time evolution of the mean separation scale of defects (vortices) in the system.This is confirmed by investigating the mean defect distance d (t), as seen in Fig. 7.The blue triangles are obtained by averaging the separation of  defects, which grows larger in time as the ensemble dilutes.One clearly sees the t 1/2 power law reflected by the separation as well as by the spectra.Interestingly, one can also observe flows of the system from the anomalous fixed point to the Gaussian fixed point (green squares in Fig. 7).Clustering of vortices leads to the initial slow evolution of the system with β 1/5, yet as the clusters decompose, they effectively behave as a randomly distributed vortex ensembles of ≈ 1500 elementary vortices, which eventually coarsen with β 1/2.
We finally emphasise that vortex-anti-vortex annihilation was studied experimentally in a quasi-two-dimensional trapping potential, following the excitation of the system by means of a laser comb pulled through the disc-shaped Bose condensate of 87 Rb atoms [46].Vortices and anti-vortices were tracked separately, and the evolution of their mean distance corroborated the predicted scalings with both, β 1/2 and β 1/5.

B. Real-time instantons
A different system exhibiting self-similar scaling far from equilibrium is the spin-1 Bose gas in d = 1 spatial dimension [48,104], modeled by the Hamiltonian where Ψ = (Ψ 1 , Ψ 0 , Ψ −1 ) T is the three-component bosonic spinor field representing the magnetic sub-levels m F = 0, ±1 of the F = 1 hyperfine manifold, and M is the atom mass.q denotes the quadratic Zeeman field strength, which shifts the energies of the m F = ±1 components relative to the m F = 0 component.The term c 0 ρ 2 encompasses density-density interactions, where ρ = Ψ † • Ψ is the total density.Spin changing collisions are described by the term ) being the generators of the so(3) Lie algebra in the three-dimensional fundamental representation.The Hamiltonian of the system is SO(3)× U(1) or, for q 0, SO(2) f z × U(1) symmetric.The mean-field phase diagram of the spinor gas spanned in the c 1 -q plane admits various distinct ground states.To prepare the system far from equilibrium, we quench q, such that the system crosses the second-order quantum phase transition line, from the polar phase (c 1 < 0, q > 2ρ|c 1 |), showing no magnetization, to the easy-plane phase (c 1 < 0, 0 < q < 2ρ|c 1 |), in which the full SO(2) f z × U(1) symmetry is broken and which, in the ground state, exhibits magnetization in the F x -F y -plane.Following the quench, the system attempts adjusting to a new ground state, and instabilities form in the Bogoliubov spin eigenmodes of the complex fields Ψ m = |Ψ m | exp(iϕ m ), which in this case excite the transverse spin degree of freedom F ⊥ ≡ F x + iF y = |F ⊥ | exp(iϕ L ), giving rise to structure formation in the so-called Larmor phase, ϕ L = ϕ 1 − ϕ −1 .During its relaxation towards equilibrium, the system develops patches of approximately constant Larmor phase, which coarsen in time (cf.Figs.8a, b).This behavior is reflected in the self-similar scaling of the transverse-spin structure factor S F ⊥ (t, k) = F ⊥ (t, k) † F ⊥ (t, k) , which takes on the form with the universal scaling function f s , reference time t ref , and universal scaling exponents α = dβ (see Fig. 8c).Our simulations find the universal function to be once more given by 0 100 200 300 400 500 , with ζ 2 and scaling exponents α β 1/4, which so far is beyond analytical predictions.
In analogy to the coarsening evolution of the vortex gas, we identify a characteristic length scale of the system by studying its topology.The extended dimensionality of the system, due to its multi-component structure, does not allow for stable topological solutions of the complex field F ⊥ as it is the case for density solitons in single-component gases in d = 1 spatial dimension.Nevertheless, the broken SO(2) f z symmetry gives rise to a non-trivial homotopy group in spin space, π 1 (S 1 ⊥ ) = Z, where S 1 ⊥ is to be understood as the unit circle in the F x -F y -plane.Hence, a length scale is introduced into the system via rare topological configurations interpolating between states of constant winding number, where L is the linear length of the system.We refer to such an event, in which the system exhibits an integer jump in Q w , as a real-time instanton.The real-time instantons manifest themselves, in the condensate, as space-time vortices, as can be seen in Figs.9a, b.Each instanton carries a charge, reflecting the integer by which the winding number jumps, as well as a topological current, j µ = ∂ µ ϕ L , which we can utilize to compute the spatio-temporal probability distribution function (PDF) P(r, t) of the instantons.The PDF decays exponentially with defect separation r, P(r, t) ∼ exp[−r/r Λ (t)], with a time varying mean separation scale r Λ (t) (cf.Fig. 9c).The mean separation scale r Λ (t) is growing algebraically, with a power law r Λ (t) ∼ t β I , where β I = 0.27(1) (see inset of Fig. 9c), which is in agreement with the self-similar scaling of the order-parameter spectrum.
V. THEORY VS.EXPERIMENT In this final section, we give a short overview of the theory development of non-thermal fixed points and briefly discuss four experiments with ultracold atomic gases, which have explored different aspects of universal dynamics close to a nonthermal fixed point.
The existence and significance of strongly non-thermal momentum power-laws, requiring a non-perturbative description reminiscent of wave turbulence, was originally proposed in the context of reheating after early-universe inflation [90,91], then later generalized to scenarios of strong matter wave turbulence in non-relativistic systems [92,125], in particular ultracold superfluids and, in their context, to the dynamics of topological defect ensembles [98-100, 108, 118, 126, 127], see also [101-103, 105, 128-131].
Universal scaling at a non-thermal fixed point in both space and time was numerically observed as algebraic time evolution of the correlation length and the condensate fraction akin to coarsening [100,108,127,129] and formalized by means of the spatio-temporal scaling form (3) for the occupationnumber distribution, for both, non-relativistic and relativistic models [52, 94, 96, 103-107, 113, 132-135], see also [112,130].It has direct applications in the context of relaxation and plasma formation in heavy-ion collisions [5,[136][137][138][139] as well as for axionic models relevant in cosmology [130,132].For previous overview articles, see [95,96,114].
The existence of non-thermal fixed points was experimentally observed in [48,49].In [48], a spinor Bose gas of 87 Rb atoms (see also Sect.IV B) was prepared in a condensate state in the polar phase, where all atoms are in the m F = 0 hyperfine component.By changing suddenly a quadratic Zeeman shift, it was quenched into the easy-plane phase, where the system wants to develop a non-vanishing angular momen-tum F in the F x -F y plane perpendicular to the direction of the Zeeman splitting.The quench thus leads to an instability, which quickly gives rise to excitations of the form of the red dashed line in Fig. 2, in the spin excitations.The subsequent self-similar evolution in the quasi one-dimensional spin excitations was characterised by α = 0.54 (6), β = 0.33 (8), and ζ ≈ 2.6.
In [49], a single-component 87 Rb Bose-condensate was quench-cooled into a quasi-one-dimensional cigar-shaped trapping potential on the surface of a microchip.As a result, strong longitudinal excitations built up in the system which gave rise to a momentum distribution resembling an ensemble of solitons [118].The self-similar scaling of the momentum distribution function was found to be anomalously strongly slowed, with α = 0.09(3), β = 0.10(3), and ζ = 2.39 (18).An extended exponential tail demonstrated the presence of a dense ensemble of solitons.
In the experiment [46], a grid of elliptical obstacles was dragged through a uniform planar 87 Rb condensate, which gave rise to the excitation of many vortices and anti-vortices.The setup served to demonstrate the buildup of so-called Onsager clusters of many elementary vortices of equal circulation.As is described in more detail in Sect.IV A, such clusters can shield vortex-anti-vortex pairs from mutually annihilating and lead to universal scaling dynamics with anomalously small exponents α and β.In the experiment, such scaling was observed in the time-evolution of the characteristic length scale measuring the mean distance between vortices.The results thus corroborated the values α = d/5 and β = 1/5 in the d = 2 dimensional dynamics, predicted in [103].
The experiment reported in [47] explored the bi-directional transport predicted in the universal dynamics as sketched in Fig. 2. The initial quench removed 77% of the atoms and 97.5% of the energy from the 39 K condensate in a cylinder trap, by turning off the interactions and lowering the trap edge for a brief amount of time.In the ensuing re-equilibration of the quench-cooled gas, both, the inverse particle, and the direct energy flow were observed.Measurements of the scaling exponents gave α = 1.08(9) and β = 0.34(4), cf.Eqs.(26a) for the IR particle flow, while the direct UV energy flow was characterized by α = −0.70(7) and β = −0.14(2),cf. the discussion in Sect.III A, in particular Eqs.(26a), (26b).Both values, β and β deviate weakly from the predictions (for η = 0), which may be explained, in the IR, by the system still being in a prescaling regime [52,140], where the exponents are slowly increasing in time.
In the experiment [50], a Bose condensate of 87 Rb atoms was driven out of equilibrium by imposing a small rotational oscillation onto the elongated quadrupole-Ioffe configuration trap.In the ensuing evolution of the momentum distribution self-similar motion towards higher momenta was observed, with α = −0.50(8), β = −0.2(4).The relation between the exponents is consistent with the prediction α/β = d for number conservation in a two-dimensional situation, which here applies to the projected distributions extracted from the data.
For the recently published results of a further experiment on a ferromagnetic spinor Bose gas, see [51].

VI. OUTLOOK
In this brief tutorial review, we have discussed the nonequilibrium phenomenon of universal scaling dynamics in strongly quenched quantum many-body systems.We introduced to the concept of non-thermal fixed points and summarized the main ideas of analytical approaches to describing the scaling behavior from first principles.
This comprises a brief outline of the 2PI formalism for obtaining a non-perturbative kinetic-theory formulation of nonthermal fixed points.Scaling exponents can be determined by power counting, assuming a pure scaling form to solve the dynamic equations for non-equilibrium two-point correlation functions such as time-evolving mode occupancies.An alternative, low-energy effective field theory description of U(N) models allows predicting the universal scaling behavior on the grounds of a perturbatively coupled Luttinger Bose liquid in the large-N regime.This entails, in particular, the scaling exponents α and β, characterizing the time evolution of the system in the vicinity of the non-thermal fixed point, and the exponent ζ defining the algebraic fall-off of the momentumspace scaling function.
At this point, let us briefly mention a novel analytical approach based on the correspondence between scaling and fixed points of the renormalization group [141,142].The crucial observation is that all the universal scaling properties can be extracted from the vicinity of a given infrared fixed point.It is therefore suggestive to try to extend this idea to the case of far-from-equilibrium self-similar dynamics and to demonstrate how non-thermal fixed points can be understood from the renormalization-group perspective.To a certain degree, this goal has been already achieved for the case of stationary (strongly) non-equilibrium configurations, see, e.g., [91,125,143].In addition, an alternative renormalization scheme involving a temporal regulator has been proposed as a suitable description of far-from-equilibrium systems even beyond the stationary case [144][145][146].However, a complete satisfactory renormalization-group description of non-thermal fixed points is still lacking.
The first attempt to implement this program, within the functional renormalization group (fRG) framework [147][148][149][150][151][152], has been made in [141,142], for the specific example of a single-component Bose gas.The employed method follows closely the works [153] and [125], in which the fRG fixed-point equations were used to analyze infrared scaling properties in Landau gauge QCD and the stochastic drivendissipative Burgers' equation, respectively.The central object in this approach is the flow equation that describes the change of correlation functions under successive application of momentum-shell integrations and thus their dependence on momentum scale [154][155][156].In the vicinity of a fixed point, the two-point functions are then parametrized in terms of the full scaling forms and of the deviations of the two-point correlators from those at vanishing cutoff scale.As we noted above, the universal scaling properties are encoded in (the asymptotic limits of) these deviation functions, determined by the fixed-point equations, which can be obtained upon integrating out the RG flow.These equations can then be solved numeri-cally in the asymptotic limits of interest allowing us to extract the universal exponents associated with far-from-equilibrium scaling dynamics at non-thermal fixed points.
Based on the experimental and numerical results as well as analytical predictions, we can conclude that universal dynamics at or close to non-thermal fixed points emerge in various settings, characterized by different symmetries of the system as well as distinguished by different initial conditions.While the examples we have discussed constitute infrared fixed points, implying that the self-similar evolution comprises transport to lower wave numbers, i.e., larger length scales, also the opposite case of ultraviolet non-thermal fixed points has been considered [53,54,96].In the former case, the respective phenomena have often be characterized as coarsening known as an ordering phenomenon in statistical physics far from equilibrium, the latter is relevant in understanding aspects of scaling in thermalization on microscopic scales such as following heavy-ion collisions.The concept of non-thermal fixed points is to provide a first-principles formulation and classification of such phenomena based on microscopic quantum field models of the respective system and their characteristic symmetry properties.
In the context of the numerical studies, we have sketched the important role of topological defects to the coarsening dynamics of the system, which typically require theoretical techniques beyond perturbative kinetic theory and nonperturbative Feynman diagrammatic methods.Moreover, such defects can show very different collective dynamic behavior and thus signal proximity of the system to different non-thermal fixed points or even a flow from one to another.While coarsening in a single-component Bose condensate in two spatial dimensions, bearing quantum vortices and antivortices, is driven by their mutual annihilation, we demonstrated that the universal infrared scaling of a one-dimensional spinor Bose gas can show related but quite different phenomena.Here, vortices appear as defects in the two-dimensional plane defined by space and evolution time, so-called real-time instantons.
Current efforts to develop a comprehensive understanding of coarsening dynamics far from equilibrium include numerical investigations into the role of disordered driven caustic dynamics, e.g., in the spinor gas showing instanton events [157], which can bear important consequences in various fields of research, e.g., for the study of cosmological structure formation.Further systems include in particular dipolar gases [158,159], which offer the possibility of exploring universal dynamics of systems with strong long-range interactions, which show a richer spectrum of phases already in equilibrium.
In summary, the physics of dynamics far from equilibrium, and in particular its possible universal characteristics, which can relate very different systems with each other, remains an exciting and rich field of fundamental research in quantum many-body physics.

FIG. 3 .
FIG. 3. Graphical representation of the re-summation scheme.(a) The two lowest-order diagrams contributing to the loop expansion of the 2PI effective action that lead to the quantum Boltzmann equation with perturbative T -matrix (13) or (14).Solid lines represent the full Green's function G(x, y), black dots the bare vertex ∼ gδ(x − y).(b) Diagram representing the re-summation approximation which replaces the diagrams in (a) within the IR regime of momenta and gives rise to the modified scaling of the T -matrix.(c) The wiggly line is the effective coupling function entering the T -matrix, which corresponds to a sum of bubble-chain diagrams, here written as an integral equation.Figure taken from [96].

2 p
FIG.4.Effective coupling g eff (p 0 , p) in d = 3 dimensions as a function of momentum p = |p|.The figure shows cuts in the p 0p-plane, with p 0 = 0.5ε p (dark solid lines) and p 0 = 1.5ε p (transparent solid lines).Different colors refer to different infrared cutoffs p Λ , see the legend.Units are set by the 'healing'-length wave number p Ξ = (2gρ nc m) 1/2 , with non-condensed particle density ρ nc .Note that p Ξ sets the scale separating the perturbative region at large p from the non-perturbative collective-scattering region within which the coupling assumes the form(16).Figure taken from[96].

FIG. 6 .
FIG. 6. Snapshots of the time evolution of the hydrodynamic velocity fields v(x, t), for a vortex-lattice initial condition (units are chosen as in Fig. 5, with M = 1/2).The color encodes the modulus |v| of the field, whereas the black flow lines indicate its orientation.The positions of (anti-)vortices are marked by (green) orange dots.Panel (a) shows the checker-board initial vortex lattice with 16×16 vortices with winding numbers w = ±6.(b)-(d) show snapshots at times t = {300, 10 3 , 10 4 }t h .Figure taken from [103].

6 d ∼ t 1/ 2 d ∼ t 1/ 5 FIG. 7 .
FIG. 7. Time evolution of the mean distance between defects d (t), starting from three different initial vortex configurations.The blue triangles show the evolution from a random distribution of N d = 2400 elementary vortices and anti-vortices at initial time t 0 = 0. Green squares (red circles) correspond to the time evolution from an initial lattice of 16 × 16 and (8 × 8) vortices with winding numbers w = ±6.Different temporal scalings d (t) ∼ t β are observed, including the flow of the system crossing over from the Gaussian non-thermal fixed point (β g 1/2) to the anomalous one (β a 1/5).Cf. the experimental results reported in[46].Figure from[103].

pFIG. 8 .
FIG. 8. Time evolution of the transvere spin F ⊥ = F x + iF y = |F ⊥ | exp(iϕ L ) in a single Truncated Wigner (TW) run.(a) Time evolution of the transverse spin length |F ⊥ |.In the initial state, the system is in the polar phase exhibiting no magnetization.After approximately one spin-changing collision time t s = 2π/ √ n|c 1 |, the system begins to reorder into the new phase and a finite spin-length emerges which fluctuates weakly about a mean value.(b) Time evolution of the Larmor phase ϕ L .Patches of approximately equal phase arise, which grow larger in time.Spin-wave excitations are seen, which propagate with the spin speed of sound c s = (ρ|c 1 |/2M) 1/2 (shown as dashed black and red lines), as well as strong phase kinks.(c) Self-similar scaling of the transverse-spin structure factor S F⊥ (t, k) = |F ⊥ (k)| 2 with universal exponents α β 1/4 and ζ 2.Figure taken from [157].
FIG. 8. Time evolution of the transvere spin F ⊥ = F x + iF y = |F ⊥ | exp(iϕ L ) in a single Truncated Wigner (TW) run.(a) Time evolution of the transverse spin length |F ⊥ |.In the initial state, the system is in the polar phase exhibiting no magnetization.After approximately one spin-changing collision time t s = 2π/ √ n|c 1 |, the system begins to reorder into the new phase and a finite spin-length emerges which fluctuates weakly about a mean value.(b) Time evolution of the Larmor phase ϕ L .Patches of approximately equal phase arise, which grow larger in time.Spin-wave excitations are seen, which propagate with the spin speed of sound c s = (ρ|c 1 |/2M) 1/2 (shown as dashed black and red lines), as well as strong phase kinks.(c) Self-similar scaling of the transverse-spin structure factor S F⊥ (t, k) = |F ⊥ (k)| 2 with universal exponents α β 1/4 and ζ 2.Figure taken from [157].

r Λ ∝ t 0. 27 FIG. 9 .
FIG. 9. Real-time instantons in the time evolution of the Larmor phase.(a) Time evolution of the winding number Q w .Integer-valued jumps are observed, which are caused by the space-time vortices seen in (b).(b) High resolution excerpt of the time evolution of the Larmor phase after a quench.A plaquette algorithm correlating phase jumps and dips in spin length locates space-time vortices, which each correspond to a winding-number jump by ±1.The winding of ϕ L by 2π around the core of the vortex is evident in the magnified section shown in the inset.(c) Probability distribution function (PDF) of defect separation.The probability decays as an exponential function exp(−r/r Λ (t)), with a time-varying mean separation scale r Λ which scales in time according to r Λ (t) ∼ t β I , with β I = 0.27(1).Figure adapted from [157].
Time evolution of the occupation number spectrum n(t, k) = h ⇤ (k) (k)i.Lengths are measured in units of the healing length ⇠ h = (2M⇢g) 1/2 , with homogeneous density ⇢ ⌘ h| (x)| 2 i, and time in units of the corresponding interaction timet h = 2M⇠ 2 h = p 2⇠ h /c s ,with speed of sound c s = (⇢g/M) 1/2 .(a) Rescaled occupation number spectrum for an initial condition with N d = 2400 randomly distributed elementary vortices with winding numbers w = ±1.The rescaling is done with ↵ g = 1.10(