Transasymptotics and hydrodynamization of the Fokker-Planck equation for gluons

We investigate the non-linear transport processes and hydrodynamization of a system of gluons undergoing longitudinal boost-invariant expansion. The dynamics is described within the framework of the Boltzmann equation in the small-angle approximation. The kinetic equations for a suitable set of moments of the one-particle distribution function are derived. By investigating the stability and asymptotic resurgent properties of this dynamical system, we demonstrate, that its solutions exhibit a rather different behavior for large (UV) and small (IR) effective Knudsen numbers. Close to the forward attractor in the IR regime the constitutive relations of each moment can be written as a multiparameter transseries. This resummation scheme allows us to extend the definition of a transport coefficient to the non-equilibrium regime naturally. Each transport coefficient is renormalized by the non-perturbative contributions of the non-hydrodynamic modes. The Knudsen number dependence of the transport coefficient is governed by the corresponding renormalization group flow equation. An interesting feature of the Yang-Mills plasma in this regime is that it exhibits transient non-Newtonian behavior while hydrodynamizing. In the UV regime the solution for the moments can be written as a power-law asymptotic series with a finite radius of convergence. We show that radius of convergence of the UV perturbative expansion grows linearly as a function of the shear viscosity to entropy density ratio. Finally, we compare the universal properties in the pullback and forward attracting regions to other kinetic models including the relaxation time approximation and the effective kinetic Arnold-Moore-Yaffe (AMY) theory.


I. INTRODUCTION
Relativistic fluid dynamics is an effective theory which describes long-wavelength phenomena. It is widely accepted that its regime of validity is restricted to systems near local thermal equilibrium. However, this traditional paradigm has recently been challenged by the overwhelming success of hydrodynamic models in describing experimental data in high energy nuclear collisions [1][2][3][4][5][6][7] as well as cold atom systems [8][9][10][11]. In these systems the initial state is far from local thermal equilibrium, and it is not fully understood how hydrodynamic behavior emerges. The search for a kinetic framework that describes far-from-equilibrium plasmas has been one of the most important research subjects in high energy nuclear collisions and condensed matter physics [12][13][14].
An important development in non-equilibrium dynamics was the discovery of emergent hydrodynamic behavior in far-from-equilibrium conditions which can be understood in terms of the mathematical theory of resurgence [15]. In this work the authors consider an extended hydrodynamic model, the Israel-Stewart equation [16], and apply it to a strongly coupled plasma undergoing Bjorken expansion. Subsequently, similar findings were obtained in many other transport models. These results show a deep connection between nonlinear relaxation towards hydrodynamic behavior, also known as hydrodynamization, and transasymptotics and transseries [17][18][19][20][21][22][23][24][25][26].
Since then a very interesting and rich physical picture has emerged: The nonlinear relaxation process towards hydrodynamic behavior, also known as hydrodynamization, is driven by the decay of non-hydrodynamic degrees of freedom. Once the non-hydrodynamic modes have died out the system enters into the hydrodynamic attractor which is entirely determined by the standard asymptotic gradient expansion. This new insight might be able to explain why hydrodynamic models work very well when applied in extreme experimental scenarios such as ultrarelativistic heavy ion collisions [27,28].
Another interesting development in the understanding of far-from-equilibrium attractors in relativistic nonequilibrium dynamics is a phase space analysis using the language of non-autonomous dynamical systems 1 [29][30][31][32]. This particular point of view allows us to characterize in simple terms the behavior of the solutions, described as flows, either at early or at late times. For instance, a global and local phase space flow analysis led to the conclusion that a large class of kinetic models under-going Bjorken expansion [33] hydrodynamize in the long time limit [21,31,32,34], whereas systems undergoing Gubser flow [35,36] never do [29,32]. Moreover, the flow structure in phase space together with the symmetries of the dynamical system constrains the asymptotic behavior of the solutions of the ODEs. For example, for weakly coupled boost invariant systems it was demonstrated [31] that the solutions of the moments equations admit power law series expansions at early times, while at late times linear perturbations of the moments decay exponentially [21,30,31,37]. This is due to the nature of the early and late-time attractors.
Attractors are understood as regions in phase space where flows accumulate either in the long or short time limit. However, in non-autonomous dynamical systems the past and future of the evolution are not the same since time translation invariance is explicitly broken. Each flow solution φ ≡ φ(φ 0 ; t, t 0 ) is written in terms of its initial value φ 0 and its initial and final times, t 0 and t respectively. In this context, it is important to differentiate the backward and forward asymptotic regions of the dynamical system. We refer to a forward attractor as an asymptotic limit of the flows at which solutions converge when t → ∞ while the initial time t 0 is fixed. In contrast, the pullback attractor is defined in the limit t 0 → 0 while keeping t fixed. It is important to emphasize that for non-autonomous dynamical systems both limits do not commute [38,39].
Most previous work on the convergence properties of kinetic equations in relativistic transport theory was based on very simple collision terms, e.g. the relaxation time approximation (RTA). In this work we study the resurgent asymptotic properties of the Boltzmann equation for a boost invariant Yang-Mills plasma governed by the weak coupling collision term in the small-angle approximation. In this limit the Boltzmann equation can be written as a nonlinear Fokker-Planck equation.
The Fokker-Planck equation (FPE) for gluons is of great interest since it captures essential aspects of the early time dynamics of QCD matter produced in ultrarelativistic heavy ion collisions [40,41]. Due to the highly nonlinear structure of the collision kernel in the small angle approximation, the FPE has been solved mostly by numerical means [42][43][44][45][46][47][48]. Very few analytical solutions for rapidly expanding systems are known in the literature [49]. One of the main achieved goals in this work is to fill this gap.
Following Grad's method [50] we map the mathematical problem of solving the FPE onto finding solutions to the kinetic equations for the Legendre moments c l [29][30][31][32]. The moment method turns out to be not only quickly convergent from a numerical point of view (see App. G of Ref. [31]), but it also provides a unique approach to understand non-trivial aspects of hydrodynamization. New solutions to the equations of motion of Legendre moments are derived by employing methods developed in the context of stability analysis of non-autonomous dynamical systems [38,39], as well as techniques from superasymp-totic and hyperasymptotic analysis [51][52][53][54]. These tools allow us to analyze the hydrodynamization process in two distinct regimes characterized by the size of dissipative corrections: Kn 1 (UV/early time) and Kn 1 (IR/late time).
The solutions of the kinetic equations for the Legendre moments enable us to extend the definition of transport coefficients to the far-from-equilibrium regime. Within our approach transport coefficients depend on the deformation history of the fluid, i.e. its rheology, and their values change as a function of the Knudsen number. Such dependence of the transport coefficients on the typical gradient size is a salient property of for non-Newtonian fluids. Our findings provide further evidence for the connection between hydrodynamization and the transient rheological behavior of the plasma [30,31].
An interesting aspect of our study is the fact that the evolution of the transport coefficients is determined by a RG flow equation where the role of the RG scale is played by the Knudsen number. Our work draws inspiration from recent arguments that any RG flow can be viewed as a dynamical system [55,56]. We show that, conversely, certain dynamical systems can be understood as RG flows, provided the existence of the slow invariant manifold. This idea was first explored in the case of boost invariant plasmas governed by a Boltzmann equation in the RTA approximation [30][31][32].
Finally, our study naturally explains the origin of the UV power series expansion considered first in [31] and later studied in [26,[57][58][59]. The power series emerges by analyzing the stability behavior of solutions close to the UV fixed point. It is shown that the finite radius of convergence of the UV power series expansion depends on the value of the shear viscosity over entropy ratio η/s. We also outline the intriguing universal properties of the FPE related to the pullback and forward attractors, and compare them with other kinetic models such as the RTA and the full leading order AMY kinetic theory [60].
The outline of this work is as follows. In the next section we review the basic ideas behind the FP equation and introduce the expansion of the distribution function in terms of its moments whose evolution equations are derived. For pedagogical purposes we study extensively a truncation of this dynamical system in Sect. III where we present in detail the transasymptotic techniques studied in this article. In Sect. IV we demonstrate that the general solutions of the Fokker-Planck equation are written in terms of multiparameter transseries after resumming the fluctuations around the UV and IR regimes, respectively. We analyze the universal properties of nonequilibrium Yang-Mills attractors of different theories in Sect. V. We conclude by giving final remarks. The technical details of our work are presented in the appendices.

II. YANG-MILLS TRANSPORT EQUATION IN THE SMALL ANGLE APPROXIMATION
We consider an interacting gluon plasma described by the Boltzmann equation in the small angle approximation [61]. We focus our discussion on the case of an expanding system with longitudinal boost invariance [33]. The dynamics is invariant under the ISO(2)⊗SO(1, 1)⊗ Z 2 symmetry group. The symmetry becomes manifest in Milne coordinates x µ = (τ, x, y, ς), where τ = √ t 2 − z 2 is the longitudinal proper time, and ς = tanh −1 (z/t) is the spatial rapidity. In this coordinate system the metric is simply g µν = (−1, 1, 1, τ 2 ). In Milne coordinates the Fokker-Planck equation for the one particle distribution f (τ, p T , p ς /τ ) ≡ f p is [40,46,47] (1) The collision kernel C[f p (τ )] in the small angle approximation takes the form [40,46] where we introduce the t'Hooft coupling λ Y M = g 2 N c ≡ 4πα s N c . The integrals J and K are given by The Coulomb logarithm l Cb in the RHS of Eq. (2) is a divergent integral of the form The IR momentum divergence originates from 2 → 2 scattering with small momentum transfer. In QCD these divergences are regularized by static and dynamic screening. The corresponding mass scale is on the order of the Debye mass. Near equilibrium, and for particles obeying Bose-Einstein statistics, we have [62] where ζ(n) and Γ(n) are the Riemann and Gamma functions respectively. The UV momentum cutoff is taken to be the mean p 2 T , which close to equilibrium is As a result, near equilibrium the Coulomb logarithm (4) is approximately given by with A = 72ζ(5)/ζ(3) ≈ 62.1. This estimate gives a Coulomb logarithm which is a constant, independent of the energy density of the medium, but logarithmically dependent on the coupling constant. This approximation has been used in many previous studies [45,46,59,63]. Replacing f eq. p with the general non-equilibrium f p in (5) and (6) and evaluating the integrals numerically at each time step allowed the authors of Refs. [42,48] to account for the time-dependence of the UV and IR momentum cutoffs.
The approximation of a constant Coulomb logarithm gives the correct dependence of the shear viscosity on the coupling constant in the near-equilibrium, weak coupling limit. However, the numerical pre-factor does not agree with calculations based on the full Hard Thermal Loop (HTL) result [64,65]. This issue was addressed in [66], where the authors propose a simple regulator that reproduces the leading order HTL result for drag and momentum diffusion in soft 2 → 2 scattering. This is a useful prescription, but it does not affect the late-time emergent hydrodynamic behavior. We shall comeback to this issue in Sect. V.

A. Ansatz for the distribution function
One of the most widely used methods to solve the Boltzmann equation was developed decades ago by Grad [50]. In this approach the problem of solving the Boltzmann equation is converted into a nonlinear set of PDEs for moments of the one-particle distribution function. This approach is very useful when analyzing the resurgence properties of the nonlinear ODEs for the moments [29][30][31][32].
In this work we will consider the following ansatz for the distribution function in a boost invariant system [30][31][32] f where p τ = p 2 x + p 2 y + (p ς /τ ) 2 , cos θ p = p ζ /(τ p τ ), and P 2l are the Legendre polynomials. In Eq. (8) the equilibrium distribution function f eq p is where T is the temperature of the system which is defined below via the Landau matching condition and ν ef f are the effective degrees of freedom which for simplicity we set to be ν ef f ≡ 1. We implicitly assume that f p is independent of spin and color. The ansatz (8) is consistent with the restrictions imposed by ISO(2) ⊗ SO(1, 1) ⊗ Z 2 .
In particular, the distribution function is invariant under the action of the Killing vectors φ i of this symmetry group, i.e., ∂f p /∂φ i = 0 [30][31][32]. We fixed µ ≡ 0 in the equilibrium distribution function (9) so that the system is not overpopulated and thus, the possible formation of a Bose-Einstein condensate [46] is not included in our approach. Furthermore, the ansatz (8) does not carry information about the non-linear relaxation of the transient  high energy tails of the distribution function which have  been studied previously within the moment method [30 -32, 46, 47, 67-70]. A more general ansatz which encodes some of the information on the high energy tails was discussed recently in Ref. [45]. The Legendre moments c are directly determined by Eq. (8), where we denote · · · X ≡ p · · · f X p and the momentum measure is p ≡ d 2 p T dp ς / (2π) 3 τ p τ . If the system reaches the thermal equilibrium state, the moments c eq. l = δ l0 . For the Bjorken flow the energy momentum tensor T µν = p µ p ν is given by [30-32, 71, 72] where we denote the time-like vector identified with the fluid velocity u µ = (1, 0, 0, 0) (with u µ u µ = −1), the space-like normal vector pointing along the ς direction is l µ = (0, 0, 0, 1) (with l µ l µ = 1) and the projection operator Ξ µν = g µν + u µ u ν + l µ l ν which is orthogonal to both u µ and l µ respectively. The energy density , longitudinal and transverse pressures, P L and P T respectively, are written in terms of the angular moments as follows It follows that = 2P T + P L as expected. The Landau matching condition for the energy density = eq. ≡ (−u · p) 2 eq. implies c 0 ≡ 1. For the Bjorken case the normalized pressure anisotropy is the ratio of the independent shear viscous tensor over the energy is written in terms of the Legendre moment c 1 , i.e., Moreover, the non-negativity property of the longitudinal and transverse pressures P L(T ) ≥ 0 implies −5/2 ≤ c 1 ≤ 5. This bounds the basin of attraction from below and above and is satisfied in general only by the exact kinetic solution to the FP equation (2). Nonetheless, it is known that any perturbative approach which aims to find an approximate solution to the Boltzmann equation does not necessarily obey this constraint [73].

B. Evolution equation of the Legendre moments evolution
We truncate the expansion in Legendre polynomials at order L and write the Legendre moments as a vector c(τ ) ≡ (c 1 (τ ), · · · , c L (τ )) . Following the procedure outlined in [30][31][32] we find that the evolution equations for the temperature and the Legendre moments are given by the following nonlinear coupled of ODEs with and The matrix elements are defined by The Lyapunov exponents and the vector Γ are defined by Finally, we have defined the matrices In the previous expressions the parameters κ, θ 0 and κ are given by respectively The derivation of the equations of motion (14) is presented in App. A. The set of ODEs (14) constitutes a nonlinear non-autonomous dynamical system [38,39] due to the explicit dependence on the proper time τ in the RHS of these equations. The strength of the effect of the collisions enters in the ODEs via the parameter θ 0 (20); we thus will vary this parameter instead of the t'Hooft coupling λ Y M when showing the numerical results. The nonlinear nature of the FP equation (1) is manifest in the mode-mode couplings among different moments, and in the temperature appearing in the RHS of Eq. (14). Finally, the time-evolution of the energy-momentum tensor (11) can be fully reconstructed from the solutions of the temperature T and the full set of Legendre moments c l . The equations of motion of the Legendre moments (10) for the conformal Boltzmann equation within the relaxation time approximation (RTA) are given by (see [30,31]) Here Θ 0 is a proportionality constant between the relaxation time and the shear viscosity over entropy ratio. This constant can be determined from relativistic kinetic theory methods which for the conformal RTA approximation gives us τ r = 5 T η s [31,[74][75][76][77][78][79][80][81][82][83][84][85], i.e. Θ 0 ≡ 5η/s. At linear order the form of the evolution equations for the Legendre moments for the RTA and FP, Eqs. (21) and (14) respectively, is very similar. At this order the main difference between the two models lies in the specific relaxation scale of the associated mode c l . For the RTA approximation all modes decay with the same relaxation scale τ r while in the Fokker-Planck case each c l has a characteristic relaxation scale λ(l) −1 which increases as a function of the l (see Eq. (20)). Therefore for the FPE there is a clear hierarchy of fast and slow modes 2 .
The nonlinearities of the FPE are encoded in the mode coupling between moments of different order and in the implicit dependence of the temperature on c 1 . These effects drive the system away from equilibrium and delay the relaxation to the equilibrium state. For instance, for the RTA Boltzmann equation a new non-hydrodynamic mode survives in the long wavelength limit due to a nonlinear coupling to the shear viscous tensor [31]. This result contradicts the common assumption of taking the typical relaxation time scale of a mode as a guide for constructing effective theories in the long time limit.

C. Dimensional reduction
There is an important simplification of the nonlinear ODEs (14) which makes the asymptotic analysis simpler. It is possible to dimensionally reduce this dynamical system from L + 1 → L by introducing the variable 3 w = τ T (τ ) ∈ R + . Since w ∼ Kn −1 [29][30][31][32], the variable w encodes the strength of the dissipative corrections. This fact will be important for re-interpreting the dynamical system of ODEs as RG flows in Sect. III B 1.
2 For hard spheres and within the kinetic theory framework similar results have been found in the past. In those models a hierarchy of scales between the fast and soft modes also emerges [68][69][70]. 3 The resurgence analysis can also be studied in terms of the original variable τ [31]. See the general Fokker-Planck case in App. C and the RTA Boltzmann case is extensively discussed in Appendixes A-D in Ref. [31].
In terms of w Eqs. (14) can be written as The dynamics of this dimensionally reduced system of ODEs depends only on the Legendre moments c l 4 and does not involve the temperature variable explicitly. It is important to emphasize that the dimensionally reduced dynamical system does not preserve the topological properties of the original nonlinear ODEs (14). For instance, a new coordinate singularity emerges when c 1 = 20 (see the denominator in the integrand of the damping function D(w, w 0 ) in Eq. (22a)) which does not exist in the τ variable. Furthermore, in the w variable, one can find a set of bounded solutions at w 0 → 0 for a fixed w, the socalled pullback attractor [31,32]; whereas the singularity at τ = 0 forbids the existence of any bounded solutions in the UV. This is of course expected on many grounds but a simple explanation is that the system undergoes a topological change under τ → w that lifts the τ = 0 singularity at the expense of introducing a new singularity at c 1 = 20 [31,32].

III. TRANSASYMPTOTIC ANALYSIS: THE TRUNCATED L=1 CASE
Before presenting the transseries solutions for the nonlinear ODEs (22) we first illustrate our techniques by considering the case where the system is truncated to a single non-hydrodynamic degree of freedom. Physically, this degree of freedom corresponds to the viscous shear tensor. This warm-up exercise illustrates the main aspects of the general asymptotic analysis to be discussed in Sect. IV.
Let us assume that the moments c l ≡ 0 for l > 1. In terms of the variable w, the truncated l = 1 case read as 4 The index i of vectors and matrices runs from i = 1, 2, ..., L.
If one performs a truncation of the non-autonomous dynamical system of Eq. (22) it is understood that c l ≡ 0 for l > L.
Before discussing the transasymptotic analysis we would like to comment on the existence of the pullback and forward attractors of the truncated ODE (23a). The extension to the full set of nonlinear ODEs (22) is similar but more difficult to visualize [31]. For now, we note that we can identify attractors by inspecting the flow diagram of the ODE (23a). We will characterize the attractors more fully below. In the top panel of Fig. 1 we observe that in the IR regime all flows converge asymptotically towards the value c 1 → 0 when w 1 regardless of their origin. This shows that there is only one IR fixed point associated to the forward attractor.
In the UV limit, on the other hand, there are two fixed points (black dots in Fig. 1) which are located along the c 1 axis in the small w 0 → 0 limit (w 0 is the initial value of w). The UV behavior of the flow lines in each fixed point is rather different. For the UV fixed point located in the positive c 1 region one observes that the flow lines emerge out of it in all the directions and thus, this fixed point is a source. In the negative c 1 region, we observe that near the UV fixed point the flow lines coming from above or below repel each other so the fixed point is a saddle. However, although the flow lines repel each other, they merge quickly in the vicinity of this fixed point. Thus, the UV fixed point located in the negative c 1 region is identified as the pullback attractor.
We will see that the behavior of the transseries solution near these fixed points is rather different. We note that that the behavior of the flows in the τ variable differs from the one in w. In the former case the temperature is a dynamical variable. It diverges in the limit τ → 0 due to the presence of the singularity at τ = 0; therefore there is no meaningful definition of a pullback attractor in the τ variable (see also [31,32]).

A. Transseries solution in the IR limit
The leading order term of each Legendre moment in Eqs. (22) is c l = O(w −l ) ∀ l > 0 in the IR regime [31]. Hence, Eq. (23a) admits the following asymptotic expansion for c 1 The coefficients u the first three coefficients are given by Herethe reality condition of c 1 was implicitly taken into account. The hydrodynamic gradient expansion of the shear viscous tensor [86] provides the following expression for c 1 [30,31] where we use Eq. (13). Comparing the previous expression and Eqs. (25) one concludes with λ(1) andκ given in Eqs. (20). We will see below that λ(1) is related to the rate at which fluctuations decay close to the IR fixed point, i.e. the Lyapunov exponent. The coefficients u 1,k that enter the IR expansion (24) are understood as transport coefficients in the IR limit. Furthermore, the set of the identities (27) turn out be of importance when generalizing the concept of a transport coefficient to the far-from-equilibrium regime [30,31] as we shall discuss in Sect. III B.
It is known that the IR expansion (24) is divergent since its coefficients grow factorially [15]. These expansions are merely formal expressions and emerge as asymptotic solutions of a certain class of differential equations [52]. In order to see this let us simply perform a linear perturbation around the perturbative expansion (24). Thus, by shifting c 1 →c 1 + δc 1 withc 1 = ∞ k=0 g k w −k while keeping terms up to O(w −1 ) one obtains the following linearized differential equation for the perturbation δc 1 whose solution reads where σ 1 is the integrating constant, S 1 is recognized as the Lyapunov exponent and b 1 is the anomalous dimension. The presence of the exponential terms is an indication that in order to capture the transient behavior of the solutions it is needed to go beyond the perturbative IR expansion (24). The leading exponential term (29) is the first type of a large set of exponentials terms that appear when summing over the fluctuations around the IR fixed point. In order to include systematically these type of terms we follow Costin's prescription [54]. First, notice that in the IR limit the ODE (23a) takes the following asymptotic form where R(c 1 , w) is a non-linear polynomial function of c 1 and w. The asymptotic limit of the differential equation (30) coincides with the prepared form of the generic class of differential equations studied by Costin [54]. Thus, given the regularity of the solutions at w → ∞ as well as the non-vanishing value for the Lyapunov exponent S 1 (29), Eq. (23a) has an exact transseries solu-tion [54] From the physical point of view, the transseries (31) describes deviations from thermal equilibrium due to exponentially damped modes multiplied by gradient terms. The exponential damped terms, usually called nonhydrodynamic modes, play a role analogous to instantons in quantum field theory and quantum mechanics. In order to determine the coefficients u (n) 1,k we simply equate (31) with Eq. (23a). As a result we get the following recursive relation where u (m) 1,k = 0 for k < 0 and we take u We conclude this section by presenting some numerical results. In Fig. 2 we show the Legendre mode c 1 as a function of w for the exact numerical solution (red line), NS (blue dotted line), NS + 2nd order (green dashed line) and truncated transseries solutions where we added two (orange dashed line) and five (magenta dashed line) transmonomials to the 1st IR perturbative order. The initial condition for the numerical solution of the ODE (23a) was c 1 (w 0 = 0.03) = {6.60119, 0, −3.02976} while fixing θ 0 = 1 5 . Note that at each order of the resummed theory we have to numerically match the transseries parameter σ. Very few cases are known where one can determine σ exactly [26]. In the present work we have determined σ using a numerical least-square fit [31]. This leads to some uncertainty, which is not unexpected given the difficulty of matching the IR and UV data with- out performing an all-order resummation. The shaded yellow area in Fig. 2 shows the variation of σ between each truncated transseries solution. Fig 2 shows that the perturbative series at first and second orders does not match the numerical results for small values of w. On the other hand, the truncated transseries (31) where only a few transmonomials (n = 2, 5 respectively) were included is very close to the exact numerical result. We note that some care has to be taken in making such comparisons. Any truncated transseries will potentially exhibit large deviations from the exact numerical result in some range of w. This is due to the growth of the inherent error associated with any type of truncation scheme, see App. B. Any resummed perturbative series with a few transmonomials has a finite radius of convergence. The full transseries solution coincides with the exact theory only when summing over all the non-perturbative sectors of the perturbative series [54].

B. Transasymptotic matching
An interesting property of the IR transseries solution (31) emerges when rearranging the terms close to the IR fixed point. This procedure is known as transasymptotic matching, and it is a well know feature of transseries solutions for a general class of differential equations [53,54]. In the small w limit there is a one to one competition between the exponentially decaying terms entering in the transseries and the inverse powers in w. However, close to the IR fixed point the instantonlike contributions are more suppressed than the leading IR perturbative terms. It is in this regime where one can indeed reshuffle and resum the small exponential terms such that with ζ 1 (w) = e −S1w w −b1 and the transasymptotic functions G k (σ 1 ζ(w)) given by The transasymptotic functions G k (σ 1 ζ(w)) effectively resum the full set of instanton-like contributions at a given order k in perturbation theory. Interestingly this matching procedure is not only valid close to the IR but it extends up to the UV, so in this sense it is transasymptotic [52]. Each G k (ζ(w)) is an analytic, Borel summable and convergent functions even if one truncates at a given order the IR perturbative expansion Eq. (33). The full transseries solution (33) is not Borel summable due to the singularity of the original differential equation (23a). This singularity can be easily determined by taking the inverse Laplace transformation of the ODE (23a) (cf. Ref. [32]). The transasymptotic matching procedure coincides with the exact solution of the ODE only when summing over all the perturbative and non-perturbative sectors, i.e. k, n → ∞ in the upper limit of the sums of Eq. (33) [52][53][54]. In general, a truncation of the perturbative expansion leads to a solution with a finite radius of convergence.
In the large w limit the transasymptotic functions G 1,k (σ 1 ζ(w)) is uniquely determine by the coefficients of the asymptotic IR expansion, namely If one uses Eqs. (25) and (27) for k = 1, 2 in the previous expression one finds the following identities which relate the transport coefficients determined from linear response theory with the asymptotic regime of the function G 1,k (σ 1 ζ(w)) in the IR regime [30,31] Thus, the transasymptotic matching procedure automatically allows us to generalize the concept of a transport coefficient beyond the linear regime by promoting the functions G 1,k to be an effective transport coefficient [30,31]. For first and second order transport coefficients we have, respectively As a result, a new physical picture for nonlinear transport emerges: Summing over all the non-perturbative sectors at each order of the IR expansion leads to an effective renormalization of the transport coefficients. Each transport coefficient exhibits a transient non-Newtonian behavior while relaxing towards its asymptotic value in the IR. This transient rheological behavior is described by a dynamical RG flow equation. This result extends earlier work in the context of the Boltzmann equation in the relaxation time approximation [30,31].
Numerical results for the renormalized transport coefficients in Eqs. (37) are shown in Fig. 3. In order to obtain these results we have to calculate the functions G 1,k with k = {1, 2} using Eq. (34). In that expression the transmonomials are given by Eqs. (31) and the coefficients u (n) 1,k are numerically determined by solving order by order the recursive relation (32). Fig. 3 shows the w-dependence of the ratios between the dynamically renormalized transport coefficient over its asymptotic value. The top panel shows the shear viscosity over entropy density ratio, and the bottom panel shows the second order transport coefficients.
The values of the transasymptotic parameters were determined by fitting the transseries solution (31)  Both panels illustrate that regardless of the initial condition each renormalized transport coefficient reaches its asymptotic value as it is expected from the properties of transasymptotic functions G 1,k in Eq. (36). On the other hand, prior to the relaxation to the asymptotic values both renormalized transport coefficients feature a transient phase which depends on the deformation history of the fluid and thus, its rheology. Clearly this transient phase is not universal since it depends on the initial condition. However, one observes that both ratios shown in Fig. 4 saturate their value to the unity around w ≈ 0.4. These results are consistent with the ones shown in Fig. 2.

Dynamical RG flow equation
When the transasymptotic matching condition is imposed over the entire domain of w, a nonlinear set of PDEs for the functions G 1,k is obtained after inserting Eq. (33) into the original ODE (23a), whereζ = ζ∂/∂ζ = ∂/∂ log ζ. Solving these PDEs is equivalent to summing over all fluctuation contributions around the IR. More importantly, these PDEs can be reinterpret as RG flow equations for the transport coefficients given the identities (36). This statement can be rigorously proven within the gradient descent approach to the RG flows of quantum field theories in the context of dynamical systems, c.f. [30,31,56]. Furthermore, a RG flow picture from a global point of view holds for our original dynamical system (22) provided the existence of of an effective potential for the ODEs, also known as Lyapunov function. Its mathematical existence is proven in App. D.
Following Refs. [30,31] the ODE for c 1 in terms of the w variable, Eq. (23a) can be rewritten as whereκ is given in Eq. (20). The previous equation is the RG equation for the moment c 1 and the function β 1 (c 1 , w) encodes the dependence on the RG time w.
The variable w ∼ Kn −1 plays the analogous role of the energy scale in QFT since it determines the size of deviations from the equilibrium state. We can think of the fugacity ζ 1 as independent of w [87], so that G 1,k (σ 1 ζ 1 ) can be regarded as a simple coefficient in the transseries solution of c 1 (31). Then, the following identity follows [30,31] where we explicitly used Eq. (39a) 6 . The Cauchy integral formula in the previous expression simply picks up the kth coefficient of the termζ 1 c 1 . In order to determine the RG equation which describes the transient behavior of the transport coefficients we consider an arbitrary observable O = O(G 1,k (σ 1 ζ 1 )). Its change along w results in the following RG equation The renormalization of the observable O is obtained by solving the previous equation. In this equation the term ζ 1 G 1,k (σ 1 ζ 1 ) is given by the identity (40) which encodes the dynamics of the original ODE via the beta function 6 The ordinary derivative respect to log w was rewritten as

C. Transseries solutions in the UV
In the previous section we showed that perturbations decay exponentially in the vicinity of the IR fixed point (29). For an arbitrarily chosen ODE linear perturbations do not necessarily decay exponentially in the neighborhood of a fixed point. For instance, consider a parameter x ∈ (0, ∞) and a function y(x) which satisfies a well defined ODE dy(x)/dx = L(x, y(x)). If this ODE has a source fixed point y 0 in the limit x → 0, i.e. the Lyapunov exponent λ close to y 0 is positive in this limit, then the solutions in a vicinity 7 around it behaves as y(x) ∼ e λx ∼ ∞ n=0 (λ x) n in the limit x 1. This approximate solution is a power law series with a finite radius of convergence which can be extended by analytic continuation.
The stability analysis of the RTA Boltzmann equation indeed showed that in the limit when the Knudsen number is large, that is in the w → 0 limit, power law series solutions for the Legendre moments emerge [31]. Similar findings were reported in [57][58][59]. The power series expansion is rather different from non-hydrodynamic modes in the gradient expansion (24). As we shall see in this section the main differences can be inferred from the behavior of linearized perturbations around the UV and IR fixed points, respectively. The case studied here illustrates the importance of the flow structure in phase space in controlling the functional form of solutions for generic dynamical systems of ODEs [32]. In the following we will explain more carefully the emergence of the power series behavior in the case of the Fokker-Planck equation.

Perturbative power series in the UV
The limit w → 0 of Eq. (23a) can be understood by changing variables w → 1/z in Eq. (23). We obtain the following differential equation for F 1 (c 1 , z) where In the limit z → ∞ the dominant term in the previous expressions is O(z −1 ). It is straightforward to see that the fixed points correspond to the roots of the polynomial that multiplies the O(1/z) term in Eq. (43), The UV fixed points are given bȳ In analogy to what we did in the IR limit one can construct a perturbative expansion in the limit z → ∞. We consider where v 1,0 =c ± 1 . A similar power series expansion was discussed in the case of kinetic models undergoing Gubser flow in [29,32]. Inserting the UV series expansion in Eq. (46) into Eq. (42) we obtain the following recursive relation for the coefficients v 1,k where v 1,k ≡ 0 iff k < 0 and v 1,0 =c ± 1 .

Stability analysis and radius of convergence
When performing asymptotic expansions it is important to check the stability of the perturbative expansion. We address this issue for the UV power series expansion (46) by analyzing the linearized perturbations around it. Linearizing Eq. (43) around the UV fixed points, i.e. c 1 →c ± 1 + δc ± 1 , one finds the following evolution equation for the linearized perturbations in the z → ∞ limit The solutions are where µ ± 1 is an integration constant associated with the UV fixed points. The power law behavior of the linearized perturbations (49) is completely different from the behavior in the IR limit (29). We note that the exponents α ± 1 do not depend on the strength of the coupling, but they are uniquely defined by the location of the UV fixed points. In the limit z → ∞ the solution δc + 1 = µ + z −2.87516 (49) follows a power law decay. The linear perturbation δc − 1 = µ − z 1.67278 (49) monotonically increases. This divergence implies that the associated fixed point only admits a power series of the form (48) which can be extended analytically as it is shown below. Alternatively one may say that the initial condition µ − 1 for the fluctuation aroundc − 1 vanishes exactly and thus, the only possible expansion around c − 1 is a power series.
The formal power series in the z variable (46) has a finite radius of convergence which can be estimated via the Cauchy-Hadamard theorem 8 [90]. We proceed to calculate the radius of convergence of the power series expansion (46) by considering first the expansion around c − 1 . The case aroundc + 1 is analyzed separately at the end of this section.
The radius of convergence R − is calculated by numerically solving the recursive relation (47). In Fig. 4 (top panel) we show the coefficients of the UV expansion |v k | −1/k vs. the order of the expansion k for different values of (4π)η/s = {2, 4, 6, 8, 10} 9 . The coefficients |v k | −1/k stabilize for k ≥ 25 − 30. This result confirms that the power series (48) has a finite radius of convergence. We analyze its dependence on the value of η/s in Fig. 4 (bottom panel) 10 . This plots shows that R − de- 8 According to the Cauchy-Hadamard theorem [90], the radius of convergence R of the formal power series of a function f (z) around the point a (with b k , a ∈ C ) is given by The asymptotic value of η/s and θ 0 can be established via Eq. (27) which leads to the following identity where we used explicitly Eqs. (20). 10 We verify the extracted numerical value of the radius of convergence by increasing k up to kmax = 100. We found that the relative difference between the saturated bound |v k | −1/k for k = 25 and k = 100 was only 0.1%. pends linearly on the value of the shear over entropy ratio η/s. The empirical relation between these two quantities extracted from this plot is where the intercept b − > 0 is very small O(10 −8 ) and consistent with zero. The linear growth of the radius of convergence as a function of η/s is intuitively understood as follows: If η/s or, equivalently, the mean free path are large then the rate of collisions is small. As a result the Yang-Mills plasma will expand freely for a longer period of time. This finding might provide an explanation for the partial success of phenomenological models at intermediate scales of momentum larger than the typical temperature where the expansion is carried out in terms of a small number of scatterings [91][92][93][94][95].
In general the radius of convergence of a power series can be extended by analytical continuation [90]. The idea is to consider a power series of the form where both v 1,k and w 0 ∈ R, v 1,0 = c 1 (w 0 ) and in this case w 0 ≥ 0. The coefficients v 1,k are determined explicitly by plugging Eq. (51) into the ODE (23a) and finding the associated recursive relation analogous to Eq. (47) which is now evaluated at w = w 0 . By construction the series (51) has a finite radius of convergence R 0 so the series will diverge outside of |w − w 0 | < R 0 . One can thus expands again the function at another point w 1 located within a distance R 0 − δ (with δ>0) centered at w 0 and determine the new coefficients of the power series. Successively applying this approach one is extending the original power series significantly beyond its original radius of convergence R 0 . We apply this method to our solution by initializing the expansion at w 0 = 0 around c − 1 and θ 0 = 1 (i.e., extremely low values of η/s → 0).
The agreement between the exact numerical result and the analytical continuation of the perturbative expansion in the UV is truly remarkable, as shown in Fig. 5.

Resummation of the fluctuations aroundc + 1
The power law decay of linearized perturbations δc + 1 in Eq. (48), together with the UV stability analysis, strongly suggests a resummation scheme around the UV fixed pointc + 1 . The form of the equation in the z variable, see Eq. (42), resembles the corresponding results obtained in the case of the RTA Boltzmann equation for Bjorken and Gubser flows [31,32]. Following these approaches we consider the following UV transseries ansatz where the UV transmonomial is with α 1 = −2.86516. We note that, unlike the IR transseries where fluctuations decay exponentially close to the fixed point, the transmonomials ξ 1 (z) (53) are nonnegligible close to the UV fixed point. For instance, the first fluctuation term z α1 ≈ z −2.87516 is comparable to terms of order O(z −3 ). This implies that the complete solution nearc + 1 requires a double summation, over both the perturbative exponent k as well as the transmonomials (53) [31,32]. We also observe that the approximate UV solution (52) has a finite radius of convergence even after resumming the fluctuations aroundc + 1 [32]. Finally, in the UV regime it is not possible to implement IR transasymptotic matching (see Sect. III B) since the solution (52) is not valid for all values of z given the finite radius of convergence.
The coefficients v (n) 1,k in Eq. (52) are determined by inserting this ansatz into Eq. (42). We obtain the recursion relation For the UV transseries (52) the coefficients v

IV. TRANSASYMPTOTIC ANALYSIS: THE GENERAL CASE
In this section we generalize the transasymptotic analysis outlined in the previous section to the case when we consider the full vector of Legendre momentsc = (c 1 , · · · , c L ) T with L ≥ 1. We briefly outline how to generate IR and UV transseries solutions for the nonautonomous dynamical system (23). The techniques presented in the section were extensively discussed in [30][31][32], and we will focus on issues specific to the Fokker-Planck equation.

A. Transseries solutions in the IR
We seek to find multiparameter transseries solutions to the ODEs (22) by following the generic procedure developed by Costin [52,54]. In the large w limit and for the nonlinear ODEs (22) the Legendre moments admit IR expansions of the form [30,31] By mathematical induction one can show that by inserting the expansion (55) .
(57) We can check that the truncated l = 1 linear equation (28) is a special case of Eq. (56). In order to solve the linearized equation (56) we perform a transformation in terms of the pseudomodes δc δc =Ṽ (w)δc,Ṽ (w) : Thus, the equation of the pseudomodes δc is where the commutator [A, B] := AB − BA was introduced. One can take diag.(V ) = (0, · · · , 0), and the other components can be chosen in such a way that one can diagonalize the matrix W as follows Therefore, the solution of Eq. (59) is The Lyapunov exponents are the diagonal components of the matrixΛ (18a). These exponents govern the rate at which each mode relaxes towards its equilibrium value, w * l = [λ(l)] −1 . We observe that for larger value of l relaxation is faster, and that there is a clear hierarchy of scales w * 1 > w * 2 > · · · > w * L . For the FP equation the slowest non-hydro mode is c 1 which is proportional to the normalized shear viscous tensor component 11 in Eq. (13). The eigenvalue w ∈ R + of the matrixŴ is the anomalous dimension of the linear perturbation δc l . Their values w depends explicitly of the coupling constant as well as the truncation order L [30][31][32]; only in the limit L → ∞ do their values coincide with the exact underlying microscopic theory. We note that w is not a universal constant. Its value depends on the kinetic model, and is sensitive to nonlinearities encoded in mode-to-mode couplings.
The general transseries solutions of the nonlinear coupled ODEs (22) are constructed by rewriting these equations in terms of the pseudomodes basis, i.e., whereṼ (w) is given in Eq. (58). Under this transforma- 11 A common misunderstanding in the literature is to assume that the existence of this hierarchy of scales determines immediately the full set of the slowest degrees of freedom of the physical system. Recent studies [31,68,[96][97][98] have shown that this assumption is incorrect since nonlinear mode-to-mode coupling among moments plays a relevant role close to the forward attractor and thus, the determination of the slow invariant manifold of the dynamical system is not uniquely determined by the mere existence of a hierarchy of scales [31,32].
tion the equation of the pseudomodesc is where 1kc k . We have denoted matrices in the pseudomode basis byM ≡Ṽ MṼ −1 and the vec-torΓ ≡Ṽ Γ. Having determined the linearized pseudomodes (61) we can write the solutions of the pseudomodes as multiparameter transseries 12 where n = (n 1 , n 2 , · · · , n L ) is a vector where each component n i is a non-negative integer which labels the nonperturbative sectors of the pseudomodes, σ i ∈ C is the integration constant, and L ∈ N is the truncation order. The term ζ j in Eq. (65b) must match the linearized solutions (61) and thus, it determines the IR data The reality condition on the distribution function implies that the Legendre moments c l are real. Although the solutions for the pseudomodesc l (65a) are complex the reconstructed transseries solutions of c l (65c) are real since the following conditions are satisfied ∀j ≥ 1 [32] 1. λ(j) ∈ R + (with λ j given by Eq. (18a)) and λ(j) = λ(j ) for any j = j .
3. If β j ∈ C then there should be a complex conjugate ,k are determined by inserting the transseries solution of the pseudomodes (65a) into its corresponding Eq. (63). As a result one gets the following recursive relation

B. Transseries solutions in the UV
In this section we derive the UV transseries solutions by generalizing the results of Sect. III C 1. We show that the linear perturbations of the Legendre modes c l around the UV fixed point, which is a source, follow a powerlaw decay. We begin by discussing some features of the stability of the dynamical system of ODEs (22c) in the UV regime.
Consider the w → 0 limit of Eq. (22c) by first changing the variable w = 1/z. We get In the z → ∞ limit the dominant terms of the previous equation are O(1/z). In this limit one can determine the fixed points by solving the following equation The solutions to this equation provides a set of vectors c = (c 1 , ...,c L ) which determine the UV fixed points of Black and red dots correspond to lmax being either an odd or an even number respectively. Blue (c1 = 5)and brown (c1 = −5/2) dashed lines correspond to the UV fixed points of the original FPE.
the dynamical system. The solutions to Eq. (69) depend on the truncation order l max , so these are not necessarily real and we need to impose the reality conditions. Furthermore, the original FPE admits two UV fixed points for the moment c 1 , which have the property that the transverse and longitudinal pressures are minimized, c 1 = 5 and c 1 = −5/2 respectively. We call these configurations maximally prolate and maximally oblate respectively. The two UV fixed points act as bounds for the basin of attraction from above and below along the c 1 ray in the infinite dimensional space of the moments c l . Thus, any truncation scheme of the distribution function can be considered a good approximate solution of the Boltzmann equation iff the UV stability properties are reproduced approximately. In this case, the truncated dynamical system of ODEs should have at least two UV fixed points for the moment c 1 .
In Fig. 6 we plot the real solutions to Eq. (69) for the momentc 1 at the fixed points as a function of the truncation order l max ∈ [1,10]. We find that if l max is odd then the real solutions forc 1 come in pairs. As l max increases (for odd values) the UV fixed points approach the expectation for the original FPE,c 1 = {5, −5/2}. On the other hand, there is only one real solution forc 1 if l max is even. This observations also hold for the UV fixed points of higher order momentsc l . Thus, the standard idea of obtaining more accurate approximate solution to the Boltzmann equation by adding more moments to the distribution function is not necessarily correct. A good truncation scheme must reproduce the flow structure of the original Boltzmann equation in the IR and UV regimes. Similar findings were reported for the RTA Boltzmann equation [31,32]. In the rest of this section we will consider the case that l max is odd.
In the limit z → ∞ we linearized Eq. (68) around a particular UV real fixed pointc, namely c →c + δc, which lead to the following equation for the linear perturbations where the matrixX and X are defined in Eq. (16a). Upon inspection we find that Eq. (48) is a particular case of this linearized equation for the perturbations. We introduce the linear pseudomodes δc = Uδc so that U diagonalizes the matrix T(c), This procedure gives the equation for the linearized pseu- Its solution is This result shows that close to the UV fixed points the linearized perturbations have a power-law behavior which is opposite to their IR counterparts (61). In general the eigenvaluet i , the integration constant µ i as well as the matrix element U ij are complex numbers. The reality condition of the linearized modes is satisfied since the coefficients of the polynomial eigenvalue problem of the matrix T are real and thus, their complex eigenvectors and eigenvalues come in complex conjugate pairs. When summing over the complete eigenbases of pseudomodes the linearized modes are real as expected. Moreover, what happens in the vicinity of any UV fixed pointc along the rayc i is encoded by the real part of the eigenvalues of the matrix T, namely R(t i ). In the limit z → ∞ the linearized perturbation δc i decays rapidly as a power law if R(t i ) < 0 while δc i increases when R(t i ) > 0.
In order to generate approximate UV solutions in the limit z → ∞ limit it is convenient to rewrite the equations of motion in terms of the pseudomode basis c i = U ij c j (with U given by Eq. (71)). In terms of the pseudomodes the nonlinear ODEs (68) read In the previous equation it must be understood that The transseries ansatz solution is built up based on the solutions of the linearized perturbations (73). These solutions can be taken as the transmonomials of the transseries which carry out the information about the non-perturbative contributions. There-fore, in the limit z → ∞ the transseries ansatz is given by [31,32] where m ∈ N L 0 , µ i ∈ C is the integration constant and L ∈ N is the truncation order. The UV data is determined by matching the transmonomial ξ j with the linearized solutions (73) and thus, the anomalous dimensions α j entering in the previous expression are the set of eigenvaluest of the linearization matrixT, Eq. (71), evaluated at the UV fixed pointsc. A similar transseries ansatz was proposed for the RTA Boltzmann equation in systems undergoing Bjorken [31] and Gubser flow [32]. The solutions of the Legendre moments in terms of the pseudomodes are given by The coefficientsũ m l,k entering in the transseries ansatz are determined entirely by inserting Eq. (75) into Eq. (74) which leads to the following recursive relation The normalization condition is chosen to beṽ 1 i,0 = 1 and the coefficientsṽ 0 i,0 = U ijcj since one expands around a given UV fixed pointc. In Eq. (78)ṽ (m) k denotes a coefficient proportional to ϕ m k . We emphasize that the transseries ansatz (75) is a convergent series with a finite radius of convergence. The rate of convergence depends on the anomalous dimension α i ≡ R(t i ) governing linearized perturbations δc i (73). The lessons learned from the l = 1 truncation studied in Sect. III C (see also Refs. [31,32]) show that the power law transseries solutions (76) can be constructed iff R(t i ) < 0. When R(t i ) > 0 the best option is to perform an analytical continuation while cancelling the transmonomial contributions in Eq. (75) by setting the integration constant µ i ≡ 0.

V. UNIVERSAL ASPECTS OF ATTRACTORS FOR DIFFERENT KINETIC MODELS
In the previous sections we analyzed the behavior of the solutions of the FPE in both the UV and IR limits. In this section we study universal aspects of the attractors in the FPE as well as the RTA Boltzmann equation and the AMY kinetic theory [60].
Following the approach of Refs. [99,100] we consider a distribution function which is squeezed along the beam direction at early time 13 . In this situation the initial longitudinal pressure vanishes exactly and thus, this configuration determines the pullback attractor of the distribution function at early times [30,32]. Under this condition the initial phase space distribution can be modeled as [100] The normalization constant is chosen such that the initial energy density per unit rapidity per transverse area is constant The initial distribution function determines the initial conditions for the Legendre moments Previous studies [86,99] have shown that universal behavior of the numerical solutions can be by analyzed in terms of observables that are less sensitive to the initial conditions. Following Refs. [99,100] we study the following observable We will analyze this observable as a function of the variablew = τ T (τ )/ [(4π)η/s]. In terms of this variable the results are insensitive to the strength of the coupling. Eq. (82) has two interesting limits in the forward and pullback attracting regions. The former is described by a few terms of the non-hydrodynamic expansion [86], while the latter is determined by expanding around the UV fixed point where the longitudinal pressure vanishes [99], 13 For the RTA Boltzmann equation the form of the initial distribution function does not play a major role. However different processes of Yang-Mills plasmas, i.e. elastic and inellastic interactions, affect the parametrization of the momentum distributions [101,102]. In our approach these effects do not play a role at the level of the moments since the free streaming expansion dominates at early times over the collision rate of the FPE as discussed in Sect. IV B.
where the constant C ∞ is determined from Eq. (51) in Ref. [100]. Alternatively one can determine this constant by fitting to the numerical data. We verified that both methods lead to the same approximate value for this constant, see also [99,100]. We numerically solved the evolution equations for the Legendre moments, Eqs. (14) and (21) respectively, together with the conservation law (14a). The solutions for the Legendre moments c l were obtained by truncating the expansion at l < l max , setting c l ≡ 0 for l > l max . This method converges rapidly and no sizable deviations are observed for l max 35 for both the RTA Boltzmann and FPE equations. In the case of the Yang-Mills AMY theory we used publicly available numerical results 14 [99,103,104].
In Fig. 7 we show numerical solutions of E as a function ofw for the RTA Boltzmann equation (dashed green line), FPE (dashed blue line) and AMY kinetic theory (red line). The general behavior of the different kinetic models is quite similar. They all exhibit a smooth transition from early-time free streaming to universal late-time hydrodynamic behavior. Deep in the UV regime the expansion (dashed gray lines) dominates over collisions regardless of the underlying collision kernel, and the pullback attractor is determined by free streaming expansion. A more surprising fact is that the interactions modifies the expansion for similar values of w ≈ 0.08, i.e. at a large Knudsen number Kn ≈ 12.5, independent of the underlying microscopic theory. We note that the existence of a pullback attractor is manifest only in the w variable. In the original proper time variable τ the limit τ 0 → 0 is not well defined because the ODEs have an essential singularity at this point [31,32].
In the UV regime the main difference among the kinetic models depicted in Fig. 7 is the values of the constant C ∞ = {0.98, 0.91, 0.89} for the AMY, FPE, and RTA models, respectively. These values allows us to quantify the difference among the predictions of RTA and FPE respect the YM-AMY model. The former differs by ∼ 10% while the latter ∼ 8%. It is somewhat surprising that a relatively simple approximation such as RTA captures many aspects of the non-equilibrium dynamics compared with YM-AMY kinetics. This observation was also noticed previously in [99,100].
In the IR regime we observe in Fig. 7 that all kinetic models reach the non-hydrodynamic behavior and thus, the forward attractor is encoded by the late-time nonhydrodynamic expansion. All the kinetic models reach the asymptotic hydrodynamic gradient expansion around w ≈ 2, and thus Kn ≈ 0.5 so the deviations from equilibrium are quite large. This finding implies that one can reach out the non-hydrodynamic behavior without having local thermal equilibrium. Given that we find forward attractors for a number of different kinetic equations, one may ask how general this feature is. In App. E we discuss a sufficient set of mathematical conditions which ensure its existence for weakly coupled boost invariant systems with highly nonlinear collision kernels beyond those studied in this work.

VI. DISCUSSION AND FINAL REMARKS
In this work we studied the hydrodynamization of a boost invariant system of gluons described by kinetic equation derived from QCD in the small angle, diffusive, approximation. We demonstrate that the physics of the Fokker-Planck equation can be recast in terms of a set of non-linear ODEs for the Legendre moments of the one-particle distribution function. We show that these kinetic equations admit transseries solutions in both the UV and IR regimes. These findings extend previous results obtained in the context of the Boltzmann equation in the relaxation time approximation [30][31][32].
We employ techniques from the theory of nonlinear dynamical systems. Applying these methods to the kinetic equations we investigate the stability properties of linearized perturbations around the IR and UV fixed points. We analyze the the emergence of transseries solutions and show that their functional form is a rather natural conse-quence of the stability properties of the nonlinear ODEs around the UV and IR fixed points, respectively.
In the IR regime we prove that the solutions of the moments equations are multiparameter transseries. The associated transmonomials are built up using the behavior of linearized perturbations around the IR fixed point. These perturbations involve the product of an integration constant, an exponentially decaying term and a power law term. The exponentially decaying term determines the rate at which linearized perturbations decay near the IR fixed point. This decay time is controlled by the product of the Lyapunov exponent and the inverse Knudsen number. These terms play an analogous role to 'instanton-like' contributions in QFT. The IR transseries effectively resums non-perturbative dissipative contributions that are absent in the usual perturbative gradient expansion. These results strongly suggest that a consistent formulation of non-hydrodynamic theories in farfrom-equilibrium regimes requires the inclusion on nonperturbative physics.
In the IR regime the constitutive relations of each Legendre moment c l contain the non-perturbative information encoded in the multiparameter transseries. The result suggests a non-perturbative dynamic renormalization scheme which goes beyond standard linear response theory. In this theoretical framework, transport coefficients can be defined in far-from-equilibrium regimes when the Knudsen number is large. At each order of the pertubative IR expansion, transport coefficients are dynamically renormalized by non-perturbative corrections in the transseries, a method known as "transasymptotic matching" in the resurgence literature. This approach has a nice interpretation in terms of an RG flow.
From the physical standpoint of view, the renormalized transport coefficients depend on the rheology of the fluid. They describe the relaxation of dissipative coefficients to their values dictated by the linear response approach. Thus, the fluid experiences transient non-Newtonian behavior prior to hydrodynamization.
On the mathematical side, we explored the relation between dynamical systems and RG flows. Indeed, it is known that the field theoretical RG flows may be discussed in the framework of autonomous dynamical systems [55,56]. Our work shows that this link can be extended to non-autonomous systems that hydrodynamize, where time plays the role of scale parameter. In other words, a non-autonomous dynamical system can be considered as an RG flow equation provided there is a slow invariant manifold, i.e. when the long-time physics is dominated by slow degrees of freedom [32]. In the case studied here and in [30][31][32] the invariant manifold is shown to exist if the IR perturbations go to zero in the long-time limit, a fact that points to the existence of a forward attractor. In particular, there are bounded solutions at w → ∞ or τ → ∞ for a fixed initial time that approach the equilibrium point. For expanding systems such as Gubser flow [29,32], the system is completely perturbative so that full hydrodynamization does not occur. Therefore, there is no slow invariant manifold, which in turn explains why the RG-flow paradigm for the transport coefficients is no longer available.
In the UV regime the nonlinear ODEs for the moments admit transseries solutions with a rather different behavior compared to their counterparts in the IR limit. These solutions are power series with a finite radius of convergence. The radius of convergence R grows linearly with the shear viscosity over entropy ratio η/s. The linear relation between R and η/s provides a rather simple explanation of previous numerical findings [72,75,76,[105][106][107][108][109][110][111] where it was noticed that initial conditions close to the maximally oblate UV fixed point converge slowly to the forward non-hydrodynamic attractor than the ones located in the vicinity of the maximally prolate UV fixed point. The latter is a saddle point for any viable truncation (e.g. if the truncation bound for l is odd) and it is found that there are a set of solutions for the Boltzmann equation that remain bounded near this point. In contrast, the maximally oblate UV fixed point is a source from which flow lines are streaming away in w 0 → 0 for a fixed w, meaning that there are no attracting regions to probe in the past of the dynamical system around this fixed point.
We show that the validity of the power series expansion can be extended by analytic continuation. As an alternative, we also introduce a new resummation scheme which accounts for the non-perturbative physics of the power-law behavior of the fluctuations around the UV fixed points. We note that the UV stability analysis also unveils a non-trivial aspect of truncating the moment expansion in relativistic kinetic theory: In order to preserve the UV structure of the original kinetic equation for Bjorken flow we have to include an odd number of Legendre moments L = 2n + 1 moments in the n → ∞ limit.
We discussed the general properties of the IR and UV regimes of the Fokker-Planck equation in comparison to the RTA Boltzmann and AMY kinetic equations. The pullback attracting region is entirely determined by the free streaming limit where the longitudinal extent of the distribution function is negligible. The attracting region in the long time limit is corresponds to hydrodynamization and can be characterized by a few terms in the gradient expansion. There are, however, minor between different models. We find that the RTA and AMY differ by at most 10% while the latter disagrees by up to 8% with the FPE result.
The techniques presented in this work are not necessarily restricted to far-from-equilibrium systems undergoing longitudinal boost invariant expansion. Indeed, these methods can be applied to describe non-equilibrium dynamics in other physical systems of relevance such as cold atoms systems, condensed matter and cosmology.
The remaining momentum integral in the RHS of Eq. (A2) reads as where By equating Eqs. (A6) and (A7) together with Eqs. (A5) into the RHS in Eq. (A2) we finally get with θ −1 0 = 5 8 π 5 λ 2 Y M l Cb . Thus, by putting together the LHS (A3) and RHS (A9) of Eq. (A2) we get the following evolution equations for the Legendre moments c l After some redefinitions of the variables and writing them in a matrix form one gets Eqs. (14).

Some useful identities of the Legendre polynomials
In the previous section we used explicitly the following identities of the Legendre polynomials P n (x) = 2 n n k=0 n k n+k−1 2 ∂ ∂ cos θ p sin 2 θ p ∂ ∂ cos θ p P 2 (cos θ p ) = −2 (2 + 1)P 2 (cos θ p ).
transseries solutions (see also Refs. [30][31][32]). Naively, one might think that a better agreement between the transseries and the numerical results can be obtained when adding more transmonomials and/or higher orders in the transseries. Nevertheless, this is not the case since in our case the radius of convergence of the transseries is finite [31,32]. Consider the truncated IR asympotic expansion of the Legendre mode c 1 , i.e., where K is the order of the truncation. The IR perturbative expansion (24) is asymptotically of Gevrais-1 class (see discussion in Sect. IIID of Ref. [31]). Thus, the asymptotic form of the coefficients u 1,k entering in Eq. (24) are where M (θ 0 ) = M 0θ β1 0 (M 0 ∈ R) is an overall factor that depends on the angleθ 0 defined in the Borel plane. The function K op (w) which optimizes the error between the asymptotic series and the exact solution can be calculated by evaluating the convergence rate as follows where • is the Floor function. Thus, the induced error R K (w) in the truncated expansion (B1) is where we used the Stirling's formula for the last line. Clearly, this estimate depends on the asymptotic behavior of |u (0) k | (B2). In that case if k < K op the error is small when |u (0) k | behaves like in Eq. (B2) and depends explicitly on w. Thus, when adding higher orders k > K op to the truncated expansion (B1) the mismatch between transseries and the exact solution increases. The situation gets more worrisome when solving the general dynamical system (14). In this general case, the leading order perturbative contribution to c (τ ) = O(τ − ) [31]. As a result, the disagreement between the multi-parameter transseries and the numerical solutions of the Legendre moments c (τ ) is larger and expected if > K op . This was precisely what some of us observed and extensively discussed in Refs. [31,32].
Appendix C: IR Transseries solutions in the proper time τ variable In Sects. III and IV we studied in detail the resurgent properties of the nonlinear ODEs (23) in terms of the w variable. However, one can also analyze the original dynamical system in terms of the variable τ as it was done for the RTA Boltzmann equation [30,31] for the detail. For completeness we present this analysis in this section.
The IR fixed points of the Legendre moments arec = 0 ∀ > 0, namely c (τ ) → 0 as τ → +∞. From Eq. (14a), one obtains the exact formal solution of T (τ ), i.e., where σ T ∈ R + is the integration constant with dimensions of energy andτ := (σ T τ ) 2/3 15 . Notice that thatτ is dimensionless and guarantees the scale invariance of c in the w-coordinate. In addition, in the second line we expanded 15 The integration and the series expansions are commutative with each other up to the value of integration constant determined from a given initial condition in general. In this sense, the value of σ T in a transseries is individually determined for each fixed points from a given initial condition.
It is straightforward to find that c (τ ) = O(τ − ) by considering the asymptotic behavior of the solutions of Eqs. (14). For instance, one obtains that c 1 andT behave asymptotically as follows As we proceeded in Sects. III and IV we rewrite the nonlinear ODEs (C2) in a matrix form, i.e., where c(τ ) = (c 1 (τ ), · · · , c L (τ )) , The matrices entering into the previous expressions are defined in Eqs. (16a) and (19). In order to build the transmonomials we linearized the dynamical system (14) around the IR fixed point, ∂F(c,τ ) ∂c δc (τ )  .
In order to solve find the behavior of the fluctuations close to the IR fixed points it is needed to solve the eigenproblem for the matrixW. Following Costin's prescription [54] we introduce the matrix V as follows where n ∈ N L 0 , σ i ∈ C is the integration constant, and L ∈ N is the truncation order of , namely 1 ≤ ≤ L.

Appendix D: On the existence of the Lyapunov functional
In this section we present a proof of the existence of a Lyapunov function for the dynamical system of ODEs (22). Roughly speaking, Lyapunov functions are positive definite functions which are monotonically decreasing along the trajectories of the flows in phase space. These functions are used to determine the stability properties of ODEs and PDEs. In the RG flow approach to dynamical systems the Lyapunov function plays an analogous role of the c function in QFT and/or the dynamical effective potential in non-newtonian mechanics. For a non-autonomous dynamical system, the existence of the Lyapunov function allows us to identify the dynamical system as a RG flow from the global point of view.
In our approach the existence of the Lyapunov function is inferred directly from the dynamical system of ODEs (22). First, we promote this non-autonomous dynamical system to an autonomous one of one dimension higher by introducing an ODE for w in terms of a new flow time ρ as follows where in the RHS of the previous expressions we introduce a positive definite differentiable function V. It is straightforward to show that the function V decreases monotonically, i.e., dV(c(ρ), w(ρ)) d log ρ = L i=1 dc i d log ρ . ∂V(c(ρ), w(ρ)) ∂ c i (ρ) + dw(ρ) d log ρ . ∂V(c(ρ), w(ρ)) ∂ w(ρ) , |β i (c(ρ), w(ρ))| 2 + |β w (c(ρ), w(ρ))| 2 ≤ 0 .

(D2)
Therefore, V satisfies the properties required to be a candidate for the global Lyapunov function of the dynamical system (22) (see Eq. B9 in Ref. [31]). Although we have proven its existence for the studied case here, it is in general extremely difficult to calculate the exact Lyapunov function (cf. [29]) and this is beyond the scope of this work.