Anisotropic fluid dynamical simulations of heavy-ion collisions
Abstract
We present VAH, a (3+1)–dimensional simulation that evolves the far-from-equilibrium quark-gluon plasma produced in ultrarelativistic heavy-ion collisions with anisotropic fluid dynamics. We solve the hydrodynamic equations on an Eulerian grid using the Kurganov–Tadmor algorithm in combination with a new adaptive Runge–Kutta method. Our numerical scheme allows us to start the simulation soon after the nuclear collision, largely avoiding the need to integrate it with a separate pre-equilibrium dynamics module. We test the code’s performance by simulating on the Eulerian grid conformal and 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 hydrodynamics in central Pb+Pb collisions and find that the former’s longitudinal flow profile responds more consistently to the fluid’s gradients along the spacetime rapidity direction.
keywords:
Ultrarelativistic heavy-ion collisions , quark-gluon plasma , relativistic hydrodynamics , computational fluid dynamicsPROGRAM 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: 1.2 GB (for a grid)
Keywords: Ultrarelativistic heavy-ion collisions, quark-gluon plasma, relativistic hydrodynamics, 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 in ultrarelativistic 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+Pb collision on a grid takes about 530 for an Intel Xeon E5-2680 v4 multi-core processor with OpenMP acceleration (see Sec. 5.3).
1 Introduction
The quark-gluon plasma is one of the most extreme phases of matter in our universe. The incredibly high temperatures and densities required to create it ( K, kg/m3) are orders of magnitude beyond those of any substance typically encountered in nature. The quark-gluon plasma also has the remarkable property of having the lowest shear viscosity to entropy density ratio of any known fluid; current phenomenological constraints put at temperatures around the deconfinement temperature MeV [1, 2]. This value is very close to the conjectured Kovtun-Son-Starinet bound , which follows from the AdS/CFT correspondence in string theory [3, 4]. Therefore, it is of great interest to the physics community to test whether or not the KSS bound holds for such strongly coupled liquids. Today, we can produce quark-gluon plasma by colliding beams of heavy ions (e.g. Pb+Pb) at very high energies. The femtoscopically small droplets of plasma formed in these ultrarelativistic nuclear collisions have incredibly short lifetimes ( s), making them difficult to probe directly. To reconstruct the medium’s transport properties, one can simulate the various stages of a heavy-ion collision with quantitative precision and predictive power and fit the theoretical model to large sets of soft-momentum ( GeV) hadronic observables [5, 6, 7, 1, 2, 8, 9].111High-energy jets, heavy quarks and electromagnetic radiation can also serve as experimental probes but their full integration into hybrid models [10, 11, 12, 13, 14, 15] is more involved than for soft-momentum hadrons. The hybrid model must be accurate enough that it can make full use of the considerable precision of the available experimental data.
The anisotropically expanding quark-gluon plasma stage is usually modeled with relativistic viscous hydrodynamics [16, 17, 18, 5, 19, 20]. Conventionally, viscous hydrodynamics only works for fluids that are both near local equilibrium (i.e. inverse Reynolds number ) and have small spacetime gradients (i.e. Knudsen number Kn ) [18, 21]. But unlike most fluids found in nature, the quark-gluon plasma is initially far from local-equilibrium () and sustains moderately large gradients (Kn ) for a good portion of its lifetime. This is primarily due to quantum fluctuations in the initial-state profile [22, 23, 24] and the rapid longitudinal expansion at early times. Under these extreme conditions, the validity of a fluid dynamical description for the quark-gluon plasma comes into question [25]. Nevertheless, second-order viscous hydrodynamic simulations, in conjunction with other multi-stage modules, have been widely successful at reproducing hadronic observables (e.g. anisotropic flow coefficients ) [26, 27]. State-of-the-art hybrid models now have enough predictive power to quantitatively constrain the shear and bulk viscosities of the quark-gluon plasma using heavy-ion experimental data and Bayesian inference [6, 7, 1, 2, 8, 9]. Hydrodynamics has also been found to be applicable in collisions between light and heavy ions (e.g. He3+Au, p+Pb) and in proton–proton collisions [12, 28], although the origin of collectivity in these small systems is still being debated [29, 30, 31, 32].
The success of second-order viscous hydrodynamics can be attributed to an underlying resummed hydrodynamic theory that also extends to large gradients [33, 34, 35]. As the associated non-hydrodynamic modes decay on microscopic time scales ,222Hydrodynamic simulations often employ relaxation-type equations from relativistic kinetic theory (e.g. DNMR [18]), but they likely do not precisely capture the transient dynamics in strongly coupled fluids [36, 37]. the system approaches a (generally non-equilibrium) hydrodynamic attractor that is well approximated by viscous hydrodynamics. This happens within the time scale333The hydrodynamization time refers to the time when viscous hydrodynamics becomes applicable, but the hydrodynamic simulation starting time can be different from this. fm/, long before the system thermalizes [38, 39, 40]. Whether or not the strongly coupled quark-gluon plasma possesses an attractor at presently experimentally accessible temperatures ( GeV) is currently unknown. A major obstacle to answering this question are the technical difficulties involved in evaluating its transport properties from first principles, including its shear and bulk viscosities [41]. Recently, it was demonstrated that resummed hydrodynamics can be obtained by expanding the microscopic Green’s function of a linearized kinetic system, an alternative to a direct resummation of the gradient expansion [42]. It might be possible to extend this concept to the non-equilibrium quark-gluon plasma by perturbing the system around local-equilibrium and expanding the resulting correlation function(s) [43, 44, 45, 46].
Despite the robustness of second-order viscous hydrodynamics, it is still susceptible to breaking down when gradients are very large (Kn ) [47, 48, 20]. The largest gradients in heavy-ion collisions occur at very early times ( fm/), especially around the edges of the fireball. If one starts the viscous hydrodynamic simulation too early, one usually encounters negative longitudinal pressures due to huge pressure anisotropies [49]; near the edges of the fireball even the transverse pressure can turn negative if the bulk viscosity peaks strongly near the quark-hadron phase transition. Not only does this cause excessive viscous heating444Viscous heating refers to internal entropy production by dissipative processes, not due to external thermal sources. but large negative pressures also redirect the matter inward, potentially resulting in cavitation. Regulation schemes can be used to tamp down the viscous pressures and prevent the simulation from crashing, but they cloud the physical predictions of the original hydrodynamic theory [5, 20].555Some regulation schemes are less extreme than others. For example, the hydrodynamic code MUSIC only regulates the dissipative currents outside the fireball region [16] while iEBE-VISHNU regulates the entire grid [5]. Because of the technical issues involved, it is more suitable to use a pre-equilibrium dynamics model before transitioning to viscous hydrodynamics at [22, 26, 50, 6, 7, 1, 2, 44, 51, 45, 52]. Still, the conformal approximation usually made in the former results in a mismatch to the latter’s non-conformal equation of state, producing (in spite of the system’s expansion) artificially positive bulk viscous pressures that can be as large as around the edges of the fireball [53]. In some situations, one cannot even switch between dynamical models instantaneously. For example, in low-energy collisions ( GeV) where the nuclear interpenetration time is comparable to the fireball lifetime, one needs to run viscous hydrodynamics in the background as soon as the participant nucleons start feeding thermal energy and net-baryon number into the newly formed fireball [54, 55, 56, 57, 58].
The shortcomings of second-order viscous hydrodynamics have motivated the development of hydrodynamic models that can better handle far-from-equilibrium situations [59, 60, 61]. The most promising candidate is anisotropic hydrodynamics, which treats the two largest dissipative effects arising in heavy-ion collisions (the pressure anisotropy and the bulk viscous pressure ) non-perturbatively [62, 63, 64, 65, 66]. While current formulations of anisotropic hydrodynamics are partially based on weakly coupled kinetic theory,666Second-order viscous hydrodynamic simulations also rely on microscopic approaches such as relativistic kinetic theory to compute the relaxation times and other second-order transport coefficients [18, 67, 68]. they are able to capture exact kinetic solutions more accurately than standard viscous hydrodynamics [47, 62, 63, 64, 69, 70, 48]. In addition, their ability to maintain positive longitudinal and transverse pressures makes them less prone to cavitation in (3+1)–dimensional simulations [49]. This means they can run at early times with little interference from viscous regulations and remain numerically stable.
The application of anisotropic hydrodynamics to heavy-ion phenomenology is still relatively new; Au+Au collisions at RHIC ( GeV) and Pb+Pb collisions at the LHC ( and 5.02 TeV) have been modeled reasonably well but so far only using smooth initial conditions and a few model parameters to adjust the fit to experimental data [71, 72, 73]. Anisotropic hydrodynamics has yet to undergo the same extensive tests and trials as viscous hydrodynamics, but it can potentially serve as an additional discrete model for the fluid dynamical stage of heavy-ion collisions. This would further increase the flexibility of the Bayesian framework developed by the JETSCAPE collaboration [1, 2], which has already incorporated a number of discrete models for the particlization stage [74]. The hope is that, by controlling large dissipative flows nonperturbatively, anisotropic hydrodynamics can constrain the transport coefficients of QCD matter more accurately.
In this paper we introduce our (3+1)–dimensional anisotropic fluid dynamical simulation for heavy-ion collisions called VAH.777The code package can be downloaded from the GitHub repository https://github.com/mjmcnelis/cpu_vah. The C++ module is based on the GPU–accelerated viscous hydrodynamic code GPU VH [20, 49], except that it implements anisotropic hydrodynamics with the (, )–matching scheme developed in Ref. [66].888At this point VAH itself has not yet been ported to GPUs. We evolve the dynamical equations on an Eulerian grid using the Kurganov–Tadmor algorithm [75], a popular method also used in other viscous hydrodynamic codes [16, 17, 20, 76]. The main improvement in our code is the ability to automatically adjust the time step after each iteration, as opposed to using a fixed value for the entire simulation. Our adaptive Runge–Kutta scheme initially uses a fine time step to resolve the rapid longitudinal expansion at early times while speeding up the evolution at later times with a coarser time step given by the Courant-Friedrichs-Lewy (CFL) condition. As a result, we can start anisotropic hydrodynamics at a very early time fm/ to model both the far-off-equilibrium dynamics stage999We assume the system is longitudinally free-streaming (i.e. and ) in the time interval before starting anisotropic hydrodynamics. This closely mirrors the situation found in other pre-hydrodynamic models [22, 50, 45]. It has also been argued [77, 78] that (at least in weakly coupled systems) for the far-off-equilibrium hydrodynamic attractor approaches the free-streaming attractor in Bjorken flow (which approximates the early evolution stage in heavy-ion collisions). with a non-conformal QCD equation of state and smoothly transition to viscous hydrodynamics.101010Anisotropic hydrodynamics reduces to second-order viscous hydrodynamics in the limit of small Knudsen and inverse Reynolds numbers (). We do not switch to a separate viscous hydrodynamics model as the Knudsen and inverse Reynolds numbers decrease, but use anisotropic hydrodynamics to evolve both the earliest far-off-equilibrium and subsequent viscous hydrodynamic stages.
The code package contains other useful features, such as the option to run either anisotropic hydrodynamics or second-order viscous hydrodynamics within the same framework. This provides the user with several discrete models to evolve the fluid dynamical stage for comparison. We further implemented user options to accelerate the simulation on a multi-core processor with OpenMP, and to use automated grid settings that optimize the overall size of the spatial grid for a given set of runtime parameters (e.g. the impact parameter and particlization switching temperature).
This paper is organized as follows: In Sec. 2 we review the anisotropic hydrodynamic equations used to evolve the energy-momentum tensor of the quark-gluon plasma. Sec. 3 discusses the numerical implementation of these dynamical equations in the code. In Sec. 4 we test our anisotropic fluid dynamical simulation for a system subject to (non)conformal Bjorken flow and conformal Gubser flow. Furthermore, we compare anisotropic hydrodynamics to second-order viscous hydrodynamics in (3+1)–dimensions. Finally, we benchmark the typical runtimes of (2+1)–d and (3+1)–d non-conformal hydrodynamic simulations in Sec. 5.
In this work we use Milne spacetime coordinates , where is the longitudinal proper time and is the spacetime rapidity. We adopt the mostly-minus convention for the metric tensor, . In the code, we solve the hydrodynamic equations in natural units ; energy, momentum and temperature have units of fm-1, energy density and pressure have units of fm-4, etc. These units are converted to physical units (e.g. energy densities in GeV/fm3) when outputting and plotting the results. The net-baryon density and baryon diffusion current are neglected in this version of the code.
2 Anisotropic fluid dynamics
In this section we provide an overview of the anisotropic hydrodynamics equations, including the equation of state and transport coefficients, to be implemented in the code that evolves the energy-momentum tensor of the quark-gluon plasma. For more details on the derivation of these hydrodynamic equations we refer the reader to Refs. [66, 65].
2.1 Energy-momentum tensor
In anisotropic fluid dynamics, it is convenient to decompose in the basis {, , }, where the fluid velocity represents the temporal direction in the local fluid rest frame (LRF) and the spatial vectors are split into the longitudinal direction and the transverse projection tensor [65]. In the Landau frame this results in the decomposition (suppressing the spacetime dependence)
| (1) |
where round parentheses denote symmetrization, i.e. . The major components of are the LRF energy density , the longitudinal pressure and the transverse pressure . Together, the pressures and capture the largest dissipative flows in heavy-ion collisions: the pressure anisotropy caused by the rapid longitudinal expansion rate at early longitudinal proper times, and the bulk viscous pressure111111Here is the average (isotropic) pressure and is the equilibrium pressure. due to critical fluctuations near the quark-hadron phase transition. The remaining components of are the longitudinal momentum diffusion current and the transverse shear stress tensor , with being the traceless double transverse projector. We refer to these components as residual shear stresses since they are typically smaller than the pressure anisotropy term in the full shear stress tensor
| (2) |
2.2 Dynamical variables
The dynamical variables that we propagate in the code are
| (3) |
Their evolution equations will be discussed below. Although and each have only two independent components, we evolve all 14 of their components independently to simplify the workflow of the algorithm.121212When propagating these extraneous components numerically, slight violations of the orthogonality and tracelessness conditions (74) can occur. We correct for these errors at each step of the simulation with the regulation scheme described in Sec. 3.5. In addition, we propagate the energy density and the fluid velocity’s spatial components131313The fluid velocity’s temporal component is . (, , ) since they appear in the hydrodynamic equations. They are inferred from the components :
| (4) |
To solve these equations, one also needs to know . From the orthogonality conditions and , there are only two nonzero components that depend on the fluid velocity:
| (5) |
where is the transverse velocity. The solution of the inferred variables (, ) from the algebraic equations (4) will be discussed in Sec. 3.3.
For (3+1)–dimensionally expanding fluids, we have a total of dynamical variables and four inferred variables to evolve on an Eulerian grid.141414Dynamical variables are evolved directly using the Kurganov–Tadmor algorithm (see Sec. 3.1). Inferred variables are determined from the dynamical variables algebraically.,151515For anisotropic fluid dynamics with a QCD equation of state [66], we further evolve the mean-field as a dynamical variable and the anisotropic variables (, , ) as additional inferred variables (see Secs. 2.4 and 2.6.3). If the system is longitudinally boost-invariant, we do not need to propagate the components , , and since they vanish by symmetry.
2.3 Conservation laws
The evolution of the components is given by the energy-momentum conservation laws
| (6) |
where the covariant derivative accounts for the curvilinear nature of the Milne coordinates. The conservation equations can be expanded as
| (7) |
where is the partial derivative and are the Christoffel symbols. In Milne spacetime, the only nonzero Christoffel symbols are
| (8) |
Thus, the set of equations (7) can be rewritten as
| (9a) | ||||
| (9b) | ||||
| (9c) | ||||
| (9d) | ||||
where the Latin indices are summed over spatial components. We can eliminate the components in (9a) by using the identity:
| (10) |
Here we introduced the three-velocity as well as the tensors and . Likewise, we can express the components in (9b-d) in terms of :
| (11) |
After some algebra one obtains [49]
| (12a) | ||||
| (12b) | ||||
| (12c) | ||||
| (12d) | ||||
For the numerical algorithm the hydrodynamic equations must be written in conservative flux form [75]:161616For a conformal fluid, the energy density’s spatial derivatives are needed to evaluate the source term .
| (13) |
Here are the currents, are the source terms and is a spatial index. Naturally, the evolution equations for already assume this form. In the next subsection we will see that the relaxation equations for the dissipative flows require additional manipulations.
2.4 Relaxation equations
The relaxation equations for the dissipative flows , , and are [66, 65]
| (14a) | ||||
| (14b) | ||||
| (14c) | ||||
| (14d) | ||||
Here a dot above any quantity denotes the co-moving time derivative , and curly brackets denote either the transverse projection of a vector, , or the traceless double transverse projection of a rank-2 tensor, . We also define the longitudinal and transverse expansion rates and , where is the LRF longitudinal derivative and is the transverse gradient. The transverse velocity-shear tensor is , while the transverse vorticity tensor is given by , where the square brackets denote anti-symmetrization. The shear and bulk viscosity to entropy density ratios and and their associated relaxation times and , along with the anisotropic transport coefficients coupled to the gradient forces, will be discussed in Sec. 2.6.
We recast the relaxation equations in conservative flux form by using the product rule identities
| (15a) | ||||
| (15b) | ||||
to rewrite the l.h.s. of Eqs. (14a-d) as
| (16a) | |||
| (16b) | |||
| (16c) | |||
| (16d) | |||
where is the fluid acceleration. Finally, we use the product rule identity
| (17) |
to rewrite Eqs. (16a-b) as
| (18a) | ||||
| (18b) | ||||
here
| (19a) | ||||
| (19b) | ||||
are the gradient source terms for and . Similarly Eqs. (16c-d) can be rewritten as
| (20a) | ||||
| (20b) | ||||
here
| (21a) | ||||
| (21b) | ||||
are the gradient source terms,
| (22a) | ||||
| (22b) | ||||
are the transverse projection source terms, and
| (23a) | ||||
| (23b) | ||||
are the geometric source terms for and . The individual components that make up the source terms in the relaxation equations (18a-b) and (20a-b) are listed in Appendices A and B.
2.5 Equation of state
In the code there are two options for the quark-gluon plasma’s equation of state : conformal and QCD. The former assumes a non-interacting gas of massless quarks and gluons:
| (26) |
where is the temperature and the degeneracy factor is
| (27) |
with colors and massless quark flavors (i.e. up, down and strange). This equation of state is primarily used to test the fluid dynamical simulation subject to either conformal Bjorken expansion or conformal Gubser expansion (see Secs. 4.1 and 4.2).
The QCD equation of state that we employ for more realistic simulations interpolates between the lattice QCD calculations provided by the HotQCD collaboration [80] and a hadron resonance gas composed of the hadrons that can be propagated in the hadronic afterburner code SMASH [81].171717In hybrid model simulations of heavy-ion collisions, the equilibrium equation of state used in the hydrodynamic module must be consistent with that in the hadronic afterburner. Otherwise, serious violations of energy and momentum conservation can occur at the hadronization phase. Fig. 1 shows the energy density, equilibrium pressure and speed of sound as a function of temperature. The data table in the code covers the temperature range GeV.
2.6 Transport coefficients
In this section we list the transport coefficients appearing in the relaxation equations in Sec. 2.4, for both QCD and conformal equations of state.
At the present moment, the transport coefficients of QCD matter (especially in the nonperturbative region GeV) are not explicitly known from first principles. Instead, we parametrize the shear and bulk viscosities and as a function of temperature. In this work, we use the best-fit parametrization models from the JETSCAPE collaboration [1, 2]. The relaxation times and anisotropic transport coefficients are computed with a quasiparticle kinetic theory model, whose equation of state is fitted to the QCD one. The kinetic model contains quasiparticles with a temperature-dependent mass and an equilibrium mean-field ; these are shown in Figure 2.
2.6.1 Shear and bulk viscosities
The shear viscosity is modeled as a linear piecewise function with a kink at temperature :
| (28) |
where is the value of at , and are the left and right slopes, respectively, and is the Heaviside step function. The bulk viscosity is parametrized as a skewed Cauchy distribution:
| (29) |
where is the normalization factor, is the peak temperature and
| (30) |
with and being the width and skewness parameters, respectively.
The best-fit values for the viscosity parameters181818They correspond to a hybrid model whose particlization phase uses the 14-moment approximation for the correction in the Cooper-Frye formula [74, 1, 2] (the viscosity parameters used for the code validation tests in Sec. 4 are slightly different because they were taken from an early draft of Ref. [2]). are , GeV, GeV-1, GeV-1, , GeV, GeV and (see Table II in Ref. [2]). The resulting temperature dependence of and is shown in Figure 3.
For conformal systems, we fix the shear viscosity to and the bulk viscosity to .
2.6.2 Shear and bulk relaxation times
In the quasiparticle kinetic model [79, 66], the shear and bulk relaxation times are proportional to the shear and bulk viscosity, respectively,
| (31) |
where the viscosity to relaxation time ratios and are
| (32a) | ||||
| (32b) | ||||
and and are evaluated with the QCD equation of state. Here, we have defined the thermodynamic integrals
| (33) |
where indicates integration with the Lorentz-invariant momentum space measure, is the quasiparticle momentum, is the square of its spatial LRF momentum, is its LRF energy, and is the local-equilibrium distribution function.
The normalized quasiparticle relaxation times and are shown in Figure 4. These are compared to the relaxation times in standard viscous hydrodynamic models, where the kinetic transport coefficients are computed in the small-mass approximation and :
| (34a) | ||||
| (34b) | ||||
One sees that the shear relaxation times are very similar to each other, except at low temperatures GeV. On the other hand, the bulk relaxation times differ by about an order of magnitude for GeV; this is due to the breakdown of the small-mass approximation, even at high temperatures GeV. As a result, the evolution of the bulk viscous pressure will be more greatly affected by critical slowing down in our anisotropic hydrodynamics model compared to standard viscous hydrodynamics.
For conformal kinetic plasmas (, ), the shear relaxation time is and the bulk relaxation time is .
2.6.3 Anisotropic transport coefficients
Finally, we list the anisotropic transport coefficients [66, 65] that are coupled to the gradient source terms in the relaxation equations for the longitudinal pressure ,
| (35a) | ||||
| (35b) | ||||
| (35c) | ||||
| (35d) | ||||
| (35e) | ||||
the transverse pressure ,
| (36a) | ||||
| (36b) | ||||
| (36c) | ||||
| (36d) | ||||
| (36e) | ||||
the longitudinal momentum diffusion current ,
| (37a) | ||||
| (37b) | ||||
| (37c) | ||||
| (37d) | ||||
| (37e) | ||||
| (37f) | ||||
| (37g) | ||||
| (37h) | ||||
and the transverse shear stress tensor :
| (38a) | ||||
| (38b) | ||||
| (38c) | ||||
| (38d) | ||||
| (38e) | ||||
| (38f) | ||||
We define the anisotropic integrals
| (39) |
where
| (40) |
and is the leading-order anisotropic distribution function.
In addition to the effective temperature this anisotropic distribution depends (through ) on two momentum anisotropy parameters and , which deform the longitudinal and transverse momentum space, respectively. In order to compute the non-conformal anisotropic transport coefficients at each time step of the simulation, we propagate as inferred variables. They are obtained from the kinetic part of , and in the quasiparticle kinetic model:
| (41a) | ||||
| (41b) | ||||
| (41c) | ||||
where is the kinetic contribution to the energy density (i.e. the total energy density minus the mean field contribution [66]), is the kinetic longitudinal pressure, and is the kinetic transverse pressure. The numerical method which we use to solve these equations will be discussed in Sec. 3.4.
The anisotropic transport coefficients in the conformal limit are listed in Appendix C.
3 Numerical scheme
In this section we discuss the numerical implementation of the hydrodynamic equations in the code. The dynamical and inferred variables are evolved on an Eulerian grid, where , and are the number of physical grid points along each spatial direction.191919For longitudinally boost-invariant systems, the number of spacetime rapidity points is set to . A grid point with cell index that corresponds to the lower left front corner of a fluid cell, has a spatial position
| (42a) | ||||
| (42b) | ||||
| (42c) | ||||
where , and are the lattice spacings. Physical fluid cells have indices , and . The numerical algorithm also requires six sets of ghost cells with depth two, which neighbor the physical grid’s faces (for an illustration, see Fig. 3 in Ref. [20]). The boundary conditions that we impose on the ghost cells neighboring the two faces are202020In conformal anisotropic hydrodynamics, the energy density also requires ghost cell boundary conditions.
| (43a) | ||||
| (43b) | ||||
| (43c) | ||||
| (43d) | ||||
and similarly for and faces after permuting the grid indices and replacing or .
The dynamical variables 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 in the RK2 scheme, we reconstruct the inferred variables from the dynamical variables. In addition, we regulate the residual shear stresses and and the mean-field . The code also provides the user with the option of using an adaptive time step to capture the fluid’s longitudinal dynamics at very early times; this will be discussed at the end of the section.
3.1 Kurganov-Tadmor algorithm
The time derivative of the dynamical variables in the partial differential equations (13) can be computed on the Eulerian grid at any time using the KT algorithm [75]:
| (44) |
where the numerical fluxes evaluated at the left and right faces of a staggered cell centered around the grid point are
| (45) |
and similarly for and after permuting the in the grid indices and the corresponding spatial component (i.e. or ).
The first two terms in Eq. (45) take the average of the currents extrapolated to the staggered cell face (or ) from the left () and right () sides. A first-order expression for the extrapolated currents can be computed using the chain rule:
| (46a) | ||||
| (46b) | ||||
| (46c) | ||||
| (46d) | ||||
where . Similarly, the extrapolated currents , , and at the remaining faces of the staggered cell are obtained by permuting the (or ) in the grid indices, the spatial components and derivatives, and the lattice spacing.
The final term in Eq. (45) takes into account the wave propagation of the discontinuities at a finite speed [75]. We define the local propagation speed component at the staggered cell faces as
| (47) |
where the extrapolated velocities are
| (48a) | ||||
| (48b) | ||||
| (48c) | ||||
| (48d) | ||||
The discontinuities propagating from the staggered cell faces depend on the extrapolated dynamical variables
| (49a) | ||||
| (49b) | ||||
| (49c) | ||||
| (49d) | ||||
The formulae for the local propagation speed components and , extrapolated velocities , , and , and extrapolated dynamical variables , , and are analogous.
The numerical spatial derivatives appearing in the extrapolated quantities (46a-d), (48a-d) and (49a-d) are computed with a minmod flux limiter [75]:
| (50a) | ||||
| (50b) | ||||
where
| (51) |
with
| (52) |
the flux limiter parameter is set to [82]. The flux limiter derivatives , , and are analogous.
In contrast, the numerical spatial derivatives appearing in the source terms are approximated with second-order central differences [20, 76]:
| (53a) | ||||
| (53b) | ||||
and similarly for , , and . The fluid velocity’s time derivative also appears in the source terms; its evaluation will be discussed the next subsection.
3.2 Two-stage Runge–Kutta scheme
Because the KT algorithm (44) admits a semi-discrete form [75] (i.e. the time derivative is continuous while the spatial derivatives are discrete), we can combine it with an RK2 ODE solver to evolve the system in time [16, 20]. Given the dynamical variables at time (discrete times are labeled with index ), we evolve the system one time step with an intermediate Euler step (omitting the spatial indices):
| (54) |
where is evaluated with r.h.s of Eq. (44). The time derivative in the source terms is approximated with a first-order backward difference [20, 17]:
| (55) |
where is the previous fluid velocity and is the previous time step.212121At the start of the hydrodynamic simulation, or , the previous time step is set to the current time step . Unless stated otherwise, we also initialize the previous fluid velocity as . From the intermediate variables (54), we reconstruct the inferred variables as well as the anisotropic variables ( , ), as described in Secs. 3.3 and 3.4 below. Afterwards, we regulate the mean field and residual shear stresses and in (see Sec. 3.5) and set the ghost cell boundary conditions for and .
Next, we evolve the system with a second intermediate Euler step
| (56) |
where the fluid velocity’s time derivative is now evaluated as
| (57) |
In the RK2 scheme, we average the two intermediate Euler steps in to update the dynamical variables at :
| (58) |
From this, we update the inferred variables (, ) and (, , ) and regulate the residual shear stresses and mean field. Finally, we set the ghost cell boundary conditions for and and proceed with the next RK2 iteration.
3.3 Reconstructing the energy density and fluid velocity
Given the hydrodynamic variables , along with , , and , we can reconstruct the energy density and fluid velocity from Eq. (4). The solution for the energy density is [66]
| (59) |
where was defined earlier and
| (60) |
The reconstruction formula (59) also requires the components , , which can be expressed in terms of hydrodynamic variables as follows:222222Writing implies that the argument in Eq. (61) is always positive.
| (61) |
Here
| (62) |
with and where .
In the cold, dilute regions surrounding the fireball (and, occasionally, in cold spots within the fluctuating fireball), the energy density can become much smaller than the freezeout energy density .232323In this work, we construct a particlization hypersurface of constant energy density GeV/fm3, which corresponds to the switching temperature GeV. The lowest value used for the switching temperature in the JETSCAPE SIMS analysis is GeV [1, 2]. While these regions are not phenomenologically important, hydrodynamic simulations are susceptible to crashing there without intervention [5, 20, 83]. To prevent this, we regulate the energy density with the formula
| (63) |
where is the minimum energy density allowed in the Eulerian grid and . For , the regulation has virtually no effect on the energy density. As , however, the energy density is smoothly regulated to . Ideally, should be the lower limit of our QCD equation of state table GeV/fm3. However, we find that the code evolving non-conformal anisotropic hydrodynamics with fluctuating initial conditions encounters fewer technical difficulties if we instead use the larger value GeV/fm3, which is still about six times smaller than our choice for .242424The constraint is not imposed for conformal systems. We checked that the regulation scheme (63) has little to no impact on the fluid’s dynamics in regions where or larger.
After regulating the energy density, we reconstruct the fluid velocity’s spatial components:
| (64a) | ||||
| (64b) | ||||
| (64c) | ||||
Because of the regulation (63), there will be inconsistencies between the inferred variables (, ) and the components but only in the dilute cold regions.
3.4 Reconstructing the anisotropic variables
To compute the non-conformal anisotropic transport coefficients in Sec. 2.6.3, we need to reconstruct the anisotropic variables
| (65) |
from , , , and , by solving the system of equations [66]
| (66) |
where
| (67) |
We solve these nonlinear equations numerically using Newton’s method. Taking the anisotropic variables prior to the intermediate Euler step (54) as our initial guess252525At we initialize the anisotropic variables (, , ) by iterating (where inside the parentheses is the temperature). , we iterate along the direction given by the Newton step
| (68) |
where is the inverse of the Jacobian
| (69) |
Specifically, we iterate the solution as
| (70) |
where is a partial step that is optimized with a line-backtracking algorithm to improve the global convergence of Newton’s method [84]. This procedure is repeated until has converged to the solution of Eq. (66).262626In our simulation tests, we found several instances when Newton’s method failed to converge near the edges of the Eulerian grid. Whenever that happens we simply set the anisotropic variables to the initial guess .
For conformal systems (, , ), the anisotropic variables and are much easier to solve. The system of equations (67) reduces to
| (71a) | ||||
| (71b) | ||||
where the functions and are listed in Appendix C. We numerically invert the longitudinal pressure to energy density ratio for :
| (72) |
From this, we can evaluate as
| (73) |
Because the solutions (72) and (73) do not require the previous values for and as input, we do not need to store them in memory during runtime.
3.5 Regulating the residual shear stresses and mean-field
After reconstructing the inferred variables, we regulate the residual shear stress components and such that they satisfy the orthogonality and traceless conditions
| (74a) | ||||
| (74b) | ||||
| (74c) | ||||
| (74d) | ||||
| (74e) | ||||
and that their overall magnitude is smaller than the longitudinal and transverse pressures:
| (75) |
The former ensures that and maintain their orthogonal and tracelessness properties within numerical accuracy during the hydrodynamic simulation. The latter prevents the residual shear stresses from overwhelming the anisotropic part of the energy-momentum tensor [49].
First, we update these components in the following order:
| (76a) | ||||
| (76b) | ||||
| (76c) | ||||
| (76d) | ||||
| (76e) | ||||
| (76f) | ||||
| (76g) | ||||
| (76h) | ||||
| (76i) | ||||
| (76j) | ||||
while leaving the components , , and unchanged. Next, we rescale all the components by the same factor :
| (77a) | ||||
| (77b) | ||||
where272727For longitudinally boost-invariant systems, we replace the second argument in Eq. (78) by .
| (78) |
with . In the validation tests discussed in Sec. 4, we did not encounter a situation where the residual shear stresses are suppressed by Eq. (78). This indicates that the residual shear stresses are naturally smaller 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, we find that the relaxation equation (24) for the mean-field can become unstable in certain spacetime regions ( GeV), causing the bulk viscous pressure to grow positive without bound. The origin of this issue is not fully understood, and we will address it in future work. For now, we remove this instability by regulating the nonequilibrium mean-field component as
| (79) |
where
| (80) |
In practice, we only find it necessary to regulate the mean-field when (i.e. in regions with ).
3.6 Adaptive time step
Most relativistic hydrodynamic codes that simulate heavy-ion collisions evolve the system with a fixed time step [16, 5]. Any choice for must satisfy the CFL condition so that the hydrodynamic simulation is at least dynamically stable [75]. However, the time step must also be small enough to resolve the fluid’s evolution rate (in particular, the large longitudinal expansion rate at early times) but, on the other hand, large enough to finish the simulation within a reasonable runtime. This balancing act places a practical limit on how early the user can start the hydrodynamic simulation (in practice, the smallest value typically used is fm/ [26]). This is not much of a concern for second-order viscous hydrodynamics since it is anyhow prone to breaking down for very early initialization times, as discussed in Sec. 1. Anisotropic hydrodynamics, on the other hand, can handle the large pressure anisotropies occurring at early times much better; to realize its full potential it should be initialized at earlier times, but for that we need to move away from a fixed time step. In this section, we introduce a new adaptive stepsize method, which automatically adjusts the successive time step to be larger or smaller than in such a way that we can push back our fluid dynamical simulation to very early times () without sacrificing numerical accuracy nor computational efficiency.282828This 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 that much longer than the hydrodynamization time (see Sec. 1).
For a system with no source terms () the KT algorithm is dynamically stable as long as the time step satisfies the CFL condition [75]
| (81) |
where () are the maximum local propagation speed components on the Eulerian grid at time . In Milne spacetime, the quark-gluon plasma’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) allows one to maintain dynamical stability with an adaptive time step that is generally larger than the fixed time step obtained by taking the limit :
| (82) |
For , however, the time step from (81) is too coarse to resolve the fluid’s gradients and relaxation rates at early times fm/ when the flow profile is very nonrelativistic. Therefore, at early times when the source terms are strongest, the adaptive time step should primarily depend on those source terms. The following implementation ensures that, as the source terms relax and the fluid velocity grows more relativistic over time, the adaptive time step naturally approaches the CFL bound (81).
First, we consider a homogeneous fluid undergoing Bjorken expansion (i.e. and ). The system of differential equations (13) for the dynamical variables reduces to
| (83) |
and we are given the initial values and time step at time . We evolve the system one time step using the RK2 scheme (58):
| (84) |
For the next iteration, we determine how much we need to adjust the time step . Adaptive stepsize methods do this by estimating the local truncation error of each time step [84]. In our method, we approximate the local truncation error of the next intermediate Euler step
| (85) |
which is
| (86) |
where denotes the -norm. If the second time derivative at is known, we can set the local truncation error to the desired error tolerance (after dropping higher-order terms)
| (87) |
where is the error tolerance parameter and is the number of dynamical variables. It is reasonable to use absolute or relative errors when is small or large, respectively [84]. In this work, we set the error tolerance parameter to .
To obtain an expression for the second time derivative, we compute the next intermediate Euler step (85) using the old time step (denoted by ):
| (88) |
This allows us to approximate with central differences:
| (89) |
Compared to the usual central difference formula there is an additional factor of that accounts for the local truncation error present in . Furthermore, the expression (89) is numerically accurate to rather than . After substituting Eqs. (85) and (89) in Eq. (87), one has for the next time step
| (90) |
where
| (91) |
and is the numerical solution to the algebraic equation
| (92) |
with .
For the general case without Bjorken symmetry, varies across the grid. We then perform the calculation (90) (replacing by ) for all spatial grid points and take the minimum value. Afterwards, we place safety bounds to prevent the time step from changing too rapidly:
| (93) |
where we set the control parameter to . Finally, we impose the CFL bound (81) on . With the new time step at hand, we resume computing the next RK2 iteration. Notice that there are no additional numerical evaluations of the flux and source terms in the adaptive RK2 scheme since we can recompute the next intermediate Euler step simply by adjusting the time step in . This allows our adaptive time step algorithm to be readily 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 .292929We do not allow the adaptive time step to become any smaller than this (i.e. we impose ). At first, the adaptive time step tends to increase with time since the fluid’s longitudinal expansion rate decreases. As the transverse flow builds up, eventually becomes bounded by the CFL condition (81).
3.7 Program summary
We close this section by summarizing the workflow of the anisotropic fluid dynamical simulation with a QCD equation of state. Figure 5 shows the flowchart of the code’s primary mode,303030The secondary (test) mode outputs the hydrodynamic evolution from the simulation and their corresponding semi-analytic solutions, if they exist. which constructs a hypersurface of constant energy density 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 the Eulerian grid at a given time . Altogether, there are three dynamical variables (q, qI, Q), two fluid velocity variables (u, up), one energy density 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 or .
-
(b)
qI holds the intermediate dynamical variables (or ).
-
(c)
Q holds the previous, updated or current dynamical variables , or .
-
(d)
u holds the current, intermediate or updated fluid velocity , or .
-
(e)
up holds the previous or current fluid velocity or .
-
(f)
e holds the current, intermediate or updated energy density , or .
-
(g)
lambda, aT and aL hold the current, intermediate or updated anisotropic variables (, , ), (, , ) or (, , )
-
(a)
-
3.
We initialize the variables q, u, up, e, lambda, aT and aL as follows:313131If the code is run with Gubser initial conditions, the hydrodynamic variables are initialized differently (see Sec. 4.2). first, we compute (or read in) the energy density profile given by an initial-state model (e.g. TENTo [24]). Then we initialize the longitudinal and transverse pressures as
(94a) (94b) where is a pressure anisotropy ratio parameter set by the user (typically a small value ). This model assumes that only the pressure anisotropy enters in the initial and while the initial bulk viscous pressure is . The initial fluid velocity is static (i.e. ) and the residual shear stresses and are set to zero. From this, we compute the components with Eq. (4). Finally, we set the mean-field to and initialize the anisotropic variables by solving Eq. (66).323232This initialization scheme is similar to how conformal free-streaming modules are initialized at [50, 6, 7, 1, 2].
-
4.
We set the ghost cell boundary conditions for q, e and u.333333The ghost cells of e are only used in conformal anisotropic hydrodynamics to approximate the source terms on the grid’s faces. We also configure the freezeout finder and set the initial time step to .
-
5.
We start the hydrodynamic evolution at (or ):
-
(a)
We call the freezeout finder every two time steps343434The user can adjust the number of time steps between freezeout finder calls. (i.e. ) to search for freezeout cells between the Eulerian grids from the current and previous calls (if , we only load the initial grid to the freezeout finder). The hypersurface volume elements and their spacetime positions are constructed with the code CORNELIUS [85]. The hydrodynamic variables at the hypersurface elements’ centroid are approximated with a 4d linear interpolation.
- (b)
-
(c)
We reconstruct the intermediate inferred variables , , , and from qI and store the results in e, u, lambda, aT and aL, respectively. We regulate the residual shear stresses and mean-field in qI and set the ghost cell boundary conditions for qI, e and u.
- (d)
-
(e)
We update the inferred variables , , , and 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 cells in the grid are below the switching temperature (i.e. ).
-
(a)
-
6.
After the hydrodynamic evolution, we deallocate the hydrodynamic variables and store the freezeout surface in memory, which can either be passed to another program or written to file.
4 Validation tests and hydrodynamic model comparisons
In this section we test the validity of our anisotropic fluid dynamical simulation for various initial-state configurations, using either the conformal or QCD equation of state. First, we run anisotropic hydrodynamics with conformal Bjorken and Gubser initial conditions [86, 87, 88], whose semi-analytic solutions can be computed accurately using a fourth-order Runge–Kutta ODE solver [69, 48]. After these two validation tests, we compare (3+1)–dimensional conformal anisotropic hydrodynamics and second-order viscous hydrodynamics353535The VAH code can also run second-order viscous hydrodynamics, where the kinetic transport coefficients are computed with either quasiparticle masses (see Fig. 2a) [79, 66] or light masses (i.e. ) [67, 20] (we use the same shear and bulk viscosities as in Fig. 3). In this work, we label the former model as quasiparticle viscous hydrodynamics (or VH) and the latter as standard viscous hydrodynamics (or VH2); for conformal systems, these two models are equivalent. For details on the numerical implementation of second-order viscous hydrodynamics, see Appendix E. for a central Pb+Pb collision with a smooth TENTo initial condition [24].
Next, we run non-conformal anisotropic hydrodynamics and viscous hydrodynamics 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 hydrodynamics in central Pb+Pb collisions with smooth or fluctuating TENTo initial conditions, with the goal of identifying situations where VAH offers definitive advantages in reliability over standard viscous hydrodynamics.
4.1 Conformal Bjorken flow test
For the first test, we evolve a conformal plasma undergoing Bjorken expansion [86]. In Bjorken flow, the system is longitudinally boost-invariant and homogeneous in the transverse plane. As a result, the fluid velocity is static and the residual shear stresses and vanish. The only two degrees of freedom are the energy density and longitudinal pressure (the transverse pressure is ). The anisotropic hydrodynamic equations simplify to [69]
| (95a) | ||||
| (95b) | ||||
where the conformal transport coefficient is given in Eq. (35a).
We start the simulation363636Bjorken flow can be simulated on either one fluid cell (i.e. or a larger grid with homogeneous initial conditions. at fm/ with an initial temperature GeV, pressure ratio and shear viscosity . We evolve the system with an adaptive time step, initially set to fm/, until the temperature drops below GeV at fm/. Figure 6 shows the evolution of the energy density normalized to its initial value and the pressure ratio . At very early times, the system is approximately free-streaming due to the rapid longitudinal expansion, resulting in a large Knudsen number and causing the energy density to decrease like . Over time, the longitudinal expansion rate decreases and the system approaches local equilibrium, i.e. and .
We compare the (0+1)–d VAH simulation to the semi-analytic solution of (95), which uses a fixed time step fm/. The simulation is in excellent agreement with the semi-analytic solution, with numerical errors staying below . We also checked the convergence of the simulation curves as we decrease the error tolerance parameter in the adaptive time step algorithm. Figure 7 shows the evolution of the adaptive time step for the parameters and . One sees that the time step increases linearly with time and is not bounded by the CFL condition (81) since the flow is stationary (i.e. . As a result, the adaptive RK2 scheme quickly reaches the switching temperature in 145 time steps, compared to about 40000 steps for the semi-analytic method.
4.2 Conformal Gubser flow test
In the next test, we evolve a conformal fluid subject to Gubser flow [88, 87], which is longitudinally boost-invariant and azimuthally symmetric. The fluid velocity profile takes the following analytic form:
| (96) |
where
| (97) |
with being the transverse radius, the azimuthal angle, and an inverse length scale parameter that determines the fireball’s transverse size at .
One can transform the polar Milne coordinates to the de Sitter coordinates where
| (98a) | ||||
| (98b) | ||||
The line element in this de Sitter space possesses symmetry [88, 87]:
| (99) |
which makes the fluid velocity stationary.373737A hat denotes a quantity in the de Sitter space . All hatted quantities are made unitless by multiplying with appropriate powers of . Hence, the hydrodynamic variables only depend on the de Sitter time , and the anisotropic hydrodynamic equations reduce to [48]
| (100a) | ||||
| (100b) | ||||
The residual shear stresses and vanish under Gubser symmetry [48]. In the semi-analytic solution (100), we start the evolution at , which corresponds to a corner in a fm fm transverse grid ( fm) at the initial time fm/, with an initial temperature and pressure ratio . We set the inverse fireball size to fm-1 and the shear viscosity to . We evolve the system with a fixed stepsize until , which corresponds to the center of our Eulerian grid ( fm) at fm/.
Next, we map the semi-analytic solution for the de Sitter energy density and longitudinal pressure back to Milne coordinates [89]:
| (101) |
We use this map to set the initial energy density and longitudinal pressure profiles and ; the initial temperature at the center of the fireball is GeV. We set the initial transverse shear stress to . Finally, we use Eq. (96) to initialize the current fluid velocity and previous fluid velocity , where the initial time step is fm/. We evolve the Gubser profile on a fm fm transverse grid with a lattice spacing of fm.
Figure 8 shows the evolution of the energy density and transverse fluid velocity in the transverse plane at various time frames.383838In order to output the hydrodynamic quantities at specific times, we readjust the adaptive time step whenever we are close to these output times (this is not shown in Fig. 10). At early times, the very hot and compact fireball expands rapidly along the longitudinal direction and quickly cools down without much change to its transverse shape. Over time, the transverse expansion overtakes the longitudinal expansion, pushing the medium radially outward as a ring of fire. One sees that the energy density and transverse fluid velocity maintain their azimuthal symmetry throughout the evolution. However, there are some numerical fluctuations around the lines at later times, especially near the peak of . This is a consequence of using a Cartesian grid with a finite spatial resolution.
Figure 9 shows the evolution of the energy density , fluid velocity component , pressure anisotropy ratio and transverse shear inverse Reynolds number / along the –axis at multiple time frames. Overall, there is very good agreement between the simulation and semi-analytic solution, except near the local extrema of , where the numerical errors are about ; these errors can be reduced by using a finer lattice spacing. The transverse shear stress should remain zero throughout the simulation. However, the transverse shear inverse Reynolds number from the VAH simulation is nonzero (albeit small, ) since the Eulerian grid does not perfectly preserve Gubser symmetry.
Figure 10 shows the evolution of the adaptive time step in the Gubser simulation. In contrast to Fig. 7, here becomes bounded by the CFL condition (81) at fm/ due to the increasing transverse expansion rate. As a result, the Gubser simulation finishes in 407 time steps, compared to 480 steps if we had used the fixed time step . Although we only gain a slight speedup, the adaptive time step is able to resolve the fluid’s early-time dynamics more accurately than the fixed time step (82) because it is initially independent of the lattice spacing. This property is especially useful when simulating more realistic nuclear collisions, where the lattice spacing required to resolve the participant nucleons’ transverse energy deposition is several times coarser than the one used in this test (e.g. fm).
4.3 (3+1)–d conformal hydrodynamic models in central Pb+Pb collisions
Here we compare (3+1)–dimensional conformal anisotropic hydrodynamics VAH to conformal second-order viscous hydrodynamics VH (see Appendix E), for a central Pb+Pb collision at zero impact parameter ( fm) with static (in Milne coordinates) and smooth initial conditions. To generate the latter we use an azimuthally symmetric TENTo energy density profile averaged over 2000 fluctuating events [24] and extended along the –direction with a smooth rapidity plateau [76] (see Appendix D for more information on the TENTo energy deposition model for high-energy nuclear collisions). The initial temperature at the center of the fireball is GeV at a starting time of which is the same for both types of hydrodynamic simulations. The fluid is evolved from an initial pressure ratio with a constant specific shear viscosity ; the runs stop once all the fluid cells are below GeV, which happens after fm/.
Figure 11 shows the evolution of , , and / along the –axis (. Initially, the pressure anisotropy in second-order viscous hydrodynamics is so large that the longitudinal pressure turns negative, especially near the grid’s edges. In comparison, the pressure ratio in anisotropic hydrodynamics is, at early times, only moderately larger in the central fireball region but quite dramatically different near its transverse edge, staying positive everywhere. The resulting differences in early-time viscous heating create a disparity between the two models for the normalization of the energy density which persists to late times even after the ratios have converged. On the other hand, the transverse flows predicted by the two models are nearly identical even at very early times since they are primarily driven by the transverse pressure gradients which, in the smooth TENTo profile, are smaller at early times than the longitudinal gradients . The transverse shear stress , which tends to counteract the fluid’s transverse acceleration, is larger in viscous hydrodynamics than in anisotropic hydrodynamics, especially along the edges of the fireball.393939The transverse shear stress is generally nonzero even for azimuthally symmetric flow profiles as long as they do not possess Gubser symmetry. Overall, however, the hydrodynamic variables are not substantially different along the transverse directions.
The longitudinal profile, however, evolves very differently in anisotropic hydrodynamics (VAH) compared to second-order viscous hydrodynamics (VH). This is shown in Fig. 12 where we plot the evolution of the dimensionless longitudinal velocity (as well as and ) along the –axis (). Similar to Fig. 11, the shear stress (2) in viscous hydrodynamics quickly overpowers the equilibrium pressure, making negative initially. Along the longitudinal direction this causes a strong reversal of the longitudinal flow at around . This unphysical feature (primarily driven by the strong longitudinal expansion rate at early times) results in a longitudinal contraction (in ) of the fireball near its forward and backward edges at the beginning of the simulation. In anisotropic hydrodynamics, the longitudinal pressure remains always positive, allowing for a stronger longitudinal flow profile whose signs are consistent with the gradients of the rapidity distribution of the energy density (142) (and thus of the thermal pressure). This gives anisotropic hydrodynamics a clear advantage over second-order viscous hydrodynamics when simulating the longitudinal dynamics of a heavy-ion collision.
In panel (d) of Fig. 12 we also plot the spacetime rapidity dependence of the inverse Reynolds number . The longitudinal momentum diffusion current is nonzero only in regions that have both longitudinal and transverse gradients. For this reason we shift the –axis transversely to fm, which corresponds to the mid-right region of the grid. As expected, the momentum diffusion current is weakest around mid-rapidity, , where the fireball profile is approximately longitudinally boost-invariant, and strongest along the sloping edges of the rapidity plateau (142). However, this longitudinal edge region is also the place where its inverse Reynolds number differs most strongly between anisotropic and standard viscous hydrodynamics. Overall, only makes up a tiny fraction of the total shear stress (2) in anisotropic hydrodynamics. The situation may change for initial conditions whose longitudinal and transverse gradients 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 hydrodynamic simulations. One sees that in VH stagnates until fm/ (see footnote 29) while the one in VAH starts increasing immediately. This indicates that at early times viscous hydrodynamics generally has a faster evolution rate and therefore requires a smaller initial time step compared to anisotropic hydrodynamics. Nevertheless, both adaptive time steps remain below their CFL bound (initially dominated by the longitudinal velocity ) until fm/. Notice that the CFL bounds for VH and VAH become virtually identical since they have almost the same transverse velocity profile (see Fig. 11b).
4.4 Non-conformal Bjorken flow test
Now we turn to testing our non-conformal hydrodynamic simulation with the QCD equation of state from Fig. 1. First, we run (3+1)–d anisotropic hydrodynamics with Bjorken initial conditions and compare it to the semi-analytic solution of the equations [66]
| (102a) | ||||
| (102b) | ||||
| (102c) | ||||
| (102d) | ||||
here we use the shear and bulk relaxation times (31)–(32), with the viscosities given by Eqs. (28)–(29) (see Fig. 3). The non-conformal anisotropic transport coefficients and are given by Eqs. (35a) and (36a), respectively. We start the simulation at fm/ with initial temperature GeV and initial pressure ratio , and evolve the system until fm/ when the temperature falls below GeV (we plot results only up to fm/). We also repeat this for quasiparticle viscous hydrodynamics (VH) and standard viscous hydrodynamics (VH2) (see footnote 35).404040VH and VH2 use the same equation of state as VAH but different transport coefficients, see Appendix E and Ref. [66]. Ideally, we would have preferred using smaller values for and to better match the longitudinally free-streaming initial condition used in the conformal hydrodynamic simulations (see Secs. 4.1–4.3 and footnote 9). However, we found that at earlier times the VAH simulation has greater difficulty reconstructing the anisotropic variables (, , ), especially directly at initialization. This indicates that our anisotropic hydrodynamic model [66] which integrates the QCD equation of state consistently with quasiparticle kinetic transport coefficients, has its limitations.
Figures 14a,b show the evolution of the normalized energy density and pressure ratio from the three hydrodynamic simulations. We also disentangle the longitudinal and transverse pressures’ viscous components and into the pressure anisotropy and bulk viscous pressure ; their evolutions relative to are shown in Figs. 14c,d. Although our initial condition for the pressure ratio is somewhat ad hoc, it quickly reaches its minimum value at fm/c, allowing for the energy density to closely follow its free-streaming limit at the beginning of the simulation. The ratio is typically larger in anisotropic hydrodynamics compared to the two viscous hydrodynamic models (until fm/). This is mainly due to differences in the higher-order transport coefficients associated with the pressure anisotropy (i.e. beyond and ) [62, 90, 63, 66]. Although in Bjorken flow the bulk viscous pressure is much smaller than the pressure anisotropy , it evolves very differently in standard viscous hydrodynamics VH2 compared to quasiparticle viscous hydrodynamics VH and anisotropic hydrodynamics VAH. Since the dimensionless bulk relaxation time in standard viscous hydrodynamics VH2 (see the dashed curve in Fig. 4b), the bulk viscous pressure quickly reaches its Navier-Stokes solution at fm/. In contrast, the bulk relaxation time used in anisotropic hydrodynamics VAH and quasiparticle viscous hydrodynamics VH is significantly larger (solid curve in Fig. 4b). As a consequence, it takes a much longer time for to relax to [79, 66].
One also sees that the simulations agree very well with their corresponding semi-analytic solutions (continuous lines in transparent color). Although we do not display them explicitly here, the numerical errors are small enough that we can unambiguously distinguish the different dynamics of the three hydrodynamic models in comparison studies discussed in the next subsections.
4.5 (3+1)–d non-conformal hydrodynamic models in central Pb+Pb collisions
4.5.1 Smooth TENTo initial conditions
Next we run non-conformal hydrodynamics for a central Pb+Pb collision with smooth, azimuthally symmetric TENTo initial conditions; we set the initial time to fm/ and the initial pressure ratio to . The initial energy density profile is almost identical to the one in Sec. 4.3, except the normalization is five times smaller, making the initial temperature at the center of the fireball GeV. We evolve the system until, at fm/, all fluid cells are below GeV.
Figure 15 shows, for the first fm/ of the simulation, the evolution of , , , /, / and given by the three hydrodynamic models along the –axis (). Here the fireball maintains higher temperatures for a longer duration than in conformal hydrodynamics since the QCD equilibrium pressure is much weaker than its Stefan–Boltzmann limit (26) (see Fig. 1a). But similar to Fig. 11a, the energy density in anisotropic hydrodynamics cools down at a slightly faster rate compared to viscous hydrodynamics. This is initially due to the moderately larger ratio in the central fireball region, which closely follows the Bjorken evolution in Fig. 14b. The viscous corrections to the longitudinal pressure increase as we move towards the edges of the fireball, but still remains positive in anisotropic hydrodynamics. The longitudinal pressure in standard viscous hydrodynamics, however, is strongly negative at early times (and also, to a lesser extent, for the quasiparticle case), although it does not significantly alter the fluid’s transverse dynamics. The ratios start to overlap at later times, mainly driven by the pressure anisotropies’ convergence in Fig. 15e.
Compared to Fig. 11b, the transverse flow in non-conformal hydrodynamics is significantly weaker due to the softer QCD equation of state. The dip 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 to the conformal case in Fig. 11d. Among the hydrodynamic models, we see that VH2 produces the smallest around the edges of the fireball since it has the largest negative bulk viscous pressure there. This is due to its ability to converge to the Navier-Stokes solution shortly after the simulation begins. In contrast, the bulk viscous pressure in VH hovers slightly above zero across the grid before falling down toward negative values at fm/. Even then, it hardly catches up to the bulk viscous pressure curves in VH2 since the differences between their bulk relaxation times grow as we move away from the center of the fireball. Overall, the bulk viscous pressure in VH is much weaker compared to VH2, which leads to a stronger transverse flow. Finally, the VAH simulation has a bulk viscous pressure that is qualitatively similar to VH (except at the edges of the grid) but slightly lags behind it. As a result, it generates the largest transverse flow out of the three simulations. The study here illustrates how the evolution of the bulk viscous pressure and transverse fluid velocity strongly depends on the choice of hydrodynamic model, specifically the temperature-dependent model for the bulk relaxation time . This is likely to have important implications for the phenomenological extraction of the bulk viscosity [1, 2].
Next, we compare the fireballs’ longitudinal evolution along the –axis () in Figure 16. When comparing the longitudinal evolution of the energy density and viscous pressures between the three hydrodynamic models, we note qualitatively similar features as already observed in Fig. 15 for their transverse evolution. While the transverse velocities shown in Fig. 15b varied between the models as a result of differences in the transverse pressures (mainly the bulk viscous pressures), here differences in their longitudinal pressures result in very different longitudinal flow profiles (Figure 16b). Most notably, the large negative in standard viscous hydrodynamics (VH2) initially causes to reverse sign very sharply in the forward and backward rapidity regions. Only by imposing strong regulations on the VH2 simulation is it able to recover the correct direction of longitudinal flow at later times as the gradients relax. The same breakdown at early times can also be seen in quasiparticle viscous hydrodynamics (VH) although there the situation is not nearly as bad. With VAH we can maintain positive longitudinal pressures so that the fireball can expand without imploding in regions where the longitudinal gradients are large. This greatly reduces the amount of regulation needed for the viscous pressures.
Finally, in Fig. 17 we plot the adaptive time step for each of the three hydrodynamic simulations. One sees that in anisotropic hydrodynamics (VAH) increases steadily until hitting the CFL bound at fm/. Its behavior is similar in quasiparticle viscous hydrodynamics (VH). The two CFL bounds are slightly separated due to differences between their transverse velocity profiles (see Fig. 15b). On the other hand, the adaptive time step in standard viscous hydrodynamics (VH2) does not pick up until fm/ (see footnote 29). Even then, it fluctuates up and down before reaching its CFL bound, which is higher than the other two bounds since its transverse flow is suppressed by a large bulk viscous pressure. Although the code runs faster per time step (see Table 1), standard viscous hydrodynamics takes considerably more time steps to evolve the fluid than the other two hydrodynamic models, suggesting that it has a harder time resolving the fireball evolution in the presence of large gradients.
4.5.2 Fluctuating TENTo initial conditions
Finally we study the differences among the three non-conformal hydrodynamic models for a central Pb+Pb collision with a fluctuating initial energy density profile, using the same model parameters as in Sec. 4.5.1. In the TENTo model used in this work, the energy density profile only has fluctuations in the transverse directions; the longitudinal profile is modeled with a finite rapidity plateau with smooth slopes at its longitudinal ends (see Appendix D).
Figure 18 shows the evolution of the fluctuating energy density profile in the transverse plane , along with the transverse slice of a particlization hypersurface of constant energy density GeV/fm3, for the three hydrodynamic simulations. Since anisotropic hydrodynamics generally has a smaller shear stress than viscous hydrodynamics (especially at early times), the transverse fluctuations across its fireball show the least dissipation or smearing. As a result, it can convert the initial-state eccentricities into anisotropic flow slightly more efficiently. The pressure anisotropy and bulk viscous pressure mainly influence the overall size of the fireball on the particlization hypersurface, affecting the final-state particle yields. The initially strong longitudinal expansion rapidly cools down the system and temporarily shrinks the fireball size at early times; a more positive longitudinal pressure helps cool the fireball even further. Compared to the two viscous hydrodynamic models, anisotropic hydrodynamics has a larger longitudinal pressure, which means its hypersurface has a slightly narrower waist. At later times, the transverse expansion overtakes the longitudinal expansion, increasing the fireball size. Although anisotropic hydrodynamics has the largest transverse flow, it also has the smallest bulk viscous pressure, enabling the fireball to cool faster and evaporate more quickly. This ultimately results in a smaller maximum fireball size. In contrast, standard viscous hydrodynamics has a much larger bulk viscous pressure. This extends the fireball’s lifetime by about fm relative to the one in anisotropic hydrodynamics, allowing it to grow larger in size.
In Figs. 19 and 20 we plot the inverse Reynolds numbers from the three hydrodynamic simulations in the plane at and in the plane at , respectively. Traditionally, the inverse Reynolds number measures the validity of the hydrodynamic expansion around local equilibrium [20, 5]. Since anisotropic hydrodynamics expands the energy-momentum tensor (1) around an anisotropic background (i.e. ), its validity is measured by the residual shear inverse Reynolds number [49]
| (103a) | ||||
as opposed to the shear and bulk inverse Reynolds numbers in second-order viscous hydrodynamics [20, 5]:
| (104) |
with . One sees that the residual shear inverse Reynolds numbers stay much smaller than unity inside the fireball, indicating that a first-order expansion in is sufficient to capture the residual shear corrections to the anisotropic hydrodynamic equations. While there is no need to regulate their strength here, the regulation scheme might be needed for collision events with sharper initial-state fluctuations.
We also plot the contribution of the nonequilibrium mean-field component to the bulk viscous pressure in anisotropic hydrodynamics:
| (105) |
We find that is moderately large and positive near the edges of the fireball but small and negative inside the fireball. However, it is necessary to regulate the mean-field in the latter region via Eq. (80) to prevent unstable growth. Specifically, the instability seems to originate near the quark-hadron phase transition, where the driving force proportional to the kinetic trace anomaly is at its strongest while the bulk relaxation rate is near a minimum. After applying the regulation (80) for a brief period, evolves freely from that point onward. While this initial instability is unfortunate, the overall impact of the non-equilibrium mean-field component on the longitudinal and transverse pressures inside the fireball region is limited.
In second-order viscous hydrodynamics, the shear inverse Reynolds number is large for a short period of time fm/ after the collision. The bulk inverse Reynolds number in standard viscous hydrodynamics is also very large near the peak of the bulk viscosity during the same time frame. To prevent the code from crashing, our regulation scheme suppresses both the shear stress and bulk viscous pressure (see Appendix E). Although it is still technically possible to run the heavy-ion simulation with standard viscous hydrodynamics, we cannot escape the effects of regulation around the base of the particlization hypersurface. In contrast, the viscous pressures in quasiparticle viscous hydrodynamics often do not require regulation414141Even in the absence of regulation, the longitudinal pressure in quasiparticle viscous hydrodynamics can still be negative in some regions at early times ( fm/) (e.g. see Figs. 15 and 16). since the bulk inverse Reynolds number is quite small.
We close this section with a comparison of the particlization hypersurfaces in the and planes. In Fig. 19 one sees that standard viscous hydrodynamics has the largest hypersurface along the –direction since its large bulk viscous pressure generates the most viscous heating. It also has the longest lifetime fm/, for the same reason. In contrast, it has a pretty narrow waist along the –direction because its longitudinal flow profile initially contracts the medium in the rapidity direction. In anisotropic hydrodynamics, the longitudinal flow begins transporting the medium away from the collision zone, following the direction of the longitudinal pressure gradients. As a result, its hypersurface has a wider waist than the one in standard viscous hydrodynamics (about larger for fm/) along the –direction.
5 Benchmarks
In this Section we benchmark the typical computational time needed to run (2+1)–d non-conformal hydrodynamic simulations of Pb+Pb collisions at LHC energies ( TeV) with fluctuating initial conditions and different values for the impact parameter and model parameters. We also benchmark the OpenMP–accelerated runtime of (3+1)–d non-conformal hydrodynamics for a central Pb+Pb collision ( fm) with the smooth TENTo initial condition and model parameters from Sec. 4.5.1.
| Mean time | Max time | Mean time per step | Mean # steps | |
| VAH | 204 (136) | 1840 (1630) | 0.639 (0.410) | 272 |
| VH | 75.0 (49.9) | 561 (320) | 0.183 (0.119) | 350 |
| VH2 | 89.7 (56.9) | 463 (358) | 0.175 (0.109) | 461 |
5.1 (2+1)–d non-conformal hydrodynamics on a fixed grid
In this test, we generate random samples for the impact parameter and Bayesian model parameters [1, 2]. The impact parameter is sampled from a piecewise linear distribution,
| (106) |
where we set fm for the lead nuclear radius. The model parameters
| (107) |
are each sampled from a distribution that is uniform within the finite intervals used in the JETSCAPE SIMS Bayesian analysis of heavy-ion collisions [1, 2].424242The TENTo initial condition model parameters (, , , , ), are defined in Appendix D [24] and the viscosity model parameters (, , , , , , , ) are discussed in Sec. 2.6.1. Finally, the switching temperature determines the particlization hypersurface of constant energy density .,434343We exclude the ratio between the dimensionless shear relaxation time and shear viscosity as a continuous model parameter [1, 2] and instead allow it to vary between discrete hydrodynamic models. The rapidity plateau model parameters and (see Appendix D) are also not considered in this benchmark test.
For each parameter sample we run (2+1)–d non-conformal hydrodynamic simulations with fluctuating TENTo initial conditions on a fm fm transverse grid ().444444All other runtime parameters (e.g. the initial time and pressure ratio ) are the same as the ones used in Sec. 4.5.,454545We found that of the anisotropic hydrodynamic simulations finish successfully (the second-order viscous hydrodynamic simulations had a success rate of ). For certain model parameter sets and/or fluctuating initial conditions, the code fails to reconstruct the anisotropic variables in the cold dilute regions, either because the maximum number of Newton iterations was exhausted or the kinetic longitudinal pressure turned negative. The resulting runtime statistics of the three hydrodynamic models are shown in Table 1. On average, it takes about to run each (2+1)–d viscous hydrodynamic simulation on the fixed grid, compared to for anisotropic hydrodynamics. The additional routine in Sec. 3.4 makes the runtime per step in anisotropic hydrodynamics about slower than in viscous hydrodynamics, but this is partially compensated by the fewer time steps required to finish each simulation.
Figure 21 shows the runtime distribution of the anisotropic hydrodynamic simulations on the fixed grid. One sees that some runs take longer than others, depending on the parameter values used. For example, a smaller impact parameter increases the participant nucleon multiplicity, which extends the fireball’s lifetime. The runtime is also very sensitive to the nucleon width parameter : halving the value for doubles the spatial resolution required464646In this work, we set the transverse lattice spacing to . to capture the spatial variations of the fluctuating TENTo profile (see Appendix 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 longer than the mean runtime (see Table 1).
5.2 (2+1)–d non-conformal hydrodynamics on an automated grid
The previous benchmark test was performed on a fixed transverse grid of fm fm. Although this grid is large enough to fit all possible fireball sizes produced in Pb+Pb collisions, evolving peripheral collisions () on it is computationally inefficient. A large grid is also unnecessary for certain Bayesian model parameter combinations that further decrease the fireball size. In this section, we present a regression algorithm that predicts the fireball radius for a given set of parameters . With this, we can automatically optimize the grid size to save computational time and memory.
To train the regression model, we generate training parameter samples arranged in the column vectors
| (108) |
to construct our feature matrix [91]. For each parameter sample, we run a single hydrodynamic simulation with a smooth, event-averaged TENTo initial condition and compute the mean fireball radius (again arranged in vector); this will be our target variable474747We define the fireball radius as the maximum transverse radius of the particlization hypersurface in a given hydrodynamic simulation.
| (109) |
A sample subset of the scattering matrix containing the five most important features is shown in Fig. 22. As expected, the fireball size strongly decreases with the impact parameter. There are also some moderately positive correlations between the fireball radius and initial condition parameters (e.g. increasing the nucleon width parameter increases the spread of the participant nucleons’ thickness function, see Appendix D, Eq. (139)). While the viscosity parameters have a much smaller effect on the fireball radius individually, the model’s fit to the train-validation data improves when we include all of them.
| VAH | VH | VH2 | |
| fm | fm | fm | |
Next, we fit a cubic polynomial Lasso regression model with standardization feature scaling [91] to the training data (108)–(109). The regularization parameter is chosen to minimize the root-mean-square error , averaged over a five-fold cross-validation [91]. The values for the cross-validated error and regularization parameter of the regression fit are listed in Table 2.
Finally, we rerun the simulations with fluctuating TENTo initial conditions, using the 200 test parameter samples from earlier. For each parameter sample, we launch the regression model to predict the mean fireball radius and set the transverse grid lengths to
| (110) |
where the margin parameter gives the fireball additional room to expand within the grid. Here we set fm not only for extra space but also to account for statistical fluctuations of the fireball radius in fluctuating TENTo events, which were not considered in the training routine.484848With enough computing resources, the user has the option to run multiple fluctuating hydrodynamic simulations per parameter sample to produce a statistical distribution for the fireball radius , which can be characterized with a mean radius and standard deviation . For more details, we refer the reader to https://github.com/mjmcnelis/fireball. We find that the automated grid has about a success rate of enclosing the fireball without touching its edges. Fig. 21 shows the anisotropic hydrodynamic runtime distribution on the automated grid. On average, the automated grid algorithm reduces the grid area by relative to the fixed grid and provides the simulation with an additional speedup (see Table 1).
While (3+1)–d hydrodynamic simulations stand to benefit the most from an optimized grid volume, we would need to implement more realistic longitudinal initial conditions than the model used in this work before retraining the regression model to predict (in addition to its transverse radius) the fireball’s average elongation along the –direction. There may also be additional parameters that characterize such fluctuating (3+1)–d initial-state models, depending on the collision system of interest. Since it takes much more computational resources to produce such a training data set, we leave the (3+1)–d automated grid for future work.
5.3 (3+1)–d non-conformal hydrodynamics with OpenMP acceleration
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)–d hydrodynamic simulations, which run typically about two orders of magnitude longer than (2+1)–d simulations [20].494949For even faster runtimes, one can parallelize relativistic hydrodynamics on a graphics processing unit [20, 76]. The present version of VAH does not yet offer this option. Fig. 23 shows the runtime of the (3+1)–d non-conformal hydrodynamic simulations from Sec. 4.5.1. The anisotropic hydrodynamic simulation takes about two and a half hours to finish on a single-core CPU; the runtime of central Pb+Pb collisions can be about an order of magnitude longer or shorter than this depending on the values used for the nucleon-width parameter and rapidity plateau parameters and (see Appendix D). On an Intel Xeon E5-2680 multi-core processor, we can significantly reduce the simulation time to about nine minutes with 28 CPU threads. However, we do not achieve a perfect speedup because some parts of the routine are not easily parallelizable (e.g. the construction of the particlization hypersurface).
6 Summary and outlook
In this work we developed a (3+1)–dimensional anisotropic fluid dynamical simulation that evolves both the pre-equilibrium and viscous hydrodynamic stages of ultrarelativistic nuclear collisions. We validated the code’s performance by reproducing the semi-analytic solutions of conformal and non-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 accurately capture the early-time dynamics of Bjorken and Gubser flow while finishing the simulation within a reasonable number of iterations. We also compared anisotropic hydrodynamics to two second-order viscous hydrodynamic models in central Pb+Pb collisions. Apart from the apparent sensitivity of the bulk viscous pressure evolution to the bulk relaxation time, the three hydrodynamic models have similar transverse dynamics, at least in the mid-rapidity region . However, the fluid’s longitudinal evolution varies greatly near the longitudinal edges of the fireball. We showed that the longitudinal pressure in (3+1)–dimensional anisotropic hydrodynamics stays positive even in the presence of large gradients at very early times, unlike in viscous hydrodynamics. This causes the fluid to initially expand outward along the spacetime rapidity direction, as expected from the outward-pointing longitudinal gradients of the thermal pressure, reducing the risk of cavitation at the beginning of the simulation. With the new development presented here, we can for the first time model even the very early pre-equilibrium dynamics stage at with a QCD equation of state, as opposed to the conformal equation of state implicit in other pre-equilibrium models [92, 93, 50, 94, 95, 96, 44, 45].
In the near term, we plan to run VAH in the JETSCAPE framework with the Maximum a Posteriori (MAP) model parameters extracted in Refs. [1, 2].505050Several model parameters associated with the conformal free-streaming module would not be used in our code, which replaces the free-streaming stage by anisotropic hydrodynamics. Since the JETSCAPE SIMS hybrid model [1, 2] uses conformal free-streaming and standard viscous hydrodynamics to evolve the pre-equilibrium and hydrodynamic stages, we will replace these two modules with our code and study the changes to the hadronic observables, such as the transverse momentum spectra and –differential anisotropic flows. A full analysis of all available data will need to wait for a Bayesian recalibration of the full evolution model using VAH as its hydrodynamic core.
Integration of VAH with the JETSCAPE framework also requires an updated version of the particlization module iS3D [74] with a new option that allows choosing the leading-order anisotropic distribution (plus residual shear corrections ) for the hadronic distribution in the Cooper–Frye formula [97, 74]. Once completed, this will allow us to investigate how different selections among a discrete set of hydrodynamic models affect both the theoretical description of heavy-ion experimental observables and the shear and bulk viscosity constraints inferred from such data–theory comparisons. In particular, it would be interesting to test the predictions from anisotropic hydrodynamics against a variety of initializations of second-order viscous hydrodynamics using different pre-hydrodynamic evolution models515151Anisotropic hydrodynamics (VAH) as presented here is not meant to be integrated with a pre-hydrodynamic module, especially one that uses a conformal approximation. This is mainly due to the technical difficulties in initializing the mean-field and anisotropic variables [66]. Available second-order viscous hydrodynamic codes can read in the energy-momentum tensor from a pre-hydrodynamic module and thus initialize the dynamical and inferred variables more easily, but this has not yet been implemented in the VAH code. (or no pre-hydrodynamic stage at all), such as the recent study done in Ref. [53].
The current code only evolves the fluid’s energy-momentum tensor components and ignores the effects of net-baryon density and diffusion. Second-order viscous hydrodynamic codes with nonzero and have already been practically implemented [98, 83], and others that initialize viscous hydrodynamics with dynamical sources [10, 54, 55, 57, 58] are currently under development. The BEShydro code [57, 58], in particular, shares a common ancestry [49] and therefore a number of similar features with VAH. The integration of both codes into a single framework, in order to utilize the adaptive time step for capturing and resolving the dynamical production of energy and baryon sources at the onset of low-energy nuclear collisions, offers itself as an interesting project. On the theoretical side, the effects of non-zero chemical potentials for net charge, baryon number and strangeness have not yet been considered in our formulation of anisotropic hydrodynamics [66]. An upgrade of 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 for assisting in the development of the automated grid algorithm during the Erdős Institute 2020 Data Science Boot Camp. Computational resources for the code validation, comparison and benchmark tests were provided by the Ohio Supercomputer Center under Project PAS0254 [99]. This work was supported by the National Science Foundation (NSF) within the framework of the JETSCAPE Collaboration under Award No. ACI-1550223. Additional partial support by the U.S. Department of Energy (DOE), Office of Science, Office for Nuclear Physics under Award No. DE-SC0004286 and within the framework 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 of the anisotropic hydrodynamic equations. The derivatives of the fluid velocity component are
| (111a) | ||||
| (111b) | ||||
| (111c) | ||||
| (111d) | ||||
The derivatives of the longitudinal basis vector are
| (112a) | ||||
| (112b) | ||||
| (112c) | ||||
| (112d) | ||||
| (112e) | ||||
| (112f) | ||||
| (112g) | ||||
| (112h) | ||||
The divergence of the spatial fluid velocity is
| (113) |
The longitudinal expansion rate is
| (114) |
The transverse expansion rate is
| (115) |
where
| (116) |
is the scalar expansion rate.
The components of the fluid acceleration
| (117) |
are
| (118a) | ||||
| (118b) | ||||
| (118c) | ||||
| (118d) | ||||
The components of the longitudinal vector’s comoving time derivative
| (119) |
are
| (120a) | ||||
| (120b) | ||||
The components of the fluid velocity’s longitudinal derivative
| (121) |
are
| (122a) | ||||
| (122b) | ||||
| (122c) | ||||
| (122d) | ||||
The transverse gradient of the fluid velocity projected along the longitudinal direction is
| (123) |
where the components of
| (124) |
are
| (125a) | ||||
| (125b) | ||||
| (125c) | ||||
| (125d) | ||||
The transverse velocity-shear tensor is525252In the code we obtain by applying the projector onto directly, rather than simplifying the expression (126).
| (126) |
where the components of
| (127) |
are
| (128a) | ||||
| (128b) | ||||
| (128c) | ||||
| (128d) | ||||
| (128e) | ||||
| (128f) | ||||
| (128g) | ||||
| (128h) | ||||
| (128i) | ||||
| (128j) | ||||
The transverse vorticity tensor is
| (129) |
where the components of
| (130) |
are
| (131a) | ||||
| (131b) | ||||
| (131c) | ||||
| (131d) | ||||
| (131e) | ||||
| (131f) | ||||
| (131g) | ||||
| (131h) | ||||
| (131i) | ||||
| (131j) | ||||
Appendix B Geometric source terms
Appendix C Conformal anisotropic transport coefficients
In the conformal limit , the anisotropic transport coefficients (35) – (38) only depend the functions , which reduce to
| (134) |
where , and are given by Eqs. (27), (72) and (73), respectively. We list the functions used in this work [66]:
| (135a) | ||||
| (135b) | ||||
| (135c) | ||||
| (135d) | ||||
| (135e) | ||||
| (135f) | ||||
| (135g) | ||||
| (135h) | ||||
| (135i) | ||||
| (135j) | ||||
| (135k) | ||||
where and .
Appendix D TENTo energy deposition model
In the TENTo model, the transverse energy deposition (GeV/fm2) of a single fluctuating nuclear collision event in the mid-rapidity region is [24]
| (136) |
where is the normalization parameter and
| (137) |
is the reduced nuclear thickness function, with being the geometric parameter. The nuclear thickness function of nucleus () is535353In this work, we only consider Pb+Pb collisions () at LHC energies TeV; the inelastic nucleon–nucleon cross section is set to fm-2.
| (138) |
where are the number of participant nucleons from nucleus (); the nucleon positions are sampled from a Woods–Saxon distribution under the constraint that each nucleon–nucleon pair in nucleus () maintains a minimum separation [24]. The participant nucleon’s thickness function is centered around its sampled transverse position (, ):
| (139) |
where is the nucleon width parameter.545454This parameter does not control the nucleon’s charge radius but rather the spread of thermal energy it deposits in the collision zone along the transverse directions. Furthermore, the multiplicity factor of each participant nucleon is sampled from the gamma distribution
| (140) |
where and is the standard deviation.
For a very brief period after the collision, we assume the system is longitudinally free-streaming and static in Milne coordinates (i.e. and ) so that . Therefore, we initialize the energy density profile of the hydrodynamic simulation at the starting time as
| (141) |
Here we also extend the transverse energy density profile along the spacetime rapidity direction with a smooth plateau distribution (unitless) [76]:
| (142) |
where is the plateau length, is the standard deviation of the half-Gaussian tails and is the Heaviside step function. The initial condition parameter values used in this work are GeV, , fm, fm, , and [76, 1, 2].
We set the lattice spacings to and to resolve the fluctuating energy density profile (141) (or event-averaged profile). For the grid size, we set the longitudinal length to to fit the rapidity plateau (142). The transverse lengths and are automatically configured by the algorithm described in Sec. 5.2.
Appendix E Numerical implementation of second-order viscous hydrodynamics
In this appendix, we summarize how second-order viscous hydrodynamics is implemented the code. We evolve viscous hydrodynamics with the same numerical algorithm discussed in Sec. 3 except the energy-momentum tensor (1) is decomposed as
| (143) |
where is the shear stress tensor and is the bulk viscous pressure. We also define the spatial projector and traceless double spatial projector . The corresponding dynamical variables
| (144) |
are propagated along with and (the mean-field and anisotropic variables (, , ) are not propagated). Although has only five independent components, we propagate all ten components in the simulation [20].555555For longitudinally boost-invariant systems, we do not propagate the shear stress components , and .
E.1 Hydrodynamic equations
Here we list the evolution equations for the dynamical variables (144) (we refer the reader to Refs. [20, 18] for details on their derivation):
| (145a) | ||||
| (145b) | ||||
| (145c) | ||||
| (145d) | ||||
| (145e) | ||||
| (145f) | ||||
where ,
| (146a) | ||||
| (146b) | ||||
are the gradient source terms for and and
| (147a) | ||||
| (147b) | ||||
are the spatial projection and geometric source terms for . We also define the velocity-shear tensor and vorticity tensor .
E.2 Transport coefficients
In quasiparticle viscous hydrodynamics [79] (i.e. from Fig. 2a), the relaxation times (, ) and first-order transport coefficients (, ) are given by Eqs. (31) – (32). The second-order transport coefficients are [66]
| (149a) | ||||
| (149b) | ||||
| (149c) | ||||
| (149d) | ||||
| (149e) | ||||
where
| (150a) | ||||
| (150b) | ||||
| (150c) | ||||
and the function are given by Eq. (33).
E.3 Reconstruction formulas for the energy density and fluid velocity
We reconstruct the energy density by solving the following nonlinear equation via Newton’s method [5, 20]:
| (153) |
where
| (154) |
with and being the QCD equation of state. Using the previous energy density for the initial guess, we iterate the solution to Eq. (153) as
| (155) |
where
| (156) |
with being the QCD speed of sound squared. We repeat the iteration (155) until we achieve sufficient convergence or the energy density falls below . If the bulk viscous pressure , we regulate it so that ; this allows us to solve for explicitly [5]:
| (157) |
Afterwards, we regulate the energy density via Eq. (63) and evaluate the fluid velocity components as
| (158a) | ||||
| (158b) | ||||
| (158c) | ||||
E.4 Regulating the shear stress and bulk viscous pressure
In this regulation scheme, we first adjust the shear stress components
| (159a) | ||||
| (159b) | ||||
| (159c) | ||||
| (159d) | ||||
| (159e) | ||||
so that satisfies the orthogonality and tracelessness conditions
| (160a) | ||||
| (160b) | ||||
Then we regulate the shear stress and bulk viscous pressure as
| (161a) | ||||
| (161b) | ||||
where
| (162) |
with . The regulation factor usually suppresses and around the edges of the fireball at early times fm/, especially in standard viscous hydrodynamics (e.g. see Figs. (19c) and (20c)).
References
- [1] D. Everett, et al., Phenomenological constraints on the transport properties of QCD matter with data-driven model averaging, arXiv:2010.03928.
- [2] D. Everett, et al., Multi-system Bayesian constraints on the transport coefficients 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 theories and 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, Applying Bayesian parameter estimation to relativistic heavy-ion collisions: simultaneous characterization of the initial state and quark-gluon plasma medium, Phys. Rev. C94 (2016) 024907. arXiv:1605.03954, doi:10.1103/PhysRevC.94.024907.
- [7] 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ürsoy, R. Snellings, A transverse momentum differential global analysis of heavy ion collisions, arXiv:2010.15130.
- [9] G. Nijs, W. van der Schee, U. Gürsoy, R. Snellings, A Bayesian analysis of heavy ion collisions with Trajectum, arXiv:2010.15134.
- [10] T. Hirano, P. Huovinen, K. Murase, Y. Nara, Integrated dynamical approach 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 and electromagnetic 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 to initializing hydrodynamic fields and mini-jet propagation in quark-gluon fluids, 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 dilepton tomography, 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üller, Coupled Boltzmann transport equations of heavy quarks and quarkonia in quark-gluon plasma, arXiv:2004.06746.
- [16] B. Schenke, S. Jeon, C. Gale, (3+1)D hydrodynamic simulation of relativistic 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 transient 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-ion collisions, Int. J. Mod. Phys. A28 (2013) 1340011. arXiv:1301.5893, doi:10.1142/S0217751X13400113.
- [20] D. Bazow, U. Heinz, M. Strickland, Massively parallel simulations of relativistic 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 University Press, London, 2013. doi:10.1093/acprof:oso/9780198528906.001.0001.
- [22] B. Schenke, P. Tribedy, R. Venugopalan, Fluctuating Glasma initial conditions 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 PHOBOS Glauber 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 to wounded nucleon and binary collision scaling in high-energy nuclear collisions, Phys. Rev. C92 (1) (2015) 011901. arXiv:1412.4708, doi:10.1103/PhysRevC.92.011901.
- [25] H. Niemi, G. Denicol, How large is the Knudsen number reached in fluid dynamical 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-Mills and 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 and electromagnetic
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: viscous hydrodynamic description of event-by-event central p+p, p+Pb and Pb+Pb collisions at 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 equilibrated quark–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 droplets of 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 , +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 hydrodynamics, arXiv:2010.11957.
- [33] M. P. Heller, M. Spalinski, Hydrodynamics Beyond the Gradient Expansion: 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 and Out of Equilibrium, Cambridge Monographs on Mathematical Physics, Cambridge University Press, 2019. arXiv:1712.05815, doi:10.1017/9781108651998.
- [35] M. Strickland, J. Noronha, G. Denicol, Anisotropic nonequilibrium hydrodynamic 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 relativistic hydrodynamics 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 simulation of 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 symmetries: 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 attractor in 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 Lattice QCD: 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 kinetic theory, 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 for hydrodynamics from weakly coupled pre-equilibrium evolution, JHEP 08 (2016) 171. arXiv:1605.04287, doi:10.1007/JHEP08(2016)171.
- [45] A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting, D. Teaney, Matching the Nonequilibrium Initial Stage of Heavy Ion Collisions to Hydrodynamics 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 kinetic theory, 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 and anisotropic 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 for Gubser 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 dynamics in 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, Thermalization in QCD: theoretical approaches, phenomenological applications, and interdisciplinary connections, arXiv:2005.12299.
- [53] T. Nunes da Silva, D. Chinellato, M. Hippert, W. Serenone, J. Takahashi, G. S. Denicol, M. Luzum, J. Noronha, Pre-hydrodynamic evolution 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 modeling 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 for heavy-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 for heavy-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 dynamics 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 collisions, Phys. Rev. C83 (2011) 034907. arXiv:1007.0130, doi:10.1103/PhysRevC.83.034907.
- [60] M. Martinez, M. Strickland, Dissipative dynamics of highly anisotropic systems, 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)-dimensional anisotropic hydrodynamics, Phys. Rev. C90 (5) (2014) 054910. arXiv:1311.6720, doi:10.1103/PhysRevC.90.054910.
- [63] D. Bazow, U. W. Heinz, M. Martinez, Nonconformal viscous anisotropic hydrodynamics, 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 expansion, 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 dissipative 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 anisotropic fluid 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 viscous pressure 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 ultrarelativistic heavy-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 motion of anisotropic fluid dynamics by a judicious choice of a moment of the Boltzmann equation, Phys. Rev. D94 (12) (2016) 125003. arXiv:1606.09019, doi:10.1103/PhysRevD.94.125003.
- [70] M. Nopoush, R. Ryblewski, M. Strickland, Anisotropic hydrodynamics for 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 quasiparticle anisotropic hydrodynamics for ultrarelativistic heavy-ion collisions, Phys. Rev. Lett. 119 (4) (2017) 042301. arXiv:1703.05808, doi:10.1103/PhysRevLett.119.042301.
- [72] M. Alqahtani, M. Nopoush, R. Ryblewski, M. Strickland, Anisotropic hydrodynamic 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 hydrodynamic modeling 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 dynamical simulations 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 nonlinear 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 and decorrelation of anisotropic flow within the open-computing-language implementation 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 solutions and attractors of higher-order viscous fluid dynamics for Bjorken flow, 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 viscous hydrodynamics 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.
- [81] J. Weil, et al., Particle production and equilibrium properties within a 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, Net baryon diffusion in fluid dynamic simulations of relativistic heavy-ion collisions, 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, Numerical Recipes in C (2nd Ed.): The Art of Scientific Computing, Cambridge University 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 central 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 de Sitter 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 exact solution of the Boltzmann equation, Phys. Rev. D 90 (12) (2014) 125026. arXiv:1408.7048, doi:10.1103/PhysRevD.90.125026.
- [90] W. Florkowski, R. Ryblewski, M. Strickland, Anisotropic hydrodynamics 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 TensorFlow: 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 shock waves in asymptotically AdS5 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 dynamics in 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 systems, 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án, C. F. Sopuerta, M. Triana, M. Zilhão, Holographic collisions in non-conformal theories, JHEP 01 (2017) 026. arXiv:1604.06439, doi:10.1007/JHEP01(2017)026.
- [96] M. P. Heller, A. Kurkela, M. Spaliński, V. Svensson, Hydrodynamization 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 hydrodynamic and 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/f5s1ph73