Stability analysis of non-thermal fixed points in longitudinally expanding kinetic theory

We use the Hamiltonian formulation of kinetic theory to perform a stability analysis of non-thermal fixed points in a non-Abelian plasma. We construct a perturbative expansion of the Fokker-Planck collision kernel in an adiabatic approximation and show that the (next-to-)leading order solutions reproduce the known non-thermal fixed point scaling exponents. Working at next-to-leading order, we derive the stability equations for scaling exponents and find the relaxation rate to the non-thermal fixed point. This approach provides the basis for an understanding of the prescaling phenomena observed in QCD kinetic theory and non-relativistic Bose gas systems.


I. INTRODUCTION
Dynamics of isolated quantum many-body systems quenched far from equilibrium has been an object of intensive study during recent years. Examples range from the dynamics of quark-gluon matter created in heavy-ion collisions [1,2] to quenches in ultracold atomic systems [3,4]. Starting from a far-from-equilibrium initial condition, these systems may exhibit a transient regime of self-similar evolution associated to a non-thermal fixed point [5,6]. As a result of such self-similarity, the nonequilibrium dynamics is fully encoded in a set of universal scaling exponents and functions. Over the last decade, the existence of non-thermal fixed points has been confirmed both experimentally [7][8][9] and in numerical studies [10][11][12][13]. On the theoretical side, progress has been made in predicting and explaining the observed scaling exponents using various techniques such as 1/N-resummed kinetic theory [14,15], low-energy effective description [16], and the functional renormalization group [17].
However, much less is known about how a general system evolves to such a self-similar regime. In [12] and [18], it was proposed that already before achieving fully developed scaling the system may exhibit a dramatic reduction in complexity such that its dynamics can be described by a few slowly evolving quantities. In particular, numerically solving the leadingorder QCD kinetic theory [19] it was observed that much before the scaling with universal exponents is established, the evolution is already governed by the fixed-point scaling function with time-dependent scaling exponents [18]. In this work, we are going to consider a toy model of an expanding Yang-Mills plasma and derive approximate equations that govern the dynamics of its scaling exponents. In particular, we derive the stability equations for scaling exponents, which can be interpreted as relaxation equations to non-thermal fixed point and demonstrate for the first time that the non-thermal fixed point is stable under small perturbations.

A. Prescaling
We begin our discussion with a quick overview of the concept of (pre)scaling, in particular in the context of heavy-ion collisions. At sufficiently high energies, where the gauge coupling is small due to asymptotic freedom [20,21], the time evolution of gluons (g) and quarks (q) is described by distribution functions f g,q (τ, p T , p z ). Since the system is longitudinally expanding, the distributions depend on transverse (p T ) and longitudinal momenta (p z ), and on proper time (τ) [22,23]. In the scaling regime, the original gluon distribution obeys with dimensionless τ → τ/τ ref and p T,z → p T,z /Q s in terms of some (arbitrary) time τ ref and characteristic momentum scale Q s . The exponents α, β, and γ are universal, and the nonthermal fixed-point distribution f S is universal up to normalizations [24], which has been established numerically using classical-statistical lattice simulations [1]. The exponents are expected to be α BMSS = −2/3, β BMSS = 0, and γ BMSS = 1/3 according to the first stage of the "bottom-up" thermalization scenario [2] based on number-conserving and small-angle scatterings, or α BD = −3/4, β BD = 0, γ BD = 1/4 in a variant of bottom-up including the effects of plasma instabilities [25]. Similarly, during prescaling the gluon distribution satisfies with non-universal time-dependent exponents α(τ), β(τ), and γ(τ). One can therefore regard prescaling as a partial fixed point at which the scaling function f S has already reached its fixed-point form, whereas the scaling exponents α, β, and γ still deviate from their asymptotic values.

B. Hamiltonian formulation of kinetic theory
In order to derive equations governing the prescaling dynamics, we are going to employ the Hamiltonian formulation arXiv:2203.02299v2 [hep-ph] 29 Jun 2022 of kinetic theory [26][27][28], the key points of which we will briefly summarize in this section. We start off with the general Boltzmann equation of a boost-invariant (in z-direction) and transversally homogeneous system [29]: Here, f is a distribution density, C[ f ] is a collision integral, τ is the longitudinal proper time, and p z and p T are longitudinal and traversal momenta, respectively. For the following discussion of prescaling, it will prove convenient to recast our problem into an (infinite) set of ordinary differential equations for moments of the occupation number f , Although for a general collision integral the expression does not have a simple form in terms of the moments n n,m , in this work, we are going to consider the kernel that is linear in f , and thus allows to reformulate the problem in the form ∂ log τ n n,m = −H n,m;n ,m n n ,m , with H n,m;n m = (2n + 1)δ n,n δ m,m − τq 2n(2n − 1)δ n−1,n δ m,m + m 2 δ n,n δ m−2,m .
Here, summation over repeated indices is implied and the momentum diffusion parameterq is parametrically given by [30,31] for SU(N c ) gauge theories in the limit of high occupancies. The Fokker-Planck-type collision integral (6) often serves as a toy model in the context of the bottom-up thermalization scenario. In the highly anisotropic limit one may, furthermore, neglect the transversal part, i.e., take C[ f ] = −q∂ 2 p z , so that the "Hamiltonian" reduces to H n,n;n ,m = (2n + 1)δ n,n − q2n(2n − 1)δ n−1,n δ m,m (10) and acquires a block diagonal structure. Here, for further convenience we have also introduced q ≡ τq.
Adopting Dirac notation we may write with n n,m ≡ n, m|ψ , H n,m;n m ≡ n, m Ĥ n , m , where the inner product is given by l 2 (Z ≥0 × Z ≥0 ) and {|n, m ≡ |n ⊗ |m } span the respective natural basis, Equation (11) with the Hamiltonian (10) will be the subject of discussion in the remaining text.

C. Adiabatic approximation
One notes thatĤ depends on y only through the parameter q(y), which immediately suggests applying the well-known adiabatic approximation from quantum mechanics. In contrast to quantum mechanics (of closed systems), however, the operatorĤ is not necessarily (anti-) Hermitian and hence the method requires some modifications. A particularly convenient generalization to the case of non-Hermitian yet diagonalizable Hamiltonians, which we summarize in App. A, was developed in [32]. The key idea is to, instead, consider with U being a transformation that diagonalizesĤ at a given instance y, The equation for |χ is given by Splitting the last term into its diagonal and off-diagonal parts, one immediately notices that, as opposed to the diagonal piecê H 0 , the off-diagonal termV is non-zero if and only if ∂ y q 0. This suggests that one may treatV as a perturbation as long as q depends on y slowly enough and thereby construct solutions to (16) in a perturbative manner: Here (see App. A), with |n and n = λ n +∂ y γ n being eigenvectors and eigenvalues ofĤ 0 , respectively, and γ n being the non-Hermitian generalization of the Berry phase, The coefficients C (l) n may be computed iteratively: where and V nm (y) = n|V(q(y)) |m .

III. PRESCALING
Now we are ready to study the time-dependent scaling exponents (prescaling) in the adiabatic approximation of the Hamiltonian formalism. To make analytical progress, we will define a small expansion parameter and will study the scaling exponents' behavior at leading and next-to-leading orders.
First, consider the l-th order contribution to the (n, m)-th moment of the distribution function: Here, we have used (25) to get the second line. For the perturbation (27) the iterative relation (22) takes a very simple form: Upon repeating the procedure (30) l times one readily obtains with C (l>k) km ≡ 0 and V(z) ≡ exp(2z) ∂ z q(z). We now recall that where T is a time-ordering operator. Since V(z i ) are ordinary numbers, the time-ordering operator drops out and we are left with where we went back to τ = τ 0 exp(y) (setting also τ 0 = 1 for brevity). Here, we have also introduced the functions and that will serve as an expansion parameter. Assembling everything together we end up with To go further, we need to truncate the series (36). To do so, let us first estimate the large-time behavior of the quantities entering the expansion (36). Near the fixed point, q ∼ τ 2α * −2β * −γ * implying the large-time behavior Hence, if 2α * −2β * −γ * +3 > 0, then we expect ∆(τ) to decay at large times τ and therefore may use it as a small parameter, at least when the scaling exponents are not too far off from their asymptotic values. Note that this condition holds both for the bottom-up [2] and for the modified [25] scaling solutions. On the contrary, for v the same analysis results in v(τ 1) ∼ const.
One may even estimate the asymptotic value as We thus conclude that the k-th term in n, m ψ (l) (τ) at large times scales as n, m ψ (l) (τ 1) with ∆(τ) playing a role of the small parameter. The leading order (LO) contribution to n, m ψ (l) (y) is then given by the l-th term in (36), the next-to-leading order (NLO) contribution is given by the (l + 1)-st term, etc. Importantly, this behavior is independent of l, which alludes to a possible need of resummation of all the terms of the same kind: n, m|ψ(τ) = n l=0 n, m ψ (l) (τ) see Fig. 1 for visualization and more details.
We are now in a position to derive equations that govern the prescaling dynamics at next-to-leading order. According to the above discussion, at this order the (n, m)-th moment of the distribution takes the form n, m|ψ(τ) NLO = (2n − 1)!!τ n−1q (τ) n C (0) Recognizing the binomial expansion we readily obtain with a nm = C (0) 0m (2n − 1)!! and b m = C (0) 1m /C (0) 0m .

B. Fixed-point equations
Up until this point, we have not assumed any particular ansatz for the time evolution of moments n n,m of the distribution function. If the prescaling assumption (2) holds, however, then we can recast the equations for the moments in terms of time-dependent scaling exponents. Following the original work [18], in order to reflect instantaneous scaling properties we redefine exponents in (2) as which for constant α reduces to the power law τ α . The rate of change of a particular moment n n,m as well as of the momentum diffusion parameterq is given by a linear combination of scaling exponents: in d = 3 spatial dimensions. This also implies Taking then the log of both sides of (45) and then the derivative with respect to log τ we end up with Since the NLO approximation is O(∆), we have to expand the log term on the right-hand side to first order in ∆ to be consistent. After some simple algebra, one then eventually arrives at First, we observe that in order for this equation to hold for any n and m (as it should during prescaling) one has to impose One immediately recognizes in (51) the scaling relation that follows from conservation of the total particle number [24]. This reflects the particle-number-conserving nature of the elastic collision kernel. It is then suggestive to also demand that the term containing m and the term containing n should individually vanish identically, too. This would result in another constraint which together with (51) indicates energy conservation [24]. The remaining equation then reads For this condition to hold b m has to be m-independent. Since see (28), the latter holds as long as n 1m (τ 0 )/n 0m (τ 0 ) does not depend on m. An important class of distributions for which this condition is always satisfied is given by separable distributions, i.e., f (τ 0 , p T , p z ) = f 1 (τ 0 , p T ) f 2 (τ 0 , p z ). We have already derived the equation (48a) governing the dynamics of v during prescaling. To obtain a similar equation for the remaining scaling exponent γ, we first take one more logarithmic derivative of both sides of (53): whereȮ ≡ ∂ log τ O. Finally, using (48a) and (48b) and imposing the constraints (51) and (52) one ends up with the system of differential equations:κ where we have introduced κ = (γ, v) and The above equations resemble flow equations describing a running of couplings in the context of renormalization group flow. The flow diagram of (57) is depicted in Fig. 2. Scaling is achieved when the flow reaches a fixed point Using (57) one recognizes the standard bottom-up scaling exponents, together with cf. (39), as a stable fixed point of the flow equations (56). Indeed, using the standard notation δκ ≡ κ − κ * one has δκ bottom-up The corresponding characteristic polynomial reads (λ + 2) (λ + 4/3) resulting in two (simple) eigenvalues with the respective eigenvectors The general solution near the fixed point is therefore given by or explicitly,

IV. CONCLUSIONS
In this work we studied the self-similar evolution phenomena in Fokker-Planck type kinetic theory. Using the Hamiltonian formalism of kinetic theory and adiabatic approximation, we were able to derive the flow equations of the timedependent scaling exponents. The fixed point of scaling exponents for the Fokker-Planck kinetic theory coincides with the scaling exponents characterizing the early stage of the bottomup thermalization scenario [2].
Working at next-to-leading order in the small expansion parameter, we found the relaxation rate for scaling exponents to the fixed point and demonstrated its stability. This analysis lays ground for the study of scaling phenomena in more complex systems, such as full QCD kinetic theory. * * * We note that an analysis of time-dependent scaling exponents in Fokker-Planck kinetic theory was performed independently by Jasmine Brewer, Bruno Scheihing-Hitschfeld and Yi Yin and made public simultaneously to the present manuscript [33].

APPENDIX Appendix A: Non-Hermitian adiabatic expansion
In this appendix, we show how one can systematically solve equations of the kind ∂ y |ψ(y) = −Ĥ(q(y)) |ψ(y) , whereĤ is diagonalizable but not necessarily (anti-)Hermitian. Following [32], we first want to find a transformation U that diagonalizesĤ at each instance y, For example, in the standard basis {|k } this transformation reads Let |ψ be a solution to the equation (A1). Define the "equivalent solution", that satisfies withĤ e (q) =Ĥ d (q) + U(q) −1 ∂ y U(q).
We now split the second term into diagonal and off-diagonal parts and introducê so thatĤ e (q) =Ĥ 0 (q) +V(q).
One can already guess that the diagonal partĤ 0 (q) governs adiabatic element of the evolution, whereas the off-diagonal pieceV(q) gives rise to non-adiabatic transitions between the quasi-energy levels. Furthermore, sinceV(q) vanishes when there is no time-dependence we anticipate that one can treat V(q) as a perturbation when q(y) depends on y slowly enough. We will therefore look for solutions in the form Here, ∂ y χ (l) (y) = −Ĥ 0 (q(y)) χ (l) (y) −V(q) χ (l−1) (y) , (A11b) for l ≥ 1. The zeroth order solution is given by with |χ(0) ≡ U(q(0)) −1 |ψ(0) . It is convenient to work in the basis of eigenvectorsĤ 0 , which we can choose to be The corresponding eigenvalues read n = λ n + ∂ y γ n , with γ n (y) = y 0 dz n| U(q(z)) −1 ∂ z U(q(z)) |n .
Expanding the l-th order solution in this basis as where we have usedĤ 0 |m = m |m . First, we notice that the last term on the left-hand side cancels the first term on the right-hand side. Multiplying both sides by n| and using the orthogonality condition n|m = δ mn we then readily obtain ∂ y C (l) n (y) = − m V nm (y) exp − y 0 ds ω nm (s) C (l−1) m (y), (A17) with ω nm (y) ≡ m (q(y)) − n (q(y)) (A18) and V nm (y) ≡ n|V(q(y)) |m .
(A20) As a final remark, we note that tedious, yet straightforward computations show that there is also no ambiguity regarding the choice of the instantaneous eigenfunctions |v k (q) . In other words, ψ (l) are invariant under reparameterizations |v k (y) → e φ k (y) |v k (y) at each order of perturbation theory.