(3+1)-dimensional anisotropic fluid dynamics with a lattice QCD equation of state

Anisotropic hydrodynamics improves upon standard dissipative fluid dynamics by treating certain large dissipative corrections non-perturbatively. Relativistic heavy-ion collisions feature two such large dissipative effects: (i) Strongly anisotropic expansion generates a large shear stress component which manifests itself in very different longitudinal and transverse pressures, especially at early times. (ii) Critical fluctuations near the quark-hadron phase transition lead to a large bulk viscous pressure on the conversion surface between hydrodynamics and a microscopic hadronic cascade description of the final collision stage. We present a new dissipative hydrodynamic formulation for non-conformal fluids where both of these effects are treated nonperturbatively. The evolution equations are derived from the Boltzmann equation in the 14-moment approximation, using an expansion around an anisotropic leading-order distribution function with two momentum-space deformation parameters, accounting for the longitudinal and transverse pressures. To obtain their evolution we impose generalized Landau matching conditions for the longitudinal and transverse pressures. We describe an approximate anisotropic equation of state that relates the anisotropy parameters with the macroscopic pressures. Residual shear stresses are smaller and are treated perturbatively, as in standard second-order dissipative fluid dynamics. The resulting optimized viscous anisotropic hydrodynamic evolution equations are derived in 3+1 dimensions and tested in a (0+1)-dimensional Bjorken expansion, using a state-of-the-art lattice equation of state. Comparisons with other viscous hydrodynamical frameworks are presented.


I. INTRODUCTION
Dissipative relativistic fluid dynamics has become the workhorse for simulations of the dynamical evolution of relativistic heavy-ion collisions [1][2][3][4][5][6][7][8]. When supplemented with realistic fluctuating initial conditions, a preequilibrium evolution module that evolves these initial conditions into starting values for the hydrodynamic evolution, and a hadronic rescattering afterburner that describes the late microscopic kinetic evolution of the collision fireball during its dilute decoupling stage, the approach has yielded impressive quantitative precision in its description of a broad set of soft hadronic observables (i.e. distributions of hadrons with momenta below about 1−2.5 GeV/c) obtained from heavy-ion collision experiments at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) [9][10][11][12], and it has demonstrated convincing predictive power when extending the calculations into new domains of collision energy [13][14][15][16][17] or for new collision systems [18][19][20][21][22][23]. Surprisingly, the phenomenological success of dissipative fluid dynamics has so far continued unabated in the description of p+Au, d+Au and 3 He+Au collisions at RHIC and p+p collisions at the LHC [20,23,24], i.e. for "small" collision systems in which the hydrodynamic model had been widely expected to break down. This finding has generated much recent work addressing two obvious questions arising from these observations: (1) What exactly are the formal criteria that ensure the applicability of relativistic dissipative fluid dynamics to small physical systems undergoing rapid collective expansion and control its eventual break-down? How far away from local thermal equilibrium can a system be and still evolve hydrodynamically? (2) Are there alternate mechanisms at work that can mimic the phenomenological signals of hydrodynamic collective flow, especially in small collision systems, without requiring strong final-state interactions among the constituents of the fireball created in the collision that lead to some degree of approximate local thermalization?
Generically, hydrodynamics is an effective macroscopic theory for the late-time, long-distance evolution of sufficiently equilibrated multiparticle systems. It is typically thought of as a gradient expansion around ideal fluid dynamics. The latter describes locally perfectly thermalized fluids in which any deviation from local thermal equilibrium is immediately erased by strong final-state interactions among the microscopic constituents, i.e. systems whose microscopic relaxation time is effectively zero. Relativistic heavy-ion collisions challenge the validity of such an expansion through extremely large density gradients in the initial state, which lead to explosive collective expansion driving the system away from local thermal equilibrium. Even worse, the ultra-relativistic collision kinematics, combined with the quantum mechanics of the initial particle production process that generates the fireball medium from the energy lost by the colliding nuclei in the collision process [25], imprints on the system a very strong initial expansion along the "longitudinal" beam direction, with an approximately boost-invariant longitudinal expansion velocity profile ("Bjorken flow" [26]), while any collective expansion in the directions transverse to the beam is initially small and only builds up later in response to transverse pressure gradients. In realistic fluids with non-zero microscopic mean free paths, the resulting large anisotropy in the collective expansion rate causes large anisotropies in the local rest frame (LRF) momentum distributions of the microscopic constituents. In such a situation an expansion around a locally isotropic thermal equilibrium distribution cannot be expected to converge well.
In addition to rapid and strongly anisotropic expansion (which is most problematic during the earliest stage of a heavy-ion collision) another large dissipative effect occurs towards the end of the evolution when the fireball matter undergoes a phase transition from a quark-gluon plasma (QGP) to a hadron resonance gas (HRG). Critical dynamics near the phase transition causes the bulk viscosity to become large and peak near the (pseudo-)critical temperature T c [27][28][29][30][31][32][33]. The resulting large bulk viscous pressure is associated with a strong deviation of the LRF momentum distribution from thermal equilibrium. Again, this provides a challenge for any expansion around a local equilibrium distribution function. Since (due to the screening of color interactions by color confinement) the constituents' mean free path in the hadron resonance gas is much larger than in the quark-gluon plasma, the Knudsen number (defined as the product of the mean free time and the scalar expansion rate) increases suddenly as the QGP turns into hadrons, to the extent that the subsequent evolution can no longer be reliably described by hydrodynamics [34]. One must therefore switch to a microscopic kinetic description basically as soon as the hadronization process is complete. At this point the critically enhanced bulk viscous pressure is still large because critical slowing down [34][35][36] prohibits it from relaxing quickly to the much smaller values expected away from the phase transition. Its effect on the hadron distribution functions in the HRG can therefore not be treated effectively as a perturbation around local thermal equilibrium and should be accounted for non-perturbatively.
In this work we develop an improved version of anisotropic hydrodynamics that accounts for large shear viscous effects caused by a strong longitudinal-transverse anisotropy of the expansion rate and for large viscous corrections caused by critical dynamics near the quarkhadron phase transition non-perturbatively. The formalism is constructed for full (3+1)-dimensional evolution and tested numerically for (0+1)-dimensional boostinvariant expansion along the beam direction (Bjorken flow [26]). Numerical results for full (3+1)-dimensional evolution of heavy-ion collisions with realistic fluctuating initial conditions [6,15,37,38] will be presented in a future publication. To derive the anisotropic hydro-dynamic evolution equations we start from an underlying kinetic theory, the relativistic Boltzmann equation in relaxation time approximation (RTA BE). While such a classical kinetic approach is known to only work for dilute and weakly coupled gases [39], it allows to derive the structure of the macroscopic hydrodynamic equations through a systematic moment expansion [40][41][42][43]. As an effective theory, the structure of these equations holds equally well for strongly and weakly coupled systems (i.e. it depends only on the separation of microscopic and macroscopic length scales), as long as one changes the material properties of the fluid (i.e. its equation of state, transport coefficients, relaxation times, etc.) that enter as input into the hydrodynamic description to the actual situation of interest.
Our approach starts from the general treatment described in [41,42,[44][45][46][47], expanding the Boltzmann equation around an anisotropic local rest frame distribution function f a which in our case is deformed around the local equilibrium distribution function by two parameters to account for one large shear stress component and a large bulk viscous pressure as described above. Our main innovation is that, following recent insights reported in [43,48], the evolution of these deformation parameters is optimized by determining them through generalized dynamical Landau matching conditions, similar to those fixing the evolution of the temperature and chemical potential. This guarantees that the leading order anisotropic distribution f a (around which the full distribution function is expanded in moments) fully accounts not only for the energy and conserved charge density, but also for the longitudinal and transverse pressures (or, equivalently, the longitudinal-transverse pressure anisotropy (which is the largest shear stress component) and the bulk viscous pressure, as described above). It also significantly simplifies the structure of the relaxation equations for the residual dissipative flows. Writing the full distribution function f as f = f a + δf , the deviation δf describes the (smaller) residual shear stress component and the charge diffusion current. In the present work we will mostly ignore conserved charges and will hence set the charge chemical potential and charge diffusion effects to zero, leaving a complete treatment to a follow-up paper.
Initially, the resulting hydrodynamic framework is formulated in terms of evolution equations for the parameters characterizing the leading-order distribution f a (temperature, chemical potential, momentumdeformation parameters), similar to the traditional approach reported in [46,47,49,50] which makes explicit reference to kinetic theory and thus requires that such a distribution function exists and is well-defined. We subsequently break this connection to an underlying kinetic theory by developing a technique that allows to express these microscopic parameters in terms of macroscopic hydrodynamic quantities so that the formalism can be completely formulated and solved as a macroscopic theory. This procedure introduces the concept of an "anisotropic equation of state" (aEOS) which we discuss at some length. We here use a weakly-interacting quasiparticle model with a temperature-dependent particle mass [51,52] to calculate this aEOS and also the required transport coefficients, keeping mind that in later applications to heavy-ion collisions these ingredients should be computed or modeled for QCD, or considered as phenomenological parameters to be determined from the experimental data.
With this work we open the door to answering the question to what extent complete second-order anisotropic fluid dynamics (which we call "viscous anisotropic hydrodynamics" or vaHydro [41]) extends the range of validity of dissipative fluid dynamics towards smaller collision systems and earlier switching times between the pre-equilibrium and hydrodynamic stages. Having a non-equilibrium hydrodynamic approach that is optimized to the particular challenges posed by ultrarelativistic collisions between nuclei as small as protons is a necessary step towards developing a quantitatively predictive dynamical framework that can set benchmarks for comparison with experimental data and with other, non-hydrodynamic explanations of the latter.
Before starting the technical part of our discussion we introduce our notation. Throughout this work we use natural units = c = k B = 1. The metric signature is taken to be "mostly minus" (+, −, −, −). The local rest frame (LRF) is defined as the Landau frame: the fluid velocity u µ is the normalized time-like eigenvector of the energy-momentum tensor, T µν u ν = Eu µ , where E is the energy density in the local rest frame. It satisfies u µ = T µν u ν u µ T µν u ν and the normalization condition u µ u µ = 1. Unless otherwise indicated we ignore conserved charges and their associated chemical potentials.
The paper is structured as follows: In Sec. II we briefly review the general structure of anisotropic hydrodynamics and its evolution equations for the energy density and longitudinal and transverse acceleration. Relaxation equations for the anisotropic dissipative flows are derived in Sec. III, starting from the relativistic Boltzmann-Vlasov equation in relaxation time approximation and implementing generalized Landau matching conditions for the evolution of the microscopic parameters characterizing the anisotropic distribution function. In Sec. IV we show how to integrate a realistic lattice QCD equation of state into the anisotropic hydrodynamic framework and reformulate its evolution equations in purely macroscopic form, i.e. without any reference to parameters that are only defined in a kinetic theory model for the microscopic dynamics. In that Section we also discuss the resolution of several technical issues arising in the process of solving the anisotropic hydrodynamic evolution equations numerically. In Sec. V we illustrate the performance of our anisotropic hydrodynamic framework in comparison with standard dissipative fluid dynamics for a simple system undergoing Bjorken flow. A summary of our findings and a brief outlook are presented in Sec. VI. Several appendices supply additional technical ingredients, including (in Appendix E) the derivation of the evolution equations for standard viscous fluid dynamics needed for the comparison shown in Sec. V.

A. Ideal fluid decomposition
In relativistic hydrodynamics, the energy-momentum tensor of a perfect fluid is best decomposed in the basis u µ u ν and ∆ µν = g µν −u µ u ν where u µ is the fluid fourvelocity (i.e. the four-velocity of the local rest frame (LRF) relative to the global frame): From here on we will mostly suppress the space-time (x) dependence for notational simplicity. The LRF energy density E and thermal equilibrium pressure P eq are recovered by projecting the energy-momentum tensor onto the temporal and spatial directions: E = u µ u ν T µν and P eq = − 1 3 ∆ µν T µν . Since the thermal pressure is isotropic none of the spatial directions in the LRF are special, and there is no advantage in decomposing the spatial projector ∆ µν further.

B. Viscous fluid decomposition
Standard dissipative fluid dynamics is formulated by decomposing T µν in the same basis, but adding dissipative corrections accounting for the bulk viscous pressure Π and the shear stress tensor π µν : This assumes that the LRF is the Landau frame, i.e. that in the LRF there is no net momentum flow. Eq. (2) implicitly assumes that in the LRF all viscous corrections are of similar order of magnitude and small relative to the equilibrium energy density and pressure, such that there still is no advantage in decomposing the locally spatial projector further into individual spatial directions.

C. Anisotropic viscous fluid decomposition
The general arguments presented in the introduction imply that in relativistic heavy-ion collisions the shear stress π µν is highly anisotropic, and the difference P L −P ⊥ between the longitudinal and transverse pressures can be quite large, due to strongly different longitudinal and transverse expansion rates. This suggests that anisotropic hydrodynamics is best formulated by using a more detailed decomposition of the energy-momentum tensor in which the spatial projector is further decomposed as ∆ µν = Ξ µν − z µ z ν , where z µ points along the beam direction and Ξ µν ≡ g µν − u µ u ν + z µ z ν projects onto the spatially transverse directions in the LRF [42]: The round parentheses around pairs of Lorentz indices indicate symmetrization: W (µ ⊥z z ν) ≡ 1 2 W µ ⊥z z ν +W ν ⊥z z µ . The anisotropic decomposition (3) clearly separates the pressures P L and P ⊥ from the other dissipative components. Given an arbitrary energy momentum tensor T µν , the anisotropic hydrodynamic quantities appearing in this decomposition can be obtained by the following projections: In the last line we introduced the symmetric traceless transverse projection tensor Ξ µν The corresponding transverse shear stress tensor π µν ⊥ describes two shear stress degrees of freedom that account for momentum diffusion currents along the transverse directions. It is traceless and orthogonal to both the fluid velocity and the direction of the pressure anisotropy: Another two shear stress degrees of freedom are encoded in the longitudinal-momentum diffusion current W µ ⊥z which is orthogonal to both u µ and z µ : The remaining fifth (and largest) shear stress component is given by P L −P ⊥ . Altogether, the 5 independent components of the standard shear stress tensor in Eq. (2) are related to those in the anisotropic decomposition (3) by while the single bulk viscous pressure degree of freedom Π in Eq. (2) is related to the thermal, longitudinal and transverse pressures in Eqs. (2,3) by Here P eq is not an independent degree of freedom but related to the energy density E by the equation of state (EOS) of the fluid, P eq (E).

D. Hydrodynamic evolution equations
Four of the ten evolution equations that control the dynamics of the energy-momentum tensor are obtained from the conservation laws for energy and momentum Projecting with u ν on the temporal direction in the LRF provides an evolution equation for the LRF energy density: Here and below a dot over or a D in front of a quantity denotes the comoving time derivative, e.g. DE ≡Ė ≡ u µ ∂ µ E. z µ D z u µ is the scalar longitudinal expansion rate, and θ ⊥ = ∇ ⊥ ·u is the scalar transverse expansion rate. The longitudinal derivative and transverse gradient in the LRF are written as D z = −z µ ∂ µ and ∇ ⊥µ = Ξ ν µ ∂ ν , respectively. σ ⊥,µν = Ξ αβ µν ∂ α u β is the transverse velocityshear tensor.

III. DISSIPATIVE RELAXATION EQUATIONS
To close the system of equations, we need six additional relaxation equations for P L , P ⊥ , W µ ⊥z and π µν ⊥ . Their dynamics is not controlled by macroscopic conservation laws but by microscopic interactions among the fluid's constituents. As discussed in Sec. I, we will here derive them by assuming a weakly-coupled dilute fluid whose microscopic physics can be described by the relativistic Boltzmann-Vlasov equation for a single particle species with a medium-dependent mass: Here f (x, p) is the single particle distribution function, C[f ] is the collision kernel, m(x) is the medium-dependent effective mass, and ∂ (p) µ is the momentum derivative.

A. Anisotropic distribution function
For anisotropic hydrodynamics we split the distribution function into a momentum-anisotropic leading-order contribution f a (x, p) and a small residual correction δf : For the leading order distribution we take the generalized Romatschke-Strickland form [53][54][55]: Here f eq (z) = 1/(e z + Θ) is the equilibrium distribution, with Θ = 1, 0, −1 for Fermi-Dirac, Boltzmann, and Bose-Einstein statistics, respectively. Λ(x) is an effective temperature andμ(x) an effective chemical potential. In this work we consider a system without conserved charges and assume, for simplicity, that particle number changing microscopic processes in the collision term C[f ] are so fast that the effective chemical potential relaxes to zero faster than any other microscopic time scale. 2 In this work the leading-order momentum anisotropy is encoded in the ellipsoidal tensor It contains two space-time dependent anisotropy parameters ξ L and ξ ⊥ . With Ξ µν p µ p ν = −p 2 ⊥,LRF (i.e. the square of the transverse momentum p ⊥,LRF in the LRF), Ω µν p µ p ν can be rewritten as The difference ξ L −ξ ⊥ can be attributed to a manifestation of shear stress (resulting in a difference between the the longitudinal and transverse pressures) while the sum ξ L +ξ ⊥ encodes a bulk viscous pressure [56,57]. Introducing the notation α L,⊥ (x) = 1 + ξ L,⊥ (x) −1/2 , Eq. (15) can be written in LRF momentum components more conveniently as To make the decomposition (14) unique one must specify the three parameters Λ(x) and α L,⊥ (x). We proceed as follows: The physical temperature T of the system is defined by the LRF energy density via the thermodynamic relation E(x) ≡ E T (x) . To relate the effective temperature Λ to T we impose the generalized Landau matching condition 3 which states that Λ(α) must be chosen such that the residual deviation δf does not contribute to the energy density. This fixes Λ(α) as a function of T ; the two agree in the limit α → 1 when the anisotropic leading-order distribution f a reduces to a locally isotropic equilibrium distribution f eq (u·p/T ). The momentum deformation parameters α L,⊥ (x) are fixed by similar generalized Landau matching conditions for the transverse and longitudinal pressures Here π zz LRF is the LRF value of the longitudinal diagonal element of the shear stress tensor π µν in the decomposition (2). Note that both the bulk viscous pressure Π and the shear stress component π zz LRF are here assumed to be "large" such that they must be accounted for already at leading order, by adjusting the parameters α L,⊥ (x) accordingly. This is done by demanding By imposing these conditions, α L,⊥ (x) are adjusted such that the longitudinal and transverse pressures P L and P ⊥ are everywhere fully accounted for by the leading-order distribution f a , with zero residual contributions from δf . This is an application of the anisotropic matching scheme proposed by Tinti in [55] and a generalization of the P Lmatching scheme proposed and studied in [43,48] to both P L and P ⊥ (or, equivalently, to the pressure anisotropy P L −P ⊥ ∼ π zz LRF and the bulk viscous pressure Π). In this matching scheme, the δf correction generates only the residual dissipative flows described by W µ ⊥z and π µν ⊥ , which break the cylindrical symmetry of the distribution function in the LRF and account for the remaining four smaller components of the shear stress tensor π µν in Eq. (7).
With the matching conditions (19), (20), and (22a,b) we have the following kinetic theory expressions for the particle number and energy densities as well as for the longitudinal and transverse pressures: 4 n (k) = u · p fa = I 1000 , (23a) The "anisotropic integrals" I nrqs over the leading-order distribution function f a that appear in these equations are defined in Eq. (A1).

B. Relaxation equations I
The relaxation equations for the dissipative flows are obtained by expressing the latter as moments of the distribution function and using the Boltzmann equation to describe its evolution, using the decomposition (14) and treating δf as a small perturbation. We start froṁ where we defined the compact notations [42] for the spatially transverse (in the LRF) components of a vector a µ or its LRF time derivativeȧ µ and the spatially transverse and traceless part of a tensor t µν or its LRF time derivativeṫ µν , as well as for the Lorentz-invariant momentum space integral, with g being a degeneracy factor counting the number of quantum states allowed for a particle with on-shell momentum p µ , and Θ(p 0 ) denoting the Heaviside step function.
After moving the time derivative D on the r.h.s. under the integral until it hits the distribution function f a or δf , we use the decomposition f = f a + δf together with to rewrite the Boltzmann-Vlasov equation (13) in the formḟ Closing this equation requires an approximation for δf . We here use the 14-moment approximation.

C. The 14-moment approximation
The 14-moment approximation derives its name from approximating δf in terms of its 14 momentum moments with p µ and p µ p ν (where the moment with the linear combination p µ p µ g µν = m 2 as weight function is equivalent with the moment taken with weight 1) [58,59]. In our case the choice of the Landau frame, together with the generalized matching conditions (19), (20), and (22a,b) and the absence of diffusion currents related to conserved charges, eliminate ten of these moments, leaving only four independent moments to construct δf . These need to be matched to the residual dissipative flows W µ ⊥z and π µν ⊥ , which each have two degrees of freedom. The 14-moment approximation for δf can thus be written as [42] δf and decoupling the resulting set of linear equations: The anisotropic integrals J nrqs appearing in these expressions are defined in Eq. (A2). As expected, the coefficients are directly proportional to the residual dissipative flows:

D. Relaxation equations II
Substituting the 14-moment approximation (32) for δf into Eqs. (24) and (28), simplifying some of the resulting terms by integrating by parts, and enforcing the generalized matching conditions, some algebra yields the following dissipative relaxation equations: 5 is the average pressure as given by kinetic theory, and ω µν The structure of Eqs. (33)- (36) is simpler than that of the corresponding equations derived in [42], not only by the absence of terms coupling to the conserved charge and diffusion currents (which only reflects the simplifying assumptions made here), but also as a result of imposing the generalized Landau matching conditions (19), (20), and (22a,b) which optimizes the evolution of the anisotropy parameters in f a and thus removes additional terms needed in [42] to correct their evolution if not chosen optimally in the first place.
The transport coefficients appearing on the right hand sides of Eqs. (33)-(36) are labeled following as much as possible the convention established in Ref. [42]. Except for the relaxation times they are given in Appendix C. 6 5 In deriving these equations one encounters terms involving the comoving time derivativeṁ of the temperature-dependent quasiparticle mass that arise from the second term on the l.h.s. of the Boltzmann-Vlasov equation (13). We eliminate them by using the chain ruleṁ = (dm/dT ) (dT /dE)Ė where we take dm/dT as external input from the quasiparticle model discussed below, evaluate the derivative dT /dE = (dE/dT ) −1 = c 2 s T /(E + Peq) from the lattice QCD EOS, and use Eq. (10) forĖ. Equations (33)-(36) (together with the transport coefficients listed in Appendix A) are found after combining the terms on the r.h.s. of Eq. (10) with other terms involving the same dissipative forces. In doing so we neglect contributions of is the residual inverse Reynolds number associated with the residual components W µ ⊥z and π µν ⊥ . 6 Please note the superscripts (k) on the transport coefficients ap-Generically they involve the "anisotropic thermodynamic integrals" over the anisotropic distribution function f a given in Appendix A. Their validity, as well as the validity of the specific relations between some of these transport coefficients listed in Appendix C, depends on the applicability of relativistic kinetic theory of a gas of weaklyinteracting quasiparticles as the underlying microscopic theory, which is not guaranteed for quark-gluon plasma. Their generalization to a realistic microscopic theory of QCD medium dynamics requires much additional work. We will here use the expressions given in the Appendices as order-of-magnitude estimates and placeholders for future more realistic sets of transport coefficients.
Equations (33)-(36) also involve two relaxation times, τ π and τ Π . τ Π controls the relaxation of the kinetic bulk viscous pressureP (k) −P (k) eq (see Eq. (8)) whereas the shear relaxation time τ π drives the relaxation of both the large shear stress component P ⊥ and the smaller ones described by W µ ⊥z and π µν ⊥ . That all shear stress components have the same relaxation time even if some of them become large is a model assumption that may be corrected in future improved calculations of the transport coefficients for strongly anisotropically expanding QGP.
Formally, the relaxation times arise from a linearization of the collision term around the local equilibrium pearing on the right hand sides of Eqs. (33) and (34). They reflect the fact that these control the evolution of the kinetic part of the longitudinal and transverse pressures. For the quasiparticle model introduced in the next Section an additional mean field enters which modifies these pressures and transport coefficients. The modified expressions will be denoted without the superscript.
distribution f eq (with temperature computed from the energy density): Literal use of this Relaxation Time Approximation (RTA) [60] gives τ π = τ Π = τ r . However, strong coupling in the quark-gluon plasma in the temperature regime just above the quark-hadron phase transition, as well as critical behavior near that phase transition, lead to very different temperature dependences of the bulk and shear viscosities and their associated relaxation times in QCD, especially around T c [28][29][30][31]61]. In particular, the bulk relaxation time τ Π is expected to be affected by "critical slowing down" [29,33,35], i.e. it should exhibit a strong peak near T c . Since large bulk viscous effects near T c are one of the main motivations for our work here, we feel compelled to account for them by introducing two different relaxation times τ π and τ Π , and tying them to phenomenologically parametrized shear and bulk viscosities η and ζ by postulating the standard kinetic theory relations [40] The (temperature dependent) isotropic thermodynamic integrals β π and β Π appearing in these relations are given further below in Eq. (82). The viscosities η and ζ are transport parameters that occur in standard "isotropic" dissipative fluid dynamics -they appear here through the relaxation times τ π and τ Π . When comparing anisotropic with standard dissipative fluid dynamics further below in Sec. V we will do so by using the same functions τ π and τ Π in both approaches.

IV. ANISOTROPIC EQUATION OF STATE
While the relaxation equations (33)-(36) were derived from the Boltzmann equation, the equations remain structurally unchanged for strongly coupled fluids. They are purely macroscopic, i.e. all terms on the r.h.s. have the form of some macroscopic driving force (proportional to the Knudsen or inverse Reynolds numbers or products thereof) multiplied by some transport coefficient. The kinetic origin of these equations is hidden in these transport coefficients. Applying the equations to strongly coupled fluids requires only that these transport coefficients, along with the equation of state relating the energy density and equilibrium pressure, are swapped out accordingly.
For the time being most of the transport coefficients of hot and dense QCD matter are still essentially unknown. While the shear and bulk viscosities will be taken as parameters whose functional forms are modeled phenomenologically and whose overall magnitudes are to be fitted to experimental observables, the remaining transport coefficients will be approximated using kinetic theory, for reasons of consistency with our derivation of the evolution equations. Their evaluation requires microscopic kinetic inputs, namely the parameters (Λ, α ⊥ , α L ) characterizing the anisotropic distribution function f a , the particle mass m, and also the temperature T . However, for the QGP equation of state, which is very precisely known from lattice QCD calculations [62,63], we want to use first-principles theoretical input.
In this Section we discuss how to consistently incorporate such direct information from QCD into a hydrodynamic framework that was originally derived from a kinetic theory with a very different EOS. We will call this procedure "integrating the lattice QCD EOS with some kinetic framework". We introduce a parametric model for an anisotropic equation of state that allows the anisotropic hydrodynamic equations, including the dissipative relaxation equations for the longitudinal and transverse pressures and the remaining shear stress components, to be solved on a purely macroscopic level. This differs from earlier implementations of the framework which relied on the solution of dynamical evolution equations of the microscopic kinetic parameters (Λ, α ⊥ , α L , m) [43,44,46,[49][50][51] (which, for the case of QCD, are not really well-defined). However, since we will need these microscopic parameters for the calculation of those transport coefficients that we compute from kinetic theory (a temporary necessity that will disappear as soon as ways have been found to calculate these transport coefficients directly from QCD), we determine them from the macroscopic hydrodynamic quantities at the end of each time evolution step, using our parametric model for the anisotropic EOS. 7 To construct this parametric model we follow Refs. [51,52,64] and parametrize the response of the pressure anisotropy and the bulk viscous pressure to anisotropic expansion using a quasiparticle EOS. The quasiparticles have a temperature-dependent mass that is chosen such that a weakly-interacting gas of these particles accurately mimics the QCD EOS. The transport coefficients are then worked out in this kinetic theory. 8 It is well known [65][66][67] that for thermodynamic consistency such an approach requires the introduction of a mean field B whose temperature dependence in equilibrium generates the temperature dependence of the quasiparticle's effective mass. It also receives additional dissipative corrections out of equilibrium [52]. 7 Note that the parametric model is not used for the equilibrium EOS Peq(E) (for which we take state-of-the-art lattice QCD results) but only to parametrize the dissipative deviations of the longitudinal and transverse pressures from Peq(E), as well as for the calculation of the remaining transport coefficients. 8 Note that an accurate description of the equation of state does not imply by any means that the kinetic theory also predicts the correct transport properties of the medium. In all likelihood it doesn't.

A. Integrating the lattice QCD EOS with a quasiparticle EOS
The key question that needs to be addressed in anisotropic hydrodynamics is how much pressure anisotropy and bulk viscous pressure is generated by a given hydrodynamic expansion rate and its anisotropy. These are the two largest and most important dissipative effects in our approach. The answer to this question depends on the microscopic properties of the medium. For QCD matter this response is presently not known. It is, however, a key ingredient in the hydrodynamic evolution model. In this subsection we model this response by that of a weakly interacting gas of quasiparticles with a medium-dependent mass m(T ). Within this model we can associate (within certain limits) any given deviations of the longitudinal and transverse pressures P L and P ⊥ from the equilibrium pressure P eq (E) with specific values for the microscopic parameters Λ, α, m T ) describing the anisotropic quasiparticle distribution function f a . These values can then be used to compute the kinetic theory values for the transport coefficients. So while the equilibrium pressure is described by the full QCD EOS from lattice QCD, the dissipative deviations of P L and P ⊥ from the equilibrium pressure are interpreted microscopically within a weakly interacting gas of massive Boltzmann particles. As we solve the hydrodynamic equations (10)- (12) together with the dissipative relaxation equations (33)-(36), we interpret the resulting deviations from local equilibrium within the quasiparticle model by writing Here the superscript (q) stands for "quasiparticle model". The zero on the l.h.s. of the first of these equations reflects the Landau matching condition E = E(T ) to the lattice QCD energy density, which also provides us with the temperature T at which the quasiparticle mass m(T ) and equilibrium mean field B eq (T ) (see below) are evaluated. In the quasiparticle model the hydrodynamic quantities on the r.h.s. of (39) consist of kinetic and mean field contributions [64] The kinetic contributions are obtained from Eqs. (23b,c,d): where T = T (E). The mean field B consists of an equilibrium part B eq and a dissipative correction δB: By Landau matching, the total quasiparticle energy density E (q) is fixed to its equilibrium value: where E This establishes a relation between the temperature and the kinetic theory parameters, provided that δB is determined. In the equilibrium limit, the quasiparticle pressure is where P (k) eq (T ) = I 2200 T, 1; m(T ) . The equilibrium terms in (43) and (45) are all functions of temperature.
For simplicity we assume that the quasiparticles have Boltzmann statistics (Θ = 0). To ensure that at asymptotically high temperature the equilibrium pressure and energy density of this Boltzmann gas approach the corresponding values of a quark-gluon gas with 2(N 2 c −1) bosonic and 4N c N f fermionic degrees of freedom, we normalize them by applying to the quasiparticle distribution function a degeneracy factor with N c = 3 colors and N f = 3 flavors, counting u, d, and s quarks only (heavier flavors are exponentially suppressed in the phenomenologically interesting temperature range and are therefore neglected). This degeneracy factor is part of the momentum integration measure p in the definition (26).
The thermal quasiparticle mass m(T ) is chosen such that the equilibrium pressure P (q) eq (E) and energy density E (q) eq (E) of the quasiparticle model agree with their lattice QCD counterparts. Technically this is done by expressing the lattice QCD entropy density S in terms of the corresponding kinetic theory expression S (q) for a gas of quasiparticles with mass m(T ) and Boltzmann statistics [64]: where K n (z) is the modified Bessel function with z = m(T )/T . For thermodynamic consistency the right hand side must satisfy S (q) = dP (q) eq /dT , which is ensured by setting whereẑ = m(T )/T . Plots of the quasiparticle mass-totemperature ratio z and equilibrium mean field B eq as functions of T , using a 2010 lattice QCD EOS obtained by the Wuppertal-Budapest Collaboration [62], can be found in [64]. Here we use the state-of-the-art QCD EOS compiled by the Beam Energy Scan Theory (BEST) Collaboration [68]. The resulting slightly modified temperature dependences of z(T ), dm(T )/dT and B eq (T ) are shown in Fig. 1 as solid lines (together with the earlier results from Ref. [64] shown as dashed lines). Equation (48) determines the mean field in equilibrium. Out of equilibrium it receives a non-equilibrium correction δB [52]. As shown in [52], thermodynamic consistency and energy-momentum conservation can be used to derive from the Boltzmann equation the following general evolution equation for B: By Landau matching the non-equilibrium correction to the quasiparticle energy density E (q) = E (k) +B must vanish, hence where δf = f −f eq . Substituting f = f a +δf and using the relaxation time approximation equation (49) takes the forṁ Note that the expression in the parentheses is the trace of the kinetic contribution to the energy-momentum tensor T µν . Since the non-equilibrium component of the mean field δB = B − B eq contributes to the bulk viscous pressure, we have replaced in Eq. (52) the relaxation time τ r by the bulk relaxation time τ Π . The time derivative of the thermal mass can be expressed in terms of the energy conservation law (10) using the chain rule. Although Eq. (50) shows that B is not an independent quantity, we find it most straightforward to use Eq. (52) to propagate the mean field B dynamically. It does not directly enter the evolution equations for the components of the energy momentum tensor as an independent variable, but is only needed for the model interpretation of the pressure anisotropy and bulk viscous pressure (which are hydrodynamic outputs) in terms of the microscopic parameters (Λ, α) needed for computing the transport coefficients in Appendix C. We use Eqs. (10), (33), (34) and (52) to evolve E, P solve these equations for the anisotropy parameters (Λ, α ⊥ , α L ), and compute the transport coefficients. Of course, the values Λ, α, m associated in this way with P L , P ⊥ and E at any point of the hydrodynamic space-time grid are model dependent, and a different parametrization of the lattice QCD EOS in terms of quasiparticles (for example, as a mixture of different types of quasiparticles with different quantum statistical properties, degeneracy factors and masses) would yield different results. For example, we have tried (and abandoned) an alternate approach where we used a weakly interacting Boltzmann gas of particles with a fixed mass to interpret the pressure anisotropy and bulk viscous pressure in terms of microscopic parameters (Λ, α; m). In that case we were unable to find solutions at early times where the strong longitudinal expansion leads to negative bulk viscous pressures large enough that no valid choice of microscopic parameters can reproduce this in the kinetic model theory. In the quasiparticle model we can partially absorb this with an out-of-equilibrium mean field contribution.
With the anisotropic EOS model (53) we can finally write down the equations of motion for the total pressures P L and P ⊥ , by combining Eqs. (33,34) with Eq. (52): Here we redefined the transport coefficients for the longitudinal and transverse pressures as detailed in Appendix D. This completes our formalism for nonconformal anisotropic hydrodynamics, where the equations of motion (10)-(12), (35), (36), (54) and (55) are purely macroscopic and structurally independent of the underlying microscopic physics, while the transport coefficients are evaluated with our specific quasiparticle kinetic model for the anisotropic equation of state.

B. Reconstructing energy density and fluid velocity
Most numerical codes developed for heavy-ion collisions (including ours) solve the energy-momentum conservation laws (9) on a fixed "Eulerian" space-time grid instead of the LRF projected conservation equations (10)- (12). The reason for this is the existence of powerful flux-corrected evolution algorithms for conservation laws of the type (9) [69][70][71]. To be able to use these algorithms one writes Eqs. (9) in conserved flux form, which for Milne coordinates x µ = (τ, x, y, η) reads The source term J µ on the r.h.s. includes both geometric (Christoffel) terms arising from the curvilinear nature of the coordinate system and the dissipative fluxes [71].
Running the evolution algorithm for one time step thus produces updated values for the first row of the energymomentum tensor and the dissipative flows encoded in P L , P ⊥ , W µ ⊥z and π µν ⊥ , as well as the mean field B. For the evaluation of the EOS and the calculation of the transport coefficients as described in the preceding subsection we need, however, the LRF energy density E, and to compute the source term J µ we need the fluid four-velocity u µ . How to obtain these from the updated output of the evolution code at the end of a time step is described in this subsection.
Next, we observe that the term 2W (τ ⊥z z µ) depends (through the components of z µ ) on u τ and u η . It cannot be subtracted from T τ µ until a relation between u τ and u η is found. To this end let us construct from known quantities the vector K µ = T τ µ − π τ µ ⊥ whose components read [72] In the last equation we used the orthogonality relation z µ W µ ⊥z = 0 to eliminate W η ⊥z . Taking the combination u η 2 K τ − u τ u η K η one obtains where as well asÃ = τ A andB = τ B are all known quantities. One can further use the normalization condition u µ u µ = 1 to rewrite Eq. (60) either as or as whereF = τ F and x = 1−F 2 . With this the components of z µ in (57) and thus 2W (τ ⊥z z µ) are now known. We next define the known vector M µ = K µ −2W (τ ⊥z z µ) , with components From Eq. (64b) one obtains immediately the transverse flow velocity components and magnitude: where M ⊥ = (M x ) 2 +(M y ) 2 . The two remaining unknown variables are u τ and E. By taking the combination u 2 ⊥ M τ − u τ u i M i one finds a relation between u τ and E: This can now be used to express u i and u ⊥ in Eqs. (65) as well as u η in Eq. (60) in terms of E and the other known quantities. With a bit of algebra, and making use of the relationF = τ M η /(M τ +P L ), the normalization condition (u τ ) 2 − u 2 ⊥ − (τ u η ) 2 = 1 then yields the following explicit reconstruction formula for the energy density: Note that, since P L and P ⊥ are evolved directly, the r.h.s. of Eq. (67) is entirely known. This is in contrast with the reconstruction formula for E in viscous hydrodynamics where one must solve some nonlinear equation F(E) = 0 numerically [72,73]. Once Eq. (67) is evaluated, the fluid velocity components u τ , u i , and u η can be determined using Eqs. (66), (65) and (60) consecutively.

C. Reconstructing the anisotropic parameters
Having reconstructed the LRF energy density E, the left hand sides of equations (53) are all known, as is the temperature T (from E(T )) and thus the particle mass m(T ). We can now use the expressions on the right hand sides of these equations to determine the anisotropic parameters (Λ, α ⊥ , α L ). We write equations (53) as We then use Newton's method in three dimensions to find the solution vector X. One starts with an initial guess vector X (0) , which is typically the value of X at the spatial grid point known from the previous time step. This vector is updated with an iteration ∆X, which is given by the matrix equation where J ij = ∂F i /∂X j is the Jacobian matrix. The analytical form of this matrix is One can simplify the computation of the matrix elements in (72) The iteration process (71) is repeated until convergence is achieved. Newton's method is a local root-finder and works well if the initial guess X (0) is sufficiently close to the solution. This may not be the case when the hydrodynamic variables are first initialized or evolving too rapidly. To better handle such situations, we include a line-backtracking algorithm, which takes partial steps of ∆X, to improve global convergence [74]. Each iteration of Newton's method requires the numerical computation of eight one-dimensional integrals, three for F and five for J. Alternatively, one may use Broyden's method, which approximates the Jacobian in terms of F , so only three integrals need to be evaluated. [74] An approximate Jacobian may complicate the linebacktracking search, which ensures a decrease of |F | only if the exact Jacobian is used. For iterations where the full Broyden step does not sufficiently decrease |F | we switch to Newton's method.
In the anisotropic equation of state model (53) there is an ambiguity between the initialization of the mean field B and of the kinetic terms E (k) , P (k) L , and P (k) ⊥ . Standard hydrodynamic initial conditions for the energymomentum tensor only provide E, P L , and P ⊥ on the initialization hypersurface. The initial energy density profile E also yields the initial temperature profile and thus the initial profile for the equilibrium part B eq (T ) of the mean field. Its comoving time derivative on the initialization surface can be obtained by taking the equilibrium limit of Eq. (52): To obtain a guess for the initial non-equilibrium deviation δB we assume that δB evolves on a time scale larger than the bulk viscous relaxation time τ Π . We can then ignore the time derivative of δB on the left hand side of Eq. (52) and obtain from the difference between Eqs. (52) and (74) the "asymptotic" initial condition As before (see footnote 5) m andṁ can be expressed in terms of the energy density E and its LRF time derivative. Having thus specified the initial profile for the mean field B = B eq + δB we can proceed to extract the initial anisotropic parameters and compute the initial values for the transport coefficients. For far-from-equilibrium initial conditions, such as those provided by the IP-Glasma model [37] where P L starts out with a very large negative value P L = −E and, after classical Yang-Mills evolution for a time of the order of the inverse saturation momentum, settles to around zero [75,76], the implied deviation P L −P eq can become so large that, with this initial choice of B, the quasiparticle model cannot accommodate it within the allowed ranges for (Λ, α ⊥ , α L ). Since typically B < 0 at high temperatures (see Fig. 1), the kinetic longitudinal pressure, P (k) L = P L +B, may in this situation be negative. Specifically, the anisotropic parameter initialization is found to fail when P L /P ⊥ 0.08. To overcome this problem, in the case of such extreme initial conditions for P L we simply adjust our initial guess for δB and increase the initial value for B until a solution for (Λ, α ⊥ , α L ) can be found. More meaningful ways of dealing with this shortcoming will be left to future work.

V. BJORKEN FLOW
In this Section we test our anisotropic hydrodynamic formalism by comparing it to standard viscous hydrodynamics for the case of (0+1)-dimensional Bjorken expansion, using the state-of-the-art lattice QCD equation of state referenced in Fig. 1. We begin by simplifying the anisotropic evolution equations (10)-(12), (35)- (36), and (54)-(55) for systems with Bjorken symmetry. In Milne coordinates, the fluid velocity is u µ = (1, 0, 0, 0), the longitudinal and transverse expansion rates are z µ D z u µ = 1/τ and θ ⊥ = 0, and the residual shear stress components W µ ⊥z and π µν ⊥ vanish by symmetry. The component T τ τ trivially reduces to E. As a result, the anisotropic hydrodynamic equations simplify tȯ where in (76d) we useḋ To fix the relaxation times τ π and τ Π we proceed as follows: we use a temperature dependent parametrization for the specific shear viscosity η/S [6], where S = S(E) is the lattice QCD entropy density and T c = 154 MeV is the pseudo-critical temperature. The model parameters (η/S) min = 0.08 and (η/S) slope = 0.85 GeV −1 were extracted from a global Bayesian analysis of RHIC and LHC heavy-ion collision data [6]. 9 Similarly, we use for the specific bulk viscosity ζ/S the parameterization from Ref. [77]: where the function f (x) is given by with A 0 = −13.45, A 1 = 27.55, A 2 = −13.77, C 1 = 0.03, C 2 = 0.001, λ 1 = 0.9, λ 2 = 0.22, λ 3 = 0.9, λ 4 = 0.25, σ 1 = 0.0025, σ 2 = 0.022, σ 3 = 0.025 and σ 4 = 0.13. For the normalization factor we choose (ζ/S) norm = 1.25 [6], and we fix the location of the peak of the specific bulk viscosity by taking T p = T c . 10 Fig. 2 shows the behavior of the specific shear and bulk viscosities as a function of temperature.
The relaxation times are then obtained from the kinetic theory relations (38), rewritten in the form 9 Note that, while some of these parameters where fitted to experimental data [6], this was done with standard viscous hydrodynamics, and slightly different values might be expected when repeating that exercise with anisotropic hydrodynamics. The precise values of the parameters in Eqs. (78)- (80) should therefore not be taken too seriously. 10 Note that this puts the peak of the bulk viscosity at a much lower temperature than assumed in most previous implementations of this parametrization (see e.g. [6,[77][78][79][80]).
using the following quasiparticle versions of the isotropic thermodynamic integrals β π and β Π [52]: Here c 2 s (E) is the squared speed of sound from lattice QCD. The system of ordinary differential equations (76) is solved using Huen's method. After each intermediate and full time step the anisotropic parameters are updated by numerically solving Eq. (68).
These anisotropic hydrodynamic results will be compared with those from second-order viscous hydrodynamics in the 14-moment approximation. The corresponding evolution equations and transport coefficients are derived in Appendix E. For Bjorken flow, the set of independent dynamical variables reduces to the energy density E, the shear stress π = −τ 2 π ηη = 2 3 (P ⊥ −P L ), and the bulk viscous pressure Π = 1 3 (2P ⊥ +P L ) − P eq . Their evolution equations simplify tȯ For the non-equilibrium mean field contribution δB we use the second-order expression [52] δB (2) where θ = ∂ µ u µ = 1/τ is the scalar expansion rate. In Eq. (84), we replaced the relaxation time τ r by τ Π . The relaxation times τ π and τ Π are obtained from Eqs. (81) and (82) while the second-order transport coefficients τ ππ , δ ππ , λ πΠ , δ ΠΠ and λ Ππ are computed from the quasiparticle model in the 14-moment approximation (after expansion around a local equilibrium distribution, see Appendix E). We will also look at how the transport coefficients, including the relaxation times, affect the viscous hydrodynamic results when computed in the small fixed mass approximation z 1 and dm/dT ≈ 0, without a mean field, which is commonly implemented in viscous hydrodynamic simulations. [78,81,82] A. Equilibrium initial conditions In Figure 3 we show the Bjorken evolution of the hydrodynamic variables in anisotropic hydrodynamics, including the total mean field and its non-equilibrium component in the quasiparticle (QP) model used to compute the transport coefficients, and compare it with that in the standard viscous hydrodynamic models. Figure 4 shows the same for the associated Knudsen and inverse Reynolds numbers. In this subsection we impose equilibrium initial conditions with initial temperature T 0 = 0.5 GeV at longitudinal proper time τ 0 = 0.25 fm/c, i.e. all non-equilibrium effects are initially zero. Figure 3a shows that all three models (anisotropic hydrodynamics with QP transport coefficients in solid red lines, stan-dard viscous hydrodynamics with QP transport coefficients in long-dashed blue lines and transport coefficients from a Boltzmann gas in a small fixed mass expansion in short-dashed green lines) produce almost identical evolutions for the energy density. The energy density decreases somewhat more slowly than for a conformal ideal fluid, indicated by the thin black line ∼ τ −4/3 . This is due to the smaller pressure of our EOS (which thus performs less longitudinal work) and to viscous heating. For reference we note that the system passes through the pseudocritical temperature T c = 154 MeV at τ c ∼ 37 fm/c, with a small spread of less than 2 fm/c between the three models.
Panel (c) shows that, if a QP model is used for the transport coefficients, the mean field B also evolves almost identically in anisotropic and standard viscous hydrodynamics. Small differences between anisotropic and standard viscous hydrodynamics with QP transport coefficients are observed in the evolution of the shear stress π (O(2%)) and the pressure ratio P L /P ⊥ (O(10%)): the effective resummation of shear viscous effects in anisotropic hydrodynamics leads to a slight reduction of the shear stress, resulting in a slightly reduced pressure anisotropy. Standard viscous hydrodynamics, with transport coefficients calculated in the small fixed mass expansion (short-dashed green lines), produces somewhat (O(15%)) larger shear stresses and stronger pressure anisotropies.
Given that the pressure anisotropy gets quite large, with P L /P ⊥ decreasing to about 30% at τ ∼ 1 fm/c, the excellent agreement between standard viscous and anisotropic hydrodynamics is somewhat unexpected. It suggests that the widely used standard viscous hydrodynamic approach is quite robust and quantitatively reliable even for large shear stresses. Similar observations were made before in [72] as well as in studies of the Bjorken dynamics of strongly coupled theories where second-order viscous hydrodynamics could be directly compared with an exact numerical solution of the underlying strong-coupling dynamics [83,84].
The largest differences between anisotropic and standard viscous hydrodynamics are seen in the evolution of the bulk viscous pressure Π (Fig. 3d) and the nonequilibrium part of the mean field δB (Fig. 3f). The two panels expose strong correlations between the evolutions of these two quantities. Both are small: (i) The bulk viscous pressure at early times is about 100 times smaller than the shear stress. While the evolution of Π is qualitatively similar (although quantitatively different by more than a factor 2 at early times) for anisotropic and standard viscous hydrodynamics with QP transport coefficients, it exhibits qualitatively different dynamics in standard viscous hydrodynamics with transport coefficients computed from the small fixed mass expansion.
(ii) Compared to the equilibrium mean field, the nonequilibrium part δB is about two orders of magnitude smaller (see panels (c) and (f) of Fig. 3). Here one observes very different trajectories for δB between the evolutions from anisotropic and standard viscous hydrody-namics, although their shapes are qualitatively similar. In addition, panel (f) shows for comparison the "asymptotic approximation" (75) (short-dashed purple curve) which should be compared to the exact numerical solution (red solid line). Obviously, the large expansion rate at early times makes the asymptotic trajectory, which is based on the assumption that δB evolves more slowly than the bulk relaxation rate, a rather crude approximation. Figure 4 shows the Knudsen and inverse Reynolds numbers associated with the shear and bulk viscous stresses. While the Knudsen and inverse Reynolds numbers associated with shear stress dominate the nonequilibrium dynamics at early times, those associated with bulk viscosity are the most relevant at late times when the system passes through the QCD phase transition. 11 In spite of the shear Knudsen number (Fig. 4a) starting out large with a value of around 2.5, the shear inverse Reynolds number (Fig. 4b) never exceeds a value of about 75-85%. This results from the delay caused by the microscopic shear relaxation time which controls the approach of the shear stress π from its zero starting point to its Navier-Stokes value and has at τ 0 = 0.25 the value τ π ≈ 0.8 fm/c. By the time R −1 π reaches its Navier-Stokes limit (shown in Fig. 4c), the shear Knudsen number has already dropped to values well below 1. We reiterate that at the peak value ∼ 3/4 of the shear inverse Reynolds number the differences between anisotropic and standard viscous hydrodynamic evolution are less than 6% as long as both are evaluated with transport coefficients computed from the same underlying kinetic theory. Fig. 4e shows the evolution of the bulk inverse Reynolds number R −1 Π , which peaks due to critical dynamics near the QCD phase transition temperature T c ; the corresponding Navier-Stokes value is shown in Fig. 4f. Because τ Π ∝ ζ (see Eq. (81)), the bulk relaxation rate slows down when the bulk viscosity peaks. This leads to "critical slowing down" of the evolution of the bulk viscous pressure Π, limiting its growth as the system cools down to T c [33,35]. Comparing the solid red and dashed blue curves in Figs. 4e and f we see that Π and thus R −1 Π never reaches much more than about half of its peak Navier-Stokes value, and it also peaks later (around τ ∼ 38 fm/c ≈ τ c , corresponding to T ≈ 0.995 T c ) than the Navier-Stokes limit which reaches its maximum already at τ ∼ 27 fm/c (corresponding to T ≈ 1.05 T c ). One observes that even near its peak at τ ∼ 38 fm/c, R −1 Π evolves almost identically in anisotropic and standard viscous hydrodynamics with QP transport coefficients. 12 A marked difference is observed, however, when the system is evolved with standard hydrodynamics using transport coefficients from a massless Boltzmann gas without a mean field (green short-dashed lines in Figs. 4d,e). It turns out that the thermodynamic integral β Π in Eq. (82) is remarkably sensitive to the degree of nonconformality of the Boltzmann gas, giving rise to a much longer bulk viscous relaxation time in the QP model than for the light Boltzmann gas without a mean field, especially in the neighborhood of T c . This is reflected in the large difference between the short-dashed green line and the other two curves for the bulk Knudsen number shown in Fig. 4d which causes the corresponding large difference in the evolution of the bulk inverse Reynolds number shown in panel (e): The much shorter relaxation time for the light Boltzmann gas allows the bulk viscous pressure to follow its Navier-Stokes limit (shown in Fig. 4f) much more closely, causing R −1 Π to rise much more steeply and to a larger peak value as the system cools towards T c than in the other two approaches where β Π is calculated from the QP model.
We have studied thermal equilibrium initial conditions with several other combinations of initial temperature T 0 and τ 0 , resulting in significantly different evolutions of the energy density and viscous pressure components (not shown here). Two features appear to be universal, however: (i) As long as we use transport coefficients computed from the same microscopic QP kinetic theory, the evolution of all components of the energy-momentum tensor, as well as of the mean field B, shows only very small differences (of the same order as shown in Fig. 3) between anisotropic and standard viscous hydrodynamics. (ii) Using instead transport coefficients for a Boltzmann gas of light fixed-mass particles, standard viscous hydrodynamics leads to significantly different evolutions for the bulk viscous pressure Π. For a meaningful comparison between anisotropic and standard viscous hydrodynamics it is therefore important that a consistent set of transport coefficients is being employed. Also, for a medium with broken conformal invariance (such as the quark-gluon plasma and other forms of QCD matter) non-conformal effects on the transport coefficients can have a large effect on the evolution of the bulk viscous pressure which may not be properly captured when using transport coefficients derived from a theory with weakly interacting degrees of freedom that have small masses.

B. Glasma-like initial conditions
In this subsection we repeat the exercise of the previous one for a different set of initial conditions, resembling those that one would get from matching the hydrodynamic evolution to a pre-equilibrium stage described by the IP-Glasma model [37]. As already described, this model predicts approximately vanishing initial longitudinal pressure P L ≈ 0 and P ⊥ ≈ E/2 [75]. (In practice, we set P L /P ⊥ = 0.01 initially). We use the same initial longitudinal proper time and temperature as before. The corresponding results are plotted in Figs. 5 and 6.
For this extreme initial condition, the default magnitude of B must be reduced by about 85% in order to be able to successfully initialize the anisotropic microscopic parameters; for B this is shown in Fig. 5c while the implications for the anisotropic microscopic parameters will be discussed in the following subsection. The highly non-equilibrium initial conditions manifest themselves in large starting values for the shear and bulk stresses and the non-equilibrium mean field. The initial shear stress (Fig. 5b) is about five times larger than its peak value for equilibrium initial conditions. The bulk viscous pressure (Fig. 5d) and non-equilibrium part of the mean field ( Fig. 5f) are for the first fm/c one to two orders of magnitude larger than for equilibrium initial conditions. In spite of this, anisotropic and standard viscous hydrodynamics still lead to almost identical evolution trajectories for the energy density (Fig. 5a) and viscous pressures (Figs. 5b,d) if QP transport coefficients are used, and if the latter are swapped out for those from a light Boltzmann gas, a significant change in the standard viscous hydrodynamic evolution is only seen for the bulk viscous pressure (short-dashed green curve in Fig. 5d). The shear stress π and the pressure ratio P L /P ⊥ emphasize the differences in the hydrodynamic models somewhat at early times (Figs. 5b,e), pushing the pressure ratio towards isotropy somewhat faster in anisotropic than in standard viscous hydrodynamics, but all three dynamical approaches converge to a common late-time behavior for π and P L /P ⊥ after about 2 fm/c (i.e. after about 3 times the initial shear relaxation time of about 0.8 fm/c). It is, however, not the case that equilibrium and Glasma-like initial conditions lead to the same temperature evolution of the system: A careful comparison of Figs. 3a and 5a shows that for the non-equilibrium initial conditions viscous heating by the large initial bulk and shear stresses causes the energy density (and therefore temperature) to drop somewhat more slowly than for equilibrium initial conditions, especially at early times. Figure 5f shows again that the "asymptotic approximation" (75) for δB (asy) (dashed purple line) is not a good approximation for the full numerical evolution of δB shown by the solid red line. Since for the Glasma-like initial conditions the non-equilibrium mean field contribution δB is initially of the same order of magnitude as the equilibrium contribution B eq , the breakdown of this approximation is visible even in the evolution of the total mean field B (solid red line) which is not at all described by B (asy) ≡ B eq +δB (asy) . Looking at the Knudsen and inverse Reynolds numbers in Fig. 6 the only striking (although obvious) difference are the large starting values for both shear and bulk inverse Reynolds numbers when using Glasma-like initial conditions. Similar to the shear stress π and pressure ratio P L /P ⊥ in Figs. 5b,e, these two observables exhibit noticeable differences at early times between anisotropic and standard viscous hydrodynamic evolution.

C. Evolution of the microscopic kinetic parameters
Although the parameters (Λ, α ⊥ , α L ) describing the slope and anisotropy of the momentum distribution of the microscopic degrees of freedom are vestiges from an underlying kinetic theory whose traces we have tried to erase as much as possible in our formulation of anisotropic hydrodynamics (hoping that eventually we can obtain the transport coefficients of QCD matter from a more fundamental approach), it is interesting to "look under the hood" and see how our parametrized anisotropic EOS works, i.e. how the QP model adjusts its microscopic parameters to accommodate the macroscopic anisotropic hydrodynamic initial conditions provided, and how it evolves them in response to the anisotropic hydrodynamic evolution of the energy-momentum tensor.
Figs. 7a,b compare, for equilibrium initial conditions, the evolution of the effective temperature parameter Λ with that of the true temperature T extracted from the energy density, and of the momentum anisotropy parameter α L with that of α ⊥ , respectively. A comparison of Figs. 7a,b with Figs. 4b,e shows that large inverse Reynolds numbers in both the shear and bulk sectors correlate with effective temperatures Λ > T and longitudinal momentum deformation parameter α L < 1. Large shear inverse Reynolds numbers correlate with α ⊥ deviating from unity in the opposite direction (i.e. with α ⊥ > 1), leading to narrower longitudinal and wider transverse momentum distributions than in the equilibrium distribution f eq , consistent with P L /P ⊥ < 1. Large bulk inverse Reynolds numbers push down both α L and α ⊥ , corresponding to negative bulk viscous pressures. At late times, when both the shear and bulk inverse Reynolds numbers approach zero, the momentum distribution approaches local equilibrium, α L,⊥ → 1 and Λ → T .
For Glasma-like initial conditions, shown in Figs. 7c,d, these generic statements for the deformation parameters α remain true but at early times the relationship between T and Λ is completely changed: the effective temperature Λ starts out much smaller than the true temperature. A low effective temperature Λ would narrow the microscopic momentum distribution in the transverse plane if it were not compensated by a very large (O(10)) initial value of α ⊥ , which upholds the kinetic energy density and transverse pressure. On the other hand, α L starts out almost at zero, reinforcing the narrowing of the longitudinal momentum distribution generated already by the small Λ value and thereby causing a very small ratio of the kinetic contributions to P L and P ⊥ . This is, of course, forced upon the system by the very anisotropic initial condition P L /P ⊥ = 0.01.
While the microscopic kinetic parameters (Λ, α) control only the kinetic contributions to energy density and pressures, the qualitative agreement of their tendencies extracted from this analysis of Figs. 7c,d with those of the total energy density and pressures shown in Fig. 5 demonstrates that the mean field B, even where large, cannot alter the sign of the pressure anisotropy (shear stress). Its value shifts the average kinetic pressure relative to the kinetic energy density and thereby has a large influence on the bulk viscous pressure.

VI. CONCLUSIONS AND OUTLOOK
In this work we presented a purely macroscopic formulation of anisotropic hydrodynamics in 3+1 space-time dimensions, parametrized with Milne coordinates. To obtain the Lorentz structure of the anisotropic hydrodynamic equations, including the relaxation equations for the dissipative flows, we started from a microscopic description in terms of a relativistic Boltzmann-Vlasov equation with a relaxation-time approximated collision term. The mean field in the Boltzmann-Vlasov equation is constructed such that the energy density and equilibrium pressure of this kinetic theory satisfy an equation of state that agrees with the lattice QCD EOS of strongly interacting matter. The macroscopic equations of motion are derived from an anisotropic moment expansion of this Boltzmann-Vlasov equation, where the distribution function is split into a momentum-anisotropic leading order term f a and a residual correction δf . To close the anisotropic moment expansion we use for the residual correction δf the 14-moment approximation. The leading-order term is constructed such that it can nonperturbatively account for the two largest dissipative effects encountered in relativistic heavy-ion collisions, a large longitudinal-transverse pressure anisotropy at early times and a large bulk viscous pressure during the phase transition of the matter from a quark-gluon plasma to color-confined hadronic matter. This requires the introduction of two momentum-anisotropy parameters α L , α ⊥ into f a whose dynamics is fixed by a novel generalization of the Landau matching conditions that ensures that there are no residual corrections from δf to the longitudinal and transverse pressures of the system. This matching scheme allows us to completely eliminate the micro-scopic parameters that define f a , and to write down, for the first time, a set of macroscopic anisotropic hydrodynamic evolution equations which make no explicit reference at all to their microscopic kinetic origin.
There are ten evolution equations for the ten components of the energy-momentum tensor. No specific assumptions are made for the equation of state relating the energy density and thermal pressure in thermal equilibrium, i.e. the equations can be used to describe any form of matter that behaves like a dissipative fluid. Two of these equations evolve the longitudinal and transverse pressures P L and P ⊥ . Instead of splitting them into a thermal equilibrium pressure, a bulk viscous pressure and a longitudinal-transverse shear stress, with the latter two quantities assumed to be small and perturbatively treatable, in our approach P L and P ⊥ themselves are evolved, with the transport coefficients controlling how far they may deviate from the thermal equilibrium pressure.
The evolution of the energy-momentum tensor components is controlled by a standard set of driving forces, such as the longitudinal and transverse expansion rates, the various components of the velocity shear tensor, the flow vorticity, etc. In addition, the dissipative flows are characterized by a set of relaxation times describing their relaxation towards their first-order Navier-Stokes limits. Consistent with the anisotropic parametrization of the leading-order distribution f a , the dissipative forces are separated into longitudinal and transverse parts using a systematic procedure involving orthogonal projection operators that was developed in Ref. [42]. They are multiplied by a set of two dozen transport coefficients. These transport coefficients, as well as the relaxation times, are material properties of the dissipative fluid to be described.
We do not know how to compute these transport coefficients for QCD matter from first principles. Therefore we use in this work a kinetic theory for weakly-interacting quasiparticles with temperature-dependent masses as a model for computing them. We write the distribution function for these quasiparticles as f = f a +δf and parametrize f a in the same way as in the kinetic theory from which we first started. From the solution of the anisotropic hydrodynamic equations we then take the energy density E, longitudinal pressure P L and transverse pressure P ⊥ , as well as the mean field B, and describe the deviations of P L and P ⊥ from the equilibrium pressure P eq (E) (which is taken from lattice QCD) in terms of the microscopic anisotropic parameters of the kinetic model. Having thus fixed the parameters of the kinetic model from the macroscopic hydrodynamic output, we can use it to compute all the transport coefficients in kinetic theory.
For the relaxation times we take previously introduced phenomenological parametrizations that were recently calibrated by a global comparison of a sophisticated dynamical model involving dissipative relativistic fluid dynamics at its core with experimental heavy-ion collision data collected at the LHC [6,85].
As a first application of this new approach we have here studied the Bjorken evolution of a longitudinally boostinvariant, transversely homogeneous system, evolving it both with anisotropic and standard viscous hydrodynamics for comparison. We found remarkable agreement between the two approaches if both used quasiparticle transport coefficients but noticeable disagreements with a standard viscous hydrodynamic simulation using transport coefficients for a weakly interacting Boltzmann gas in the small fixed-mass limit. This suggests an unexpected robustness of the standard viscous hydrodynamic approach even in the presence of large shear and bulk viscous effects. A final assessment of the relative strengths and weaknesses of anisotropic vs. standard viscous hydrodynamics will, however, have to await the availability of full (3+1)-dimensional numerical evolution comparisons which are presently being pursued.
A key motivation for anisotropic hydrodynamics is that by accounting for the large dissipative components already at leading order, by parametrizing them into the leading-order distribution function f a , the remaining dissipative flows arising from the residual deviation δf in the decomposition f = f a +δf should be smaller than the dissipative flows in standard viscous hydrodynamics where they arise from δf in the decomposition f = f eq +δf . Bjorken flow does not allow to test this expectation because for systems with Bjorken symmetry the residual dissipative flows arising from δf vanish anyhow exactly by symmetry. Full (3+1)-dimensional simulations will allow to answer this question. Taking the results reported in the last chapter of Ref. [72] (based on a version of the present framework that did not treat the bulk viscous pressure non-perturbatively) as guidance for what to expect, anisotropic hydrodynamics as formulated here should indeed make the residual shear stress components significantly smaller in the center of the fireball, where the largest shear stresses are generated at early times by longitudinal expansion. However, the same may not necessarily hold for cells near the transverse edge of the fireball where the transverse expansion rate can exceed the longitudinal one and where accounting non-perturbatively for large effects associated with anisotropies relative to the beam axis (as we do here) may not offer significant advantages. We hope to be able to soon present numerical results that show how these expectations bear out in practice.
We close by noting that the observed sensitivity of the Bjorken evolution to the chosen model for computing the transport coefficients puts some urgency to the question how to compute the transport coefficients of anisotropic hydrodynamics from first principles for a theory such as hot and dense QCD. We have to leave this as a challenge for future work. withp = p/Λ, the angular integrals in Eqs. (A1,A2) can be evaluated analytically, yielding I nrqs = α 2q+2 ⊥ α r+1 L Λ n+s+2 4π 2 (2q)!! ∞ 0 dpp n+s+1 × R nrq (α ⊥ , α L ;m/p)f eq ( p 2 +m 2 ), (A5) wherem = m/Λ and the functions R nrq are defined as R nrq (α ⊥ , α L ;m/p) = w n−r−2q−1 × 1 −1 d cos θ sin 2q θ cos r θ(1 + z sin 2 θ) (n−r−2q−1)/2 , with w = α 2 L + (m/p) 2 and z = α 2 ⊥ − α 2 L w 2 . The radial momentum integral can be computed numerically with generalized Gauss-Laguerre quadrature. For reference, we list the functions R nrq that are needed in this paper: and write the term f afa as Integrating by parts with respect to cos θ, where the boundary term vanishes for q ≥ 1, gives the following relation: (r−1) I n−2,r−2,q,s+1 −I n−2,r,q−1,s+1 + Λ (s+1) I n,r,q,s−1 , which for (n, r, q, s) = (4, 2, 1, −1) yields Eq. (73c). Those controlling the evolution of the total transverse pressure P ⊥ arē Appendix E: Viscous hydrodynamic equations Here we derive the viscous hydrodynamic equations (83) and their transport coefficients. We start with the quasiparticle case (long-dashed blue) and derive the relaxation equation for δB and its second-order solution (84). The general evolution equation for δB is given by [52] where theṁ/m term can be written aṡ We replace the time derivativeĖ with the energy conservation law in viscous hydrodynamicṡ E + (E + P eq + Π)θ − π µν σ µν = 0 , where θ = ∂ µ u µ is the scalar expansion rate and σ µν = ∆ αβ µν ∂ β u α is the velocity shear tensor. For the secondorder relaxation equation, we only need the first-order approximationĖ ≈ −(E + P eq )θ, thuṡ The equation of motion for δB then reduces to To first order in deviations from equilibrium δB = 0. The second-order solution is given by (84) [52]. In Eq. (E5), we truncate the third-order term ∝ δB θ to arrive at the second-order relaxation equation