An improved anisotropic hydrodynamics ansatz

We introduce an improved form for the anisotropic hydrodynamics distribution function which explicitly takes into account the free-streaming and equilibrating contributions separately. We demonstrate that with this improvement one can better reproduce exact results available in the literature for the evolution of moments of the distribution function, in particular, for moments which contain no powers of the longitudinal momentum in their definition (m=0 moments). Using the resulting dynamical equations, we extract the non-equilibrium attractor associated with our improved aHydro ansatz and demonstrate that the improvement also allows one to better reproduce the exact dynamical attractor obtained using kinetic theory in the relaxation time approximation, particularly at early rescaled times and for m=0 moments.


I. INTRODUCTION
In the very early universe (a few microseconds after the Big Bang), the quark-gluon plasma (QGP) is believed to have existed where the density can reach values ten times higher than those of ordinary nuclei. It was speculated theoretically that one can reach these extreme conditions by colliding two heavy nuclei with ultrarelativistic energies. In this collision, the temperatures can be million times hotter than the core of the sun, and a fraction of the kinetic energies of the two colliding nuclei transform to heat the QCD vacuum within an extremely small volume. Because of the appearance of modern accelerator facilities, ultrarelativistic heavy-ion collisions (URHICs) be able to provide an opportunity to systematically create and study different phases of the bulk nuclear matter. In heavyion collision experiments at Relativistic Heavy-Ion Collider (RHIC) located at Brookhaven National Laboratory, USA, and Large Hadron Collider (LHC) at European Organization for Nuclear Research (CERN), Geneva, the new state of matter (the QGP) is widely believed created. Results obtained at RHIC energies and recently at LHC energies strongly suggested the formation of a quark-gluon plasma (QGP) which may be close to (local) thermodynamic equilibrium, albeit in a tiny volume (∼ 100 − 1000 fm 3 ). After the QGP is generated, it is expected to expand, cool, and then hadronize in the final stage of its evolution, with a QGP lifetime on the order of 10 fm/c in central collisions [1][2][3].
Heavy-ion collisions such as those at RHIC provide a primary tool to study the thermodynamic and transport properties of the QGP. Of remarkable importance is knowledge of time evolution of the rapidly expanding the QGP that produced in these URHICs. For this purpose, one can use a basic theoretical approach called relativistic hydrodynamics to describe the QGP. The resulting models describe the collective behavior of the soft hadrons with P T 2 GeV quite well. In early studies, it was found that the QGP created at RHIC energies was well described by models which assume ideal hydrodynamic behavior from very early times τ 1 fm/c [4][5][6][7][8]. Traditionally, one can apply ideal hydrodynamics if the system is in perfect isotropic local thermal equilibrium. Based on these early studies, it was expected that the QGP would isotropize on a timescale τ ∼ 0.5 fm/c. In practice, however, when one includes viscous corrections to the hydrodynamical models  one observes that at times τ 2 fm/c there can still be sizable differences between the transverse pressure, P T , and longitudinal pressure, P L which is associated with the existence of a non-equilibrium hy-drodynamic attractor . In addition, as one moves closer the transverse/longitudinal edges of the QGP, the size of the pressure anisotropies increases at all times [79][80][81]. Faced with this, researchers suggested to find another method to formulate hydrodynamics in a momentum-space anisotropic QGP. Recently, there have been theoretical and phenomenological studies that try to better account for large deviations from isotropy by relaxing the assumption that the QGP is close to local isotropic thermal equilibrium. To address this issue, they introduced a framework called anisotropic hydrodynamics (aHydro) in order to describe the non-equilibrium dynamics of relativistic systems, without breaking important physics constraints such as the positivity of the one-particle distribution function [82][83][84][85][86][87].
In a prior paper [88], comparisons between three hydrodynamic models and exact solutions of the RTA Boltzmann equation [89][90][91] were presented. It was found that linearized viscous hydrodynamics performed more poorly than the canonical formulation of aHydro in reproducing the exact attractor for all moments. However, although the canonical aHydro formulation [82,83] did a reasonable job in describing moments with m > 0, Ref. [88] found that it did not provide a good approximation for moments with m = 0. The failure of the canonical formulation was postulated to be due to the fact that the exact solutions to the RTA Boltzmann equation have an explicit two-component nature and cannot be accurately described by a single ellipsoidal form. As a result, it would be interesting to implement aHydro with a two-component ansatz for the distribution function to see if a better description of moments with m = 0 can be achieved. Additionally, it would be interesting to see if this also results in a quantitative improvement for higher-order moments.
In this document, we report on our progress in obtaining improved dynamical equations for anisotropic hydrodynamics through the use of an improved ansatz for the form of the underlying aHydro distribution function which explicitly includes a free streaming contribution. We demonstrate that with this improvement one can better reproduce exact results available in the literature for the evolution of moments of the distribution function, in particular, for moments which contain no powers of the longitudinal momentum in their definition (m = 0 moments). Using the resulting dynamical equations, we extract the nonequilibrium attractor associated with our improved aHydro ansatz and demonstrate that the improvement also allows one to better reproduce the exact dynamical attractor obtained using kinetic theory in the relaxation time approximation, particularly at early rescaled times and for m = 0 moments. We will focus our attention in this first work on a conformal system undergoing boost-invariant and transversally homogeneous Bjorken expansion, however, the method introduced herein is easily extended to full 3+1d.
The paper is organized as follows. In Sec. II we present the basic setup and assumptions used for the system and introduce our improved aHydro distribution function ansatz. We then use the first and second moments of the Boltzmann equation to obtain equations of motion for the dynamical parameters appearing in the new ansatz. We do this explicitly for a system undergoing boost-invariant 0+1d Bjorken expansion. Using the resulting dynamical equations we obtain the time evolution of all moments of the distribution function. In Sec.
III, we present our numerical results and discuss. In Sec. IV we present our conclusions and an outlook for the future.

II. SETUP
For the current work, we assume a system of massless particles. Furthermore, we assume that the system is undergoing boost invariant longitudinal expansion (v z = z/t) and expands only along the beam-line axis, ignoring the effects of transverse dynamics. Accordingly one can assume a homogeneous distribution in the transverse directions and set v x,y = 0. By taking into account these assumptions only proper-time derivatives remain and the dynamics reduces to 0+1d dimensional evolution [92]. In order to better describe free streaming contributions to the evolution of the one-particle distribution function, we propose an improved aHydro one-particle distribution function of the form where the first term is the free-streaming contribution and the second term is the equilibrating contribution. In Eq. (1), Λ 0 is the initial momentum scale, f 0 is the initial particle distribution with where ξ 0 is the initial momentum-space anisotropy, τ 0 is the initial proper time, and where RS indicates the anisotropic Romatschke-Strickland form [93]. The equilibrium distribution function f eq may be taken to be a Bose-Einstein, Fermi-Dirac, or Boltzmann distribution. Here we will assume that f eq is given by a Boltzmann distribution. For free-streaming distribution function f 0 is also of RS form but with ξ = ξ FS and Λ = Λ 0 , i.e.
Here we use the label '0' to emphasize that this contribution is constrained by the initial condition for the distribution function. Additionally, −1 < ξ < ∞ is a parameter that indicates the strength and type of momentum-space anisotropy. By stretching (−1 < ξ < 0) or squeezing (ξ > 0) the underlying isotropic distribution function f eq along one direction in momentum-space, one can obtain an anisotropic distribution function.

A. Moments of the improved distribution function
To calculate the energy density and pressures in the local rest frame (LRF), one can integrate the distribution function (1) times p µ p ν using the Lorentz-invariant integration After performing this operation one finds that all moments of the distribution function can be decomposed into two terms, e.g.
where in Eqs. (8) and (9) the left hand sides are the non-equilibrium energy density and longitudinal pressure, respectively. For a conformal system, one can use = 2P T + P L to determine the transverse pressure In general, one can compute a large set of moments of the one-particle distribution function Taking a general moment of Eq. (1), one finds Note that certain moments map to familiar hydrodynamics variables, e.g. taking n = 1 and m = 0, one obtains the number density Taking n = 2 and m = 0, one can evaluate the energy density via and taking n = 0 and m = 1, one obtains the longitudinal pressure Since both the free-streaming and equilibrating contributions are of RS form, we can compute the moments for both of these contributions using [61] M nm with where 2 F 1 is a hypergeometric function, y = 1/ √ 1 + ξ, and it has been assumed that the underlying isotropic distribution function is a Boltzmann distribution function. In practice, we will scale these moments by their equilibrium limit, which assuming Boltzmann statistics, Using the improved aHydro ansatz (1) one obtains where we have introduce the scaled moments .
Note that one has M nm aHydro (τ )=1 if the system is in equilibrium.

First moment
Our starting point is the Boltzmann equation for massless particles where the collisional kernel is taken to be the relaxation-time approximation (RTA) collisional kernel and u µ is the four-velocity associated with the local rest frame. Herein we will focus our attention on a system that is transversally homogenous and subject to boost-invariant Bjoken flow (0+1d). In order to preserve conformal invariance, the equilibration time must be inversely proportional to the local temperature and, for RTA, is given by [94,95] τ eq (T ) = 5η/T , whereη = η/s is the ratio of the shear viscosity η to entropy density s.
The first moment of the left-hand side of the Boltzmann equation reduces to ∂ µ T µν ; however, in the relaxation time approximation the first moment of the collisional kernel the right hand side results in a constraint that must be satisfied in order to conserve energy and momentum, i.e. dP p µ C[f ] = 0. This constraint is referred to as the matching condition and allows one to compute the local effective temperature of the system. In RTA, it results in the following constraint equation where the effective temperature T appears in eq .
As a result of this constraint, computing the first moment gives Expanding this equation out in terms of the non-vanishing components of the energymomentum tensor, for a 0+1d system, one obtainṡ withΛ = Λ/Λ 0 and R(ξ) = 1 2 which satisfy 3R = 2R T + R L .

Matching condition
At any time, we define the local effective temperature T (τ ) of the fluid using the canonical matching condition which results from the vanishing of the right-hand-side of the first moment of the Boltzmann equation. Using the improved form one finds with

Second moment
To close the system of equations, we use the uzz projection of the second-moment minus the 1/3 of the sum of uxx, uyy, and uzz. That brings us to the second moment equations.
For the second moment equation of motion, we will perform a similar manipulation by starting from the relaxation-time approximation (RTA) Boltzmann equation We then encounter a rank three tensor which is defined as where N dof is the number of degrees of freedom. One obtains the following equation of motion from the second moment of the RTA Boltzmann equation [96] ∂ µ I µνλ = 1 τ eq (u µ I µνλ eq − u µ I µνλ ) .
Note that I µνλ is symmetric with respect to interchanges of µ, ν, and λ and traceless in any pair of indices (massless particles/conformal invariance).
In an isotropic system, one finds I x = I y = I z = I 0 with Using the canonical aHydro form one finds with S u (ξ) = 3 + 2ξ (1 + ξ) 3/2 , which satisfy 2S T + S L = S u .
The i = {x, y, z} equations result from with the co-moving derivative D = u µ ∂ µ , the expansion scalar θ = ∂ µ u µ , and θ i ≡ −u µ D i X µ i . For the case of 0+1d Bjorken expansion one has D = ∂ τ , θ = ∂ µ u µ = 1/τ , θ x = θ y = 0 and Based on this, one has ∂ τ I y + 1 τ I y = 1 τ eq (I eq − I y ) , The first two equations (xx and yy projections) both give The third equation (zz projection) gives Taking the zz projection minus one-third of the sum of the xx, yy, and zz projections gives We note that all of the free streaming contributions vanish. Solving forξ using (37) we As mentioned previously, in the limit τ → τ 0 , one has D = 1 and hence the second term on the right-hand-side of (38) will diverge at τ = τ 0 unless either ξ = 0 or ξ = −1. The latter condition makes the entire right hand side vanish and hence does not allow for dynamical evolution of ξ. For this reason we will use lim τ →τ 0 ξ(τ ) = 0.

Cross check (D = 0)
As cross-check on our results, one can recompute the second-moment equation with D = 0 to see if it agrees with results available in the literature. In this case on finds that the zz projection gives and the xx and yy projections both give (log S T ) ξ + 5∂ τ log Λ + 1 τ = 1 τ eq where, in both cases, we used T = R 1/4 (ξ)Λ.
Finally, with D = 0, taking the zz projection minus one-third of the sum of the xx, yy, and zz projections gives One can verify explicitly that Eq. (37) reduces to this in the limit D → 0.

ANISOTROPIC ATTRACTOR
In this section we present some representative numerical solutions using different initial conditions along with the attractor solution to which they flow. For this purpose, we solve the first and second differential equations corresponding to Eq. (24) and (37) for the evolution of ξ(τ ) and Λ(τ ). However to evolve these equations we need to know the damping function.
Herein, we solve the integral equation by using an iterative method. In the first iteration, we assume that the temperature evolution contained within the integral defining D(τ, τ 0 ) is given by ideal hydrodynamics, i.e. T guess (τ ) = T 0 (τ 0 /τ ) 1/3 . We then solve the dynamical equations (24) and (37). From this we obtain the approximate dependence of the effective temperature T on proper time using Eq. (26). The resulting effective temperature T (τ ) is then used to load the damping function for the next iteration. We repeat this process until the effective temperature and longitudinal pressure converge to a part in 10 8 . In practice, this can be achieved with only five iterations. Once converged, the solutions for ξ(τ ) and Λ(τ ) can be used to compute the full distribution function using Eq. (1) and all moments of the distribution function using Eq. (16). 1 In Fig. 1 we present a contour plot of the one-particle distribution function at the proper time at which the contribution from the free streaming part and equilibrating part contribute equally. 2 Generically, the exact solution for the one-particle distribution function contains two independent components [88][89][90][91]. The first component is an anisotropic part which has been squeezed in the longitudinal direction and is exponentially damped at late times. This contribution represents the subset particles that never had any interaction at all. Statistically, there is always such a population of particles. As a function of time, 1 One can substantially reduce the number of iterations required by initializing instead with the canonical aHydro evolution equations. 2 This occurs when D(τ, τ 0 ) = 1/2. FIG. 1. Visualization of the one-particle distribution function at a given moment in proper time.
A bimodal structure can be seen, with the two contributions corresponding to a highly squeezed free-streaming component (inner ellipsoid) and a less anisotropic equilibrating contribution (outer ellipsoid).
this contribution becomes compressed along the longitudinal direction in momentum space resulting in P FS L → 0 as the system evolves. This contribution comes from the first term in the exact solution Eq. (1) which corresponds to the free streaming contribution. Note that, because of the damping function D(τ, τ 0 ) in the first term in Eq. (1), the amplitude of this very narrow ridge will decrease in time exponentially. The second visible component in Fig. 1 is an isotropizing part which dominates at late times. This contribution comes from the second term in Eq. (1).
In Fig. 2, we present the evolution of the scaled moments of the distribution function as a function the scaled time and we compare to the exact RTA solution (black solid lines) obtained in Refs. [88][89][90].
Results from the new aHydro and old aHydro ansatze are shown as red dashed and blue dot-dashed lines, respectively. In all cases shown, the new aHydro ansatz provides a better ansatz is able to reproduce the dynamics of low-order moments much better than the old ansatz. This is particularly striking for moments with m = 0 for which we see that the new aHydro ansatz is very close to the exact results for all n shown. 3 We note, however, for higher moments, e.g. M 33 , we see that the new aHydro ansatz interpolates between the exact solution at early times and the old aHydro result at lates times. As a result, one sees larger deviations from the exact solution in these moments.
In order to provide more quantitative comparison of the two methods, in Fig. 3 we 3 We have checked that this holds true for larger n than shown in Fig. 2. One finds that the new aHydro ansatz fails to accurately describe the evolution of the scaled moments with m = 3 which are dominated by isotropizing contribution at late times. As a consequence, the new aHydro does not provide reliable approximations for these moments and the problem becomes more severe as one increases for m > 1. Turning to the old aHydro ansatz, for m = 1, one sees that, although the old aHydro ansatz does a reasonable job in describing the m = 1 moments, as n and m are increased or decreased, the results become significantly worse. Note that, even given the caveats mentioned above, comparing the old and new aHydro ansatze, we see that the new approach dramatically improves agreement with the exact RTA attractor and, in particular, can be used to fix the problem encountered with moments with m = 0.

IV. CONCLUSIONS
In this study, our goal was to find an improved set of anisotropic hydrodynamic evolution equations that can more faithfully describe the non-equilibrium dynamics of the quark-gluon plasma created in relativistic heavy ion collisions at RHIC and LHC. We introduced a new version of anisotropic hydrodynamics that includes separate free-streaming and equilibrating contributions which allows for a better description of exact solutions to the Boltzmann equation available in the literature. We computed explicit expressions for the first and second moments of the one-particle distribution function in the new aHydro approach and used these to obtain the new 0+1d conformal equations of motion given by Eqs. (24) and (37). We presented comparisons of the numerical solution of the conformal 0+1d equations of motion for both the old and new aHydro schemes with the exact RTA solution. Our results demonstrated that the new aHydro form allows one to have a bimodal distribution function similar to what is seen in the exact RTA solution for the one-particle distribution function. We then computed the evolution of the scaled moments as a function the scaled time, w, and demonstrated that the new aHydro ansatz provides a better approximation to the exact solution than the original aHydro ansatz. Finally, we determined the nonequilibrium attractor associated with the new aHydro scheme and demonstrated that it provides much better agreement with the exact RTA attractor than the original aHydro scheme, in particular for moments with m = 0. In the future, it would be interesting to apply the ansatz obtained here to full 3+1d anisotropic hydrodynamics, including temperaturedependent masses for the particles similar to 'canonical' quasiparticle aHydro [85,97].