Recent progress in anisotropic hydrodynamics

The quark-gluon plasma created in a relativistic heavy-ion collisions possesses a sizable pressure anisotropy in the local rest frame at very early times after the initial nuclear impact and this anisotropy only slowly relaxes as the system evolves. In a kinetic theory picture, this translates into the existence of sizable momentum-space anisotropies in the underlying partonic distribution functions,<<. In such cases, it is better to reorganize the hydrodynamical expansion by taking into account momentum-space anisotropies at leading-order in the expansion instead of as a perturbative correction to an isotropic distribution. The resulting anisotropic hydrodynamics framework has been shown to more accurately describe the dynamics of rapidly expanding systems such as the quark-gluon plasma. In this proceedings contribution, I review the basic ideas of anisotropic hydrodynamics, recent progress, and present a few preliminary phenomenological predictions for identified particle spectra and elliptic flow.


Introduction
The phenomenological application of viscous hydrodynamics to the dynamics of the quark-gluon plasma (QGP) created in heavy-ion collisions has been tremendously successful [1][2][3]. Despite this success, there are spacetime regions where standard viscous hydrodynamics is being pushed beyond its limits, such as at very early times after the initial heavy-ion collision τ < 1 fm/c and near the (semi)-dilute edges of the system. In these spacetime regions, viscous hydrodynamics itself tells you that there may be trouble, since the shear Knudsen number, Kn π = τ π ∂ µ u µ , and the inverse Reynolds number, R −1 π = π µν π µν /P eq , can become quite large. The situation gets worse as the beam energy is decreased or one considers small collision systems such as pA and pp. Additionally, when one wants to compute observables other than soft hadron production, such as heavy quarkonium suppression, photon emission, dilepton emission, etc. one traditionally employs a kinetic description which requires knowledge of the full momentum dependence of the underlying parton distribution function(s), f (x, p). In standard viscous hydrodynamics approaches, f (x, p) is expressed as power series with the leading term being the isotropic thermal contribution and the shear and bulk corrections to isotropic equilibrium being expressed as polynomials in the momenta. As a result, one cannot necessarily trust the high-momentum limit of production rates since this maps to regions in which the partonic distributions functions can become negative.
In order to address these problems, in Refs. [4,5] it was first shown that, in the context of relativistic transport theory, one could change the ansatz for the expansion point for the distribution function and use this reorganized expansion to derive so-called anisotropic hydrodynamics (aHydro). In Ref. [5], in particular, it was shown that, for a 0+1d boost invariant system, aHydro could reproduce arXiv:1611.05056v1 [nucl-th] 15 Nov 2016 EPJ Web of Conferences both the ideal and free streaming limits of the dynamics and that, in the limit of small anisotropies, the dynamical equations reduced identically to those of Israel-Stewart viscous hydrodynamics. Since these two original papers, there has been a great deal of progress in anisotropic hydrodynamics [4][5][6][7][8][9][10][11][12][13][14][15][16] including applications to cold atomic gases near the unitary limit [17,18]. In parallel, there have been efforts to construct exact solutions to the Boltzmann equation in some simple cases which can be used to test the efficacy of various dissipative hydrodynamics approaches and it has been shown that anisotropic hydrodynamics most accurately reproduces all known exact solutions even in the limit of very large η/s and/or initial momentum-space anisotropy [15,[19][20][21][22][23][24]. The recent focus has been on turning aHydro into a practical phenomenological tool which has a realistic equation of state and self-consistent anisotropic hadronic freeze-out in order to make comparisons to experimental data. In this proceedings, I review the moment formulation of aHydro which begins with an ansatz for the leading-order one-particle distribution that contains a symmetric anisotropy tensor which allows for multiple anisotropy parameters. I present some preliminary comparisons with experimental data from LHC 2.76 TeV Pb-Pb collisions and provide an outlook for the future.

Distribution function ansatz
Typically, when one derives dissipative hydrodynamics from transport, the starting assumption is that the distribution function can be expanded around an assumed state of isotropic equilibrium with all non-equilibrium corrections collected into "perturbative corrections" which are typically expressed as an orthogonal polynomial expansion in the momentum contracted with the viscous stress tensor [25,26]. In aHydro, one instead replaces the leading-order term in this series by an anisotropically deformed distribution function of the form [6,8,9] f (x, p) ≡ f eq p µ Ξ µν (x)p ν /λ(x) , with λ being a local temperature-like scale and Ξ µν ≡ u µ u ν + ξ µν − Φ∆ µν being the anisotropy tensor, which parametrizes the anisotropic form of distribution function. In the definition of Ξ µν , u µ is fluid four-velocity, ξ µν is a traceless symmetric anisotropy tensor, Φ is related to the bulk-viscous pressure correction, and ∆ µν ≡ g µν − u µ u ν projects out components of a four-vector which are transverse to u µ . The quantities entering the ansatz obey u µ u µ = 1, ξ µ µ = 0, and u µ ξ µν = u µ ∆ µν = 0. Since ξ µν is transverse to u µ it can be expanded in terms of spacelike basis vectors X µ , Y µ and Z µ which satisfy u µ X µ = u µ Y µ = u µ Z µ = 0. The ansatz (1) contains information about both shear and bulk corrections at leading order [9]. In what follows, we assume the distribution to be of Boltzmann form, i.e. f iso (p) ≡ exp (−p) and that the anisotropy tensor ξ µν is diagonal in the local rest frame (LRF) Since ξ µν is traceless, there are only two independent anisotropies. In practice, it is convenient to parameterize these two independent degrees of freedom and the Φ parameter in terms of three momentum-space ellipticities in which case, under the assumption that the anisotropy tensor is diagonal in the LRF, allows us to write simply

Equations of motion
In anisotropic hydrodynamics, the equations of motion are obtained from moments of the Boltzmann equation with the collisional kernel taken in the relaxation time approximation (RTA) In what follows, we have used the first and second moments, which result in equations of the form where T µν is the energy-momentum tensor consistent with Eq. (4) and I µ 1 ···µ n is defined by where N dof is the number of degrees of freedom. The general equilibrium tensor I µ 1 ···µ n eq is the same with f (x, p) → f eq (|p|, T (x)). In terms of the general moments, I µ 1 ···µ n , one has N µ = I µ , T µν = I µν , etc. With the assumption that the anisotropy tensor is diagonal, we need equations of motion for the three ellipticities, three independent components of u µ , and λ. We obtain these from the four equations for energy-momentum conservation and three diagonal projections of the second moment equation. The local effective temperature T (x) is determined by Landau-matching, which requires that the equilibrium and non-equilibrium energy densities are the same at all points in spacetime, i.e. ε( α, λ) = ε eq (T ).

Equation of State
If the QGP is anisotropic in momentum-space it is not obvious how one should implement a realistic lattice-based equation of state since this is an implicitly equilibrium concept. There are currently two approaches being followed in the literature. In the first, dubbed the standard approach, one obtains the dynamical equations necessary in the conformal limit, m → 0. In this limit, the components of T µν multiplicatively factorize and one can then connect the isotropic parts of energy density and pressures using a realistic EoS [27]. In the second approach, dubbed the quasiparticle approach, one assumes that the system is comprised of quasiparticles with a temperature-dependent mass, m → m(T ), which is fit to lattice QCD data [13,16]. The quasiparticle method is more self-consistent in the way breaking of conformal invariance is implemented, however, at the moment it is infeasible to use for 3+1d simulations. For this reason, in the results presented herein we will use the standard method for implementing the equation of state. For the underlying realistic equilibrium EoS we used the Krakow parametrization of lattice data which is matched onto a hadron resonance gas at low temperatures [28].

Freezeout
As the system cools, the QGP undergoes a crossover from quark and gluons to hadronic degrees of freedom which subsequently freeze out kinetically when their mean free path becomes large. In order to compare the result of hydrodynamic models to experimental data, one needs to calculate the differential particle spectra at freezeout. In practice, we construct a constant energy density hypersurface, defined through the effective temperature T FO = ε −1 (ε FO ). Then, by computing the number of particles that cross this hypersurface, one can determine the number of hadrons produced in heavy-ion collisions at the freezeout using where i labels the hadronic species, N i ≡ (2s i +1)(2g i +1) is the degeneracy factor with s i and g i being the spin and isospin of the state in question, and f i is the distribution function for particle species i taking into account the appropriate quantum statistics. For each species, we assume that the same anisotropic distribution function form can be used [16,27].

Initial Conditions
In order to solve the aHydro dynamical equations one has to choose initial conditions at τ = τ 0 , i.e. for 3+1d aHydro with only diagonal anisotropies in the LRF, one needs seven three-dimensional profiles: where ς is the spatial rapidity. In this work, we assume that the initial transverse profile for the effective temperature (determined via Landau matching) is given by the optical Glauber model. We assume that the initial energy density is proportional to the scaled initial density of the sources such that the initial effective temperature is where b is the impact parameter and the proportionality constant ε 0 is chosen in such a way as to correspond to a given central temperature, T 0 = ε −1 (ε 0 ). The density of sources is constructed using the following mixed model where ρ ± WN is the density of wounded nucleons from the left/right-moving nuclei and ρ BC is the density of binary collisions, both of which are obtained using the optical limit of the Glauber model The longitudinal profile is taken to be For the LHC case studied here, we use κ = 0.145 for the mixing factor and an inelastic crosssection of σ in = 62 mb. The parameters in the longitudinal profile (13) were fitted to reproduce the pseudorapidity distribution of charged particles, with the results being ∆ς = 2.5 and σ ς = 1.4 and the shift in rapidity was calculated according to the formula [29]  where all functions are understood to be evaluated at a particular value of b and x ⊥ . The participant velocity is defined as v P ≡ ( √ s/2) 2 − (m N /2) 2 /( √ s/2), m N is the nucleon mass, and √ s is the center-of-momentum collision energy. In Eq. (12) we have made use of the thickness function where the nuclear density is given by the Woods-Saxon profile For Pb-Pb collisions, we use ρ 0 = 0.17 fm −3 for the nuclear saturation density, R = 6.48 fm for the nuclear radius, and a = 0.535 fm for the surface diffuseness of the nucleus. In the calculations presented herein, we assumed that the produced matter has no initial transverse flow, i.e. u x (τ 0 , x ⊥ , ς) = u x (τ 0 , x ⊥ , ς) = 0, while the initial longitudinal flow is of Bjorken form u z (τ 0 , x ⊥ , ς) = z/t. For simplicity, in this work the initial anisotropy parameters are assumed to be homogeneous and isotropic, α(τ 0 , x ⊥ , ς) = 1.

Numerical methods
We solve the 3+1d aHydro equations on a lattice with spacing a ⊥ = 0.2 fm and a ζ = 0.2 with N x = N y = N ς = 128. For temporal updates we used 4th-order Runge-Kutta with a temporal step size of = 0.01 fm/c. A weighted LAX scheme with a rather small weighting factor of λ = 0.01 was used to regulate possible numerical instabilities associated with shock-wave formation [6]. We varied the centrality of the collision by setting the impact parameter to the average value expected in each centrality class according to the optical Glauber model as indicated in Table 1. Each of the configurations was evolved until the maximum effective temperature in the entire volume was below 120 MeV. From each of these results, a freeze-out hypersurface corresponding to a fixed effective temperature was extracted. The microscopic parameters α, u, and λ on the hypersurface were then fed to a version of THERMINATOR 2 which had been modified to account for an ellipsoidal distribution function [30]. THERMINATOR 2 produces sampled event-by-event hadronic production from the exported freezeout hypersurface and also performs hadronic feed down (resonance decays) for each sampled event.  Figure 1. Comparison of aHydro model predictions with experimental data for π + (blue), K + (red), and p production (green). The left panel shows the particle spectra in the 0-5% centrality class and the right panel shows v 2 in the 20-30% centrality class, both as a function of transverse momentum p T . The best fit corresponded to T 0 = 0.56 GeV, T FO = 0.13 GeV, and 4πη/s = 3. Data shown are from the ALICE collaboration [31,32]. Experimental error bars shown are statistical only.

Results
Herein we consider LHC 2.76 TeV Pb-Pb collisions and compare to experimental data available from the ALICE collaboration [31,32]. Using the initial conditions specified previously, we varied 4πη/s ∈ {1, 2, 3} and, in each case, varied the initial central temperature T 0 and the freeze-out temperature T FO and compared with ALICE data. From our preliminary analysis, the values of T 0 560 MeV and T FO 130 MeV gave the best fit to the spectrum and v 2 .
In Fig. 1 we show the resulting π + , K + , and p transverse momentum spectrum in the 0-5% centrality class (left) and v 2 in the 20-30% centrality class (right). As can be seen from this figure, the model does a very good job reproducing the identified particle v 2 in the 20-30% centrality class; however, we see significant deviations in the low-p T π + , K + , and p spectra. This discrepancy is most likely a result of the "standard method" for implementing a realistic equation of state. As was shown in Ref. [16], one finds that the standard method significantly underestimates the number of low p T hadrons, whereas the quasiparticle method for implementing the equation of state is in better agreement with standard second-order viscous hydrodynamics at low p T .
In Fig. 2 we show two results. In the left panel, we show the charged particle v 2 as a function of p T in the 30-40% centrality class and on the right we show the p T -integrated charged particle v 2 as a function of centrality. As can be seen from the left panel, there is reasonable agreement between the aHydro model employed herein and the charged particle v 2 . In the right panel, we show three results corresponding to 4πη/s ∈ {1, 2, 3} which are indicated in red, blue, and green, respectively, and are ordered from top to bottom at 35% centrality, respectively. As we can see from the right panel there is disagreement between the model and the data for very central events. This is to be expected since we used smooth optical Glauber initial conditions. All v 2 generated in central collisions is due to fluctuations and, in the first two centrality classes, there are significant contributions to v 2 from initial state fluctuations. In the peripheral centrality classes we also note that the model prediction falls more quickly to zero than the data seem to indicate; however, we note there is a large systematic error associated with the reported v 2 in the most peripheral centrality class. Finally, we note that the maximum model sensitivity to the assumed values of η/s occurs at around 40% centrality.  Figure 2. Comparison of aHydro model predictions with experimentally measured charged particle v 2 as a function of p T in the 30-40% centrality class (left) and the p T -integrated charged particle v 2 as a function of centrality (right). Data shown are from the ALICE collaboration [31,32]. Experimental error bars shown are statistical only.

Conclusions
In this proceedings, I have reviewed recent progress in anisotropic hydrodynamics. The anisotropic hydrodynamics framework has been checked against exact solution of the Boltzmann equation in some non-trivial but simple cases in which exact solution is possible and one finds that, in all cases considered, the anisotropic hydrodynamics approach most accurately reproduces the evolution of the system. Building on this success, the current focus of the aHydro program is to turn the formalism into a practical phenomenological tool that can be used to model the non-equilibrium dynamics of the QGP created in AB nuclear collisions including pA and, down the road, pp collisions. As demonstrated herein, significant progress in this direction has been made. We now have functioning 3+1d versions of leading-order aHydro which include multiple anisotropy parameters that are associated with the shear and bulk corrections to the distribution function.
The code now includes "anisotropic freeze-out" which is tightly integrated with a customized version of THERMINATOR 2, which takes care of the final hadronic production and resonance decays. Comparisons with ALICE data indicate that, even given the fact that we are using smooth Glauber initial conditions, the model is able to reproduce the qualitative features of the data overall, but in the case of the p T -differential v 2 we find very promising agreement between the 3+1d aHydro results and experimental data. The key standout phenomenologically lies in our inability to properly fit the low p T part of the hadronic spectra. We have suggested that this is due to an inconsistency inherent in the standard method for implementing the equation of state in aHydro. Looking to the future, we are currently working on optimizing the quasiparticle aHydro approach in order to make it feasible to use in 3+1d simulations. Initial testing shows that this will be possible and offers hope that one can have a self-consistent way to implement conformal symmetry breaking while at the same time having an efficient code. This will become critical as we begin explorations of fluctuating initial conditions.