The non-equilibrium attractor for kinetic theory in relaxation time approximation

I demonstrate that the concept of a non-equilibrium attractor can be extended beyond the lowest-order moments typically considered in hydrodynamic treatments. Using a previously obtained exact solution to the relaxation-time approximation Boltzmann equation for a transversally homogeneous and boost-invariant system subject to Bjorken flow, I derive an equation obeyed by all moments of the one-particle distribution function. Using numerical solutions, I show that, similar to the pressure anisotropy, all moments of the distribution function exhibit attractor-like behavior wherein all initial conditions converge to a universal solution after a short time with the exception of moments which are sensitive to modes with zero longitudinal momentum and high transverse momentum. In addition, I compute the exact solution for the distribution function itself on very fine lattices in momentum space and demonstrate that (a) an attractor for the full distribution function exists and (b) solutions with generic initial conditions relax to this solution, first at low momentum and later at high momentum.

The hydrodynamization of the QGP can be understood by studying the approach of different initial conditions to a non-equilibrium attractor [38]. For example, one finds in practice that the pressure anisotropy (P L /P T ) developed from a set of different initial conditions all approach a universal attractor solution in a time on the order of a few fm/c with the time scale increasing as one decreases the effective coupling or initial temperature. Once the solutions have collapsed onto the attractor, information about the details of the initial condition used are wiped out, suggesting a kind of pseudo-thermalization of the system [47]. One shortcoming of prior works is that, when discussing the attractor, past authors (including this author) have restricted their attention to the amplitude, ϕ, which can be related to a particular ratio of a linear combination of low-order moments of the oneparticle distribution function. A natural question is whether or not attractor-like behavior is observed in the evolution of all moments of the one-particle distribution function and, if an attractor for higher-moments exists, what is the timescale for a generic initial condition to approach the attractor for a given moment.
In this paper, I address this question by considering the behavior of general moments of the one-particle distribution function in the case of a 0+1d conformal system subject to boost-invariant longitudinal Bjorken flow and a relaxation-time approximation (RTA) collisional kernel. For this purpose, I make use of a semi-analytic exact solution of the conformal RTA Boltzmann equation which has been obtained previously [51,52]. This exact solution takes the form of a one-dimensional integral equation for the energy density which can be solved iteratively and mapped to a solution for the effective temperature of the conformal system. Once the effective temperature is known, one can solve for all moments of the one-particle distribution function and the distribution function itself. I will present results for a subset of moments which are indicative of the general behavior seen.
I find that there exists attractor-like behavior in all moments, however, low-order moments which contain only powers E and not p z exhibit a slower approach to their respective attractors. This is explained in terms of a spectrum of modes with nearly zero p z and an average p T which increases with time. These modes have a small impact on the overall dynamics but dominate these particular moments. For all other moments, I find that the timescale for approach to the attractor only weakly depends on the order of the moment considered. The convergence of all modes to their respective attractor on approximately the same time scale suggests that the evolution of the full one-particle distribution function exhibits attractor-like behavior. I demonstrate that it does by solving for the full oneparticle distribution function associated with the attractor and then I study the flow of the distribution function towards the distribution function's attractor.
The structure of this paper is as follows. In Sec. 2, I review the basic setup for finding the exact solution to the 0+1d RTA Boltzmann equation. In Sec. 3, I present the integral equation obeyed by a general moment M nm of the one-particle distribution function and demonstrate that, for low-order moments, it reduces to results obtained previously in the literature. In Sec. 4, I calculate the moments using anisotropic hydrodynamics and secondorder viscous hydrodynamics. In Sec. 5, I present my numerical results and discuss. In Sec. 6, I present my conclusions and an outlook for the future.

Setup and review
In this section I review how to obtain the exact solutions to the 0+1d RTA Boltzmann equation. The method is based on the original work presented in Refs. [51,52]. 2

RTA Boltzmann equation
My starting point is the RTA Boltzmann equation with The quantity τ eq = 5η/T is the relaxation time whereη = η/s is the shear viscosity to entropy density ratio and T is the local spacetime-dependent effective temperature, which is proportional to the fourth root of the local energy density. For a conformal system, the equilibrium distribution function f eq may be taken to be a Bose-Einstein, Fermi-Dirac, or Boltzmann distribution. Herein, I will assume that f is given by a Boltzmann distribution 3) The effective temperature T can be obtained via the Landau matching condition which demands that the energy density calculated from the distribution function f is equal to the energy density determined from an equilibrium distribution, f eq . The quantity u µ is the four-velocity of the local rest frame of the matter (fluid four velocity). I will assume Bjorken flow, in which case the Minkowski-space components of the four-flow are u µ = (t/τ, 0, 0, z/τ ), where τ = √ t 2 − z 2 is the longitudinal proper-time. In Milne coordinates, Bjorken flow is static, i.e. u τ = 1 and u x,y,ς = 0.
The use of this simple form of the kinetic equation given by Eqs. (2.1) and (2.2) is motivated by the fact that there are many results obtained within this approximation, allowing one to make comparisons with other approaches. In particular, there exist consistent second and third-order calculations of the kinetic coefficients in RTA, for example, see Refs. [58][59][60][61][62][63]. Finally, and perhaps most importantly, in this simple case it is possible to solve the kinetic equation exactly using straightforward numerical algorithms.

Thermodynamic functions
For a single species of scalar massless particles obeying equilibrium classical statistics the particle density, entropy density, energy density, and pressure are In what follows I make use of the relation ε eq = 3P eq when a specification of the equilibrium equation of state is required.

Boost-invariant variables
In the case of one-dimensional boost-invariant expansion (0+1d), all scalar functions of spacetime depend only on the longitudinal proper time τ . To proceed, in addition to the timelike flow four-vector u µ one can introduce a spacelike vector that is orthogonal in all frames and corresponds to the z-direction in the local rest frame of the matter, z µ LAB = (z/τ, 0, 0, t/τ ) and z µ LRF = (0, 0, 0, 1). The one-particle distribution function f (x, p) is a scalar under Lorentz transformations. The requirement of boost invariance implies that in this case f (x, p) may depend only on three variables: τ , w, and p T [64,65]. The boost-invariant variable w is defined by (2.5) Note that, in the prior equation, z is the spatial coordinate, not to be confused with the basis vector z µ . With the help of w and p T one defines Using (2.5) and (2.6) one can easily find the energy and the longitudinal momentum of a particle The momentum integration measure is In the following I will consider massless partons, m = 0.

Boost-invariant form of the kinetic equation
Using the boost-invariant variables introduced in the previous section, one finds [51,52] p Using Eqs. (2.9) in Eq. (2.1) and simplifying, the 0+1d Boltzmann equation takes a particularly simple form [51,52] ∂f where the equilibrium distribution function can be written as (2.11) In the results section, I make use of the fact that symmetries require that f (τ, w, p T ) be an even function of w and that it depends only on the magnitude of the transverse momentum p T .

Exact solution
Once cast in the form (2.11), the solution becomes straightforward, with the result being [51,52,66] f (τ, w, where D is the damping function (2.13) At τ = τ 0 the distribution function f reduces to the initial distribution function, f 0 . For the conformal RTA Boltzmann equation, one has [58,59] τ eq (τ ) = 5η T (τ ) , (2.14) whereη ≡ η/s is the ratio of the shear viscosity to entropy density. In order to solve for f using Eq. (2.12) one needs to know T (τ ). An efficient way to do this is by considering the integral equation for the energy density obtained by integrating both sides of Eq. (2.12) times (p·u) 2 , i.e. ε = dP (p·u) 2 f . After performing this operation on Eq. (2.12), the left hand side becomes the non-equilibrium energy density, ε(τ ), while the first term on the right hand side encodes the free-streaming contribution to the energy density evolution, and the second term on the right hand side dominates the late-time behavior for any τ eq = ∞. To close the equation, one implements the Landau matching condition ε(τ ) = ε eq (T ) on the left-hand-side so that one obtains an integral equation for ε eq (T ) ∝ T 4 (τ ). Before presenting this, however, one needs to specify the form of the initial distribution function f 0 appearing in Eq. (2.12) in order to fix the free streaming contribution.

Initial distribution function
In this paper, I consider initial conditions of Romatschke-Strickland form [3] with a Boltzmann distribution as the underlying isotropic distribution This reduces to an isotropic Boltzmann distribution if the anisotropy parameter ξ 0 = ξ(τ 0 ) vanishes. In this case, the transverse momentum scale Λ 0 is equal to the system's initial temperature T 0 . By direct calculation one obtains where H(y) = y π 0 dφ sin φ y 2 cos 2 φ + sin 2 φ , and α 0 = α(τ 0 ) denotes the initial value of the anisotropy parameter with 0 ≤ α ≤ ∞.

Integral equation for the effective temperature
The resulting integral equation for the effective temperature is [51,52] This equation can be solved iteratively by constructing a discrete lattice in proper-time, making an initial guess for T (τ i ), and using numerical quadratures to compute the required integrals on the right hand side. I will provide details of the lattice size, spacing, tolerance, etc. in the results section. Once the solution for T (τ ) is obtained by iterating this procedure, it can then be used in Eq. (2.12) to reconstruct the full one-particle distribution function point-by-point in momentum space. In principle, one could compute f and then determine moments from that, however, this proves to be quite inefficient and computationally demanding. Instead, one can derive integral equations which can be used to determine the evolution of all moments of the distribution function.

Evolution equation for a general moment
To extend the treatment presented in the previous section to general moments, one introduces In principle, powers of p 2 T could also appear in a general moment, however, such moments can be expressed as a linear combination of the two-index moment appearing above using . For n = 1 and m = 0 one obtains the number density For n = 2 and m = 0 one obtains the energy density and for n = 0 and m = 1 one obtains the longitudinal pressure For a conformal system, one can use ε = 2P T + P L to determine the transverse pressure. This follows by using the mass shell condition, . In the general case, using the boost-invariant variables introduced earlier, the moment integral necessary becomes Taking a general moment of Eq. (2.12) one obtains Using the results obtained above, one can write this as Eq. (3.6) is one of the main results obtained herein.

Cross checks
As a check of the general moment equation, one can verify that it reduces to known equations for the low-order moments of the distribution function available in the literature. Taking the n = 2 and m = 0, and relabeling H 20 → H, one obtains Using the fact that M 20 (τ ) = ε(τ ) = ε eq (τ ) = 3T 4 (τ )/π 2 , this reduces to the known integral equation for the effective temperature [51,52] Taking n = 0 and m = 1, and relabeling H 01 → H L , one obtains the longitudinal pressure [51,52]

Evolution of the moments in dissipative hydrodynamics
In the results section I will compare results obtained using Eq. (3.6) with those obtained using anisotropic hydrodynamics (aHydro) [67][68][69][70] and second-order viscous hydrodynamics (vHydro) [58,59]. In anticipation of this, I now list the results for a general moment in both schemes. For aHydro the moments depend only on the anisotropy parameter α and for vHydro they depend only on the ratio of the shear viscous correction π = π ς ς to the energy density.

Evolution of the moments in anisotropic hydrodynamics
In aHydro, one uses the first and second moments of the Boltzmann equation to solve for the evolution of α and Λ with the one-particle distribution function assumed to be of the form [67,68] (4.1) Using this form, one finds with α(τ ) = 1/ 1 + ξ(τ ).
When comparing results it is useful to rescale each moment by its equilibrium value For aHydro, this gives where I have introduced the scaled moments .

(4.5)
It is straightforward to verify that in the isotropic (equilibrium) limit, one has M nm aHydro (τ ) = 1.

Evolution of the moments in second-order viscous hydrodynamics
In vHydro, the 0+1d one-particle distribution function takes the form [58,59] whereπ = π/ε. Computing M nm using this form, one finds withπ ≡ π/ε. Computing the scaled moment, one obtains Note that the Navier-Stokes (NS) result for a general moment can be obtained from the above expression by takingπ = 16η/(9τ T ). Also note that from the above equation one finds that M n0 vHydro (τ ) = 1.

Navier-Stokes thermalization time
Before proceeding, I mention that using the above relation and taking the Navier-Stokes limit, one can solve for the scaled time As one can see from this result, the Navier-Stokes limit predicts that large m and n moments will equilibrate at a later time than small m and n moments. In the results section, I will compare this result to the thermalization time extracted from the exact solution to the Boltzmann equation. 3

Numerical results
I will now present some representative numerical solutions using different initial conditions and the attractor solution to which they flow. For a given value ofη and set of initial conditions specified by α 0 and T 0 , I solve the integral equation (2.19) numerically. For this purpose I wrote a CUDA-based GPU code which allows one to efficiently solve the integral equation (2.19) efficiently on very large lattices. The code uses a logarithmically-spaced grid in proper-time in order to more accurately account for the effects of large early time gradients. The code is included in the arXiv bundle for this paper and is also publicly available for download using the link provided in Ref. [71]. Once the solution for the effective temperature is obtained, one can use Eq. (3.6) to obtain any moment required. One can also use Eq. (2.12) to reconstruct the full one-particle distribution in a grid in momentum-space. For all results presented here, I iterated the integral equation for the temperature (2.19) until the result converged to sixteen digits at all values of τ . I analyze the (pseudo-)thermalization of the system by considering the scaled moments (4.3) as a function of the scaled time w ≡ τ /τ eq = τ T /5η. 4,5 The rate at which the scaled  moments approach unity provides information about the thermalization of the system. 6 Some familiar quantities like the scaled number density and longitudinal pressure are given by M 10 (τ ) = n(τ )/n eq (τ ) and M 01 (τ ) = P L (τ )/P eq (τ ), respectively.
In Fig. 1, I compare the attractor (black solid line) with a set of representative solutions (dashed/dotted colored lines) with differing levels of initial momentum-space anisotropy (0.1 ≤ α 0 ≤ 1.5) and fixed initial temperature T 0 = 1 GeV at τ 0 = 0.1 fm/c. For the solutions with different initial conditions (dashed/dotted colored lines) I used 2048 points spaced logarithmically between τ = 0.1 and 100 fm/c. For the attractor solution, I used 4096 points spaced logarithmically between τ = 0.001 and 1000 fm/c and tuned the initial anisotropy to α 0 0.0025 following a method similar to the one outlined in the appendix of Ref. [40]. As can be seen from this figure, all solutions approach the attractor solution in a finite time. The slowest approach appears to be for moments with m = 0 which appear in 6 One exception is M 20 which equals unity at all times due to energy conservation. Figure 2. Visualization of the one-particle distribution function obtained using a typical (nonattractor) anisotropic initial condition. I will make these statements more quantitative shortly, but before that I would like to discuss the source of the slow thermalization and approach to the attractor for the moments with m = 0. This subset of moments contains no power of p z in the integrand and so are sensitive to the behavior of the one-particle distribution in the vicinity of p z = 0. In order to understand this behavior better, one can obtain the full one-particle distribution from the exact solution using Eq. (2.12). For this purpose, I first obtain the solution to the integral equation for the temperature (2.19) and then evaluate Eq. (2.12) in a grid in p x,y /T and p z /T . For the results show herein, I used a 500 × 500 grid in the scaled transverse and longitudinal momentum with |p T,z |/T ≤ 50.
In Fig. 2 I show snapshots of the one-particle distribution function at three different scaled times. For this figure, I used a typical anisotropic initial condition from the set shown in Fig. 1. As can be seen from this figure, the exact solution for the one-particle distribution function generically contains two visually identifiable components. The first is an anisotropic piece which becomes increasingly more compressed into the region with . The second, more isotropic, component which can be seen in Fig. 2 comes from the second term in (2.12). This contribution dominates at late times. Turning now to Fig. 3, in this figure I present similar snapshots of the one-particle distribution function, however, this time using the temperature evolution obtained from the attractor solution to Eq. (2.19). In this case, one sees similar features to the generic solution presented previously, however, one sees that in this case the attractor initial condition used and the free streaming term result in a set highly-squeezed modes that have nearly p z = 0. It is precisely this set of squeezed modes which cause the approach to equilibrium to proceed more slowly for these moments than for moments with m = 0. For m = 0, the powers of p z naturally reduce the impact of the highly squeezed modes. For m = 0 the moments will be dominated by these squeezed modes at early times and, as n is increased, the magnitude of M n0 will increase for all n > 2 at early times. As an example, in Fig. 4, I plot the integrand of the M 30 moment as a function of p T /T for p z = 0. As this figure demonstrates, at early times, the integrand for this moment is dominated by modes with high transverse momentum. The peak of the integrand moves to higher p T and has an increasing magnitude at early times, however, at late times, the amplitude of the peak is diminished due to the damping function which appears in the free streaming contribution to the one-particle distribution function. Despite the presence of these highly squeezed modes, the dynamics of the longitudinal pressure and other low-order moments with m > 0 are not significantly affected.
Returning to the behavior of the general moments obtained from the exact solution, in Figs Fig. 6, one sees that the prediction that M n0 should increase in magnitude at early times as n is increased due to squeezed modes is borne out by the results. The other panels demonstrate that, for fixed m, increasing n results in slower thermalization. Again this is consistent with the slower relaxation of high momentum modes.
To make this more quantitative, one can extract a "thermalization time" w therm by solving for the scaled time at which a given moment falls within 10% of its equilibrium value, i.e. when the M nm = 0.9. The results obtained are shown in Fig. 7. For this plot I exclude M 20 , since this is always equal to unity by energy conservation. In addition, if the thermalization time was longer that the maximum rescaled time (in this case w max 112), then I don't plot a value. As Fig. 7 demonstrates, the thermalization time increases as both n and m are increased. Again, this is indicative of the slower thermalization of high-   Fig. 8. In Fig. 8, the top left panel shows w c as a function of m with the set of lines corresponding to different values of n. As can be seen from this figure, all moments with m ≥ 1 have short pseudo-thermalization times that decrease with increasing m in the range shown. This is surprising since in Fig. 7 one sees that increasing m for fixed n results in larger thermalization times. Turning to the top right panel of Fig. 8, one sees once again that, for the moments shown, one has w c < w therm , with their separation increasing as one increases m and n for m > 2. The moments with m = 0, 1, 2 are found to have pseudo-thermalization times that increase with increasing n. Because of the limited range shown, it is possible that other moments might start increasing at very large n. Finally, because of their slow relaxation, in the bottom panel of Fig. 8 I plot w c as a function of n for moments with m = 0, separately. 7 From this panel one sees evidence that, for m = 0, the pseudo-thermalization time increases approximately linearly large n. Comparing with Fig. 7 one sees that, for the m = 0 moments, it is possible to have w c ∼ w therm or to even invert their order.
Summarizing, one sees that in the range 0 ≤ n ≤ 8, 1 ≤ m ≤ 8 one has w c < w therm with the separation increasing with m and n. From the results, one also sees indications that for 1 ≤ m ≤ 2, w c increases with n for large n, but slower than w therm increases. For moments with m = 0, one sees large w c which can become on the order of or exceed w therm . Based on this, it's hard to draw firm conclusions about the nature of the (psuedo)thermalization of the system for all possible moments. To shed some light on this question, one can compute the full distribution function for the exact solution using the attractor initial conditions. As mentioned previously, once the solution to the integral equation for the temperature (2.19) is obtained, one can evaluate Eq. (2.12) in a grid in p x,y /T and p z /T .
The result for the distribution function associated with the attractor is visualized in Fig. 3, however, it's hard to draw quantitative conclusions from such plots. In Figs. 9 and 10, I present the distribution function along the line with p T = 0 as a function of p z (top row) and along the line with p z = 0 as a function of p T (bottom row). Fig. 9 shows the lowmomentum region p x,y,z /T ≤ 3 and Fig. 9 shows the high-momentum region p x,y,z /T ≤ 40. Focusing on Fig. 9, firstly one sees that the two typical solutions (red short-dashed and blue long-dashed lines, respectively) converge towards the attractor solution (black solid line), however, the approach is not uniform in the sense that in some regions of momentum the solutions approach the attractor from below while in other regions they approach it from above. Secondly, one sees that the approach to attractor and thermal solutions (black solid and green dot-dashed lines, respectively) is slower along the transverse line (p z = 0) than along the longitudinal line (p x,y = 0). This is consistent with there being a subset high 7 For M 20 , I chose wc to be the minimum time since M 20 = 1 at all times by Landau matching. momentum modes with high p T and nearly p z = 0 which equilibrate more slowly than the rest of the modes.
This can be seen more clearly if one zooms out and investigates the behavior of the distribution function at higher momentum, as shown in Fig. 10. From this figure, one sees that indeed the distribution in the longitudinal direction (top row) approaches the attractor solution and thermal approximation very quickly, whereas along the transverse direction (bottom row) one sees a two-component distribution function, with a subset of low-momentum modes which quickly approach the attractor being visually distinct from modes with high-momentum that have a much more shallow exponential slope. In addition, one sees a sharp transition from the low-momentum to the high-momentum behavior and that this transition point, p transition , moves to higher momentum as a function of w.
To investigate the behavior of the transition point p transition more quantitatively one can take the ratio of the distribution function evaluated along the line with p x,y = 0 (longitudinal direction) to the same along the line with p z = 0 (transverse direction), i.e. f (0, 0, p)/f (p, 0, 0). I present plots of this ratio in Fig. 11. As can be see from this figure, there is an isotropization "front" that moves to the right and which can be associated with the position of the transition point p transition in the distribution function visible in the bottom row of Fig. 10. If one defines the transition point as the point at which the longitudinal to transverse ratio along the transverse directions equals 1/2, one finds that, for the attractor solution, at w 2 the hydrodynamization front moves out in temperature-scaled momentum at a constant speed of 1.5, i.e. p transition /T 1.5 w. A similar hydrodynamization front can be seen in the typical solutions (dashed lines), however, in this case one finds that the hydrodynamization front speed is faster than 1.5 and approaches 1.5 at asymptotically late times as a given typical solution converges to the attractor solution.
Finally, I would like to compare the exact results obtained for a subset of moments with results obtained from the aHydro and vHydro dissipative hydrodynamics approximations. For vHydro, I determine the attractor using the complete second-order viscous hydrodynamics equations of Denicol, Niemi, Molnar, and Rischke (DNMR) [58,59]. For aHydro, I use the method introduced originally by Florkowski and Tinti [69] which utilizes the first and second moments of the Boltzmann equation. For details concerning the determination of the attractor for both aHydro and vHydro, I refer to the reader to Ref. [47]. In both cases, the attractor is determined from the solution of a one-dimensional ordinary one-dimensional ordinary differential equation subject to the appropriate initial condition. For vHydro, one extractsπ = π/ and, using this, one can reconstruct the solution for any moment using Eq. (4.8). From the vHydro result, one can also obtain the Navier-Stokes (NS) result by takingπ = 16η/(9τ T ). For aHydro, one extracts α associated with the attractor solution and uses Eq. (4.4) to compute the moments. In Fig. 12, I compare the exact attractor (black solid lines) with the aHydro attractor (red dashed lines), DNMR attractor (blue long dashed lines), and the NS limit for each moment (green dot-dashed lines). In all cases shown, aHydro provides a better approximation to the exact moments than the vHydro or NS solutions.
In the case of aHydro, one sees that the worst agreement is for the m = 0 moments (leftmost column of Fig. 12). This can be understood from the fact that within aHydro one assumes that the one-particle distribution function is given by a single spheroid in momentum space. As a result, aHydro fails to accurately describe the evolution of the modes with m = 0 which are dominated by the free-streaming part of the evolution and remain highly anisotropic at all times. For m > 0 one sees that, as m and n are increased, the aHydro results differ more and more from the exact solutions in the transition region around w ∼ 1 − 10, however, the aHydro solution does not appear to "break" in the sense of giving unphysical (e.g. negative) results for any of the moments considered. 8 Turning to the vHydro and NS results, firstly one sees that for m = 0 both schemes predict M n0 = 1 at all times. This is by construction in vHydro since the δf correction to the vHydro distribution function vanishes in order to guarantee that the number density and energy density are unaffected by the viscous correction. As a consequence, vHydro and NS do not provide reliable approximations for these moments for n = 2. sees that, although the vHydro and NS results do a reasonable job in describing the low n and m moments, as n and m are increased these approximations become significantly worse. In fact, one sees that for all moments with m > 0, both vHydro and NS predict that the various moments become negative at early times. This is due to the breakdown of the near-equilibrium assumption and violates positivity of the moments from Eq. (3.1). The problem becomes more severe as one increases m and n due to the increasing powers of p z in the moment integrands. Relatedly, in Fig. 13, I compare the thermalization time w therm obtained from the exact solution (see Fig. 7 for the values of n for each curve) with the NS solution (left panel) and the aHydro solution (right panel). As can be seen from this figure, as m and n increase standard viscous hydrodynamics diverges from the exact solution, whereas aHydro continues to provide a reasonable approximation for the thermalization time determined from the exact solution for each moment. have dubbed the "pseudo-thermalization time", w c . I demonstrated that, within RTA, the pseudo-thermalization time depends primarily on the transverse momentum cut considered and that the longitudinal momenta converge to their attractor on a much faster time scale. The fact that the transverse modes isotropize more slowly was shown to originate from the free streaming contribution to the exact solution, which contributes a narrow (highly momentum-space anisotropic) band of modes in the vicinity of p z ∼ 0 which dominate moments of the one-particle distribution function which do not contain powers of p z , i.e. M n0 .
Contrarily, I demonstrated that moments with m > 0 have pseudo-thermalization times that are parametrically shorter than the corresponding thermalization time for 1 ≤ m ≤ 8 and 0 ≤ n ≤ 8. I presented a detailed comparison of the moments in this range and provided numerical estimates for both the thermalization and pseudo-thermalization times as a function of m and n.
Due to the fact that evolution of a large set of moments is sometimes difficult to interpret, I also presented solutions for the distribution function obtained from the exact RTA attractor solutions. Using these solutions, I demonstrated that there is an hydrodynamization front that propagates out along the p z = 0 plane from low-momenta to high-momenta with a constant speed, p T /T 1.5 w. With this in hand, the behavior of the various moments becomes easier to understand, in particular, the behavior of the M n0 moments.
Based on these findings, one can determine the time scale for momenta up to p T ∼ 2 GeV to thermalize. Assuming a hadronic freeze-out temperature of T FO = 150 MeV, then requires that momenta p T /T 13 are approximately thermalized. Using the propagation speed for the hydrodynamization front then gives w 9. To convert this to a physical time, one can use the exact solution with the following (typical LHC) initial conditions τ 0 = 0.25 fm/c, T 0 = 0.6 GeV, and η/s = 0.2, giving τ 8 fm/c. One can compare this to the time it takes the system to reach the freeze-out temperature, which using the same initial conditions, is τ FO 23 fm/c. Relatedly, one can determine the maximum transverse momentum which is thermalized at freeze-out. Using the same initial conditions, one finds w FO 17.3 giving p T 3.9 GeV. Finally, I presented a comparison of the attractors for moments of the one-particle distribution function using the aHydro, DNMR second-order vHydro, and Navier-Stokes approximations to the exact attractor solution. For all three hydrodynamic schemes, I was able to find compact analytical expressions that could be used to compute any moment based on the exact solution for the pressure anisotropy (or equivalently the amplitude ϕ) in the corresponding scheme. I found that the aHydro attractor for a general moment was the best approximation to exact attractor and that it also provided reliable estimates for the thermalization time of higher-order moments. One obvious shortcoming of aHydro is evident from this work, namely the use of a single spheroid. The exact solution for the distribution function has a two-component form with a highly-anisotropic free-streaming contribution and isotropizing component. It would be interesting to see if one could generalize aHydro to include such a two-component distribution. Turning, in the end, to firstand second-order viscous hydrodynamics, I found that, while they provide good approximations for low-order moments, higher-order moments are poorly reproduced and the resulting estimates for the pseudo-thermalization time for the higher-moments do not agree with the exact solutions. For this reason, when considering higher moments, one can no longer associate the pseudo-thermalization time with the "hydrodynamization time".
Looking to the future, it would be interesting to understand more about the numerically observed hydrodynamization front. I presented evidence that the propagation speed for the attractor solution is 3/2 in dimensionless units. It would be nice to have an analytical understanding of this result. Along these lines, Heller, Kurkela and Spalinski [72] demonstrated that the solutions to the exact RTA kinetic integral equation (2.19) possess infinitely many transient modes that carry the majority of information about the distribution function entering the energy density at late times. Since the evolution of the distribution function (2.12) can be written entirely in terms of the solution to (2.19), the modes identified by Heller et al should also govern the evolution of the full distribution function and, hence, all moments. This seems to provide a natural explanation of the emergence of an attractor for all moments. It would be interesting to study the evolution of higher moments and the distribution function itself using their approach. Other natural extensions of this work include, for example, studying the dependence on the underlying statistics of the isotropic distribution function similar to Ref. [54,63] and studying the effect of enforcing number conservation for a single component scalar field or a multi-component quark-gluon system, as done in Ref. [45].

Corrections
In this version (arXiv v3), I have corrected an error which affected Figs. 2, 3, 9, 10, and 11. This stemmed from a bug in the original reconstruction of the full distribution function. All other figures remain unchanged. I thank the authors of Ref. [73], in particular C. Chattopadhyay, for calling this to my attention. Luckily, all discussions and conclusions are unaffected by the change in the content of these figures. For accuracy of the description, I have taken this opportunity to change the nomenclature "isotropization front" to "hydrodynamization front".