Dynamical systems and nonlinear transient rheology of the far-from-equilibrium Bjorken flow

In relativistic kinetic theory, the one-particle distribution function is approximated by an asymptotic perturbative power series in Knudsen number which is divergent. For the Bjorken flow, we expand the distribution function in terms of its moments and study their nonlinear evolution equations. The resulting coupled dynamical system can be solved for each moment consistently using a multi-parameter transseries which makes the constitutive relations inherit the same structure. A new non-perturbative dynamical renormalization scheme is born out of this formalism that goes beyond the linear response theory. We show that there is a Lyapunov function, aka dynamical potential, which is, in general, a function of the moments and time satisfying Lyapunov stability conditions along RG flows connected to the asymptotic hydrodynamic fixed point. As a result, the transport coefficients get dynamically renormalized at every order in the time-dependent perturbative expansion by receiving non-perturbative corrections present in the transseries. The connection between the integration constants and the UV data is discussed using the language of dynamical systems. Furthermore, we show that the first dissipative correction in the Knudsen number to the distribution function is not only determined by the known effective shear viscous term but also a new high energy non-hydrodynamic mode. It is demonstrated that the survival of this new mode is intrinsically related to the nonlinear mode-to-mode coupling with the shear viscous term. Finally, we comment on some possible phenomenological applications of the proposed non-hydrodynamic transport theory.


I. INTRODUCTION
Recent measurements of multiparticle cumulants in d +Au , p+Pb, p+Au, and 3 He+Au systems have provided a compelling evidence for collective flow in heavy-light ions collisions [1][2][3][4][5]. Similar observations had been made previously in nucleusnucleus collisions (cf. Refs. [6,7] and references therein). Altogether, these measurements seem to imply that the collective evolution of a nuclear fireball created in these experiments can be understood in terms of hydrodynamics.
The phenomenological success of hydrodynamic models in producing an accurate description of these extreme experimental situations have raised some questions as per the validity of the assumptions that hydrodynamics relies on. For instance, the local thermal equilibrium hypothesis has been usually thought of as a necessary and sufficient criterion for the applicability of the fluid-dynamical equations of motion. Different theoretical studies [8][9][10][11][12][13][14][15][16][17][18][19][20] have shown that the equilibrium condition seems to be too restrictive, and that surprisingly no inconsistency arises when hydrodynamics is used for interpreting a system sitting far from equilibrium. These findings have repeatedly pointed out to the possibility of generalizing relativistic hydrodynamic theories for systems with the local thermal equilibrium removed [21][22][23]. Although the idea of employing fluid dynamics in non-thermal equilibrium physics is not entirely new [24,25], little progress has been made to formulate it from first principles.
In recent years we have learned more about certain generic properties of far-from-equilibrium hydrodynamics that follows suitably from the nonlinear nature of fluid dynamics. The breakdown of the naive perturbative asymptotic expansion is, therefore, imminent, meaning that a series ansatz fails to converge and accordingly it cannot be a full-blown solution to the underlying nonlinear differential equations. This has been long known in mathematics [26]. An example of this fact in hydrodynamics was given in [27] where the gradient expansion of the energy-momentum tensor of a conformal fluid was shown to be divergent. The divergences of the energy-momentum tensor ought to beg for the existence of additional degrees of freedom called nonhydrodynamic modes. Since then more examples have appeared in generic setups in both strong and weakly coupled regimes [21,[27][28][29][30][31][32][33][34][35][36][37][38][39][40].
A more relevant and important question to ask in this context has to do with how non-hydrodynamic modes can affect the transport properties of the system and what their possible phenomenological significance is. We addressed this question in our previous work [41] by investigating the nonlinear dynamics of a far-from-equilibrium weakly coupled plasma undergoing Bjorken expansion. Our work first treated the equations governing the time evolution of momentum moments of the one-particle distribution function as a coupled nonlinear dynamical system, the solutions of which are found to be multi-parameter transseries with real integration constants σ . These transseries carry exponentially small (non-perturbative) information in terms of the Knudsen Kn and inverse Reynolds Re −1 numbers which both play the role of expansion parameters in the transseries, whereas the relaxation time τ r characterizes the size of exponential corrections. A linear combination of these moments defines the constitutive relations at every order in the perturbative expansion whose coefficients are indeed related to the transport coefficients. Then summing over all exponentially small factors and other monomials of the expansion parameters in the transseris renders an effective dynamical renormalization of the transport coefficients away from the values obtained at equilibrium using linear response theory. It should be noted that in this approach, 'dynamical renormalization' means that there exists a positive-definite monotonically decreasing functional also known as Lyapunov functional in the stability theory, along every flow line in the space of moments, that can be easily derived from our dynamical system.
One advantage that comes with studying dynamical systems is the ability to have full control over the late-time (or IR) as well as early-time (or UV) behavior of the flow lines. We should keep in mind that the the alternative techniques such as Borel resummation heavily rely on the asymptotics of the divergent series. However, in a dynamical system of evolving moments, the resurgent multi-parameter transseries is directly found as a formal solution without the need for going to Borel plane where UV information is completely lost. The Borel singularities in this sense only encode IR data such as the the size of exponential corrections which in our approach corresponds to the eigenvalues of the linearization matrix commonly known as Lyapunov exponents around the asymptotic hydrodynamic fixed point. Therefore, the theory after Borel resummation is only asymptotically equivalent to the original theory and drawing any conclusions on the fate of solutions of the original theory in the UV from a resumed series is not in general possible.
The organization of this paper is as follows. We extend and generalize our previous findings [41] to include higher nonhydrodynamic modes in Sects. II-III while providing technical details of our derivations in the Appendices. In Sect. IV, we explain further the deep relationship between dynamical systems and the resurgence theory in the time-dependent (nonautonomous) system at hand, which was first discussed in the context of relativistic hydrodynamics in Ref. [42]. Finally, in Sect. V it will be reported that a new high energy non-hydrodynamic mode exists whose first-order perturbative asymptotic term goes like ∼ (τT ) −1 , which obviously displays the same decay as the Navier-Stokes (NS) shear viscous tensor component. The survival of this new non-hydrodynamic mode is due to the nonlinear mode-to-mode coupling with the effective shear viscous tensor. Some of the quantitative and qualitative properties of this new mode will also be discussed. As a last couple of remarks, we should remind the readers that for the sake of keeping the paper self-contained, we have included a small number of mathematical definitions from dynamical systems in App. B that will make the read easier. Also, as mentioned above, we will adopt the field theory language when talking about the early-time and late-time regimes, which then can be interchangeable with UV and IR, respectively.

II. KINETIC MODEL
The statistical and transport properties of the weakly coupled systems are usually described by kinetic theory. Within this approach, it becomes important to understand the evolution of one-particle distribution function from the Boltzmann equation [43][44][45]. Fluid-dynamical equations of motion arise as a coarse-graining process of the distribution function where the fastest degrees of freedom are integrated out. The approximations made in the coarse-graining procedure do not necessarily lead to the same macroscopic evolution equations [46,47], and nor do they in general draw an accurate picture of the physics when compared against exact numerical solutions [10-13, 48, 49]. The common approaches to finding approximate solutions of the Boltzmann equation are the Chapman-Enskog approximation [50] and Grad's moment method [43,51]. In the Chapman-Enskog method, the distribution function is expanded as a power series in the Knudsen number around its equilibrium value, which eventually leads to a divergent series [52,53]. In Grad's moment method, the distribution function is expanded about the equilibrium state in terms of an orthogonal set of polynomials, and the coefficients of such expansion are average momentum moments [51]. The second approach, however, turns out to be more convenient when dealing with far-from-equilibrium backgrounds [20,[54][55][56][57].
In this work, we also adopt the latter method and study the full nonlinear aspects of the moments' evolution equations. We carry out our research plan by considering a weakly coupled relativistic system of massless particles which undergoes Bjorken expansion [58] with vanishing chemical potential. We simplify our problem by assuming the Boltzmann equation within the relaxation time approximation (RTA-BE).
For the Bjorken flow we use the Milne coordinates x µ = (τ, x, y, ς ) (with the longitudinal proper time τ = √ t 2 − z 2 and pseudorapidity ς = arctanh (z/t)), and the metric is g µν = diag. −1, 1, 1, τ 2 . The time-like normal vector identified with the fluid velocity is taken to be u µ = (1, 0, 0, 0) (with u µ u µ = −1) and the spatial-like normal vector pointing along the ς -direction is l µ = (0, 0, 0, 1) (with l µ l µ = 1). The Bjorken flow symmetries reduce the RTA-BE to the following relaxation type equation [59] ∂ τ f τ, p T , p ς = − 1 τ r (τ) f τ, p T , p ς − f eq. (−u · p/T ) , where p T is the transverse momentum, p ς denotes the momentum component along ς -direction, and τ r represents the relaxation time scale. Also, f eq. is the local equilibrium Jüttner distribution which without loss of generality is chosen to be of Maxwell-Boltzmann type, i.e., f eq (x) = e −x . We recall that τ r sets the time scale at which the system relaxes to its thermal equilibrium. We shall consider models whose relaxation time is a power law in the effective temperature, namely For ∆ = 1 means that the constant θ 0 is dimensionful, while for ∆ = 0 the theory is conformally invariant. 1 For pedagogical purposes we perform explicit calculations for the case of ∆ = 0 in the main body of the paper, while the gist of results for the general case are discussed in Appendices D, E and F. The RTA-BE (1) can be solved exactly [59]. In what follows, we will take a distinct approach in which the mathematical problem of solving the Boltzmann equation is recast into seeking solutions to a set of nonlinear ODEs for the moments. For the Bjorken flow we propose the following ansatz for the single-particle distribution function [41] f τ, p T , where p τ = p 2 T + (p ς /τ) 2 is the energy of the particle in the comoving frame, L n and P 2l denote the generalized Laguerre and the Legendre polynomials, respectively. The ansatz (3) allows us to study hydrodynamization processes [41]. Furthermore, the nonlinear relaxation of the low energy (n = 0) as well as the high energy tails (n > 0) of the distribution function are better understood when f (τ, p T , p ς ) is expanded in terms of orthogonal polynomials [61,62].
The moments c nl are read directly from Eq. (3) as follows 2 where Γ(n) is the Gamma function and the momentum average of any observable O(x µ , p µ ) weighted by an arbitrary distribution then have that the hydrodynamic equilibrium (asymptotic IR fixed point) is given by c eq. nl = δ n0 δ l0 . For the Bjorken flow the energy momentum tensor T µν = p µ p ν is [55, 65-69] 3 T µν = ε u µ u ν + P L l µ l µ + P T Ξ µν , with the projector operator Ξ µν = g µν + u µ u ν − l µ l ν which is orthogonal to both u µ and l µ . The energy density ε, the transverse and longitudinal pressures (P L and P T respectively) can be written in terms of the moments c nl [41] ε = (−u · p) 2 = 3 It follows that ε = 2P T + P L from these expressions. We also mention the nonnegativity of pressure components, P T , P L ≥ 0, sets the physical range for c 01 as −2.5 ≤ c 01 ≤ 5.
These bounds are satisfied by the exact solution of the RTA-BE (1) but it is not expected to be satisfied by a particular truncation scheme for the distribution function. The energy-momentum conservation together with the Landau matching condition for the energy density imply c 00 ≡ 1. The only independent (normalized) shear viscous componentπ ≡ τ 2 π ς ς /ε for the Bjorken flow is proportional to the moment c 01 [41] In a previous work of some of us [41], it was shown that the time evolution of the temperature T and the moments c nl with n ≥ 0 and l ≥ 1 is described by an infinite number of coupled nonlinear ODEs as follows where the coefficients are given by , , The hierarchy of equations in (9) constitute a dynamical system where moments of different degree n and l mix amongselves in a nontrivial way. The nonlinear nature of the RTA-BE is manifest in the set of ODEs (9b) by the mode-to-mode coupling term ∼ c nl c 01 . Nonetheless, one observes that the low energy modes c 0l decouple entirely from the high energy ones c nl with n > 0. As a result, the time evolution of the energy momentum tensor is fully reconstructed from the solutions of the temperature and the modes c 0l . Notice that one cannot deduce the same thing about the high energy modes whose evolution receive major contribution from the lower energy modes. Furthermore, the high energy moments shall play a role in the stability of the system and its convergence to its asymptotic thermal state, which is the subject of discussion in Sect. V.

III. TRANSSERIES SOLUTIONS AND COSTIN'S FORMULA
In this section, we construct the transseries solution to the dynamical system in (9). Once one reduces to the Boltzmann equation to a dynamical system, next immediate step would be to construct an exact form of the transseries around each asymptotic fixed point and obtain all the coefficients recursively from the evolution equations. In what follows, we will attempt to build the general form of the exact transseries solution by slightly modifying the original set of ODEs. For technical reasons, we also start with a truncation of the dynamical system at 0 ≤ l ≤ L and 0 ≤ n ≤ N, hence the original Boltzmann distribution will be reproduced by taking the limit N, L → ∞. Since we are interested in building the solutions starting at late times, we will assume the convergence at infinity, i.e., c(τ) → 0 for τ → ∞.
Let us prepare the set of ODEs in the dynamical system (9) in the following asymptotically linearized form: whereΛ and B are constant matrices, c is an ((N + 1) × (L + 1) − 1)-dimensional vector given by c = (c 01 , . . . , c 0L , c 10 , c 11 , . . . , c 1L , . . . , c N0 , . . . , c NL ), and A is a constant vector. Here, we have defined a new time coordinate as w = τT (τ) which behaves asymptotically like w ∼ τ 2/3 at late times. We remind that rank(Λ) = rank(B) = (N + 1) × (L + 1) − 1 =: I. Furthermore, we assumeΛ is a diagonal matrix proportial to the unit matrix. To construct the transseries solution, we suitably diagonalize B by defining an invertible matrix U such that These transformations cast Eqs. (11) in the following form We hereby call the components of the vectorsc pseudomodes due to being generally complex-valued, thus not physical. This makes the components of the matrix U complex-valued as well. But the inverse transformations ofc achieved by the action of U −1 always yield real vectors that eventually contribute to the observables in our theory. The classical asymptotics beyond naive asymptotic power series expansion can be carried out with the help of a transseries ansatz first put forward by O. Costin in Ref. [26]. In that seminal work, the author proves that the set of ODEs written in the prepared form (14) has the exact transseries solutioñ where I = dim (c), m ∈ N I 0 is an integer vector, and the dot denotes the inner product between any two vectors. The real numbers σ i are going be referred to as "integration constants" throughout this work as they would really symbolize the constants to be obtained if we were to integrate the ODEs in Eqs. (14). Here, for simplicity we have defined E (m) k (w) that stand for the basis of transmonomials (i.e., the exponential factors and the fractional powers w −1 ).
The transseries data such as the coefficientsũ (n) i,k , the Lyapunov exponents S i and the anomalous dimensionsb i can be recursively determined by the evolution equations up to normalization of σ i . Without loss of generality, we pick the normalization fixed by u (m) j,0 = δ i j for m j = δ i j . This leaves no room for ambiguity in determining the coefficientsũ (m) i,k . It is noteworthy that the transseries solution of c i (w) can be reproduced by the inverse transformations of (13), and one can find that c i (w) has essentially the same transseries asc i (w) due to the fact that the matrix U acts on the index i (mode number) inũ (m) i,k only. Although the transseries ofc i (w) given by (15) is generally complex-valued due tob i andũ (0) ik taking values in complex numbers, we can still recover a real transseries solution for the moments c i (w). The important fact is that ifb i for some i is complex, it is always accompanied by its complex conjugate counterpart, i.e.,b i =b * j for some j, where the associated coefficients satisfyũ j,k under a certain normalization of the eigenvectors ofB in such a way that U −1 Ii = 1 for every i where again I is the rank of matrixB. Therefore, the reality of c i (w) is guaranteed if the following conditions are satisfied: A. Evolution equations: N = 0 case In terms of w, Eqs. (9) can be reduced to the following non-autonomous dynamical system: The second equation is valid for any l > 0. Since temperature T is now sort of washed away from Eq. (22b), we shall solve these ODEs for c nl by the mathematical techniques discussed in Refs. [26,70]. These tools are well known in the context of resurgence theory. A curious reader is invited to check e.g., Ref. [71] for technical details, Ref. [72] for a nice introduction to the subject, and Ref. [73] for the summary of its applications to quantum mechanics and quantum field theories). We will afterwards substitute c 01 in Eq. (22a) to solve for T . It is convenient to cast Eqs. (22a)-(22b) in the following familiar form where Here, the index of vectors and matrices runs over i = 1, . . . , L. To directly apply Costin's formula, one has to diagonalize B in Eq. (23) using the transformations in (13) to get the equation We can easily see that c 0i can be reproduced byc and U as As is seen in Eq. (31), from now on whenever N = 0, we may drop the index 0 in c 0l for ease of notation. Notice that the asymptotic expansion ofc i , c i both start at order O(w −1 ) asymptotically, and due to this we haveũ i,0 = 0. Before trying to solve the dynamical system in (23), it would be instructive to first linearize the r.h.s. of it to obtain the so-called IR data S i ,b i in Eq. (18) of the transseries ansatz. These are going to be the input data obtained uniquely from the profile of the solution of the linarized equation around the asymptotic fixed point of the dynamical system at w → ∞. To do so, we expand the factor where c 1 (w) = (c 1 (w), 0, . . . , 0) . The coefficients of the lowest power of w −1 , i.e., u i,k , can be found from this equation, and one can readily see that c i=l = O(w −l ). In particular, Now, we may proceed to linearize the equation (33) by substituting c i →c i + δ c i withc i = ∑ k u (0) i,k w −k , and expanding to first order in δ c i Using the matrix U, the solution to the linearized equation is found to be given by where σ i is the integration constant and δc = Uδ c. Therefore, we have the IR data (Lyapunov exponent and anomalous dimension) as Finally, we substitute the ansatz (15) in Eq. (32) and then insert the resulting transseries for c in Eqs.(30)-(31) to find the recursive relation for the transseries coefficients as where u (m ) i,0 = 1 for m i = δ i,i , we can solve this equation order by order to get the transseries data u (n) 1,k , S i andb i without any (imaginary) ambiguity. Note that the IR data should match the numbers found using linearization around the asymptotic fixed point in Eq. (37). The technique advocated in the N = 0 case above can be extended to include higher order moments c nl . The energy dependence of hydrodynamization is only materialized by c nl (n = 1, . . . , N) even if c 0l have been for the most part sidelined in the literature. To get a complete and correct picture, however, all the moments ought to be taken into consideration. In principle, this means that one has to eventually solve an infinite-dimensional dynamical system for getting the full hydrodynamization process in both IR and UV regimes, which is obviously ambitious. As a result, a practical approach is to consider a truncated dynamical system whose proposed transseries solution has of course the ability to be generalized to the exact result. So our starting point in this section will be Eq. (22b), α nl c nl+1 + β nl c nl + γ nl c nl−1 − n (ρ l c n−1l+1 + ψ l c n−1l + φ l c n−1l−1 ) + 3c nl 2θ 0 .
The equation (39) can be written in a concise matrix form, where the quantities c,Λ, A are explicitly given by and we have defined the block matrices B and D as Here, the blocksB 00 ,D 00 are both L × L matrices,B 10 ,D 10 are (L + 1) × L matrices, and the remaining blocks are (L + 1) × (L + 1) matrices. As in the N = 0 case, we will be defining the matrix U to diagonalize B, We finally employ the transseries ansatz in the diagonalized form of Eqs. (40)-(41) as before to obtain the recursive relations involving the transseries data, namely , and m is the index of each non-perturbative sector of the pseudomodes, c i (w). Eq. (51) can be solved recursively. The higher non-perturbative sectors are connected to the lower ones through modeto-mode couplings of the form c 01 c nl , whereas operations in charge of promoting the asymptotic order are differentiation and multiplication by a 1/w factor. Once the perturbative order is fixed, the 1st-order non-perturbative sector can be constructed explicitly.
For pedagogic purposes, we show how to compute the leading order contribution in the perturbative sector of the transseries corresponding to m = 0 (also known as the leading order in the asymptotic expansion). As an example, first, we determine which moments behave like O(1/w). This is effectively done by taking k = 0 in (51). Since the asymptotic fixed point has all the moments other than c 00 go to 0 at w → ∞, the coefficientsũ 0 i,0 have to vanish. By default, all the other coefficients promoting the asymptotic order also vanish. For the lowest asymptotic behavior for each n and l, Eq. (40) gives Here, we relabeled the index i by n, l for convenience. Since c nl ∼ O(1/w) and A nl = 0 only when (n, l) ∈ {(0, 1), (1, 1)}, we find that both c 01 and c 11 have a non-zero coefficient for the 1st-order asymptotics given by Since A is constructed by γ 01 and φ 1 , only these two moments survive at the 1st asymptotic order, that is, c 01 (cf. Eq. (34) aka Navier-Stokes) and c 11 ∼ −φ 1 θ 0 /w. For the other moments, the leading asymptotic order is dictated by their neighbors c n l alongside mode-to-mode coupling. The equation responsible for identifying the leading order of c nl ∼ O(1/w k ) (k > 1) reads All the other operations promoting the asymptotic order vanish, such as differentiation. For example, the leading asymptotic order of c 2l (l = 0, 1, 2) is identified by As can be seen, they both have the lowest asymptotic order 1/w 2 because of c 11 . By virtue of Eq. (51), the leading order asymptotics of higher moments is determined to be c nl ∼ O(1/w k ) where k = max(n, l) for n, l ≥ 2. The higher modes, therefore, decay faster in general as all the Lyapunov exponents are equal due to the RTA-BE having a single scale which in this case is θ 0 , and the suppression in the IR is just determined by the leading asymptotic order. All the modes satisfying n, l ≤ 1, however, are exempt from this rule as c 10 ∼ 1/w 2 , and as mentioned earlier, the two remaining lower modes c 01 and c 11 are the slowest modes of all in the IR.
Once the coefficients of the bare asymptotic series for each c nl are obtained, the transseries coefficients can be achieved by solving the equation including higher m > 0 non-perturbative sectors. As an example, we show the transseries and leading order bare asymptotics of five moments c 01 , . . . , c 21 of the truncated dynamical system N = 2, L = 1 in Fig. 1 4 . Each moment is renormalized by the inclusion of non-perturbative/perturbative sectors available using the transasymptotics as C nl,k /w k + O(1/w k ), where the transasymptotic matching is approximated up to 15th-order transmonomials. On the numerical front, the integration constants σ are determined by employing simultaneous least-squares method, which aims to minimize the difference of the exact and transseries solutions of all the moments involved in a given truncated dynamical system simultaneously. So the overall average error will be reduced across all the moments. In this figure, five integration constants σ 1 , . . . , σ 5 are given by their optimized values along with individual standard deviations reported in the plots. Also, the light blue dashed lines stand for the transseries solution with the best optimized σ i while the blue shaded areas highlighting the possible variation of transseries due to standard deviation. Furthermore, bare asymptotic expansion at 1/w, 1/w 2 , and 1/w 3 orders for each moment is plotted as a comparison. These plots show that continuing to a higher asymptotic order will not necessarily lead to a better result, as opposed to adding more transmonomials at each order, which results in a better approximation of the transasymptotic matching overall (see next section).

C. Transasymptotic matching
In this section, we construct the transasymptotic matching condition responsible for the full-blown form of the time evolution of the transport coefficients including all the effects of non-perturbative sectors, which is a 1st-order PDE. In doing so, we first 4 In this work the initial condition for the distribution function is given by the RS ansatz of the distribution function [74], i.e., where Λ 0 and −1 < ξ 0 < ∞ are the initial effective temperature and momentum anisotropy along the ς direction. Thus the set of initial conditions for the moments c nl depend on the initial time τ 0 and initial anisotropy parameter ξ 0 . See App. E for further details. T 0~5 00 Mev, ξ 0 =1000 and sum over m in Eq. (51). Then the transasymptotic matching condition turns out to have the following formal form Here,C i,k is supposed to be only a function of σ i ζ i . We again note that Eq. (60) is a 1st-order PDE, whose solution yields the time evolution of the transport coefficients with all the σ i ζ i (w) included. Moreover, the integration constants may be determined from the transseries coefficientsũ (m) i,k . Due to its partial differential nature, it is technically not easy to solve (60) for the general case even though for L = 1, N = 0, it can be exactly solved. For instance, the first three coefficients are given by solving (60) in this case as [41] where For a system consisting of a larger number of moments, the above equation becomes a PDE, which renders the matching condition hard to solve. However, the transasymptotic matching, in this case, is at least effectively consolidated by summing over a large number of transmonomials just in the same manner as in the renormalization group equation for a running coupling constant, where loop diagram contributions are added to the r.h.s. in order to capture a more realistic RG flow. As a special example, we guide the readers to take a look at App. D 1 in which the transasymptotic matching for N = 0, L = 1, ∆ = 1 system yields polynomials of finite degree, which are the exact solutions to the dynamical system. In Fig. 1, the transmonomials are taken up to 15th order to get as close as possible to the exact solution of the truncated system. We again mention that the renormalized moment is meant to be a solution compatible with the full-fledged transasymptotics at its leading asymptotic order.

D. Comment on resurgence and cancellation of the imaginary ambiguities
Once one has an ODE of the form (11), we can realize imaginary ambiguity cancellation in a systematic way. In general cases, it is hard to prove the cancellation of imaginary ambiguities. However, for ODEs of the type (11) there are mathematically rigorous proofs for the ambiguity cancellations due to O. Costin (cf. [26] and references within for more technical details). In this section, we consider an easy example involving only one mode i = 1, and therefore we omit the index for simplicity.
Let us begin with a formal transseries ansatz given by Suppose that F(w) solves (11) and let us only focus on the 0th non-perturbative sector a (0) k . By calculating the radius of convergence one can see that the power series (64) is divergent. The approximate form of the growing upper bound is with some positive real constants, M and S. This means thatΦ (0) (w) is asymptotically of Gevrey-1 class 5 , hence the Borel transform is just enough for convergence purposes, defined by Note that Φ (0) (ξ ) has a finite convergence radius and r c = S. SinceΦ (0) (w) is a divergent series, there should exist a branchcut on the positive real axis with a branch-point S on the Borel plane. The position of the first singularity corresponds to the Lyapunov exponent of the flow lines approaching the asymptotically stable fixed point in the dynamical system. Notice that the appearance of a branch-cut is a reminder that the series expansions are coupled to higher order transmonomials such as exp(−Sw)/w which independently satisfy the same ODE, i.e. Eq. (11).
The original divergent series (64) can be reproduced through the Laplace integral, and subsequently taking the asymptotic limit It is noteworthy that the imaginary ambiguity depending on the contour of integration is nothing but an artifact of going to the Borel plane and it is exponentially suppressed in the large-w limit. However, when one observes a singularity on the Borel plane, one can build a relationship among non-perturbative sectors by taking the Hankel contour going around the singularity in the Laplace integral, The relations obtained in this way are called resurgence relations.
We now are in position to introduce the Stokes automorphism S θ . Formally speaking, the integration constants in the transseries may jump once the singularity (Stokes) rays on the Borel plane are crossed. This is associated with the existence of an automorphism that takes one integration constant to another plus some new contribution, defined by where {ρ θ } denotes the set of singular points along the Stokes ray with the contour enclosing the singularities being controlled by an angle parameter θ . Here, ∆ ρ represents an abstract derivative operator known as alien derivative which can be understood as an infinitesimal change in the asymptotic behavior due to a singularity in the direction parametrized by θ . Note that is not a divergent series along θ direction. However, when a singularity appears for a particular angle θ , we can make the resurgence relations among different sectors. In our case, the set of singularities on the real positive axis can be given by {ρ 0 } = {nS | n ∈ N}, where again n = 1 shows the Lyapunov exponent of the flow lines flowing to the asymptotic stable fixed point of the dynamical system. By taking the Hankel contour with a singularity at nS, one can obtain the so-called bridge equation∆ and the relationships among different sectors are given by where A ∈ i R the famous Stokes constant. This equation means that information of some sector is carried over to higher sectors once Stokes rays are crossed.
Although the construction of resurgence relations is generally an independent issue from the imaginary ambiguity cancellation, we can systematically find such a cancellation mechanism for the transseries solution of the ODE in (11) without the need to resort to the discussion of this section. The significantly important fact is that a nonzero r.h.s. Eq. (76), satisfies the same ODE. Therefore, the Borel-summed transseries along the θ = 0 ray may well be expressed by the arbitrariness in the integration constant, and the fact that nothing prevents one from extending σ ∈ R to σ ∈ C. The relationship between two integration constants across a Stokes ray is then given by Therefore, the imaginary ambiguity can be cancelled by shifting the integration constant in such a way that Let us demonstrate the imaginary ambiguity cancellation for simple cases. For L = 1 and N = 0, the leading large-order growth of the asymptotic expansion is given by [75] where S 1 = 3 2θ 0 ,b 1 = −18/35, and C(θ 0 ) is an overall factor depending on θ 0 . Here, we have usedũ (1) 1,0 = 1 in the second line. Since w always appears with a constant factor 1/θ 0 inc l (w), one finds that where C 0 is a constant. We measured C 0 by numerical fitting with the a prepared form (79) and obtained Therefore, the Stokes constant A 1 is related to C(θ ) in the following way leading to the imaginary ambiguity cancellation once we set For general L and N,ũ l,k are complex-valued so are the integration constants σ l for the reality condition to hold. For L = 2, the overall factor in (79) is found to be in whichb 1 = (75 + i √ 10655)/110 =b * 2 . By fitting numerically we can estimate the complex constant C 0 as Therefore, the Stokes constant as a function of θ 0 can be fully evaluated from Eq. (82). Finally, the imaginary ambiguity cancellation follows by just shifting the integration constant as

IV. GLOBAL DYNAMICS AND ASYMPTOTIC STABILITY ANALYSIS
In this section we will explore the generic properties of the Bjorken flow from the perspective of dynamical systems with some useful definitions provided in App. B for readability.

A. Phase space and asymptotic UV/IR fixed points
The explicit time-dependence of (9) means that the system is non-autonomous and the phase space is larger than that of timeindependent (autonomous) system with the same number of variables. Let us assume that the total number of Legendre and Laguerre terms in the expansion of the distribution function is L + 1 and N + 1, respectively. Then one can write down the phase space of this truncated dynamical system as where each subspace M i consists of L dimensions, and we excluded the constant component c 00 = 1 from the phase space. So (87) is just the hyperspace c 00 = 1 of the full phase space. Eventually, the limit N, L → ∞ has to be taken, meaning that the true physical dynamical system at its full glory is infinite-dimensional 6 . The tuple (M 0 , . . . ) is just short for the product of all the entries, and the time manifold t is taken to be the set of positive real numbers R + . This choice is made due to the fact that Eq. (9a) is not regular at τ = 0 so to keep the continuity, τ < 0 is not allowed. Since temperature T has its own equation in (9), we also add its manifold T to the crowd and let it take values over T = R + . In total we are then left with (L + 1)(N + 1) + 1 + 1 − 1 = (L + 1)(N + 1) + 1 dimensions for M . By looking at the case where N = 0, the dynamical system reduces to L + 2 dimensional subspace (M 0 , T, t) which is denoted by M 0 ⊂ M . The special facet of this subspace is that it indeed encompasses the invariant manifold 7 because its equations are completely decoupled from the rest (n > 0) in case we started with a larger system involving more moments, i.e., N > 0. Therefore, any solution of this sub-system entirely lies inside M 0 at all times τ > 0. The second observation is that the mode-tomode coupling terms in (9) for n > 0 all depend on the elements of the subspace M 0 . So it is crucial to understand the dynamics happening in this subspace by first locating the asymptotic fixed points and analyzing their stability and finally obtaining the flow lines connecting these fixed points to find the qualitative shape of the invariant manifold.
In terms of the variable w = τT (τ), the non-autonomous dynamical system is described by where w = R + := R + ∪ {0}. Note that w = 0 has been allowed since the combination of τT (τ) does essentially have a welldefined zero limit. The invariant manifold now is inside W 0 parametrized by c 0l . The structure of this invariant manifold will be the subject of next subsection. There are however three major differences between the two parametrizations of the phase space of the Bjorken flow. Namely, 1-In terms of τ, the early-time τ < 1 or UV limit of the theory does have a pair of fixed points for the truncation order L, N = ∞; one being the maximally oblate point at which the longitudinal pressure p L vanishes and transverse pressure p T is maximized, and one more fixed point, at which p L is maximum and transverse pressure is vanishing, to be called maximally prolate point. 3-In the original version of the dynamical system, Eqs. (9), the temperature does not converge in the UV so practically speaking it is impossible to connect continuously to the maximally oblate fixed point at T = 0 by running the evolution equations backward in time, say, from somewhere in the IR. It is also not possible to find a complete flow connecting to another UV fixed point present in the Bjorken flow phase space M [41]. In terms of w, the UV structure of the theory is altered through a change in the stability of fixed points at w → 0 even though the IR will remain intact as mentioned above. It is now possible to search for flow lines starting at the maximally prolate point (see Fig 3 for instance), but still it will not be feasible to have a flow that begins its journey at exactly the maximally oblate fixed point. This latter solution, if existed, would be a critical line due to the extreme fine-tuning required around a saddle point to initiate a structurally unstable flow on the boundary of the invariant manifold that connects to the IR fixed point c nl = δ n0 δ l0 at w → ∞. This is explained below. If we consider the full theory with L = ∞ and N = 0 in the M phase space of the Bjorken flow, there are two UV fixed points in general whose coordinates are given by maximally oblate point : where we have adopted the notation (τ, T, {c 0l , c 02 , . . . }). In Fig. 2 the stability analysis of these two fixed points in the truncated system N = 0, L = 31 as well as N = 30, L = 31 is shown. At a general truncation order, the stability analysis is summarized as follows for any relaxation time τ r = θ 0 T 1−∆ (0 ≤ ∆ ≤ 1): (a) N = 0, L = odd: The system has 2 fixed points in the UV. The maximally prolate fixed point is a source while the maximally oblate fixed point is a saddle point with one repelling direction, i.e., index 1 fixed point); (b) N = 0, L = even: The system has only an index 1 maximally oblate point. The total index is still conserved since the other fixed point becomes a complex saddle; (c) N > 0, L = odd: The maximally prolate fixed point becomes an index(I − N + 1) saddle point whereas the maximally oblate fixed point is an index(N + 1) saddle point; (d) N > 0, L = even: The maximally oblate fixed point is a saddle of index(N + 1).
Here, we have again used the notation I = (L + 1)(N + 1) − 1. In L = odd case the total UV index is always I + 2.
Considering that an asymptotically stable hydrodynamic fixed point has zero index, then the total Morse index is I + 2. A comparison between the numerical solutions of the flow lines in L = odd system and the exact RTA-BE shows that I + 2 is the correct topological invariant of the real phase portrait of the RTA-BE. Note that the temperature is always an attractive direction in both UV and IR, and therefore it does not contribute to the index. Furthermore, due to the lack of one real fixed point (maximally prolate point) in the UV when L = even, this kind of truncation will not be physically sensible and will be hereby omitted from our discussions. 8 The phase portraits of different 3d sectors of the dynamical system in (39) at odd truncation orders are plotted in Figs. 3 and 4.
An interesting question that comes to mind concerns the existence of a flow line that connects the maximally oblate fixed point to the asymptotic hydrodynamic fixed point in the phase portrait of the Bjorken dynamical system using w coordinate. There are two things that one might keep in mind here. First, the invariant manifold is in the subspace W 0 ⊂ W or simply the subspace of the phase space with coordinates c 0l . Since the maximally oblate point for N > 0 is a saddle point of index N + 1 (see (c) above), it must be located on the boundary of the invariant manifold. So at least in the N > 0, L = odd system, there cannot be a critical line since it will always flow outside of the basin of attraction of the invariant manifold due to the repelling directions being in the subspace N n=0 W n . Second, for N = 0, the maximally oblate point is essentially a sink (time is the only repelling direction). As a result, the c 0l set initially to be the coordinates of this fixed point at w = 0 will change its position with time, while being a fixed point at any moment w > 0 until eventually it becomes the hydrodynamic fixed point at w → ∞. So there will not be a flow line (process) to begin with, proving that the critical line does not exist in the w parametrization of the Bjorken flow just like it could not exist in the τ coordinate but for a completely different reason. 9 FIG. 2: Eigenvalues of the Jacobian matrix around asymptotic UV fixed points in the truncated system. Red and blue points denote eigenvalues of the maximally oblate point and the spiral source, respectively. In (a), maximally oblate point behaves like a spiral sink because ℜ(b i ) are all negative. For the maximally prolate point, however, ℜ(b i ) are all positive, a property that classifies it as a spiral source. In (b), both of these points are shown to be able to have positive and negative ℜ(b i ) at the same time, hence they are in general spiral saddle points. From the perspective of transseries, upon choosing those integration constants σ i = 0 associated with the monomials with b i satisfying ℜ(b i ) > 0, we can make the stability of the UV fixed points change to a source in the general case N > 0.
FIG. 3: The phase portraits for the system truncated at N = 0, L = 3. The initial value of the flow is taken to be c nl (w 0 ) = c nl (w 0 ; ξ 0 ). In the UV (early time), the flows are coming out of a spiral source as seen on the left where anistropy parameters are close to −1 that corresponds to c 01 ∼ 6.6 [41]. In the range closer to the maximally oblate point, that is 10 < ξ 0 < 1000, the flow lines do not converge since they are initiated outside of the basin of attraction (of the invariant manifold). The maximally oblate point is exactly located on the boundary of the basin. Note that the flow lines are shown to be spiraling in the UV as w → 0, a symptom of being in the vicinity of a spiral saddle. We can see the forward attractor in the IR for the flows initiated in the basin of attraction.

B. Initial value problem for transseries and invariant manifold
Although the transseries contains a large amount of information about a given dynamical system, fixing the integration constants is an important, yet challenging, problem. On the one hand, since our (truncated) dynamical system has I variables-thus I first-order ODEs, it consequently has I integration constants needed to be tuned in order for the basin of attraction (of the invariant space) to be determined. On the other hand, the initial condition of the exact integral solution of the Boltzmann equation has only three parameters, namely τ 0 , T 0 , and ξ 0 .
By redefining the ODEs using w as a time coordinate, T 0 can be completely separated from the rest of the dynamical system which now has the form (22a). The new initial time w 0 determines a point on a flow (more precisely a process), hence σ i do only depend on the anisotropy parameter ξ 0 . We could behold from the numerical results of the integral solution of the Boltzmann equation that every flow is closed on an open subset of the W as long as ξ 0 is chosen from the domain [−1, +∞). Moreover, we recall that if a flow is initiated in the subspace W 0 ⊂ W parametrized by c 0l , it always stays in that subspace, meaning that the invariant manifold has to be an embedding in the space (W 0 , w). But how does transseries know about the invariant manifold? Since ξ 0 determines the basin of attraction of the isolated invariant space of the Boltzmann equation, there has to be a map from the ξ 0 -space to W 0 such that where w = [0, +∞) is the time manifold and F is the space of integral curves given by the n = 0 ODEs in the dynamical system. Because of the uniqueness theorem, the map in (91) is a bijection only when L, N = ∞, and consequently one can find that the invariant manifold is a two dimensional surface embedded in the space F . It is a quite challenging problem to figure out the correspondence defined by (92), but we can narrow R I down to a subspace formed by those σ(ξ 0 ) which yield flow lines starting in the basin of attraction of the invariant manifold. This means that σ(ξ 0 ) is able to give the invariant manifold once the stability analysis of each fixed point in the UV is performed and the approximate boundary of this space is determined by finding possible critical lines. To facilitate the search for the invariant manifold, one can construct the transseries around every fixed point individually at w → 0, which would be something for the form where the anomalous dimensionsb i of the pseudomodes are now the eigenvalues of the linearization matrix about the corresponding fixed point. 10 Choosing σ i = 0 for every i for which ℜ(b i ) < 0, we get a maximally prolate point that is a source (that is otherwise a spiral) fixed point. We will discuss the stability of the UV fixed points in Sect. IV A. An immediate question that comes to mind is whether one can construct a complete flow line and/or a critical line from the transseries solutions around individual fixed points. One part of the answer involves topological arguments to be discussed in an upcoming work from the perspective of dynamical systems and equivariant cohomology Conley index [76]. Another piece comes from transferring the information from one transseries to another by tuning σ i to a particular set of numbers from both UV and IR sides, provided that the flow line achieved this way remains on the invariant manifold at all times.
As a final remark, let us mention that one can analytically continue the transseries by rotating in the time axis w → e iα w such that transseries solution around the asymptotic hydrodynamic fixed point is now exactly a Fourier decomposition at α = π/2. Neglecting an overall imaginary factor at each non-perturbative order, it is easy to see that σ i serve as amplitudes of the individual wave components in the Fourier decomposition. The significant advantage of analytic continuation is that the attractor entails oscillatory components which were exponentially suppressed in real time, thus fitting the transseries to the complex numerical solutions of the analytically continued Boltzmann equation is much more tractable.

C. Dynamical renormalization of 2nd order transport coefficients
In this section, we aim to construct an RG equation for the transport coefficients in analogy with the gradient descent approach to the RG flows in the context of quantum field theories using the language of dynamical systems (cf. Ref. [77]).
The variable w so far has been playing the role of flow 'time' in our system of ODEs, i.e., Eq. (22b). But in this section, we want to interpret it as playing the role of 'energy scale' for the renormalization scheme arising from resuming all the nonperturbative fluctuations around the asymptotic expansion of non-hydrodynamic modes c := {c nl }. The reason is simple: the Knudsen number Kn = τ r (τ)|D µ u µ | ∼ w −1 in Bjorken flow provides the only parameter by which the system evolves from the UV regime all the way to IR. Since the transseries coefficients of individual moments include transport coefficients with exponentially suppressed factors in w following [41], the far-from-equilibrium dynamics of the moments can then be associated to the running of these coefficients as w changes. Hence, preparing a time-dependent (non-autonomous) dynamical system for c automatically amounts to having a renormalization group equation on the phase space of moments c and w, together denoted by W , the so-called phase space of the Bjorken flow in w parametrization as in (88).
A flow (process) on the phase space is described by a continuous map such that φ w,w 0 (c 0 ) = c(w), and φ w,w 0 (φ s,w 0 (c)) = φ w+s,w 0 (c) where w ∈ R + is the RG time. The IR regime of the theory is captured by the behavior of the flow lines approaching to an asymptotic stable fixed point aka hydrodynamic fixed point satisfying c eq l = δ 0l ∈ M for which φ w,w 0 (c eq l ) = c eq l for all w ∈ R + . We alternatively refer to c eq l as the asymptotic IR fixed point since it is reached at w → ∞.
Let us formulate our dynamical system as where µ := w/θ is the RG scale and log(µ) ∈ R. As µ varies, the moments c are mapped to themselves in a self-similar way, and any initial RG scale µ 0 can initiate a process to a later µ by the renormalization group action. Since the r.h.s. only depends on c and µ, one can write This is nothing but an RG equation in the space of moments c and the vector β(µ, c) consists of β -functions, each encoding the dependence of every non-hydrodynamic mode on the renormalization scale, µ, in the process of equilibration. In this expression, it is convenient to redefine our notations asũ l,k →ũ l,k θ k , F l,k → F l,k θ k , Note that the ordinary derivative along an RG flow can be expressed by two partial derivatives as In the formalism we have been seeking to build so far, one prominent assumption is that an observable O is able to be expressed in terms of c i (w) (orc i (w)), namely O = O(c(w)). Therefore by solving the RG equation via transseries, we can achieve a renormalized form of the observable, say a transport coefficient, up to an arbitrary order of our choosing.
It is straightforward to define the RG equation by starting with the ODE in (40). To make things more tractable, we again go back to the RG time w and compute the derivative of an observable O(c(w)) w.r.t. log w dO(c(w)) d log w = wheref = Uf andβ := wf. Here, use was made of d/d log w = (b − Sw) ·ζ + ∂ /∂ log w and the fact thatc i (w) = ∑ ∞ k=1C i,k (σζ(w))w −k . We should also point out thatζ i := ∂ /∂ log(σ i ζ i ). The first line describes the scaling behavior of O in terms of w, and the second line contains information of the dynamical system through β. Hence, solving the RG equation determines a renormalization of the observable O. Since we want to consider the RG equation for the transport coefficients, we are slightly better off with the definition of O changed to O = O(C k (w)). Roughly speaking, this implies that O is a function of the transport coefficients and depends only on ζ(w). Consequently, the scaling behavior could be derived in a fashion similar to what we did to obtain (99), so However, the link between the scaling behavior and the dynamical information hidden in β is very non-trivial. Yet, one can proceed to calculate the RG equation for the L = 1, N = 0 system. To do so, we employ the scaling behavior in (100) which gives Suppose now that ζ 1 and w are independent of each other, i.e.,c 1 =c 1 (σ 1 ζ 1 , w) andC 1,k =C 1,k (σ 1 ζ 1 ). In this case we have To remove w from the l.h.s. of Eq. (102), we take a contour integration around the origin after multiplying by w k−1 , which then givesζ Hence, Eq. (100) for L = 1 and N = 0 reads whereζ 1C1,k (σ 1 ζ 1 ) is given by Eq. (103). The definition in (104) has a direct connection with the dynamical information encoded in β 1 , so that (104) is regarded as the RG equation of the transport coefficients by setting O(C 1,k (w)) = C 1,k (w). In the L = 1, N = 0 system, choosing k = 1 in Eq. (104) and solving for the renormalized quantity C 1,1 gives that it is proportional to the shear-to-entropy ratio at the equilibrium [41]. Plugging C 1,1 back in the expansion of the non-hydrodynamic moment c 1 obtained using the linear response theory (see App. A) automatically promotes the 1st-order transport coefficient to the non-equilibrium case compatible with the transasymptotic matching along an RG flow approaching the IR fixed point asymptotically. This essentially gives a dynamically renormalized η/s [41] as follows where − 3 40 C 1,1 (w → ∞) ≡ (η/s) 0 [41]. Likewise, the 2nd order renormalized transport coefficients are given by where again 9 40 C 1,2 (w → ∞) ≡ T τ π η s 0 − λ 1 s 0 . We conclude this section by discussing whether there is some sort of a Lyapunov function that would identify the RG flows from a global dynamical standpoint. It is well-known that one can formally define an object analogous to a dynamical potential from the Newtonian mechanics. For simplicity, we assume that β is independent of w (a scenario that would hold once the limit θ 0 → ∞ is taken). 11 Now, we define a positive-definite differentiable function V dc i (w) We can easily show that V is a monotonically decreasing function in terms of w as thus a candidate for a global Lyapunov function satisfying the conditions in Eq. (B9). 12

V. HYDRODYNAMIZATION OF SOFT AND HARD MODES
In previous sections, we demonstrated that the moments c nl solving the ODEs in the dynamical system of (22b) turn out to be of the multi-parameter transseries form. In this section, we continue our quest for what type of information one can get from these solutions on a more physical ground by analyzing the late-time behavior of different momentum and energy sectors of the distribution function. In doing so, we stumble upon something interesting: a flow line in the space of moments initiated at an arbitrary state in the basin of attraction of the asymptotically stable fixed point, is controlled at late times by not only the known non-hydrodynamic moment c 01 but also the mode c 11 . The latter also happens to possess the same perturbative decay as c 01 . We recall that this immediately leads to the statement that the distribution function and any observable projecting onto/or involving at least the l = 1 sector of the distribution function would receive two major IR contributions. In this section, we shall analyze in great details the impact of the slowest (= O(1/(τT ))) non-hydrodynamic modes in the IR.
Let us start rewrite our ansatz as where the soft (s), semi-hard (sh), and hard (h) sectors of the distribution function are respectively defined by (113c) 11 When we keep w in the β -function, we have to promote the non-autonomous system to an autonomous system of one dimension higher by introducing an ODE for w in terms of a new flow time ρ as follows dc(ρ) d log ρ = β(c(ρ), w(ρ)), dw(ρ) d log ρ = β w (w(ρ)).
In order to assess which sector of the distribution function hydrodynamizes faster, we focus on the decay of the following normalized moments [53,78] 13M M nm are then certainly sensitive to the energy and momentum tails of the distribution function. To this explicitly, first note that M nm ≥ 0 if and only if f (τ, p T , p ς ) ≥ 0. This condition holds for the exact solution of the RTA-BE, but not necessarily for the approximate distribution function. Second, the Landau matching condition for the energy density impliesM 20 ≡ 1. Last, some of the moments in (114) such asM 10 = n/n eq. (, where n and n eq. are the non-equilibrium and thermodynamic particle densities, respectively) andM 01 = P L /P 0 have a direct interpretation in terms of the usual macroscopic variables. After taking the timeM nm in Eq. (114), we end up with the following equation The solutions to this system indicate that different momentum tails of the distribution function reach values at equilibrium state asymptotically, that is the asymptotic hydrodynamic fixed point. In order to put the problem into the framework of dynamical systems, we writeM nm (114) as a linear combination of the moments c nl using the ansatz (3) of the distribution function, i.e. 14 M nm (τ) = 1 + 2m + 1 2 ∑ k,l=0 k+l>0 where Γ(x) is the Gamma function and a b is the binomial coefficient. Therefore the solutions to the dynamical system in (115) are also written in the form of multi-parameter transseries. 15 These solutions flow to the stable fixed point at late times once the initial state is chosen in the basin of attraction of the invariant manifold (of the truncated system), which indicates that each momentum sector is guaranteed to equilibrate eventually. In other words,M nm → 1 since c eq. nl → δ n0 δ l0 asymptotically. In Sect. III B, (cf. Eq. (53)), we showed that the slowest non-hydrodynamic moments decay perturbatively like c 01 , c 11 ∼ 1/w and c 10 , c 20 , c 21 , c 22 ∼ 1/w 2 at large w. Thus, the asymptotic behavior ofM nm is dominated by the slowest moments appearing in each sector of the distribution function.
We study first the normalized momentsM nm for n, l > 0 where there are contributions coming only from both soft and semihard sectors the distribution function whose slowest non-hydrodynamic modes are c 01 and c 11 , respectively. In this case, the impact of these non-hydrodynamic modes is determined by taking the following asymptotic limits in Eq. (117): • Soft (s) regime in which there is only the leading order asymptotic contribution of c 01 ∼ 1/w, i.e., 16 13 In this work we follow the notation of Strickland [78]. 14 In order to get this expression, we utilize the following identities [79,80] ∞ 0 dx e −x x γ−1 L (µ) 15 Alternatively one can also find a compact form of the momentM nm from the exact RTA-BE solution of the distribution function (see Eq. (A7) in App. F). 16 In Eq. (4.6) of Ref. [78], the author assumes explicitly the 14-moment approximation when truncating the distribution function, that is By substituting this expression in our definition (114), one obtains the result derived by Strickland, Eq.(4.8) [78], which obviously differs from ours (118). This discrepancy arises from the truncation procedure. Here, we do not truncate the momentum and energy dependence of the distribution function and rather keep these in their exact form. Instead, we truncate the number of moments entering the distribution function at our discretion. We remind the reader of the work of Denicol et al. [47], where it was shown that the 14-moment approximation does not provide a unique set of equations of motion for the dissipative currents, an ambiguity which was resolved in Ref. [46]. • Soft + semi-hard (s + sh) regime which incorporates the leading order asymptotic contribution of both c 01 and c 11 up to O(w −1 ), namelȳ In Fig. 5, a number ofM nm are plotted versus the variable w, obtained from the exact RTA-BE solution (E8) (black line), the s (118) (green dot-dashed line), and the s + sh (119) (blue dashed line) regimes. We verify numerically that the IR behavior is not affected by the choice of initial conditions as long as they agree with the bounds set by the invariant manifold. In order to get a hold of how fast the normalized moments equilibrate, we set an arbitrary saturation bound at which a given moment falls within 5% of its equilibrium value, that is to say M nm = 1 ± δ with δ = 0.05, and the + or − is taken depending on whether the moment monotonically increases or decreases, respectively. In this setting, as seen in Fig. 5, it can be verified that independently of the truncation scheme, the normalized moments will asymptote to one, as expected. In the UV, however, the s and s + sh regimes disagree with the exact RTA-BE result. This should not come as a surprise since both limits are only valid when the distance from the hydrodynamic fixed point is small. But in the IR region, the momentsM nm s first reach the saturation bound whilē M nm s+sh andM nm exact approach this value much later on, and above all else, almost at the same time. In the same limit, we behold a remarkable agreement between the exact and s + sh solutions. Thus, the asymptotic behavior of the normalized momentsM nm (n, l > 0) can only be given entirely by a linear combination of both non-hydrodynamic modes c 01 and c 11 . In Fig. 6, the value of w = w c is numerically evaluated. As is shown, if the indices for energy, n, and longitudinal momentum, m, are increased, w s for the exact and s + sh moments marks a later time. The convergence ofM nm s depends only on m, which is evident from Eq. (118). The numerical results in this figure attests that in fact the exact and s + sh moments are in good agreement, and therefore s + sh is the correct description of the IR regime of different energy and momentum sectors of the distribution function solely based on a linear combination of both c 01 and c 11 as written in Eq. (118) 17 . Given this agreement, one concludes that the larger the index n > 0, the more relevant the inclusion of the new non-hydrodynamic mode, c 11 . We want to mention that this is solely due to the power of dynamical systems since the 1/w decay of both c 01 and c 11 in the perturbation theory around the hydrodynamic fixed point is directly confirmed from the form of ODEs in (115).
We are now in position to analyze the normalized momentsM nm for n > 0, l = 0. In this case, only the hard sector of the distribution function contributes and the slowest modes playing a role in Eq. (114) are c 10 and c 20 , both asymptotically decaying as 1/w 2 . As a double check, consider for instance the modeM 10 = n/n eq. . For the Bjorken flow, the Chapman-Enskog expansion in the RTA approximation [81] shows that the NS corrections to the particle density n vanishes as chemical potential µ is turned off, and thus the leading order dissipative corrections are of the second order in the Knudsen number ∼ 1/w 2 . In general, it is known that second-order fluid dynamics theories, e.g. Israel-Stewart theory, do not reproduce properly the behavior of the heat flow and particle density as well as the dominant corrections arising from couplings between the particle density and the shear stress tensor [81,82]. Now, as done before, we study the impact of the non-hydrodynamic modes entering Eq. (114) by taking two truncation limits 17 A careful reader would notice a minor disagreement (less than 5%) between the values of w exact c and w s+sh c . This small mismatch is due to the difficulty in extracting numerically the large-time behavior of the exact RTA-BE solver used here. While concluding the draft of this paper, we were informed of an optimized algorithm that uses logarithmically-spaced grid in proper time and reduces the numerical error at late times. This new code is available in the ancillary files included in the arXiv version of Ref. [78]. In Fig. 7 the evolution of several normalized momentsM n0 are shown. In the IRM n0 c 10 +c 20 (h) andM n0 exact reach the asymptotic hydrodynamic fixed point faster than the momentsM n0 c 10 (h) , supporting the compatibility of the exact normalized moments with the former. Note that the asymptotic behavior of the normalized moments M n0 is entirely determined by the linear combination of the non-hydrodynamic modes c 10 and c 20 . This conclusion does not hold when n ≥ 3, where only the mode c 10 contributes.
In Fig. 8 the values of w c vs. n for the moments M n0 are plotted. It is seen that the saturation bound is approached as n increases if n ≥ 3, as opposed to when n = 0, 1 for which w c decreases. Now, when comparing with the results plotted in Fig. 6, we observe thatM nm (n > 0, m ≥ 1) asymptote to their values at the equilibrium state later than the momentsM n0 do. In terms of our orthogonal polynomial basis expansion, it is easy to understand this as the slowest non-hydrodynamic modes in the soft and semi-hard sectors decay like ∼ 1/w, whereas in the hard sector, the slowest modes die away like 1/w 2 . On the other hand, for the given bound δ = 0.05, there is a disagreement between all the hard asymptotic limits (120) and the exact results due to the truncation limitations. As we emphasized previously, close to the IR fixed point, the inclusion of both moodes c 10 and c 20 (Eq. (120b)) is a must, which is also seen in Fig. 8. The rest of this section is dedicated to showing how the presence of both c 01 and c 11 is a crucial factor in unveiling the asymptotic behavior of the exact distribution function in the IR. We calculate the deviations from the thermal equilibrium distribution function defined as following the procedure outlined in Ref. [83]. In the s and s + sh regimes, δ f reads where T NS denotes the NS solution of the conservation law (9a). The Fig. 9 outlines the time evolution of δ f for the exact and asymptotic limits. As a double check, other valid initial conditions are considered by varying the anisotropy parameter ξ 0 . Two remarks are due here. First, no momentum and energy sector of the distribution function does in fact thermalize homogeneously, and second, the low energy particles (p τ /T < 1) thermalize faster than the highly energetic (p τ /T > 1) particles. This is not surprising since it is known that for weakly coupled systems, and in the absence of external fields, soft particles equilibrate faster [61,62,84,85]. On the one hand, the IR limit of s truncation (122a) does not capture the relaxation of the energy and momentum tails when compared against the exact RTA-BE solution. On the other hand, the exact RTA-BE result is perfectly compatible with the s + sh limit (122b), which in turn confirms the necessity of including the new mode c 11 . Therefore, the proper description of the IR in the Bjorken flow has to incorporate the dynamics of this non-hydrodynamic mode as well.
In sum, the numerical results presented in this section together with the previous analysis of the multi-parameter transseries solutions for the moments c nl lead to the true mechanism behind the nonlinear relaxation processes of different momentum sectors of the distribution function, which is completely governed by various mode-to-mode couplings among the moments in the framework of dynamical systems. Furthermore, the IR regime of the distribution function is determined uniquely by two non-hydrodynamic modes c 01 and c 11 . The conclusions derived in this section would still hold if one used other relaxation time models (cf. App. F).

VI. CONCLUSIONS
In this work, we have proposed a new dynamical renormalization scheme which allows us to study the nonlinear transient relaxation processes of a weakly-coupled plasma undergoing Bjorken expansion. The distribution function was expanded in terms of orthogonal polynomials. The nonlinear dynamics of the RTA-BE is then investigated through the kinetic equations of the coefficients entering this expansion, i.e. the average momentum moments c nl of the distribution function. We not only developed further our earlier findings [41] to include higher modes, but we also found new interesting results summarized in the following • Based on the seminal works of Costin [26,86] we show that the coupled system of nonlinear ODEs for the moments c nl admit analytic multi-parameter transseries solutions. At every given order in the time-dependent perturbative asymptotic expansion of each mode, the summation over all the transient non-perturbative sectors appearing in the trans-series leads to renormalized transport coefficients. This presents a new description of the transport coefficients in the regimes far from equilibrium with an associated renormalization group equation, going beyond the usual linear response theory.
• As long as 0 ≤ ∆ < 1, the T = 0 limit of the Boltzmann equation is explicitly time-independent, meaning that the distribution function f does not evolve with time. In the language of our dynamical system, this limit is equivalent to saying that the system is autonomous. The stability of the maximally oblate point shows in this case that it is a sink and consequently, there is no flow line that could connect it to the asymptotic hydrodynamic fixed point at τ → ∞. As soon as T > 0, we have that the dynamical system explicitly depends on τ > 0 (i.e., it is non-autonomous) and continuously connected to the IR theory (hydrodynamics) at τ → ∞. But now the maximally oblate point is out of reach, hence a flow line (process) cannot be found with τ 0 = 0, T 0 = 0 which leads to hydrodynamical behavior at late times. In other words, the actual phase space M of the Bjorken flow is a disjoint union of hypersurfaces T = 0 and T = 0.
• The above conclusion does not hold true in the w = τT parametrization. The phase space W is different in the UV in the sense that the system stays always non-autonomous for all w ≥ 0. Also, notice that the temperature has been washed away from W at the cost of introducing a new singularity hypersurface c 01 = 20. However, the maximally prolate point is the only fixed point in the UV that can connect to the IR as shown on the l.h.s. of Fig. 3. The maximally oblate point that in principle should correspond to the free-streaming limit, is a saddle point of index N + 1. On the one hand, if N = 0, it is basically a sink from the perspective of an observer sitting in the subspace W 0 ⊂ W parametrized by c 0l . Therefore, if the initial conditions are set to be the coordinates of this fixed point, the observer only sees the maximally oblate fixed point moving to the IR as time passes. On the other hand, there is no bifurcation or change of stability along the way. Rather, the same fixed point at w = 0 is just dislocated to become the fixed point at w → ∞. As a result, there is no complete flow line connecting this fixed point to the hydrodynamic fixed point.
Also, if N > 0, there appear to be repelling directions in the subspace N n=1 W n . However, the invariant manifold is in W 0 , and hence the search for a solution hits a snag again. If there was such a line, then we would have an example of a critical line. This critical flow line is often mistaken in the literature as a global attractor as if the system was autonomous. In non-autonomous systems such as the Bjorken flow, the attractor is a concept encoding information of the past or future of the system of equations governing its dynamics independently of each other, as shown in Fig. 3. The relevant attractor is then either a forward or pullback attractor. In sum, given our study of the dynamical system for the Bjorken flow, we conclude that there is no critical line connecting the maximally oblate point to the hydrodynamic fixed point, which is otherwise known as "attractor solution" in the hydrodynamics literature.
• For the Bjorken flow the asymptotic behavior of the kinetic equations unveils the existence of a new non-hydrodynamic mode whose decay rate is exactly the same as the usual NS shear viscous component, i.e., ∼ (τT ) −1 . The origin of the new mode's decay is due to the nonlinear mode-to-mode coupling c 11 c 01 which dominates close to the asymptotic hydrodynamic fixed point. This non-hydrodynamic mode is the slowest high energy mode and it determines quantitatively the late-time behavior of the transient high energy tails of the distribution function.
In the present paper, we left out many interesting questions which certainly require appropriate answers. For instance, within our approach it is important to understand the analytic behavior of the retarded correlators together with the subtle interplay between branch cuts and poles as a function of the coupling [87,88]. It will be also important to investigate the possible phenomenological consequences of the new non-hydrodynamic mode c 11 . For instance, one might wonder about the impact of this high energy mode in the modeling of a jet crossing an expanding QGP (quark-gluon plasma) using kinetic theory. Another possibility is to understand if the origin of azimuthal anisotropies at large momenta is connected to c 11 . It has been argued by different authors that the azimuthal anisotropies at intermediate p T are related to non-hydrodynamic transport [89][90][91]. On a more theoretical ground, it would be very interesting to come up with extending our analysis to the challenging case of nonlinear PDEs such as Israel-Stewart hydrodynamic equations for either 2 + 1 or 3 + 1 dimensions. Finally, our work can be generalized to studying nonlinear aspects of time-dependent systems of ODEs whose description is cast in the form of the properties of a dynamical system. We have left these exciting topics for future research projects.
with f eq. = d f eq. (x)/dx, θ = ∇ · u,p := p/T and tilde indicates that the quantity is dimensionless. The transport coefficients in a E * (A6e) will be given below. For the Bjorken flow the gradient expansion of the distribution function reduces to the expression [63,64] δ f = χ p C p −p 2 4 9T 4 τ 2 −p 20 9T 2 τ 2 −p 2 4 9T 2 τ 2 −p 3 4 9T 2 τ 2 −p 5 64 315T 2 τ 2 +p 6 16 315T 2 τ 2 + · · · P 0 (cos θ ) 18 Here, we have used the Cauchy-Stokes decomposition of the fluid velocity together with the (conformal) conservation laws to first order in gradients The use of conservation laws at first order in the gradients ensures that the value of δ f does not change the equilibrium energy density [101].
whereχ p = −C p f eq . One immediately recognizes that the Chapman-Enskog expansion of the distribution function is an asymptotic series in terms of 1/(τT ). Close to the thermal equilibrium 1/(τT ) ≈ τ −2/3 , which in the case of Bjorken flow, it becomes proportional to the Knudsen and inverse Reynolds numbers. For the Bjorken flow using the second order Chapman-Enskog expansion for δ f (A7) in Eq. (4) results in the asymptotic series expansion of the moments c 01 and c 02 as follows The transport coefficients appearing in Eq. (A8) are given by [92] (A9c) Appendix B: A self-contained dictionary for the dynamical systems In this section, we list a dictionary of words from the terminology of the mathematical field of dynamical systems used in Sect. IV. This is a self-contained and compact dictionary that could help the reader with some terminology of time-dependent dynamical systems. 1-Dynamical system. A set of differential equations in terms of a state vector whose components are real numbers determined by a set of points in some suitable state space. Any small variation in the state of the system leads in turn to a change in the numbers. In short, the evolution of a dynamical system is based on a fixed deterministic rule that describes what future states follow from the current state. When explicitly time-dependent, however, equations of the system also evolve in form, and thus the rule changes as time passes.
2-Non-autonomous dynamical system. A differential equation of the form is explicitly time-dependent through the function g. This defines a non-autonomous dynamical system if the state function x(τ) is dynamically evolving with a flow time τ. T (τ), and equivalent of g in (B1) is the r.h.s. of (9).
3-(Asymptotic) fixed point. The time-dependent solutions x(τ) of g(x(τ), τ) = 0 are the moving fixed points of the system (B1). For τ → ∞ these define a point(s) x * at which the state vector is in a steady equilibrium state for all times t > 0 if τ = ∞ + t, referred to as fixed points of the non-autonomous system in (B1). Note that here the time τ takes values over t = R + 19 . If we reparametrize flow time to ρ such that t = R, the fixed points can now be equally defined at ρ → −∞.
Hence τ(ρ) turns into a new independent variable in the dynamical system such that is autonomous. The fixed points are then defined as usual by finding those (x * , τ * ) for which g(x * , τ * ) = 0 = h(τ * ). 19 In math literature, a non-autonomous system is known to not admit any fixed point by promoting it to an autonomous system of one dimension higher and observing that the r.h.s. of this secondary system never vanishes. Nonetheless, our definition of a fixed point for a non-autonomous system is not for all times as in the autonomous case, rather over just enough period of time beyond τ = ∞ in order to distinguish it as an equilibrium state. Studying topology of the flow lines around such fixed points is thus more complicated and requires a lot of additional mathematical input [76].

5-Phase portrait.
It is a geometric representation of the flow lines for a dynamical system in its phase space given a set of initial conditions represented by (x i , τ i ) with i ∈ I where I is simply an index set.
6-Invariant manifold. As defined above, it is a subset of the phase portrait which contains all the flow lines initiated inside it (or at its boundary) which keep being restricted in there at all times.

9-
The Hartman-Grobman (HG) theorem. It shows that near a hyperbolic fixed point (x * , ∞), the nonlinear system (B1) has the same qualitative structure as the linear system Therefore, HG implies that two hyperbolic systems are locally topologically flow equivalent iff their unstable manifolds have equal dimensions.
10-Index of (x * , ∞). The number of components of λ ∞ whose real parts are positive is called the (Morse) index of that fixed point, which is technically the dimension of the unstable manifold around every hyperbolic fixed point. So any nonzero index means that fixed point is basically unstable. This instability is severe for a larger index. For example, in a d dimensional system (x 1 , . . . , x d−1 , τ), an index n fixed point is an unstable saddle whose stable manifold is d − ndimensional. The instability of an index n fixed point known as a source is far more pronounced than a saddle. Finally, a stable fixed point (so-called sink) has always a vanishing index 20 .
• Invariant manifold. Suppose φ is a flow on the phase space (X, t) equipped with a spatial metric dist. A family A := τ∈t A τ of nonempty subsets of X is invariant with respect to φ if for all τ ≥ τ 0 A goes by the name of invariant manifold.
12(a)-Forward attraction. Let φ be a flow (process). A nonempty, compact and invariant manifold A := τ A τ for compact I is said to forward attract if for all x 0 ∈ X and τ 0 ∈ t. The distance function between two sets A, B is defined as usual on the metric space X as dist(A, B) = inf x∈A,y∈B dist(x, y).

12(b)-Pullback attraction. A nonempty, compact and invariant non-autonomous set
for all x 0 ∈ X and τ ∈ t. 20 The imaginary part of λ ∞ if present, means that the flow line is going to spiral its way in or out; thus the name stable or unstable spiral, respectively.
13-Basin of attraction. A neighborhood of an attractor is called the basin of attraction B if it consists of all points that evolve to the attractor in the limit τ → ∞. In other words, B is the set of all points b ∈ B in the phase space provided that for any open neighborhood N of the attractor, there is a positive constant s such that φ τ,τ 0 (b) ∈ N for all real τ > s.
14-Lyapunov function. Let us assume that the point (x, τ) at τ → ∞ is (0, ∞) which defines an asymptotically stable fixed point of a non-autonomous dynamical system. Then if there is a function V (x, τ) such that the Lyapunov stability conditions are satisfied, V (x, τ) is called a Lyapunov function.
It may be appropriate here to mention that in general, 11(a)-(b) are very different and independent concepts. Pullback attractor (B8) contains information about the past (early-time) of a non-autonomous dynamical system, whereas forward attractor (B7) makes use of information about the future. It well may be the case that one exists while the other one does not in a nonautonomous system. However, in the case of an autonomous system, they are equivalent and the resulting attraction is global with the set A referred to as global attractor.
Appendix C: Constructing the transseries for T (τ)

Algebraic properties of transseries
Let us suppose the following transseries is given where φ (n) k are complex-valued coefficients. We define the basis of the transseries aka E (n) The transseries can be regarded as an element of a vector space defined as The operations on this vector space are briefly explained below. Note that the transseries is always closed under these operations.

Product
The product operation × between two transseries is an additive operation on both transmonomial and asymptotic orders given by By assuming the ansatz for c 01 (τ) as we can get from Eq. (C23) Hence, it is easy to read off the coefficients for any n and k = 0 for n = 0 and k ≥ 1 Therefore, inserting Eq. (C27) in (C22) results in where Ψ(τ) converges to zero at late times, that is limτ →∞ Ψ(τ) = 0.

The c-T dynamical system
In the analysis of the previous section, we used w as a time variable, but something else that may be done is to consider the same analysis for the dynamical system including explicitly the temperature, which will require constructing transseries for T (τ) as well. In this subsection, we seek to prove that this dynamical system entails the same structure as (11), which is compatible with Costin's prepared form.
The dynamical system in the τ-coordinate takes the following form where B is a constant matrix, D is a diagonal constant matrix, and A is a constant vector. From Eq. (C32), one can obtain the formal solution of T (τ) as where C 0 is an integration constant, and C T (τ) is defined as In the IR, or the late-time limit τ → ∞, T (τ) is observed to behave like with a constant T 0 . Hence, by redefining this constant, the temperature T (τ) takes the form Now, by working (C36) into (C31) and changing the parametersτ andθ 0 aŝ the dynamical system becomes whereT Note that the temperature T (τ) can be reproduced by As discussed in App. C, we find thatT (τ) will be expressed bŷ where Ψ(τ) is a continuous function with a transseries formula identical to that ofc, which converges to zero asymptotically. Therefore, Eq. (C38) bears the familiar prepared form as in Eq. (11) with only one significant difference from a physical standpoint: the temperature now is a dimension in the phase space of the dynamical system. In the UV limit, the divergence of T (τ) as τ → 0 in general (see Eq. (C32)) suggests that T = 0 and T = 0 theories are disconnected from each other, and there cannot be any flow line (process) that connects a UV fixed point to the hydrodynamic fixed point achieved at τ → ∞ in the original kinetic model of the Bjorken flow.
Appendix D: Bjorken dynamical system and its trasseries solutions for general τ r = θ 0 /T 1−∆ In this appendix, we extend the dynamical system in Eq. (9) to include general relaxation time, τ r = θ 0 /T 1−∆ with 0 ≤ ∆ ≤ 1, which is very straightforward. The generalized system has the same number of asymptotic fixed points with their stability being identical to the case ∆ = 0 as discussed in Sect. IV A. So except for some small quantitative alterations in the flow lines around UV and IR, the qualitative picture of the phase portrait remains the same. In other words, there are no changes in global dynamics. The generalized dynamical system is described by As in the main bulk of the paper, we change to the coordinate w = τT 1−∆ and consider the dynamical system to be of the form where the matrices and vectors involved here have the same definitions as in (41). Being again compatible with Costin's prepared form, we conclude that the transseries can be explained with the same form used to explore the solutions of the ∆ = 0 model defined by (15). By defining an invertible matrix U for the diagonalization purposes again, one finds the recursive relation involving the transseries data as with the IR data S i andb i in the transseries being where b i are the eigenvalues of the matrix B.
In addition, the transaymptotic mathing is obtained in the following, whereζ i = ∂ /∂ log ζ i . Summing over m yields the transasymptotic matching for the general ∆ as 20 (1 + ∆/2)(b ·ζ − k) + b i C i,k − (1 + ∆/2)S ·ζC i,k+1 + 3 2θ 0C i,k+1 + 20Ã i δ k,0 C 1,k−k b ·ζ − k C i,k + 20 The transseries and leading order asymptotics of five moments c 01 , . . . , c 21 for ∆ = 1 system are depicted in Fig. 10. To avoid repeating ourselves here, we refer the interested reader to the explanations given in the text prior to Fig. 1. 1. Transasymptotic matching for c 01 of the N = 0, L = 1, ∆ = 1 dynamical system As a result of Eq. (D6), the mode-to-mode coupling term c 01 dc nl dw vanishes, and therefore the transasymptotic matching effectively becomes a finite expansion of the transmonomials packaged into ζ at each order. The first five transseries coefficients are then found to be exactly given bỹ In Eq. (E5) the reader must bear in mind that the Landau matching condition for the energy density at τ = τ 0 when using the RS distribution function (E2) implies T 0 = R 1/4 (ξ 0 )Λ 0 , where the function R(ξ ) is defined by [20] R(ξ ) = 1 2 The momentsM nl (114) can now be evaluated from the exact solution (E1) as follows where F nm 0 (τ, τ ) = F nm eq (τ, τ ) = The exact solutions for both c nl andM nl are determined numerically by having the expression for the temperature from Eq. (E3). by varying its dimension from 1 to L and compare against exact result (4) by using the same set of initial values of the moments c l . The initial conditions for the exact c nl (E6) are chosen to be τ 0 = 0.25 fm/c, T 0 = 0.6 GeV and ξ 0 = 1000. In Fig. 16 we present the numerical results for our findings by considering the initial conditions c 0l (w 0 ) = c 0l (w 0 ; ξ 0 ). When the dimension of the truncated dynamical system increases, they get closer to the exact solutions of the full RTA-BE. We notice that the larger the index l, the more moments are needed in order to match the exact solution of RTA-BE. In Fig. 17, the exact moments c nl (n = 1, 3, 5) are compared to the solutions of the truncated system. This figure encodes the information of the non-linear structure of the original dynamical system: the evolution equation of c nl couples only to the moments of the same or lower order in n, and at the same time it couples to the next moment c n,l+1 . So c nl in Fig. 17 become closer to the exact RTA result if L is increased for a fixed value of n. We comment that the numerical results presented in Fig. 16 and 17 do not depend on the particular choice of initial conditions as long as we are bound to the basin of attraction of the invariant manifold of the truncated system.