Prescaling relaxation to nonthermal attractors

We study how isotropic and homogeneous far-from-equilibrium quantum systems relax to nonthermal attractors, which are of interest for cold atoms and nuclear collisions. We demonstrate that a first-order ordinary differential equation governs the self-similar approach to nonthermal attractors, i.e., the prescaling. We also show that certain natural scaling-breaking terms induce logarithmically slow corrections that prevent the scaling exponents from reaching the constant values during the system's lifetime. We propose that, analogously to hydrodynamic attractors, the appropriate mathematical structure to describe such dynamics is the transseries. We verify our analytic predictions with state-of-the-art 2PI simulations of the large-N vector model and QCD kinetic theory.

Introduction.-Thermalization of isolated quantum many-body systems is an important contemporary research problem of a broad scope.Its relevance ranges from cold atom systems, through QCD in ultrarelativistic nuclear collisions all the way to gravity and black hole physics [1].Given the complexity of modeling quantum many-body dynamics and the richness of non-equilibrium phenomena, emergent regularities that form a basis for a quantitative understanding are of particular interest.
In this work we are concerned with an important instance of such an emergent regularity: far-fromequilibrium self-similar time evolution of nonthermal attractors, also known as nonthermal fixed points.These phenomena are transient stages in the thermalization dynamics, whose defining feature is self-similar scaling behavior in time.Consider a momentum distribution function f (t, p) of a homogeneous and isotropic system, where t is time and p spatial momentum.The system reaches a nonthermal attractor, when f scales with time in a characteristic momentum range with constant scaling exponents α ∞ and β ∞ .Indeed, such behavior corresponds to a vast reduction in the complexity, as the knowledge of the distribution function at some time allows one to determine the distribution function at a different time by a simple rescaling.Nonthermal attractors appear in the studies of isolated quantum systems across a wide range of energy scales: ultracold quantum gases [2][3][4][5][6][7][8][9][10][11], ultrarelativistic nuclear collisions [12][13][14] and early universe cosmology [15,16].Despite significant interest in nonthermal attractors, a quantitative understanding of how a system approaches a nonthermal fixed point remains elusive [17][18][19][20][21].
In [17] it was proposed that even prior to reaching the nonthermal attractor (1) the system can exhibit prescaling, where f has already assumed the fixed-point shape f S but continues to evolve with time dependent scaling exponents α(t) and β(t) f (t, p) = A(t) f S (B(t)|p|). ( The prescaling factor A(t) = exp[ t t0 dt ′ α(t ′ )/t ′ ] reduces to the fixed-point scaling of Eq. (1) when α(t) approaches α ∞ .The same holds for B(t) in terms of β(t).
Given that scaling (1) is an asymptotic late time statement known to be reached slowly, the systems of interest might in fact spend a much greater fraction of their lifetime prescaling (2) rather than scaling (1).Therefore, a quantitative understanding of prescaling is as important as understanding scaling itself.
In our work, we develop a simple theoretical description of prescaling dynamics that uses the same assumptions as the ones used to derive scaling.We test our predictions using strongly-correlated large-N vector model and weak coupling QCD kinetic theory simulations.
Scaling implies prescaling.-Understandingprescaling requires identifying laws governing time evolution of A(t) and B(t) (or, alternatively, α(t) and β(t)).As we show, these laws have a surprisingly simple origin and form.
The key role in deriving scaling (1) is played by conserved quantities: particle number density n = d d p f /(2π) d or energy density E = d d p ω p f /(2π) d , where d is the number of spatial dimensions and ω p is the dispersion relation of particles.We focus on ω p ∼ |p| z .Applying conservation of n or E in the momentum regime of interest imposes the relation α ∞ = σ β ∞ between the scaling exponents [22].When n = const, then σ = d, while E = const gives σ = (d + z).The conserved quantities are local in time, which means that they in fact constrain also prescaling exponents in exactly the same way: α(t) = σ β(t).Equivalently, A(t) = B(t) σ .This implies that there is only one independent degree of freedom in the isotropic and homogeneous prescaling, which we will choose to be B(t).
The time evolution for the independent prescaling factor B(t) is still subject to the equation of motion for f .In the case of a kinetic theory it is given by the Boltzmann equation with collision kernel In the present section, we assume the collision kernel to be a homogeneous functional of particles momenta, i.e., to simply scale under Eq.( 2) by A(t) µα B(t) µ β ≡ B(t) σµα+µ β for some real numbers µ α,β .This assumption applies to many (but not all) collision kernels describing nonthermal attractors (see the Supplemental Material [23] for explicit examples).Typically overoccupation singles out terms with the highest power of the distribution function and associated matrix elements often happen to scale homogeneously under rescalings of momenta.For such collision kernels, we can separate timedependent and (rescaled) momentum-dependent contributions by substituting the prescaling ansatz (2) in the Boltzmann equation (3) and using where p = B(t) p is the rescaled momentum, 1/β ∞ = (1 − µ α )σ − µ β , and D 1 is a separation of variables constant.The intrinsic time dependence of our setup implies nonzero D 1 , which can be fixed at any time in the prescaling evolution and we choose D 1 = β(t 0 )/t 0 .The idea of separation of variables in the context of nonthermal attractors appeared already in [24], but for homogeneous C[f ] only solutions with constant scaling exponents were considered.Our key observation here is that prescaling is encapsulated by the general solution From Eq. ( 5) it is clear that prescaling induces power-law corrections to scaling.The prescaling originates from the presence of nonzero t * = t 0 (1 − β ∞ /β 0 ), where t 0 and β 0 ≡ β(t 0 ) correspond to the initial data.Its appearance comes as no surprise: the dynamics of the system in question is time translationally invariant and therefore t appearing in formulas needs to be measured with respect to some time, t * .The reason why it does not appear in Eq. ( 1) is because the exact scaling is an asymptotic late time statement and dependence on t * drops.Note that while β ∞ and α ∞ are theory specific and independent of initial conditions, t * will depend on a chosen initial state.Before we move to testing Eq. ( 5) using ab initio solutions of quantum dynamics, let us reiterate that this result originates from the Boltzmann equation and pertinent conservation laws.These are exactly the same constraints as used in a conventional scaling analysis [22].Prescaling in the case of collision kernels being homogeneous functionals of momenta can therefore be understood as a direct consequence of the existence of scaling.
Prescaling in large-N vector model.-Webegin by testing Eq. ( 5) against the full quantum dynamics of a nextto-leading order (NLO) large-N vector model at small coupling λ.This model features a dual cascade [22,25] consisting of an inverse particle cascade in the IR and  a direct energy cascade to the UV.We focus on the IR scaling, which is characterized by effective particle number conservation with (α ∞ , β ∞ ) = (d/2, 1/2) and was also realized in cold quantum gases [2,4].This fixed point arises for large occupations f ∼ 1/λ ≫ 1 such that its description requires going beyond a standard kinetic theory analysis.Large-N kinetic theory addresses this regime due to inclusion of relevant resummations [22,26].The corresponding collision kernel scales homogeneously with µ α = 1 and µ β = −2 [22].We perform ab initio studies of this fixed point in 3 + 1 dimensions using the 2PI formalism following [21] to validate Eq. ( 5) and the underlying assumptions.Below the mass gap the equal-time statistical function F (t, t, |p|) reduces to f (t, |p|) [22], where the system is initialized with f (t 0 , |p|) = f 0 θ(Q − |p|) with f 0 = 100/λ, λ = 0.01 and Q is the characteristic hard scale far from equilibrium.Following [17] we extract the prescaling factors A(t) and B(t) from the time evolution of integral moments of F (t, t, |p|) [23], e.g., B(t) = n(t)E(t 0 )/(E(t)n(t 0 )) [27].In the upper right panel of Fig. 1 we show how rescaling with A(t) and B(t) leads to an early collapse of distributions at different times, while a considerable spread remains when rescaling with the fixed point exponents (left panel).The evolution of the extracted B(t) is shown in the lower panel to be well described by Eq. (5a) already at early times (dashed black line), and only asymptotes to the corresponding fixed point scaling behavior (1) (solid green line).The full NLO quantum dynamics are shown to be captured remarkably well by our effective kinetic description of Eq. ( 5) already from times close to initialization.
Prescaling in isotropic QCD kinetic theory.-Wemove now to studying prescaling dynamics in QCD, whose nonthermal fixed point plays an important role in our understanding of thermalization dynamics in weaklycoupled models of ultrarelativistic nuclear collisions [1].We use QCD kinetic theory, where the evolution of the color and polarization averaged gluon distribution function f (t, p) is described by 2 ↔ 2 and 1 ↔ 2 processes [28]: see [23,[28][29][30].Nonthermal fixed points can be reached from a wide range of initial conditions including large occupation numbers [31][32][33][34][35][36][37], which we implement via Here g 2 is the square of the coupling and n 0 is the initial occupation.We consider n 0 = 1 and g 2 = 10 −8 .To obtain the precise late time behavior we initialize at t i Q = 0 and evolve for very long times until t f Q = 10 8 .Results will be given in units of the characteristic energy scale Q as given by the maximum of |p| 2 f (t i , p).We discuss explicitly only the pure gluon simulations where the scaling phenomenon is encountered after checking that our results do not change under the inclusion of quark/anti-quark dynamics.
A scaling analysis for the vacuum QCD collision kernel together with energy conservation σ = d + z = 4 reveals the direct energy cascade fixed point (α ∞ , β ∞ ) = (−4/7, −1/7), see [38][39][40].However, the overall scaling of the elastic collision kernel is broken by the presence of the Debye mass m the Debye mass decreases over time such that the violation of scaling only leads to a delay in the approach to the fixed point.This (diminishing) breaking of scaling for QCD kinetic theory is demonstrated in Fig. 2, where we extract β(t) (solid color lines) from different moments of the distribution function [23].We do not display α(t) explicitly as we find the scaling relation α(t)/β(t) = 4 realized to a very good accuracy.β(t) is observed to approach the fixed point value β ∞ = −1/7 (dashed gray line) but a finite deviation remains due to m D (t) even after eight orders of magnitude in simulation time.The decreasing but finite spread between different moments further demonstrates that different momenta of the distribution function approach the fixed point on different timescales.A fast convergence of the moments and collapse to β ∞ is in contrast found for simulations without the elastic collision kernel C 2↔2 [23].In the following we will study these deviations analytically in a small-angle scattering approximation of C 2↔2 .
Effect of scaling breaking terms in the Fokker-Planck approximation.-Thebreaking of scaling inhibits prescaling exponents extracted from different moments to share the same universal prescaling dynamics.Nevertheless, at qualitative level the scaling dynamics can be reasonably modeled via the Fokker-Planck (FP) approximation [41][42][43].This approach assumes the dominance of small angle scatterings and has previously been used in the context of nonthermal attractors [44] and prescaling [19,20].We will compare our analytical results from the FP approximation to simulations using the full QCD collision kernel.The corresponding FP collision kernel allows us to factorize the scaling-breaking Coulomb logarithm, which involves the ratio of the UV scale, the characteristic gluon energy ⟨p⟩, and the IR scale, the Debye mass where we only display terms relevant for the scaling analysis [23].We note that the time-dependence due to the Coulomb log can be identified with the background term in the formalism of [24].The FP-kernel does not scale homogeneously and the solution (5) does not apply.However, separation of variables can still be performed as the necessary property is factorization in t and pdependencies.Upon separating variables with associated constant D 2 , see Eq. ( 4), and relating α(t) = σβ(t) as before, we now obtain Equation ( 7) can be directly integrated, but then B(t) appears as an argument of a nontrivial transcendental function.A more useful approach is to derive from Eq. ( 7) a second order differential equation for β(t) which is directly solved by β(t) = β ∞ = −1/7.The solution to Eq. (8) (dashed black line) is shown in Fig. 2 to capture well the evolution of β(t) obtained from solving the full QCD kinetic theory collision kernel from times shortly after initialization of the system over more than eight orders of magnitude.We obtain this result by solving Eq. ( 8) with initial conditions for β0 determined consistently from the full QCD kinetic theory data at t 0 [23].
In the inset we show that solving Eq. ( 8) shortly after initialization describes qualitatively well the late-time dynamics.
With the applicability of Eq. ( 8) demonstrated in Fig. 2, we now use it to study the prescaling dynamics in the vicinity of the fixed point.We can linearize Eq. ( 8) in perturbations δβ(t) around the fixed point value β ∞ = −1/7, which yields a power-law decay from below δβ(t) ∼ −1/t.The full solution to Eq. ( 8) at late times however turns out to be governed by logarithmic corrections not captured by this linearization procedure.We find that a consistent late-time solution to Eq. ( 8) is given by where we used Q as the reference scale but emphasize that the choice of a constant does not matter at late enough times.Similar late-time power-law [20] corrections from linearization and late-time logarithmic [19] corrections induced by the temporal evolution of the Coulomb logarithm were found in the FP approximation for the Baier-Mueller-Schiff-Son [12] fixed point in longitudinal expanding plasma.The simple power-law approach to the fixed point found in the absence of scaling breaking terms in Eq. ( 5b) is therefore enriched to involve both fast (power-law) and slow (inverse powers of logarithms and slower) behavior.This discussion is reminiscent of the transseries form [45,46] for late time dynamics of the energymomentum tensor of matter undergoing longitudinal boost-invariant expansion [47][48][49][50].There slow modes came from relativistic hydrodynamics and exponentially faster modes from transient excitations.Here slow modes come from the Debye mass breaking the homogeneity of the collision kernel with respect to rescalings of momenta and fast modes are the original prescaling excitations encountered already in (5b).Similarly to [47][48][49][50], it is not difficult to gather finite order indications that the series containing slow modes (9a) has a vanishing radius of convergence with β m+1,0 /β m,0 ∼ m [23].Curiously, the leading (at each m) doubly logarithmic term behaves geometrically: β m+1,m /β m,m−1 = −1.It would be interesting to develop systematic understanding of this behavior, including resummations of the resulting transseries.A good starting point might be analysis of the Painlevé I equation in [51,52], which is also second order and exhibits expansion in three building blocks analogous to ours t, log(t Q) and log (log (t Q)).
In Fig. 3 we visualize the attractive nature of the prescaling dynamics by extracting prescaling exponents from different initial conditions.All simulations are initialized with variations in parameters of the class of initial condition used in this work apart from the data rep- resented by light blue, which uses box initial conditions The prescaling exponents extracted from different simulations are all found to converge to a universal late-time behavior, which we additionally show is well described by Eq. (9b) (solid black line).Furthermore, we want to emphasize the similarity between the behavior shown in Fig. 3 and hydrodynamic attractors, where different solutions converge to a single universal curve which at sufficiently late times is described by relativistic hydrodynamics [47,53,54].
The above analysis has an important bearing on the appearance of scaling.The regime when the highest order terms in the collision kernel dominate parametrically ends when the typical occupancy becomes of O(1).This is realized if t f Q ≥ α −7/4 S [13].At that time, we have a deviation of δβ(t)/β ∼ 1/ log(α −7/4 S ) ≲ 0.03 with g 2 = 10 −8 .As a consequence of this, QCD kinetic theory will therefore still show percent deviations from the fixed point values when the direct energy cascade ceases and ultraviolet modes |p|/Q ≥ 1 start to thermalize.
Conclusions.-Westudied the approach of isotropic and spatially homogeneous quantum many-body systems to nonthermal attractors.Our results demonstrate that the prescaling is governed by a simple first-order ordinary differential equation obtained from the underlying dynamics via emergent conservation laws.
Our analytical prediction implies that prescaling entails infinitely many power-law corrections to constant scaling exponents.They conspire to a simple time off-set in the fixed point scaling.We have successfully tested our simple formula for prescaling against ab initio simulations of a relativistic vector model QFT using 2PI formalism and QCD kinetic theory simulations.Our QCD kinetic theory simulations span eight orders of magnitude in time and provide the most accurate extraction of scaling exponents to date.
The exact scaling associated with nonthermal attractors requires the collision term to be a homogeneous functional of particle momenta at large occupations.For QCD kinetic theory this property is violated by the Debye mass term that regulates the Coulomb divergence in the elastic scattering matrix element.We demonstrate that exact scaling exponents are not reached during the lifetime of the system.Using the Fokker-Planck approximation to QCD kinetic theory we show that the scalingbreaking Couloumb logarithm significantly enriches the prescaling dynamics.The late-time behavior is given by a factorial divergent series that includes inverse powers of logarithms and positive powers of double logarithms of time.This constitutes a striking structural similarity with theoretical descriptions of hydrodynamic attractors in the boost-invariant models of nuclear collisions.
Our work shows that prescaling is an unavoidable consequence of nonthermal attractors.Therefore our analytical predictions for prescaling can be verified experimentally in cold atom systems.Furthermore, we uncovered that scaling breaking terms generate rich prescaling dynamics that bares similarities to transseries in the context of hydrodynamic attractors.It would be fascinating to utilize the enormous degree of control in cold atom systems to induce scaling breaking terms and experimentally discover the phenomenology of transseries.
Note added: Several months after the arXiv submission of the present work, Ref. [55] appeared and reported the observation of the prescaling solutions that we derived in Eq. ( 5) in an ultra-cold atom experiment.
We now show that the inelastic collision kernel (A5b) scales with µ α = 3 and µ β = −1 under prescaling in the non-expanding system (2).We consider only the first term since both have the same scaling behavior.Crucially, we need to know the scaling behavior of the splitting rate, which can be extracted from its two prevalent limiting regimes for soft gluon radition z = p ′ /p ≪ 1: the Bethe-Heitler (BH) limit in which interferences between successive scatterings are negligible and the Landau-Pomeranchuk-Migdal (LPM) limit in which successive scattering events by the medium interfere destructively.The rate in the respective limit reads where we only included those terms relevant for a scaling analysis with diffusion coefficient in this overoccupied scenario.For q(µ)| µ=em D in the BH limit, we thus have q(µ) = A(t) 2 B(t) −3 q and For the LPM limit, µ is to next-to-leading-logarithmic order given self-consistently via which shows that µ scales self-consistently as µ 2 ∼ A(t)B(t) −2 such that the time dependence drops out of log(µ 2 /(2m 2 D )).Accordingly, for the LPM limit we thus find again The rate therefore scales like the Debye mass in both limits and we will adopt this scaling for the complete rate γ g gg ∼ A(t)B(t) −2 γg gg .This leads to the overall scaling prediction gg (p; p′ , k′ )B(t)δ (1) (p − p′ − k′ ) × A(t) 2 f S, pp(f S, p′ p + fk′ p) − f S, p′ pf S, k′ p (A14) such that we can identify µ α = 3 and µ β = −1 as anticipated.C 1↔2 will therefore lead to the direct energy cascade fixed point (α ∞ , β ∞ ) = (−4/7, −1/7).
The scaling analysis for the elastic collision kernel in the absence of the Debye mass has been performed in [39] and can simply be generalized to the prescaling case with again µ α = 3 and µ β = −1.This scaling analysis however needs to be augmented by the inclusion of the Debye mass as we discussed above, which leads to an effective regulation of the divergent soft contributions to the elastic collision kernel and prevents one from extracting an overall scaling behavior thereof.

Fokker-Planck approximation
We assume the gluons to interact via elastic small-angle scatterings such that the collision kernel takes a FP form [1,19,41,43] where we highlighted the contributions relevant for a scaling analysis and here ⟨. . .⟩ ≡ ).The advantage of the FP approximation for a prescaling analysis becomes apparent here, as the contribution due to the Debye mass in the elastic QCD scattering matrix element is simply factorized into a logarithm of the characteristic UV and IR scale.Plugging in the prescaling ansatz, we then find directly for the overoccupied system that

Breaking of scaling
The breaking of scaling by the Debye mass in the elastic collision kernel discussed in the main text is visualized in the left and middle panel of Fig. A1, where we show |p| 2 f (t, |p|) rescaled with fixed point exponents.The rescaling for simulations with only inelastic scatterings (middle) shows a very clear collapse of all curves from shortly after initialization over an evolution of six orders of magnitude, but a spread cannot be removed for all times even by time-dependent rescalings if one includes elastic scatterings (top left) due to the presence of the Debye mass.In the right panel we extract prescaling exponents for only inelastic scatterings C 1↔2 from different moments of the distribution function (see Eq. (A3)).The resulting approach of β(t) to the fixed point value β ∞ = −1/7 (dashed gray line) is demonstrated to be well-described over more than seven orders of magnitude in time by Eq. ( 5) of the main text (dashed black line) with µ α = 3 and µ β = −1, see Eq. (A15), even at surprisingly early times shortly after initialization.In the inset, we show that this solution initialized at a very early time t 0 Q ∼ 2 with t * obtained from β(t 0 ) ≡ β 0 captures the evolution at intermediate times t Q ∼ 10 4 only qualitatively.If we obtain t * in the same way at a later time t 0 Q ∼ 10 4 , Eq. ( 5) describes the late-time evolution quantitatively well as given by the solid black line.

Derivation of Eq. (8)
Here we give a brief account of how one derives Eq. ( 7) and Eq. ( 8) of the main text.As we demonstrated above, the FP collision kernel does not scale homogeneously under the prescaling ansatz and the corresponding prescaling solutions are therefore not captured by Eq. ( 5) of the main text.The t-and p-dependent contributions of the FP Large order behavior of transseries We iteratively solved for the coefficients in Eq. (9a) of the main text in Mathematica.This allowed us to generate exactly within several hours the lowest 38 orders of the expansion in inverse powers of log (Q t).In Fig. A2 we show that the generated coefficients β m,0 ∼ m! for m already as small as 10.Furthermore, we also found that the leading late time coefficient at each order obey β m,m−1 = (−1) m 1 7 , i.e. they form a geometric series allowing to simply resum their contribution.

1 n 0 FIG. 3 .
FIG.3.Comparison of prescaling trajectories in QCD kinetic theory simulations.Solid lines correspond to prescaling exponents extracted from simulations with different initial conditions, where for visibility only results from the first moment are displayed.
FIG. A2.Ratio test indication that the series (9a) has a vanishing radius of convergence