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 $(\eta_{s}=0)$, along with the transverse slice of a
particlization hypersurface of constant energy density
$\mathcal{E}_{\text{sw}}=0.116$ 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
$\mathcal{P}_{L}{-}\mathcal{P}_{\perp}$ and bulk viscous pressure $\Pi$ 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 $1$
fm$/c$ relative to the one in anisotropic hydrodynamics, allowing it to grow
larger in size.

Figure 19: (Color online) The $\tau{-}x$ slice at $y{\,=\,}\eta_{s}{\,=\,}0$
of the residual shear and mean-field inverse Reynolds numbers in anisotropic
hydrodynamics (VAH, top row) and the shear and bulk inverse Reynolds numbers
in quasiparticle viscous hydrodynamics (VH, middle row) and standard viscous
hydrodynamics (VH2, bottom row), for the fluctuating (3+1)–d T${}_{\text{\sc
R}}$ENTo event. The white contour lines are $\tau{-}x$ slices ($y=\eta_{s}=0$)
of the same particlization hypersurface as Fig. 18. Spacetime regions that are
regulated according to Sec. 3.5 and App. E are circled in black.  Figure 20:
(Color online) The same as Fig. 19 but for the $\tau{-}\eta_{s}$ plane at
$x{\,=\,}y{\,=\,}0$. The $\eta_{s}$–axis in the upper left panel is shifted
transversely to $(x,y)=(7.04,0)$ fm; the corresponding hypersurface slice
(dashed white) is smaller than the one at $x{\,=\,}y{\,=\,}0$ (solid white).

In Figs. 19 and 20 we plot the inverse Reynolds numbers from the three
hydrodynamic simulations in the $\tau{-}x$ plane at $y{\,=\,}\eta_{s}{\,=\,}0$
and in the $\tau{-}\eta_{s}$ plane at $x{\,=\,}y{\,=\,}0$, 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.
$T^{\mu\nu}_{a}=\mathcal{E}u^{\nu}u^{\nu}+\mathcal{P}_{L}z^{\mu}z^{\nu}-\frac{1}{2}\mathcal{P}_{\perp}\Xi^{\mu\nu}$),
its validity is measured by the residual shear inverse Reynolds number [49]

$\displaystyle\text{Re}^{-1}_{\pi W}=$
$\displaystyle\,\sqrt{\frac{\pi_{\perp}{\cdot\,}\pi_{\perp}-2W_{\perp
z}{\,\cdot\,}W_{\perp z}}{\mathcal{P}_{L}^{2}+2\mathcal{P}_{\perp}^{2}}}\,,$
(103a)

as opposed to the shear and bulk inverse Reynolds numbers in second-order
viscous hydrodynamics [20, 5]:

$\text{Re}^{-1}_{\pi}=\,\frac{\sqrt{\pi\cdot\pi}}{\mathcal{P}_{\text{eq}}\sqrt{3}}\,,\qquad\text{Re}^{-1}_{\Pi}=\,\frac{|\Pi|}{\mathcal{P}_{\text{eq}}}\,,$
(104)

with $\pi{\,\cdot\,}\pi=\pi_{\mu\nu}\pi^{\mu\nu}$. One sees that the residual
shear inverse Reynolds numbers stay much smaller than unity inside the
fireball, indicating that a first-order expansion in $O(\text{Re}^{-1}_{\pi
W})$ 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
$\delta B$ to the bulk viscous pressure in anisotropic hydrodynamics:

$\text{Re}_{\delta B}^{-1}=\frac{|\delta B|}{\mathcal{P}_{\text{eq}}}\,.$
(105)

We find that $\delta B$ 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
$\mathcal{E}^{(k)}-2\mathcal{P}_{\perp}^{(k)}-\mathcal{P}_{L}^{(k)}=\mathcal{E}-2\mathcal{P}_{\perp}-\mathcal{P}_{L}-4B$
is at its strongest while the bulk relaxation rate $\tau_{\Pi}^{-1}$ is near a
minimum. After applying the regulation (80) for a brief period, $\delta B$
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 $\Delta\tau\sim 1$ fm/$c$ after the
collision. The bulk inverse Reynolds number in standard viscous hydrodynamics
is also very large near the peak of the bulk viscosity $\zeta/\mathcal{S}$
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 ($\tau<0.6$ fm/$c$) (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 $\tau{-}x$ and $\tau{-}\eta_{s}$ planes. In Fig. 19 one sees that standard
viscous hydrodynamics has the largest hypersurface along the $x$–direction
since its large bulk viscous pressure generates the most viscous heating. It
also has the longest lifetime $\tau_{f}\sim 17$ fm/$c$, for the same reason.
In contrast, it has a pretty narrow waist along the $\eta_{s}$–direction
because its longitudinal flow profile $u^{\eta}$ 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 $\Delta\eta_{s}\sim 0.5$ larger for $\tau\sim 1-2$
fm/$c$) along the $\eta_{s}$–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 ($\sqrt{s_{\text{NN}}}=2.76$ TeV) with fluctuating initial conditions
and different values for the impact parameter $b$ and model parameters. We
also benchmark the OpenMP–accelerated runtime of (3+1)–d non-conformal
hydrodynamics for a central Pb+Pb collision ($b=0$ fm) with the smooth
T${}_{\text{\sc R}}$ENTo initial condition and model parameters from Sec.
4.5.1.

|  Mean time $(s)$ |  Max time $(s)$ |  Mean time per step $(s)$ | 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
Table 1: The mean runtime, maximum runtime, mean runtime per step and mean
number of steps of the (2+1)–d non-conformal hydrodynamic simulations of Pb+Pb
collisions on the fixed grid (auto–grid). (The number of time steps in the
last column does not change for the auto–grid.) These statistics were
generated by running a total of $10000$ simulations ($200$ test parameter
samples times $50$ fluctuating T${}_{\text{\sc R}}$ENTo events) on single-core
Intel Xeon E5-2680 v4 CPUs for each of the three hydrodynamic models.

###  5.1 (2+1)–d non-conformal hydrodynamics on a fixed grid

In this test, we generate $200$ random samples for the impact parameter $b$
and Bayesian model parameters [1, 2]. The impact parameter is sampled from a
piecewise linear distribution,

$P(b)=\Bigg{\\{}\begin{array}[]{ll}b/(2R_{A}^{2})&(0\leq b\leq 2R_{A}),\\\
0&(b{\,>\,}2R_{A}),\end{array}$ (106)

where we set $R_{A}=7$ fm for the lead nuclear radius. The model parameters

$P_{B}=\left[N,p,w,d_{\text{min}},\sigma_{k},T_{\text{sw}},(\eta/\mathcal{S})_{\text{kink}},T_{\eta},a_{\text{low}},a_{\text{high}},(\zeta/\mathcal{S})_{\text{max}},T_{\zeta},w_{\zeta},\lambda_{\zeta}\right]$
(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 T${}_{\text{\sc R}}$ENTo initial condition model parameters
($N$, $p$, $w$, $d_{\text{min}}$, $\sigma_{k}$), are defined in Appendix D
[24] and the viscosity model parameters ($(\eta/\mathcal{S})_{\text{kink}}$,
$T_{\eta}$, $a_{\text{low}}$, $a_{\text{high}}$,
$(\zeta/\mathcal{S})_{\text{max}}$, $T_{\zeta}$, $w_{\zeta}$,
$\lambda_{\zeta}$) are discussed in Sec. 2.6.1. Finally, the switching
temperature $T_{\text{sw}}$ determines the particlization hypersurface of
constant energy density
$\mathcal{E}_{\text{sw}}=\mathcal{E}(T_{\text{sw}})$.,434343We exclude the
ratio between the dimensionless shear relaxation time and shear viscosity
$b_{\pi}=\tau_{\pi}T/(\eta/\mathcal{S})$ as a continuous model parameter [1,
2] and instead allow it to vary between discrete hydrodynamic models. The
rapidity plateau model parameters $\eta_{\text{flat}}$ and $\sigma_{\eta}$
(see Appendix D) are also not considered in this benchmark test.

Figure 21: (Color online) The runtime distribution of the (2+1)–d non-
conformal anisotropic hydrodynamic simulations on the fixed grid (blue) and
auto-grid (red).

For each parameter sample $\\{b,P_{B}\\}$ we run $50$ (2+1)–d non-conformal
hydrodynamic simulations with fluctuating T${}_{\text{\sc R}}$ENTo initial
conditions on a $30$ fm $\times$ $30$ fm transverse grid
($\eta_{s}=0$).444444All other runtime parameters (e.g. the initial time
$\tau_{0}$ and pressure ratio $R$) are the same as the ones used in Sec.
4.5.,454545We found that $93\%$ of the anisotropic hydrodynamic simulations
finish successfully (the second-order viscous hydrodynamic simulations had a
success rate of $99{+}\%$). 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
$\mathcal{P}_{L}^{(k)}$ turned negative. The resulting runtime statistics of
the three hydrodynamic models are shown in Table 1. On average, it takes about
$85s$ to run each (2+1)–d viscous hydrodynamic simulation on the fixed grid,
compared to $200s$ for anisotropic hydrodynamics. The additional routine in
Sec. 3.4 makes the runtime per step in anisotropic hydrodynamics about
$3.5\times$ 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 $b$ increases the participant nucleon multiplicity, which extends
the fireball’s lifetime. The runtime is also very sensitive to the nucleon
width parameter $w$: halving the value for $w$ doubles the spatial resolution
required464646In this work, we set the transverse lattice spacing to $\Delta
x=\Delta y=\frac{1}{5}w$. to capture the spatial variations of the fluctuating
T${}_{\text{\sc R}}$ENTo 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 $30$
fm $\times$ $30$ fm. Although this grid is large enough to fit all possible
fireball sizes produced in Pb+Pb collisions, evolving peripheral collisions
($b\lesssim 2R_{A}$) 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
$\\{b,P_{B}\\}$. With this, we can automatically optimize the grid size to
save computational time and memory.

Figure 22: (Color online) A subset of the scattering matrix used to train the
automated grid for (2+1)–d non-conformal anisotropic hydrodynamic simulations.

To train the regression model, we generate $1000$ training parameter samples
arranged in the column vectors

$\bm{A}=\left[\bm{b},\bm{P}_{B}\right],$ (108)

to construct our feature matrix [91]. For each parameter sample, we run a
single hydrodynamic simulation with a smooth, event-averaged T${}_{\text{\sc
R}}$ENTo initial condition and compute the mean fireball radius $\bar{\bm{r}}$
(again arranged in vector); this will be our target variable474747We define
the fireball radius $r$ as the maximum transverse radius of the particlization
hypersurface in a given hydrodynamic simulation.

$\bm{y}=\bm{\bar{r}}\,.$ (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 $w$ 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
---|---|---|---
$\delta\bar{r}_{\text{RMSE}}$ |  $0.29$ fm |  $0.19$ fm |  $0.22$ fm
$\alpha$ | $0.013$ | $0.009$ | $0.007$
Table 2:  The cross-validated root-mean-square error
$\delta\bar{r}_{\text{RMSE}}$ and regularization parameter $\alpha$ of the
Lasso regression fit for each hydrodynamic model.

Next, we fit a cubic polynomial Lasso regression model with standardization
feature scaling [91] to the training data (108)–(109). The regularization
parameter $\alpha$ is chosen to minimize the root-mean-square error
$\delta\bar{r}_{\text{RMSE}}$, 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 T${}_{\text{\sc R}}$ENTo
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 $\bar{r}_{\text{pred}}$ and set the transverse grid lengths to

$L_{x}=L_{y}=2\left(\bar{r}_{\text{pred}}+\delta\bar{r}_{\text{RMSE}}+\ell\right)\,,$
(110)

where the margin parameter $\ell$ gives the fireball additional room to expand
within the grid. Here we set $\ell=2.5$ fm not only for extra space but also
to account for statistical fluctuations of the fireball radius in fluctuating
T${}_{\text{\sc R}}$ENTo 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 $r$, which can be
characterized with a mean radius and standard deviation
$\bm{Y}=[\bm{\bar{r}},\bm{\sigma_{r}}]$. For more details, we refer the reader
to https://github.com/mjmcnelis/fireball. We find that the automated grid has
about a $99.6\%$ 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 $32-34\,\%$ relative to the fixed grid and provides the simulation with an
additional $1.5-1.6\times$ 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 $\eta_{s}$–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.

Figure 23: (Color online) The simulation runtime of (3+1)–d non-conformal
anisotropic hydrodynamics (red dots), quasiparticle viscous hydrodynamics
(blue dots) and standard viscous hydrodynamics (green dots) on an Intel Xeon
E5-2680 v4 multi-core processor for the smooth T${}_{\text{\sc R}}$ENTo
initial condition from Sec. 4.5.1. The ideal runtime is inversely proportional
to the number of parallel CPU threads (same-colored continuous lines).

###  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 $w$ and rapidity plateau
parameters $\eta_{\text{flat}}$ and $\sigma_{\eta}$ (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
$\eta_{s}=0$. 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 $\tau_{0}{\,\ll\,}\tau_{\text{hydro}}$ 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 $p_{T}$–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 $f_{a}$ (plus residual
shear corrections $\delta\tilde{f}$) 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 $n_{B}$ and $V_{B}^{\mu}$ 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 $u^{\tau}$ are

$\displaystyle\partial_{\tau}u^{\tau}$
$\displaystyle=v^{x}\partial_{\tau}u^{x}+v^{y}\partial_{\tau}u^{y}+\tau^{2}v^{\eta}\partial_{\tau}u^{\eta}+\tau
v^{\eta}u^{\eta}\,,$ (111a) $\displaystyle\partial_{x}u^{\tau}$
$\displaystyle=v^{x}\partial_{x}u^{x}+v^{y}\partial_{x}u^{y}+\tau^{2}v^{\eta}\partial_{x}u^{\eta}\,,$
(111b) $\displaystyle\partial_{y}u^{\tau}$
$\displaystyle=v^{x}\partial_{y}u^{x}+v^{y}\partial_{y}u^{y}+\tau^{2}v^{\eta}\partial_{y}u^{\eta}\,,$
(111c) $\displaystyle\partial_{\eta}u^{\tau}$
$\displaystyle=v^{x}\partial_{\eta}u^{x}+v^{y}\partial_{\eta}u^{y}+\tau^{2}v^{\eta}\partial_{\eta}u^{\eta}\,.$
(111d)

The derivatives of the longitudinal basis vector $z^{\mu}$ are

$\displaystyle\partial_{\tau}z^{\tau}$
$\displaystyle=\frac{\tau}{\sqrt{1{+}u_{\perp}^{2}}}\left(\partial_{\tau}u^{\eta}-\frac{u^{\eta}\left(u^{x}\partial_{\tau}u^{x}{+}u^{y}\partial_{\tau}u^{y}\right)}{1{+}u_{\perp}^{2}}\right)+\frac{z^{\tau}}{\tau}\,,$
(112a) $\displaystyle\partial_{x}z^{\tau}$
$\displaystyle=\frac{\tau}{\sqrt{1{+}u_{\perp}^{2}}}\left(\partial_{x}u^{\eta}-\frac{u^{\eta}\left(u^{x}\partial_{x}u^{x}{+}u^{y}\partial_{x}u^{y}\right)}{1{+}u_{\perp}^{2}}\right)\,,$
(112b) $\displaystyle\partial_{y}z^{\tau}$
$\displaystyle=\frac{\tau}{\sqrt{1{+}u_{\perp}^{2}}}\left(\partial_{y}u^{\eta}-\frac{u^{\eta}\left(u^{x}\partial_{y}u^{x}{+}u^{y}\partial_{y}u^{y}\right)}{1{+}u_{\perp}^{2}}\right)\,,$
(112c) $\displaystyle\partial_{\eta}z^{\tau}$
$\displaystyle=\frac{\tau}{\sqrt{1{+}u_{\perp}^{2}}}\left(\partial_{\eta}u^{\eta}-\frac{u^{\eta}\left(u^{x}\partial_{\eta}u^{x}{+}u^{y}\partial_{\eta}u^{y}\right)}{1{+}u_{\perp}^{2}}\right)\,,$
(112d) $\displaystyle\partial_{\tau}z^{\eta}$
$\displaystyle=\frac{1}{\tau\sqrt{1{+}u_{\perp}^{2}}}\left(\partial_{\tau}u^{\tau}-\frac{u^{\tau}\left(u^{x}\partial_{\tau}u^{x}{+}u^{y}\partial_{\tau}u^{y}\right)}{1{+}u_{\perp}^{2}}\right)-\frac{z^{\eta}}{\tau}\,,$
(112e) $\displaystyle\partial_{x}z^{\eta}$
$\displaystyle=\frac{1}{\tau\sqrt{1{+}u_{\perp}^{2}}}\left(\partial_{x}u^{\tau}-\frac{u^{\tau}\left(u^{x}\partial_{x}u^{x}{+}u^{y}\partial_{x}u^{y}\right)}{1{+}u_{\perp}^{2}}\right)\,,$
(112f) $\displaystyle\partial_{y}z^{\eta}$
$\displaystyle=\frac{1}{\tau\sqrt{1{+}u_{\perp}^{2}}}\left(\partial_{y}u^{\tau}-\frac{u^{\tau}\left(u^{x}\partial_{y}u^{x}{+}u^{y}\partial_{y}u^{y}\right)}{1{+}u_{\perp}^{2}}\right)\,,$
(112g) $\displaystyle\partial_{\eta}z^{\eta}$
$\displaystyle=\frac{1}{\tau\sqrt{1{+}u_{\perp}^{2}}}\left(\partial_{\eta}u^{\tau}-\frac{u^{\tau}\left(u^{x}\partial_{\eta}u^{x}{+}u^{y}\partial_{\eta}u^{y}\right)}{1{+}u_{\perp}^{2}}\right)\,.$
(112h)

The divergence of the spatial fluid velocity $v^{i}=u^{i}/u^{\tau}$ is

$\partial_{i}v^{i}=\frac{\partial_{x}u^{x}-v^{x}\partial_{x}u^{\tau}+\partial_{y}u^{y}-v^{y}\partial_{y}u^{\tau}+\partial_{\eta}u^{\eta}-v^{\eta}\partial_{\eta}u^{\tau}}{u^{\tau}}\,.$
(113)

The longitudinal expansion rate is

$\begin{split}\theta_{L}&=z_{\mu}D_{z}u^{\mu}=-z_{\mu}z^{\nu}\partial_{\nu}u^{\mu}-z_{\mu}z^{\nu}\Gamma^{\mu}_{\nu\lambda}u^{\lambda}\\\
&=-(z^{\tau})^{2}\partial_{\tau}u^{\tau}+z^{\tau}z^{\eta}(\tau^{2}\partial_{\tau}u^{\eta}-\partial_{\eta}u^{\tau})+(\tau
z^{\eta})^{2}\partial_{\eta}u^{\eta}+\tau(z^{\eta})^{2}u^{\tau}\,.\end{split}$
(114)

The transverse expansion rate is

$\theta_{\perp}=\nabla_{\perp\mu}u^{\mu}=\theta-\theta_{L}\,,$ (115)

where

$\begin{split}\theta&=D_{\mu}u^{\mu}=\partial_{\mu}u^{\mu}+\Gamma^{\mu}_{\mu\nu}u^{\nu}\\\
&=\partial_{\tau}u^{\tau}+\partial_{x}u^{x}+\partial_{y}u^{y}+\partial_{\eta}u^{\eta}+\frac{u^{\tau}}{\tau}\end{split}$
(116)

is the scalar expansion rate.

The components of the fluid acceleration

$a^{\mu}=Du^{\mu}=u^{\nu}\partial_{\nu}u^{\mu}+u^{\nu}\Gamma^{\mu}_{\nu\lambda}u^{\lambda}$
(117)

are

$\displaystyle a^{\tau}$
$\displaystyle=u^{\tau}\partial_{\tau}u^{\tau}+u^{x}\partial_{x}u^{\tau}+u^{y}\partial_{y}u^{\tau}+u^{\eta}\partial_{\eta}u^{\tau}+\tau(u^{\eta})^{2}\,,$
(118a) $\displaystyle a^{x}$
$\displaystyle=u^{\tau}\partial_{\tau}u^{x}+u^{x}\partial_{x}u^{x}+u^{y}\partial_{y}u^{x}+u^{\eta}\partial_{\eta}u^{x}\,,$
(118b) $\displaystyle a^{y}$
$\displaystyle=u^{\tau}\partial_{\tau}u^{y}+u^{x}\partial_{x}u^{y}+u^{y}\partial_{y}u^{y}+u^{\eta}\partial_{\eta}u^{y}\,,$
(118c) $\displaystyle a^{\eta}$
$\displaystyle=u^{\tau}\partial_{\tau}u^{\eta}+u^{x}\partial_{x}u^{\eta}+u^{y}\partial_{y}u^{\eta}+u^{\eta}\partial_{\eta}u^{\eta}+\frac{2u^{\tau}u^{\eta}}{\tau}\,.$
(118d)

The components of the longitudinal vector’s comoving time derivative

$\dot{z}^{\mu}=Dz^{\mu}=u^{\nu}\partial_{\nu}z^{\mu}+u^{\nu}\Gamma^{\mu}_{\nu\lambda}z^{\lambda}$
(119)

are

$\displaystyle\dot{z}^{\tau}$
$\displaystyle=u^{\tau}\partial_{\tau}z^{\tau}+u^{x}\partial_{x}z^{\tau}+u^{y}\partial_{y}z^{\tau}+u^{\eta}\partial_{\eta}z^{\tau}+\tau
u^{\eta}z^{\eta}\,,$ (120a) $\displaystyle\dot{z}^{\eta}$
$\displaystyle=u^{\tau}\partial_{\tau}z^{\eta}+u^{x}\partial_{x}z^{\eta}+u^{y}\partial_{y}z^{\eta}+u^{\eta}\partial_{\eta}z^{\eta}+\frac{u^{\tau}z^{\eta}{+}u^{\eta}z^{\tau}}{\tau}\,.$
(120b)

The components of the fluid velocity’s longitudinal derivative

$D_{z}u^{\mu}=-z^{\nu}\partial_{\nu}u^{\mu}-z^{\nu}\Gamma^{\mu}_{\nu\lambda}u^{\lambda}$
(121)

are

$\displaystyle D_{z}u^{\tau}$
$\displaystyle=-z^{\tau}\partial_{\tau}u^{\tau}-z^{\eta}\partial_{\eta}u^{\tau}-\tau
u^{\eta}z^{\eta}\,,$ (122a) $\displaystyle D_{z}u^{x}$
$\displaystyle=-z^{\tau}\partial_{\tau}u^{x}-z^{\eta}\partial_{\eta}u^{x}\,,$
(122b) $\displaystyle D_{z}u^{y}$
$\displaystyle=-z^{\tau}\partial_{\tau}u^{y}-z^{\eta}\partial_{\eta}u^{y}\,,$
(122c) $\displaystyle D_{z}u^{\eta}$
$\displaystyle=-z^{\tau}\partial_{\tau}u^{\eta}-z^{\eta}\partial_{\eta}u^{\eta}-\frac{u^{\tau}z^{\eta}{+}u^{\eta}z^{\tau}}{\tau}\,.$
(122d)

The transverse gradient of the fluid velocity projected along the longitudinal
direction is

$z_{\nu}\nabla_{\perp}^{\mu}u^{\nu}=\Xi^{\mu}_{\alpha}z_{\nu}D^{\alpha}u^{\nu}\,,$
(123)

where the components of

$z_{\nu}D^{\alpha}u^{\nu}=g^{\alpha\beta}z_{\nu}\partial_{\beta}u^{\nu}+g^{\alpha\beta}z_{\nu}\Gamma^{\nu}_{\beta\lambda}u^{\lambda}$
(124)

are

$\displaystyle z_{\nu}D^{\tau}u^{\nu}$
$\displaystyle=z^{\tau}\partial_{\tau}u^{\tau}-\tau^{2}z^{\eta}\partial_{\tau}u^{\eta}-\tau
u^{\eta}z^{\eta}\,,$ (125a) $\displaystyle z_{\nu}D^{x}u^{\nu}$
$\displaystyle=-z^{\tau}\partial_{x}u^{\tau}+\tau^{2}z^{\eta}\partial_{x}u^{\eta}\,,$
(125b) $\displaystyle z_{\nu}D^{y}u^{\nu}$
$\displaystyle=-z^{\tau}\partial_{y}u^{\tau}+\tau^{2}z^{\eta}\partial_{y}u^{\eta}\,,$
(125c) $\displaystyle z_{\nu}D^{\eta}u^{\nu}$
$\displaystyle=-\frac{1}{\tau^{2}}\left(z^{\tau}\partial_{\eta}u^{\tau}-\tau^{2}z^{\eta}\partial_{\eta}u^{\eta}+\tau\left(u^{\eta}z^{\tau}{-}u^{\tau}z^{\eta}\right)\right)$
(125d)

The transverse velocity-shear tensor is525252In the code we obtain
$\sigma_{\perp}^{\mu\nu}$ by applying the projector
$\Xi^{\mu\nu}_{\alpha\beta}$ onto $D^{(\alpha}u^{\beta)}$ directly, rather
than simplifying the expression (126).

$\sigma_{\perp}^{\mu\nu}=\Xi^{\mu\nu}_{\alpha\beta}D^{(\alpha}u^{\beta)}\,,$
(126)

where the components of

$D^{(\alpha}u^{\beta)}=\frac{1}{2}\left(g^{\alpha\rho}\partial_{\rho}u^{\beta}+g^{\beta\rho}\partial_{\rho}u^{\alpha}+g^{\alpha\rho}\Gamma^{\beta}_{\rho\lambda}u^{\lambda}+g^{\beta\rho}\Gamma^{\alpha}_{\rho\lambda}u^{\lambda}\right)$
(127)

are

$\displaystyle D^{(\tau}u^{\tau)}$ $\displaystyle=\partial_{\tau}u^{\tau}\,,$
(128a) $\displaystyle D^{(\tau}u^{x)}$
$\displaystyle=\frac{1}{2}\left(\partial_{\tau}u^{x}-\partial_{x}u^{\tau}\right)\,,$
(128b) $\displaystyle D^{(\tau}u^{y)}$
$\displaystyle=\frac{1}{2}\left(\partial_{\tau}u^{y}-\partial_{y}u^{\tau}\right)\,,$
(128c) $\displaystyle D^{(\tau}u^{\eta)}$
$\displaystyle=\frac{1}{2}\left(\partial_{\tau}u^{\eta}-\frac{\partial_{\eta}u^{\tau}}{\tau^{2}}\right)\,,$
(128d) $\displaystyle D^{(x}u^{x)}$ $\displaystyle=-\partial_{x}u^{x}\,,$
(128e) $\displaystyle D^{(x}u^{y)}$
$\displaystyle=-\frac{1}{2}\left(\partial_{x}u^{y}+\partial_{y}u^{x}\right)\,,$
(128f) $\displaystyle D^{(x}u^{\eta)}$
$\displaystyle=-\frac{1}{2}\left(\partial_{x}u^{\eta}+\frac{\partial_{\eta}u^{x}}{\tau^{2}}\right)\,,$
(128g) $\displaystyle D^{(y}u^{y)}$ $\displaystyle=-\partial_{y}u^{y}\,,$
(128h) $\displaystyle D^{(y}u^{\eta)}$
$\displaystyle=-\frac{1}{2}\left(\partial_{y}u^{\eta}+\frac{\partial_{\eta}u^{y}}{\tau^{2}}\right)\,,$
(128i) $\displaystyle D^{(\eta}u^{\eta)}$
$\displaystyle=-\frac{1}{\tau^{2}}\left(\partial_{\eta}u^{\eta}+\frac{u^{\tau}}{\tau}\right)\,.$
(128j)

The transverse vorticity tensor is

$\begin{split}\omega_{\perp}^{\mu\nu}&=\Xi^{\mu}_{\alpha}\Xi^{\nu}_{\beta}D^{[\alpha}u^{\beta]}\\\
&=D^{[\mu}u^{\nu]}-\frac{u^{\mu}a^{\nu}{-}u^{\nu}a^{\mu}{+}z^{\mu}(D_{z}u^{\nu}{+}z_{\alpha}D^{\nu}u^{\alpha}){-}z^{\nu}(D_{z}u^{\mu}{+}z_{\alpha}D^{\mu}u^{\alpha})}{2}\end{split}$
(129)

where the components of

$D^{[\alpha}u^{\beta]}=\frac{1}{2}\left(g^{\alpha\rho}\partial_{\rho}u^{\beta}-g^{\beta\rho}\partial_{\rho}u^{\alpha}+g^{\alpha\rho}\Gamma^{\beta}_{\rho\lambda}u^{\lambda}-g^{\beta\rho}\Gamma^{\alpha}_{\rho\lambda}u^{\lambda}\right)$
(130)

are

$\displaystyle D^{[\tau}u^{\tau]}$ $\displaystyle=0\,,$ (131a) $\displaystyle
D^{[\tau}u^{x]}$
$\displaystyle=\frac{1}{2}\left(\partial_{\tau}u^{x}+\partial_{x}u^{\tau}\right)\,,$
(131b) $\displaystyle D^{[\tau}u^{y]}$
$\displaystyle=\frac{1}{2}\left(\partial_{\tau}u^{y}+\partial_{y}u^{\tau}\right)\,,$
(131c) $\displaystyle D^{[\tau}u^{\eta]}$
$\displaystyle=\frac{1}{2}\left(\partial_{\tau}u^{\eta}+\frac{\partial_{\eta}u^{\tau}}{\tau^{2}}\right)+\frac{u^{\eta}}{\tau}\,,$
(131d) $\displaystyle D^{[x}u^{x]}$ $\displaystyle=0\,,$ (131e) $\displaystyle
D^{[x}u^{y]}$
$\displaystyle=-\frac{1}{2}\left(\partial_{x}u^{y}-\partial_{y}u^{x}\right)\,,$
(131f) $\displaystyle D^{[x}u^{\eta]}$
$\displaystyle=-\frac{1}{2}\left(\partial_{x}u^{\eta}-\frac{\partial_{\eta}u^{x}}{\tau^{2}}\right)\,,$
(131g) $\displaystyle D^{[y}u^{y]}$ $\displaystyle=0\,,$ (131h) $\displaystyle
D^{[y}u^{\eta]}$
$\displaystyle=-\frac{1}{2}\left(\partial_{y}u^{\eta}+\frac{\partial_{\eta}u^{y}}{\tau^{2}}\right)\,,$
(131i) $\displaystyle D^{[\eta}u^{\eta]}$ $\displaystyle=0\,.$ (131j)

##  Appendix B Geometric source terms

Here we list the components of the geometric source terms $G_{W}^{\mu}$ and
$G_{\pi}^{\mu\nu}$ (Eq. (23)) that appear in the relaxation equations (20) for
$W_{\perp z}^{\mu}$ and $\pi_{\perp}^{{\mu\nu}}$, respectively:

$\displaystyle G_{W}^{\tau}$ $\displaystyle=\tau u^{\eta}W_{\perp
z}^{\eta}\,,$ (132a) $\displaystyle G_{W}^{x}$ $\displaystyle=0\,,$ (132b)
$\displaystyle G_{W}^{y}$ $\displaystyle=0\,,$ (132c) $\displaystyle
G_{W}^{\eta}$ $\displaystyle=\frac{u^{\tau}W_{\perp z}^{\eta}+u^{\eta}W_{\perp
z}^{\tau}}{\tau}\,.$ (132d) $\displaystyle G_{\pi}^{\tau\tau}$
$\displaystyle=2\tau u^{\eta}\pi_{\perp}^{\tau\eta}\,,$ (133a) $\displaystyle
G_{\pi}^{\tau x}$ $\displaystyle=\tau u^{\eta}\pi_{\perp}^{x\eta}\,,$ (133b)
$\displaystyle G_{\pi}^{\tau y}$ $\displaystyle=\tau
u^{\eta}\pi_{\perp}^{y\eta}\,,$ (133c) $\displaystyle G_{\pi}^{\tau\eta}$
$\displaystyle=\tau
u^{\eta}\pi_{\perp}^{\eta\eta}+\frac{u^{\tau}\pi_{\perp}^{\tau\eta}{+}u^{\eta}\pi_{\perp}^{\tau\tau}}{\tau}\,,$
(133d) $\displaystyle G_{\pi}^{xx}$ $\displaystyle=0\,,$ (133e) $\displaystyle
G_{\pi}^{xy}$ $\displaystyle=0\,,$ (133f) $\displaystyle G_{\pi}^{x\eta}$
$\displaystyle=\frac{u^{\tau}\pi_{\perp}^{x\eta}+u^{\eta}\pi_{\perp}^{\tau
x}}{\tau}\,,$ (133g) $\displaystyle G_{\pi}^{yy}$ $\displaystyle=0\,,$ (133h)
$\displaystyle G_{\pi}^{y\eta}$
$\displaystyle=\frac{u^{\tau}\pi_{\perp}^{y\eta}+u^{\eta}\pi_{\perp}^{\tau
y}}{\tau}\,,$ (133i) $\displaystyle G_{\pi}^{\eta\eta}$
$\displaystyle=\frac{2\left(u^{\tau}\pi_{\perp}^{\eta\eta}+u^{\eta}\pi_{\perp}^{\tau\eta}\right)}{\tau}\,.$
(133j)

##  Appendix C Conformal anisotropic transport coefficients

In the conformal limit $m=B=0$, the anisotropic transport coefficients (35) –
(38) only depend the functions $\mathcal{I}_{nrqs}$, which reduce to

$\mathcal{I}_{nrqs}=\frac{g(n{+}s{+}1)!\,\alpha_{L}^{r+1}\Lambda^{n+s+2}\mathcal{R}_{nrq}}{4\pi^{2}(2q)!!}\,,$
(134)

where $g$, $\alpha_{L}$ and $\Lambda$ are given by Eqs. (27), (72) and (73),
respectively. We list the functions $\mathcal{R}_{nrq}$ used in this work
[66]:

$\displaystyle\mathcal{R}_{200}$
$\displaystyle=\alpha_{L}\big{(}1+(1+\xi_{L})t_{L}\big{)}\,,$ (135a)
$\displaystyle\mathcal{R}_{220}$
$\displaystyle=\frac{-1+(1+\xi_{L})t_{L}}{\xi_{L}\alpha_{L}}\,,$ (135b)
$\displaystyle\mathcal{R}_{201}$
$\displaystyle=\frac{1+(\xi_{L}-1)t_{L}}{\xi_{L}\alpha_{L}}\,,$ (135c)
$\displaystyle\mathcal{R}_{240}$
$\displaystyle=\frac{3+2\xi_{L}-3(1+\xi_{L})t_{L})}{\xi_{L}^{2}\alpha_{L}^{3}}\,,$
(135d) $\displaystyle\mathcal{R}_{202}$
$\displaystyle=\frac{3+\xi_{L}+(1+\xi_{L})(\xi_{L}-3)t_{L}}{\xi_{L}^{2}(1+\xi_{L})\alpha_{L}^{3}}\,,$
(135e) $\displaystyle\mathcal{R}_{221}$
$\displaystyle=\frac{-3+(3+\xi_{L})t_{L})}{\xi_{L}^{2}\alpha_{L}^{3}}\,,$
(135f) $\displaystyle\mathcal{R}_{441}$
$\displaystyle=\frac{-15+13\xi_{L}+3(1+\xi_{L})(5+\xi_{L})t_{L}}{4\xi_{L}^{3}\alpha_{L}^{3}}\,,$
(135g) $\displaystyle\mathcal{R}_{402}$
$\displaystyle=\frac{3(\xi_{L}-1)+(\xi_{L}(3\xi_{L}-2)+3)t_{L}}{4\xi_{L}^{2}\alpha_{L}}\,,$
(135h) $\displaystyle\mathcal{R}_{421}$
$\displaystyle=\frac{3+\xi_{L}+(1+\xi_{L})(\xi_{L}-3)t_{L}}{4\xi_{L}^{2}\alpha_{L}}\,,$
(135i) $\displaystyle\mathcal{R}_{422}$
$\displaystyle=\frac{15+\xi_{L}+(\xi_{L}(\xi_{L}-6)-15)t_{L}}{4\xi_{L}^{3}\alpha_{L}^{3}}\,,$
(135j) $\displaystyle\mathcal{R}_{403}$
$\displaystyle=\frac{(\xi_{L}-3)(5+3\xi_{L})+3(1+\xi_{L})(\xi_{L}(\xi_{L}-2)+5)t_{L}}{4\xi_{L}^{3}(1+\xi_{L})\alpha_{L}^{3}}\,,$
(135k)

where $\xi_{L}=\alpha_{L}^{-2}-1$ and
$t_{L}=\mathrm{arctan}\sqrt{\xi_{L}}/\sqrt{\xi_{L}}$.

##  Appendix D T${}_{\text{\sc R}}$ENTo energy deposition model

In the T${}_{\text{\sc R}}$ENTo model, the transverse energy deposition
(GeV/fm2) of a single fluctuating nuclear collision event in the mid-rapidity
region is [24]

$\frac{dE_{T}}{dxdyd\eta_{s}}_{|{\eta_{s}=0}}=N\times T_{R}(x,y)\,,$ (136)

where $N$ is the normalization parameter and

$T_{R}(x,y)=\left(\frac{T^{p}_{A}(x,y)+T^{p}_{B}(x,y)}{2}\right)^{1/p}$ (137)

is the reduced nuclear thickness function, with $p$ being the geometric
parameter. The nuclear thickness function of nucleus $A$ ($B$) is535353In this
work, we only consider Pb+Pb collisions ($A=B=208$) at LHC energies
$\sqrt{s_{\text{NN}}}=2.76$ TeV; the inelastic nucleon–nucleon cross section
is set to $\sigma_{\text{NN}}=6.4$ fm-2.

$T_{A,B}(x,y)=\sum_{n=1}^{N_{\text{part},A,B}}\gamma_{n}\,T_{p}(x-x_{n},y-y_{n})\,,$
(138)

where $N_{\text{part},A,B}$ are the number of participant nucleons from
nucleus $A$ ($B$); the nucleon positions are sampled from a Woods–Saxon
distribution under the constraint that each nucleon–nucleon pair in nucleus
$A$ ($B$) maintains a minimum separation $d_{\text{min}}$ [24]. The
participant nucleon’s thickness function is centered around its sampled
transverse position ($x_{n}$, $y_{n}$):

$T_{p}(x-x_{n},y-y_{n})=\frac{1}{2\pi
w^{2}}\times\exp\left[-\frac{(x{-}x_{n})^{2}+(y{-}y_{n})^{2}}{2w^{2}}\right]\,,$
(139)

where $w$ 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 $\gamma_{n}$ of each participant nucleon is sampled
from the gamma distribution

$P(\gamma)=\frac{k^{k}\gamma^{k-1}\exp\left[-k\gamma\right]}{\Gamma[k]}\,,$
(140)

where $k=\sigma_{k}^{-2}$ and $\sigma_{k}$ is the standard deviation.

For a very brief period $\tau_{0}$ after the collision, we assume the system
is longitudinally free-streaming and static in Milne coordinates (i.e.
$\mathcal{P}_{L}/\mathcal{P}_{\text{eq}}\sim 0$ and $\bm{u}\sim\bm{0}$) so
that $\mathcal{E}(x)\propto 1/\tau_{0}$. Therefore, we initialize the energy
density profile of the hydrodynamic simulation at the starting time $\tau_{0}$
as

$\mathcal{E}(\tau_{0},x,y,\eta_{s})=\frac{1}{\tau_{0}}\times\frac{dE_{T}}{dxdyd\eta_{s}}_{|\eta_{s}=0}\times
f_{L}(\eta_{s})\,.$ (141)

Here we also extend the transverse energy density profile along the spacetime
rapidity direction with a smooth plateau distribution (unitless) [76]:

$f_{L}(\eta_{s})=\exp\left[-\frac{\big{(}|\eta_{s}|-\frac{1}{2}\eta_{\text{flat}}\big{)}^{2}\,\Theta\big{(}|\eta_{s}|-\frac{1}{2}\eta_{\text{flat}}\big{)}}{2\sigma_{\eta}^{2}}\right]\,,$
(142)

where $\eta_{\text{flat}}$ is the plateau length, $\sigma_{\eta}$ is the
standard deviation of the half-Gaussian tails and $\Theta$ is the Heaviside
step function. The initial condition parameter values used in this work are
$N=14.19$ GeV, $p=0.06$, $w=1.11$ fm, $d_{\text{min}}=1.45$ fm,
$\sigma_{k}=1.03$, $\eta_{\text{flat}}=4.0$ and $\sigma_{\eta}=1.8$ [76, 1,
2].

We set the lattice spacings to $\Delta x=\Delta y=\frac{1}{5}w$ and
$\Delta\eta_{s}=\frac{1}{5}\sigma_{\eta}$ to resolve the fluctuating energy
density profile (141) (or event-averaged profile). For the grid size, we set
the longitudinal length to $L_{\eta}=\eta_{\text{flat}}+10\sigma_{\eta}$ to
fit the rapidity plateau (142). The transverse lengths $L_{x}$ and $L_{y}$ 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

$T^{\mu\nu}=\mathcal{E}u^{\mu}u^{\nu}-(\mathcal{P}_{\text{eq}}{+}\Pi)\Delta^{\mu\nu}+\pi^{\mu\nu}\,,$
(143)

where $\pi^{\mu\nu}=\Delta^{\mu\nu}_{\alpha\beta}T^{\alpha\beta}$ is the shear
stress tensor and
$\Pi=-\frac{1}{3}\Delta_{\mu\nu}T^{\mu\nu}-\mathcal{P}_{\text{eq}}$ is the
bulk viscous pressure. We also define the spatial projector
$\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$ and traceless double spatial
projector
$\Delta^{\mu\nu}_{\alpha\beta}=\frac{1}{2}(\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}+\Delta^{\nu}_{\beta}\Delta^{\mu}_{\alpha}-\frac{2}{3}\Delta^{\mu\nu}\Delta_{\alpha\beta})$.
The corresponding dynamical variables

$\bm{q}=(T^{\tau\mu},\pi^{\mu\nu},\Pi)$ (144)

are propagated along with $\mathcal{E}$ and $\bm{u}$ (the mean-field $B$ and
anisotropic variables ($\Lambda$, $\alpha_{\perp}$, $\alpha_{L}$) are not
propagated). Although $\pi^{\mu\nu}$ 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
$\pi^{\tau\eta}$, $\pi^{x\eta}$ and $\pi^{y\eta}$.

###  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):

$\displaystyle\partial_{\tau}T^{\tau\tau}+\partial_{i}(v^{i}T^{\tau\tau})=$
$\displaystyle-\frac{T^{\tau\tau}{+}\tau^{2}T^{\eta\eta}}{\tau}+(\pi^{\tau\tau}{-}\mathcal{P}_{\text{eq}}{-}\Pi)\partial_{i}v^{i}$
(145a)
$\displaystyle+v^{i}\partial_{i}(\pi^{\tau\tau}{-}\mathcal{P}_{\text{eq}}{-}\Pi)-\partial_{i}\pi^{\tau
i}\,,$ $\displaystyle\partial_{\tau}T^{\tau x}+\partial_{i}(v^{i}T^{\tau x})=$
$\displaystyle-\frac{T^{\tau
x}}{\tau}-\partial_{x}(\mathcal{P}_{\text{eq}}{+}\Pi)+\pi^{\tau
x}\partial_{i}v^{i}+v^{i}\partial_{i}\pi^{\tau x}$ (145b)
$\displaystyle-\partial_{i}\pi^{xi}\,,$ $\displaystyle\partial_{\tau}T^{\tau
y}+\partial_{i}(v^{i}T^{\tau y})=$ $\displaystyle-\frac{T^{\tau
y}}{\tau}-\partial_{y}(\mathcal{P}_{\text{eq}}{+}\Pi)+\pi^{\tau
y}\partial_{i}v^{i}+v^{i}\partial_{i}\pi^{\tau y}$ (145c)
$\displaystyle-\partial_{i}\pi^{yi}\,,$
$\displaystyle\partial_{\tau}T^{\tau\eta}+\partial_{i}(v^{i}T^{\tau\eta})=$
$\displaystyle-\frac{3T^{\tau\eta}}{\tau}-\frac{\partial_{\eta}(\mathcal{P}_{\text{eq}}{+}\Pi)}{\tau^{2}}+\pi^{\tau\eta}\partial_{i}v^{i}+v^{i}\partial_{i}\pi^{\tau\eta}$
(145d) $\displaystyle-\partial_{i}\pi^{\eta i}\,,$
$\displaystyle\partial_{\tau}\pi^{\mu\nu}+\partial_{i}(v^{i}\pi^{\mu\nu})=$
$\displaystyle\,\pi^{\mu\nu}\partial_{i}v^{i}+\frac{1}{u^{\tau}}\left[-\frac{\pi^{\mu\nu}}{\tau_{\pi}}+\mathcal{I}_{\pi^{\prime}}^{\mu\nu}-\mathcal{P}_{\pi^{\prime}}^{\mu\nu}-\mathcal{G}_{\pi^{\prime}}^{\mu\nu}\right]\,,$
(145e) $\displaystyle\partial_{\tau}\Pi+\partial_{i}(v^{i}\Pi)=$
$\displaystyle\,\Pi\partial_{i}v^{i}+\frac{1}{u^{\tau}}\left[-\frac{\Pi}{\tau_{\Pi}}+\mathcal{I}_{\Pi}\right]\,,$
(145f)

where
$T^{\eta\eta}=(\mathcal{E}{+}\mathcal{P}_{\text{eq}}{+}\Pi)(u^{\eta})^{2}+(\mathcal{P}_{\text{eq}}{+}\Pi)/\tau^{2}+\pi^{\eta\eta}$,

$\displaystyle\mathcal{I}^{\mu\nu}_{\pi^{\prime}}=$
$\displaystyle\,2\beta_{\pi}\sigma^{\mu\nu}+\Delta^{\mu\nu}_{\alpha\beta}\big{(}2\pi^{\lambda(\alpha}\omega^{\beta)}_{\,\,\,\lambda}-\bar{\tau}_{\pi\pi}\pi^{\lambda(\alpha}\sigma^{\beta)}_{\,\,\,\lambda}\big{)}-\bar{\delta}_{\pi\pi}\pi^{\mu\nu}\theta$
(146a) $\displaystyle+\bar{\lambda}_{\pi\Pi}\Pi\sigma^{\mu\nu}\,,$
$\displaystyle\mathcal{I}_{\Pi}=$
$\displaystyle\,-\beta_{\Pi}\theta-\bar{\delta}_{\Pi\Pi}\Pi\theta+\bar{\lambda}_{\Pi\pi}\pi^{\mu\nu}\sigma_{\mu\nu}\,,$
(146b)

are the gradient source terms for $\pi^{\mu\nu}$ and $\Pi$ and

$\displaystyle\mathcal{P}^{\mu\nu}_{\pi^{\prime}}=$
$\displaystyle\,\left(\pi^{\mu\alpha}u^{\nu}+\pi^{\nu\alpha}u^{\mu}\right)a_{\alpha}\,,$
(147a) $\displaystyle\mathcal{G}^{\mu\nu}_{\pi^{\prime}}=$
$\displaystyle\,u^{\gamma}\Gamma^{\mu}_{\gamma\lambda}\pi^{\nu\lambda}+u^{\gamma}\Gamma^{\nu}_{\gamma\lambda}\pi^{\mu\lambda}\,,$
(147b)

are the spatial projection and geometric source terms for $\pi^{\mu\nu}$. We
also define the velocity-shear tensor
$\sigma^{\mu\nu}=\Delta^{\mu\nu}_{\alpha\beta}D^{(\alpha}u^{\beta)}$ and
vorticity tensor
$\omega^{\mu\nu}=\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}D^{[\alpha}u^{\beta]}$.

The components of $\sigma^{\mu\nu}$ and $\mathcal{G}^{\mu\nu}_{\pi^{\prime}}$
are the same as $\sigma_{\perp}^{\mu\nu}$ and $\mathcal{G}^{\mu\nu}_{\pi}$
after replacing $\Xi^{\mu\nu}_{\alpha\beta}\to\Delta^{\mu\nu}_{\alpha\beta}$
and $\pi_{\perp}^{{\mu\nu}}\to\pi^{\mu\nu}$ in Eqs. (126) and (133),
respectively. The components of $\omega^{\mu\nu}$ are

$\omega^{\mu\nu}=D^{[\mu}u^{\nu]}-\frac{u^{\mu}a^{\nu}{-}u^{\nu}a^{\mu}}{2}\,,$
(148)

where $D^{[\mu}u^{\nu]}$ is given by Eq. (130).

###  E.2 Transport coefficients

In quasiparticle viscous hydrodynamics [79] (i.e. $m=m(T)$ from Fig. 2a), the
relaxation times ($\tau_{\pi}$, $\tau_{\Pi}$) and first-order transport
coefficients ($\beta_{\pi}$, $\beta_{\Pi}$) are given by Eqs. (31) – (32). The
second-order transport coefficients are [66]

$\displaystyle\tau_{\pi\pi}=$
$\displaystyle\,\frac{10+4\bar{c}_{\pi}m^{2}\mathcal{I}_{22}}{7}\,,$ (149a)
$\displaystyle\delta_{\pi\pi}=$
$\displaystyle\,\frac{4}{3}+\bar{c}_{\pi}\mathcal{I}_{22}\Big{(}\frac{m^{2}}{3}-m\frac{dm}{d\mathcal{E}}(\mathcal{E}{+}\mathcal{P}_{\text{eq}})\Big{)}\,,$
(149b) $\displaystyle\lambda_{\pi\Pi}=$
$\displaystyle\,\frac{6}{5}-\frac{2m^{4}}{15}\big{(}\bar{c}_{\mathcal{E}}\mathcal{I}_{00}+\bar{c}_{\Pi}\mathcal{I}_{01}\big{)}\,,$
(149c) $\displaystyle\delta_{\Pi\Pi}=$
$\displaystyle\,1-c_{s}^{2}-\frac{m^{4}}{9}\big{(}\bar{c}_{\mathcal{E}}\mathcal{I}_{00}+\bar{c}_{\Pi}\mathcal{I}_{01}\big{)}$
(149d)
$\displaystyle-m\frac{dm}{d\mathcal{E}}(\mathcal{E}{+}\mathcal{P}_{\text{eq}})\Big{(}\bar{c}_{\mathcal{E}}\mathcal{I}_{21}+\frac{5}{3}\bar{c}_{\Pi}\mathcal{I}_{22}+\frac{3}{m^{2}}\Big{)}\,,$
$\displaystyle\lambda_{\Pi\pi}=$
$\displaystyle\,\frac{1}{3}-c_{s}^{2}+\frac{\bar{c}_{\pi}m^{2}\mathcal{I}_{22}}{3}\,,$
(149e)

where

$\displaystyle\bar{c}_{\pi}=$ $\displaystyle\,\frac{1}{2\mathcal{I}_{42}}\,,$
(150a) $\displaystyle\bar{c}_{\mathcal{E}}=$
$\displaystyle\,-\frac{\mathcal{I}_{41}}{\frac{5}{3}\mathcal{I}_{40}\mathcal{I}_{42}-\mathcal{I}_{41}^{2}}\,,$
(150b) $\displaystyle\bar{c}_{\Pi}=$
$\displaystyle\,\frac{\mathcal{I}_{40}}{\frac{5}{3}\mathcal{I}_{40}\mathcal{I}_{42}-\mathcal{I}_{41}^{2}}\,,$
(150c)

and the function $\mathcal{I}_{nq}$ are given by Eq. (33).

In standard viscous hydrodynamics [67, 20] (i.e. $\bar{m}=m/T\ll 1$), the
relaxation times are given by Eq. (34) and

$\displaystyle\beta_{\pi}\approx$
$\displaystyle\,\frac{\mathcal{E}{+}\mathcal{P}_{\text{eq}}}{5}+O(\bar{m}^{2})\,,$
(151a) $\displaystyle\beta_{\Pi}\approx$
$\displaystyle\,15\Big{(}\frac{1}{3}-c_{s}^{2}\Big{)}^{2}(\mathcal{E}{+}\mathcal{P}_{\text{eq}})+O(\bar{m}^{5})\,.$
(151b)

The second-order transport coefficients (149) reduce to

$\displaystyle\tau_{\pi\pi}\approx$
$\displaystyle\,\frac{10}{7}+O(\bar{m}^{2})\,,$ (152a)
$\displaystyle\delta_{\pi\pi}\approx$
$\displaystyle\,\frac{4}{3}+O(\bar{m}^{2})\,,$ (152b)
$\displaystyle\lambda_{\pi\Pi}\approx$
$\displaystyle\,\frac{6}{5}+O(\bar{m}^{2}\ln\bar{m})\,,$ (152c)
$\displaystyle\delta_{\Pi\Pi}\approx$
$\displaystyle\,\frac{2}{3}+O(\bar{m}^{2}\ln\bar{m})\,,$ (152d)
$\displaystyle\lambda_{\Pi\pi}\approx$
$\displaystyle\,\frac{8}{5}\Big{(}\frac{1}{3}-c_{s}^{2}\Big{)}^{2}+O(\bar{m}^{4})\,.$
(152e)

Both models use the same shear and bulk viscosities as anisotropic
hydrodynamics (e.g. Fig. 3).

###  E.3 Reconstruction formulas for the energy density and fluid velocity

We reconstruct the energy density by solving the following nonlinear equation
via Newton’s method [5, 20]:

$f(\mathcal{E})=0\,,$ (153)

where

$f(\mathcal{E})=\left(\bar{M}^{\tau}{-}\mathcal{E}\right)\left(\bar{M}^{\tau}{+}\mathcal{P}_{\text{eq}}{+}\Pi\right)-(\bar{M}^{x})^{2}-(\bar{M}^{y})^{2}-(\tau\bar{M}^{\eta})^{2}\,,$
(154)

with $\bar{M}^{\mu}=T^{\tau\mu}-\pi^{\tau\mu}$ and
$\mathcal{P}_{\text{eq}}=\mathcal{P}_{\text{eq}}(\mathcal{E})$ being the QCD
equation of state. Using the previous energy density for the initial guess, we
iterate the solution to Eq. (153) as

$\mathcal{E}\leftarrow\mathcal{E}-\frac{f(\mathcal{E})}{df/d\mathcal{E}}\,,$
(155)

where

$\frac{df}{d\mathcal{E}}=c_{s}^{2}(\bar{M}^{\tau}{-}\mathcal{E})-\bar{M}^{\tau}-\mathcal{P}_{\text{eq}}-\Pi\,,$
(156)

with $c_{s}^{2}=d\mathcal{P}_{\text{eq}}/d\mathcal{E}$ being the QCD speed of
sound squared. We repeat the iteration (155) until we achieve sufficient
convergence or the energy density falls below $\mathcal{E}_{\text{min}}$. If
the bulk viscous pressure $\Pi<-\mathcal{P}_{\text{eq}}$, we regulate it so
that $\mathcal{P}_{\text{eq}}+\Pi=0$; this allows us to solve for
$\mathcal{E}$ explicitly [5]:

$\mathcal{E}=\bar{M}^{\tau}-\frac{(\bar{M}^{x})^{2}{+}(\bar{M}^{y})^{2}{+}(\tau\bar{M}^{\eta})^{2}}{{\bar{M}}^{\tau}}\,.$
(157)

Afterwards, we regulate the energy density via Eq. (63) and evaluate the fluid
velocity components as

$\displaystyle u^{x}$
$\displaystyle=\frac{\bar{M}^{x}}{\sqrt{\left(\mathcal{E}{+}\mathcal{P}_{\text{eq}}{+}\Pi\right)\left(\bar{M}^{\tau}{+}\mathcal{P}_{\text{eq}}{+}\Pi\right)}}\,,$
(158a) $\displaystyle u^{y}$
$\displaystyle=\frac{\bar{M}^{y}}{\sqrt{\left(\mathcal{E}{+}\mathcal{P}_{\text{eq}}{+}\Pi\right)\left(\bar{M}^{\tau}{+}\mathcal{P}_{\text{eq}}{+}\Pi\right)}}\,,$
(158b) $\displaystyle u^{\eta}$
$\displaystyle=\frac{\bar{M}^{\eta}}{\sqrt{\left(\mathcal{E}{+}\mathcal{P}_{\text{eq}}{+}\Pi\right)\left(\bar{M}^{\tau}{+}\mathcal{P}_{\text{eq}}{+}\Pi\right)}}\,.$
(158c)

###  E.4 Regulating the shear stress and bulk viscous pressure

In this regulation scheme, we first adjust the shear stress components

$\displaystyle\pi^{\eta\eta}\leftarrow$
$\displaystyle\,\frac{1}{\tau^{2}\left(1{+}u_{\perp}^{2}\right)}\Big{[}\pi^{xx}\big{(}(u^{x})^{2}{-}(u^{\tau})^{2}\big{)}+\pi^{yy}\big{(}(u^{y})^{2}{-}(u^{\tau})^{2}\big{)}$
$\displaystyle+2\big{(}\pi^{xy}u^{x}u^{y}+\tau^{2}(\pi^{x\eta}u^{x}{+}\pi^{y\eta}u^{y})u^{\eta}\big{)}\Big{]}\,,$
(159a) $\displaystyle\pi^{\tau x}\leftarrow$
$\displaystyle\,\frac{\pi^{xx}u^{x}+\pi^{xy}u^{y}+\tau^{2}\pi^{x\eta}u^{\eta}}{u^{\tau}}\,,$
(159b) $\displaystyle\pi^{\tau y}\leftarrow$
$\displaystyle\,\frac{\pi^{xy}u^{x}+\pi^{yy}u^{y}+\tau^{2}\pi^{y\eta}u^{\eta}}{u^{\tau}}\,,$
(159c) $\displaystyle\pi^{\tau\eta}\leftarrow$
$\displaystyle\,\frac{\pi^{x\eta}u^{x}+\pi^{y\eta}u^{y}+\tau^{2}\pi^{\eta\eta}u^{\eta}}{u^{\tau}}\,,$
(159d) $\displaystyle\pi^{\tau\tau}\leftarrow$ $\displaystyle\,\frac{\pi^{\tau
x}u^{x}+\pi^{\tau y}u^{y}+\tau^{2}\pi^{\tau\eta}u^{\eta}}{u^{\tau}}\,,$ (159e)

so that $\pi^{\mu\nu}$ satisfies the orthogonality and tracelessness
conditions

$\displaystyle\pi^{\mu\nu}u_{\nu}$ $\displaystyle=0\,,$ (160a)
$\displaystyle\pi^{\mu}_{\mu}$ $\displaystyle=0\,.$ (160b)

Then we regulate the shear stress and bulk viscous pressure as

$\displaystyle\pi^{\mu\nu}$
$\displaystyle\leftarrow\gamma_{\text{reg}}\pi^{\mu\nu}\,,$ (161a)
$\displaystyle\Pi$ $\displaystyle\leftarrow\gamma_{\text{reg}}\Pi\,,$ (161b)

where

$\gamma_{\text{reg}}=\min\Bigg{(}1,\sqrt{\frac{3\mathcal{P}_{\text{eq}}^{2}}{\pi{\,\cdot\,}\pi+3\Pi^{2}}}\Bigg{)}\,,$
(162)

with $\pi{\,\cdot\,}\pi=\pi_{{\mu\nu}}\pi^{\mu\nu}$. The regulation factor
$\gamma_{\text{reg}}$ usually suppresses $\pi^{\mu\nu}$ and $\Pi$ around the
edges of the fireball at early times $\tau<1$ fm/$c$, 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 $\sqrt{s}=5.02$ TeV, Phys. Lett. B774 (2017) 351–356.  arXiv:1701.07145, doi:10.1016/j.physletb.2017.09.077.
  * [29] P. Romatschke, Do nuclear collisions create a locally 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 $p+p$, $p$+Pb, and Pb+Pb collisions, Phys. Rev. C 102 (5) (2020) 054908.  arXiv:2008.01709, doi:10.1103/PhysRevC.102.054908.
  * [32] C. Plumberg, The multiplicity dependence of pion interferometry in 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