# Comprehensive simulation of heavy-ion collisions at non-zero baryon chemical
potential

A. De<EMAIL_ADDRESS>School of Physics & Astronomy, University of
Minnesota, Minneapolis, MN 55455, USA  J. I. Kapusta<EMAIL_ADDRESS>School
of Physics & Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
M. Singh<EMAIL_ADDRESS>School of Physics & Astronomy, University of
Minnesota, Minneapolis, MN 55455, USA  T. Welle<EMAIL_ADDRESS>School of
Physics & Astronomy, University of Minnesota, Minneapolis, MN 55455, USA

###### Abstract

We present results of hydrodynamic modelling of Au-Au collisions from
$\sqrt{s_{NN}}$ = 7.7 to 200 GeV. Our simulations have three novel components.
Firstly, we use a Linear EXtrapolation of Ultrarelativistic nucleon-nucleon
Scattering to nucleus-nucleus collisions (LEXUS) inspired Monte-Carlo initial
state model. Secondly, we use a crossover equation of state at finite baryon
densities without a critical point. Finally, we use departure functions
derived from the quasiparticle theory of transport coefficients for hadronic
matter at non-zero baryon densities.

††preprint: APS/123-QED

##  I Introduction

The Beam Energy Scan (BES) at the Relativistic Heavy-Ion Collider (RHIC) made
the non-zero baryon chemical potential region of the quantum chromodynamics
(QCD) phase diagram accessible to experiments. BES collides gold nuclei at a
range of beam-energies and consequently scans different regions of the QCD
phase diagram. Experimental programs at the Facility for Antiproton and Ion
Research (FAIR) in Darmstadt and at the Nuclotron-based Ion Collider fAcility
(NICA) in Dubna also collide ions at similar energies. These experimental
programs complement each other and ensure that a broad region of the QCD phase
diagram is covered.

Lattice QCD calculations have shown that the deconfined partons in a cooling
quark-gluon plasma (QGP) change phase to a hadron gas by a smooth crossover
when the net baryon density is zero [1, 2, 3]. This is expected for Pb+Pb
collisions at the Large Hadron Collider (LHC) and for $\sqrt{s_{NN}}=200$ GeV
Au+Au collisions at RHIC. This is in contrast to the first order phase
transition expected between the two phases at high baryon densities [4, 5, 6].
This first-order phase transition line in the QCD phase diagram is expected to
end in a critical point where there will be a second-order phase transition.
One of the main goals of the BES program is to quantify the location and
dynamics of the first and second order phase transitions in QCD. To achieve
this goal, we would need phenomenological modelling of these experiments to
extract relevant physics insights from these experiments.

Significant progress has been made towards phenomenological description of
high-energy heavy ion-collisions where net baryon densities are close to zero.
References [7, 8, 9, 10] provide recent reviews on modelling ultrarelativistic
heavy-ion collisions. The standard model of heavy-ion collisions is comprised
of a pre-hydrodynamic phase, a hydrodynamic phase, and a post hydrodynamic
phase. The pre-hydrodynamic phase is usually modeled by approximating the
colliding nuclei as being infinitely thin in the beam direction [11, 12]. The
is justified because of significant Lorentz contraction. A variety of initial
state models are in popular use for these high energy collisions [13, 14, 15,
16, 17, 18, 19, 20]. The energy density profile from the pre-hydrodynamic
phase provides the initial conditions for the hydrodynamic phase. The
expanding matter is simulated by relativistic viscous fluid dynamics. Apart
from the energy and momentum conservation equations of the hydrodynamics, one
also needs to know the equation of state (EOS) of the nuclear matter.
Typically, in the high energy regime, the EOS is obtained by matching the high
temperature EOS from lattice QCD with low temperature EOS from hadron
resonance gas models using a parameterization [21, 22] or a switching
functions [23]. As the fluid expands, the mean free paths become large and
hydrodynamics is no longer an appropriate theory for describing the system. At
this stage we switch to the kinetic theory description. Hadron distributions
are obtained from the fluid by matching the energy-momentum tensor of the
fluid with the energy-momentum distribution function of hadrons using the
Cooper-Frye prescription [24]. Viscous corrections to the fluid energy-
momentum tensor are matched to the non-equilibrium corrections to the particle
distribution functions using departure functions. The hadrons produced collide
with each other and resonances decay, eventually leading to the chemical
freezeout followed by the kinetic freezeout, at which point all the particles
free stream to detectors.

Extending the standard model of heavy-ion collisions to BES energies with non-
zero baryon densities poses challenges. A recent report by the BEST
collaboration discusses these issues in detail [25]. We will deal with some of
the issues here.

The assumption of ultrarelativistic colliding nuclei being infinitely thin no
longer holds as the collision energies are much lower. Consequently, the
initial state dynamics changes from being approximately two-dimensional to
being fully three-dimensional. This necessitates dynamical initialization [26,
27, 28, 29] of hydrodynamics. For this, we need a 3+1D space-time and momentum
distribution of initial energy and charge sources. The assumption of
longitudinal boost-invariance is no longer a good approximation. There are
some initial state models which deal with these issues [30, 26, 31, 28, 27].

Here we propose a new LEXUS-inspired 3D initial state. LEXUS stands for Linear
EXtrapolation of Ultrarelativistic nucleon-nucleon Scattering to nucleus-
nucleus collisions; the model was introduced in Ref. [32]. LEXUS treats
nucleus-nucleus collisions as superposition of individual nucleon-nucleon
collisions. It uses parameterized data from nucleon-nucleon collisions and
extrapolates it to describe nucleus-nucleus collisions. LEXUS was originally
formulated in momentum space. Here we use a Monte-Carlo sampling of nucleon
positions and treat nuclear collisions as a sequence of binary nucleon-nucleon
collisions to obtain spatial information. Energy loss in each binary collision
is given by a distribution from LEXUS which is fit to nucleon-nucleon
collision data. The idea is to use known results from nucleon-nucleon
collisions to fix the free parameters of the model.

As the baryon charge densities are non-zero, one needs to keep track of an
additional $U(1)$ charge conservation along with the usual energy-momentum
conservation of fluid dynamics. Also, an equation of state is required where
the thermal quantities are also a function of baryon potential in addition to
being functions of temperature. Some progress has been made in this direction
recently. These calculations usually extrapolate lattice calculations to
Taylor expansion coefficients in some order of chemical potential over
temperature [33, 34, 35, 36, 37, 38]. In this work we employ a crossover EOS
[23] which matches parameterized perturbative QCD EOS at high temperatures
without a critical point to hadron resonance gas EOS at low temperatures using
a switching function.

Departure functions need modifications to account for the additional $U(1)$
charge. They have largely been calculated at zero baryon densities [39, 40].
We calculate and implement departure functions calculated using quasiparticle
theory [41] at finite baryon chemical potentials and within the relaxation
time approximation.

This is a comprehensive heavy ion collision simulation at finite baryon
chemical potential meant to be seen as a baseline study. The organization of
this paper is as follows. The initial state model is given in Sec II. We begin
with a short summary of the original LEXUS paper. This helps in understanding
the various model choices we make in the LEXUS inspired 3D initial state. We
discuss the EOS used in Sec. III. We calculate the departure functions at
finite baryon chemical potentials in Sec. IV. Finally, we give our comparisons
with data in Sec. V and our conclusions in Sec. VI.

##  II LEXUS inspired 3D initial state

###  II.1 Summary of the original LEXUS model

We begin by giving a brief summary of the original LEXUS model [32]. Knowledge
of the original LEXUS paper is not required to understand this paper but it
helps to understand the various motivations and inspirations of the LEXUS
based 3D model proposed in this paper.

In the absence of an ab-initio QCD calculation for nuclear collisions, LEXUS
served as a baseline study compared to other phenomenological models of
initial states. When LEXUS was first published, it described the rapidity and
transverse momentum distributions of baryons in central sulfur-sulfur and
lead-lead collisions at the SPS pretty well. LEXUS relies only on the data
from nucleon-nucleon collisions which is extrapolated to form a model of
nucleus-nucleus collisions.

The salient features of LEXUS are as follows

  * •

In LEXUS, nucleons are arranged in rows. When the projectile and the target
rows pass through each other, all the nucleons in one row collide with all the
nucleons in the other row. The main object of interest in LEXUS is the two
particle rapidity distribution $W_{mn}$ for the $m$th projectile nucleon and
$n$th target nucleon immediately after their collision. The $W_{mn}$ is a
result of a collision of a projectile nucleon, which has already undergone
$n-1$ previous collisions and a target nucleon, which has undergone $m-1$
previous collisions. For a projectile nucleon with rapidity $y^{\prime}_{P}$
colliding with a target nucleon of rapidity $y^{\prime}_{T}$ resulting in two
nucleons with rapidities $y_{P}$ and $y_{T}$, we have

$\displaystyle\hskip 20.0ptW_{mn}(y_{P},y_{T})$ $\displaystyle=$
$\displaystyle\int
dy^{\prime}_{P}dy^{\prime}_{T}W^{P}_{mn-1}(y^{\prime}_{P})W^{T}_{m-1n}(y^{\prime}_{T})$
(1) $\displaystyle\times$ $\displaystyle
K(y^{\prime}_{P}+y^{\prime}_{T}\rightarrow y_{P}+y_{T}).$

Here $y$ is the momentum rapidity. The collision kernel $K$ is chosen to be
Markovian. In the original LEXUS model, this kernel was chosen to be

$\displaystyle K(y^{\prime}_{P}+y^{\prime}_{T}\rightarrow y_{P}+y_{T})$
$\displaystyle=$ $\displaystyle\lambda K_{\text{inelastic}}$ (2)
$\displaystyle+$ $\displaystyle(1-\lambda)K_{\text{elastic}},$

where

$\displaystyle
K_{\text{inelastic}}=\frac{\cosh(y_{P}-y^{\prime}_{T})}{\sinh(y^{\prime}_{P}-y^{\prime}_{T})}\frac{\cosh(y^{\prime}_{P}-y_{T})}{\sinh(y^{\prime}_{P}-y^{\prime}_{T})}$
(3)

and

$\displaystyle
K_{\text{elastic}}=\delta(y^{\prime}_{P}-y_{P})\delta(y^{\prime}_{T}-y_{T}).$
(4)

The $\cosh$ function is chosen for the distribution because, in a high-energy
nucleon-nucleon collision, the distribution of outgoing nucleons is flat in
longitudinal momentum (or, in other words, hyperbolic cosine in rapidity). The
coefficient $\lambda$ is the fraction of nucleon-nucleon collisions that are
inelastic and non-diffractive.

  * •

LEXUS arrives at the following expression for the final baryon rapidity
distribution arising from the projectile participants

$\displaystyle\frac{dN_{P}}{dy}(y,\mathbf{b})$ $\displaystyle=$
$\displaystyle\sum_{\bar{m}=1}^{A_{P}}\sum^{\bar{m}}_{m=1}\sum_{n=1}^{A_{T}}W^{P}_{mn}(y)$
(5)
$\displaystyle\times\int\frac{d^{2}s_{P}}{\sigma_{NN}}\mathcal{P}_{n}^{T}(\mathbf{s}_{T})\mathcal{P}_{\bar{m}}^{P}(\mathbf{s}_{P})$

Here $\sigma_{NN}$ is the nucleon-nucleon cross-section which was chosen to be
a constant, and the $s_{P}$ and $s_{T}$ refer to the transverse positions.
There is a symmetrical contribution from the target participants. The
$\mathcal{P}$ was chosen to be a binomial distribution, subject to the
condition that the nucleons follow the Woods-Saxon distribution.

  * •

Transverse momentum is put into LEXUS in the form of a random walk in the
transverse momentum space. The average transverse momentum squared of a baryon
after $k$ collisions is $k\langle p_{T}^{2}\rangle_{NN}$, where $\langle
p_{T}^{2}\rangle_{NN}$ is the average in a nucleon-nucleon collision.

  * •

There was no hydrodynamic evolution used in the original LEXUS model. The
rapidity and transverse momentum distributions of baryons obtained from this
model compared favorably to experimental data available at the time. There
were comparisons made for multiplicity of negatively charged hadrons and their
rapidity and transverse momentum distributions.

###  II.2 LEXUS-inspired 3D initial state

As in LEXUS, we assume that the nucleons follow straight-line trajectories,
striking nucleons from the other nucleus that lie in their path and
interacting with them like two nucleons will interact in free space. The model
does not distinguish between neutrons and protons. The original LEXUS model
was formulated in momentum space but for our Monte Carlo model, we need to
specify the coordinate space information.

We sample the nucleonic positions randomly from a Woods Saxon distribution.

$\rho(r)=\frac{\rho_{0}}{1+\mathrm{exp}[(r-R)/a]}$ (6)

Here $R$ is the nuclear radius, $\rho_{0}$ is the density at the center of the
nucleus, and $a$ is the nuclear skin-thickness. Once the nucleonic positions
are specified, Lorentz contraction is applied in the beam direction. We
initialize all the nucleons with the beam rapidity and no transverse motion.
In the future, this can be easily generalized to include transverse flow. If
the two nuclei are moving with opposite velocities of $v_{z}$ and the radius
of the nucleus is $R$, the overlap time of the two nuclei is

$\tau_{\text{overlap}}=\frac{2R}{\gamma
v_{z}}=\frac{2R}{\sinh(y_{\text{beam}})}.$ (7)

Here $\gamma$ is the Lorentz factor and
$y_{\text{beam}}=\text{arccosh}(\sqrt{s_{NN}}/(2m_{N}))$ is the beam rapidity.
The collision energy per nucleon pair is $\sqrt{s_{NN}}$. The mass of the
proton is $m_{N}=0.938$ GeV. Longitudinal thickness in the overlap leads to a
considerable overlapping time of $\tau_{\text{overlap}}\sim 2-3$ fm for lower
BES energies. The target and the projectile nucleons are initialized with
velocities

$\displaystyle v_{z}^{\text{Projectile}}=\tanh(y_{\text{beam}}),$ (8)
$\displaystyle v_{z}^{\text{Target}}=-\tanh(y_{\text{beam}}).$ (9)

For the target moving to the left and the projectile moving to the right, we
set
$\mathrm{max}\\{z_{i}^{\text{Projectile}}\\}=\mathrm{min}\\{z_{j}^{\text{Target}}\\}=0$
at time $t=0$. In other words, we set the zero of time to be $t=0$ at the
moment when the right-most projectile nucleon crosses the left-most target
nucleon. The longitudinal position of this nucleon pair is defined to be
$z=0$. This particular nucleon pair may or may not collide depending on their
positions in the transverse reaction plane. The nucleons travel in a straight
line and collide with other nucleons within a fixed non-diffractive inelastic
scattering cross section. We use the geometric interpretation of the nucleon-
nucleon cross section and the transverse positions of the nucleons to
determine whether a collision takes place. The total nucleonic scattering
cross section is chosen to be a constant 42 mb, which is approximately the
average of the total nucleon cross section for the collision energies
considered in this paper. We only accept $60\%$ of the collisions because that
is the percentage of total collisions that are inelastic and non-diffractive
[42]. Only the non-diffractive inelastic collisions source into hydrodynamic
evolution.

The model assumes that the nucleons from one nucleus strike nucleons from
other nucleus that lie in their path and interact with them exactly as they
would in free space. The longitudinal coordinate and time of binary collision
is determined by the space-time location of the nucleons crossing each other
with $z=0$ and $t=0$ defined above. The transverse coordinates are defined as

$\displaystyle x_{\text{binary-collision}}=(x_{P}+x_{T})/2,$ (10)
$\displaystyle y_{\text{binary-collision}}=(y_{T}+y_{T})/2.$ (11)

Here $(x_{T},y_{T})$ and $(x_{P},y_{P})$ are the transverse positions of the
participating target and projectile nucleons, respectively.

Participating nucleons continue on their trajectory with reduced momentum from
the point of binary collision. Nucleons can undergo multiple collisions and
can even reverse direction if it they have undergone significant momentum
change from a collision. These binary collision positions can be turned into
space-time rapidity coordinates. Figure 1 shows the binary collision positions
for a $\sqrt{s_{NN}}=200$ GeV and a $\sqrt{s_{NN}}=11.5$ GeV Au+Au collision,
at zero impact parameter. One can see the effect of Lorentz transformation in
the spread in space-time rapidity $\eta=\tanh^{-1}(z/t)$ of the binary
collision locations. There is greater Lorentz contraction at
$\sqrt{s_{NN}}=200$ GeV, consequently binary collision locations occupy a
narrower region in space-time rapidity $\eta$.

Figure 1: Binary collision locations for a Au + Au collision at
$\sqrt{s_{NN}}=200$ GeV (a) and at $\sqrt{s_{NN}}=11.5$ GeV (b).

Let us now describe the rapidity and proper time dependence of the energy and
net baryon density deposition. The energy lost in a single binary collision is
sampled from the probability distribution

$P(y_{\text{loss}})=\frac{\cosh(2y^{\text{total}}_{\text{rest-
frame}}-y_{\text{loss}})}{\sinh(2y^{\text{total}}_{\text{rest-
frame}})-\sinh(y^{\text{total}}_{\text{rest-frame}})}.$ (12)

The absolute value of the incoming nucleon’s rapidity in the pair rest frame
is $y^{\text{total}}_{\text{rest-frame}}$ and the rapidity loss in a given
binary collision is $y_{\text{loss}}$. The probability distribution is chosen
to be a $\cosh$ distribution as explained in the previous section and is
inspired from Eq. (3). The relation between Eq. (3) and the above distribution
is explained in appendix A. The distribution is normalized to 1. This
distribution function is approximately the same distribution function that was
used in the original LEXUS paper. Such distributions have also been used in
Refs. [26, 43, 44]. We sample the rapidity lost $y_{\text{loss}}$ in the range
$[0,y^{\text{total}}_{\text{rest-frame}}]$, where

$y^{\text{total}}_{\text{rest-
frame}}=\left|y_{\text{target}}\right|+\left|y_{\text{projectile}}\right|.$
(13)

|
---|---
|
Figure 2: Binary collision energy deposition in the Cartesian (a) and Milne (b) coordinates. Participant nucleon locations in the Cartesian (c) and Milne (d) coordinates. The plots are for a zero impact parameter $\sqrt{s_{NN}}=200$ GeV Au+Au collision. |
---|---
|
Figure 3: Same as fig. 2 but for a $\sqrt{s_{NN}}=11.5$ GeV Au+Au collision.

The rapidity lost in a single binary collision $y_{\text{loss}}$ determines
the amount of energy sourced into the hydrodynamic equations. The energy that
is deposited in the binary collision is distributed as a Gaussian in rapidity
space. We make this choice because the rapidity distribution of particles
produced in nucleon-nucleon collisions is approximately Gaussian as evidenced
in the Landau model [45]. The width is

$\sigma_{\eta}^{2}=\frac{8}{3}\frac{c_{0}^{2}}{1-c_{0}^{4}}\ln\left(\frac{\sqrt{s_{b}}}{2m_{N}}\right)=\ln\left(\frac{\sqrt{s_{b}}}{2m_{N}}\right)$
(14)

Here $\sqrt{s_{b}}$ signifies the center of mass energy in a particular binary
collision. The speed of sound squared $c_{0}^{2}=1/3$ is a good fit for
nucleon-nucleon collision data [46]. The energy deposition from the collisions
in the hydrodynamic phase occurs after a fixed time
$t_{\text{thermalize}}=0.5$ fm. The energy deposition at
$\\{\tau_{c},\eta_{c},x_{c},y_{c}\\}$ with a width of $\sigma_{\eta}$ is
instantaneous in the $\tau,\eta$ coordinate. The original LEXUS model has no
free parameters. Here, in the LEXUS inspired 3-D initial state, the
thermalization time for deposited energies and participants is the only free
parameter. Importantly, there is no normalization factor associated with the
initial distribution of energy. Once the two-particle $y_{\text{loss}}$ is
sampled, it is equally subtracted from the rapidity of the colliding nucleons
before the collision. That means that even if the target nucleus is moving to
the left and projectile nucleus is moving to the right initially, depending on
the magnitude of $y_{\text{loss}}$ sampled in any given collision, the
projectile might up end up moving to the left or the target to the right.

Once a nucleon has undergone a collision and has deposited some of its energy,
it is called a participant. Participants can undergo further collisions. We
assume that all participants, after they have undergone their last binary
collision, will deposit their energy and baryon charge into the quark-gluon
plasma. The nucleons that have not taken part in any binary collisions, called
spectators, are allowed to pass through. The energy and baryon charge from
participants are deposited after they are propagated for $\Delta\tau=0.5$ fm
after their last collision. This means that the participants will have
propagated to different positions (based on their individual energies) before
being absorbed into the produced matter.

Figures 2 and 3 show the distribution of collision energy deposition and
participant locations in laboratory $t-z$ and Milne $\tau-\eta$ coordinates
for Au+Au collisions at 200 GeV and 11.5 GeV, respectively, at zero impact
parameter. The left-hand plots in both Figures 2 and 3 represent the positions
of the binary collision energy distributions for the two collision energies.
The right-hand plots of the same figures represent the positions of
participants for the two beam energies. One can see the effects of Lorentz
contraction in the spread in the space-time rapidity and $z$ coordinate of the
binary collision energy depositions and participants. For the 200 GeV
collision, there is a wider interval in rapidity around midrapidity, where
binary collision energy deposition is independent of rapidity. This is a check
that at higher energies, the Bjorken model is a more accurate approximation.
This is also visible when one compares the $t-z$ plots for the two energies.
The extent of the deposition regions depicted in the plots is $1\sigma$ of the
aforementioned Gaussian distribution mentioned in Eq. (14). One can see the
fading of the plots at the edges of these deposition regions, which just
signifies the tailing of energy depositions. The $\tau-\eta$ plots start from
$\tau=0.5\,\mathrm{fm}/c$, which is a consequence of the parameter choice of
$\tau_{\text{thermalize}}$ we explained earlier.

The different participant position profiles for the two energies is due to the
difference in energies. For 200 GeV collisions, even after binary collisions,
nucleons possess considerable energy that could be one or two orders of
magnitude greater than their rest mass. Consequently, they still traverse at
speeds close to the speed of light and graze the light cone as is evident in
the figure. On the other hand, for 11.5 GeV collisions, participant nucleons
possess energy of the order of their rest mass and therefore travel at much
lower velocities. This is why they are more uniformly spread in the rapidity
direction.

###  II.3 Initializing Hydrodynamics

We use the publicly available hydrodynamic solver MUSIC [47]. In this section,
we describe how we initialize our hydrodynamical equations.

Each binary collision energy deposition carries the following information in
order to fully quantify the source terms for the hydrodynamic evolution.

$\tau_{c},x_{c},y_{c},\eta_{c},\sigma_{\eta},E,p_{z}$ (15)

Here $E$ is the energy deposited and $p_{z}$ is the momentum of the center of
mass of the binary collision. Participants have the following information.

$\tau,x,y,\eta,E,p_{z}$ (16)

Here $E$ is the remaining energy of the nucleon and $p_{z}$ is the momentum
with which it was traveling at the time it was dissolved into the dense
matter. Each of them has a baryonic charge of $+1$.

These baryon energy and participant depositions act as sources for the
hydrodynamic equations

$\displaystyle\partial_{\nu}T^{\mu\nu}$ $\displaystyle=$ $\displaystyle
S^{\mu}_{\mathrm{source}}(\tau,\mathbf{x}),$ (17)
$\displaystyle\partial_{\mu}J^{\mu}_{B}$ $\displaystyle=$
$\displaystyle\rho_{B,\mathrm{source}}(\tau,\mathbf{x}).$ (18)

The source terms for these equations are provided by our initial state model.

The energy-momentum source $S^{\mu}$ is obtained from all the energy-momentum
depositions as

$S^{\mu}_{\mathrm{source}}=\sum_{i}p_{i}^{\mu}\,f_{\mathrm{smear}}.$ (19)

Here $p^{\mu}_{i}=(E_{i},0,0,p_{z,i})$ is the energy-momentum of binary
collisions events and of the participants at the local space-time point. We
use a Gaussian smearing profile $f_{\mathrm{smear}}$ given by

$f_{\mathrm{smear}}=\frac{1}{N}\exp\left[-\frac{(x_{i}-x_{c})^{2}+(y_{i}-y_{c})^{2}}{\sigma_{\perp}^{2}}-\frac{(\eta_{i}-\eta_{c})^{2}}{\sigma_{\eta}^{2}}\right],$
(20)

with

$N=\Delta\tau(2\pi)^{(3/2)}\sigma_{\perp}^{2}\tau_{i}\sigma_{\eta}.$ (21)

Here $(\tau_{i},x_{i},y_{i},\eta_{i})$ are local the space-time coordinates on
the hydro grid, $\Delta\tau$ is the time-step size on the hydrodynamic solver,
and the transverse smearing width $\sigma_{\perp}$ is chosen to be 0.5 fm. The
energy-momentum source from the participants can be similarly obtained using
Eq. (19), with the space-time positions of the collisions being replaced by
the final positions of the participants. For the participants, the transverse
smearing width $\sigma_{\perp}$ is kept the same as 0.5 fm while the
longitudinal smearing width $\sigma_{\eta}$ is chosen to be 0.2.

The baryon current source $\rho_{B}$ has contributions from every participant
given by

$\displaystyle\rho_{B,\mathrm{source}}$ $\displaystyle=$
$\displaystyle\sum_{i}b_{i}\frac{u_{\mu}p_{i}^{\mu}}{p^{\tau}_{i}}f_{\mathrm{smear}}$
(22) $\displaystyle=$
$\displaystyle\sum_{i}b_{i}\left[u^{0}+\frac{u^{3}p_{z,i}}{E_{i}}\right]f_{\mathrm{smear}}$

The local fluid four-velocity is given by $u^{\mu}$ and $b_{i}$ is the
baryonic charge.

##  III The equation of state

The equation of state used in the hydrodynamic stage of the model is taken
from Ref. [23]. The pressure is expressed as the sum of contributions from a
hadron resonance gas (HRG) model with an excluded volume correction, and a
model obtained from perturbative QCD as

$\displaystyle P(T,\mu_{B})$ $\displaystyle=$
$\displaystyle(1-S(T,\mu_{B}))P_{\mathrm{HRG}}(T,\mu_{B})$ (23)
$\displaystyle+S(T,\mu_{B})P_{\mathrm{pQCD}}(T,\mu_{B})$

The function $S(T,\mu_{B})$ takes values between $0$ and $1$ and parameterizes
the degree to which the hadronized and dehadronized phases contribute. Its
functional form is given by:

$S(T,\mu_{B})=\exp\left[-\left(\frac{T^{2}}{T_{0}^{2}}+\frac{\mu_{B}^{2}}{\mu_{0}^{2}}\right)^{-\frac{r}{2}}\right]$
(24)

The constants $T_{0}=177.12\,\textrm{MeV}$ and $r=5$, which is restricted to
integers, are determined by fitting to lattice values at $T=0$. We take
$\mu_{0}=3\pi T_{0}$. This function is constructed to be strictly increasing
in $T$ and $\mu_{B}$. It is smooth and infinitely differentiable so as not to
introduce discontinuities which would cause a phase transition.

The function $P_{\mathrm{HRG}}(T,\mu_{B})$ is the pressure for an ideal gas of
hadron resonances with an excluded volume from each particle proportional to
its total energy. The proportionality constant $\epsilon_{0}=1.15$ GeV/fm3 is
determined from a fit to lattice values. We include all hadrons composed of u,
d, and s quarks as listed by Ref. [48].

The function $P_{\mathrm{pQCD}}(T,\mu_{B})$ is the standard EOS determined
from perturbative QCD involving 3 flavors of massless quarks [49] with slight
modifications. The first is that the renormalization scale for the running
coupling is given by $M=C_{M}\sqrt{(\pi T)^{2}+(\mu_{B}/3)^{2}}$, with the
constant $C_{M}=3.352$ determined via fit to lattice data. Second, since the
running coupling depends upon $M$ through the quantity
$t=\ln(M^{2}/\Lambda_{\overline{MS}}^{2})$, there will be a divergence at low
temperatures, where the pressure should only depend on the HRG model anyway.
To regulate this divergence, $t$ is replaced with
$t=\ln(C_{S}^{2}+M^{2}/\Lambda_{\overline{MS}}^{2})$, with the constant
$C_{S}=4.28$ again being determined from fitting. The
$\Lambda_{\overline{MS}}$ is taken to be $290\,\textrm{MeV}$.

From $P(T,\mu_{B})$, the baryon number, entropy and energy densities can be
determined from

$s=\left(\frac{\partial P}{\partial T}\right)_{\mu_{B}}\;,\quad
n_{B}=\left(\frac{\partial P}{\partial\mu_{B}}\right)_{T}$ (25)
$\epsilon=Ts+\mu_{B}n_{B}-P$ (26)

For the purposes of our model, $T,\mu_{B},P$ and $s$ are tabulated as
functions of $n_{B}$ and $\epsilon$ through numerical root finding.

##  IV Departure Functions

To first order in a departure from equilibrium the quasiparticle distribution
function for species $a$ is

$f_{a}=f_{a}^{\rm eq}\left(1+\phi_{a}\right)$ (27)

where $f_{a}^{\rm eq}$ is the distribution function in thermal and chemical
equilibrium. In what follows we use a relativistic Boltzmann distribution

$f_{a}^{\rm eq}=\exp[-(E_{a}-\mu_{a})/T]$ (28)

where $\mu_{a}=b_{a}\mu_{B}$ and $b_{a}$ is the baryon number of species $a$.
To match the viscous and thermal conduction contributions to the energy-
momentum tensor and baryon current the $\phi_{a}$ must have the form

$\displaystyle\phi_{a}$ $\displaystyle=$ $\displaystyle-
A_{a}\partial_{\rho}u^{\rho}-B_{a}p_{a}^{\nu}D_{\nu}\left(\frac{\mu_{B}}{T}\right)$
(29) $\displaystyle+$ $\displaystyle
C_{a}p_{a}^{\mu}p_{a}^{\nu}\left(D_{\mu}u_{\nu}+D_{\nu}u_{\mu}+{\textstyle{\frac{2}{3}}}\Delta_{\mu\nu}\partial_{\rho}u^{\rho}\right)$

The functions $A_{a}$, $B_{a}$ and $C_{a}$ only depend on momentum $p$ while
the 4-velocity $u^{\mu}$ only depends on space-time coordinate $x$. Detailed
studies were carried out at zero baryon density [50] and later extended to
nonzero baryon density [41]. One can express the departure from equilibrium as

$\displaystyle f_{a}(E_{a},T,\mu_{B})$ $\displaystyle=$ $\displaystyle
f_{a}^{\rm eq}(E_{a}^{0},T^{0},\mu_{B}^{0})+\delta f_{a}$ (30)
$\displaystyle=$ $\displaystyle f_{a}^{\rm
eq}(E_{a},T^{0},\mu_{B}^{0})+\delta\tilde{f}_{a}$

Here $E_{a}^{0}$ denotes the equilibrium single particle energy and $E_{a}$
the total nonequilibrium energy; it is the latter which is conserved in the
particle collisions. It is the $\delta\tilde{f}_{a}$ which determine the
transport coefficients, not the $\delta f_{a}$. They are related by

$\delta f_{a}=\left[1-\frac{T(\partial E_{a}/\partial
T)_{\sigma}}{E_{a}-\mu_{a}+T\left(\partial\mu_{a}/\partial
T\right)_{\sigma}}\right]\delta\tilde{f}_{a}$ (31)

where $\sigma$ is the entropy per baryon. If there are no mean fields (scalar
or vector) present then $(\partial E_{a}/\partial T)_{\sigma}=0$ and $\delta
f_{a}=\delta\tilde{f}_{a}$. This is the situation we assume here.

It is useful to know the contributions to the pressure, energy density, baryon
density, entropy density, and heat capacity from a single species of particle.

$\displaystyle P_{a}$ $\displaystyle=$ $\displaystyle T\int
d\Gamma_{a}f_{a}^{\rm eq}=\int d\Gamma_{a}\frac{p^{2}}{3E_{a}}f_{a}^{\rm eq}$
$\displaystyle\epsilon_{a}$ $\displaystyle=$ $\displaystyle\int
d\Gamma_{a}E_{a}f_{a}^{\rm eq}$ $\displaystyle n_{Ba}$ $\displaystyle=$
$\displaystyle b_{a}n_{a}=\frac{b_{a}P_{a}}{T}$ $\displaystyle Ts_{a}$
$\displaystyle=$ $\displaystyle\frac{1}{3T}\int d\Gamma_{a}p^{2}f_{a}^{\rm
eq}-\mu_{B}n_{Ba}$

where

$d\Gamma_{a}=(2s_{a}+1)\frac{d^{3}p_{a}}{(2\pi)^{3}}$ (33)

###  IV.1 Shear Viscosity

The shear viscosity is

$\eta=\frac{2}{15}\sum_{a}\int d\Gamma_{a}\frac{p^{4}}{E_{a}}f_{a}^{\rm
eq}C_{a}$ (34)

The usual simplifying assumption when computing the departure distribution
during particlization is that $C_{a}$ is independent of energy and of particle
species [51, 39]. Integration by parts gives

$C_{a}=\frac{\eta/w}{2T^{2}}$ (35)

where enthalpy density is $w=Ts+\mu_{B}n_{B}=\epsilon+P$. This generalizes the
oft-used formula

$C_{a}=\frac{\eta/s}{2T^{3}}$ (36)

to nonzero baryon density.

In the relaxation time approximation

$C_{a}=\frac{\tau_{a}(E_{a})}{2TE_{a}}$ (37)

Numerical results from the linear $\sigma$ model suggest that
$\tau_{a}(E_{a})\propto E_{a}$ [52]. This makes physical sense since higher
momentum or higher mass particles should take longer to reach kinetic
equilibrium. If this is the case, and if one takes
$\tau_{a}(E_{a})=\tau^{\prime}E_{a}$ with $\tau^{\prime}$ a constant, then

$\displaystyle C_{a}$ $\displaystyle=$ $\displaystyle\frac{\tau^{\prime}}{2T}$
$\displaystyle\tau^{\prime}$ $\displaystyle=$ $\displaystyle\frac{\eta}{2Tw}$
(38)

where $\tau^{\prime}$ has units of 1/energy2 and is independent of species.
This may be used to inform the thermal conductivity and bulk viscosity.

###  IV.2 Thermal Conductivity

The formula for the thermal conductivity associated with the baryon current is

$\lambda=\frac{1}{3}\left(\frac{w}{n_{B}T}\right)^{2}\sum_{a}b_{a}\int
d\Gamma_{a}\frac{p^{2}}{E_{a}}f_{a}^{\rm eq}B_{a}\,.$ (39)

Due to energy-momentum conservation, if we have a particular solution
$B_{a}^{\rm par}$ to the integro-differential equation arising from the
Boltzmann equation we can generate another solution as $B_{a}=B_{a}^{\rm
par}-b$, where $b$ is a constant independent of particle species $a$. This
freedom is resolved by the Landau-Lifshitz condition of fit which requires
that $\delta T^{0j}=0$ in the local rest frame. The result is that

$b=\frac{1}{3Tw}\sum_{a}\int d\Gamma_{a}p^{2}f_{a}^{\rm eq}B_{a}^{\rm par}$
(40)

where

$3Tw=\sum_{a}d\Gamma_{a}p^{2}f_{a}^{\rm eq}$ (41)

Substitution into expression (39) gives

$\lambda=\frac{1}{3}\left(\frac{w}{n_{B}T}\right)^{2}\sum_{a}\int
d\Gamma_{a}\frac{p^{2}}{E_{a}}\left(b_{a}-\frac{n_{B}E_{a}}{w}\right)f_{a}^{\rm
eq}B_{a}^{\rm par}$ (42)

The simplest approximation is to take $B_{a}^{\rm par}$ = constant, as in the
case of the shear viscosity. However, it is easily shown that this results in
$\lambda=0$. Alternatively, with this approximation one gets $b=B_{a}^{\rm
par}$, and so $B_{a}=0$ and again $\lambda=0$.

In the relaxation time approximation

$B_{a}^{\rm
par}=\frac{\tau_{a}(E_{a})}{E_{a}}\left(b_{a}-\frac{n_{B}}{w}E_{a}\right)$
(43)

To be consistent with expression (35) one uses
$\tau_{a}(E_{a})=\tau^{\prime}E_{a}$. Then

$b=\frac{\tau^{\prime}}{3Tw}\left[\sum_{a}b_{a}\int d\Gamma_{a}p^{2}f_{a}^{\rm
eq}-\frac{n_{B}}{w}\sum_{a}\int d\Gamma_{a}p^{2}E_{a}f_{a}^{\rm eq}\right]$
(44)

After some manipulation of integrals and using thermodynamic identities this
can be written as

$b=\frac{\tau^{\prime}T}{w^{2}}\left[s\sum_{a}b_{a}w_{a}-n_{B}w-n_{B}T\frac{\partial
w}{\partial T}\right]$ (45)

and then as

$\displaystyle b$ $\displaystyle=$
$\displaystyle\frac{\tau^{\prime}T}{w^{2}}\Big{[}Ts(T\chi_{T\mu}+\mu_{B}\chi_{\mu\mu})$
(46) $\displaystyle-$ $\displaystyle
Tn_{B}(T\chi_{TT}+\mu_{B}\chi_{T\mu})-n_{B}w\Big{]}$

Here the susceptibilities are

$\chi_{xy}=\frac{\partial^{2}P(T,\mu)}{\partial x\partial y}$ (47)

Finally

$B_{a}=\tau^{\prime}\left(b_{a}-\frac{n_{B}}{w}E_{a}\right)-b(T,\mu_{B})$ (48)

Unfortunately the expression for $B_{a}$ is not a simple pocket formula as it
is for $C_{a}$.

Using Boltzmann statistics, the integrals we need are

$\displaystyle\sum_{a}b_{a}^{2}\int d\Gamma_{a}\frac{p^{2}}{E_{a}}f_{a}^{\rm
eq}=3T^{2}\chi_{\mu\mu}$ $\displaystyle\sum_{a}b_{a}\int
d\Gamma_{a}p^{2}f_{a}^{\rm
eq}=3T^{2}(T\chi_{T\mu}+\mu_{B}\chi_{\mu\mu}+n_{B})$
$\displaystyle\sum_{a}\int d\Gamma_{a}p^{2}E_{a}f_{a}^{\rm eq}=6T^{2}w$
$\displaystyle+3T^{2}(T^{2}\chi_{TT}+2T\mu_{B}\chi_{T\mu}+\mu_{B}^{2}\chi_{\mu\mu})$
(49)

which results in

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\tau^{\prime}T^{2}\left[\left(\frac{s}{n_{B}}\right)^{2}\chi_{\mu\mu}-2\left(\frac{s}{n_{B}}\right)\chi_{T\mu}+\chi_{TT}\right]$
(50) $\displaystyle=$
$\displaystyle\tau^{\prime}\left(\frac{T}{n_{B}}\right)^{2}w\det\chi
v_{\sigma}^{2}=\frac{\eta T}{2}\frac{\det\chi}{n_{B}^{2}}v_{\sigma}^{2}$

where $v_{\sigma}$ is the speed of sound (see below) and

$\det\chi=\chi_{TT}\chi_{\mu\mu}-\chi^{2}_{T\mu}$ (51)

###  IV.3 Bulk Viscosity

The bulk viscosity is

$\zeta=\frac{1}{3}\sum_{a}\int d\Gamma_{a}\frac{p^{2}}{E_{a}}f_{a}^{\rm
eq}A_{a}\,.$ (52)

In the relaxation time approximation the particular solution is

$A_{a}^{\rm
par}=\frac{\tau_{a}}{3T}\left\\{\frac{p^{2}}{E_{a}}-3\left[v_{n}^{2}E_{a}+(v_{s}^{2}-v_{n}^{2})\mu_{a}\right]\right\\}$
(53)

where, for brevity of notation, we have defined

$\displaystyle v_{n}^{2}$ $\displaystyle=$ $\displaystyle\left(\frac{\partial
P}{\partial\epsilon}\right)_{n}=\frac{s\chi_{\mu\mu}-n_{B}\chi_{\mu
T}}{T\det\chi}$ $\displaystyle v_{s}^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial
P}{\partial\epsilon}\right)_{s}=\frac{n_{B}\chi_{TT}-s\chi_{\mu
T}}{\mu_{B}\det\chi}$ $\displaystyle v_{\sigma}^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial
P}{\partial\epsilon}\right)_{\sigma}=\frac{v_{n}^{2}Ts+v_{s}^{2}\mu_{B}n_{B}}{w}$
(54)

Of course waves do not physically propagate at constant $n$ or $s$, only at
constant $\sigma$.

In order to satisfy the Landau-Lifshitz condition of fit we must allow for the
functional form

$A_{a}=A_{a}^{\rm par}-a_{E}E_{a}-a_{B}b_{a}$ (55)

Here

$\displaystyle a_{E}$ $\displaystyle=$
$\displaystyle\frac{X_{B}Z_{B}-Y_{B}Z_{E}}{Y_{E}X_{B}-X_{E}Y_{B}}$
$\displaystyle a_{B}$ $\displaystyle=$
$\displaystyle\frac{Y_{E}Z_{E}-X_{E}Z_{B}}{Y_{E}X_{B}-X_{E}Y_{B}}$ (56)

where

$\displaystyle X_{E}$ $\displaystyle=$ $\displaystyle
T(T^{2}\chi_{TT}+2\mu_{B}T\chi_{T\mu}+\mu_{B}^{2}\chi_{\mu\mu})$
$\displaystyle X_{B}$ $\displaystyle=$ $\displaystyle
T(T\chi_{T\mu}+\mu_{B}\chi_{\mu\mu})$ $\displaystyle Y_{E}$ $\displaystyle=$
$\displaystyle T(T\chi_{T\mu}+\mu_{B}\chi_{\mu\mu})$ $\displaystyle Y_{B}$
$\displaystyle=$ $\displaystyle T\chi_{\mu\mu}$ (57)

and

$\displaystyle Z_{E}$ $\displaystyle=$ $\displaystyle\sum_{a}\int
d\Gamma_{a}E_{a}A_{a}^{\rm par}f_{a}^{\rm eq}$ $\displaystyle Z_{B}$
$\displaystyle=$ $\displaystyle\sum_{a}b_{a}\int d\Gamma_{a}A_{a}^{\rm
par}f_{a}^{\rm eq}$ (58)

As in the case of thermal conductivity, one should take
$\tau_{a}=\tau^{\prime}E_{a}$ to be consistent with the standard result for
the shear viscosity. This leads to complicated formulas for $Z_{E}$ and
$Z_{B}$. To evaluate the lengthy integrals it is useful to have the identities

$\displaystyle T\frac{\partial}{\partial\mu_{B}}f_{a}^{\rm eq}$
$\displaystyle=$ $\displaystyle b_{a}f_{a}^{\rm eq}$ $\displaystyle
T^{2}\frac{\partial}{\partial T}f_{a}^{\rm eq}$ $\displaystyle=$
$\displaystyle(E_{a}-\mu_{a})f_{a}^{\rm eq}$ $\displaystyle
T\left(T\frac{\partial}{\partial
T}+\mu_{B}\frac{\partial}{\partial\mu_{B}}\right)f_{a}^{\rm eq}$
$\displaystyle=$ $\displaystyle E_{a}f_{a}^{\rm eq}$ (59)

which are valid for Boltzmann statistics. The results are

$\displaystyle
Z_{B}/T\tau^{\prime}=n_{B}+(1-v_{n}^{2}-v_{s}^{2})\mu_{B}\chi_{\mu\mu}+(1-2v_{n}^{2})T\chi_{T\mu}$
$\displaystyle-(v_{n}^{2}+v_{s}^{2})T\mu_{B}\chi_{T\mu\mu}-v_{n}^{2}T^{2}\chi_{TT\mu}-v_{s}^{2}\mu_{B}^{2}\chi_{\mu\mu\mu}$
(60)

and

$\displaystyle Z_{E}/T\tau^{\prime}=2w+(1-3v_{n}^{2})T^{2}\chi_{TT}$
$\displaystyle+2(1-2v_{n}^{2}-v_{s}^{2})T\mu_{B}\chi_{T\mu}+(1-v_{n}^{2}-v_{s}^{2})\mu_{B}^{2}\chi_{\mu\mu}$
$\displaystyle-
v_{n}^{2}T^{3}\chi_{TTT}-(2v_{n}^{2}+v_{s}^{2})T^{2}\mu_{B}\chi_{TT\mu}$
$\displaystyle-(v_{n}^{2}+2v_{s}^{2})T\mu_{B}^{2}\chi_{T\mu\mu}-v_{s}^{2}\mu_{B}^{3}\chi_{\mu\mu\mu}$
(61)

The bulk viscosity can be expressed in terms of the equation of state and
susceptibilities, but the expression is very long and not displayed here.

###  IV.4 Relaxation time for baryon current

Now the question is how to choose the relaxation time for the baryon current
$\tau_{B}$. This is the time constant that appears in the Cattaneo equation
[53, 54]. Properly speaking the baryon diffusion constant is $D_{B}$ and has
units of length (in natural units with $c=1$).

$D_{B}=\frac{\lambda
T}{\chi_{\mu\mu}}\left(\frac{n_{B}}{w}\right)^{2}=\frac{\kappa_{B}}{T\chi_{\mu\mu}}$
(62)

Collective baryon fluctuations travel with speed $v_{B}^{2}=D_{B}/\tau_{B}$
[53, 54]. This means that $\tau_{B}>D_{B}$ in order that signals not travel
faster than the speed of light. Using the above formulas for $\lambda$ and
$\tau^{\prime}$ we have

$D_{B}=\frac{\eta}{2}\left(\frac{T}{w}\right)^{2}\frac{\det\chi}{\chi_{\mu\mu}}v_{\sigma}^{2}$
(63)

This expression is simplified if all the particles are assumed to have a
constant relaxation time [55]. The two expressions are compared in appendix B.

The two speeds $v_{B}$ and $v_{\sigma}$ cannot be too different, so we choose
them equal to obtain

$\tau_{B}=\frac{\eta}{2}\left(\frac{T}{w}\right)^{2}\frac{\det\chi}{\chi_{\mu\mu}}$
(64)

Figure 4: The $v_{2}$ for $\pi^{+}$ for the centrality bin $0-10\%$. Shaded
areas represent statistical uncertainties. Experimental data are from [56,
57]. Figure 5: The $v_{2}$ for proton for the centrality bin $0-10\%$. Shaded
areas represent statistical uncertainties. Experimental data are from [56,
57].

##  V Results

Our goal is to conduct realistic heavy-ion collision simulations at finite
baryon chemical potential with a LEXUS-inspired initial-state model, a
crossover EOS without a critical point, and departure functions derived using
quasiparticle theory at non-zero baryon chemical potential. In many regards,
it is a straightforward modeling of the heavy ion simulations to test our
physics understanding of the initial state at finite baryon chemical potential
and incorporating departure functions that are specifically calculated for
non-zero baryon chemical potential. In this section we compare our simulations
with experimental data from RHIC. We calculate the single particle spectra and
the flow harmonics.

Our LEXUS-inspired model provides the initial conditions for MUSIC. The
produced matter is then hydrodynamically evolved until freezeout which is
assumed to occur on the constant energy density surface 0.205 GeV/fm3. We use
the baryon diffusion constant and baryon relaxation times given in equations
(63) and (64). The particle spectra are computed using the Cooper-Frye
procedure and our new departure functions. The Cooper-Frye procedure is
performed using mode 3 in MUSIC and the obtained particle distributions
undergo resonance decays using mode 4 in MUSIC. All hadronic resonances with
mass less than 1.8 GeV in the PDG table [48] were included in our study. No
hadronic after-burner was employed in this work. Effects of an hadronic
afterburner is left for future investigation.

Impact parameters are randomly sampled between $0$ and $20$ fm. The centrality
is then determined by categorizing the $5\%$ of events with the highest total
energy as the $0-5\%$ central events, the next $5\%$ events in terms of total
energy as the $5-10\%$ central events, and so on. This process is close to
what is done in experiments, which bin events on charged particle
multiplicity. Thus, the total energy in the initial state is used as a proxy
for the final multiplicity.

###  V.1 Hydrodynamic flow

One of the key signatures of flow in heavy ion collisions is the second
harmonic coefficient of particle multiplicity $v_{2}$.

$\displaystyle\frac{dN}{p_{T}dp_{T}dyd\phi}=\frac{dN}{2\pi
p_{T}dp_{T}dy}\left[1+\sum_{n}2v_{n}(y,p_{T})\right.$
$\displaystyle\times\left.\;\cos(n\phi-n\Psi_{n}(p_{T}))\vphantom{\sum_{n}2v_{n}(y,p_{T})}\right]$
(65)

The event-plane angles $\Psi_{n}(p_{T})$ are determined event by event and
given as

$\Psi_{n}(p_{T})=\left(\tan^{-1}\frac{\sum_{i}\sin(n\phi_{i})}{\sum_{i}\cos(n\phi_{i})}\right)/n.$
(66)

We compute the $v_{2}$ flow harmonic as a function of $p_{T}$. The shear
viscosity to entropy ratio $\eta/s$ that we use is a constant for a given
collision energy and not T and $\mu_{B}$ dependent. The ratio $\eta/s$ was
adjusted to match the $p_{T}$-differential $v_{2}$ data for various collision
energies in the $0-10\%$ centrality class as a function of $p_{T}$, when
departure functions are included. Calculations without departure functions are
done with same value of $\eta/s$. Table 1 lists the $\eta/s$ for various
collision energies.

Collision energy | $\eta/s$
---|---
200 GeV | 0.08
62.4 GeV | 0.08
14.5 GeV | 0.08
11.5 GeV | 0.08
7.7 GeV | 0.14
Table 1: $\eta/s$ for various collision energies.

Figure 4 shows the $p_{T}$-differential $v_{2}$ for $\pi^{+}$ at midrapidity
for the above mentioned energies compared to STAR data. Theoretical
calculations are shown with the statistical error band. There is good
agreement between the simulation and the experimental data. The agreement with
STAR data for protons for the $0-10\%$ centrality class is also reasonable as
shown in Fig. 5.

One can see that the agreement with the data ceases to be very good beyond
$p_{T}$ of about 1.5 GeV/c. That is to be expected because hydrodynamics is a
long-wavelength theory and is not particularly effective for large $p_{T}$.

The $v_{2}$ is mainly driven by the initial geometry of the system. Higher
harmonics are more sensitive to small-scale initial fluctuations. The $v_{3}$
and $v_{4}$ measurements for 200 GeV collisions are available and are compared
to the model predictions in Figs. 6 and 7, respectively. Our results are in
good agreement with the experimental data for $v_{3}$ and $v_{4}$ as well.

Departure functions have a sizeable effect on $p_{T}$ differential $v_{n}$s.
This effect is more pronounced as we go to lower collision energies where
baryon densities are higher. Inclusion of departure functions requires lower
values of $\eta/s$ to explain the data.

Figure 6: The $v_{3}$ for $\sqrt{s_{NN}}=200$ GeV Au+Au collisions for
centrality $0-10\%$. Shaded areas represent statistical uncertainties.
Experimental data are from from [58]. Figure 7: The $v_{4}$ for
$\sqrt{s_{NN}}=200$ GeV Au+Au collisions for centrality $0-10\%$. Shaded areas
represent statistical uncertainties. Experimental data are from [58].

We present our predictions for the $\pi^{+}$ $v_{3}$ and $v_{4}$ for BES
energies in Fig. 8. These should be compared to BES measurements when the data
becomes available [59].

Figure 8: Predictions for $v_{3}$ and $v_{4}$ for $\pi^{+}$ for the centrality
bin $0-10\%$. Shaded areas represent statistical uncertainties. Figure 9: Pion
transverse momentum distributions for Au+Au collision for the centrality bin
$0-5\%$. Experimental data are from [60, 61, 62]. Figure 10: Proton transverse
momentum distributions for Au+Au for the centrality bin $0-5\%$. Experimental
data are from [60, 61, 62].

###  V.2 Effect of the departure functions on hadron rapidity and transverse
momentum distributions

We are working under the assumption of an overlapping period of applicability
between hydrodynamic and kinetic theories which lets us match the energy-
momentum tensor and baryon current from fluid dynamics to particle
distributions in the kinetic theory. We calculate $p_{T}$ spectra for pions
and protons with and without the use of departure functions. The results are
compared to experimental data. These comparisons are shown in Figs.9 and 10
for five collision energies. There is an enhancement when we include the
$\delta f$ terms, which is to be expected, since there is an increase in the
hadrons being sampled because of the out-of-equilibrium corrections.
Corrections due to the $\delta f$ terms are more prominent for the proton
yield than for the pion yield. That makes sense because the out-of-equilibrium
corrections depend on the relaxation times. We have taken the relaxation times
to be proportional to the mass of the hadronic species; hence, the
contribution of the $\delta f$ terms to the particle multiplicities are
greater in the proton yields compared to the pion yields.

In general, we expect our LEXUS based initial state model to be more accurate
for low-energy heavy-ion collisions compared to high-energy heavy-ion
collisions. The reason is that we are treating nucleus-nucleus collisions as a
sequence of nucleon-nucleon collisions. One would surmise that as one goes to
higher and higher collision energies, there will be additional physics from
partonic degrees of freedom that are missing from a LEXUS based model. Our aim
is to give a 3D model for lower-energy collisions, so we should not be too
concerned about discrepancies at higher-energy collisions.

There are various other factors that might contribute to the discrepancies
with the data. We are not using an hadronic afterburner in this work; this is
because we want to focus on the effects of the departure functions $\delta f$.
Using an hadronic afterburner will require particle sampling and will increase
the computational cost to achieve similar statistics. The hadronic
rescatterings among light mesons and baryons largely blue-shift the
distributions and shift them to higher $p_{T}$. This is the pion wind effect
that pushes the heavier particles to the high $p_{T}$ region [63]. The effect
should be more pronounced in the case of protons than pions. If one compares
Figs. 9 and 10, the difference with data is greater for protons. Hence, the
addition of hadronic rescatterings will have the desired effect. The net
proton rapidity distribution is expected to be widened by scatterings with
other hadrons.

We have not included contributions from weak decays of baryons which STAR has
included. This will also enhance the proton multiplicity in our results.

We do not consider the possibility of transverse flow in the initial state.
The presence of collective flow in the initial state might necessitate
increasing the $\eta/s$ at lower energies. This would lead to more entropy
production, which will enhance the multiplicities and bring them closer to the
experimental data.

We are not considering bulk viscosity for the sake of simplicity. The presence
of bulk viscosity would cause an increase in the multiplicities of all the
hadronic species. The inclusion of bulk viscosity will involve more tuning of
this work and is deferred to future investigations.

This paper includes the effects of baryon number but not electric charge. This
calculation assumes that the ratio of electric charge to baryon number is one-
half. Consequently, if it were not for the small difference in proton and
neutron masses, the Cooper-Frye procedure would produce the same exact
distribution for the two particles. As this mass difference is very small, the
proton to neutron ratio is almost 1 before the resonance decays. In reality,
before collisions, gold nuclei have more neutrons than protons. During the
collision, there is a net conversion of neutrons to protons accompanied by
more negatively charged pions than positively charged ones, but the proton to
neutron ratio will be less than one. So, conserving both baryon number and
electric charge requires that the proton to neutron ratio will be somewhere
between $Z/A$ and 1, though in this calculation it is almost 1. This was
studied long ago [64]. Doing this better will require the use of chemical
potentials for electric charge and strangeness in addition to baryon number,
which is left for a future investigation.

We want to emphasize that, although there is some physics missing in our
initial state, our initial-state model is an extrapolation of nucleon-nucleon
collision data. There is no additional normalization factor that has been
introduced to match the predicted multiplicities of hadrons with the
experimental data. There is a common practice of introducing normalization
factors in initial-state models to match experimental data which represent
unknown physics in the initial-state model, which is not necessary in this
model. For better or worse,it is absolutely normalized.

Let us now turn our attention to the rapidity distributions. The invariant
yields of pions and protons as functions of pseudorapidity are plotted in
Figs. 11 and 12. Comparisons have been made to show how the distributions
differ with and without the inclusion of departure functions. How they vary
with pseudorapidity is more relevant for low-energy collisions where the
Bjorken model is no longer a good approximation. Detailed experimental data is
not yet available. If one compares the $\pi^{+}$ pseudorapidity distribution
at 7.7 GeV and at 200 GeV, it is flatter at higher energies, which is a
consequence of the initial binary collision energy depositions being flatter
for higher energies. This is another affirmation that the Bjorken model is a
good model at higher energies.

The influence of departure functions is more pronounced for proton
distributions, as is evident in Fig. 12. The consequences of the initial
space-time conditions for participants are shown in Figs. 2 and 3. The
separation of participant space-time deposition at higher collision energies
leads to separate baryon pseudorapidity peaks present at higher energies,
whereas a more uniform space-time deposition of participants at lower
collision energies leads to the broader pseudorapidity distribution of protons
at lower energies.

Figure 11: Pion pseudorapidity distributions for Au+Au for the centrality bin
$0-5\%$. Shaded areas represent statistical uncertainties. Figure 12: Proton
pseudorapidity distributions for Au+Au for the centrality bin $0-5\%$. Shaded
areas represent statistical uncertainties.

##  VI Conclusions

The BES program at RHIC has conducted Au+Au heavy ion collisions at center-of-
momentum collision energies from 7.7 GeV to 200 GeV with the estimated baryon
chemical potential ranging from 20 MeV to 420 MeV [65]. In order to understand
the experimental results and in order to ‘detect’ the presence of the QCD
critical point, we need high precision simulations at lower collision
energies. This work is a baseline simulation to such studies.

We used an initial state inspired by LEXUS. The input is the measured binary
nucleon collision cross-section. We employ dynamical initialization of the
hydrodynamic solver. The hydrodynamic solver used in this paper is MUSIC [47].
We use departure functions calculated at finite baryon chemical potential
within the relaxation time approximation. The EOS used is a crossover equation
of state without a critical point.

We compared the transverse momentum dependent flow coefficients $v_{2}$ and
the single particle transverse momentum distributions with the STAR data.
Although we find reasonable agreement with $v_{2}$ for five collision energies
between 7.7 GeV and 200 GeV, our model underestimates the transverse momentum
distributions with respect to experimental data. We discuss the possible
sources of these discrepancies. We believe that including an hadronic
afterburner, and accounting for a non-zero bulk viscosity, considering initial
flow in the initial state , the hadronic multiplicities predicted with our
framework will have a reasonable agreement with experimental data.

Obviously, a future direction is to incorporate hadronic scatterings and a
bulk viscosity. It will be interesting to see the application of our initial
state model to asymmetric collision systems like Cu+Au, 3He+Au and d+Au and
see how the boost invariance is broken in this framework. One should also
investigate how to include flow in the initial state and try to probe the
rapidity dynamics at lower energies. We also present our predictions for
higher flow harmonics which can be compared to the data when they become
available [59]. Of course, we as a community also need to find a way to
simulate hydrodynamics across a critical point and across a first order phase
transition, and only then we will have a full physics understanding of heavy
ion collisions at finite baryon densities. We will report on progress in this
direction in forthcoming publications.

## Acknowledgments

This work was supported by the U.S. DOE Grant No. DE-FG02-87ER40328.

##  Appendix A Energy loss distribution in 3D initial state

The normalized collision kernel in LEXUS is given in Eq. (3). The probability
of two nucleons with rapidities $y^{\prime}_{p}$ and $y^{\prime}_{T}$
colliding to give nucleons with rapidities $y_{P}$ and $y_{T}$

$\displaystyle P(y^{\prime}_{P}+y^{\prime}_{T}\rightarrow y_{P}+y_{T})$
$\displaystyle\propto$
$\displaystyle\cosh(y_{P}-y^{\prime}_{T})\cosh(y^{\prime}_{P}-y_{T})$
$\displaystyle=$
$\displaystyle\frac{1}{2}(\cosh(y_{P}-y^{\prime}_{T}+y^{\prime}_{P}-y_{T})$
$\displaystyle+\cosh(y_{P}-y^{\prime}_{T}-y^{\prime}_{P}+y_{T}))$

Recall that

$\displaystyle y^{\text{total}}_{\text{rest-frame}}$ $\displaystyle=$
$\displaystyle y^{\prime}_{P}-y^{\prime}_{T}$ (68) $\displaystyle
y_{\text{loss}}$ $\displaystyle=$ $\displaystyle
y^{\prime}_{P}-y^{\prime}_{T}-(y_{P}-y_{T})$ (69)

Also, in pair rest frame, $y^{\prime}_{P}+y^{\prime}_{T}=0$ and as both the
nucleons lose $y_{\text{loss}}/2$ units in rapidity, $y_{P}+y_{T}=0$. So,

$\displaystyle P(y_{\text{loss}})$ $\displaystyle\propto$
$\displaystyle\cosh(y_{P}-y^{\prime}_{T})\cosh(y^{\prime}_{P}-y_{T})$ (70)
$\displaystyle=$
$\displaystyle\frac{1}{2}\left(\cosh(2y^{\text{total}}_{\text{rest-
frame}}-y_{\text{loss}})+1\right)$ $\displaystyle\approx$
$\displaystyle\frac{1}{2}\cosh(2y^{\text{total}}_{\text{rest-
frame}}-y_{\text{loss}})$

Normalizing this gives us the distribution in Eq. (12).

##  Appendix B Comparing Thermal Conductivities

Reference [55] uses the baryon diffusion constant

$\kappa_{B}={\textstyle{\frac{1}{3}}}\tau
n_{B}\left[\coth\left(\frac{\mu_{B}}{T}\right)-\frac{3Tn_{B}}{w}\right]$ (71)

where $\tau$ is an energy independent relaxation time common to all particles.
For numerical studies they take $\tau=C_{B}/T$ and vary the dimensionless
parameter $C_{B}$. Where does this expression for $\kappa_{B}$ come from and
how does it compare to ours?

First note the relationship

$\kappa_{B}=\left(\frac{Tn_{B}}{w}\right)^{2}\lambda$ (72)

Expression (71) begins with quantum statistics. In that case [41]

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left(\frac{w}{n_{B}T}\right)^{2}\sum_{a}\int
d\Gamma_{a}\frac{p^{2}}{E_{a}}\frac{\tau_{a}(E_{a})}{E_{a}}$ (73)
$\displaystyle\times$
$\displaystyle\left(b_{a}-\frac{n_{B}}{w}E_{a}\right)^{2}f_{a}^{\rm
eq}(1-f_{a}^{\rm eq})$

This expression already enforces the condition of fit. When the relaxation
time is the same constant for all particles the particular solution Eq. (43)
automatically satisfies the condition of fit without any need for a nonzero
additive constant $b$ even with quantum statistics. To see that rewrite

$\displaystyle\lambda=\frac{\tau}{3}\left(\frac{w}{n_{B}T}\right)^{2}\sum_{a}b_{a}\int
d\Gamma_{a}\frac{p^{2}}{E_{a}^{2}}\left(b_{a}-\frac{n_{B}}{w}E_{a}\right)f_{a}^{\rm
eq}(1-f_{a}^{\rm eq})$
$\displaystyle-\frac{\tau}{3}\left(\frac{w}{n_{B}T^{2}}\right)\sum_{a}\int
d\Gamma_{a}\frac{p^{2}}{E_{a}}\left(b_{a}-\frac{n_{B}}{w}E_{a}\right)f_{a}^{\rm
eq}(1-f_{a}^{\rm eq})$ (74)

For an equilibrium Fermi-Dirac distribution function

$T\frac{\partial}{\partial\mu_{B}}f_{a}^{\rm eq}=b_{a}f_{a}^{\rm
eq}(1-f_{a}^{\rm eq})$ (75)

and

$T\left(T\frac{\partial}{\partial
T}+\mu_{B}\frac{\partial}{\partial\mu_{B}}\right)f_{a}^{\rm
eq}=E_{a}f_{a}^{\rm eq}(1-f_{a}^{\rm eq})$ (76)

so that the second contribution to $\lambda$ above is zero. Hence one can use
the simpler expression

$\lambda=\frac{\tau}{3}\left(\frac{w}{n_{B}T}\right)^{2}\sum_{a}b_{a}\int
d\Gamma_{a}\frac{p^{2}}{E_{a}^{2}}\left(b_{a}-\frac{n_{B}}{w}E_{a}\right)f_{a}^{\rm
eq}(1-f_{a}^{\rm eq})$ (77)

Now

$\sum_{a}b_{a}\int d\Gamma_{a}\frac{p^{2}}{E_{a}}f_{a}^{\rm eq}(1-f_{a}^{\rm
eq})=3Tn_{B}$ (78)

whereas the other integral

$\sum_{a}b_{a}^{2}\int d\Gamma_{a}\frac{p^{2}}{E_{a}^{2}}f_{a}^{\rm
eq}(1-f_{a}^{\rm eq})$

cannot readily be expressed in terms of thermodynamic functions. Therefore
Ref. [55] assumed massless particles to express the integral as

$\sum_{a}b_{a}^{2}\int d\Gamma_{a}f_{a}^{\rm eq}(1-f_{a}^{\rm
eq})=T\chi_{\mu\mu}$ (79)

which results in

$\kappa_{B}={\textstyle{\frac{1}{3}}}\tau
n_{B}\left[\frac{T\chi_{\mu\mu}}{n_{B}}-\frac{3Tn_{B}}{w}\right]$ (80)

The above formula is also true for massless particles obeying Boltzmann
statistics. For massless quarks obeying quantum statistics which all have the
same chemical potential

$\frac{T\chi_{\mu\mu}}{n_{B}}=\left(\frac{1+\mu_{B}^{2}/3\pi^{2}T^{2}}{1+\mu_{B}^{2}/9\pi^{2}T^{2}}\right)\frac{T}{\mu_{B}}$
(81)

whereas for massless quarks obeying Boltzmann statistics

$\frac{T\chi_{\mu\mu}}{n_{B}}=\coth\left(\frac{\mu_{B}}{T}\right)$ (82)

It is the latter which is used in [55], resulting in Eq. (71). In contrast,
our expression is

$\kappa_{B}=\tau^{\prime}\frac{T^{4}}{w}\det\chi\;v_{\sigma}^{2}=\frac{\eta}{2}\frac{T^{3}}{w^{2}}\det\chi\;v_{\sigma}^{2}$
(83)

## References

  * Aoki  _et al._ [2006] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, The Order of the quantum chromodynamics transition predicted by the standard model of particle physics, Nature 443, 675 (2006), arXiv:hep-lat/0611014 .
  * Ding  _et al._ [2015] H.-T. Ding, F. Karsch, and S. Mukherjee, Thermodynamics of strong-interaction matter from Lattice QCD, Int. J. Mod. Phys. E 24, 1530007 (2015), arXiv:1504.05274 [hep-lat] .
  * Bazavov  _et al._ [2019] A. Bazavov, F. Karsch, S. Mukherjee, and P. Petreczky (USQCD), Hot-dense Lattice QCD: USQCD whitepaper 2018, Eur. Phys. J. A 55, 194 (2019), arXiv:1904.09951 [hep-lat] .
  * Fukushima and Hatsuda [2011] K. Fukushima and T. Hatsuda, The phase diagram of dense QCD, Rept. Prog. Phys. 74, 014001 (2011), arXiv:1005.4814 [hep-ph] .
  * Fukushima and Sasaki [2013] K. Fukushima and C. Sasaki, The phase diagram of nuclear and quark matter at high baryon density, Prog. Part. Nucl. Phys. 72, 99 (2013), arXiv:1301.6377 [hep-ph] .
  * Fischer [2019] C. S. Fischer, QCD at finite temperature and chemical potential from Dyson–Schwinger equations, Prog. Part. Nucl. Phys. 105, 1 (2019), arXiv:1810.12938 [hep-ph] .
  * Heinz and Snellings [2013] U. Heinz and R. Snellings, Collective flow and viscosity in relativistic heavy-ion collisions, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826 [nucl-th] .
  * Braun-Munzinger  _et al._ [2016] P. Braun-Munzinger, V. Koch, T. Schäfer, and J. Stachel, Properties of hot and dense matter from relativistic heavy ion collisions, Phys. Rept. 621, 76 (2016), arXiv:1510.00442 [nucl-th] .
  * Romatschke and Romatschke [2019] P. Romatschke and U. Romatschke,  _Relativistic Fluid Dynamics In and Out of Equilibrium_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2019) arXiv:1712.05815 [nucl-th] .
  * Shen and Yan [2020] C. Shen and L. Yan, Recent development of hydrodynamic modeling in heavy-ion collisions, Nucl. Sci. Tech. 31, 122 (2020), arXiv:2010.12377 [nucl-th] .
  * Romatschke [2010] P. Romatschke, New Developments in Relativistic Viscous Hydrodynamics, Int. J. Mod. Phys. E 19, 1 (2010), arXiv:0902.3663 [hep-ph] .
  * Ollitrault [2008] J.-Y. Ollitrault, Relativistic hydrodynamics for heavy-ion collisions, Eur. J. Phys. 29, 275 (2008), arXiv:0708.2433 [nucl-th] .
  * Miller  _et al._ [2007] M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Glauber modeling in high energy nuclear collisions, Ann. Rev. Nucl. Part. Sci. 57, 205 (2007), arXiv:nucl-ex/0701025 .
  * Drescher and Nara [2007a] H. J. Drescher and Y. Nara, Effects of fluctuations on the initial eccentricity from the Color Glass Condensate in heavy ion collisions, Phys. Rev. C 75, 034905 (2007a), arXiv:nucl-th/0611017 .
  * Drescher and Nara [2007b] H.-J. Drescher and Y. Nara, Eccentricity fluctuations from the color glass condensate at RHIC and LHC, Phys. Rev. C 76, 041903 (2007b), arXiv:0707.0249 [nucl-th] .
  * Krasnitz and Venugopalan [2000] A. Krasnitz and R. Venugopalan, The Initial energy density of gluons produced in very high-energy nuclear collisions, Phys. Rev. Lett. 84, 4309 (2000), arXiv:hep-ph/9909203 .
  * Krasnitz and Venugopalan [2001] A. Krasnitz and R. Venugopalan, The Initial gluon multiplicity in heavy ion collisions, Phys. Rev. Lett. 86, 1717 (2001), arXiv:hep-ph/0007108 .
  * Pang  _et al._ [2012] L. Pang, Q. Wang, and X.-N. Wang, Effects of initial flow velocity fluctuation in event-by-event (3+1)D hydrodynamics, Phys. Rev. C 86, 024911 (2012), arXiv:1205.5019 [nucl-th] .
  * Schenke  _et al._ [2012] B. Schenke, P. Tribedy, and R. Venugopalan, Fluctuating Glasma initial conditions and flow in heavy ion collisions, Phys. Rev. Lett. 108, 252301 (2012), arXiv:1202.6646 [nucl-th] .
  * Moreland  _et al._ [2015] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions, Phys. Rev. C 92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Huovinen and Petreczky [2010] P. Huovinen and P. Petreczky, QCD Equation of State and Hadron Resonance Gas, Nucl. Phys. A 837, 26 (2010), arXiv:0912.2541 [hep-ph] .
  * Moreland and Soltz [2016] J. S. Moreland and R. A. Soltz, Hydrodynamic simulations of relativistic heavy-ion collisions with different lattice quantum chromodynamics calculations of the equation of state, Phys. Rev. C 93, 044913 (2016), arXiv:1512.02189 [nucl-th] .
  * Albright  _et al._ [2014] M. Albright, J. Kapusta, and C. Young, Matching Excluded Volume Hadron Resonance Gas Models and Perturbative QCD to Lattice Calculations, Phys. Rev. C 90, 024915 (2014), arXiv:1404.7540 [nucl-th] .
  * Cooper and Frye [1974] F. Cooper and G. Frye, Comment on the Single Particle Distribution in the Hydrodynamic and Statistical Thermodynamic Models of Multiparticle Production, Phys. Rev. D 10, 186 (1974).
  * An  _et al._ [2022] X. An _et al._ , The BEST framework for the search for the QCD critical point and the chiral magnetic effect, Nucl. Phys. A 1017, 122343 (2022), arXiv:2108.13867 [nucl-th] .
  * Shen and Schenke [2018] C. Shen and B. Schenke, Dynamical initial state model for relativistic heavy-ion collisions, Phys. Rev. C 97, 024907 (2018), arXiv:1710.00881 [nucl-th] .
  * Du  _et al._ [2019] L. Du, U. Heinz, and G. Vujanovic, Hybrid model with dynamical sources for heavy-ion collisions at BES energies, Nucl. Phys. A 982, 407 (2019), arXiv:1807.04721 [nucl-th] .
  * Shen  _et al._ [2017] C. Shen, G. Denicol, C. Gale, S. Jeon, A. Monnai, and B. Schenke, A hybrid approach to relativistic heavy-ion collisions at the RHIC BES energies, Nucl. Phys. A 967, 796 (2017), arXiv:1704.04109 [nucl-th] .
  * Akamatsu  _et al._ [2018] Y. Akamatsu, M. Asakawa, T. Hirano, M. Kitazawa, K. Morita, K. Murase, Y. Nara, C. Nonaka, and A. Ohnishi, Dynamically integrated transport approach for heavy-ion collisions at high baryon density, Phys. Rev. C 98, 024909 (2018), arXiv:1805.09024 [nucl-th] .
  * Shen and Alzhrani [2020] C. Shen and S. Alzhrani, Collision-geometry-based 3D initial condition for relativistic heavy-ion collisions, Phys. Rev. C 102, 014909 (2020), arXiv:2003.05852 [nucl-th] .
  * Okai  _et al._ [2017] M. Okai, K. Kawaguchi, Y. Tachibana, and T. Hirano, New approach to initializing hydrodynamic fields and mini-jet propagation in quark-gluon fluids, Phys. Rev. C 95, 054914 (2017), arXiv:1702.07541 [nucl-th] .
  * Jeon and Kapusta [1997] S. Jeon and J. I. Kapusta, Linear extrapolation of ultrarelativistic nucleon-nucleon scattering to nucleus-nucleus collisions, Phys. Rev. C 56, 468 (1997), arXiv:nucl-th/9703033 .
  * Hegde [2014] P. Hegde (BNL–Bielefeld–CCNU), The QCD equation of state to $\mathcal{O}(\mu_{B}^{4})$ from lattice QCD, Nucl. Phys. A 931, 851 (2014), arXiv:1408.6305 [hep-lat] .
  * Guenther  _et al._ [2017] J. N. Guenther, R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz, A. Pasztor, C. Ratti, and K. K. Szabó, The QCD equation of state at finite density from analytical continuation, Nucl. Phys. A 967, 720 (2017), arXiv:1607.02493 [hep-lat] .
  * Bazavov  _et al._ [2017] A. Bazavov _et al._ , The QCD Equation of State to $\mathcal{O}(\mu_{B}^{6})$ from Lattice QCD, Phys. Rev. D 95, 054504 (2017), arXiv:1701.04325 [hep-lat] .
  * Monnai  _et al._ [2019] A. Monnai, B. Schenke, and C. Shen, Equation of state at finite densities for QCD matter in nuclear collisions, Phys. Rev. C 100, 024907 (2019), arXiv:1902.05095 [nucl-th] .
  * Noronha-Hostler  _et al._ [2019] J. Noronha-Hostler, P. Parotto, C. Ratti, and J. M. Stafford, Lattice-based equation of state at finite baryon number, electric charge and strangeness chemical potentials, Phys. Rev. C 100, 064910 (2019), arXiv:1902.06723 [hep-ph] .
  * Borsányi  _et al._ [2021] S. Borsányi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pásztor, C. Ratti, and K. K. Szabó, Lattice QCD equation of state at finite chemical potential from an alternative expansion scheme, Phys. Rev. Lett. 126, 232001 (2021), arXiv:2102.06660 [hep-lat] .
  * Teaney [2003] D. Teaney, The Effects of viscosity on spectra, elliptic flow, and HBT radii, Phys. Rev. C 68, 034913 (2003), arXiv:nucl-th/0301099 .
  * Bozek [2010] P. Bozek, Bulk and shear viscosities of matter created in relativistic heavy-ion collisions, Phys. Rev. C 81, 034909 (2010), arXiv:0911.2397 [nucl-th] .
  * Albright and Kapusta [2016] M. Albright and J. I. Kapusta, Quasiparticle Theory of Transport Coefficients for Hadronic Matter at Finite Temperature and Baryon Density, Phys. Rev. C 93, 014903 (2016), arXiv:1508.02696 [nucl-th] .
  * Videbaek and Hansen [1995] F. Videbaek and O. Hansen, Baryon rapidity loss and mid-rapidity stacking in high-energy nucleus-nucleus collisions, Phys. Rev. C 52, 2684 (1995).
  * Hwa [1984] R. C. Hwa, Degradation of Proton Momentum Through Nuclei, Phys. Rev. Lett. 52, 492 (1984).
  * Csernai and Kapusta [1984] L. P. Csernai and J. I. Kapusta, PROTON STOPPING POWER OF HEAVY NUCLEI, Phys. Rev. D 29, 2664 (1984).
  * Florkowski [2010] W. Florkowski,  _Phenomenology of Ultra-Relativistic Heavy-Ion Collisions_ (2010).
  * Carruthers and Doung-van [1973] P. Carruthers and M. Doung-van, Rapidity and angular distributions of charged secondaries according to the hydrodynamical model of particle production, Phys. Rev. D 8, 859 (1973).
  * Schenke  _et al._ [2010] B. Schenke, S. Jeon, and C. Gale, (3+1)D hydrodynamic simulation of relativistic heavy-ion collisions, Phys. Rev. C 82, 014903 (2010), arXiv:1004.1408 [hep-ph] .
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2020, 083C01 (2020).
  * Vuorinen [2003] A. Vuorinen, The Pressure of QCD at finite temperatures and chemical potentials, Phys. Rev. D 68, 054017 (2003), arXiv:hep-ph/0305183 .
  * Chakraborty and Kapusta [2011] P. Chakraborty and J. I. Kapusta, Quasiparticle theory of shear and bulk viscosities of hadronic matter, Phys. Rev. C 83, 014906 (2011).
  * S. R. de Groot, W. A. van Leeuwen and Ch. G. van Weert [1980] S. R. de Groot, W. A. van Leeuwen and Ch. G. van Weert,  _Relativistic Kinetic Theory: Principles and Applications_ (North-Holland, Amsterdam, 1980).
  * Chakraborty and Kapusta [2017] P. Chakraborty and J. I. Kapusta, Departure from equilibrium of the quasiparticle distribution functions in high energy nuclear collisions, Phys. Rev. C 95, 014907 (2017).
  * Kapusta and Young [2014] J. I. Kapusta and C. Young, Causal baryon diffusion and colored noise, Phys. Rev. C 90, 044902 (2014).
  * Kapusta and Plumberg [2018] J. I. Kapusta and C. Plumberg, Causal electric charge diffusion and balance functions in relativistic heavy ion collisions, Phys. Rev. C 97, 014906 (2018), [Erratum: Phys.Rev.C 102, 019901 (2020)].
  * G. S. Denicol, C. Gale, S. Jeon, A. Monnai, B. Schenke, and C. Shen [2018] G. S. Denicol, C. Gale, S. Jeon, A. Monnai, B. Schenke, and C. Shen, Net baryon diffusion in fluid dynamic simulations of relativistic heavy-ion collisions, Phys. Rev. C 98, 034916 (2018).
  * Adare  _et al._ [2016a] A. Adare _et al._ (PHENIX), Transverse energy production and charged-particle multiplicity at midrapidity in various systems from $\sqrt{s_{NN}}=7.7$ to 200 GeV, Phys. Rev. C 93, 024901 (2016a), arXiv:1509.06727 [nucl-ex] .
  * Adamczyk  _et al._ [2016] L. Adamczyk _et al._ (STAR), Centrality dependence of identified particle elliptic flow in relativistic heavy ion collisions at $\sqrt{s_{NN}}$=7.7–62.4 GeV, Phys. Rev. C 93, 014907 (2016), arXiv:1509.08397 [nucl-ex] .
  * Adare  _et al._ [2016b] A. Adare _et al._ (PHENIX), Measurement of the higher-order anisotropic flow coefficients for identified hadrons in Au$+$Au collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV, Phys. Rev. C 93, 051902 (2016b), arXiv:1412.1038 [nucl-ex] .
  * Parfenov [2020] P. Parfenov (STAR), Elliptic ($v_{2}$) and triangular ($v_{3}$) anisotropic flow of identified hadrons from the STAR Beam Energy Scan program, J. Phys. Conf. Ser. 1690, 012128 (2020), arXiv:2012.06759 [hep-ex] .
  * Abelev  _et al._ [2009] B. I. Abelev _et al._ (STAR), Systematic Measurements of Identified Particle Spectra in $pp,d^{+}$ Au and Au+Au Collisions from STAR, Phys. Rev. C 79, 034909 (2009), arXiv:0808.2041 [nucl-ex] .
  * Adamczyk  _et al._ [2017] L. Adamczyk _et al._ (STAR), Bulk Properties of the Medium Produced in Relativistic Heavy-Ion Collisions from the Beam Energy Scan Program, Phys. Rev. C 96, 044904 (2017), arXiv:1701.07065 [nucl-ex] .
  * Adam  _et al._ [2020] J. Adam _et al._ (STAR), Bulk properties of the system formed in $Au+Au$ collisions at $\sqrt{s_{\mathrm{NN}}}$ =14.5 GeV at the BNL STAR detector, Phys. Rev. C 101, 024905 (2020), arXiv:1908.03585 [nucl-ex] .
  * Ryu  _et al._ [2015] S. Ryu, J. F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon, and C. Gale, Importance of the Bulk Viscosity of QCD in Ultrarelativistic Heavy-Ion Collisions, Phys. Rev. Lett. 115, 132301 (2015), arXiv:1502.01675 [nucl-th] .
  * Kapusta [1977] J. I. Kapusta, Particle Production in the Nuclear Fireball Model, Phys. Rev. C 16, 1493 (1977).
  * Li and Kapusta [2019] M. Li and J. I. Kapusta, Large Baryon Densities Achievable in High Energy Heavy Ion Collisions Outside the Central Rapidity Region, Phys. Rev. C 99, 014906 (2019), arXiv:1808.05751 [nucl-th] .