Anisotropic fluid dynamical simulations of heavy-ion collisions
AAnisotropic fluid dynamical simulations ofheavy-ion collisions
Mike McNelis a, ∗ , Dennis Bazow a , Ulrich Heinz a,b a Department of Physics, The Ohio State University, Columbus, OH 43210-1117, USA b Institut f¨ur Theoretische Physik, J. W. Goethe Universit¨at, Max-von-Laue-Str. 1,D-60438 Frankfurt am Main, Germany
Abstract
We present
VAH , a (3+1)–dimensional simulation that evolves the far-from-equilibrium quark-gluon plasma produced in ultrarelativistic heavy-ion colli-sions with anisotropic fluid dynamics. We solve the hydrodynamic equationson an Eulerian grid using the Kurganov–Tadmor algorithm in combinationwith a new adaptive Runge–Kutta method. Our numerical scheme allows usto start the simulation soon after the nuclear collision, largely avoiding theneed to integrate it with a separate pre-equilibrium dynamics module. Wetest the code’s performance by simulating on the Eulerian grid conformaland non-conformal Bjorken flow as well as conformal Gubser flow, whose(0+1)–dimensional solutions are precisely known. Finally, we compare non-conformal anisotropic hydrodynamics to second-order viscous hydrodynamicsin central Pb+Pb collisions and find that the former’s longitudinal flow pro-file responds more consistently to the fluid’s gradients along the spacetimerapidity direction.
Keywords:
Ultrarelativistic heavy-ion collisions, quark-gluon plasma,relativistic hydrodynamics, computational fluid dynamics ∗ Corresponding author.Email addresses: [email protected] (M. McNelis), [email protected] (U. Heinz)
Preprint submitted to Computer Physics Communications January 11, 2021 a r X i v : . [ nu c l - t h ] J a n ROGRAM SUMMARY
Manuscript Title:
Anisotropic fluid dynamical simulations of heavy-ion collisions
Authors:
Mike McNelis, Dennis Bazow, Ulrich Heinz
Program Title:
VAH
Licensing provisions:
GPLv3
Programming Language:
C++
Computer:
Laptop, desktop, cluster
Operating System:
GNU/Linux distributions, Mac OS X
Global memory usage: × ×
63 grid)
Keywords:
Ultrarelativistic heavy-ion collisions, quark-gluon plasma, relativistichydrodynamics, computational fluid dynamics
Classification:
12 Gases and Fluids, 17 Nuclear physics
External routines/libraries:
GNU Scientific Library (GSL)
Nature of problem:
Modeling the far-from-equilibrium dynamics of quark-gluon plasma produced inultrarelativistic heavy-ion collisions.
Solution method:
Kurganov–Tadmor algorithm, adaptive stepsize method
Running time:
A (3+1)–d non-conformal anisotropic fluid dynamical simulation of a central Pb+Pbcollision on a 129 × ×
63 grid takes about 530 s for an Intel Xeon E5-2680 v4multi-core processor with OpenMP acceleration (see Sec. 5.3). ontents1 Introduction 52 Anisotropic fluid dynamics 10 T R ENTo initial conditions . . . . . . . . . . . 554.5.2 Fluctuating T R ENTo initial conditions . . . . . . . . 59 T R ENTo energy deposition model 79Appendix E Numerical implementation of second-order viscoushydrodynamics 81
E.1 Hydrodynamic equations . . . . . . . . . . . . . . . . . . . . . 81E.2 Transport coefficients . . . . . . . . . . . . . . . . . . . . . . . 83E.3 Reconstruction formulas for the energy density and fluid velocity 84E.4 Regulating the shear stress and bulk viscous pressure . . . . . 854 . Introduction
The quark-gluon plasma is one of the most extreme phases of matter inour universe. The incredibly high temperatures and densities required tocreate it ( T ∼ K, ρ ∼ kg/m ) are orders of magnitude beyondthose of any substance typically encountered in nature. The quark-gluonplasma also has the remarkable property of having the lowest shear viscos-ity to entropy density ratio of any known fluid; current phenomenologicalconstraints put η/ S ∼ . T c = 154 MeV [1, 2]. This value is very close to the conjecturedKovtun-Son-Starinet bound η/ S ≥ / π , which follows from the AdS/CFTcorrespondence in string theory [3, 4]. Therefore, it is of great interest tothe physics community to test whether or not the KSS bound holds for suchstrongly coupled liquids. Today, we can produce quark-gluon plasma bycolliding beams of heavy ions (e.g. Pb+Pb) at very high energies. The fem-toscopically small droplets of plasma formed in these ultrarelativistic nuclearcollisions have incredibly short lifetimes ( τ f ∼ − s), making them diffi-cult to probe directly. To reconstruct the medium’s transport properties,one can simulate the various stages of a heavy-ion collision with quantitativeprecision and predictive power and fit the theoretical model to large sets ofsoft-momentum ( p T < The hybridmodel must be accurate enough that it can make full use of the considerableprecision of the available experimental data.The anisotropically expanding quark-gluon plasma stage is usually mod-eled with relativistic viscous hydrodynamics [5, 16–20]. Conventionally, vis-cous hydrodynamics only works for fluids that are both near local equilibrium(i.e. inverse Reynolds number Re − (cid:28)
1) and have small spacetime gradients(i.e. Knudsen number Kn (cid:28)
1) [18, 21]. But unlike most fluids found in na-ture, the quark-gluon plasma is initially far from local-equilibrium (Re − ∼ ∼
1) for a good portion of itslifetime. This is primarily due to quantum fluctuations in the initial-stateprofile [22–24] and the rapid longitudinal expansion at early times. Underthese extreme conditions, the validity of a fluid dynamical description for thequark-gluon plasma comes into question [25]. Nevertheless, second-order vis- High-energy jets, heavy quarks and electromagnetic radiation can also serve as exper-imental probes but their full integration into hybrid models [10–15] is more involved thanfor soft-momentum hadrons. v n ) [26, 27]. State-of-the-art hybrid models nowhave enough predictive power to quantitatively constrain the shear and bulkviscosities of the quark-gluon plasma using heavy-ion experimental data andBayesian inference [1, 2, 6–9]. Hydrodynamics has also been found to beapplicable in collisions between light and heavy ions (e.g. He +Au, p+Pb)and in proton–proton collisions [12, 28], although the origin of collectivity inthese small systems is still being debated [29–32].The success of second-order viscous hydrodynamics can be attributed toan underlying resummed hydrodynamic theory that also extends to largegradients [33–35]. As the associated non-hydrodynamic modes decay on mi-croscopic time scales τ r (cid:28) τ hydro , the system approaches a (generally non-equilibrium) hydrodynamic attractor that is well approximated by viscoushydrodynamics. This happens within the time scale τ hydro ∼ c , longbefore the system thermalizes [38–40]. Whether or not the strongly coupledquark-gluon plasma possesses an attractor at presently experimentally ac-cessible temperatures ( T ∼ . − . (cid:29)
1) [20, 47,48]. The largest gradients in heavy-ion collisions occur at very early times( τ < . c ), especially around the edges of the fireball. If one starts the vis-cous hydrodynamic simulation too early, one usually encounters negative lon- Hydrodynamic simulations often employ relaxation-type equations from relativistickinetic theory (e.g. DNMR [18]), but they likely do not precisely capture the transientdynamics in strongly coupled fluids [36, 37]. The hydrodynamization time τ hydro refers to the time when viscous hydrodynamicsbecomes applicable, but the hydrodynamic simulation starting time τ can be differentfrom this. P L due to huge pressure anisotropies [49]; near the edgesof the fireball even the transverse pressure P ⊥ can turn negative if the bulkviscosity ζ/ S peaks strongly near the quark-hadron phase transition. Notonly does this cause excessive viscous heating but large negative pressuresalso redirect the matter inward, potentially resulting in cavitation. Regula-tion schemes can be used to tamp down the viscous pressures and prevent thesimulation from crashing, but they cloud the physical predictions of the orig-inal hydrodynamic theory [5, 20]. Because of the technical issues involved, itis more suitable to use a pre-equilibrium dynamics model before transitioningto viscous hydrodynamics at τ = τ hydro [1, 2, 6, 7, 22, 26, 44, 45, 50–52].Still, the conformal approximation usually made in the former results in amismatch to the latter’s non-conformal equation of state, producing (in spiteof the system’s expansion ) artificially positive bulk viscous pressures that canbe as large as Π ∼ P eq around the edges of the fireball [53]. In some situ-ations, one cannot even switch between dynamical models instantaneously.For example, in low-energy collisions ( √ s NN ∼ −
50 GeV) where the nu-clear interpenetration time is comparable to the fireball lifetime, one needsto run viscous hydrodynamics in the background as soon as the participantnucleons start feeding thermal energy and net-baryon number into the newlyformed fireball [54–58].The shortcomings of second-order viscous hydrodynamics have motivatedthe development of hydrodynamic models that can better handle far-from-equilibrium situations [59–61]. The most promising candidate is anisotropichydrodynamics, which treats the two largest dissipative effects arising inheavy-ion collisions (the pressure anisotropy P L − P ⊥ and the bulk vis-cous pressure Π) non-perturbatively [62–66]. While current formulationsof anisotropic hydrodynamics are partially based on weakly coupled kinetictheory, they are able to capture exact kinetic solutions more accurately thanstandard viscous hydrodynamics [47, 48, 62–64, 69, 70]. In addition, their Viscous heating refers to internal entropy production by dissipative processes, not dueto external thermal sources. Some regulation schemes are less extreme than others. For example, the hydrodynamiccode
MUSIC only regulates the dissipative currents outside the fireball region [16] while iEBE-VISHNU regulates the entire grid [5]. Second-order viscous hydrodynamic simulations also rely on microscopic approachessuch as relativistic kinetic theory to compute the relaxation times and other second-ordertransport coefficients [18, 67, 68]. √ s NN = 200 GeV) andPb+Pb collisions at the LHC ( √ s NN = 2 .
76 and 5.02 TeV) have been mod-eled reasonably well but so far only using smooth initial conditions and a fewmodel parameters to adjust the fit to experimental data [71–73]. Anisotropichydrodynamics has yet to undergo the same extensive tests and trials asviscous hydrodynamics, but it can potentially serve as an additional dis-crete model for the fluid dynamical stage of heavy-ion collisions. This wouldfurther increase the flexibility of the Bayesian framework developed by theJETSCAPE collaboration [1, 2], which has already incorporated a numberof discrete models for the particlization stage [74]. The hope is that, by con-trolling large dissipative flows nonperturbatively, anisotropic hydrodynamicscan constrain the transport coefficients of QCD matter more accurately.In this paper we introduce our (3+1)–dimensional anisotropic fluid dy-namical simulation for heavy-ion collisions called
VAH . The C++ mod-ule is based on the GPU–accelerated viscous hydrodynamic code
GPU VH [20, 49], except that it implements anisotropic hydrodynamics with the ( P L , P ⊥ )–matching scheme developed in Ref. [66]. We evolve the dynamical equa-tions on an Eulerian grid using the Kurganov–Tadmor algorithm [75], a pop-ular method also used in other viscous hydrodynamic codes [16, 17, 20, 76].The main improvement in our code is the ability to automatically adjustthe time step ∆ τ n after each iteration, as opposed to using a fixed valuefor the entire simulation. Our adaptive Runge–Kutta scheme initially usesa fine time step to resolve the rapid longitudinal expansion at early timeswhile speeding up the evolution at later times with a coarser time step givenby the Courant-Friedrichs-Lewy (CFL) condition. As a result, we can startanisotropic hydrodynamics at a very early time τ = 0 .
05 fm/ c to model boththe far-off-equilibrium dynamics stage with a non-conformal QCD equation The code package can be downloaded from the GitHub repository https://github.com/mjmcnelis/cpu_vah . At this point
VAH itself has not yet been ported to GPUs. We assume the system is longitudinally free-streaming (i.e. P L / P eq ≈ E ∝ /τ )in the time interval 0 < τ ≤ τ before starting anisotropic hydrodynamics. This closely
8f state and smoothly transition to viscous hydrodynamics. The code package contains other useful features, such as the option torun either anisotropic hydrodynamics or second-order viscous hydrodynam-ics within the same framework. This provides the user with several discretemodels to evolve the fluid dynamical stage for comparison. We further im-plemented user options to accelerate the simulation on a multi-core processorwith OpenMP, and to use automated grid settings that optimize the overallsize of the spatial grid for a given set of runtime parameters (e.g. the impactparameter and particlization switching temperature).This paper is organized as follows: In Sec. 2 we review the anisotropichydrodynamic equations used to evolve the energy-momentum tensor of thequark-gluon plasma. Sec. 3 discusses the numerical implementation of thesedynamical equations in the code. In Sec. 4 we test our anisotropic fluid dy-namical simulation for a system subject to (non)conformal Bjorken flow andconformal Gubser flow. Furthermore, we compare anisotropic hydrodynam-ics to second-order viscous hydrodynamics in (3+1)–dimensions. Finally,we benchmark the typical runtimes of (2+1)–d and (3+1)–d non-conformalhydrodynamic simulations in Sec. 5.In this work we use Milne spacetime coordinates x µ = ( τ, x, y, η s ), where τ = √ t − z is the longitudinal proper time and η s = tanh − ( z/t ) is thespacetime rapidity. We adopt the mostly-minus convention for the metrictensor, g µν = diag(1 , − , − , − /τ ). In the code, we solve the hydrody-namic equations in natural units (cid:126) = c = k B = 1; energy, momentum andtemperature have units of fm − , energy density and pressure have units offm − , etc. These units are converted to physical units (e.g. energy densitiesin GeV/fm ) when outputting and plotting the results. The net-baryon den-sity n B and baryon diffusion current V µB are neglected in this version of thecode. mirrors the situation found in other pre-hydrodynamic models [22, 45, 50]. It has also beenargued [77, 78] that (at least in weakly coupled systems) for τ → Anisotropic hydrodynamics reduces to second-order viscous hydrodynamics in thelimit of small Knudsen and inverse Reynolds numbers (Kn , Re − (cid:28) . Anisotropic fluid dynamics In this section we provide an overview of the anisotropic hydrodynamicsequations, including the equation of state and transport coefficients, to beimplemented in the code that evolves the energy-momentum tensor T µν ( x )of the quark-gluon plasma. For more details on the derivation of these hy-drodynamic equations we refer the reader to Refs. [65, 66]. In anisotropic fluid dynamics, it is convenient to decompose T µν ( x ) inthe basis { u µ ( x ), z µ ( x ), Ξ µν ( x ) } , where the fluid velocity u µ represents thetemporal direction in the local fluid rest frame (LRF) and the spatial vec-tors are split into the longitudinal direction z µ and the transverse projectiontensor Ξ µν = g µν − u µ u ν + z µ z ν [65]. In the Landau frame this results in thedecomposition (suppressing the spacetime dependence) T µν = E u µ u ν + P L z µ z ν − P ⊥ Ξ µν + 2 W ( µ ⊥ z z ν ) + π µν ⊥ , (1)where round parentheses denote symmetrization, i.e. W ( µ ⊥ z z ν ) = ( W µ ⊥ z z ν + W ν ⊥ z z µ ). The major components of T µν are the LRF energy density E = u µ u ν T µν , the longitudinal pressure P L = z µ z ν T µν and the transverse pres-sure P ⊥ = − Ξ µν T µν . Together, the pressures P L and P ⊥ capture thelargest dissipative flows in heavy-ion collisions: the pressure anisotropy ∆ P = P L − P ⊥ caused by the rapid longitudinal expansion rate at early longitu-dinal proper times, and the bulk viscous pressure Π = ¯
P − P eq due tocritical fluctuations near the quark-hadron phase transition. The remain-ing components of T µν are the longitudinal momentum diffusion current W µ ⊥ z = − Ξ µα z ν T αν and the transverse shear stress tensor π µν ⊥ = Ξ µναβ T αβ ,with Ξ µναβ = (Ξ µα Ξ νβ + Ξ νβ Ξ µα − Ξ µν Ξ αβ ) being the traceless double transverseprojector. We refer to these components as residual shear stresses since theyare typically smaller than the pressure anisotropy term in the full shear stresstensor π µν = 13 ( P L −P ⊥ )(2 z µ z ν +Ξ µν ) + 2 W ( µ ⊥ z z ν ) + π µν ⊥ . (2) Here ¯ P = ( P L + 2 P ⊥ ) is the average (isotropic) pressure and P eq = P eq ( E ) is theequilibrium pressure. .2. Dynamical variables The dynamical variables that we propagate in the code are q = ( T τµ , P L , P ⊥ , W µ ⊥ z , π µν ⊥ ) . (3)Their evolution equations will be discussed below. Although W µ ⊥ z and π µν ⊥ each have only two independent components, we evolve all 14 of their compo-nents independently to simplify the workflow of the algorithm. In addition,we propagate the energy density E and the fluid velocity’s spatial compo-nents u = ( u x , u y , u η ) since they appear in the hydrodynamic equations.They are inferred from the components T τµ : T τµ = E u τ u µ + P L z τ z µ − P ⊥ Ξ τµ + 2 W ( τ ⊥ z z µ ) + π τµ ⊥ . (4)To solve these equations, one also needs to know z µ . From the orthogonalityconditions z µ z µ = − z µ u µ = 0, there are only two nonzero componentsthat depend on the fluid velocity: z µ = 1 (cid:112) u ⊥ (cid:18) τ u η , , , u τ τ (cid:19) , (5)where u ⊥ = (cid:112) ( u x ) + ( u y ) is the transverse velocity. The solution of theinferred variables ( E , u ) from the algebraic equations (4) will be discussedin Sec. 3.3.For (3+1)–dimensionally expanding fluids, we have a total of 20 dynam-ical variables and four inferred variables to evolve on an Eulerian grid. , If the system is longitudinally boost-invariant, we do not need to propagatethe components T τη , W µ ⊥ z , π µη ⊥ and u η since they vanish by symmetry. When propagating these extraneous components numerically, slight violations of theorthogonality and tracelessness conditions (74) can occur. We correct for these errors ateach step of the simulation with the regulation scheme described in Sec. 3.5. The fluid velocity’s temporal component is u τ = (cid:112) u x ) + ( u y ) + ( τ u η ) . Dynamical variables are evolved directly using the Kurganov–Tadmor algorithm (seeSec. 3.1). Inferred variables are determined from the dynamical variables algebraically. For anisotropic fluid dynamics with a QCD equation of state [66], we further evolvethe mean-field B as a dynamical variable and the anisotropic variables (Λ, α L , α ⊥ ) asadditional inferred variables (see Secs. 2.4 and 2.6.3). .3. Conservation laws The evolution of the components T τµ is given by the energy-momentumconservation laws D µ T µν = 0 , (6)where the covariant derivative D µ accounts for the curvilinear nature of theMilne coordinates. The conservation equations can be expanded as ∂ µ T µν + Γ µµλ T λν + Γ νµλ T λµ = 0 , (7)where ∂ µ is the partial derivative and Γ µνλ are the Christoffel symbols. InMilne spacetime, the only nonzero Christoffel symbols areΓ τηη = τ, Γ ητη = Γ ηητ = 1 τ . (8)Thus, the set of equations (7) can be rewritten as ∂ τ T ττ + ∂ i T τi = − T ττ + τ T ηη τ , (9a) ∂ τ T τx + ∂ j T xj = − T τx τ , (9b) ∂ τ T τy + ∂ j T yj = − T τy τ , (9c) ∂ τ T τη + ∂ j T ηj = − T τη τ , (9d)where the Latin indices ( i, j ) ∈ ( x, y, η ) are summed over spatial components.We can eliminate the components T τi in (9a) by using the identity: T τi = ( E + P ⊥ ) u τ u i + L τi + W τi + π τi ⊥ = v i T ττ + v i ( P ⊥ −L ττ −W ττ − π ττ ⊥ ) + L τi + W τi + π τi ⊥ . (10)Here we introduced the three-velocity v i = u i /u τ as well as the tensors L µν =∆ P z µ z ν and W µν = 2 W ( µ ⊥ z z ν ) . Likewise, we can express the components T ij in (9b-d) in terms of T τi : T ij = ( E + P ⊥ ) u i u j − P ⊥ g ij + L ij + W ij + π ij ⊥ = v j T τi − v j ( L τi + W τi + π τi ⊥ ) − P ⊥ g ij + L ij + W ij + π ij ⊥ . (11)12fter some algebra one obtains [49] ∂ τ T ττ + ∂ i ( v i T ττ ) = − T ττ + τ T ηη τ + ( L ττ + W ττ + π ττ ⊥ −P ⊥ ) ∂ i v i (12a)+ v i ∂ i ( L ττ + W ττ + π ττ ⊥ −P ⊥ ) − ∂ η L τη − ∂ i W τi − ∂ i π τi ⊥ ,∂ τ T τx + ∂ i ( v i T τx ) = − T τx τ − ∂ x P ⊥ + ( W τx + π τx ⊥ ) ∂ i v i (12b)+ v i ∂ i ( W τx + π τx ⊥ ) − ∂ η W xη − ∂ i π xi ⊥ ,∂ τ T τy + ∂ i ( v i T τy ) = − T τy τ − ∂ y P ⊥ + ( W τy + π τy ⊥ ) ∂ i v i (12c)+ v i ∂ i ( W τy + π τy ⊥ ) − ∂ η W yη − ∂ i π yi ⊥ ,∂ τ T τη + ∂ i ( v i T τη ) = − T τη τ − ∂ η P ⊥ τ + ( L τη + W τη + π τη ⊥ ) ∂ i v i (12d)+ v i ∂ i ( L τη + W τη + π τη ⊥ ) − ∂ η L ηη − ∂ i W ηi − ∂ i π ηi ⊥ . For the numerical algorithm the hydrodynamic equations must be writtenin conservative flux form [75]: ∂ τ q ( x ) + ∂ i F i ( x ) = S ( τ, q ( x ) , u ( x ) , E ( x ) , ∂ m q ( x ) , ∂ µ u ( x )) . (13)Here F i = v i q are the currents, S are the source terms and m ∈ ( x, y, η ) isa spatial index. Naturally, the evolution equations for T τµ already assumethis form. In the next subsection we will see that the relaxation equationsfor the dissipative flows require additional manipulations. The relaxation equations for the dissipative flows P L , P ⊥ , W µ ⊥ z and π µν ⊥ are [65, 66]˙ P L = P eq − ¯ P τ Π − P L −P ⊥ τ π / ζ Lz θ L + ¯ ζ L ⊥ θ ⊥ − W µ ⊥ z ˙ z µ (14a) For a conformal fluid, the energy density’s spatial derivatives are needed to evaluatethe source term ∂ i P ⊥ = ( ∂ i E − ∂ i P L ).
13 ¯ λ LW u W µ ⊥ z D z u µ + ¯ λ LW ⊥ W µ ⊥ z z ν ∇ ⊥ ,µ u ν − ¯ λ Lπ π µν ⊥ σ ⊥ ,µν , ˙ P ⊥ = P eq − ¯ P τ Π + P L −P ⊥ τ π + ¯ ζ ⊥ z θ L + ¯ ζ ⊥⊥ θ ⊥ + W µ ⊥ z ˙ z µ (14b)+ ¯ λ ⊥ W u W µ ⊥ z D z u µ − ¯ λ ⊥ W ⊥ W µ ⊥ z z ν ∇ ⊥ ,µ u ν + ¯ λ ⊥ π π µν ⊥ σ ⊥ ,µν , ˙ W { µ }⊥ z = − W µ ⊥ z τ π + 2¯ η Wu Ξ µν D z u ν − η W ⊥ z ν ∇ µ ⊥ u ν − (cid:0) ¯ τ Wz Ξ µν + π µν ⊥ (cid:1) ˙ z ν (14c) − ¯ λ WW u W µ ⊥ z θ L + ¯ δ WW W µ ⊥ z θ ⊥ + ¯ λ WW ⊥ σ µν ⊥ W ⊥ z,ν + ω µν ⊥ W ⊥ z,ν + ¯ λ Wπu π µν ⊥ D z u ν − ¯ λ Wπ ⊥ π µν ⊥ z λ ∇ ⊥ ,ν u λ , ˙ π { µν }⊥ = − π µν ⊥ τ π + 2¯ η ⊥ σ µν ⊥ − W { µ ⊥ z ˙ z ν } + ¯ λ ππ π µν ⊥ θ L − ¯ δ ππ π µν ⊥ θ ⊥ (14d) − ¯ τ ππ π λ { µ ⊥ σ ν }⊥ ,λ + 2 π λ { µ ⊥ ω ν }⊥ ,λ − ¯ λ πW u W { µ ⊥ z D z u ν } . Here a dot above any quantity denotes the co-moving time derivative u γ D γ ,and curly brackets denote either the transverse projection of a vector, t { µ } =Ξ µα t α , or the traceless double transverse projection of a rank-2 tensor, t { µν } =Ξ µναβ t αβ . We also define the longitudinal and transverse expansion rates θ L = z µ D z u µ and θ ⊥ = ∇ ⊥ ,µ u µ , where D z = − z ν D ν is the LRF longitu-dinal derivative and ∇ µ ⊥ = Ξ µν D ν is the transverse gradient. The trans-verse velocity-shear tensor is σ µν ⊥ = Ξ µναβ D ( α u β ) , while the transverse vor-ticity tensor is given by ω µν ⊥ = Ξ µα Ξ νβ D [ α u β ] , where the square brackets D [ α u β ] = ( D α u β − D β u α ) denote anti-symmetrization. The shear and bulkviscosity to entropy density ratios η/ S and ζ/ S and their associated re-laxation times τ π and τ Π , along with the anisotropic transport coefficientscoupled to the gradient forces, will be discussed in Sec. 2.6.We recast the relaxation equations in conservative flux form by using theproduct rule identities˙ W { µ }⊥ z = Ξ µα u γ D γ W α ⊥ z = u γ D γ W µ ⊥ z − W α ⊥ z u γ D γ Ξ µα , (15a)˙ π { µν }⊥ = Ξ µναβ u γ D γ π αβ ⊥ = u γ D γ π µν ⊥ − π αβ ⊥ u γ D γ Ξ µναβ (15b)to rewrite the l.h.s. of Eqs. (14a-d) as˙ P L = u γ ∂ γ P L , (16a)14 P ⊥ = u γ ∂ γ P ⊥ , (16b)˙ W { µ }⊥ z = u γ ∂ γ W µ ⊥ z + u γ Γ µγλ W λ ⊥ z + W α ⊥ z ( u µ a α − z µ ˙ z α ) , (16c)˙ π { µν }⊥ = u γ ∂ γ π µν ⊥ + u γ Γ µγλ π νλ ⊥ + u γ Γ νγλ π µλ ⊥ (16d)+ π µα ⊥ ( u ν a α − z ν ˙ z α ) + π να ⊥ ( u µ a α − z µ ˙ z α ) , where a µ = ˙ u µ is the fluid acceleration. Finally, we use the product ruleidentity u γ ∂ γ P L, ⊥ = u τ (cid:104) ∂ τ P L, ⊥ + ∂ i ( v i P L, ⊥ ) − P L, ⊥ ∂ i v i (cid:105) (17)to rewrite Eqs. (16a-b) as ∂ τ P L + ∂ i ( v i P L ) = P L ∂ i v i + 1 u τ (cid:34) P eq − ¯ P τ Π − P L −P ⊥ τ π / I L (cid:35) , (18a) ∂ τ P ⊥ + ∂ i ( v i P ⊥ ) = P ⊥ ∂ i v i + 1 u τ (cid:34) P eq − ¯ P τ Π + P L −P ⊥ τ π + I ⊥ (cid:35) ; (18b)here I L = ¯ ζ Lz θ L + ¯ ζ L ⊥ θ ⊥ − W µ ⊥ z ˙ z µ + ¯ λ LW u W µ ⊥ z D z u µ + ¯ λ LW ⊥ W µ ⊥ z z ν ∇ ⊥ ,µ u ν (19a) − ¯ λ Lπ π µν ⊥ σ ⊥ ,µν , I ⊥ = ¯ ζ ⊥ z θ L + ¯ ζ ⊥⊥ θ ⊥ + W µ ⊥ z ˙ z µ + ¯ λ ⊥ W u W µ ⊥ z D z u µ − ¯ λ ⊥ W ⊥ W µ ⊥ z z ν ∇ ⊥ ,µ u ν (19b)+ ¯ λ ⊥ π π µν ⊥ σ ⊥ ,µν are the gradient source terms for P L and P ⊥ . Similarly Eqs. (16c-d) can berewritten as ∂ τ W µ ⊥ z + ∂ i ( v i W µ ⊥ z ) = W µ ⊥ z ∂ i v i + 1 u τ (cid:20) − W µ ⊥ z τ π + I µW − P µW − G µW (cid:21) , (20a) ∂ τ π µν ⊥ + ∂ i ( v i π µν ⊥ ) = π µν ⊥ ∂ i v i + 1 u τ (cid:20) − π µν ⊥ τ π + I µνπ − P µνπ − G µνπ (cid:21) ; (20b)here I µW = Ξ µν (cid:0) η Wu D z u ν − ¯ τ Wz ˙ z ν (cid:1) − η W ⊥ z ν ∇ µ ⊥ u ν − π µν ⊥ ˙ z ν − ¯ λ WW u W µ ⊥ z θ L (21a)+ ¯ δ WW W µ ⊥ z θ ⊥ + ¯ λ WW ⊥ σ µν ⊥ W ⊥ z,ν + ω µν ⊥ W ⊥ z,ν + ¯ λ Wπu π µν ⊥ D z u ν ¯ λ Wπ ⊥ π µν ⊥ z λ ∇ ⊥ ,ν u λ , I µνπ = Ξ µναβ (cid:0) π λ ( α ⊥ ω β ) ⊥ ,λ − ¯ τ ππ π λ ( α ⊥ σ β ) ⊥ ,λ − W ( α ⊥ z ˙ z β ) − ¯ λ πW u W ( α ⊥ z D z u β ) (cid:1) (21b)+ 2¯ η ⊥ σ µν ⊥ + ¯ λ ππ π µν ⊥ θ L − ¯ δ ππ π µν ⊥ θ ⊥ are the gradient source terms, P µW = W α ⊥ z ( u µ a α − z µ ˙ z α ) , (22a) P µνπ = π µα ⊥ ( u ν a α − z ν ˙ z α ) + π να ⊥ ( u µ a α − z µ ˙ z α ) (22b)are the transverse projection source terms, and G µW = u γ Γ µγλ W λ ⊥ z , (23a) G µνπ = u γ Γ µγλ π νλ ⊥ + u γ Γ νγλ π µλ ⊥ (23b)are the geometric source terms for W µ ⊥ z and π µν ⊥ . The individual componentsthat make up the source terms in the relaxation equations (18a-b) and (20a-b) are listed in Appendices A and B.If we evolve anisotropic fluid dynamics with a QCD equation of state (seeSecs. 2.5 and 2.6) we also propagate a mean field B as a dynamical variable.Its relaxation equation is [66, 79]˙ B = B eq − Bτ Π − ˙ mm ( E − P ⊥ −P L − B ) (24)or, in conservative flux form, ∂ τ B + ∂ i ( v i B ) = B ∂ i v i + 1 u τ (cid:20) B eq − Bτ Π − ˙ mm ( E − P ⊥ −P L − B ) (cid:21) , (25)where B eq ( E ) is the equilibrium mean field and m ( E ) is the quasiparticlemass. In the code there are two options for the quark-gluon plasma’s equationof state P eq ( E ): conformal and QCD. The former assumes a non-interactinggas of massless quarks and gluons: P eq = E gT π , (26)16 .0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.0 ℰℰ SB P eq P SB ( a ) c s ( b ) T ( GeV ) Figure 1: (Color online)
Left:
The QCD energy density (blue) and equilibrium pressure(red), normalized by their Stefan–Boltzmann limits (26).
Right:
The squared speed ofsound (purple) as a function of temperature; the gray line indicates the conformal limit c s = . where T is the temperature and the degeneracy factor is g = π (cid:20) (cid:0) N c − (cid:1) + 72 N c N f (cid:21) , (27)with N c = 3 colors and N f = 3 massless quark flavors (i.e. up, down andstrange). This equation of state is primarily used to test the fluid dynam-ical simulation subject to either conformal Bjorken expansion or conformalGubser expansion (see Secs. 4.1 and 4.2).The QCD equation of state that we employ for more realistic simulationsinterpolates between the lattice QCD calculations provided by the HotQCDcollaboration [80] and a hadron resonance gas composed of the hadrons thatcan be propagated in the hadronic afterburner code SMASH [81]. Fig. 1shows the energy density, equilibrium pressure and speed of sound as a func-tion of temperature. The data table in the code covers the temperature range T ∈ [0 . , .
0] GeV. In hybrid model simulations of heavy-ion collisions, the equilibrium equation of stateused in the hydrodynamic module must be consistent with that in the hadronic afterburner.Otherwise, serious violations of energy and momentum conservation can occur at thehadronization phase. .0 0.2 0.4 0.6 0.8 1.00.02.04.06.08.0 mT ( a ) - - B eq P eq ( b ) T ( GeV ) Figure 2: (Color online) The quasiparticle mass to temperature ratio (blue, left) andequilibrium mean field normalized to the QCD equilibrium pressure (orange, right) asfunctions of temperature (for conformal systems, m = 0 and B eq = 0). In this section we list the transport coefficients appearing in the relaxationequations in Sec. 2.4, for both QCD and conformal equations of state.At the present moment, the transport coefficients of QCD matter (espe-cially in the nonperturbative region T ∈ [0 . , .
5] GeV) are not explicitlyknown from first principles. Instead, we parametrize the shear and bulk vis-cosities ( η/ S )( T ) and ( ζ/ S )( T ) as a function of temperature. In this work,we use the best-fit parametrization models from the JETSCAPE collabora-tion [1, 2]. The relaxation times and anisotropic transport coefficients arecomputed with a quasiparticle kinetic theory model, whose equation of stateis fitted to the QCD one. The kinetic model contains quasiparticles witha temperature-dependent mass m ( T ) and an equilibrium mean-field B eq ( T );these are shown in Figure 2. The shear viscosity is modeled as a linear piecewise function with a kinkat temperature T η :( η/ S )( T ) = ( η/ S ) kink + ( T − T η ) (cid:0) a low Θ( T η − T ) + a high Θ( T − T η ) (cid:1) , (28)where ( η/ S ) kink is the value of η/ S at T η , a low and a high are the left and rightslopes, respectively, and Θ is the Heaviside step function. The bulk viscosity18 .0 0.2 0.4 0.6 0.8 1.00.00.10.20.30.4 η / ( a ) ζ / ( b ) T ( GeV ) Figure 3: (Color online) The temperature parametrization of ( η/ S )( T ) and ( ζ/ S )( T ) usedin this work. (For conformal systems, we set η/ S = 0 . ζ/ S = 0.) is parametrized as a skewed Cauchy distribution:( ζ/ S )( T ) = ( ζ/ S ) max Λ ζ ( T ) Λ ζ ( T ) + ( T − T ζ ) , (29)where ( ζ/ S ) max is the normalization factor, T ζ is the peak temperature andΛ ζ ( T ) = w ζ (cid:0) λ ζ sgn( T − T ζ ) (cid:1) , (30)with w ζ and λ ζ being the width and skewness parameters, respectively.The best-fit values for the viscosity parameters are ( η/ S ) kink = 0 . T η = 0 .
223 GeV, a low = − .
776 GeV − , a high = 0 .
37 GeV − , ( ζ/ S ) max =0 . T ζ = 0 .
12 GeV, w ζ = 0 .
072 GeV and λ ζ = − .
122 (see Table II inRef. [2]). The resulting temperature dependence of η/ S and ζ/ S is shown inFigure 3.For conformal systems, we fix the shear viscosity to η/ S = 0 . ζ/ S = 0. 19 .0 0.2 0.4 0.6 0.8 1.00.00.51.01.52.02.53.0 τ π T ( a ) τ Π T ( b ) T ( GeV ) Figure 4: (Color online) The dimensionless shear and bulk relaxation times computedwith the quasiparticle kinetic model (solid color) and small-mass approximation (dashedcolor). (For conformal systems, τ π T = 5 η/ S and τ Π T = 0.) In the quasiparticle kinetic model [66, 79], the shear and bulk relaxationtimes are proportional to the shear and bulk viscosity, respectively, τ π = ηβ π , τ Π = ζβ Π , (31)where the viscosity to relaxation time ratios β π and β Π are β π ( T ) = I T , (32a) β Π ( T ) = 53 β π + c s (cid:0) m dmdT I − ( E + P eq ) (cid:1) , (32b)and c s and P eq are evaluated with the QCD equation of state. Here, we havedefined the thermodynamic integrals I nq = g (cid:90) P ( u · p ) n − q ( − p · ∆ · p ) q f eq (2 q +1)!! , (33) They correspond to a hybrid model whose particlization phase uses the 14-momentapproximation for the δf correction in the Cooper-Frye formula [1, 2, 74] (the viscosityparameters used for the code validation tests in Sec. 4 are slightly different because theywere taken from an early draft of Ref. [2]). (cid:82) P · · · ≡ (cid:82) d p/E p . . . indicates integration with the Lorentz-invariantmomentum space measure, p µ is the quasiparticle momentum, − p · ∆ · p isthe square of its spatial LRF momentum, u · p = (cid:112) m ( T ) − p · ∆ · p is itsLRF energy, and f eq = exp (cid:2) − u · p/T (cid:3) is the local-equilibrium distributionfunction.The normalized quasiparticle relaxation times τ π T and τ Π T are shown inFigure 4. These are compared to the relaxation times in standard viscoushydrodynamic models, where the kinetic transport coefficients are computedin the small-mass approximation ¯ m = m/T (cid:28) dm/dT = 0: τ π T ≈ η/ S + O (cid:0) ¯ m (cid:1) , (34a) τ Π T ≈ ζ T E + P eq ) (cid:0) − c s (cid:1) + O (cid:0) ¯ m (cid:1) . (34b)One sees that the shear relaxation times are very similar to each other, exceptat low temperatures T < . T < . T ∼ m = 0, dm/dT = 0), the shear relaxationtime is τ π = 5 η/ ( S T ) and the bulk relaxation time is τ Π = 0. Finally, we list the anisotropic transport coefficients [65, 66] that arecoupled to the gradient source terms in the relaxation equations for the lon-gitudinal pressure P L ,¯ ζ Lz = I − P L + B ) + m dmd E ( E + P L ) I , (35a)¯ ζ L ⊥ = I − P L − B + m dmd E ( E + P ⊥ ) I , (35b)¯ λ LW u = I I + m dmd E I , (35c)¯ λ LW ⊥ = 1 − ¯ λ LW u , (35d)¯ λ Lπ = I I + m dmd E I , (35e)21he transverse pressure P ⊥ ,¯ ζ ⊥ z = I − P ⊥ − B + m dmd E ( E + P L ) I , (36a)¯ ζ ⊥⊥ = 2( I −P ⊥ − B ) + m dmd E ( E + P ⊥ ) I , (36b)¯ λ ⊥ W ⊥ = 2 I I + m dmd E I , (36c)¯ λ ⊥ W u = ¯ λ ⊥ W ⊥ − , (36d)¯ λ ⊥ π = 1 − I I − m dmd E I , (36e)the longitudinal momentum diffusion current W µ ⊥ z ,¯ η Wu = 12 (cid:0) P L + B − I ) , (37a)¯ η W ⊥ = 12 ( P ⊥ + B − I ) , (37b)¯ τ Wz = P L − P ⊥ , (37c)¯ δ WW = ¯ λ WW ⊥ −
12 + m dmd E ( E + P ⊥ ) (cid:18) I I (cid:19) , (37d)¯ λ WW u = 2 − I I − m dmd E ( E + P L ) (cid:18) I I (cid:19) , (37e)¯ λ WW ⊥ = 2 I I − , (37f)¯ λ Wπu = I I , (37g)¯ λ Wπ ⊥ = ¯ λ Wπu − , (37h)and the transverse shear stress tensor π µν ⊥ :¯ η ⊥ = P ( k ) ⊥ − I , (38a)¯ δ ππ = 34 ¯ τ ππ + 12 − m dmd E ( E + P ⊥ ) (cid:18) I I (cid:19) , (38b)¯ τ ππ = 2 − I I , (38c)¯ λ ππ = ¯ λ Wπu − m dmd E ( E + P L ) (cid:18) I I (cid:19) , (38d)22 λ πW u = ¯ λ WW ⊥ − , (38e)¯ λ πW ⊥ = ¯ λ πW u + 2 . (38f)We define the anisotropic integrals I nrqs = g (2 π ) (cid:90) P ( u · p ) n − r − q ( − z · p ) r ( − p · Ξ · p ) q E sa f a (2 q )!! , (39)where E a = (cid:113) m ( T ) − ( p · Ξ · p ) /α ⊥ + ( − z · p ) /α L (40)and f a = e − E a / Λ is the leading-order anisotropic distribution function.In addition to the effective temperature Λ this anisotropic distribution f a depends (through E a ) on two momentum anisotropy parameters α L and α ⊥ ,which deform the longitudinal and transverse momentum space, respectively.In order to compute the non-conformal anisotropic transport coefficients ateach time step of the simulation, we propagate (Λ , α ⊥ , α L ) as inferred vari-ables. They are obtained from the kinetic part of E , P L and P ⊥ in thequasiparticle kinetic model: E ( k ) = E −
B , (41a) P ( k ) L = P L + B , (41b) P ( k ) ⊥ = P ⊥ + B , (41c)where E ( k ) = I is the kinetic contribution to the energy density (i.e. thetotal energy density minus the mean field contribution [66]), P ( k ) L = I isthe kinetic longitudinal pressure, and P ( k ) ⊥ = I is the kinetic transversepressure. The numerical method which we use to solve these equations willbe discussed in Sec. 3.4.The anisotropic transport coefficients in the conformal limit are listed inAppendix C.
3. Numerical scheme
In this section we discuss the numerical implementation of the hydro-dynamic equations in the code. The dynamical and inferred variables areevolved on an ( N x +4) × ( N y +4) × ( N η +4) Eulerian grid, where N x , N y and23 η are the number of physical grid points along each spatial direction. Agrid point with cell index ( i, j, k ) that corresponds to the lower left frontcorner of a fluid cell, has a spatial position x i = (cid:2) i − −
12 ( N x − (cid:3) ∆ x , (42a) y j = (cid:2) j − −
12 ( N y − (cid:3) ∆ y , (42b) η s,k = (cid:2) k − −
12 ( N η − (cid:3) ∆ η s , (42c)where ∆ x , ∆ y and ∆ η s are the lattice spacings. Physical fluid cells haveindices i ∈ [2 , N x +1], j ∈ [2 , N y +1] and k ∈ [2 , N η +1]. The numericalalgorithm also requires six sets of ghost cells with depth two, which neighborthe physical grid’s faces (for an illustration, see Fig. 3 in Ref. [20]). Theboundary conditions that we impose on the ghost cells neighboring the two( y, η s ) faces are q ,j,k = q ,j,k = q ,j,k , (43a) u ,j,k = u ,j,k = u ,j,k , (43b) q N x +2 ,j,k = q N x +3 ,j,k = q N x +1 ,j,k , (43c) u N x +2 ,j,k = u N x +3 ,j,k = u N x +1 ,j,k , (43d)and similarly for ( x, η s ) and ( x, y ) faces after permuting the grid indices andreplacing N x → N y or N η .The dynamical variables q are updated using a two-stage Runge–Kutta(RK2) scheme, where the time derivatives are evaluated with the Kurganov–Tadmor (KT) algorithm [16, 20, 75]. After each intermediate Euler step inthe RK2 scheme, we reconstruct the inferred variables from the dynamicalvariables. In addition, we regulate the residual shear stresses W µ ⊥ z and π µν ⊥ and the mean-field B . The code also provides the user with the option ofusing an adaptive time step to capture the fluid’s longitudinal dynamics atvery early times; this will be discussed at the end of the section. For longitudinally boost-invariant systems, the number of spacetime rapidity pointsis set to N η = 1. In conformal anisotropic hydrodynamics, the energy density E also requires ghost cellboundary conditions. .1. Kurganov-Tadmor algorithm The time derivative of the dynamical variables ∂ τ q in the partial differ-ential equations (13) can be computed on the Eulerian grid at any time τ using the KT algorithm [75]:( ∂ τ q ) ijk = − H xi + ,j,k − H xi − ,j,k ∆ x − H yi,j + ,k − H yi,j − ,k ∆ y − H ηi,j,k + − H ηi,j,k − ∆ η s + S ijk (cid:0) τ, q ijk , u ijk , E ijk , ( ∂ m q ) ijk , ( ∂ µ u ) ijk (cid:1) , (44)where the numerical fluxes evaluated at the left and right faces of a staggeredcell centered around the grid point ( i, j, k ) are H xi ± ,j,k = 12 (cid:104) F x + i ± ,j,k + F x − i ± ,j,k − s xi ± ,j,k (cid:0) q + i ± ,j,k − q − i ± ,j,k (cid:1)(cid:105) , (45)and similarly for H yi,j ± ,k and H ηi,j,k ± after permuting the ± in the gridindices and the corresponding spatial component (i.e. x → y or η ).The first two terms in Eq. (45) take the average of the currents extrap-olated to the staggered cell face ( i + , j, k ) (or ( i − , j, k )) from the left ( − )and right (+) sides. A first-order expression for the extrapolated currentscan be computed using the chain rule: F x − i + ,j,k = F xijk + ∆ x (cid:104) ( ∂ x v x ) ijk q ijk + v xijk ( ∂ x q ) ijk (cid:105) , (46a) F x + i + ,j,k = F xi +1 ,j,k − ∆ x (cid:104) ( ∂ x v x ) i +1 ,j,k q i +1 ,j,k + v xi +1 ,j,k ( ∂ x q ) i +1 ,j,k (cid:105) , (46b) F x − i − ,j,k = F xi − ,j,k + ∆ x (cid:104) ( ∂ x v x ) i − ,j,k q i − ,j,k + v xi − ,j,k ( ∂ x q ) i − ,j,k (cid:105) , (46c) F x + i − ,j,k = F xijk − ∆ x (cid:104) ( ∂ x v x ) ijk q ijk + v xijk ( ∂ x q ) ijk (cid:105) , (46d)where F xijk = v xijk q ijk . Similarly, the extrapolated currents F y + i,j ± ,k , F y − i,j ± ,k , F η + i,j,k ± and F η − i,j,k ± at the remaining faces of the staggered cell are obtainedby permuting the ± (or ±
1) in the grid indices, the spatial components andderivatives, and the lattice spacing.The final term in Eq. (45) takes into account the wave propagation of thediscontinuities q + i ± ,j,k − q − i ± ,j,k at a finite speed [75]. We define the local25ropagation speed component s xi ± ,j,k at the staggered cell faces ( i ± , j, k )as s xi ± ,j,k = max (cid:0) | v x − i ± ,j,k | , | v x + i ± ,j,k | (cid:1) , (47)where the extrapolated velocities are v x − i + ,j,k = v xijk + ∆ x ∂ x v x ) ijk , (48a) v x + i + ,j,k = v xi +1 ,j,k − ∆ x ∂ x v x ) i +1 ,j,k , (48b) v x − i − ,j,k = v xi − ,j,k + ∆ x ∂ x v x ) i − ,j,k , (48c) v x + i − ,j,k = v xijk − ∆ x ∂ x v x ) ijk . (48d)The discontinuities q + i ± ,j,k − q − i ± ,j,k propagating from the staggered cell faces( i ± , j, k ) depend on the extrapolated dynamical variables q − i + ,j,k = q ijk + ∆ x ∂ x q ) ijk , (49a) q + i + ,j,k = q i +1 ,j,k − ∆ x ∂ x q ) i +1 ,j,k , (49b) q − i − ,j,k = q i − ,j,k + ∆ x ∂ x q ) i − ,j,k , (49c) q + i − ,j,k = q ijk − ∆ x ∂ x q ) ijk . (49d)The formulae for the local propagation speed components s yi,j ± ,k and s ηi,j,k ± ,extrapolated velocities v y + i,j ± ,k , v y − i,j ± ,k , v η + i,j,k ± and v η − i,j,k ± , and extrapolateddynamical variables q + i,j ± ,k , q − i,j ± ,k , q + i,j,k ± and q − i,j,k ± are analogous.The numerical spatial derivatives appearing in the extrapolated quan-tities (46a-d), (48a-d) and (49a-d) are computed with a minmod flux lim-iter [75]:( ∂ x q ) ijk = M (cid:16) Θ q ijk − q i − ,j,k ∆ x , q i +1 ,j,k − q i − ,j,k x , Θ q i +1 ,j,k − q ijk ∆ x (cid:17) , (50a)( ∂ x v x ) ijk = M (cid:16) Θ v xijk − v xi − ,j,k ∆ x , v xi +1 ,j,k − v xi − ,j,k x , Θ v xi +1 ,j,k − v xijk ∆ x (cid:17) , (50b)where M ( a, b, c ) = minmod( a, minmod( b, c )) , (51)26ith minmod( a, b ) = sgn( a ) + sgn( b )2 × min( | a | , | b | ) ; (52)the flux limiter parameter is set to Θ = 1 . ∂ y q ) ijk , ( ∂ y v y ) ijk , ( ∂ η q ) ijk and ( ∂ η v η ) ijk are analogous.In contrast, the numerical spatial derivatives appearing in the sourceterms S ijk are approximated with second-order central differences [20, 76]:( ∂ x q ) ijk = q i +1 ,j,k − q i − ,j,k x , (53a)( ∂ x u ) ijk = u i +1 ,j,k − u i − ,j,k x , (53b)and similarly for ( ∂ y q ) ijk , ( ∂ y u ) ijk , ( ∂ η q ) ijk and ( ∂ η u ) ijk . The fluid velocity’stime derivative ( ∂ τ u ) ijk also appears in the source terms; its evaluation willbe discussed the next subsection. Because the KT algorithm (44) admits a semi-discrete form [75] (i.e. thetime derivative is continuous while the spatial derivatives are discrete), wecan combine it with an RK2 ODE solver to evolve the system in time [16, 20].Given the dynamical variables q n,ijk ≡ q ijk ( τ n ) at time τ = τ n (discrete timesare labeled with index n ), we evolve the system one time step ∆ τ n with anintermediate Euler step (omitting the spatial indices): q I ,n +1 = q n + ∆ τ n E ( τ n , q n , u n , E n ; u n − , ∆ τ n − ) , (54)where E = ∂ τ q is evaluated with r.h.s of Eq. (44). The time derivative ∂ τ u inthe source terms is approximated with a first-order backward difference [17,20]: ( ∂ τ u ) n = u n − u n − ∆ τ n − , (55)where u n − is the previous fluid velocity and ∆ τ n − is the previous timestep. From the intermediate variables (54), we reconstruct the inferredvariables ( E I ,n +1 , u I ,n +1 ) as well as the anisotropic variables (Λ I ,n +1 α ⊥ , I ,n +1 , At the start of the hydrodynamic simulation, n = 0 or τ = τ , the previous time step∆ τ n − is set to the current time step ∆ τ n . Unless stated otherwise, we also initialize theprevious fluid velocity as u n − = u n . L, I ,n +1 ), as described in Secs. 3.3 and 3.4 below. Afterwards, we regulate themean field B and residual shear stresses W µ ⊥ z and π µν ⊥ in q I ,n +1 (see Sec. 3.5)and set the ghost cell boundary conditions for q I ,n +1 and u I ,n +1 .Next, we evolve the system with a second intermediate Euler step Q n +2 = q I ,n +1 + ∆ τ n E ( τ n +∆ τ n , q I ,n +1 , u I ,n +1 , E I ,n +1 ; u n , ∆ τ n ) , (56)where the fluid velocity’s time derivative is now evaluated as( ∂ τ u ) I ,n +1 = u I ,n +1 − u n ∆ τ n . (57)In the RK2 scheme, we average the two intermediate Euler steps in Q n +2 toupdate the dynamical variables at τ = τ n + ∆ τ n : q n +1 = q n + Q n +2 . (58)From this, we update the inferred variables ( E n +1 , u n +1 ) and (Λ n +1 , α ⊥ ,n +1 , α L,n +1 ) and regulate the residual shear stresses and mean field. Finally, weset the ghost cell boundary conditions for q n +1 and u n +1 and proceed withthe next RK2 iteration. Given the hydrodynamic variables T τµ , along with P L , P ⊥ , W µ ⊥ z and π τµ ⊥ ,we can reconstruct the energy density and fluid velocity from Eq. (4). Thesolution for the energy density is [66] E = M τ − L ττ − ( M x ) +( M y ) M τ + P ⊥ −L ττ − ( τ M η ) ( M τ + P ⊥ −L ττ )( M τ + P L ) , (59)where L ττ = ∆ P ( z τ ) was defined earlier and M µ = T τµ − W ( τ ⊥ z z µ ) − π τµ ⊥ . (60)The reconstruction formula (59) also requires the components z τ , z η , whichcan be expressed in terms of hydrodynamic variables as follows: z τ = τ F (cid:112) − ( τ F ) , z η = 1 τ (cid:112) − ( τ F ) . (61) Writing F = u η /u τ implies that the argument 1 − ( τ F ) in Eq. (61) is always positive. F = A − B (cid:112) τ ( B − A )1 + ( τ B ) , (62)with A = K η / ( K τ + P L ) and B = W τ ⊥ z / ( τ ( K τ + P L )) where K µ = T τµ − π τµ ⊥ .In the cold, dilute regions surrounding the fireball (and, occasionally, incold spots within the fluctuating fireball), the energy density can becomemuch smaller than the freezeout energy density E sw . While these regionsare not phenomenologically important, hydrodynamic simulations are sus-ceptible to crashing there without intervention [5, 20, 83]. To prevent this,we regulate the energy density with the formula
E ← E + + E min e −E + / E min , (63)where E min is the minimum energy density allowed in the Eulerian grid and E + = max(0 , E ). For E (cid:29) E min , the regulation has virtually no effect on theenergy density. As
E →
0, however, the energy density is smoothly regulatedto E min . Ideally, E min should be the lower limit of our QCD equation ofstate table E low = 3 . × − GeV/fm . However, we find that the codeevolving non-conformal anisotropic hydrodynamics with fluctuating initialconditions encounters fewer technical difficulties if we instead use the largervalue E min = 0 .
02 GeV/fm , which is still about six times smaller than ourchoice for E sw . We checked that the regulation scheme (63) has little to noimpact on the fluid’s dynamics in regions where
E ∼ E sw or larger.After regulating the energy density, we reconstruct the fluid velocity’sspatial components: u x = M x (cid:112) ( E + P ⊥ )( M τ + P ⊥ − L ττ ) , (64a) u y = M y (cid:112) ( E + P ⊥ )( M τ + P ⊥ − L ττ ) , (64b) u η = F (cid:114) M τ + P ⊥ − L ττ E + P ⊥ . (64c) In this work, we construct a particlization hypersurface of constant energy density E sw = 0 .
116 GeV/fm , which corresponds to the switching temperature T sw = 0 .
136 GeV.The lowest value used for the switching temperature in the JETSCAPE SIMS analysis is T sw = 0 .
135 GeV [1, 2]. The constraint
E ≥ E min is not imposed for conformal systems. E , u ) and the components T τµ but only in the dilute coldregions. To compute the non-conformal anisotropic transport coefficients in Sec. 2.6.3,we need to reconstruct the anisotropic variables X = Λ α ⊥ α L (65)from E , P L , P ⊥ , B and m ( E ), by solving the system of equations [66] f ( X ) = , (66)where f ( X ) = I ( X ) − E + B I ( X ) − P L − B I ( X ) − P ⊥ − B . (67)We solve these nonlinear equations numerically using Newton’s method. Tak-ing the anisotropic variables prior to the intermediate Euler step (54) as ourinitial guess X g = (Λ n , α ⊥ ,n , α L,n ) T , we iterate X along the direction givenby the Newton step ∆ X = − J − f , (68)where J − is the inverse of the Jacobian J = ∂ f ∂ X = I ( X )Λ I − ( X )Λ α ⊥ I − ( X )Λ α L I ( X )Λ I − ( X )Λ α ⊥ I − ( X )Λ α L I ( X )Λ I − ( X )Λ α ⊥ I − ( X )Λ α L . (69) At τ = τ we initialize the anisotropic variables (Λ , α ⊥ , , α L, ) by iterating X g =( T, , T (where T inside the parentheses is the temperature). X ← X + λ ∆ X , (70)where λ ∈ [0 ,
1] is a partial step that is optimized with a line-backtrackingalgorithm to improve the global convergence of Newton’s method [84]. Thisprocedure is repeated until X has converged to the solution of Eq. (66). For conformal systems ( B = 0, m = 0, α ⊥ = 1), the anisotropic variablesΛ and α L are much easier to solve. The system of equations (67) reduces to E = 3 g Λ R ( α L )2 π , (71a) P L = g Λ R ( α L )2 π , (71b)where the functions R and R are listed in Appendix C. We numericallyinvert the longitudinal pressure to energy density ratio for α L : R ( α L ) R ( α L ) = P L E . (72)From this, we can evaluate Λ asΛ = (cid:18) π E g R ( α L ) (cid:19) / . (73)Because the solutions (72) and (73) do not require the previous values for Λand α L as input, we do not need to store them in memory during runtime. After reconstructing the inferred variables, we regulate the residual shearstress components W µ ⊥ z and π µν ⊥ such that they satisfy the orthogonality andtraceless conditions u µ W µ ⊥ z = 0 , (74a) z µ W µ ⊥ z = 0 , (74b) In our simulation tests, we found several instances when Newton’s method failed toconverge near the edges of the Eulerian grid. Whenever that happens we simply set theanisotropic variables to the initial guess X g . µ π µν ⊥ = 0 , (74c) z µ π µν ⊥ = 0 , (74d) π µ ⊥ ,µ = 0 , (74e)and that their overall magnitude is smaller than the longitudinal and trans-verse pressures: (cid:113) π ⊥ ,µν π µν ⊥ − W ⊥ z,µ W µ ⊥ z ≤ (cid:113) P L + 2 P ⊥ . (75)The former ensures that W µ ⊥ z and π µν ⊥ maintain their orthogonal and trace-lessness properties within numerical accuracy during the hydrodynamic sim-ulation. The latter prevents the residual shear stresses from overwhelmingthe anisotropic part of the energy-momentum tensor [49].First, we update these components in the following order: W τ ⊥ z ← u τ ( W x ⊥ z u x + W y ⊥ z u y )1 + u ⊥ , (76a) W η ⊥ z ← W τ ⊥ z u η u τ , (76b) π yy ⊥ ← π xy ⊥ u x u y − π xx ⊥ (cid:0) u y ) (cid:1) u x ) , (76c) π τx ⊥ ← u τ ( π xx ⊥ u x + π xy ⊥ u y )1 + u ⊥ , (76d) π τy ⊥ ← u τ ( π xy ⊥ u x + π yy ⊥ u y )1 + u ⊥ , (76e) π xη ⊥ ← π τx ⊥ u η u τ , (76f) π yη ⊥ ← π τy ⊥ u η u τ , (76g) π τη ⊥ ← u τ ( π xη ⊥ u x + π yη ⊥ u y )1 + u ⊥ , (76h) π ηη ⊥ ← π τη ⊥ u η u τ , (76i) π ττ ⊥ ← π τx ⊥ u x + π τy ⊥ u y + τ π τη ⊥ u η u τ , (76j)while leaving the components W x ⊥ z , W y ⊥ z , π xx ⊥ and π xy ⊥ unchanged. Next, we32escale all the components by the same factor ρ reg ∈ [0 , W µ ⊥ z ← ρ reg W µ ⊥ z , (77a) π µν ⊥ ← ρ reg π µν ⊥ , (77b)where ρ reg = min (cid:32) , (cid:115) P L + 2 P ⊥ π ⊥ · π ⊥ − W ⊥ z · W ⊥ z (cid:33) , (78)with π ⊥ · π ⊥ = π ⊥ ,µν π µν ⊥ . In the validation tests discussed in Sec. 4, we didnot encounter a situation where the residual shear stresses are suppressedby Eq. (78). This indicates that the residual shear stresses are naturallysmaller than the leading-order anisotropic pressures during the simulation.Nevertheless, we keep this procedure in place as a precaution.For hydrodynamic simulations with two or three spatial dimensions, wefind that the relaxation equation (24) for the mean-field B can become un-stable in certain spacetime regions ( T ∼ . − .
16 GeV), causing the bulkviscous pressure Π to grow positive without bound. The origin of this issue isnot fully understood, and we will address it in future work. For now, we re-move this instability by regulating the nonequilibrium mean-field component δB = B − B eq as δB ← κ reg δB, (79)where κ reg = min (cid:16) , − | B eq | δB (cid:17) ∀ δB < . (80)In practice, we only find it necessary to regulate the mean-field when δB < > Most relativistic hydrodynamic codes that simulate heavy-ion collisionsevolve the system with a fixed time step ∆ τ n = ∆ τ [5, 16]. Any choice for∆ τ must satisfy the CFL condition so that the hydrodynamic simulationis at least dynamically stable [75]. However, the time step must also besmall enough to resolve the fluid’s evolution rate (in particular, the largelongitudinal expansion rate θ L ∼ /τ at early times) but, on the other hand, For longitudinally boost-invariant systems, we replace the second argument in Eq. (78)by (cid:112) P ⊥ /π ⊥ · π ⊥ . τ ∼ . c [26]). This is not much of a concern for second-order viscoushydrodynamics since it is anyhow prone to breaking down for very earlyinitialization times, as discussed in Sec. 1. Anisotropic hydrodynamics, onthe other hand, can handle the large pressure anisotropies occurring at earlytimes much better; to realize its full potential it should be initialized at earliertimes, but for that we need to move away from a fixed time step. In thissection, we introduce a new adaptive stepsize method, which automaticallyadjusts the successive time step ∆ τ n +1 to be larger or smaller than ∆ τ n insuch a way that we can push back our fluid dynamical simulation to veryearly times ( τ ∼ . − .
05 fm /c ) without sacrificing numerical accuracynor computational efficiency. For a system with no source terms ( S = ) the KT algorithm is dy-namically stable as long as the time step satisfies the CFL condition [75]∆ τ n ≤ ∆ τ CFL = 18 min (cid:18) ∆ xs x max ( τ n ) , ∆ ys y max ( τ n ) , ∆ η s s η max ( τ n ) (cid:19) , (81)where s i max ( τ n ) ( i = x, y, η ) are the maximum local propagation speed com-ponents on the Eulerian grid at time τ n . In Milne spacetime, the quark-gluonplasma’s flow profile is not ultrarelativistic throughout most of its evolution(as long as one uses a QCD equation of state). Thus, the criterium (81) allowsone to maintain dynamical stability with an adaptive time step that is gen-erally larger than the fixed time step obtained by taking the limit s i max → τ n = 18 min (∆ x, ∆ y, ∆ η s ) . (82)For S (cid:54) = , however, the time step ∆ τ CFL from (81) is too coarse to resolvethe fluid’s gradients and relaxation rates at early times τ < . c whenthe flow profile is very nonrelativistic. Therefore, at early times when thesource terms are strongest, the adaptive time step should primarily dependon those source terms. The following implementation ensures that, as the This is also useful for constructing the particlization hypersurface in peripheral heavy-ion collisions or small collision systems (e.g. p+p), whose fireball lifetimes are not thatmuch longer than the hydrodynamization time τ hydro (see Sec. 1). E = T ττ and u = ). The system of differential equations (13) for thedynamical variables q ( τ ) reduces to ∂ τ q ( τ ) = S ( τ, q ) , (83)and we are given the initial values q n and time step ∆ τ n at time τ n . Weevolve the system one time step ∆ τ n using the RK2 scheme (58): q n +1 = q n + ∆ τ n (cid:0) S ( τ n , q n ) + S ( τ n +∆ τ n , q n +∆ τ n S ( τ n , q n ) (cid:1) . (84)For the next iteration, we determine how much we need to adjust the timestep ∆ τ n +1 . Adaptive stepsize methods do this by estimating the local trun-cation error of each time step [84]. In our method, we approximate the localtruncation error of the next intermediate Euler step q I ,n +2 = q n +1 + ∆ τ n +1 S ( τ n +∆ τ n , q n +1 ) , (85)which is (cid:15) n +2 = 12 (cid:13)(cid:13) ( ∂ τ q ) n +1 (cid:13)(cid:13) ∆ τ n +1 + O (∆ τ n +1 ) , (86)where (cid:107) ... (cid:107) denotes the (cid:96) -norm. If the second time derivative ( ∂ τ q ) n +1 at τ = τ n + ∆ τ n is known, we can set the local truncation error to the desirederror tolerance (after dropping higher-order terms)12 (cid:13)(cid:13) ( ∂ τ q ) n +1 (cid:13)(cid:13) ∆ τ n +1 = δ × max( N / q , (cid:13)(cid:13) q I ,n +2 (cid:13)(cid:13) ) , (87)where δ is the error tolerance parameter and N q is the number of dynamicalvariables. It is reasonable to use absolute or relative errors when (cid:13)(cid:13) q I ,n +2 (cid:13)(cid:13) is small or large, respectively [84]. In this work, we set the error toleranceparameter to δ = 0 . (cid:63) ): q (cid:63) I ,n +2 = q n +1 + ∆ τ n S ( τ n +∆ τ n , q n +1 ) . (88)This allows us to approximate ( ∂ τ q ) n +1 with central differences:( ∂ τ q ) n +1 = 2( q (cid:63) I ,n +2 − q n +1 + q n )∆ τ n + O (∆ τ n ) . (89)35ompared to the usual central difference formula there is an additional factorof 2 that accounts for the local truncation error present in q (cid:63) I ,n +2 . Further-more, the expression (89) is numerically accurate to O (∆ τ n ) rather than O (∆ τ n ). After substituting Eqs. (85) and (89) in Eq. (87), one has for thenext time step ∆ τ n +1 = max (cid:0) ∆ τ (abs) n +1 , ∆ τ (rel) n +1 (cid:1) , (90)where ∆ τ (abs) n +1 = ∆ τ n (cid:118)(cid:117)(cid:117)(cid:116) δ N / q (cid:13)(cid:13)(cid:13) q (cid:63) I ,n +2 − q n +1 + q n (cid:13)(cid:13)(cid:13) (91)and ∆ τ (rel) n +1 is the numerical solution to the algebraic equation N / q (cid:0) ∆ τ (rel) n +1 (cid:1) (cid:0) ∆ τ (abs) n +1 (cid:1) = (cid:113) (cid:107) q n +1 (cid:107) + 2 q n +1 · S n +1 ∆ τ (rel) n +1 + (cid:107) S n +1 (cid:107) (cid:0) ∆ τ (rel) n +1 (cid:1) , (92)with S n +1 = S ( τ n +∆ τ n , q n +1 ).For the general case without Bjorken symmetry, q ( x ) varies across thegrid. We then perform the calculation (90) (replacing S by E ) for all spatialgrid points and take the minimum value. Afterwards, we place safety boundsto prevent the time step from changing too rapidly:(1 − α ) ∆ τ n ≤ ∆ τ n +1 ≤ (1+ α )∆ τ n , (93)where we set the control parameter to α = 0 .
5. Finally, we impose the CFLbound (81) on ∆ τ n +1 . With the new time step at hand, we resume comput-ing the next RK2 iteration. Notice that there are no additional numericalevaluations of the flux and source terms in the adaptive RK2 scheme since wecan recompute the next intermediate Euler step q I ,n +2 simply by adjustingthe time step in q (cid:63) I ,n +2 . This allows our adaptive time step algorithm to bereadily integrated into our numerical scheme.When we start the hydrodynamic simulation, we initialize the time step∆ τ to be 20 times smaller than the initial time τ . At first, the adap-tive time step ∆ τ n tends to increase with time since the fluid’s longitudinalexpansion rate decreases. As the transverse flow builds up, ∆ τ n eventuallybecomes bounded by the CFL condition (81). We do not allow the adaptive time step to become any smaller than this (i.e. weimpose ∆ τ n ≥ . τ ). .7. Program summary We close this section by summarizing the workflow of the anisotropicfluid dynamical simulation with a QCD equation of state. Figure 5 showsthe flowchart of the code’s primary mode, which constructs a hypersurfaceof constant energy density E sw for a particlization module:1. We read in the runtime parameters and configure the Eulerian grid.2. We allocate memory to store the dynamical and inferred variables in theEulerian grid at a given time τ n . Altogether, there are three dynamicalvariables ( q , qI , Q ), two fluid velocity variables ( u , up ), one energydensity variable e and three anisotropic variables ( lambda , aT , aL ).They store one of the following variables during the RK2 iteration:(a) q holds the current or updated dynamical variables q n or q n +1 .(b) qI holds the intermediate dynamical variables q I ,n +1 (or q (cid:63) I ,n +1 ).(c) Q holds the previous, updated or current dynamical variables q n − , q n +1 or q n .(d) u holds the current, intermediate or updated fluid velocity u n , u I ,n +1 or u n +1 .(e) up holds the previous or current fluid velocity u n − or u n .(f) e holds the current, intermediate or updated energy density E n , E I ,n +1 or E n +1 .(g) lambda , aT and aL hold the current, intermediate or updatedanisotropic variables (Λ n , α ⊥ ,n , α L,n ), (Λ I ,n +1 , α ⊥ , I ,n +1 , α L, I ,n +1 )or (Λ n +1 , α ⊥ ,n +1 , α L,n +1 )3. We initialize the variables q , u , up , e , lambda , aT and aL as follows: first, we compute (or read in) the energy density profile given by aninitial-state model (e.g. T R ENTo [24]). Then we initialize the longi-tudinal and transverse pressures as P L = 3 R P eq R , (94a) P ⊥ = 3 P eq R , (94b) The secondary (test) mode outputs the hydrodynamic evolution from the simulationand their corresponding semi-analytic solutions, if they exist. If the code is run with Gubser initial conditions, the hydrodynamic variables areinitialized differently (see Sec. 4.2). tart programRead in parameters and configure grid No Allocate memory and initialize dynamicaland inferred variables. Set ghost cells.Start hydrodynamicevolution (n = 0) n mod 2 = 0
Compute adaptive time step Δτ n Find freezeoutcells min(T ijk ) < T sw No Evolve system one time step Δτ n Deallocate memory and savefreezeout surface
Yes Yes
End program n++
Figure 5: (Color online) Program flowchart of
VAH . R ∈ [0 ,
1] is a pressure anisotropy ratio parameter set by theuser (typically a small value R ≤ . P L − P ⊥ enters in the initial P L and P ⊥ whilethe initial bulk viscous pressure is Π = 0. The initial fluid velocity isstatic (i.e. u = up = ) and the residual shear stresses W µ ⊥ z and π µν ⊥ areset to zero. From this, we compute the components T τµ with Eq. (4).Finally, we set the mean-field to B = B eq and initialize the anisotropicvariables by solving Eq. (66).
4. We set the ghost cell boundary conditions for q , e and u . We also con-figure the freezeout finder and set the initial time step to ∆ τ = 0 . τ .5. We start the hydrodynamic evolution at τ = τ (or n = 0):(a) We call the freezeout finder every two time steps (i.e. n mod 2 =0) to search for freezeout cells between the Eulerian grids from thecurrent and previous calls (if n = 0, we only load the initial gridto the freezeout finder). The hypersurface volume elements d σ µ and their spacetime positions are constructed with the code COR-NELIUS [85]. The hydrodynamic variables at the hypersurfaceelements’ centroid are approximated with a 4d linear interpola-tion.(b) We compute the intermediate Euler step (54) and store the results q I ,n +1 in qI (if n >
0, we compute q (cid:63) I ,n +1 and the adaptive timestep ∆ τ n using the algorithm in Sec. 3.6 before evaluating q I ,n +1 ).Afterwards, we swap the variables u ↔ up so that up ← u n .(c) We reconstruct the intermediate inferred variables E I ,n +1 , u I ,n +1 ,Λ I ,n +1 , α ⊥ , I ,n +1 and α L, I ,n +1 from qI and store the results in e , u , lambda , aT and aL , respectively. We regulate the residual shearstresses and mean-field in qI and set the ghost cell boundary con-ditions for qI , e and u .(d) We compute the second intermediate Euler step (56) and updatethe dynamical variables q n +1 via Eq. (58), which is stored in Q . This initialization scheme is similar to how conformal free-streaming modules areinitialized at τ → The ghost cells of e are only used in conformal anisotropic hydrodynamics to approx-imate the source terms ∝ ∂ i E on the grid’s faces. The user can adjust the number of time steps between freezeout finder calls. q ↔ Q so that q ← q n +1 and Q ← q n .(e) We update the inferred variables E n +1 , u n +1 , Λ n +1 , α ⊥ ,n +1 and α L,n +1 from q and store the results in e , u , lambda , aT and aL ,respectively. We regulate the residual shear stresses and mean-field in q and set the ghost cell boundary conditions for q , e and u .(f) Steps (a) – (e) are repeated until the temperature of all fluid cellsin the grid are below the switching temperature (i.e. min( T ijk ) 4. Validation tests and hydrodynamic model comparisons In this section we test the validity of our anisotropic fluid dynamical sim-ulation for various initial-state configurations, using either the conformal orQCD equation of state. First, we run anisotropic hydrodynamics with con-formal Bjorken and Gubser initial conditions [86–88], whose semi-analyticsolutions can be computed accurately using a fourth-order Runge–KuttaODE solver [48, 69]. After these two validation tests, we compare (3+1)–dimensional conformal anisotropic hydrodynamics and second-order viscoushydrodynamics for a central Pb+Pb collision with a smooth T R ENTo initial condition [24].Next, we run non-conformal anisotropic hydrodynamics and viscous hy-drodynamics with Bjorken initial conditions and compare them to their semi-analytic solutions [66]. Finally, we study the differences between (3+1)–dimensional non-conformal anisotropic hydrodynamics and viscous hydrody-namics in central Pb+Pb collisions with smooth or fluctuating T R ENTo The VAH code can also run second-order viscous hydrodynamics, where the kinetictransport coefficients are computed with either quasiparticle masses m ( T ) (see Fig. 2a) [66,79] or light masses (i.e. m/T (cid:28) 1) [20, 67] (we use the same shear and bulk viscosities as inFig. 3). In this work, we label the former model as quasiparticle viscous hydrodynamics (or VH ) and the latter as standard viscous hydrodynamics (or VH VAH offersdefinitive advantages in reliability over standard viscous hydrodynamics. For the first test, we evolve a conformal plasma undergoing Bjorken ex-pansion [86]. In Bjorken flow, the system is longitudinally boost-invariantand homogeneous in the transverse plane. As a result, the fluid velocity u µ = (1 , , , 0) is static and the residual shear stresses W µ ⊥ z and π µν ⊥ van-ish. The only two degrees of freedom are the energy density E = T ττ andlongitudinal pressure P L (the transverse pressure is P ⊥ = ( E −P L )). Theanisotropic hydrodynamic equations simplify to [69] ∂ τ E = − E + P L τ , (95a) ∂ τ P L = E − P L τ π + ¯ ζ Lz τ , (95b)where the conformal transport coefficient ¯ ζ Lz is given in Eq. (35a). We startthe simulation at τ = 0 . 01 fm/ c with an initial temperature T = 1 . 05 GeV,pressure ratio R = 10 − and shear viscosity η/ S = 0 . 2. We evolve the systemwith an adaptive time step, initially set to ∆ τ = 5 × − fm/ c , until thetemperature drops below T sw = 0 . 136 GeV at τ f ∼ 20 fm/ c . Figure 6 showsthe evolution of the energy density normalized to its initial value E andthe pressure ratio P L / P ⊥ . At very early times, the system is approximatelyfree-streaming due to the rapid longitudinal expansion, resulting in a largeKnudsen number and causing the energy density to decrease like E ≈ E τ /τ .Over time, the longitudinal expansion rate θ L = 1 /τ decreases and the systemapproaches local equilibrium, i.e. E ∝ τ − / and P L / P ⊥ → VAH simulation to the semi-analytic solutionof (95), which uses a fixed time step ∆ τ = 5 × − fm/ c . The simulation isin excellent agreement with the semi-analytic solution, with numerical errorsstaying below 0 . δ in the adaptive time stepalgorithm. Figure 7 shows the evolution of the adaptive time step ∆ τ n forthe parameters δ = 0 . 004 and α = 0 . 5. One sees that the time step increases Bjorken flow can be simulated on either one fluid cell (i.e. N x = N y = N η = 1) or alarger grid with homogeneous initial conditions. - - ℰℰ ( a ) τ / τ ( τ / τ ) / τ = / c T = = - P L P ⟂ ( b ) VAHSemi - analytic η / = ℰ / ℰ exact ( P L / P ⟂ )/( P L / P ⟂ ) exact τ ( fm / c ) Figure 6: (Color online) Conformal Bjorken evolution of the normalized energy density E / E (a) and pressure ratio P L / P ⊥ (b) from the anisotropic hydrodynamic simulation(dashed red) and semi-analytic solution (transparent red). The gray lines in (a) show twodifferent power laws for comparison. The subpanels at the bottom show the ratio betweenthe numerical simulation and the semi-analytic solutions (solid red). linearly with time and is not bounded by the CFL condition (81) since theflow is stationary (i.e. ∆ τ CFL → ∞ ). As a result, the adaptive RK2 schemequickly reaches the switching temperature in 145 time steps, compared toabout 40000 steps for the semi-analytic method. In the next test, we evolve a conformal fluid subject to Gubser flow [87,88], which is longitudinally boost-invariant and azimuthally symmetric. Thefluid velocity profile takes the following analytic form: u τ = cosh κ , u x = sinh κ cos φ , u y = sinh κ sin φ , u η = 0 , (96)where κ ( τ, r ) = tanh − (cid:20) q τ r q ( τ + r ) (cid:21) , (97)42 .01 0.1 1.0 1010 - - τ ( fm / c ) Δ τ n ( f m / c ) VAHSemi - analytic δ = ( a ) α = Figure 7: (Color online) The evolution of the adaptive time step ∆ τ n (dashed red) in theconformal Bjorken simulation. The semi-analytic method (transparent red) uses a fixedtime step of ∆ τ n = 5 × − fm/ c . with r = (cid:112) x + y being the transverse radius, φ = tan − ( y/x ) the azimuthalangle, and q an inverse length scale parameter that determines the fireball’stransverse size R ⊥ ∼ /q at τ = 0 + .One can transform the polar Milne coordinates ˜ x µ = ( τ, r, φ, η s ) to the deSitter coordinates ˆ x µ = ( ρ, θ, φ, η s ) where ρ ( τ, r ) = − sinh − (cid:20) − q (cid:0) τ − r (cid:1) qτ (cid:21) , (98a) θ ( τ, r ) = tan − (cid:20) qr q ( τ − r ) (cid:21) . (98b)The line element in this de Sitter space possesses SO (3) q ⊗ SO (1 , ⊗ Z symmetry [87, 88]: d ˆ s = − dρ + cosh ρ (cid:0) dθ +sin θdφ (cid:1) + dη s , (99)which makes the fluid velocity ˆ u µ = (1 , , , 0) stationary. Hence, the hy- A hat denotes a quantity in the de Sitter space ˆ x µ = ( ρ, θ, φ, η s ). All hatted quantitiesare made unitless by multiplying with appropriate powers of τ . ρ , and the anisotropichydrodynamic equations reduce to [48] ∂ ρ ˆ E = ( ˆ P L − E ) tanh ρ , (100a) ∂ ρ ˆ P L = ˆ E − P L τ π − (cid:0) P L + ˆ¯ ζ Lz (cid:1) tanh ρ . (100b)The residual shear stresses ˆ W µ ⊥ z and ˆ π µν ⊥ vanish under Gubser symmetry [48].In the semi-analytic solution (100), we start the evolution at ρ = − . 2, whichcorresponds to a corner in a 14 fm × 14 fm transverse grid ( r = 7 √ τ = 0 . 01 fm/ c , with an initial temperature ˆ T = 0 . R = 10 − . We set the inverse fireball size to q = 1 . − andthe shear viscosity to η/ S = 0 . 2. We evolve the system with a fixed stepsize∆ ρ = 10 − until ρ f = 1 . 1, which corresponds to the center of our Euleriangrid ( r = 0 fm) at τ f = 3 . 01 fm/ c .Next, we map the semi-analytic solution for the de Sitter energy densityˆ E ( ρ ) and longitudinal pressure ˆ P L ( ρ ) back to Milne coordinates [89]: E ( τ, r ) = ˆ E (cid:0) ρ ( τ, r ) (cid:1) τ , P L ( τ, r ) = ˆ P L (cid:0) ρ ( τ, r ) (cid:1) τ . (101)We use this map to set the initial energy density and longitudinal pressureprofiles E ( τ , x, y ) and P L ( τ , x, y ); the initial temperature at the center ofthe fireball is T , center = 1 . 05 GeV. We set the initial transverse shear stressto π µν ⊥ = 0. Finally, we use Eq. (96) to initialize the current fluid velocity u ← u µ ( τ , x, y ) and previous fluid velocity up ← u µ ( τ − ∆ τ , x, y ), wherethe initial time step is ∆ τ = 5 × − fm/ c . We evolve the Gubser profile ona 14 fm × 14 fm transverse grid with a lattice spacing of ∆ x = ∆ y = 0 . 05 fm.Figure 8 shows the evolution of the energy density and transverse fluid ve-locity u ⊥ = (cid:112) ( u x ) + ( u y ) in the transverse plane at various time frames. At early times, the very hot and compact fireball expands rapidly along thelongitudinal direction and quickly cools down without much change to itstransverse shape. Over time, the transverse expansion overtakes the longi-tudinal expansion, pushing the medium radially outward as a ring of fire . In order to output the hydrodynamic quantities at specific times, we readjust theadaptive time step whenever we are close to these output times (this is not shown inFig. 10). igure 8: (Color online) Conformal Gubser evolution of the energy density E (GeV/fm )(left column) and transverse fluid rapidity u ⊥ (right column) in the anisotropic hydrody-namic simulation. - ℰ η / = τ = / c ρ = - T = R = - Δ x = Δ y = ( a ) - - - VAHSemi - analytic u x τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c ( b ) ℰ / ℰ exact u x / u exact x - ( c ) P L P ⟂ ( d ) π ⟂ · π ⟂ P ⟂ - - - - - ( P L / P ⟂ )/( P L / P ⟂ ) exact - - - - - x ( fm ) Figure 9: (Color online) Conformal Gubser evolution of (a) E (GeV/fm ), (b) u x , (c) P L / P ⊥ , and (d) √ π ⊥ · π ⊥ / ( P ⊥ √ 2) along the x –axis ( y = 0), given by the anisotropichydrodynamic simulation (dashed color) and semi-analytic solution (transparent color).The subpanels below panels (a-c) show the ratio between numerical simulation and semi-analytic solution (solid color). One sees that the energy density and transverse fluid velocity maintain theirazimuthal symmetry throughout the evolution. However, there are some nu-merical fluctuations around the lines y = ± x at later times, especially nearthe peak of u ⊥ . This is a consequence of using a Cartesian grid with a finitespatial resolution. 46 .01 0.1 1.010 - - τ ( fm / c ) Δ τ ( f m / c ) Δτ n Δτ CFL δ = α = Δ x = Δ y = ( a ) Figure 10: (Color online) The evolution of the adaptive time step ∆ τ n (dashed red) andCFL bound ∆ τ CFL (transparent red) in the conformal Gubser simulation. Figure 9 shows the evolution of the energy density E , fluid velocity com-ponent u x , pressure anisotropy ratio P L / P ⊥ and transverse shear inverseReynolds number √ π ⊥ · π ⊥ / ( P ⊥ √ 2) along the x –axis at multiple time frames.Overall, there is very good agreement between the simulation and semi-analytic solution, except near the local extrema of u x , where the numericalerrors are about 1 − π µν ⊥ should remain zero throughout thesimulation. However, the transverse shear inverse Reynolds number from the VAH simulation is nonzero (albeit small, (cid:46) τ n becomes bounded by the CFLcondition (81) at τ ∼ . c due to the increasing transverse expansionrate. As a result, the Gubser simulation finishes in 407 time steps, comparedto 480 steps if we had used the fixed time step ∆ τ = ∆ x/ 8. Although weonly gain a slight 1 . × speedup, the adaptive time step is able to resolvethe fluid’s early-time dynamics more accurately than the fixed time step (82)because it is initially independent of the lattice spacing. This property isespecially useful when simulating more realistic nuclear collisions, where the47attice spacing required to resolve the participant nucleons’ transverse energydeposition is several times coarser than the one used in this test (e.g. ∆ x =∆ y ∼ . − . Here we compare (3+1)–dimensional conformal anisotropic hydrodynam-ics VAH to conformal second-order viscous hydrodynamics VH (see Ap-pendix E), for a central Pb+Pb collision at zero impact parameter ( b = 0fm) with static (in Milne coordinates) and smooth initial conditions. Togenerate the latter we use an azimuthally symmetric T R ENTo energy den-sity profile averaged over 2000 fluctuating events [24] and extended alongthe η s –direction with a smooth rapidity plateau [76] (see Appendix D formore information on the T R ENTo energy deposition model for high-energynuclear collisions). The initial temperature at the center of the fireball is T , center = 1 . 05 GeV at a starting time of τ = 0 . 01 fm /c which is the same forboth types of hydrodynamic simulations. The fluid is evolved from an initialpressure ratio R = 10 − with a constant specific shear viscosity η/ S = 0 . T sw = 0 . 136 GeV, which hap-pens after τ f ∼ − c .Figure 11 shows the evolution of E , u x , P L / P ⊥ and √ π ⊥ · π ⊥ / ( P ⊥ √ x –axis ( y = η s = 0). Initially, the pressure anisotropy P L −P ⊥ insecond-order viscous hydrodynamics is so large that the longitudinal pressureturns negative, especially near the grid’s edges. In comparison, the pressureratio P L / P ⊥ in anisotropic hydrodynamics is, at early times, only moderatelylarger in the central fireball region but quite dramatically different near itstransverse edge, staying positive everywhere. The resulting differences inearly-time viscous heating create a disparity between the two models for thenormalization of the energy density which persists to late times even afterthe P L / P ⊥ ratios have converged. On the other hand, the transverse flows u x predicted by the two models are nearly identical even at very early times sincethey are primarily driven by the transverse pressure gradients ∂ x P ⊥ which, inthe smooth T R ENTo profile, are smaller at early times than the longitudinalgradients ∼ /τ . The transverse shear stress π µν ⊥ , which tends to counteractthe fluid’s transverse acceleration, is larger in viscous hydrodynamics thanin anisotropic hydrodynamics, especially along the edges of the fireball. The transverse shear stress π µν ⊥ is generally nonzero even for azimuthally symmetric - - ( a ) ℰ η / = τ = / c T = = - - - ( b ) u x VAHVH Δ x = Δ y = Δη s = τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c - - ( c ) P L P ⟂ - - ( d ) π ⟂ · π ⟂ P ⟂ x ( fm ) Figure 11: (Color online) The evolution of (a) E (GeV/fm ), (b) u x , (c) P L / P ⊥ and (d) √ π ⊥ · π ⊥ / ( P ⊥ √ 2) along the x –axis ( y = η s = 0), given by conformal anisotropic hydrody-namics ( VAH , dashed color) and second-order viscous hydrodynamics ( VH , transparentcolor), for the smooth (3+1)–dimensional T R ENTo initial condition described in the text. Overall, however, the hydrodynamic variables are not substantially differentalong the transverse directions.The longitudinal profile, however, evolves very differently in anisotropichydrodynamics ( VAH ) compared to second-order viscous hydrodynamics( VH ). This is shown in Fig. 12 where we plot the evolution of the dimen-sionless longitudinal velocity τ u η (as well as E and P L / P ⊥ ) along the η s –axis flow profiles as long as they do not possess Gubser symmetry. - - ( a ) ℰ η / = τ = / c T = = - - - ( b ) τ u η VAHVH Δ x = Δ y = Δη s = τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c - ( c ) P L P ⟂ - ( d ) - W ⟂ z · W ⟂ z P L + P ⟂ x ⟂ = ( ) fm η s Figure 12: (Color online) The evolution of (a) E (GeV/fm ), (b) τ u η , (c) P L / P ⊥ and(d) √ W ⊥ z · W ⊥ z / (cid:112) P L +2 P ⊥ along the η s –axis ( x = y = 0), computed with conformalanisotropic hydrodynamics ( VAH , dashed color) and second-order viscous hydrodynamics( VH , transparent color), for the smooth (3+1)–d T R ENTo initial condition described inthe text. Note that in panel (d) the η s –axis is shifted transversely to ( x, y ) = (7 . , 0) fm. ( x = y = 0). Similar to Fig. 11, the shear stress (2) in viscous hydrodynam-ics quickly overpowers the equilibrium pressure, making P L negative initially.Along the longitudinal direction this causes a strong reversal of the longitudi-nal flow τ u η at around | η s | ∼ 7. This unphysical feature (primarily driven bythe strong longitudinal expansion rate at early times) results in a longitudinalcontraction (in η s ) of the fireball near its forward and backward edges at thebeginning of the simulation. In anisotropic hydrodynamics, the longitudinalpressure remains always positive, allowing for a stronger longitudinal flow50 .01 0.1 1.010 - - τ ( fm / c ) Δ τ ( f m / c ) Δτ n Δτ CFL δ = α = Δη s = Δ x = Δ y = ( a ) Figure 13: (Color online) The evolution of the adaptive time step ∆ τ n (dashed color)and CFL bound ∆ τ CFL (transparent color) in conformal anisotropic hydrodynamics (red)and second-order viscous hydrodynamics (blue) for the smooth (3+1)–d T R ENTo initialcondition. profile whose signs are consistent with the gradients of the rapidity distri-bution of the energy density (142) (and thus of the thermal pressure). Thisgives anisotropic hydrodynamics a clear advantage over second-order viscoushydrodynamics when simulating the longitudinal dynamics of a heavy-ioncollision.In panel (d) of Fig. 12 we also plot the spacetime rapidity dependence ofthe inverse Reynolds number (cid:112) − W ⊥ z · W ⊥ z / ( P L + 2 P ⊥ ). The longitudinalmomentum diffusion current W µ ⊥ z is nonzero only in regions that have bothlongitudinal and transverse gradients. For this reason we shift the η s –axistransversely to ( x, y ) = (7 . , 0) fm, which corresponds to the mid-right re-gion of the grid. As expected, the momentum diffusion current is weakestaround mid-rapidity, η s ∼ 0, where the fireball profile is approximately lon-gitudinally boost-invariant, and strongest along the sloping edges of the ra-pidity plateau (142). However, this longitudinal edge region is also the placewhere its inverse Reynolds number differs most strongly between anisotropicand standard viscous hydrodynamics. Overall, W µ ⊥ z only makes up a tinyfraction of the total shear stress (2) in anisotropic hydrodynamics. The sit-51ation may change for initial conditions whose longitudinal and transversegradients are larger than the ones used in this test comparison.Finally, we plot in Fig. 13 the adaptive time step for each of the two hydro-dynamic simulations. One sees that in VH ∆ τ n stagnates until τ ∼ . 02 fm/ c (see footnote 29) while the one in VAH starts increasing immediately. Thisindicates that at early times viscous hydrodynamics generally has a fasterevolution rate and therefore requires a smaller initial time step compared toanisotropic hydrodynamics. Nevertheless, both adaptive time steps remainbelow their CFL bound (initially dominated by the longitudinal velocity u η )until τ ∼ . − c . Notice that the CFL bounds for VH and VAH becomevirtually identical since they have almost the same transverse velocity profile(see Fig. 11b). Now we turn to testing our non-conformal hydrodynamic simulation withthe QCD equation of state from Fig. 1. First, we run (3+1)–d anisotropichydrodynamics with Bjorken initial conditions and compare it to the semi-analytic solution of the equations [66] ∂ τ E = − E + P L τ , (102a) ∂ τ P L = P eq − ¯ P τ Π − P L −P ⊥ τ π / ζ Lz τ , (102b) ∂ τ P ⊥ = P eq − ¯ P τ Π + P L −P ⊥ τ π + ¯ ζ ⊥ z τ , (102c) ∂ τ B = B eq − Bτ Π + E + P L mτ dmd E ( E − P ⊥ −P L − B ) ; (102d)here we use the shear and bulk relaxation times (31)–(32), with the viscosi-ties given by Eqs. (28)–(29) (see Fig. 3). The non-conformal anisotropictransport coefficients ¯ ζ Lz and ¯ ζ ⊥ z are given by Eqs. (35a) and (36a), respec-tively. We start the simulation at τ = 0 . 05 fm/ c with initial temperature T = 0 . 718 GeV and initial pressure ratio R = 0 . 3, and evolve the systemuntil τ f ∼ 80 fm/ c when the temperature falls below T sw = 0 . 136 GeV (weplot results only up to τ = 20 fm/ c ). We also repeat this for quasiparticleviscous hydrodynamics ( VH ) and standard viscous hydrodynamics ( VH - ℰℰ ( a ) τ / τ ℰ ideal ℰ τ = / c T = = ( b ) P L P ⟂ VAHVHVH2 Semi - analytic - - - - - ( c ) P L - P ⟂ P eq - ( d ) Π P eq τ ( fm / c ) Figure 14: (Color online) Non-conformal Bjorken evolution of (a) E / E , (b) the pressureratio P L / P ⊥ , (c) the pressure anisotropy ( P L −P ⊥ ), and (d) the bulk viscous pressure Π(the latter two normalized to the equilibrium pressure P eq ), given by anisotropic hydrody-namics ( VAH , dashed red), quasiparticle viscous hydrodynamics ( VH , dashed blue) andstandard viscous hydrodynamics ( VH 2, dashed green), along with their semi-analytic solu-tions (transparent color). The lower gray curve in panel (a) is from an ideal hydrodynamiccalculation (i.e. P L = P ⊥ = P eq ( E ) and η/ S = ζ/ S = 0). (see footnote 35). Ideally, we would have preferred using smaller valuesfor τ and R to better match the longitudinally free-streaming initial con-dition used in the conformal hydrodynamic simulations (see Secs. 4.1–4.3 VH and VH P eq ( E ) as VAH but different transportcoefficients, see Appendix E and Ref. [66]. VAH simu-lation has greater difficulty reconstructing the anisotropic variables (Λ, α ⊥ , α L ), especially directly at initialization. This indicates that our anisotropichydrodynamic model [66] which integrates the QCD equation of state con-sistently with quasiparticle kinetic transport coefficients, has its limitations.Figures 14a,b show the evolution of the normalized energy density E / E and pressure ratio P L / P ⊥ from the three hydrodynamic simulations. We alsodisentangle the longitudinal and transverse pressures’ viscous components P L −P eq = ∆ P +Π and P ⊥ −P eq = − ∆ P +Π into the pressure anisotropy∆ P = P L −P ⊥ and bulk viscous pressure Π = ( P L +2 P ⊥ ) − P eq ; their evolu-tions relative to P eq are shown in Figs. 14c,d. Although our initial conditionfor the pressure ratio is somewhat ad hoc, it quickly reaches its minimumvalue at τ ∼ . P L / P ⊥ ratio istypically larger in anisotropic hydrodynamics compared to the two viscoushydrodynamic models (until τ ∼ c ). This is mainly due to differ-ences in the higher-order transport coefficients associated with the pressureanisotropy (i.e. beyond η/ S and τ π ) [62, 63, 66, 90]. Although in Bjorken flowthe bulk viscous pressure Π is much smaller than the pressure anisotropy ∆ P ,it evolves very differently in standard viscous hydrodynamics VH VH and anisotropic hydrodynamics VAH . Since the dimensionless bulk relaxation time τ Π T < VH NS = − ζ/τ at τ ∼ c . In contrast, the bulk relaxation time used in anisotropic hydro-dynamics VAH and quasiparticle viscous hydrodynamics VH is significantlylarger (solid curve in Fig. 4b). As a consequence, it takes a much longer timefor Π to relax to Π NS [66, 79].One also sees that the simulations agree very well with their correspond-ing semi-analytic solutions (continuous lines in transparent color). Althoughwe do not display them explicitly here, the numerical errors are small enoughthat we can unambiguously distinguish the different dynamics of the threehydrodynamic models in comparison studies discussed in the next subsec-tions. 54 .5. (3+1)–d non-conformal hydrodynamic models in central Pb+Pb colli-sions4.5.1. Smooth T R ENTo initial conditions Next we run non-conformal hydrodynamics for a central Pb+Pb collisionwith smooth, azimuthally symmetric T R ENTo initial conditions; we set theinitial time to τ = 0 . 05 fm/ c and the initial pressure ratio to R = 0 . 3. Theinitial energy density profile is almost identical to the one in Sec. 4.3, exceptthe normalization ∝ /τ is five times smaller, making the initial temperatureat the center of the fireball T , center = 0 . 718 GeV. We evolve the system until,at τ f ∼ − 16 fm/ c , all fluid cells are below T sw = 0 . 136 GeV.Figure 15 shows, for the first ∆ τ = 6 fm/ c of the simulation, the evolutionof E , u x , P L / P ⊥ , √ π ⊥ · π ⊥ /( P ⊥ √ ( P L − P ⊥ )/ P eq and Π / P eq given bythe three hydrodynamic models along the x –axis ( y = η s = 0). Here thefireball maintains higher temperatures for a longer duration than in conformalhydrodynamics since the QCD equilibrium pressure is much weaker than itsStefan–Boltzmann limit (26) (see Fig. 1a). But similar to Fig. 11a, the energydensity in anisotropic hydrodynamics cools down at a slightly faster ratecompared to viscous hydrodynamics. This is initially due to the moderatelylarger P L / P ⊥ ratio in the central fireball region, which closely follows theBjorken evolution in Fig. 14b. The viscous corrections to the longitudinalpressure increase as we move towards the edges of the fireball, but P L stillremains positive in anisotropic hydrodynamics. The longitudinal pressure instandard viscous hydrodynamics, however, is strongly negative at early times(and also, to a lesser extent, for the quasiparticle case), although it doesnot significantly alter the fluid’s transverse dynamics. The P L / P ⊥ ratiosstart to overlap at later times, mainly driven by the pressure anisotropies’convergence in Fig. 15e.Compared to Fig. 11b, the transverse flow u x in non-conformal hydrody-namics is significantly weaker due to the softer QCD equation of state. Thedip in the speed of sound at the quark-hadron phase transition (see Fig. 1b)also flattens the transverse velocity gradients near the edges of the fireball.This causes a significant decrease in the transverse shear stress relative tothe conformal case in Fig. 11d. Among the hydrodynamic models, we seethat VH u x around the edges of the fireball sinceit has the largest negative bulk viscous pressure there. This is due to itsability to converge to the Navier-Stokes solution shortly after the simulationbegins. In contrast, the bulk viscous pressure in VH hovers slightly above55 .010.11.010100 ( a ) ℰ sw = / fm ℰ τ = / c T = = - - ( b ) u x VAHVHVH2 Δ x = Δ y = Δη s = τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c ( c ) P L P ⟂ ( d ) π ⟂ · π ⟂ P ⟂ - - - - - - - ( e ) P L - P ⟂ P eq - - - - - - - ( f ) Π P eq x ( fm ) Figure 15: (Color online) The evolution of (a) E (GeV/fm ), (b) u x , (c) P L / P ⊥ , (d) √ π ⊥ · π ⊥ / ( P ⊥ √ ( P L −P ⊥ ) / P eq and (f) Π / P eq along the x –axis ( y = η s = 0), givenby non-conformal anisotropic hydrodynamics ( VAH , dashed color), quasiparticle viscoushydrodynamics ( VH , transparent color) and standard viscous hydrodynamics ( VH 2, dot-ted transparent color), for the smooth (3+1)–d T R ENTo initial condition. The minimumenergy density parameter is set to E min = 0 . 02 GeV/fm . τ ∼ c . Even then, it hardly catches up to the bulk viscous pressure curvesin VH VH is much weaker compared to VH 2, which leads to a stronger trans-verse flow. Finally, the VAH simulation has a bulk viscous pressure thatis qualitatively similar to VH (except at the edges of the grid) but slightlylags behind it. As a result, it generates the largest transverse flow out of thethree simulations. The study here illustrates how the evolution of the bulkviscous pressure and transverse fluid velocity strongly depends on the choiceof hydrodynamic model, specifically the temperature-dependent model forthe bulk relaxation time τ Π . This is likely to have important implications forthe phenomenological extraction of the bulk viscosity ζ/ S [1, 2].Next, we compare the fireballs’ longitudinal evolution along the η s –axis( x = y = 0) in Figure 16. When comparing the longitudinal evolution of theenergy density and viscous pressures between the three hydrodynamic mod-els, we note qualitatively similar features as already observed in Fig. 15 fortheir transverse evolution. While the transverse velocities shown in Fig. 15bvaried between the models as a result of differences in the transverse pressures(mainly the bulk viscous pressures), here differences in their longitudinalpressures result in very different longitudinal flow profiles (Figure 16b). Mostnotably, the large negative P L in standard viscous hydrodynamics ( VH 2) ini-tially causes τ u η to reverse sign very sharply in the forward and backwardrapidity regions. Only by imposing strong regulations on the VH VH ) although there the situa-tion is not nearly as bad. With VAH we can maintain positive longitudinalpressures so that the fireball can expand without imploding in regions wherethe longitudinal gradients are large. This greatly reduces the amount ofregulation needed for the viscous pressures.Finally, in Fig. 17 we plot the adaptive time step for each of the three hy-drodynamic simulations. One sees that in anisotropic hydrodynamics ( VAH )∆ τ n increases steadily until hitting the CFL bound at τ ∼ c . Its be-havior is similar in quasiparticle viscous hydrodynamics ( VH ). The two CFLbounds are slightly separated due to differences between their transverse ve-locity profiles (see Fig. 15b). On the other hand, the adaptive time step instandard viscous hydrodynamics ( VH 2) does not pick up until τ ∼ . c .010.11.010100 ( a ) ℰ sw = / fm ℰ τ = / c T = = - - ( b ) τ u η VAHVHVH2 Δ x = Δ y = Δη s = τ = / c τ = / c τ = / c τ = / c τ = / c τ = / c ( c ) P L P ⟂ ( d ) - W ⟂ z · W ⟂ z P L + P ⟂ x ⟂ = ( ) fm - - - - - - ( e ) P L - P ⟂ P eq - - - - - - ( f ) Π P eq η s Figure 16: (Color online) The evolution of (a) E (GeV/fm ), (b) τ u η , (c) P L / P ⊥ ,(d) √ W ⊥ z · W ⊥ z / (cid:112) P L +2 P ⊥ , (e) ( P L −P ⊥ ) / P eq , and (f) Π / P eq along the η s –axis( x = y = 0), given by non-conformal anisotropic hydrodynamics ( VAH , dashed color),quasiparticle viscous hydrodynamics ( VH , transparent color) and standard viscous hy-drodynamics ( VH 2, dotted transparent color), for the smooth (3+1)–d T R ENTo initialcondition. Note that in panel (d) the η s –axis is shifted transversely to ( x, y ) = (7 . , 0) fm. .1 1.0 1010 - τ ( fm / c ) Δ τ ( f m / c ) Δτ n Δτ CFL δ = α = Δη s = Δ x = Δ y = ( a ) Figure 17: (Color online) The evolution of the adaptive time step ∆ τ n (dashed color)and CFL bound ∆ τ CFL (transparent color) in non-conformal anisotropic hydrodynam-ics ( VAH , red), quasiparticle viscous hydrodynamics ( VH , blue) and standard viscoushydrodynamics ( VH 2, green), for the smooth (3+1)–d T R ENTo initial condition. (see footnote 29). Even then, it fluctuates up and down before reaching itsCFL bound, which is higher than the other two bounds since its transverseflow is suppressed by a large bulk viscous pressure. Although the code runsfaster per time step (see Table 1), standard viscous hydrodynamics takesconsiderably more time steps to evolve the fluid than the other two hydro-dynamic models, suggesting that it has a harder time resolving the fireballevolution in the presence of large gradients. T R ENTo initial conditions Finally we study the differences among the three non-conformal hydrody-namic models for a central Pb+Pb collision with a fluctuating initial energydensity profile, using the same model parameters as in Sec. 4.5.1. In the T R ENTo model used in this work, the energy density profile only has fluc-tuations in the transverse directions; the longitudinal profile is modeled witha finite rapidity plateau with smooth slopes at its longitudinal ends (seeAppendix D).Figure 18 shows the evolution of the fluctuating energy density profile inthe transverse plane ( η s = 0), along with the transverse slice of a particliza-59 igure 18: (Color online) The evolution of the QCD energy density profile (GeV/fm ) inthe central transverse plane ( η s = 0), given by non-conformal anisotropic hydrodynamics( VAH , left column), quasiparticle viscous hydrodynamics ( VH , middle column) and stan-dard viscous hydrodynamics ( VH 2, right column) for the fluctuating (3+1)–d T R ENTo event described in the text. The white contour lines are η s = 0 slices at the shown timeframes of a particlization hypersurface of constant energy density E sw = 0 . 116 GeV/fm . E sw = 0 . 116 GeV/fm , for thethree hydrodynamic simulations. Since anisotropic hydrodynamics generallyhas a smaller shear stress than viscous hydrodynamics (especially at earlytimes), the transverse fluctuations across its fireball show the least dissipationor smearing. As a result, it can convert the initial-state eccentricities intoanisotropic flow slightly more efficiently. The pressure anisotropy P L −P ⊥ and bulk viscous pressure Π mainly influence the overall size of the fireballon the particlization hypersurface, affecting the final-state particle yields.The initially strong longitudinal expansion rapidly cools down the systemand temporarily shrinks the fireball size at early times; a more positive lon-gitudinal pressure helps cool the fireball even further. Compared to thetwo viscous hydrodynamic models, anisotropic hydrodynamics has a largerlongitudinal pressure, which means its hypersurface has a slightly narrowerwaist. At later times, the transverse expansion overtakes the longitudinalexpansion, increasing the fireball size. Although anisotropic hydrodynamicshas the largest transverse flow, it also has the smallest bulk viscous pres-sure, enabling the fireball to cool faster and evaporate more quickly. Thisultimately results in a smaller maximum fireball size. In contrast, standardviscous hydrodynamics has a much larger bulk viscous pressure. This ex-tends the fireball’s lifetime by about 1 fm /c relative to the one in anisotropichydrodynamics, allowing it to grow larger in size.In Figs. 19 and 20 we plot the inverse Reynolds numbers from the three hy-drodynamic simulations in the τ − x plane at y = η s = 0 and in the τ − η s planeat x = y = 0, respectively. Traditionally, the inverse Reynolds number mea-sures the validity of the hydrodynamic expansion around local equilibrium[5, 20]. Since anisotropic hydrodynamics expands the energy-momentumtensor (1) around an anisotropic background (i.e. T µνa = E u ν u ν + P L z µ z ν − P ⊥ Ξ µν ), its validity is measured by the residual shear inverse Reynoldsnumber [49] Re − πW = (cid:115) π ⊥ · π ⊥ − W ⊥ z · W ⊥ z P L + 2 P ⊥ , (103a)as opposed to the shear and bulk inverse Reynolds numbers in second-orderviscous hydrodynamics [5, 20]:Re − π = √ π · π P eq √ , Re − = | Π |P eq , (104)61 igure 19: (Color online) The τ − x slice at y = η s = 0 of the residual shear and mean-fieldinverse Reynolds numbers in anisotropic hydrodynamics ( VAH , top row) and the shearand bulk inverse Reynolds numbers in quasiparticle viscous hydrodynamics ( VH , middlerow) and standard viscous hydrodynamics ( VH 2, bottom row), for the fluctuating (3+1)–d T R ENTo event. The white contour lines are τ − x slices ( y = η s = 0) of the sameparticlization hypersurface as Fig. 18. Spacetime regions that are regulated according toSec. 3.5 and App. E are circled in black. igure 20: (Color online) The same as Fig. 19 but for the τ − η s plane at x = y = 0. The η s –axis in the upper left panel is shifted transversely to ( x, y ) = (7 . , 0) fm; the correspondinghypersurface slice (dashed white) is smaller than the one at x = y = 0 (solid white). π · π = π µν π µν . One sees that the residual shear inverse Reynoldsnumbers stay much smaller than unity inside the fireball, indicating thata first-order expansion in O (Re − πW ) is sufficient to capture the residual shearcorrections to the anisotropic hydrodynamic equations. While there is noneed to regulate their strength here, the regulation scheme might be neededfor collision events with sharper initial-state fluctuations.We also plot the contribution of the nonequilibrium mean-field component δB to the bulk viscous pressure in anisotropic hydrodynamics:Re − δB = | δB |P eq . (105)We find that δB is moderately large and positive near the edges of the fireballbut small and negative inside the fireball. However, it is necessary to regulatethe mean-field in the latter region via Eq. (80) to prevent unstable growth.Specifically, the instability seems to originate near the quark-hadron phasetransition, where the driving force proportional to the kinetic trace anomaly E ( k ) − P ( k ) ⊥ − P ( k ) L = E − P ⊥ − P L − B is at its strongest while the bulkrelaxation rate τ − is near a minimum. After applying the regulation (80)for a brief period, δB evolves freely from that point onward. While thisinitial instability is unfortunate, the overall impact of the non-equilibriummean-field component on the longitudinal and transverse pressures insidethe fireball region is limited.In second-order viscous hydrodynamics, the shear inverse Reynolds num-ber is large for a short period of time ∆ τ ∼ c after the collision. Thebulk inverse Reynolds number in standard viscous hydrodynamics is also verylarge near the peak of the bulk viscosity ζ/ S during the same time frame.To prevent the code from crashing, our regulation scheme suppresses boththe shear stress and bulk viscous pressure (see Appendix E). Although it isstill technically possible to run the heavy-ion simulation with standard vis-cous hydrodynamics, we cannot escape the effects of regulation around thebase of the particlization hypersurface. In contrast, the viscous pressures inquasiparticle viscous hydrodynamics often do not require regulation sincethe bulk inverse Reynolds number is quite small. Even in the absence of regulation, the longitudinal pressure in quasiparticle viscoushydrodynamics can still be negative in some regions at early times ( τ < . c ) (e.g.see Figs. 15 and 16). ean time ( s ) Max time ( s ) Mean time per step ( s ) Mean T R ENTo events) on single-core IntelXeon E5-2680 v4 CPUs for each of the three hydrodynamic models. We close this section with a comparison of the particlization hypersur-faces in the τ − x and τ − η s planes. In Fig. 19 one sees that standard viscoushydrodynamics has the largest hypersurface along the x –direction since itslarge bulk viscous pressure generates the most viscous heating. It also hasthe longest lifetime τ f ∼ 17 fm/ c , for the same reason. In contrast, it has apretty narrow waist along the η s –direction because its longitudinal flow pro-file u η initially contracts the medium in the rapidity direction. In anisotropichydrodynamics, the longitudinal flow begins transporting the medium awayfrom the collision zone, following the direction of the longitudinal pressuregradients. As a result, its hypersurface has a wider waist than the one instandard viscous hydrodynamics (about ∆ η s ∼ . τ ∼ − c )along the η s –direction. 5. Benchmarks In this Section we benchmark the typical computational time neededto run (2+1)–d non-conformal hydrodynamic simulations of Pb+Pb colli-sions at LHC energies ( √ s NN = 2 . 76 TeV) with fluctuating initial conditionsand different values for the impact parameter b and model parameters. Wealso benchmark the OpenMP–accelerated runtime of (3+1)–d non-conformalhydrodynamics for a central Pb+Pb collision ( b = 0 fm) with the smooth T R ENTo initial condition and model parameters from Sec. 4.5.1. In this test, we generate 200 random samples for the impact parameter b and Bayesian model parameters [1, 2]. The impact parameter is sampled65 500 1000 1500 2000 Time (s) E v e n t s Fixed gridAuto grid Figure 21: (Color online) The runtime distribution of the (2+1)–d non-conformalanisotropic hydrodynamic simulations on the fixed grid (blue) and auto-grid (red). from a piecewise linear distribution, P ( b ) = (cid:40) b/ (2 R A ) (0 ≤ b ≤ R A ) , b > R A ) , (106)where we set R A = 7 fm for the lead nuclear radius. The model parameters P B = (cid:2) N, p, w, d min , σ k , T sw , ( η/ S ) kink , T η , a low , a high , ( ζ/ S ) max , T ζ , w ζ , λ ζ (cid:3) (107)are each sampled from a distribution that is uniform within the finite intervalsused in the JETSCAPE SIMS Bayesian analysis of heavy-ion collisions [1,2]. , For each parameter sample { b, P B } we run 50 (2+1)–d non-conformalhydrodynamic simulations with fluctuating T R ENTo initial conditions on The T R ENTo initial condition model parameters ( N , p , w , d min , σ k ), are defined inAppendix D [24] and the viscosity model parameters (( η/ S ) kink , T η , a low , a high , ( ζ/ S ) max , T ζ , w ζ , λ ζ ) are discussed in Sec. 2.6.1. Finally, the switching temperature T sw determinesthe particlization hypersurface of constant energy density E sw = E ( T sw ). We exclude the ratio between the dimensionless shear relaxation time and shear vis-cosity b π = τ π T / ( η/ S ) as a continuous model parameter [1, 2] and instead allow it to varybetween discrete hydrodynamic models. The rapidity plateau model parameters η flat and σ η (see Appendix D) are also not considered in this benchmark test. 66 30 fm × 30 fm transverse grid ( η s = 0). , The resulting runtime statis-tics of the three hydrodynamic models are shown in Table 1. On average, ittakes about 85 s to run each (2+1)–d viscous hydrodynamic simulation on thefixed grid, compared to 200 s for anisotropic hydrodynamics. The additionalroutine in Sec. 3.4 makes the runtime per step in anisotropic hydrodynam-ics about 3 . × slower than in viscous hydrodynamics, but this is partiallycompensated by the fewer time steps required to finish each simulation.Figure 21 shows the runtime distribution of the anisotropic hydrodynamicsimulations on the fixed grid. One sees that some runs take longer than oth-ers, depending on the parameter values used. For example, a smaller impactparameter b increases the participant nucleon multiplicity, which extends thefireball’s lifetime. The runtime is also very sensitive to the nucleon width pa-rameter w : halving the value for w doubles the spatial resolution required to capture the spatial variations of the fluctuating T R ENTo profile (see Ap-pendix D) and therefore increases the runtime by roughly a factor of eight.As a result, the maximum runtime can be about an order of magnitude longerthan the mean runtime (see Table 1). The previous benchmark test was performed on a fixed transverse grid of30 fm × 30 fm. Although this grid is large enough to fit all possible fireballsizes produced in Pb+Pb collisions, evolving peripheral collisions ( b (cid:46) R A )on it is computationally inefficient. A large grid is also unnecessary forcertain Bayesian model parameter combinations that further decrease thefireball size. In this section, we present a regression algorithm that predictsthe fireball radius for a given set of parameters { b, P B } . With this, we canautomatically optimize the grid size to save computational time and memory.To train the regression model, we generate 1000 training parameter sam- All other runtime parameters (e.g. the initial time τ and pressure ratio R ) are thesame as the ones used in Sec. 4.5. We found that 93% of the anisotropic hydrodynamic simulations finish successfully(the second-order viscous hydrodynamic simulations had a success rate of 99+%). Forcertain model parameter sets and/or fluctuating initial conditions, the code fails to recon-struct the anisotropic variables in the cold dilute regions, either because the maximumnumber of Newton iterations was exhausted or the kinetic longitudinal pressure P ( k ) L turnednegative. In this work, we set the transverse lattice spacing to ∆ x = ∆ y = w . b [ f m ] -0.600.6 p w [ f m ] d m i n [ f m ] T s w [ G e V ] b [fm] r [ f m ] -0.6 0 0.6 p w [fm] d min [fm] T sw [GeV] r [fm] Figure 22: (Color online) A subset of the scattering matrix used to train the automatedgrid for (2+1)–d non-conformal anisotropic hydrodynamic simulations. ples arranged in the column vectors A = [ b , P B ] , (108)to construct our feature matrix [91]. For each parameter sample, we runa single hydrodynamic simulation with a smooth, event-averaged T R ENTo initial condition and compute the mean fireball radius ¯ r (again arranged in68 AH VH VH2 δ ¯ r RMSE . 29 fm 0 . 19 fm 0 . 22 fm α . 013 0 . 009 0 . Table 2: The cross-validated root-mean-square error δ ¯ r RMSE and regularization parameter α of the Lasso regression fit for each hydrodynamic model. vector); this will be our target variable y = ¯ r . (109)A sample subset of the scattering matrix containing the five most importantfeatures is shown in Fig. 22. As expected, the fireball size strongly decreaseswith the impact parameter. There are also some moderately positive cor-relations between the fireball radius and initial condition parameters (e.g.increasing the nucleon width parameter w increases the spread of the partic-ipant nucleons’ thickness function, see Appendix D, Eq. (139)). While theviscosity parameters have a much smaller effect on the fireball radius individ-ually, the model’s fit to the train-validation data improves when we includeall of them.Next, we fit a cubic polynomial Lasso regression model with standardiza-tion feature scaling [91] to the training data (108)–(109). The regularizationparameter α is chosen to minimize the root-mean-square error δ ¯ r RMSE , aver-aged over a five-fold cross-validation [91]. The values for the cross-validatederror and regularization parameter of the regression fit are listed in Table 2.Finally, we rerun the simulations with fluctuating T R ENTo initial condi-tions, using the 200 test parameter samples from earlier. For each parametersample, we launch the regression model to predict the mean fireball radius¯ r pred and set the transverse grid lengths to L x = L y = 2 (cid:0) ¯ r pred + δ ¯ r RMSE + (cid:96) (cid:1) , (110)where the margin parameter (cid:96) gives the fireball additional room to expandwithin the grid. Here we set (cid:96) = 2 . We define the fireball radius r as the maximum transverse radius of the particlizationhypersurface in a given hydrodynamic simulation. R ENTo events, which were not considered in the training routine. Wefind that the automated grid has about a 99 . 6% success rate of enclosingthe fireball without touching its edges. Fig. 21 shows the anisotropic hy-drodynamic runtime distribution on the automated grid. On average, theautomated grid algorithm reduces the grid area by 32 − 34 % relative to thefixed grid and provides the simulation with an additional 1 . − . × speedup(see Table 1).While (3+1)–d hydrodynamic simulations stand to benefit the most froman optimized grid volume, we would need to implement more realistic longi-tudinal initial conditions than the model used in this work before retrainingthe regression model to predict (in addition to its transverse radius) thefireball’s average elongation along the η s –direction. There may also be ad-ditional parameters that characterize such fluctuating (3+1)–d initial-statemodels, depending on the collision system of interest. Since it takes muchmore computational resources to produce such a training data set, we leavethe (3+1)–d automated grid for future work. The code also includes the option to use OpenMP acceleration on a multi-core processor. This feature is especially useful for speeding up (3+1)–dhydrodynamic simulations, which run typically about two orders of magni-tude longer than (2+1)–d simulations [20]. Fig. 23 shows the runtime ofthe (3+1)–d non-conformal hydrodynamic simulations from Sec. 4.5.1. Theanisotropic hydrodynamic simulation takes about two and a half hours tofinish on a single-core CPU; the runtime of central Pb+Pb collisions can beabout an order of magnitude longer or shorter than this depending on the val-ues used for the nucleon-width parameter w and rapidity plateau parameters η flat and σ η (see Appendix D). On an Intel Xeon E5-2680 multi-core proces-sor, we can significantly reduce the simulation time to about nine minuteswith 28 CPU threads. However, we do not achieve a perfect speedup because With enough computing resources, the user has the option to run multiple fluctuatinghydrodynamic simulations per parameter sample to produce a statistical distribution forthe fireball radius r , which can be characterized with a mean radius and standard deviation Y = [ ¯ r , σ r ]. For more details, we refer the reader to https://github.com/mjmcnelis/fireball . For even faster runtimes, one can parallelize relativistic hydrodynamics on a graphicsprocessing unit [20, 76]. The present version of VAH does not yet offer this option. ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● Threads R un t i m e ( s ) VAHVHVH2 Figure 23: (Color online) The simulation runtime of (3+1)–d non-conformal anisotropichydrodynamics (red dots), quasiparticle viscous hydrodynamics (blue dots) and standardviscous hydrodynamics (green dots) on an Intel Xeon E5-2680 v4 multi-core processor forthe smooth T R ENTo initial condition from Sec. 4.5.1. The ideal runtime is inverselyproportional to the number of parallel CPU threads (same-colored continuous lines). some parts of the routine are not easily parallelizable (e.g. the constructionof the particlization hypersurface). 6. Summary and outlook In this work we developed a (3+1)–dimensional anisotropic fluid dynam-ical simulation that evolves both the pre-equilibrium and viscous hydrody-namic stages of ultrarelativistic nuclear collisions. We validated the code’sperformance by reproducing the semi-analytic solutions of conformal andnon-conformal Bjorken flow and conformal Gubser flow on a Cartesian grid.Thanks to the adaptive time step algorithm derived in Sec. 3.6, we can ac-curately capture the early-time dynamics of Bjorken and Gubser flow whilefinishing the simulation within a reasonable number of iterations. We alsocompared anisotropic hydrodynamics to two second-order viscous hydrody-namic models in central Pb+Pb collisions. Apart from the apparent sen-sitivity of the bulk viscous pressure evolution to the bulk relaxation time,the three hydrodynamic models have similar transverse dynamics, at least71n the mid-rapidity region η s = 0. However, the fluid’s longitudinal evolu-tion varies greatly near the longitudinal edges of the fireball. We showedthat the longitudinal pressure in (3+1)–dimensional anisotropic hydrody-namics stays positive even in the presence of large gradients at very earlytimes, unlike in viscous hydrodynamics. This causes the fluid to initiallyexpand outward along the spacetime rapidity direction, as expected fromthe outward-pointing longitudinal gradients of the thermal pressure, reduc-ing the risk of cavitation at the beginning of the simulation. With the newdevelopment presented here, we can for the first time model even the veryearly pre-equilibrium dynamics stage at τ (cid:28) τ hydro with a QCD equation ofstate, as opposed to the conformal equation of state implicit in other pre-equilibrium models [44, 45, 50, 92–96].In the near term, we plan to run VAH in the JETSCAPE framework withthe Maximum a Posteriori (MAP) model parameters extracted in Refs. [1,2]. Since the JETSCAPE SIMS hybrid model [1, 2] uses conformal free-streaming and standard viscous hydrodynamics to evolve the pre-equilibriumand hydrodynamic stages, we will replace these two modules with our codeand study the changes to the hadronic observables, such as the transversemomentum spectra and p T –differential anisotropic flows. A full analysis ofall available data will need to wait for a Bayesian recalibration of the fullevolution model using VAH as its hydrodynamic core.Integration of VAH with the JETSCAPE framework also requires an up-dated version of the particlization module iS3D [74] with a new option thatallows choosing the leading-order anisotropic distribution f a (plus residualshear corrections δ ˜ f ) for the hadronic distribution in the Cooper–Frye for-mula [74, 97]. Once completed, this will allow us to investigate how differentselections among a discrete set of hydrodynamic models affect both the the-oretical description of heavy-ion experimental observables and the shear andbulk viscosity constraints inferred from such data–theory comparisons. Inparticular, it would be interesting to test the predictions from anisotropichydrodynamics against a variety of initializations of second-order viscoushydrodynamics using different pre-hydrodynamic evolution models (or no Several model parameters associated with the conformal free-streaming module wouldnot be used in our code, which replaces the free-streaming stage by anisotropic hydrody-namics. Anisotropic hydrodynamics ( VAH ) as presented here is not meant to be integratedwith a pre-hydrodynamic module, especially one that uses a conformal approximation. n B and V µB have alreadybeen practically implemented [83, 98], and others that initialize viscous hy-drodynamics with dynamical sources [10, 54, 55, 57, 58] are currently underdevelopment. The BEShydro code [57, 58], in particular, shares a commonancestry [49] and therefore a number of similar features with VAH . The inte-gration of both codes into a single framework, in order to utilize the adaptivetime step for capturing and resolving the dynamical production of energy andbaryon sources at the onset of low-energy nuclear collisions, offers itself as aninteresting project. On the theoretical side, the effects of non-zero chemicalpotentials for net charge, baryon number and strangeness have not yet beenconsidered in our formulation of anisotropic hydrodynamics [66]. An upgradeof our code package that will include them is planned for the future. 7. Acknowledgements M.M. would like to thank Seyed Sabok-Sayr from Rutgers University forassisting in the development of the automated grid algorithm during theErd˝os Institute 2020 Data Science Boot Camp. Computational resources forthe code validation, comparison and benchmark tests were provided by theOhio Supercomputer Center under Project PAS0254 [99]. This work wassupported by the National Science Foundation (NSF) within the frameworkof the JETSCAPE Collaboration under Award No. ACI-1550223. Additionalpartial support by the U.S. Department of Energy (DOE), Office of Science,Office for Nuclear Physics under Award No. DE-SC0004286 and within theframework of the BEST and JET Collaborations is also acknowledged. Appendix A Gradient source terms In this appendix, we list the gradients that appear in the source terms ofthe anisotropic hydrodynamic equations. The derivatives of the fluid velocity This is mainly due to the technical difficulties in initializing the mean-field and anisotropicvariables [66]. Available second-order viscous hydrodynamic codes can read in the energy-momentum tensor from a pre-hydrodynamic module and thus initialize the dynamical andinferred variables more easily, but this has not yet been implemented in the VAH code. u τ are ∂ τ u τ = v x ∂ τ u x + v y ∂ τ u y + τ v η ∂ τ u η + τ v η u η , (111a) ∂ x u τ = v x ∂ x u x + v y ∂ x u y + τ v η ∂ x u η , (111b) ∂ y u τ = v x ∂ y u x + v y ∂ y u y + τ v η ∂ y u η , (111c) ∂ η u τ = v x ∂ η u x + v y ∂ η u y + τ v η ∂ η u η . (111d)The derivatives of the longitudinal basis vector z µ are ∂ τ z τ = τ (cid:112) u ⊥ (cid:32) ∂ τ u η − u η ( u x ∂ τ u x + u y ∂ τ u y )1+ u ⊥ (cid:33) + z τ τ , (112a) ∂ x z τ = τ (cid:112) u ⊥ (cid:32) ∂ x u η − u η ( u x ∂ x u x + u y ∂ x u y )1+ u ⊥ (cid:33) , (112b) ∂ y z τ = τ (cid:112) u ⊥ (cid:32) ∂ y u η − u η (cid:0) u x ∂ y u x + u y ∂ y u y (cid:1) u ⊥ (cid:33) , (112c) ∂ η z τ = τ (cid:112) u ⊥ (cid:32) ∂ η u η − u η (cid:0) u x ∂ η u x + u y ∂ η u y (cid:1) u ⊥ (cid:33) , (112d) ∂ τ z η = 1 τ (cid:112) u ⊥ (cid:32) ∂ τ u τ − u τ ( u x ∂ τ u x + u y ∂ τ u y )1+ u ⊥ (cid:33) − z η τ , (112e) ∂ x z η = 1 τ (cid:112) u ⊥ (cid:32) ∂ x u τ − u τ ( u x ∂ x u x + u y ∂ x u y )1+ u ⊥ (cid:33) , (112f) ∂ y z η = 1 τ (cid:112) u ⊥ (cid:32) ∂ y u τ − u τ (cid:0) u x ∂ y u x + u y ∂ y u y (cid:1) u ⊥ (cid:33) , (112g) ∂ η z η = 1 τ (cid:112) u ⊥ (cid:32) ∂ η u τ − u τ (cid:0) u x ∂ η u x + u y ∂ η u y (cid:1) u ⊥ (cid:33) . (112h)The divergence of the spatial fluid velocity v i = u i /u τ is ∂ i v i = ∂ x u x − v x ∂ x u τ + ∂ y u y − v y ∂ y u τ + ∂ η u η − v η ∂ η u τ u τ . (113)The longitudinal expansion rate is θ L = z µ D z u µ = − z µ z ν ∂ ν u µ − z µ z ν Γ µνλ u λ = − ( z τ ) ∂ τ u τ + z τ z η ( τ ∂ τ u η − ∂ η u τ ) + ( τ z η ) ∂ η u η + τ ( z η ) u τ . (114)74he transverse expansion rate is θ ⊥ = ∇ ⊥ µ u µ = θ − θ L , (115)where θ = D µ u µ = ∂ µ u µ + Γ µµν u ν = ∂ τ u τ + ∂ x u x + ∂ y u y + ∂ η u η + u τ τ (116)is the scalar expansion rate.The components of the fluid acceleration a µ = Du µ = u ν ∂ ν u µ + u ν Γ µνλ u λ (117)are a τ = u τ ∂ τ u τ + u x ∂ x u τ + u y ∂ y u τ + u η ∂ η u τ + τ ( u η ) , (118a) a x = u τ ∂ τ u x + u x ∂ x u x + u y ∂ y u x + u η ∂ η u x , (118b) a y = u τ ∂ τ u y + u x ∂ x u y + u y ∂ y u y + u η ∂ η u y , (118c) a η = u τ ∂ τ u η + u x ∂ x u η + u y ∂ y u η + u η ∂ η u η + 2 u τ u η τ . (118d)The components of the longitudinal vector’s comoving time derivative˙ z µ = Dz µ = u ν ∂ ν z µ + u ν Γ µνλ z λ (119)are ˙ z τ = u τ ∂ τ z τ + u x ∂ x z τ + u y ∂ y z τ + u η ∂ η z τ + τ u η z η , (120a)˙ z η = u τ ∂ τ z η + u x ∂ x z η + u y ∂ y z η + u η ∂ η z η + u τ z η + u η z τ τ . (120b)The components of the fluid velocity’s longitudinal derivative D z u µ = − z ν ∂ ν u µ − z ν Γ µνλ u λ (121)are D z u τ = − z τ ∂ τ u τ − z η ∂ η u τ − τ u η z η , (122a) D z u x = − z τ ∂ τ u x − z η ∂ η u x , (122b)75 z u y = − z τ ∂ τ u y − z η ∂ η u y , (122c) D z u η = − z τ ∂ τ u η − z η ∂ η u η − u τ z η + u η z τ τ . (122d)The transverse gradient of the fluid velocity projected along the longitudinaldirection is z ν ∇ µ ⊥ u ν = Ξ µα z ν D α u ν , (123)where the components of z ν D α u ν = g αβ z ν ∂ β u ν + g αβ z ν Γ νβλ u λ (124)are z ν D τ u ν = z τ ∂ τ u τ − τ z η ∂ τ u η − τ u η z η , (125a) z ν D x u ν = − z τ ∂ x u τ + τ z η ∂ x u η , (125b) z ν D y u ν = − z τ ∂ y u τ + τ z η ∂ y u η , (125c) z ν D η u ν = − τ (cid:0) z τ ∂ η u τ − τ z η ∂ η u η + τ ( u η z τ − u τ z η ) (cid:1) (125d)The transverse velocity-shear tensor is σ µν ⊥ = Ξ µναβ D ( α u β ) , (126)where the components of D ( α u β ) = 12 (cid:16) g αρ ∂ ρ u β + g βρ ∂ ρ u α + g αρ Γ βρλ u λ + g βρ Γ αρλ u λ (cid:17) (127)are D ( τ u τ ) = ∂ τ u τ , (128a) D ( τ u x ) = 12 ( ∂ τ u x − ∂ x u τ ) , (128b) D ( τ u y ) = 12 (cid:0) ∂ τ u y − ∂ y u τ (cid:1) , (128c) D ( τ u η ) = 12 (cid:18) ∂ τ u η − ∂ η u τ τ (cid:19) , (128d) In the code we obtain σ µν ⊥ by applying the projector Ξ µναβ onto D ( α u β ) directly, ratherthan simplifying the expression (126). ( x u x ) = − ∂ x u x , (128e) D ( x u y ) = − (cid:0) ∂ x u y + ∂ y u x (cid:1) , (128f) D ( x u η ) = − (cid:18) ∂ x u η + ∂ η u x τ (cid:19) , (128g) D ( y u y ) = − ∂ y u y , (128h) D ( y u η ) = − (cid:18) ∂ y u η + ∂ η u y τ (cid:19) , (128i) D ( η u η ) = − τ (cid:18) ∂ η u η + u τ τ (cid:19) . (128j)The transverse vorticity tensor is ω µν ⊥ = Ξ µα Ξ νβ D [ α u β ] = D [ µ u ν ] − u µ a ν − u ν a µ + z µ ( D z u ν + z α D ν u α ) − z ν ( D z u µ + z α D µ u α )2 (129)where the components of D [ α u β ] = 12 (cid:16) g αρ ∂ ρ u β − g βρ ∂ ρ u α + g αρ Γ βρλ u λ − g βρ Γ αρλ u λ (cid:17) (130)are D [ τ u τ ] = 0 , (131a) D [ τ u x ] = 12 ( ∂ τ u x + ∂ x u τ ) , (131b) D [ τ u y ] = 12 (cid:0) ∂ τ u y + ∂ y u τ (cid:1) , (131c) D [ τ u η ] = 12 (cid:18) ∂ τ u η + ∂ η u τ τ (cid:19) + u η τ , (131d) D [ x u x ] = 0 , (131e) D [ x u y ] = − (cid:0) ∂ x u y − ∂ y u x (cid:1) , (131f) D [ x u η ] = − (cid:18) ∂ x u η − ∂ η u x τ (cid:19) , (131g) D [ y u y ] = 0 , (131h) D [ y u η ] = − (cid:18) ∂ y u η + ∂ η u y τ (cid:19) , (131i)77 [ η u η ] = 0 . (131j) Appendix B Geometric source terms Here we list the components of the geometric source terms G µW and G µνπ (Eq. (23)) that appear in the relaxation equations (20) for W µ ⊥ z and π µν ⊥ ,respectively: G τW = τ u η W η ⊥ z , (132a) G xW = 0 , (132b) G yW = 0 , (132c) G ηW = u τ W η ⊥ z + u η W τ ⊥ z τ . (132d) G ττπ = 2 τ u η π τη ⊥ , (133a) G τxπ = τ u η π xη ⊥ , (133b) G τyπ = τ u η π yη ⊥ , (133c) G τηπ = τ u η π ηη ⊥ + u τ π τη ⊥ + u η π ττ ⊥ τ , (133d) G xxπ = 0 , (133e) G xyπ = 0 , (133f) G xηπ = u τ π xη ⊥ + u η π τx ⊥ τ , (133g) G yyπ = 0 , (133h) G yηπ = u τ π yη ⊥ + u η π τy ⊥ τ , (133i) G ηηπ = 2 (cid:0) u τ π ηη ⊥ + u η π τη ⊥ (cid:1) τ . (133j) Appendix C Conformal anisotropic transport coefficients In the conformal limit m = B = 0, the anisotropic transport coeffi-cients (35) – (38) only depend the functions I nrqs , which reduce to I nrqs = g ( n + s +1)! α r +1 L Λ n + s +2 R nrq π (2 q )!! , (134)78here g , α L and Λ are given by Eqs. (27), (72) and (73), respectively. Welist the functions R nrq used in this work [66]: R = α L (cid:0) ξ L ) t L (cid:1) , (135a) R = − ξ L ) t L ξ L α L , (135b) R = 1 + ( ξ L − t L ξ L α L , (135c) R = 3 + 2 ξ L − ξ L ) t L ) ξ L α L , (135d) R = 3 + ξ L + (1 + ξ L )( ξ L − t L ξ L (1 + ξ L ) α L , (135e) R = − ξ L ) t L ) ξ L α L , (135f) R = − 15 + 13 ξ L + 3(1 + ξ L )(5 + ξ L ) t L ξ L α L , (135g) R = 3( ξ L − 1) + ( ξ L (3 ξ L − 2) + 3) t L ξ L α L , (135h) R = 3 + ξ L + (1 + ξ L )( ξ L − t L ξ L α L , (135i) R = 15 + ξ L + ( ξ L ( ξ L − − t L ξ L α L , (135j) R = ( ξ L − ξ L ) + 3(1 + ξ L )( ξ L ( ξ L − 2) + 5) t L ξ L (1 + ξ L ) α L , (135k)where ξ L = α − L − t L = arctan √ ξ L / √ ξ L . Appendix D T R ENTo energy deposition model In the T R ENTo model, the transverse energy deposition (GeV/fm ) ofa single fluctuating nuclear collision event in the mid-rapidity region is [24] dE T dxdydη s | η s =0 = N × T R ( x, y ) , (136)where N is the normalization parameter and T R ( x, y ) = (cid:18) T pA ( x, y ) + T pB ( x, y )2 (cid:19) /p (137)79s the reduced nuclear thickness function, with p being the geometric param-eter. The nuclear thickness function of nucleus A ( B ) is T A,B ( x, y ) = N part ,A,B (cid:88) n =1 γ n T p ( x − x n , y − y n ) , (138)where N part ,A,B are the number of participant nucleons from nucleus A ( B );the nucleon positions are sampled from a Woods–Saxon distribution underthe constraint that each nucleon–nucleon pair in nucleus A ( B ) maintains aminimum separation d min [24]. The participant nucleon’s thickness functionis centered around its sampled transverse position ( x n , y n ): T p ( x − x n , y − y n ) = 12 πw × exp (cid:34) − ( x − x n ) + ( y − y n ) w (cid:35) , (139)where w is the nucleon width parameter. Furthermore, the multiplicityfactor γ n of each participant nucleon is sampled from the gamma distribution P ( γ ) = k k γ k − exp [ − kγ ]Γ[ k ] , (140)where k = σ − k and σ k is the standard deviation.For a very brief period τ after the collision, we assume the system islongitudinally free-streaming and static in Milne coordinates (i.e. P L / P eq ∼ u ∼ ) so that E ( x ) ∝ /τ . Therefore, we initialize the energy densityprofile of the hydrodynamic simulation at the starting time τ as E ( τ , x, y, η s ) = 1 τ × dE T dxdydη s | η s =0 × f L ( η s ) . (141)Here we also extend the transverse energy density profile along the spacetimerapidity direction with a smooth plateau distribution (unitless) [76]: f L ( η s ) = exp − (cid:0) | η s | − η flat (cid:1) Θ (cid:0) | η s | − η flat (cid:1) σ η , (142) In this work, we only consider Pb+Pb collisions ( A = B = 208) at LHC energies √ s NN = 2 . 76 TeV; the inelastic nucleon–nucleon cross section is set to σ NN = 6 . − . This parameter does not control the nucleon’s charge radius but rather the spread ofthermal energy it deposits in the collision zone along the transverse directions. η flat is the plateau length, σ η is the standard deviation of the half-Gaussian tails and Θ is the Heaviside step function. The initial conditionparameter values used in this work are N = 14 . 19 GeV, p = 0 . w = 1 . d min = 1 . 45 fm, σ k = 1 . η flat = 4 . σ η = 1 . x = ∆ y = w and ∆ η s = σ η to resolvethe fluctuating energy density profile (141) (or event-averaged profile). Forthe grid size, we set the longitudinal length to L η = η flat + 10 σ η to fit therapidity plateau (142). The transverse lengths L x and L y are automaticallyconfigured by the algorithm described in Sec. 5.2. Appendix E Numerical implementation of second-order viscoushydrodynamics In this appendix, we summarize how second-order viscous hydrodynam-ics is implemented the code. We evolve viscous hydrodynamics with thesame numerical algorithm discussed in Sec. 3 except the energy-momentumtensor (1) is decomposed as T µν = E u µ u ν − ( P eq +Π)∆ µν + π µν , (143)where π µν = ∆ µναβ T αβ is the shear stress tensor and Π = − ∆ µν T µν −P eq is thebulk viscous pressure. We also define the spatial projector ∆ µν = g µν − u µ u ν and traceless double spatial projector ∆ µναβ = (∆ µα ∆ νβ + ∆ νβ ∆ µα − ∆ µν ∆ αβ ).The corresponding dynamical variables q = ( T τµ , π µν , Π) (144)are propagated along with E and u (the mean-field B and anisotropic vari-ables (Λ, α ⊥ , α L ) are not propagated). Although π µν has only five indepen-dent components, we propagate all ten components in the simulation [20]. E.1 Hydrodynamic equations Here we list the evolution equations for the dynamical variables (144) (werefer the reader to Refs. [18, 20] for details on their derivation): ∂ τ T ττ + ∂ i ( v i T ττ ) = − T ττ + τ T ηη τ + ( π ττ −P eq − Π) ∂ i v i (145a) For longitudinally boost-invariant systems, we do not propagate the shear stress com-ponents π τη , π xη and π yη . v i ∂ i ( π ττ −P eq − Π) − ∂ i π τi ,∂ τ T τx + ∂ i ( v i T τx ) = − T τx τ − ∂ x ( P eq +Π) + π τx ∂ i v i + v i ∂ i π τx (145b) − ∂ i π xi ,∂ τ T τy + ∂ i ( v i T τy ) = − T τy τ − ∂ y ( P eq +Π) + π τy ∂ i v i + v i ∂ i π τy (145c) − ∂ i π yi ,∂ τ T τη + ∂ i ( v i T τη ) = − T τη τ − ∂ η ( P eq +Π) τ + π τη ∂ i v i + v i ∂ i π τη (145d) − ∂ i π ηi ,∂ τ π µν + ∂ i ( v i π µν ) = π µν ∂ i v i + 1 u τ (cid:20) − π µν τ π + I µνπ (cid:48) − P µνπ (cid:48) − G µνπ (cid:48) (cid:21) , (145e) ∂ τ Π + ∂ i ( v i Π) = Π ∂ i v i + 1 u τ (cid:20) − Π τ Π + I Π (cid:21) , (145f)where T ηη = ( E + P eq +Π)( u η ) + ( P eq +Π) /τ + π ηη , I µνπ (cid:48) = 2 β π σ µν + ∆ µναβ (cid:0) π λ ( α ω β ) λ − ¯ τ ππ π λ ( α σ β ) λ (cid:1) − ¯ δ ππ π µν θ (146a)+ ¯ λ π Π Π σ µν , I Π = − β Π θ − ¯ δ ΠΠ Π θ + ¯ λ Π π π µν σ µν , (146b)are the gradient source terms for π µν and Π and P µνπ (cid:48) = ( π µα u ν + π να u µ ) a α , (147a) G µνπ (cid:48) = u γ Γ µγλ π νλ + u γ Γ νγλ π µλ , (147b)are the spatial projection and geometric source terms for π µν . We also de-fine the velocity-shear tensor σ µν = ∆ µναβ D ( α u β ) and vorticity tensor ω µν =∆ µα ∆ νβ D [ α u β ] . 82he components of σ µν and G µνπ (cid:48) are the same as σ µν ⊥ and G µνπ after replac-ing Ξ µναβ → ∆ µναβ and π µν ⊥ → π µν in Eqs. (126) and (133), respectively. Thecomponents of ω µν are ω µν = D [ µ u ν ] − u µ a ν − u ν a µ , (148)where D [ µ u ν ] is given by Eq. (130). E.2 Transport coefficients In quasiparticle viscous hydrodynamics [79] (i.e. m = m ( T ) from Fig. 2a),the relaxation times ( τ π , τ Π ) and first-order transport coefficients ( β π , β Π )are given by Eqs. (31) – (32). The second-order transport coefficients are [66] τ ππ = 10 + 4¯ c π m I , (149a) δ ππ = 43 + ¯ c π I (cid:16) m − m dmd E ( E + P eq ) (cid:17) , (149b) λ π Π = 65 − m (cid:0) ¯ c E I + ¯ c Π I (cid:1) , (149c) δ ΠΠ = 1 − c s − m (cid:0) ¯ c E I + ¯ c Π I (cid:1) (149d) − m dmd E ( E + P eq ) (cid:16) ¯ c E I + 53 ¯ c Π I + 3 m (cid:17) ,λ Π π = 13 − c s + ¯ c π m I , (149e)where ¯ c π = 12 I , (150a)¯ c E = − I I I − I , (150b)¯ c Π = I I I − I , (150c)and the function I nq are given by Eq. (33).83n standard viscous hydrodynamics [20, 67] (i.e. ¯ m = m/T (cid:28) β π ≈ E + P eq O ( ¯ m ) , (151a) β Π ≈ (cid:16) − c s (cid:17) ( E + P eq ) + O ( ¯ m ) . (151b)The second-order transport coefficients (149) reduce to τ ππ ≈ 107 + O ( ¯ m ) , (152a) δ ππ ≈ 43 + O ( ¯ m ) , (152b) λ π Π ≈ 65 + O ( ¯ m ln ¯ m ) , (152c) δ ΠΠ ≈ 23 + O ( ¯ m ln ¯ m ) , (152d) λ Π π ≈ (cid:16) − c s (cid:17) + O ( ¯ m ) . (152e)Both models use the same shear and bulk viscosities as anisotropic hydrody-namics (e.g. Fig. 3). E.3 Reconstruction formulas for the energy density and fluid velocity We reconstruct the energy density by solving the following nonlinear equa-tion via Newton’s method [5, 20]: f ( E ) = 0 , (153)where f ( E ) = (cid:0) ¯ M τ −E (cid:1) (cid:0) ¯ M τ + P eq +Π (cid:1) − ( ¯ M x ) − ( ¯ M y ) − ( τ ¯ M η ) , (154)with ¯ M µ = T τµ − π τµ and P eq = P eq ( E ) being the QCD equation of state.Using the previous energy density for the initial guess, we iterate the solutionto Eq. (153) as E ← E − f ( E ) df /d E , (155)where dfd E = c s ( ¯ M τ −E ) − ¯ M τ − P eq − Π , (156)84ith c s = d P eq /d E being the QCD speed of sound squared. We repeat theiteration (155) until we achieve sufficient convergence or the energy densityfalls below E min . If the bulk viscous pressure Π < −P eq , we regulate it sothat P eq + Π = 0; this allows us to solve for E explicitly [5]: E = ¯ M τ − ( ¯ M x ) +( ¯ M y ) +( τ ¯ M η ) ¯ M τ . (157)Afterwards, we regulate the energy density via Eq. (63) and evaluate thefluid velocity components as u x = ¯ M x (cid:113)(cid:0) E + P eq +Π (cid:1) (cid:0) ¯ M τ + P eq +Π (cid:1) , (158a) u y = ¯ M y (cid:113)(cid:0) E + P eq +Π (cid:1) (cid:0) ¯ M τ + P eq +Π (cid:1) , (158b) u η = ¯ M η (cid:113)(cid:0) E + P eq +Π (cid:1) (cid:0) ¯ M τ + P eq +Π (cid:1) . (158c) E.4 Regulating the shear stress and bulk viscous pressure In this regulation scheme, we first adjust the shear stress components π ηη ← τ (cid:0) u ⊥ (cid:1) (cid:104) π xx (cid:0) ( u x ) − ( u τ ) (cid:1) + π yy (cid:0) ( u y ) − ( u τ ) (cid:1) + 2 (cid:0) π xy u x u y + τ ( π xη u x + π yη u y ) u η (cid:1)(cid:105) , (159a) π τx ← π xx u x + π xy u y + τ π xη u η u τ , (159b) π τy ← π xy u x + π yy u y + τ π yη u η u τ , (159c) π τη ← π xη u x + π yη u y + τ π ηη u η u τ , (159d) π ττ ← π τx u x + π τy u y + τ π τη u η u τ , (159e)so that π µν satisfies the orthogonality and tracelessness conditions π µν u ν = 0 , (160a) π µµ = 0 . (160b)85hen we regulate the shear stress and bulk viscous pressure as π µν ← γ reg π µν , (161a)Π ← γ reg Π , (161b)where γ reg = min (cid:32) , (cid:115) P π · π + 3Π (cid:33) , (162)with π · π = π µν π µν . The regulation factor γ reg usually suppresses π µν andΠ around the edges of the fireball at early times τ < c , especially instandard viscous hydrodynamics (e.g. see Figs. (19c) and (20c)). References [1] D. Everett, et al., Phenomenological constraints on the transport prop-erties of QCD matter with data-driven model averaging, arXiv:2010.03928 .[2] D. Everett, et al., Multi-system Bayesian constraints on the transportcoefficients of QCD matter, arXiv:2011.01430 .[3] P. Kovtun, D. T. Son, A. O. Starinets, Holography and hydrodynamics:Diffusion on stretched horizons, JHEP 10 (2003) 064. arXiv:hep-th/0309213 , doi:10.1088/1126-6708/2003/10/064 .[4] J. M. Maldacena, The Large N limit of superconformal field theoriesand supergravity, Int. J. Theor. Phys. 38 (1999) 1113–1133. arXiv:hep-th/9711200 , doi:10.1023/A:1026654312961 .[5] C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, U. Heinz, The iEBE-VISHNU code package for relativistic heavy-ion collisions, Comput.Phys. Commun. 199 (2016) 61–85. arXiv:1409.8164 , doi:10.1016/j.cpc.2015.08.039 .[6] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, U. Heinz, Ap-plying Bayesian parameter estimation to relativistic heavy-ion colli-sions: simultaneous characterization of the initial state and quark-gluonplasma medium, Phys. Rev. C94 (2016) 024907. arXiv:1605.03954 , doi:10.1103/PhysRevC.94.024907 .867] J. E. Bernhard, Bayesian parameter estimation for relativistic heavy-ion collisions, Ph.D. thesis, Duke University (2018-04-19). arXiv:1804.06469 .[8] G. Nijs, W. van der Schee, U. G¨ursoy, R. Snellings, A transversemomentum differential global analysis of heavy ion collisions, arXiv:2010.15130 .[9] G. Nijs, W. van der Schee, U. G¨ursoy, R. Snellings, A Bayesian analysisof heavy ion collisions with Trajectum, arXiv:2010.15134 .[10] T. Hirano, P. Huovinen, K. Murase, Y. Nara, Integrated dynamicalapproach to relativistic heavy ion collisions, Prog. Part. Nucl. Phys.70 (2013) 108–158. arXiv:1204.5814 , doi:10.1016/j.ppnp.2013.02.002 .[11] J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon,C. Gale, Production of photons in relativistic heavy-ion collisions,Phys. Rev. C93 (4) (2016) 044906. arXiv:1509.06738 , doi:10.1103/PhysRevC.93.044906 .[12] C. Shen, J.-F. Paquet, G. S. Denicol, S. Jeon, C. Gale, Collectivity andelectromagnetic radiation in small systems, Phys. Rev. C95 (1) (2017)014906. arXiv:1609.02590 , doi:10.1103/PhysRevC.95.014906 .[13] M. Okai, K. Kawaguchi, Y. Tachibana, T. Hirano, New approach toinitializing hydrodynamic fields and mini-jet propagation in quark-gluonfluids, Phys. Rev. C 95 (5) (2017) 054914. arXiv:1702.07541 , doi:10.1103/PhysRevC.95.054914 .[14] G. Vujanovic, J.-F. Paquet, C. Shen, G. S. Denicol, S. Jeon, C. Gale,U. Heinz, Exploring the influence of bulk viscosity of QCD on dileptontomography, Phys. Rev. C 101 (2020) 044904. arXiv:1903.05078 , doi:10.1103/PhysRevC.101.044904 .[15] X. Yao, W. Ke, Y. Xu, S. A. Bass, B. M¨uller, Coupled Boltzmann trans-port equations of heavy quarks and quarkonia in quark-gluon plasma, arXiv:2004.06746 . 8716] B. Schenke, S. Jeon, C. Gale, (3+1)D hydrodynamic simulation of rel-ativistic heavy-ion collisions, Phys. Rev. C82 (2010) 014903. arXiv:1004.1408 , doi:10.1103/PhysRevC.82.014903 .[17] B. Schenke, S. Jeon, C. Gale, Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics, Phys. Rev. Lett. 106 (2011)042301. arXiv:1009.3244 , doi:10.1103/PhysRevLett.106.042301 .[18] G. S. Denicol, H. Niemi, E. Molnar, D. H. Rischke, Derivation of tran-sient relativistic fluid dynamics from the Boltzmann equation, Phys.Rev. D85 (2012) 114047, [Erratum: Phys. Rev.D 91 (2015) 039902]. arXiv:1202.4551 , doi:10.1103/PhysRevD.85.114047 .[19] C. Gale, S. Jeon, B. Schenke, Hydrodynamic modeling of heavy-ioncollisions, Int. J. Mod. Phys. A28 (2013) 1340011. arXiv:1301.5893 , doi:10.1142/S0217751X13400113 .[20] D. Bazow, U. Heinz, M. Strickland, Massively parallel simulations ofrelativistic fluid dynamics on graphics processing units with CUDA,Comput. Phys. Commun. 225 (2018) 92–113. arXiv:1608.06577 , doi:10.1016/j.cpc.2017.01.015 .[21] L. Rezzolla, O. Zanotti, Relativistic Hydrodynamics, Oxford UniversityPress, London, 2013. doi:10.1093/acprof:oso/9780198528906.001.0001 .[22] B. Schenke, P. Tribedy, R. Venugopalan, Fluctuating Glasma initialconditions and flow in heavy ion collisions, Phys. Rev. Lett. 108 (2012)252301. arXiv:1202.6646 , doi:10.1103/PhysRevLett.108.252301 .[23] C. Loizides, J. Nagle, P. Steinberg, Improved version of the PHOBOSGlauber Monte Carlo, SoftwareX 1-2 (2015) 13–18. arXiv:1408.2549 , doi:10.1016/j.softx.2015.05.001 .[24] J. S. Moreland, J. E. Bernhard, S. A. Bass, Alternative ansatz towounded nucleon and binary collision scaling in high-energy nuclearcollisions, Phys. Rev. C92 (1) (2015) 011901. arXiv:1412.4708 , doi:10.1103/PhysRevC.92.011901 . 8825] H. Niemi, G. Denicol, How large is the Knudsen number reached in fluiddynamical simulations of ultrarelativistic heavy ion collisions?, arXiv:1404.7327 .[26] C. Gale, S. Jeon, B. Schenke, P. Tribedy, R. Venugopalan, Event-by-event anisotropic flow in heavy-ion collisions from combined Yang-Millsand viscous fluid dynamics, Phys. Rev. Lett. 110 (1) (2013) 012302. arXiv:1209.6330 , doi:10.1103/PhysRevLett.110.012302 .[27] C. Shen, The standard model for relativistic heavy-ion collisions andelectromagnetic tomography, Ph.D. thesis, Ohio State U. (2014-07-25).URL http://rave.ohiolink.edu/etdc/view?acc_num=osu1405931790 [28] R. D. Weller, P. Romatschke, One fluid to rule them all: viscoushydrodynamic description of event-by-event central p+p, p+Pb andPb+Pb collisions at √ s = 5 . 02 TeV, Phys. Lett. B774 (2017) 351–356. arXiv:1701.07145 , doi:10.1016/j.physletb.2017.09.077 .[29] P. Romatschke, Do nuclear collisions create a locally equilibratedquark–gluon plasma?, Eur. Phys. J. C 77 (1) (2017) 21. arXiv:1609.02820 , doi:10.1140/epjc/s10052-016-4567-x .[30] W. Zhao, Y. Zhou, K. Murase, H. Song, Searching for small dropletsof hydrodynamic fluid in proton–proton collisions at the LHC, Eur.Phys. J. C 80 (9) (2020) 846. arXiv:2001.06742 , doi:10.1140/epjc/s10052-020-8376-x .[31] C. Plumberg, Hanbury-Brown–Twiss interferometry and collectivity in p + p , p +Pb, and Pb+Pb collisions, Phys. Rev. C 102 (5) (2020) 054908. arXiv:2008.01709 , doi:10.1103/PhysRevC.102.054908 .[32] C. Plumberg, The multiplicity dependence of pion interferometry in hy-drodynamics, arXiv:2010.11957 .[33] M. P. Heller, M. Spalinski, Hydrodynamics Beyond the Gradient Ex-pansion: Resurgence and Resummation, Phys. Rev. Lett. 115 (7) (2015)072501. arXiv:1503.07514 , doi:10.1103/PhysRevLett.115.072501 .[34] P. Romatschke, U. Romatschke, Relativistic Fluid Dynamics In andOut of Equilibrium, Cambridge Monographs on Mathematical Physics,89ambridge University Press, 2019. arXiv:1712.05815 , doi:10.1017/9781108651998 .[35] M. Strickland, J. Noronha, G. Denicol, Anisotropic nonequilibrium hy-drodynamic attractor, Phys. Rev. D 97 (3) (2018) 036020. arXiv:1709.06644 , doi:10.1103/PhysRevD.97.036020 .[36] W. Florkowski, M. P. Heller, M. Spalinski, New theories of relativistichydrodynamics in the LHC era, Rept. Prog. Phys. 81 (4) (2018) 046001. arXiv:1707.02282 , doi:10.1088/1361-6633/aaa091 .[37] W. van der Schee, P. Romatschke, S. Pratt, Fully dynamical simulationof central nuclear collisions, Phys. Rev. Lett. 111 (22) (2013) 222302. arXiv:1307.2539 , doi:10.1103/PhysRevLett.111.222302 .[38] P. Romatschke, Relativistic fluid dynamics far from local equilibrium,Phys. Rev. Lett. 120 (1) (2018) 012301. arXiv:1704.08699 , doi:10.1103/PhysRevLett.120.012301 .[39] P. Romatschke, Relativistic hydrodynamic attractors with broken sym-metries: non-conformal and non-homogeneous, JHEP 12 (2017) 079. arXiv:1710.03234 , doi:10.1007/JHEP12(2017)079 .[40] D. Almaalol, A. Kurkela, M. Strickland, Non-equilibrium attractorin high-temperature QCD plasmas, Phys. Rev. Lett. 125 (12) (2020)122302. arXiv:2004.05195 , doi:10.1103/PhysRevLett.125.122302 .[41] A. Bazavov, F. Karsch, S. Mukherjee, P. Petreczky, Hot-dense LatticeQCD: USQCD whitepaper 2018, Eur. Phys. J. A 55 (11) (2019) 194. arXiv:1904.09951 , doi:10.1140/epja/i2019-12922-0 .[42] M. McNelis, U. Heinz, Hydrodynamic generators in relativistic kinetictheory, Phys. Rev. C 101 (5) (2020) 054901. arXiv:2001.09125 , doi:10.1103/PhysRevC.101.054901 .[43] S. Jeon, U. Heinz, Introduction to hydrodynamics, Int. J. Mod.Phys. E24 (10) (2015) 1530010. arXiv:1503.03931 , doi:10.1142/S0218301315300106 .[44] L. Keegan, A. Kurkela, A. Mazeliauskas, D. Teaney, Initial conditions forhydrodynamics from weakly coupled pre-equilibrium evolution, JHEP08 (2016) 171. arXiv:1605.04287 , doi:10.1007/JHEP08(2016)171 .9045] A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting, D. Teaney,Matching the Nonequilibrium Initial Stage of Heavy Ion Collisions toHydrodynamics with QCD Kinetic Theory, Phys. Rev. Lett. 122 (12)(2019) 122302. arXiv:1805.01604 , doi:10.1103/PhysRevLett.122.122302 .[46] S. Kamata, M. Martinez, P. Plaschke, S. Ochsenfeld, S. Schlichting,Hydrodynamization and nonequilibrium Green’s functions in kinetictheory, Phys. Rev. D 102 (5) (2020) 056003. arXiv:2004.06751 , doi:10.1103/PhysRevD.102.056003 .[47] W. Florkowski, R. Ryblewski, M. Strickland, Testing viscous andanisotropic hydrodynamics in an exactly solvable case, Phys. Rev. C 88(2013) 024903. arXiv:1305.7234 , doi:10.1103/PhysRevC.88.024903 .[48] M. Martinez, M. McNelis, U. Heinz, Anisotropic fluid dynamics forGubser flow, Phys. Rev. C95 (5) (2017) 054907. arXiv:1703.10955 , doi:10.1103/PhysRevC.95.054907 .[49] D. P. Bazow, Fluid dynamics for the anisotropically expanding quark-gluon plasma, Ph.D. thesis, Ohio State U. (2017).URL http://rave.ohiolink.edu/etdc/view?acc_num=osu1493056954338129 [50] J. Liu, C. Shen, U. Heinz, Pre-equilibrium evolution effects on heavy-ion collision observables, Phys. Rev. C91 (6) (2015) 064906, [Erratum:Phys. Rev.C92,no.4,049904(2015)]. arXiv:1504.02160 , doi:10.1103/PhysRevC.92.049904,10.1103/PhysRevC.91.064906 .[51] A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting, D. Teaney,Effective kinetic description of event-by-event pre-equilibrium dynamicsin high-energy heavy-ion collisions, Phys. Rev. C 99 (3) (2019) 034910. arXiv:1805.00961 , doi:10.1103/PhysRevC.99.034910 .[52] J. Berges, M. P. Heller, A. Mazeliauskas, R. Venugopalan, Thermal-ization in QCD: theoretical approaches, phenomenological applications,and interdisciplinary connections, arXiv:2005.12299 .[53] T. Nunes da Silva, D. Chinellato, M. Hippert, W. Serenone, J. Taka-hashi, G. S. Denicol, M. Luzum, J. Noronha, Pre-hydrodynamic evo-91ution and its signatures in final-state heavy-ion observables, arXiv:2006.02324 .[54] C. Shen, B. Schenke, Dynamical initial state model for relativistic heavy-ion collisions, Phys. Rev. C97 (2018) 024907. arXiv:1710.00881 , doi:10.1103/PhysRevC.97.024907 .[55] C. Shen, B. Schenke, Dynamical initialization and hydrodynamic model-ing of relativistic heavy-ion collisions, Nucl. Phys. A 982 (2019) 411–414. arXiv:1807.05141 , doi:10.1016/j.nuclphysa.2018.08.007 .[56] Y. Akamatsu, et al., Dynamically integrated transport approach forheavy-ion collisions at high baryon density, Phys. Rev. C98 (2) (2018)024909. arXiv:1805.09024 , doi:10.1103/PhysRevC.98.024909 .[57] L. Du, U. Heinz, G. Vujanovic, Hybrid model with dynamical sources forheavy-ion collisions at BES energies, Nucl. Phys. A 982 (2019) 407–410. arXiv:1807.04721 , doi:10.1016/j.nuclphysa.2018.09.015 .[58] L. Du, U. Heinz, (3+1)-dimensional dissipative relativistic fluid dynam-ics at non-zero net baryon density, Comput. Phys. Commun. 251 (2020)107090. arXiv:1906.11181 , doi:10.1016/j.cpc.2019.107090 .[59] W. Florkowski, R. Ryblewski, Highly-anisotropic and strongly-dissipative hydrodynamics for early stages of relativistic heavy-ion colli-sions, Phys. Rev. C83 (2011) 034907. arXiv:1007.0130 , doi:10.1103/PhysRevC.83.034907 .[60] M. Martinez, M. Strickland, Dissipative dynamics of highly anisotropicsystems, Nucl. Phys. A848 (2010) 183–197. arXiv:1007.0889 , doi:10.1016/j.nuclphysa.2010.08.011 .[61] M. Martinez, R. Ryblewski, M. Strickland, Boost-invariant (2+1)-dimensional anisotropic hydrodynamics, Phys. Rev. C85 (2012) 064913. arXiv:1204.1473 , doi:10.1103/PhysRevC.85.064913 .[62] D. Bazow, U. Heinz, M. Strickland, Second-order (2+1)-dimensionalanisotropic hydrodynamics, Phys. Rev. C90 (5) (2014) 054910. arXiv:1311.6720 , doi:10.1103/PhysRevC.90.054910 .9263] D. Bazow, U. W. Heinz, M. Martinez, Nonconformal viscous anisotropichydrodynamics, Phys. Rev. C 91 (6) (2015) 064903. arXiv:1503.07443 , doi:10.1103/PhysRevC.91.064903 .[64] L. Tinti, Anisotropic matching principle for the hydrodynamic expan-sion, Phys. Rev. C94 (4) (2016) 044902. arXiv:1506.07164 , doi:10.1103/PhysRevC.94.044902 .[65] E. Molnar, H. Niemi, D. H. Rischke, Derivation of anisotropic dis-sipative fluid dynamics from the Boltzmann equation, Phys. Rev.D93 (11) (2016) 114025. arXiv:1602.00573 , doi:10.1103/PhysRevD.93.114025 .[66] M. McNelis, D. Bazow, U. Heinz, (3+1)-dimensional anisotropicfluid dynamics with a lattice QCD equation of state, Phys. Rev.C97 (5) (2018) 054912. arXiv:1803.01810 , doi:10.1103/PhysRevC.97.054912 .[67] G. S. Denicol, S. Jeon, C. Gale, Transport coefficients of bulk viscouspressure in the 14-moment approximation, Phys. Rev. C90 (2) (2014)024912. arXiv:1403.0962 , doi:10.1103/PhysRevC.90.024912 .[68] S. Ryu, J. F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon,C. Gale, Importance of the bulk viscosity of QCD in ultrarelativisticheavy-ion collisions, Phys. Rev. Lett. 115 (13) (2015) 132301. arXiv:1502.01675 , doi:10.1103/PhysRevLett.115.132301 .[69] E. Molnar, H. Niemi, D. H. Rischke, Closing the equations of motionof anisotropic fluid dynamics by a judicious choice of a moment of theBoltzmann equation, Phys. Rev. D94 (12) (2016) 125003. arXiv:1606.09019 , doi:10.1103/PhysRevD.94.125003 .[70] M. Nopoush, R. Ryblewski, M. Strickland, Anisotropic hydrodynamicsfor conformal Gubser flow, Phys. Rev. D 91 (4) (2015) 045007. arXiv:1410.6790 , doi:10.1103/PhysRevD.91.045007 .[71] M. Alqahtani, M. Nopoush, R. Ryblewski, M. Strickland, (3+1)D quasi-particle anisotropic hydrodynamics for ultrarelativistic heavy-ion col-lisions, Phys. Rev. Lett. 119 (4) (2017) 042301. arXiv:1703.05808 , doi:10.1103/PhysRevLett.119.042301 .9372] M. Alqahtani, M. Nopoush, R. Ryblewski, M. Strickland, Anisotropichydrodynamic modeling of 2.76 TeV Pb-Pb collisions, Phys. Rev.C96 (4) (2017) 044910. arXiv:1705.10191 , doi:10.1103/PhysRevC.96.044910 .[73] D. Almaalol, M. Alqahtani, M. Strickland, Anisotropic hydrodynamicmodeling of 200 GeV Au-Au collisions, Phys. Rev. C 99 (4) (2019)044902. arXiv:1807.04337 , doi:10.1103/PhysRevC.99.044902 .[74] M. McNelis, D. Everett, U. Heinz, Particlization in fluid dynamicalsimulations of heavy-ion collisions: The iS3D module, Comput. Phys.Commun. 258 (2021) 107604. arXiv:1912.08271 , doi:10.1016/j.cpc.2020.107604 .[75] A. Kurganov, E. Tadmor, New high-resolution central schemes for non-linear conservation laws and convection–diffusion equations, J. Comput.Phys. 160 (1) (2000) 241–282. doi:10.1006/jcph.2000.6459 .URL http://dx.doi.org/10.1006/jcph.2000.6459 [76] L.-G. Pang, H. Petersen, X.-N. Wang, Pseudorapidity distribution anddecorrelation of anisotropic flow within the open-computing-languageimplementation CLVisc hydrodynamics, Phys. Rev. C97 (2018) 064918. arXiv:1802.04449 , doi:10.1103/PhysRevC.97.064918 .[77] S. Jaiswal, C. Chattopadhyay, A. Jaiswal, S. Pal, U. Heinz, Exact solu-tions and attractors of higher-order viscous fluid dynamics for Bjorkenflow, Phys. Rev. C 100 (3) (2019) 034901. arXiv:1907.07965 , doi:10.1103/PhysRevC.100.034901 .[78] A. Kurkela, W. van der Schee, U. A. Wiedemann, B. Wu, Early-and Late-Time Behavior of Attractors in Heavy-Ion Collisions, Phys.Rev. Lett. 124 (10) (2020) 102301. arXiv:1907.08101 , doi:10.1103/PhysRevLett.124.102301 .[79] L. Tinti, A. Jaiswal, R. Ryblewski, Quasiparticle second-order viscoushydrodynamics from kinetic theory, Phys. Rev. D95 (5) (2017) 054007. arXiv:1612.07329 , doi:10.1103/PhysRevD.95.054007 .[80] A. Bazavov, et al., Equation of state in (2+1)-flavor QCD, Phys. Rev.D90 (2014) 094503. arXiv:1407.6387 , doi:10.1103/PhysRevD.90.094503 . 9481] J. Weil, et al., Particle production and equilibrium properties withina new hadron transport approach for heavy-ion collisions, Phys. Rev.C94 (2016) 054905. arXiv:1606.06642 , doi:10.1103/PhysRevC.94.054905 .[82] H. Marrochio, J. Noronha, G. S. Denicol, M. Luzum, S. Jeon, C. Gale,Solutions of conformal Israel-Stewart relativistic viscous fluid dynamics,Phys. Rev. C 91 (1) (2015) 014903. arXiv:1307.6130 , doi:10.1103/PhysRevC.91.014903 .[83] G. S. Denicol, C. Gale, S. Jeon, A. Monnai, B. Schenke, C. Shen, Netbaryon diffusion in fluid dynamic simulations of relativistic heavy-ioncollisions, Phys. Rev. C 98 (3) (2018) 034916. arXiv:1804.10557 , doi:10.1103/PhysRevC.98.034916 .[84] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numeri-cal Recipes in C (2nd Ed.): The Art of Scientific Computing, CambridgeUniversity Press, New York, NY, USA, 1992.[85] P. Huovinen, H. Petersen, Particlization in hybrid models, Eur.Phys. J. A48 (2012) 171. arXiv:1206.3371 , doi:10.1140/epja/i2012-12171-9 .[86] J. D. Bjorken, Highly relativistic nucleus-nucleus collisions: the cen-tral rapidity region, Phys. Rev. D27 (1983) 140–151. doi:10.1103/PhysRevD.27.140 .[87] S. S. Gubser, A. Yarom, Conformal hydrodynamics in Minkowski and deSitter spacetimes, Nucl. Phys. B846 (2011) 469–511. arXiv:1012.1314 , doi:10.1016/j.nuclphysb.2011.01.012 .[88] S. S. Gubser, Symmetry constraints on generalizations of Bjorken flow,Phys. Rev. D82 (2010) 085027. arXiv:1006.0006 , doi:10.1103/PhysRevD.82.085027 .[89] G. S. Denicol, U. W. Heinz, M. Martinez, J. Noronha, M. Strickland,Studying the validity of relativistic hydrodynamics with a new exactsolution of the Boltzmann equation, Phys. Rev. D 90 (12) (2014) 125026. arXiv:1408.7048 , doi:10.1103/PhysRevD.90.125026 .9590] W. Florkowski, R. Ryblewski, M. Strickland, Anisotropic hydrodynam-ics for rapidly expanding systems, Nucl. Phys. A 916 (2013) 249–259. arXiv:1304.0665 , doi:10.1016/j.nuclphysa.2013.08.004 .[91] A. Gron, Hands-On Machine Learning with Scikit-Learn and Tensor-Flow: Concepts, Tools, and Techniques to Build Intelligent Systems,1st Edition, O’Reilly Media, Inc., 2017.[92] P. M. Chesler, L. G. Yaffe, Holography and colliding gravitational shockwaves in asymptotically AdS spacetime, Phys. Rev. Lett. 106 (2011)021601. arXiv:1011.3562 , doi:10.1103/PhysRevLett.106.021601 .[93] P. M. Chesler, L. G. Yaffe, Numerical solution of gravitational dynamicsin asymptotically anti-de Sitter spacetimes, JHEP 07 (2014) 086. arXiv:1309.1439 , doi:10.1007/JHEP07(2014)086 .[94] P. M. Chesler, Colliding shock waves and hydrodynamics in small sys-tems, Phys. Rev. Lett. 115 (24) (2015) 241602. arXiv:1506.02209 , doi:10.1103/PhysRevLett.115.241602 .[95] M. Attems, J. Casalderrey-Solana, D. Mateos, D. Santos-Oliv´an, C. F.Sopuerta, M. Triana, M. Zilh˜ao, Holographic collisions in non-conformaltheories, JHEP 01 (2017) 026. arXiv:1604.06439 , doi:10.1007/JHEP01(2017)026 .[96] M. P. Heller, A. Kurkela, M. Spali´nski, V. Svensson, Hydrodynamiza-tion in kinetic theory: Transient modes and the gradient expansion,Phys. Rev. D 97 (9) (2018) 091503. arXiv:1609.04803 , doi:10.1103/PhysRevD.97.091503 .[97] F. Cooper, G. Frye, Single-particle distribution in the hydrodynamicand statistical thermodynamic models of multiparticle production, Phys.Rev. D 10 (1974) 186–189. doi:10.1103/PhysRevD.10.186 .[98] A. Monnai, Dissipative hydrodynamic effects on baryon stopping, Phys.Rev. C 86 (2012) 014908. arXiv:1204.4713 , doi:10.1103/PhysRevC.86.014908 .[99] Ohio Supercomputer Center (1987).URL http://osc.edu/ark:/19495/f5s1ph73http://osc.edu/ark:/19495/f5s1ph73