# Sensitivity of Au+Au collisions to the symmetric nuclear matter equation of
state
at 2–5 nuclear saturation densities

Dmytro Oliinychenko<EMAIL_ADDRESS>Institute for Nuclear
Theory, University of Washington, Box 351550, Seattle, Washington 98195, USA
Agnieszka Sorensen<EMAIL_ADDRESS>Institute for Nuclear Theory,
University of Washington, Box 351550, Seattle, Washington 98195, USA  Volker
Koch  Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley,
California 94720, USA  Larry McLerran  Institute for Nuclear Theory,
University of Washington, Box 351550, Seattle, Washington 98195, USA

###### Abstract

We demonstrate that proton and pion flow measurements in heavy-ion collisions
at incident energies ranging from 1 to 20 GeV per nucleon in the fixed target
frame can be used for an accurate determination of the symmetric nuclear
matter equation of state at baryon densities equal 2–4 times nuclear
saturation density $n_{0}$. We simulate Au+Au collisions at these energies
using a hadronic transport model with an adjustable vector mean-field
potential dependent on baryon density $n_{B}$. We show that the mean field can
be parametrized to reproduce a given density dependence of the speed of sound
at zero temperature $c_{s}^{2}(n_{B},T=0)$, which we vary independently in
multiple density intervals to probe the differential sensitivity of heavy-ion
observables to the equation of state at these specific densities. Recent flow
data from the STAR experiment at the center-of-mass energies
$\sqrt{s_{\textrm{NN}}}=\\{3.0,4.5\\}\ \rm{GeV}$ can be described by our
model, and a Bayesian analysis of these data indicates a hard equation of
state at $n_{B}\in(2,3)n_{0}$ and a possible phase transition at
$n_{B}\in(3,4)n_{0}$. More data at
$\sqrt{s_{\textrm{NN}}}=2$–$5~{}\textrm{GeV}$, as well as a more thorough
analysis of the model systematic uncertainties will be necessary for a more
precise conclusion.

##  I Introduction

The baryon density is approximately $n_{B}\approx 0.16~{}\textrm{fm}^{-3}$ at
the center of a nucleus and averages to about $0.12~{}\textrm{fm}^{-3}$ over
its entire volume, almost independently of the size of the given nucleus Wong
(1998). Neglecting the finite size effects in nuclei as well as Coulomb
interactions, one arrives at an idealized theoretical concept of nuclear
matter, which is in equilibrium at density $n_{0}\approx 0.16$ fm-3, often
called the nuclear matter saturation density Wong (1998); Bethe (1971). The
only way of obtaining a substantially more dense nuclear matter in a
laboratory is to collide heavy nuclei at relativistic incident velocities.
Such collisions of two nuclei result in a rapid (timescales on the order of a
few fm/$c$) compression and heating, followed by expansion and cooling of the
produced fireball. The outcomes of both the compression and the expansion
phase depend on the strong interactions between the constituents of the
fireball. In equilibrium, the influence of these interactions on the
properties of the medium is described by the equation of state (EOS), that is
by the dependence of the equilibrium pressure $P$ on temperature $T$, net
baryon density $n_{B}$, net strangeness density $n_{S}$, and net charge
density $n_{Q}$111Depending on the problem at hand, it is sometimes more
convenient to represent the EOS in the form $P(T,\mu_{B},\mu_{S},\mu_{Q})$,
where $\mu_{i}$ denotes the chemical potential associated with the conserved
charge $n_{i}$ and $i\in\\{B,S,Q\\}$; in the form
$P(\mathcal{E},n_{B},n_{S},n_{Q})$, where $\mathcal{E}$ is the energy density;
or in the form $P(\mathcal{E},n_{B},n_{S},n_{I3})$, where instead of the
charge density $n_{Q}$ one considers the isospin projection density $n_{I3}$,
with the two densities closely related through the relation between the
electric charge and the isospin projection of a hadron,
$Q=I_{3}+\frac{1}{2}(B+S)$..

We note that while the EOS is only defined in equilibrium, the systems created
in heavy-ion collisions may well be out-of-equilibrium. To incorporate non-
equilibrium effects in simulations, one can either introduce viscous
corrections to a fluid-dynamic description, which works well for small
deviations from equilibrium, or one can resort to transport theory, which
accounts for the full non-equilibrium dynamics governed by the Boltzmann
equation. Within transport-theoretical approaches, the model parameters (in
particular, the mean-field interaction) can be adjusted to reproduce relevant
observables measured in heavy-ion reactions such as flow. The EOS is then
obtained by using these model parameters in the transport equation, which
under the assumption that the system is in equilibrium can be used to
calculate, e.g., the pressure.

Ultimately, one of the major goals of heavy-ion collision experiments is to
extract the EOS within the experimentally accessible domain. This domain is,
admittedly, limited; for example, $n_{S}\approx 0$ in heavy-ion collisions.
Even more importantly, compression is always accompanied by heating, and as a
result regions characterized by high $n_{B}$ and low $T$ are not probed at any
collision energy Arsene  _et al._ (2007). This can be seen in Fig. 1, which
shows phase trajectories of the central region of a heavy-ion collision at
different values of the incident kinetic energy per nucleon (excluding the
rest mass) $E_{\textrm{lab}}$, obtained from simulations using the hadronic
transport code SMASH Weil  _et al._ (2016) (version 2.1 sma ).

Figure 1: Trajectories of a central region (i.e., a sphere of radius
$r=2~{}\textrm{fm}$ located at the event-averaged center-of-mass of the
collision system) in Au+Au collisions at impact parameter $b=6~{}\rm{fm}$,
obtained from simulations in SMASH. The contours demonstrate qualitatively
which ranges of densities are probed at which collision energies. The
simulations are performed in the center-of-mass frame and the time $t$ is also
given in that frame; the time difference between the shown data points is 1
fm/$c$, and the first points shown are for $t=-2$ fm/$c$ when nuclei are not
yet touching (by the convention used in SMASH, $t=0$ is the time at which
nuclei would touch in a central collision; note that for mid-central
collisions, the center of mass of the system at $t=0$ is between the nuclei,
which explains the corresponding low values of density). For collisions at
$E_{\textrm{lab}}=1$ and 2 GeV/nucleon (red and blue marks, respectively) the
trajectories are traversed in the counterclockwise direction with increasing
time, while for collisions at 4 GeV/nucleon and above, the trajectories are
traversed in the clockwise direction. The energy density $\mathcal{E}$ is
calculated from the $00$-component of the energy-momentum tensor $T^{\mu\nu}$
in the Landau frame; note that mean-field contributions are not taken into
account in this calculation of $T^{\mu\nu}$, that is
$\mathcal{E}=\mathcal{E}_{\textrm{kin}}$, to enable a clearer inference of the
different temperatures reached in collisions at different energies. The
simulations used a mean-field potential parametrized to reproduce the standard
Skyrme EOS at $n_{B}\in[0,2]n_{0}$ and requiring that
$c_{s}^{2}(n_{B},T=0)=0.3$ at higher densities (see Sec. II.2 for details of
the EOS construction). Changes in the trajectories due to employing different
EOSs are shown in Fig. 2.  Figure 2: Phase trajectories in a Au+Au collisions
at $E_{\textrm{lab}}=4~{}\textrm{GeV/nucleon}$, obtained using the same
simulation setup as in Fig. 1, demonstrating that varying the EOS only in the
density region $n_{B}\in(2,3)n_{0}$ leads to a higher (lower) reached values
of the density and energy per nucleon in the case of the softer (harder) EOS.
Here, in contrast to previous works where the stiffness of the EOS is
controlled solely by the value of nuclear matter incompressibility at
saturation $K_{0}$ (see, e.g., Danielewicz  _et al._ (2002); Fuchs  _et al._
(2001)), the stiffness of the EOS is changed by varying the value of the speed
of sound $c_{s}^{2}$ at $T=0$ for $n_{B}\in(2,3)n_{0}$, with small (large)
values of $c_{s}^{2}$ corresponding to a relatively soft (hard) EOS in that
density range. The trajectories are traversed in the clockwise direction with
increasing time, and the time difference between the shown data points is 1
fm/$c$. Below $n_{B}=2n_{0}$, the mean-field potential is parametrized to
reproduce the standard Skyrme EOS, while above $n_{B}=2n_{0}$, the mean-field
is parametrized to reproduce values of $c_{s}^{2}$ from a representative set,
$c_{s}^{2}(n_{B}>2n_{0},T=0)=\\{0.1,~{}0.5,~{}0.9\\}$ (see Sec. II for details
of the EOS construction).

The role of the EOS in heavy-ion collisions is intuitively understandable: a
stiffer EOS (meaning an EOS with a relatively large
$\partial{P}/\partial{n_{B}}$ or $\partial{P}/\partial{\mathcal{E}}$) leads to
less compression, less heating, and a faster transverse expansion than a
softer EOS (see, e.g., Fuchs  _et al._ (2001)). In Fig. 2, we not only
demonstrate that this is indeed the case in our simulations, but we also show
that there is a large sensitivity to changing the stiffness of the EOS only in
a specific density region (we will elaborate on this point in the following
sections). In terms of observables, the slower expansion characterizing a
softer EOS is expected to produce lower mean transverse momenta (see, e.g.,
Steinheimer  _et al._ (2022)), while the increased heating in this case leads
one to expect a larger thermal dilepton and photon yield (see, e.g., Huovinen
_et al._ (1999); Savchuk  _et al._ (2022); Li  _et al._ (2023)). Moreover, as
the fireball lifetime is longer for a softer EoS, one would also expect that
the combination of the femtoscopic radii
$R_{\textrm{out}}^{2}-R_{\textrm{side}}^{2}$, which can be shown to be
proportional to the duration of the emission of detected particles Bertsch
(1989); Pratt  _et al._ (1990), will be larger in this case.

Nevertheless, the most stringent currently available constraints on the
symmetric nuclear matter EOS come from angular distributions in the transverse
plane ${dN}/{d\phi}$, where $\phi$ denotes the azimuthal angle, which are
highly sensitive to the EOS and, at the same time, measurable with high
precision Danielewicz and Odyniec (1985). Indeed, the sensitivity of the flow
observables has been demonstrated by multiple hydrodynamic Stocker  _et al._
(1980); Ollitrault (1992); Rischke  _et al._ (1995); Stoecker (2005);
Brachmann  _et al._ (2000); Csernai and Rohrich (1999); Ivanov and Soldatov
(2015) and hadronic transport Hartnack  _et al._ (1994); Li and Ko (1998);
Danielewicz  _et al._ (2002); Le Fèvre  _et al._ (2016); Wang  _et al._
(2018); Colonna (2020); Nara and Ohnishi (2022) models. We notice, however,
that most of these works compare only several EOSs (often just an EOS with a
phase transition to a quark-gluon plasma to an EOS without such transition),
and do not attempt to quantify the sensitivity by parametrizing the EOS
continuously and constraining the parameters. The works that do explore a
range of possible EOSs Le Fèvre  _et al._ (2016); Wang  _et al._ (2018);
Danielewicz  _et al._ (2002) parametrize the EOS with a single parameter: the
incompressibility of isospin-symmetric nuclear matter at $n_{B}=n_{0}$ and
$T=0$, defined as $K_{0}\equiv 9\left({\partial P}/{\partial
n_{B}}\right)|_{n_{B}=n_{0}}$ or equivalently as
$K_{0}=9n_{B}^{2}\left({\partial^{2}(E/A)}/{\partial
n_{B}^{2}}\right)|_{n_{B}=n_{0}}$, where $E/A$ is energy per baryon. In this
work, we aim at exploring the sensitivity of the flow observables to the EOS
parametrized by varying both the incompressibility $K_{0}$ and values of the
speed of sound at different ranges of baryon density:
$c^{2}_{[2,3]n_{0}}\equiv c_{s}^{2}\big{[}n_{B}\in(2,3)n_{0}\big{]}$,
$c^{2}_{[3,4]n_{0}}\equiv c_{s}^{2}\big{[}n_{B}\in(3,4)n_{0}\big{]}$, and
$c^{2}_{[4,5]n_{0}}\equiv c_{s}^{2}\big{[}n_{B}\in(4,5)n_{0}\big{]}$. We
consider this work as the first step towards a full Bayesian analysis of all
available observables with a flexible and suitably parametrized EOS. Here,
however, we use only a limited set of measurements, that is flow measurements
as provided by the E895 Collaboration Pinkenburg  _et al._ (1999); Liu  _et
al._ (2000); Chung  _et al._ (2002); Klay  _et al._ (2002) and by recent
results from the STAR Collaboration Adam  _et al._ (2021); Abdallah  _et al._
(2022) based on Phase II of the Beam Energy Scan (BES) program at the
Relativistic Heavy Ion Collider (RHIC).

The effects of the EOS on ${dN}/{d\phi}$ are known to be substantial in
noncentral collisions at energies for which the speed of the spectators is
comparable to that of the fireball expansion Sorge (1997); Danielewicz  _et
al._ (2002). The ${dN}/{d\phi}$ distribution of protons around midrapidity
$y^{\prime}=0$ (where $y^{\prime}=y/y_{\textrm{beam}}$ is the center-of-mass
rapidity $y$ scaled by the beam rapidity in the center-of-mass frame
222Rapidity is defined as
$y(\beta)=\frac{1}{2}\ln\big{[}{(1+\beta)}/{(1-\beta)}\big{]}$, where $\beta$
is the velocity, and a Lorentz boost corresponding to the relative velocity of
two inertial frames $\beta_{0}$ is equivalent to a shift of $y$ by
$y(\beta_{0})$. Therefore, the center-of-mass rapidity $y$ for collisions in
the fixed-target mode can be obtained by taking $y\to y-y_{\textrm{beam}}$. In
the center-of-mass frame, rapidity distributions of measured particles are
centered around $y=0$, while the beam rapidities are at $\pm
y_{\textrm{beam}}$. To make results at different collision energies, and
therefore different spreads in rapidity, easier to compare, one uses the
scaled rapidity $y^{\prime}=y/y_{\textrm{beam}}$, since for the colliding
beams one then always has $y^{\prime}_{\textrm{beam}}=\pm 1$ in the center-of-
mass frame. Notice that this property of $y^{\prime}$ is satisfied only in the
center-of-mass frame. ) has maxima coincident with the reaction plane in the
case where the spectators move out of the way of the expanding fireball fast
enough. In the opposite case, where the spectators block the in-plane fireball
expansion, the preferential emission occurs in the out-of-plane direction (due
to the role of the spectators, this phenomenon is often called a “squeeze-
out”). This behavior is captured by the second Fourier coefficient of
${dN}/{d\phi}$,

$\displaystyle v_{2}=\langle\cos 2\phi\rangle=\frac{\int
d\phi~{}\cos(2\phi)~{}\frac{dN}{d\phi}}{\int d\phi~{}\frac{dN}{d\phi}}~{},$
(1)

known as the elliptic flow, which quantifies the difference between the in-
plane and out-of-plane emission. A positive $v_{2}$ indicates a preferential
emission toward angles $\phi\approx 0$ and $\phi\approx\pi$, while a negative
$v_{2}$ indicates a preferential emission toward $\phi\approx\pi/2$ and
$\phi\approx 3\pi/2$. For collisions at $E_{\textrm{lab}}\approx
1$–$10~{}\textrm{GeV/nucleon}$, for which the spectators still occupy the
vicinity of the collision region during the expansion phase, a faster fireball
expansion due to a stiffer EOS will be correspondingly more forcefully blocked
by the spectators, resulting in a larger squeeze-out effect and a more
negative $v_{2}$.

Another observable used to constrain the EOS is the slope of the directed flow
at mid-rapidity,

$\displaystyle\frac{dv_{1}}{dy^{\prime}}\bigg{|}_{y^{\prime}=0}=\frac{d\langle\cos\phi\rangle}{dy^{\prime}}\bigg{|}_{y^{\prime}=0}~{},$
(2)

where the directed flow $v_{1}$ is defined as the first Fourier coefficient of
${dN}/{d\phi}$, so that a positive $v_{1}=\langle\cos\phi\rangle$
characterizes a preferential emission in the $\phi\in[-\pi/2,\pi/2]$ direction
compared to the $\phi\in[\pi/2,3\pi/2]$ direction. Intuitively,
$dv_{1}/dy^{\prime}$ quantifies the strength of spectator deflection Voloshin
_et al._ (2010). To visualize the physical meaning of $dv_{1}/dy^{\prime}$,
one can imagine that spectators are two masses pushed apart by a spring
oriented along the transverse plane, where the spring represents the
compressed fireball. In the center-of-mass frame, the spectators moving in the
positive $y$ direction will be deflected towards $\phi\in[-\pi/2,\pi/2]$,
while the spectators moving in the negative $y$ direction will be deflected
towards $\phi\in[\pi/2,3\pi/2]$, resulting in a positive $v_{1}$ for $y>0$ and
a negative $v_{1}$ for $y<0$ 333Note that this picture has a certain
dependence on how one defines the coordinate axes in the problem. Here, we
define the impact parameter vector $\bm{b}$ (considered in the center-of-mass
frame, at the moment of the closest approach of the centers of the two nuclei)
as starting at the center of the heavy ion moving in the negative $y$
direction and ending at the center of the heavy-ion moving in the positive $y$
direction, and we define the $x$ axis of the transverse plane to be parallel
to the so-defined $\bm{b}$. The azimuthal angle $\phi$ is then measured from
the $x$ axis. A realization of this system of coordinates, in which the cross
product of the basis vectors of the transverse plane is parallel to the $z$
axis, is depicted in Fig. 1 of Danielewicz  _et al._ (2002). In general, by
convention, one defines the coordinate axes such that $dv_{1}/dy$ is positive
for spectators deflected away from the collision zone.. As can be easily
visualized using the spring analogy, a softer EOS results in a weaker
deflection and a smaller $dv_{1}/dy^{\prime}$. Moreover, it has been shown,
both in hydrodynamic and transport simulations, that a sufficiently soft EOS
(including EOSs with phase transitions) can lead to a negative
$dv_{1}/dy^{\prime}$ Stoecker (2005); Zhang  _et al._ (2018a), which within
our analogy corresponds to a situation in which the spring pulls the two
masses closer together.

One of the most prominent constraints on the EOS using the flow analysis of
heavy-ion collisions Danielewicz  _et al._ (2002) is based on comparing a
hadronic transport model Danielewicz and Bertsch (1991); Danielewicz (2000)
and flow data measured in the E895 experiment Pinkenburg  _et al._ (1999); Liu
_et al._ (2000). The main result of the work Danielewicz  _et al._ (2002) is
constraining the incompressibility $K_{0}$ of symmetric nuclear matter, the
only parameter of the EOS varied in that study, to be between 210 and 380 MeV.
(We note here that in studies applying such simple parametrizations of the EOS
to heavy-ion collisions at relativistic energies, which primarily probe large
baryon densities, the incompressibility is treated as a parameter which
specifies the behavior of the EOS at densities above the saturation density
$n_{0}$.) Such a spread in the extracted values of $K_{0}$ originates not from
data uncertainties, but from the fact that the used model was not able to
describe the $v_{2}$ and $dv_{1}/dy^{\prime}$ data simultaneously. While
reproducing the $v_{2}$ measurements required larger values of $K_{0}$ (and
therefore a harder EOS), the $dv_{1}/dy^{\prime}$ data required smaller values
of $K_{0}$ (corresponding to a softer EOS). In this work, we test whether a
relaxation of the EOS model that allows one to independently vary the
stiffness of the EOS in different density regions can lead to a more
consistent description of the flow data and, consequently, to a more precise
constraint on the EOS. For example, if densities to which $dv_{1}/dy^{\prime}$
is most sensitive are higher than the corresponding densities for $v_{2}$,
then we may make the EOS soft at high densities to reproduce the measured
$dv_{1}/dy^{\prime}$, and stiff at lower densities to reproduce the measured
$v_{2}$. Such idea was, in fact, suggested in Danielewicz  _et al._ (2002),
but the chosen parametrization of the EOS did not allow to implement it. In
this work, we construct a more flexible EOS parametrization and use it to fit
the E895 data available at the time of Danielewicz  _et al._ (2002) as well as
the newest STAR data for the same energy range Adam  _et al._ (2021); Abdallah
_et al._ (2022).

The structure of this work is as follows: Section II gives a short overview of
the dynamical evolution in hadronic transport (Sec. II.1), introduces a
parametrization of the mean-field potential reproducing a given behavior of
the speed of sound as a function of baryon density (Sec. II.2), and provides a
summary of simulation and implementation choices relevant to the study at hand
(Sec. II.3). Section III uses the developed model to explore the sensitivity
of the flow observables to the stiffness of the EOS in separate density
intervals. Section IV.1 presents our exploration of describing the available
flow data from E895 and STAR experiments, where we conclude that there is a
significant discrepancy between the results from the two data sets. Finally,
Sec. IV.2 discusses the results of our Bayesian analysis of the STAR flow
data, including a discussion of the slight tension with the EOS inferences
from neutron star observations. We summarize and provide an outlook in Section
V.

##  II Methodology

In the following we use natural units in which $\hbar=c=k_{B}=1$.

###  II.1 Simulation framework

As already discussed above, at projectile kinetic energies $E_{\textrm{lab}}$
below 20 GeV per nucleon in the fixed target frame, corresponding to
$\sqrt{s}_{\textrm{NN}}\lesssim 6.4~{}\textrm{GeV}$, spectators play a very
important role both for $dv_{1}/dy^{\prime}$ and for $v_{2}(y^{\prime}=0)$,
which was also explicitly demonstrated in hadronic transport simulations Zhang
_et al._ (2018a). This is in contrast to the highest RHIC
($\sqrt{s_{\textrm{NN}}}$ = 200 GeV) and LHC ($\sqrt{s_{\textrm{NN}}}$ = 2.76
and 5.02 TeV) energies, where mid-rapidity observables are unaffected by the
spectators. In consequence, most of the state-of-the-art hydrodynamic codes,
intended for very high energies and neglecting spectators, are not applicable
at $E_{\textrm{lab}}=2$–$20$ GeV/nucleon without modifications. Therefore, we
choose to employ a hadronic transport simulation in which spectators are
naturally included throughout the evolution. Our particular code of choice is
the transport code SMASH Weil  _et al._ (2016) (version 2.1 sma ), modified
according to the prescription given in the next subsection. SMASH is a
relativistic Boltzmann–Ueling–Uhlenbeck (BUU) type of hadronic transport with
vector-density dependent mean-field potentials, which means that it is a Monte
Carlo solver of the following kinetic equations

$\Pi_{\mu}\partial^{\mu}_{x}f_{i}(x,p)+\Pi_{\nu}\big{(}\partial_{\mu}^{x}A^{\nu}\big{)}\partial_{p}^{\mu}f_{i}(x,p)=I_{\textrm{coll}}^{(i)}\,.$
(3)

Here, $x^{\mu}$ and $p^{\mu}$ are position and momentum four-vectors (with
indices suppressed where it is beneficial for clarity), $\Pi^{\mu}(x,p)\equiv
p^{\mu}-A^{\mu}(x)$ is the kinetic momentum and $A^{\nu}$ is a vector field
dependent on the baryon current as

$\displaystyle A^{\mu}=\alpha(n_{B})j_{B}^{\mu}\,,$ (4)

where $\alpha(n_{B})$ is a chosen function of the rest frame density
$n_{B}\equiv\sqrt{(j_{B})_{\mu}(j_{B})^{\mu}}$, $f_{i}$ is the distribution
function for the particle species $i$, and $I_{\textrm{coll}}^{(i)}$ is the
collision integral for the $i$th particle species. The distribution function
$f_{i}$ depends on position as well as energy and momentum, and can be written
in a form that accounts for the mass-shell condition explicitly:

$\displaystyle
f_{i}(x^{\mu},p^{\mu})=2(2\pi)\Theta(\Pi^{0})\delta\big{(}\Pi^{\mu}\Pi_{\mu}-m_{i}^{2}\big{)}~{}\tilde{f}_{i}(x^{\mu},\bm{p})~{}.$
(5)

In the simulation, the continuous distribution function
$\tilde{f}_{i}(x,\bm{p})$ describing a system of $N_{i}$ particles of species
$i$ is approximated by using a standard test-particle ansatz, that is by
taking $\tilde{f}(x,\bm{p})$ to be of the form

$\displaystyle\tilde{f}_{i}(x,\bm{p})=\frac{1}{N_{T}}\sum_{j=1}^{N_{T}N_{i}}\delta^{(3)}\big{(}\bm{x}-\bm{x}_{j}(t)\big{)}\delta^{(3)}\big{(}\bm{p}-\bm{p}_{j}(t)\big{)}~{}.$
(6)

The above equation means that the continuous distribution function
$\tilde{f}_{i}(x,\bm{p})$, describing a physical system composed of $N_{i}$
particles, is sampled by $N=N_{i}N_{T}$ discrete points, or test particles,
where the factor $N_{T}\gg 1$ is known as the number of test particles per
particle or the oversampling factor. Intuitively, the larger the number of
samples (i.e., the larger $N_{T}$ is used) is, the better is the description
of $\tilde{f}_{i}(x,\bm{p})$. To preserve the essential properties of a system
of $N_{i}$ particles in a simulation involving $N_{T}N_{i}$ test particles
(such as the collision rate per test particle or the local density), all
collision cross sections are reduced by a factor of $N_{T}$, and a single
test-particle contributes to the energy density and charge densities with a
factor of $1/N_{T}$.

By substituting Eqs. (5) and (6) into the Vlasov equation (that is, the
Boltzmann equation, Eq. (3), with the collision term set to zero,
$I_{\textrm{coll}}^{(i)}$ = 0), one obtains the following equations of motion
for a single test particle in a system described by $f_{i}(x,p)$:

$\displaystyle\frac{dx^{\mu}}{dt}$
$\displaystyle=\frac{\Pi^{\mu}}{\Pi^{0}}~{},$ (7)
$\displaystyle\frac{d\Pi^{\mu}}{dt}$
$\displaystyle=\frac{\Pi_{\nu}}{\Pi^{0}}F^{\mu\nu}~{},$ (8)

where $F^{\mu\nu}=\partial^{\mu}A^{\nu}-\partial^{\nu}A^{\mu}$ and each test
particle satisfies the mass-shell condition $\Pi^{\mu}\Pi_{\mu}=m_{i}^{2}$;
details of this derivation can be found in Appendix A. Propagating the test
particles according to equations of motion, Eqs. (7)-(8), together with
performing decays and particle-particle collisions effectively solves Eq. (3),
as at each time $t$ the distribution function is given by (6). Importantly,
while Eqs. (7)-(8) have the same form as standard relativistic equations of
motion in an external vector field, in our case the vector field is directly
dependent on local baryon density, that is, the equations governing the system
are self-consistent.

An alternative derivation of Eqs. (7)-(8) within a relativistic vector density
functional (VDF) model of dense nuclear matter is shown in Sorensen and Koch
(2021), where the equilibrium thermodynamics and transport equations are
derived within a single approach. Here, in contrast to the approach taken in
Sorensen and Koch (2021) (where the contribution of the vector interactions to
the energy density is assumed to be a polynomial in baryon density), we
derived Eqs. (7)-(8) using the method of test-particles (see Appendix A),
which allows us to chose a completely arbitrary form of the vector single-
particle potential $A^{\mu}$; this property will be utilized in Sec. II.2.

By taking moments of the Boltzmann equation, Eq. (3), one can derive the
conservation laws in the standard way; details of this derivation are given in
Appendix B. In particular, the conserved current density of the $i$-th species
is given by

$\displaystyle j_{i}^{\mu}$
$\displaystyle=g_{i}\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{\Pi^{\mu}}{\Pi^{0}}~{}\tilde{f}_{i}~{},$
(9)

where $g_{i}$ is the degeneracy, and the net baryon current is defined as

$\displaystyle j_{B}^{\mu}\equiv\sum_{i}B_{i}j_{i}^{\mu}~{},$ (10)

where $B_{i}$ is the baryon number of the $i$th species. Meanwhile, the
energy-momentum tensor assumes the following form:

$\displaystyle T^{\mu\nu}$
$\displaystyle=\sum_{i}g_{i}\int\frac{d^{3}\Pi}{(2\pi)^{3}}~{}\frac{\Pi^{\mu}\Pi^{\nu}}{\Pi^{0}}~{}\tilde{f}_{i}+A^{\mu}j_{B}^{\nu}$
$\displaystyle\hskip
14.22636pt-~{}g^{\mu\nu}~{}\bigg{(}n_{B}U(n_{B})-\int_{0}^{n_{B}}dn^{\prime}~{}U(n^{\prime})\bigg{)}~{},$
(11)

where we introduce the notation $U(n_{B})$ to distinguish the zeroth component
vector potential $A^{\mu}$ calculated in the rest frame of the baryon current,
$U(n_{B})\equiv A^{0}(n_{B})\big{|}_{\begin{subarray}{c}\textrm{rest}\\\
\textrm{frame}\end{subarray}}=\alpha(n_{B})n_{B}$. Note, in particular, that
in equilibrium $T^{\mu\nu}=\textrm{diag}(\mathcal{E},P,P,P)$, giving us access
to the thermodynamic properties of the system.

The basis of our approach is the fact that by parametrizing the function
$\alpha(n_{B})$, introduced in Eq. (4), one can reproduce given properties of
nuclear matter in equilibrium and thus control the EOS of baryons, while at
the same time the mean-field interactions which lead to a given EOS also enter
the equations of motion, Eqs. (7)-(8), which continue to be well-defined for
an out-of-equilibrium evolution and are used in the simulations.

###  II.2 Parametrization of the mean-field potential by the speed of sound

In this section, we will use the formalism described above to introduce a
parametrization of the mean-field potential $A^{\mu}$ that reproduces a given
behavior of the speed of sound squared $c_{s}^{2}$ as a function of the baryon
density $n_{B}$.

In equilibrium, the baryon density of the $i$th baryonic species is given by

$\displaystyle
n_{B,i}=g_{i}\int\frac{d^{3}\Pi}{(2\pi)^{3}}~{}\left[e^{\beta\big{(}\Pi^{0}+\alpha(n_{B})n_{B}-\mu_{B}\big{)}}+1\right]^{-1}~{},$
(12)

where $\beta=1/T$ and $\Pi^{0}=\sqrt{\bm{\Pi}^{2}+m_{i}^{2}}$. At zero
temperature and in the absence of scalar interactions (which would lead to
effective masses $m_{i}^{*}<m_{i}$ and, consequently, smaller energies
required to excite more massive species), the only baryon species present in
the range of baryon densities relevant to heavy-ion collisions,
$n_{B}\in[0,5]n_{0}$, are protons and neutrons (nucleons); therefore, to
simplify the notation, in the following derivation we will drop the index $i$
and take $g=4$, $m=m_{N}=938~{}\textrm{MeV}$. Note that assuming only nucleons
to be present at $T=0$ does not preclude exciting other baryon states at
$T>0$; moreover, while within our formalism the EOS is essentially fixed for
nucleons at $T=0$, it still displays nontrivial behavior as a function of $T$,
and is applicable to complex systems of many baryonic species that inevitably
arise in considerations at finite temperature as well as in heavy-ion
collisions.

Taking the $T\to 0$ limit and integrating Eq. (12) over the Fermi sphere leads
to

$\displaystyle\mu_{B}(n_{B},T=0)$
$\displaystyle=\alpha(n_{B})n_{B}+\left[m_{N}^{2}+\left(\frac{6\pi^{2}n_{B}}{g}\right)^{2/3}\right]^{1/2}\,.$
(13)

At the same time, at $T=0$ the expression for the speed of sound squared is

$\displaystyle c_{s}^{2}(n_{B},T=0)$
$\displaystyle=\frac{n_{B}}{\mu_{B}\left(\frac{\partial{n_{B}}}{\partial{\mu_{B}}}\right)}~{}.$
(14)

Solving the above differential equation for $\mu_{B}$ yields

$\displaystyle\mu_{B}(n_{B},T=0)$
$\displaystyle=\mu_{B}\big{(}n_{B}^{(0)}\big{)}\exp\Bigg{[}\int_{n_{B}^{(0)}}^{n_{B}}dn^{\prime}~{}\frac{c_{s}^{2}(n^{\prime})}{n^{\prime}}\Bigg{]}~{},$
(15)

where $n_{B}^{(0)}$ is some density at which we know the corresponding value
of the chemical potential $\mu_{B}\big{(}n_{B}^{(0)}\big{)}$. Equating the
left-hand sides of Eqs. (13) and (15), we obtain the single-particle rest
frame potential $U(n_{B})\equiv
A^{0}(n_{B})\big{|}_{\begin{subarray}{c}\textrm{rest}\\\
\textrm{frame}\end{subarray}}=\alpha(n_{B})n_{B}$,

$\displaystyle U(n_{B})$
$\displaystyle=\mu_{B}\big{(}n_{B}^{(0)}\big{)}\exp\Bigg{[}\int_{n_{B}^{(0)}}^{n_{B}}dn^{\prime}~{}\frac{c_{s}^{2}(n^{\prime})}{n^{\prime}}\Bigg{]}$
$\displaystyle\hskip
14.22636pt-~{}\left[m_{N}^{2}+\left(\frac{6\pi^{2}n_{B}}{g}\right)^{2/3}\right]^{1/2}~{}.$
(16)

This form of $U(n_{B})$ allows one to parametrize the vector potential in an
intuitive way, that is via the density dependence of the speed of sound at
zero temperature.

For simplicity, in our study we choose the following piecewise functional form
of $c_{s}^{2}(n_{B})$:

$\displaystyle
c_{s}^{2}(n_{B})=\left.\begin{cases}c_{s}^{2}(\text{Skyrme}),&n_{B}<n_{1}=2n_{0}\\\
c_{1}^{2},&n_{1}<n_{B}<n_{2}\\\ c_{2}^{2},&n_{2}<n_{B}<n_{3}\\\ \dots&\\\
c_{m}^{2},&n_{m}<n_{B}\end{cases}\right.$ (17)

Such parametrization can be considered as a zeroth order interpolation of an
arbitrary function $c_{s}^{2}(n_{B})$. In the above, we take into account
that, for $n_{B}<2n_{0}$, an arbitrary parametrization would be inadequate as
in this region of baryon density the potential is already considerably well
constrained. Therefore, in the density range $n_{B}\in[0,2]n_{0}$, we adopt a
polynomial parametrization of the potential often referred to as the Skyrme
potential,

$\displaystyle
U_{\textrm{Sk}}(n_{B})=C_{1}\left(\frac{n_{B}}{\tilde{n}}\right)^{b_{1}-1}+C_{2}\left(\frac{n_{B}}{\tilde{n}}\right)^{b_{2}-1}\,,$
(18)

where $\tilde{n}=0.168$ fm-3, $C_{1}=-209.2$ MeV, $C_{2}=156.5$ MeV,
$b_{1}=2$, and $b_{2}=2.35$. The Skyrme potential with these values of the
parameters is the default mean field in SMASH, producing a nuclear matter
ground state at $n_{0}=0.166~{}\textrm{fm}^{-3}$ with binding energy per
nucleon $E_{\textrm{bin}}=-15.65$ MeV and a moderate incompressibility of
$K_{0}=236.73$ MeV (the small discrepancy between $n_{0}$ and $\tilde{n}$ is
an inaccuracy of the default SMASH parametrization, which, however, does not
affect any of our results). We note here that given $n_{0}$ and
$E_{\textrm{bin}}$, which do not vary significantly between different
approaches, whether the Skyrme EOS is soft or hard is entirely controlled by
the value of $K_{0}$, with small values of $K_{0}$ for soft EOSs and large
values of $K_{0}$ for hard EOSs. Unless stated otherwise, we employ this
default parametrization of the potential at $n_{B}<2n_{0}$. Above $2n_{0}$,
the potential is controlled by the speed of sound as can be seen in Eq. (17),
and by substituting this piecewise functional form of $c_{s}^{2}$ into Eq.
(16), we obtain

$\displaystyle
U(n_{B})=\left.\begin{cases}U_{\mathrm{Sk}}(n_{B})~{},&~{}~{}~{}n_{B}<n_{1}=2n_{0}\\\
\Big{[}U_{\mathrm{Sk}}(n_{1})+\mu^{*}(n_{1})\Big{]}\left(\frac{n_{B}}{n_{1}}\right)^{c_{1}^{2}}-\mu^{*}(n_{B})~{},&~{}~{}~{}n_{1}<n_{B}<n_{2}\\\
\Big{[}U_{\mathrm{Sk}}(n_{1})+\mu^{*}(n_{1})\Big{]}\left(\frac{n_{B}}{n_{k}}\right)^{c_{k}^{2}}\prod^{k}_{i=2}\left(\frac{n_{i}}{n_{i-1}}\right)^{c_{i-1}^{2}}-\mu^{*}(n_{B})~{},&~{}~{}~{}n_{k}<n_{B}<n_{k+1}\end{cases}\right.$
(19)

where we denote

$\displaystyle\mu^{*}(n_{B})\equiv\left[m_{N}^{2}+\left(\frac{6\pi^{2}n_{B}}{g}\right)^{2/3}\right]^{1/2}~{}.$
(20)

From Eq. (19), it is straightforward to obtain $\alpha(n_{B})$ which enters
the expression for the vector field, Eq. (4), with the latter in turn entering
the equations of motion, Eqs. (7)-(8).

###  II.3 Details of the simulations

We implement the vector mean field $A^{\mu}$, parametrized by the speed of
sound according to Eq. (19), by modifying the existing SMASH implementation of
the VDF model Sorensen and Koch (2021) which our model generalizes. The VDF
approach, where the EOS is parametrized using a polynomial form of
$\alpha(n_{B})$, is fully relativistically covariant and, in particular, leads
to fully covariant equations of motion, Eqs. (7)-(8) 444We note that, while
the mean-field equations are fully covariant, the implementation of the
Boltzmann equation with a finite number of test particles breaks the
covariance (for example, collisions between test particles occur at a finite
distance, the mean field is obtained on a grid with a finite spacing, etc.).
However, in the the limit of an infinite number of test particles these
effects disappear, and full covariance is restored.. The generalization of the
VDF model presented in this paper is, likewise, fully relativistically
covariant. In particular, the VDF vector mean-field used to reproduce the
Skyrme potential at $n_{B}<2n_{0}$, as described above, is fully
relativistically covariant as well. We stress this because usually Skyrme
potentials used in hadronic transport codes are, in contrast to our case,
nonrelativistic. One important advantage of using covariant mean fields is
that a Lorentz-contacted, boosted nucleus is still in the ground state. In a
nonrelativistic treatment, on the other hand, one has to either not allow for
Lorentz contraction of the incoming nuclei, which results in overestimating
the penetration time, or the Lorentz-contracted colliding nuclei are not in
the ground state, thus possibly underestimating the repulsion needed to
describe, e.g., flow data.

The calculation of the baryon density, necessary for solving the equations of
motion, is performed on a Cartesian lattice with lattice spacing 0.8 fm in
$x$, $y$, and $z$ directions. We have tested two different methods of weighing
(also called smearing) a test particle’s contribution to density at each node
of the lattice: Gaussian smearing and triangular smearing, described in
Appendix D.1. In the case of Gaussian smearing, the weight for a given test
particle’s contribution to density on the lattice is Lorentz-contracted in the
direction of that test particle’s motion. In the case of triangular smearing,
the lattice itself is contracted in the $z$ direction by a factor
$\gamma={\sqrt{s_{NN}}}/{2m_{N}}$. In Appendix D.2, we show that the flow
observables are independent of the type of smearing used for the density
calculation (provided that certain technical details of the density
calculation are done correctly; see the appendix for more details). To obtain
the results presented in the following, we have used the Gaussian smearing,
which had a slightly lower computational cost due to the smaller required
number of lattice nodes.

Throughout the simulations we apply a mixed ensembles method, within which we
simultaneously simulate $N_{\textrm{ens}}$ events (in this context also called
ensembles), each with $N_{T}$ test particles per particle. The test particles
can collide with each other only within one ensemble, and their cross sections
are reduced by a factor of $N_{T}$ so that the scattering rate per test
particle is preserved. However, the baryon density, the vector mean field, and
the estimated distribution function for Pauli blocking are computed based on
test particles in all $N_{\textrm{ens}}$ ensembles. We note that for a
reliable density calculation, the product $N_{T}N_{\textrm{ens}}$ should be
sufficiently large; see Appendix D.2.

It is important to stress here that we simulate multiple hadronic species and
their reactions. In Eq. (3), the index $i$ goes over all particle species
implemented in SMASH, that is over around 120 hadronic species (this number
does not account for isospin states, i.e., $\pi^{+}$, $\pi^{-}$, and $\pi^{0}$
are counted as one species; antibaryons and baryons are also counted as one
species). Most of these species are short-lived hadronic resonances. Masses,
width, decay channels, and branching ratios of the resonances are taken from
the Particle Data Group compendium Olive  _et al._ (2014) whenever available.
Unknown or poorly known properties were adjusted to better fit the known cross
sections Steinberg  _et al._ (2019). The collision integral on the right-hand
side of Eq. (3) includes elastic and inelastic $2\leftrightarrow 2$ reactions,
$1\leftrightarrow 2$ resonance formations and decays, as well as $2\to n$
reactions realized as a string excitation and breakup Mohs  _et al._ (2020).
Details about resonance properties and reaction cross sections can be found in
Weil  _et al._ (2016); Steinberg  _et al._ (2019).

We note that we do not utilize recently implemented stochastic collision rates
Staudenmaier  _et al._ (2021); Garcia-Montero  _et al._ (2022), resorting
instead to a standard geometric collision criterion within which collisions
are performed if $d_{ij}<\sqrt{\sigma_{ij}/\pi}$, where $d_{ij}$ is the
transverse distance between two particles in their center-of-mass frame,
computed at the time of the closest approach of the two particles, and
$\sigma_{ij}$ is the cross section. A consequence of this choice is that
collisions occur at a non-zero distance between particles. This causes
spurious effects described in detail in Cheng  _et al._ (2002): there, it was
shown that the finite interaction range effectively increases the viscosity
and, moreover, that it makes hadronic flow dependent on the technical details
of the collision implementation, which differ between various hadronic
transport codes. This problem is remedied by using a large number of test
particles per particle, $N_{T}$. As shown in Cheng  _et al._ (2002), different
transport codes can yield significantly different results at $N_{T}=1$, but
already at $N_{T}=32$ they were all in agreement regardless of implementation
details. We performed similar tests in our simulations, focusing on the effect
of different values of $N_{T}$ and $N_{\textrm{ens}}$ on flow observables; a
detailed description of our results can be found in Appendix D.2. Based on
these tests, we use $N_{T}=20$ and $N_{\textrm{ens}}=50$, which we find to be
large enough to yield simulation results that are independent of both the
specific density calculation scheme and specific collision criterion used.

At collision energies in the range
$E_{\textrm{lab}}\in[1,10]~{}\textrm{GeV/nucleon}$, proton production can be
substantially influenced by the light nuclei production mechanism employed in
the simulation. One can estimate the relative importance of the light nuclei
by comparing their yields at midrapidity to proton yields. While the role of
the light nuclei diminishes as one goes from $E_{\textrm{lab}}=2$ to
$8~{}\textrm{GeV/nucleon}$ (for example, at
$E_{\textrm{lab}}=2~{}\textrm{GeV/nucleon}$ the proton yield at mid-rapidity
is around 52 and the deuteron yield is around 22, while at
$E_{\textrm{lab}}=8~{}\textrm{GeV/nucleon}$ the proton yield remains
approximately the same and the deuteron yield reduces to around around 14 Witt
_et al._ (1999)), it is nevertheless a very substantial effect in the whole
energy range explored in our study.

The situation is different for the flow observables, for which the effects of
the light nuclei on the proton flow are smaller. However, they are still non-
negligible, see Mohs  _et al._ (2022) for a dedicated study. In this work, to
account for these effects, we employ a deuteron production model introduced in
Oliinychenko  _et al._ (2019), where deuterons are produced dynamically in a
chain of reactions $NN\leftrightarrow d^{\prime}$, $Nd^{\prime}\leftrightarrow
Nd$, effectively reducing to nucleon catalysis $NNN\leftrightarrow Nd$; or in
$NN\leftrightarrow d^{\prime}$, $\pi d^{\prime}\leftrightarrow\pi d$,
effectively reducing to pion catalysis $\pi NN\leftrightarrow\pi d$; or in the
reaction $NN\leftrightarrow\pi d$. The implementation and the cross sections
are a part of the publicly available SMASH 2.1 sma , and are described in
detail in Oliinychenko  _et al._ (2019). The most important deuteron
production reaction at the considered collision energies, due to the high
relative abundance of nucleons compared to pions, is nucleon catalysis. While
in our approach deuterons are dynamically produced and influence the entire
duration of the collisions, no light nuclei other than deuterons are produced
in our simulations. Both this fact and the uncertainty of the deuteron
production model contribute to the total systematic uncertainties in our
results due to light nuclei production (for other approaches to simulating
light nuclei, see, e.g., Le Fèvre  _et al._ (2019); Ikeno  _et al._ (2016)).

Above, we have highlighted some of the known unknowns, i.e., features of the
simulation that contribute to its systematic uncertainty: the implementation
of the density calculation (smearing), the collision scheme, and the light
nuclei production model. In addition, in Appendix D we further discuss
possible effects on the simulation results due to the nucleus initialization
model (Appendix D.3), Pauli blocking (Appendix D.4), the absence of isospin-
dependent, Coulomb, and momentum-dependent interactions (Appendix D.5), and
the absence of in-medium cross sections (Appendix D.6). Overall, some of these
effects are substantial and require further systematic study. A comparative
analysis taking into account all these possible effects is beyond the scope of
the current study. Within the effects explored in the present work, we find
the effects due to the EOS to be of the leading order.

##  III Sensitivity of the flow measurements to the EOS

Figure 3: Sensitivity of flow observables to the EOS in SMASH simulations of
Au+Au collisions at impact parameter $b=6$ fm. The top panels (a1)-(c1) show
the independent parametrizations of the EOS in three separate density regions;
panels (a2)-(c2) show $dv_{1}/dy^{\prime}$ for protons (full circles),
deuterons (open circles), $\pi^{+}$ (full squares), and $\pi^{-}$ (open
squares); panels (a3)-(c3) show $v_{2}(y^{\prime}=0)$ for protons and
deuterons; panels (a4)-(c4) show $v_{2}(y^{\prime}=0)$ for pions. Note that
the deuteron flow is divided by 2 for both $dv_{1}/dy^{\prime}$ and
$v_{2}(y^{\prime}=0)$, and the points for deuterons and pions are slightly
shifted to the right along the $E_{\textrm{lab}}$ axis for clarity. Results
shown within the columns (a), (b), (c) correspond to different regions in
which $c_{s}^{2}$ is varied, as can be seen in the top panels; within each
column, each EOS (top panels) and the corresponding simulation results (lower
panels) are shown using a separate color (purple, blue, green, orange, and red
for $c_{s}^{2}=0.1$, 0.3, 0.5, 0.7, and 0.9, respectively). Note that we show
$dv_{1}/dy^{\prime}$ instead of $dv_{1}/dy$, where $y^{\prime}$ is defined in
the center-of-mass frame as $y^{\prime}=y/y_{\textrm{beam}}$, so that
$y^{\prime}_{\textrm{beam}}=\pm 1$ in the center-of-mass frame. No $p_{T}$ cut
is imposed.

Our ultimate goal is to constrain the high-density EOS of symmetric nuclear
matter using flow measurements in heavy-ion collisions555We note that while Au
nuclei have a finite isospin asymmetry,
$\delta\equiv(N_{n}-N_{p})/(N_{n}+N_{p})\approx 0.198$, the asymmetric
contribution to the EOS is proportional to $\delta^{2}\approx 0.04$ and
therefore substantially suppressed. Moreover, while the leading contribution
to the symmetry energy is linear in density, the leading contribution to the
symmetric EOS is quadratic in density, which makes the symmetry energy
contribution less and less important as density increases. Consequently, not
only is extracting the EOS of asymmetric nuclear matter considerably
challenging at beam energies above $1~{}\textrm{GeV/nucleon}$, but also, to a
good approximation, at these energies one can neglect the symmetry energy
contribution and assume that the extracted EOS is that of symmetric nuclear
matter. We make such an assumption in this work.. In order to do this, in this
work we parametrize the EOS using the behavior of $c_{s}^{2}(T=0,n_{B})$ as a
function of the baryon density (as described in Sec. II.2), implement it into
the SMASH hadronic transport approach with a vector mean field
$A^{\mu}=\alpha(n_{B})j^{\mu}$ (with details of the implementation discussed
in Sec. II.3), simulate heavy-ion collisions, and analyze particle flow from
the obtained simulation data. In this section, we set out to investigate two
properties of flow observables as used for the extraction of the nuclear
matter EOS:

  * •

Sensitivity: How much is flow in heavy-ion collisions sensitive to changes of
the EOS in a specific density range?

  * •

Specificity: How much is flow in heavy-ion collisions sensitive to other
parameters of our simulations (beyond the EOS)?

An ideal observable would be sensitive to the EOS at a given density range and
insensitive to any other parameters. Additionally, a necessary condition to
extract the EOS from experimental data is that variations in simulation
results due to the EOS are larger than the experimental errors on the
measurements. Our results show (see, e.g., Fig. 3) that the latter condition
is easily satisfied for flow observables, as the change of the flow due to
varying the EOS is much larger than the experimental precision. For example,
recent STAR measurements of the proton flow at $\sqrt{s_{\mathrm{NN}}}=3.0$
and 4.5 GeV Adam  _et al._ (2021); Abdallah  _et al._ (2022) provide
$dv_{1}/dy^{\prime}$ with an absolute error around 0.01 and
$v_{2}(y^{\prime}=0)$ with an absolute error around 0.003.

Figure 4: Interpolated dependence of the proton flow in mid-central (with
impact parameter $b=6$ fm) Au+Au collisions on the speed of sound at different
densities for a series of collision energies
$E_{\textrm{lab}}=\\{2,~{}4,~{}6,~{}8\\}~{}\textrm{GeV/nucleon}$; see text for
more details.

To qualitatively test the sensitivity of flow observables to the EOS in
specific density ranges, we perform simulations using a parametrization of the
EOS with varying values of $c_{s}^{2}(n_{B},T=0)$ only in the density range of
interest, keeping a fixed dependence on $n_{B}$ everywhere else, and study how
much the flow results depend on the used EOS. In Fig. 3, we show proton,
deuteron, and pion $dv_{1}/dy^{\prime}$ and $v_{2}(y^{\prime}=0)$ obtained in
simulations where $c_{s}^{2}(n_{B},T=0)$ was varied only in the density range
$n_{B}\in(2,3)n_{0}$ (first column), only in the density range
$n_{B}\in(3,4)n_{0}$ (second column), and only in the density range
$n_{B}\in(4,5)n_{0}$ (third column); for $n_{B}<2n_{0}$,
$c_{s}^{2}(n_{B},T=0)$ took on the density-dependence coming from the
underlying default Skyrme EOS (see Sec. II.2 for more details), while in the
remaining, non-varied regions $c_{s}^{2}(n_{B},T=0)$ took on an arbitrary
constant value (see top panels in Fig. 3).

Figure 3 clearly shows that in Au+Au mid-central collisions the proton,
deuteron, and pion $dv_{1}/dy^{\prime}$ and $v_{2}(y^{\prime}=0)$ are very
sensitive to the EOS at $n_{B}\in(2,3)n_{0}$, with the maximal sensitivity to
this density region at the kinetic beam energy around
$E_{\textrm{lab}}=4~{}\textrm{GeV/nucleon}$. It is evident that a stiffer EOS
produces a larger $dv_{1}/dy^{\prime}$ and a smaller $v_{2}(y^{\prime}=0)$ for
protons, deuterons, and pions. One can see that the proton and deuteron flow
is also sensitive to the density region $n_{B}\in(3,4)n_{0}$, where the
maximum sensitivity is reached at the collision energy around
$E_{\textrm{lab}}=8~{}\textrm{GeV/nucleon}$. This sensitivity, however, is
smaller than sensitivity to the $n_{B}\in(2,3)n_{0}$ region. The flow of all
particles at any collision energy becomes rather insensitive to the EOS for
$n_{B}\in(4,5)n_{0}$, which limits the experimental opportunities to constrain
the EOS at $n_{B}>4n_{0}$ from heavy-ion collisions (we note, however, that
there is a possibility to use the deuteron flow at
$E_{\textrm{lab}}=16~{}\textrm{GeV/nucleon}$, where some sensitivity is still
present; see Fig. 3). While one could argue that if the EOS is very soft at
lower densities, then the fireball may spend more time at higher densities
and, consequently, flow observables might be sensitive to those higher
densities, in Fig. 3 we already take a soft EOS at low values of $n_{B}$ when
varying the EOS at $n_{B}\in(4,5)n_{0}$, and the sensitivity at the
$n_{B}\in(4,5)n_{0}$ region is still rather small.

Notice that the deuteron flow (for our deuteron production model see Sec.
II.3), which in Fig. 3 is divided by 2, can be measured with the same
experimental precision as the proton flow. Therefore, effectively, deuteron
flow is twice more sensitive to the EOS. Pion flow exhibits a moderate
dependence on the EOS in the $n_{B}\in(2,3)n_{0}$ region, but at higher
densities it becomes rather insensitive to the EOS at any collision energy.
Naturally, the most precise constraints can potentially be achieved if one
combines experimental data about proton, deuteron, and pion flow.

Let us now assess these sensitivities quantitatively. For this purpose, we
simulated mid-central (with impact parameter $b=6$ fm) Au+Au collisions at
collision energies
$E_{\textrm{lab}}=\\{2,~{}4,~{}6,~{}8\\}~{}\textrm{GeV/nucleon}$. The
simulations were run for 50 sets of values of $c^{2}_{[2,3]n_{0}}$ and
$c^{2}_{[3,4]n_{0}}$ (design points), sampled randomly in the
$(c^{2}_{[2,3]n_{0}},c^{2}_{[3,4]n_{0}})$ plane; $c^{2}_{[4,5]n_{0}}$ was not
varied and instead was set to equal 0.3, based on our observation (see Fig. 3)
that it does not appreciably influence the flow at the considered energies.
Next, a Gaussian emulator was used to interpolate between the design points
and obtain the flow dependence on the values of $c^{2}_{[2,3]n_{0}}$ and
$c^{2}_{[3,4]n_{0}}$, shown in Fig. 4. In this figure, we observe that at the
considered energies the dependence of $v_{2}(y^{\prime}=0)$ and especially of
$dv_{1}/dy^{\prime}$ on $c_{s}^{2}$ is very well described by a linear
approximation $a_{1}+a_{2}c^{2}_{[2,3]n_{0}}+a_{3}c^{2}_{[3,4]n_{0}}$, where
$a_{1},a_{2},a_{3}$ are regression coefficients specific to each energy, as
long as $c_{s}^{2}$ is positive. Based on this, we can summarize the
regression coefficients for protons in Au+Au collisions at $b=6$ fm as
follows:

$\displaystyle\begin{pmatrix}\frac{d{v_{1}}}{d{y^{\prime}}}|_{2\textrm{
GeV}}\\\ \frac{d{v_{1}}}{d{y^{\prime}}}|_{4\textrm{ GeV}}\\\
\frac{d{v_{1}}}{d{y^{\prime}}}|_{8\textrm{ GeV}}\\\ v_{2}|_{2\textrm{ GeV}}\\\
v_{2}|_{4\textrm{ GeV}}\\\ v_{2}|_{8\textrm{ GeV}}\end{pmatrix}\approx$
$\displaystyle\begin{pmatrix}0.40\\\ 0.04\\\ -0.04\\\ 0.00\\\ 0.04\\\
0.04\end{pmatrix}+\begin{pmatrix}0.35&0.00\\\ 0.50&0.14\\\ 0.36&0.24\\\
-0.06&0.00\\\ -0.06&-0.03\\\
-0.03&-0.02\end{pmatrix}\begin{pmatrix}c^{2}_{[2,3]n_{0}}\vspace{2mm}\\\
c^{2}_{[3,4]n_{0}}\end{pmatrix}\,.$ (21)

We reiterate that the above relations are rather accurate (with the
coefficient of determination $R^{2}>0.95$) for positive values of $c_{s}^{2}$.
Note that the numbers in the second matrix on the right-hand side are
effective measures of the sensitivities, i.e., they indicate how well a
measurement at a given collision energy constrains the EOS, assuming that
there are no other parameters influencing the flow except the EOS. For
example, measuring proton $dv_{1}/dy^{\prime}$ at
$E_{\textrm{lab}}=4~{}\textrm{GeV/nucleon}$ with an error bar of 0.01, one
obtains $c^{2}_{[2,3]n_{0}}$ with the lowest possible error of around
$0.01/0.50=0.02$, and $c^{2}_{[3,4n_{0}]}$ with the lowest possible error of
around $0.01/0.14\approx 0.07$. From the same matrix one can also see that
$dv_{1}/dy^{\prime}$ measurements are more constraining for the EOS, even
though $v_{2}(y^{\prime}=0)$ is measured with a better precision, because the
sensitivity of $v_{2}(y^{\prime}=0)$ is lower. Based on both Figs. 3 and 4 as
well as the sensitivity matrix in Eq. (21), one can see that the EOS for
$n_{B}\in(2,3)n_{0}$ is mainly constrained by measurements at lower energies,
around $E_{\textrm{lab}}=4~{}\textrm{GeV/nucleon}$, while the EOS for
$n_{B}\in(3,4)n_{0}$ is mainly constrained by measurements at higher energies,
around $E_{\textrm{lab}}=8~{}\textrm{GeV/nucleon}$.

However, the flow is influenced by parameters other than the EOS parameters.
We can evaluate the importance of varying these other parameters (often called
nuisance parameters) by finding an equivalent change in the EOS parameters, by
which we mean finding how big of an adjustment in, e.g., $c^{2}_{[2,3]n_{0}}$
is needed to compensate for the change in flow due to the variation of
nuisance parameters. For example, from Fig. 18 one can compute that, at
$E_{\textrm{lab}}=4~{}\textrm{GeV/nucleon}$, scaling all cross sections down
by a factor of 0.6 has an effect on $dv_{1}/dy^{\prime}$ that can be entirely
compensated by increasing $c^{2}_{[2,3]n_{0}}$ by around 0.2, as can be seen
from Eq. (21). Such considerations pave the way to a full Bayesian analysis,
where both the EOS parameters and the nuisance parameters are varied. However,
here we do not perform such an elaborate analysis, given that a computational
effort of such a scale should be motivated by a preliminary exploratory study
which we provide in this work. Thus we restrict ourselves to varying the EOS
parameters and only roughly estimate systematic errors due to the nuisance
parameters.

Figure 5: Impact parameter dependence of the proton flow in Au+Au collisions
as simulated in SMASH. At $n_{B}\in(0,2)n_{0}$, the mean-field is parametrized
to reproduce the default Skyrme EOS, while at higher densities the mean-field
is parametrized to yield $c^{2}_{s}[n_{B}\in(2,3)n_{0}]=0.1$ and
$c^{2}_{s}[n_{B}>3n_{0}]=0.3$. Notice that the decrease of
$|v_{2}|(y^{\prime}=0)$ in peripheral collisions is typical for potentials
without explicit momentum-dependent terms, such as the one we employ in this
work. With momentum-dependent terms, $|v_{2}|(y^{\prime}=0)$ exhibits a
monotonic growth against $b$ and attains large values in peripheral
collisions; see Fig. 4 of Danielewicz (2000). Which of these scenarios is
realized in nature remains to be tested, which will provide strong constraints
for the momentum-dependent potential terms.

##  IV Results

###  IV.1 Using both E895 and STAR data

Let us start the discussion of our results by reviewing the experimental data
that we want to fit to constrain the EOS. We use the flow data from the E895
Collaboration, measured at
$E_{\textrm{lab}}=\\{2,~{}4,~{}6,~{}8\\}~{}\textrm{GeV/nucleon}$
($\sqrt{s_{\textrm{NN}}}=\\{2.7,~{}3.3,~{}3.8,~{}4.3\\}~{}\textrm{GeV}$)
Pinkenburg  _et al._ (1999); Liu  _et al._ (2000), as well as the recent data
from the STAR Collaboration at
$\sqrt{s_{\textrm{NN}}}=\\{3,~{}4.5\\}~{}\textrm{GeV}$, spanning the collision
energy region where the flow is most sensitive to the variation of the EOS
above $n_{B}=2n_{0}$ (as already shown in Fig. 3). Let us briefly summarize
possible systematic differences between measurements from these two
experiments. The centrality selection in the E895 experiment was done by
considering events with charged particle multiplicities $M$ belonging to a
given range of fractions of $M_{\textrm{max}}$, where $M_{\textrm{max}}$ is
the value near the upper limit of the $M$ spectrum at which the height of the
distribution has fallen to half of its plateau value. The data presented by
the E895 Collaboration are for events with charged particle multiplicity
between $0.5M_{\textrm{max}}$ and $0.75M_{\textrm{max}}$, which by comparisons
to models was found to correspond to an impact parameter
$b=5$–$7~{}\textrm{fm}$ Pinkenburg  _et al._ (1999); Liu  _et al._ (2000).
This impact parameter estimation may be imprecise, but we can estimate the
error due to a shift in $b$ of $\pm 1~{}\textrm{fm}$ by using simulation
results presented in Fig. 5: for example, at
$E_{\textrm{lab}}=4~{}\textrm{GeV/nucleon}$, varying $b$ from 5 to 7 fm is
equivalent to a decrease in the proton $dv_{1}/dy^{\prime}$ of around 0.05,
which can be compensated by increasing $c^{2}_{[2,3]n_{0}}$ by around 0.1, as
can be seen from Eq. (21). The STAR centrality selection at
$\sqrt{s_{\textrm{NN}}}=3$ GeV is 10–40%, which is found Adamczewski-Musch
_et al._ (2018) to correspond to $b=4.7$–$9.3~{}\textrm{fm}$. At
$\sqrt{s_{\textrm{NN}}}=4.5~{}\textrm{GeV}$, the centrality selection for
calculating $v_{1}$ is 10–25%, corresponding to $b=4.7$–$7.4~{}\textrm{fm}$,
while for $v_{2}$ it is 0–30%, corresponding to $b=0.0$–$8.1~{}\textrm{fm}$.
These differences in centrality selection, along with differences in cuts on
the transverse momentum $p_{T}$, are summarized in Table 1.

Figure 6: Comparison of E895 Pinkenburg  _et al._ (1999); Liu  _et al._ (2000)
and STAR measurements Adam  _et al._ (2021); Abdallah  _et al._ (2022) for
$dv_{1}/dy^{\prime}$ (a) and $v_{2}(y^{\prime}=0)$ (b) of protons. While the
centrality selection and the used $p_{T}$ cuts are not identical (see Table
1), we find that correcting for these factors cannot explain the apparent
discrepancy in the $dv_{1}/dy^{\prime}$ measurements. Figure 7: Influence of
the lower $p_{T}$ cut on proton $dv_{1}/dy^{\prime}$ (a) and
$v_{2}(y^{\prime}=0)$ (b) at collision energy
$\sqrt{s_{\textrm{NN}}}=3~{}\textrm{GeV}$; we note that in the simulations, we
used an EOS parametrized with the maximum a posteriori probability parameters
$K_{0}$, $c^{2}_{[2,3]n_{0}}$, $c^{2}_{[3,4]n_{0}}$, which we will discuss in
the following. The shown dependency suggests that the difference between E895
and STAR measurements at $\sqrt{s_{\textrm{NN}}}\approx 3~{}\textrm{GeV}$
cannot be ascribed to different $p_{T}$ cuts applied by the experiments; see
Table 1.

Looking at the flow measurements shown in Fig. 6, one can see an apparent
disagreement between the experiments. It is valid to ask whether they can be
explained by a difference in the centrality or $p_{T}$ selections. We find
that they cannot, although to make this conclusion we have to rely on
simulations. Our analysis shows that changing the lower $p_{T}$ bound from 0.4
to 0.1 GeV (as used in the E895 experiment) would lower the STAR
$dv_{1}/dy^{\prime}$ at $\sqrt{s_{NN}}=3~{}\textrm{GeV}$, making it closer to
the E895 data. However, the magnitude of this change does not exceed 0.05 (see
Fig. 7), which is insufficient to explain the discrepancy. Decreasing the
upper bound on the centrality selection for STAR at
$\sqrt{s_{\textrm{NN}}}=3~{}\textrm{GeV}$ can only increase
$dv_{1}/dy^{\prime}$ (see Fig. 5), and therefore increase the discrepancy.
This makes us conclude that there is a substantial disagreement between the
E895 and the STAR flow data, even after taking into the account that they are
measured at somewhat different centrality and $p_{T}$ cuts.

Data |  $E_{\textrm{lab}}$ [GeV] |  $\sqrt{s_{\textrm{NN}}}$ [GeV] |  $b$ [fm] |  $p_{T}$-cut [GeV]
---|---|---|---|---
E895 | 2 | 2.7 | 5 – 7 | [0.1, 2.0]
E895 | 4 | 3.4 | 5 – 7 | [0.1, 2.0]
E895 | 6 | 3.8 | 5 – 7 | [0.2, 2.0]
E895 | 8 | 4.3 | 5 – 7 | [0.4, 2.0]
STAR | 2.9 | 3.0 | 4.7 – 9.3 | [0.4, 2.0]
STAR | 8.9 | 4.5 |  $v_{1}$: 4.7 – 7.4 | [0.4, 2.0]
|  |  |  $v_{2}$: 0.0 – 8.10 |
Table 1: Comparison of different experimental conditions for flow measurements
Pinkenburg  _et al._ (1999); Liu  _et al._ (2000); Adam  _et al._ (2021);
Abdallah  _et al._ (2022).

Figure 8: Priors (top) and posteriors (bottom) from analysis employing both
E895 and STAR data, where the posteriors were selected by restricting the
spread in simulation results for $dv_{1}/dy^{\prime}(E_{\textrm{lab}})$ and
$v_{2}(E_{\textrm{lab}})$ to encompass most of available data. The plot of
$c_{s}^{2}$ against density in the (a1) panel, showing all EOSs used in the
current analysis, has three characteristic regions: for $n_{B}<2n_{0}$, the
speed of sound squared follows the behavior given by the default Skyrme EOS;
for $n_{B}>4n_{0}$, $c_{s}^{2}$ assumes a constant value of 0.3; for
$2n_{0}<n_{B}<3n_{0}$ and $3n_{0}<n_{B}<4n_{0}$, thin horizontal grey lines
indicate values of $c_{s}^{2}$ in the corresponding region used in the current
analysis. The colorful lines joining the horizontal lines for
$2n_{0}<n_{B}<3n_{0}$ and $3n_{0}<n_{B}<4n_{0}$ indicate particular
combinations of the values of $c^{2}_{[2,3]n_{0}}$ and $c^{2}_{[3,4]n_{0}}$
used, with the colors continually changing depending on the values connected.
Lines of the same color on all panels correspond to the same combination of
$c^{2}_{[2,3]n_{0}}$ and $c^{2}_{[3,4]n_{0}}$; panels (a1) and (a2) serve as a
color legend. By looking at the thick posterior lines, panels (b2) and (c2),
one can see that it is not possible, within our model, to fit all experimental
data with a single EOS: $v_{2}$ needs a harder EOS and $dv_{1}/dy^{\prime}$
from E895 needs a softer one.

Given this disagreement in the experimental data, we decide to perform two
analyses. First, we try to find a range of the EOS parameters that allows one
to roughly encompass both the E895 and the STAR proton flow data. For this
purpose we simulate Au+Au collisions at $b=6$ fm and vary only
$c^{2}_{[2,3]n_{0}}$ and $c^{2}_{[3,4]n_{0}}$, keeping fixed the default
Skyrme parametrization at lower densities (see Sec. II.2 for more details).
Both varied parameters are allowed to take any values in the range
$c_{s}^{2}\in[-1,1]$.

We note here that because a negative speed of sound squared cannot occur in
equilibrated matter, regions of the phase diagram where it appears in
calculations are often “corrected” by performing the Maxwell construction. At
the same time, $c_{s}^{2}(n_{B},T=0)<0$ simply corresponds to a negative slope
of the pressure as a function of density, which occurs when the potential
$U(n_{B})$ becomes attractive within a certain density range. This in turn
indicates that the system is unstable and that, given sufficient time, a phase
separation would eventually take place; in fact, it has been shown that such a
phase separation indeed occurs in hadronic transport models employing the
corresponding mean-field interactions (for a recent study, see Sorensen and
Koch (2021) or Sorensen (2021)). Therefore let us stress, especially in the
context of hadronic transport simulations which are well suited for evolving
systems out of equilibrium, that there is nothing inherently wrong in an EOS
leading to a negative $c_{s}^{2}(n_{B},T=0)$, as it simply indicates that the
particle-particle interactions are on average attractive. Whether a phase
separation indeed occurs in that case within heavy-ion collision simulations
depends on multiple factors, including the average duration of the collision,
the characteristic time for spinodal decomposition, and the average
temperature reached in the collision region; nevertheless, regardless of a
phase transition occurring or not, the interactions affect the evolution of
the system and, consequently, the obtained values of $dv_{1}/dy^{\prime}$ and
$v_{2}(y^{\prime}=0)$.

We start our analysis by assuming prior distributions of the EOS parameters
$c^{2}_{[2,3]n_{0}}$ and $c^{2}_{[3,4]n_{0}}$ to be uniform, and use the
degree of agreement between the experimental data and simulations utilizing
sampled sets of parameters to infer posterior distributions, that is
conditional distributions of parameters given the measurements. In Fig. 8, we
show priors (top) and posteriors (bottom) from analysis employing both E895
and STAR data, where the posteriors were selected by restricting the spread in
simulation results for $dv_{1}/dy^{\prime}(E_{\textrm{lab}})$ and
$v_{2}(E_{\textrm{lab}})$ to encompass most of the available data. The plot of
$c_{s}^{2}$ against density, showing all EOSs used in the current analysis
(top left panel), has three characteristic regions: (i) for $n_{B}<2n_{0}$,
the speed of sound squared follows the behavior given by the default Skyrme
EOS for all used EOSs; (ii) for $n_{B}>4n_{0}$, $c_{s}^{2}$ assumes a constant
value of 0.3 for all used EOSs; (iii) for $n_{B}\in(2,3)n_{0}$ and
$n_{B}\in(3,4)n_{0}$, thin horizontal grey lines indicate values of
$c_{s}^{2}$ in the corresponding regions that were used in the current
analysis, having been sampled uniformly in the respective regions. The
colorful straight lines joining the horizontal lines for $n_{B}\in(2,3)n_{0}$
and $n_{B}\in(3,4)n_{0}$ indicate particular combinations of the values of
$c^{2}_{[2,3]n_{0}}$ and $c^{2}_{[3,4]n_{0}}$ used in a given EOS, with the
colors continually changing depending on the values connected; for example,
blue lines represent EOSs with considerably negative values of
$c^{2}_{[2,3]n_{0}}$ and large positive values of $c^{2}_{[3,4]n_{0}}$. The
same color coding is used in the plots showing flow results for all considered
EOSs (top middle and right panels). In the posterior plot for $c_{s}^{2}$
(bottom left panel) we emphasize two of the shown EOSs: An orange thick line
represents an EOS that is very stiff at $n_{B}\in(2,3)n_{0}$ and very soft at
$n_{B}\in(3,4)n_{0}$, and that can describe $v_{2}$ from both experiments as
well as $dv_{1}/dy^{\prime}$ from STAR, but not $dv_{1}/dy^{\prime}$ from E895
(see bottom middle and right panels). In contrast, a thick grey line
represents an EOS that is moderately soft both at $n_{B}\in(2,3)n_{0}$ and at
$n_{B}\in(3,4)n_{0}$, and that can describe the $dv_{1}/dy^{\prime}$ from
E895. (We also note here that based on Fig. 8, one can see that a negative
$c^{2}_{[2,3]n_{0}}$ tends to generate $dv_{1}/dy^{\prime}<0$, which is in
agreement with a known result from hydrodynamics that a phase transition in
the EOS can lead to a negative $dv_{1}/dy^{\prime}$ Stoecker (2005). As we
mentioned in the introduction, one can intuitively interpret this as the
fireball acting like an attractive spring on the spectators, causing them to
be “deflected inwards”.)

Our broad prior, shown in the upper left panel of Fig. 8, covers a large range
of proton $dv_{1}/dy^{\prime}$ and $v_{2}(y^{\prime}=0)$. Even though the
substantial disagreement between the E895 and the STAR data makes a precise
fit difficult, our posterior distribution restricts possible EOSs to ones with
$c^{2}_{[2,3]n_{0}}>0.1$. Nevertheless, the most important conclusion from
this analysis is that we cannot simultaneously describe the
$dv_{1}/dy^{\prime}$ and $v_{2}(y^{\prime}=0)$ measurements obtained by the
E895 experiment. The E895 $dv_{1}/dy^{\prime}$ data prefer softer potentials,
while the E895 $v_{2}$ data prefer harder ones. This result is in agreement
with the conclusion from the seminal work Danielewicz  _et al._ (2002), where
an attempt was made to describe the same data using a different hadronic
transport code (known as pBUU, Danielewicz and Bertsch (1991); Danielewicz
(2000)) and varying the stiffness of the EOS by changing the incompressibility
$K_{0}$.

Figure 9: Priors (red dashed lines) and posteriors (blue solid lines) obtained
using the STAR proton flow measurements.

###  IV.2 Using only STAR data

Given the difficulty to simultaneously describe the E895 $dv_{1}/dy^{\prime}$
and $v_{2}(y^{\prime}=0)$ results, as well as the discrepancy between the E895
and the STAR data, a more rigorous Bayesian analysis with more parameters is
unlikely to be helpful if one analyzes both the E895 and the STAR data
together. In that case, the final constraint will depend strongly on the
systematic errors of the E895 experiment, which are not shown in Pinkenburg
_et al._ (1999); Liu  _et al._ (2000). Therefore we decided to try discarding
the older E895 data from our fit and focus only on the recent STAR proton flow
data. The centrality selection and the $p_{T}$ cut in our simulations are the
same as those for the STAR experiment, listed in Table 1. This time, we vary
three parameters: the incompressibility $K_{0}$, which controls the curvature
of $\mathcal{E}/n_{B}$ (and, therefore, $U(n_{B})$) at $n_{B}=n_{0}$ (for the
other two constraints of the Skyrme EOS we take $n_{0}=0.166$ fm-3 and
$E_{\textrm{bin}}$ = -15.65 MeV), $c^{2}_{[2,3]n_{0}}$, and
$c^{2}_{[3,4]n_{0}}$. For $n_{B}>4n_{0}$, $c_{s}^{2}$ is fixed at 0.3. To
infer the probability distributions of the parameter values, we use the
JETSCAPE statistical framework for Bayesian analysis introduced and applied in
Everett  _et al._ (2021).

Figure 10: Posterior parameter distribution obtained using the STAR proton
flow measurements. The maximum a posteriori probability (MAP) parameters are
$K_{0}=285\pm 67~{}\textrm{MeV}$, $c^{2}_{[2,3]n_{0}}=0.49\pm 0.13$,
$c^{2}_{[3,4]n_{0}}=-0.03\pm 0.15$.

In Fig. 9 one can see that, in contrast to the case including the E895
experiment data, our model can fit both the $dv_{1}/dy^{\prime}$ and the
$v_{2}(y^{\prime}=0)$ measurements from the STAR experiment. The posterior
distribution of the parameters is shown in Fig. 10. The incompressibility
$K_{0}$ is not very well constrained by the data, which was expected as the
STAR collision energies (and similarly the E895 collision energies) mainly
probe densities $n_{B}>2n_{0}$, while a much stronger constraint on $K_{0}$
can be obtained from lower energy experiments Le Fèvre  _et al._ (2016) (an
extensive discussion of current constraints on $K_{0}$ can be found in Dutra
_et al._ (2012); however, see also Stone  _et al._ (2014)). Nevertheless, we
can see in Fig. 10 that there is a clear anticorrelation between $K_{0}$ and
$c^{2}_{[2,3]n_{0}}$, i.e., that a stiffer EOS at lower densities may be
compensated by a somewhat softer EOS at higher densities. Moreover,
$c^{2}_{[2,3]n_{0}}$ is rather strongly constrained, mainly by the STAR
$\sqrt{s_{\textrm{NN}}}=3~{}\textrm{GeV}$ data, while $c^{2}_{[3,4]n_{0}}$ is
mainly constrained by the STAR $\sqrt{s_{\textrm{NN}}}=4.5~{}\textrm{GeV}$
data (see Fig. 3 and the accompanying description in the text).

The maximum a posteriori probability (MAP) parameters are $K_{0}=285\pm
67~{}\textrm{MeV}$, $c^{2}_{[2,3]n_{0}}=0.49\pm 0.13$,
$c^{2}_{[3,4]n_{0}}=-0.03\pm 0.15$. Given this result, we conclude that the
recent STAR proton flow measurements indicate a very hard potential at
$n_{B}\in(2,3)n_{0}$, followed by a substantial softening at
$n_{B}\in(3,4)n_{0}$. Indeed, in Fig. 11, showing a scatter plot of the
pressure as a function of baryon density for 50000 EOSs obtained by sampling
$K_{0}$, $c^{2}_{[2,3]n_{0}}$, and $c^{2}_{[3,4]n_{0}}$ from the posterior
distribution and putting $c^{2}_{[4,+\infty]n_{0}}=0.3$ (we stress here that
it was not possible to constrain the values of $c_{s}^{2}$ for $n_{B}>4n_{0}$
within our study, see Sec. III for more details), one can see that the
pressure prefers a substantial softening for $n_{B}\in(3,4)n_{0}$.

Figure 11: Scatter plot of the pressure as a function of baryon density,
obtained with 50000 EOSs with $K_{0}$, $c^{2}_{[2,3]n_{0}}$, and
$c^{2}_{[3,4]n_{0}}$ sampled from the posterior distribution. The solid line
corresponds to the mean of the pressures, while the thick and thin dashed line
corresponds to the $\pm 1\sigma$ and $\pm 3\sigma$ contour around the mean,
respectively.

A negative $c^{2}_{[3,4]n_{0}}$ (or, equivalently, a negative slope of
pressure for $n_{B}\in(3,4)n_{0}$) is a sign that the EOS exhibits a first-
order phase transition in that region. Using the thermodynamics of the model
(see Appendix C), one can identify the critical temperature for the EOS with
the MAP parameters at $T_{c}\approx 70\ \textrm{MeV}$. At the same time, an
EOS with $c^{2}_{[3,4]n_{0}}=-0.18$, which is a value smaller by 1$\sigma$
than the central value, corresponds to a critical temperature of $T_{c}\approx
260\ \textrm{MeV}$, while an EOS with the value of $c^{2}_{[3,4]n_{0}}$ larger
by 1$\sigma$ than the central value, $c^{2}_{[3,4]n_{0}}=0.12$, does not lead
to a first-order phase transition at all, although it still displays a
significant softening. This is, of course, a very large spread of possible
thermodynamic properties of dense nuclear matter in this region. Drawing firm
conclusions on the EOS using our study is further complicated by the lack of
momentum-dependent interactions in our model. In particular, since at higher
energies momentum-dependent interactions are repulsive, describing the
experimental results using a model without momentum-dependence leads to a
spuriously stiff EOS. If this effect is substantial for $n_{B}\in(2,3)n_{0}$,
then we can see from Fig. 10 that a very stiff EOS at $n_{B}\in(2,3)n_{0}$ can
induce a very soft EOS at $n_{B}\in(3,4)n_{0}$. Since it is unclear to what
extent the stiffness of the EOS for densities $n_{B}\in(2,3)n_{0}$ is caused
by the the lack of momentum-dependent interactions, we remain cautious about
making strong conclusions on the EOS.

Figure 12: Top panel: The speed of sound $c_{s}^{2}$ calculated for the
central values of the MAP parameters for both exactly symmetric nuclear matter
(solid black line) and pure neutron matter, the latter obtained using three
values of the slope parameter of the symmetry energy $L$ (yellow long-dashed,
green dash-dotted, and blue short-dashed lines); a smaller (bigger) $L$
results in a smaller (bigger) $c_{s}^{2}$ in pure neutron matter. Note that,
for all values of $L$ used, the transformation to pure neutron matter results
in a disappearance of the first-order phase transition. Bottom panel: The
speed of sound $c_{s}^{2}$ for exactly symmetric nuclear matter at three
values of $c^{2}_{[3,4]n_{0}}=\\{0.12,-0.03,-0.18\\}$ (thick black long-
dashed, medium solid, and thin short-dashed line, respectively), corresponding
to the central MAP value and values at $\pm 1\sigma$, and $c_{s}^{2}$ in pure
neutron matter obtained using a high value of $L$ for SYM1 (ASYM1, thick
yellow long-dashed line) and a low value of $L$ for SYM3 (ASYM3, thin short-
dashed blue line). The spread between the ASYM1 and ASYM3 curves illustrates
the uncertainty in $c_{s}^{2}$ in pure neutron matter at $n_{B}\in[3,4]n_{0}$
given both the uncertainty in our results and the uncertainty in the value of
$L$.

Nevertheless, taking our results at their face value, we see that they imply a
substantial softening of the EOS at large baryon densities. One may ask
whether this result is consistent with the current knowledge of the EOS based
on neutron star data Bedaque and Steiner (2015); Tews  _et al._ (2018);
Fujimoto  _et al._ (2020); Marczenko  _et al._ (2023), which indicates a very
stiff EOS at moderate densities, including a local maximum in the speed of
sound at which $c_{s}^{2}>1/3$, and significantly disfavors large softening of
the EOS that inevitably arises when a first-order phase transition is present.
Importantly, the differences between symmetric (heavy-ion collisions) and
asymmetric nuclear matter (neutron stars) should not be overlooked. To test
the influence of the symmetry energy on the behavior of $c_{s}^{2}$ as a
function of density, we take a simple expansion of the symmetry energy around
$n=n_{0}$ (see Appendix E for more details) and use it to transform the speed
of sound squared extracted from our analysis, calculated for exactly symmetric
nuclear matter, to the speed of sound squared in pure neutron matter. Given
the uncertainty in estimates for the value of the symmetry energy slope
parameter, we test three values of $L$ corresponding to the central and
central $\pm 1\sigma$ values, which, taking $L=58.7\pm 28.1~{}\textrm{MeV}$
Oertel  _et al._ (2017) (see also Li and Han (2013)), means that we use
$L=\\{30.6,~{}58.7,~{}86.6\\}~{}\textrm{MeV}$. In the top panel of Fig. 12,
the $c_{s}^{2}$ in symmetric nuclear matter, calculated using the central
values of the MAP parameters (that is for $K_{0}=285~{}\textrm{MeV}$,
$c^{2}_{[2,3]n_{0}}=0.49$, $c^{2}_{[3,4]n_{0}}=-0.03$, and setting
$c^{2}_{[4,\infty]n_{0}}=0.3$), is shown with the black solid line. With
yellow long-dashed, green dash-dotted, and short-dashed blue lines we show the
corresponding results for pure neutron matter, using the three values of the
slope parameter of the symmetry energy $L$. It is evident that a smaller
(bigger) $L$ results in a smaller (bigger) $c_{s}^{2}$ in pure neutron matter.
Note that for all values of $L$ used, the transformation to pure neutron
matter results in a disappearance of the first-order phase transition. In the
bottom panel, we show $c_{s}^{2}$ for exactly symmetric nuclear matter at
three values of $c^{2}_{[3,4]n_{0}}=\\{0.12,-0.03,-0.18\\}$ (thick black long-
dashed, medium solid, and thin short-dashed line, respectively), corresponding
to the central MAP value and boundary values within 1$\sigma$. For the pure
neutron matter, we show two curves of $c_{s}^{2}$, obtained using a high value
of $L$ for SYM1 (ASYM1, thick long-dashed yellow line) and a low value of $L$
for SYM3 (ASYM3, thin short-dashed blue line). The spread between the ASYM1
and ASYM3 curves illustrates the uncertainty in the speed of sound in pure
neutron matter at $n_{B}\in[3,4]n_{0}$ given both the uncertainty in our
results and the uncertainty in the value of the slope parameter $L$. We note
here that this spread might be even larger given the large values of the
symmetry energy slope parameter reported by the PREXII experiment, $L=106\pm
37~{}\textrm{MeV}$ Adhikari  _et al._ (2021). We also point out that recently,
an extensive study was performed in which the influence of the symmetry energy
expansion parameters on the conversion between neutron matter EOS and
symmetric matter EOS was studied in detail Yao  _et al._ .

Figure 13: Pressure in pure neutron matter as a function of baryon density for
50 EOSs sampled from the posterior distribution, calculated using two limiting
values of L: $L=30.6~{}\textrm{MeV}$ (blue solid lines) and
$L=86.8~{}\textrm{MeV}$ (orange solid lines). Also shown are the previous
constraints based on heavy-ion data from Danielewicz  _et al._ (2002),
likewise calculated for a soft (grey shaded region) and a stiff (yellow shaded
region) symmetry energy, as well as a constraint based on the recent NICER
analysis of the J0740 pulsar measurement Miller  _et al._ (2021) (area with
grey backward stripes) and a constraint based on a global analysis of
available neutron star data Essick  _et al._ (2020) (area with green forward
stripes).

Based on the above discussion, we see that in general a lack of a first-order
phase transition in pure neutron matter does not exclude a first-order phase
transition in exactly symmetric nuclear matter, and in particular we conclude
that while there is some tension between our results and the neutron star
data, the discrepancy is not significant. Given both the obtained statistical
uncertainty of our results and the fact that the study we present is minimal,
at this point we do not claim that there indeed is a first-order phase
transition in exactly symmetric nuclear matter around $n_{B}\in[3,4]n_{0}$,
but rather, we point out that such a possibility is consistent with the range
if possible EOSs indicated by our study.

Figure 14: Pressure in pure neutron matter as a function of energy density for
50 EOSs sampled from the posterior distribution, calculated using two limiting
values of L: $L=30.6~{}\textrm{MeV}$ (blue solid lines) and
$L=86.8~{}\textrm{MeV}$ (orange solid lines). Also shown are constraints
obtained in two analyses using an interpolation of the EOS at intermediate
densities together with perturbative QCD constraints: Annala  _et al._ (2020)
(area with green forward stripes) and the more recent Gorda  _et al._ (2023)
(maroon shaded area).

With the symmetry energy expansion, we can explore not only the behavior of
the speed of sound, but also of the pressure in pure neutron matter. In Fig.
13, we show pressure in pure neutron matter as a function of baryon density
for 50 EOSs sampled from the posterior distribution, calculated using two
limiting values of the slope parameter, $L=30.6~{}\textrm{MeV}$ (blue solid
lines) and $L=86.8~{}\textrm{MeV}$ (orange solid lines). Also shown are some
of the constraints on the pressure from other works: the seminal constraint
from Danielewicz  _et al._ (2002) for soft (grey shaded area) and stiff
(yellow shaded area) asymmetry energy, a constraint based on an analysis of
the NICER measurements of the J0740 pulsar Miller  _et al._ (2021) (area with
grey backward stripes), and a constraint from a theoretical analysis utilizing
chiral effective field theory coupled with quantum Monte Carlo methods and
constrained by observational data on massive pulsars, gravitational waves, and
the recent NICER measurement Essick  _et al._ (2020) (area with green forward
stripes). It is evident that compared to Danielewicz  _et al._ (2002), our
results favor a much stiffer EOS at lower densities and a somewhat softer EOS
at higher densities. The idea that such an EOS would be successful in
describing the data was already suggested, although not explored, in
Danielewicz  _et al._ (2002). In particular, we see that the posterior EOSs
obtained using a stiff slope parameter are largely consistent with the EOS as
inferred by Essick  _et al._ (2020).

In Fig. 14, we show pressure in pure neutron matter as a function of energy
density, where the curves showing the constraint obtained in this work have
been obtained in the same fashion as for Fig. 13. Additionally, we show limits
obtained from two analyses using an interpolation of the EOS at intermediate
densities together with perturbative QCD constraints, Annala  _et al._ (2020)
(area with green forward stripes) and the more recent Gorda  _et al._ (2023)
(maroon shaded area). Interestingly, here the EOSs from our softer set are
more consistent with the tighter constraint of Gorda  _et al._ (2023) for
energy densities $\mathcal{E}\lesssim
300~{}\textrm{MeV}\cdot\textrm{fm}^{-3}$, or equivalently for densities
$n_{B}\lesssim 2n_{0}$, but conversely, the EOSs from the tighter set are more
consistent with that constraint at $\mathcal{E}\gtrsim
500~{}\textrm{MeV}\cdot\textrm{fm}^{-3}$, or equivalently for $n_{B}\gtrsim
3n_{0}$. This suggests an interesting possibility that the symmetry energy
could be softer at moderate densities and harder at large densities.

Finally, let us stress that our results for $n_{B}\geq 4n_{0}$ in Figs. 11,
12, 13, and 14 originate from using $c^{2}_{[4,\infty]n_{0}}=0.3$, a value
that we set ad hoc due to the lack of a constraint from heavy-ion collisions.

##  V Summary and discussion

In this work, we performed a Bayesian analysis of the flow data from heavy-ion
collision experiments using hadronic transport simulations. Within our
framework, the mean-field potentials can be freely parametrized by the
density-dependence of the speed of sound squared at $T=0$ and the
incompressibility of nuclear matter, $K_{0}$. To constrain the EOS, we choose
a piecewise parametrization of the speed sound $c_{s}^{2}$ in which
$c_{s}^{2}[n_{B}<2n_{0}]$ is that of a Skyrme model (with stiffness controlled
by varied values of $K_{0}$), while $c_{s}^{2}[n_{B}\in(2,3)n_{0}]$ and
$c_{s}^{2}[n_{B}\in(3,4)n_{0}]$ assume constant values sampled from the
interval $[-1,1]$; for $n_{B}>4n_{0}$, we assume $c_{s}^{2}[n_{B}>4n_{0}]=0.3$
(see Sec. II.2 for details). After assessing the sensitivity of the flow
measurements to the EOS (Sec. III), we put a constraint on the dense nuclear
matter EOS by performing a Bayesian analysis of simulation results and
experimental data (Sec. IV). Our study leads us to make the following
conclusions:

  * •

Flow observables at $\sqrt{s_{\textrm{NN}}}=2.5$–$5~{}\textrm{GeV}$ are very
sensitive to the dense nuclear matter EOS at $n_{B}\in(2,4)n_{0}$. While lower
densities can be studied by means of collisions at lower energies, there is
almost no possibility to constrain the EOS at baryon densities $n_{B}>4n_{0}$
from AA collisions, at least based on the analysis of flow observables. A
similar conclusion was obtained in a concurrent UrQMD study Steinheimer  _et
al._ (2022).

  * •

In particular, we find that the proton flow can yield a very tight constraint
on the EOS, and an even better constraint can be obtained from the deuteron
flow for which the sensitivity to the EOS is twice as large as that for
protons. To a lesser extent, pion flow can also be used to help constrain the
EOS.

  * •

Even given a large freedom to vary the EOS differentially in different density
regions, we cannot describe the E895 proton flow data: $dv_{1}/dy^{\prime}$
prefers a relatively soft EOS while $v_{2}(y^{\prime}=0)$ prefers a harder
one. However, the more recent STAR data, which additionally seems to be in
disagreement with the E895 data, can be described within our model, yielding
$K_{0}=285\pm 67~{}\textrm{MeV}$, $c^{2}_{[2,3]n_{0}}=0.49\pm 0.13$,
$c^{2}_{[3,4]n_{0}}=-0.03\pm 0.15$ and thus indicating a very hard EOS at
$n_{B}\in(2,3)n_{0}$ and a possible phase transition at $n_{B}\in(3,4)n_{0}$.

  * •

While the obtained set of EOS parameters results in simulations that reproduce
experimental measurements, the extraction of the EOS performed in this work is
not definitive given the substantial systematic uncertainties of the model
originating from the absence of momentum-dependent potentials. Additional
contributions to the uncertainty may come from the choice of the deuteron
production model, absence of light nuclei other than deuterons in the model,
absence of isospin-dependent and Coulomb interactions, and other sources,
including technical details of the simulation.

The possibility of a very stiff EOS at relatively low densities above the
nuclear saturation density, which our study suggests, is supported by analyses
of neutron star EOSs Bedaque and Steiner (2015); Tews  _et al._ (2018);
Fujimoto  _et al._ (2020); Marczenko  _et al._ (2023), where the strongest
experimental astrophysical constraint comes from the existence of neutron
stars with masses above two solar masses Demorest  _et al._ (2010); Antoniadis
_et al._ (2013); Riley  _et al._ (2021); Fonseca  _et al._ (2021). On the
other hand, our finding of a relatively soft EOS at moderate densities,
$n_{B}\in(3,4)n_{0}$, seems to create some tension with neutron star studies.
However, we stress that we only constrain a symmetric nuclear matter EOS. To
obtain the neutron star EOS, knowledge of the isospin-dependence of the EOS is
necessary; while to some extent one can use the symmetry energy expansion to
perform this transformation, at this time the coefficients of the expansion
still carry significant uncertainties, and additionally the expansion might
only be valid at moderate densities, leaving the behavior of the symmetry
energy at high baryon densities largely unknown. Given this uncertainty and
the statistical uncertainty from our Bayesian analysis (see Fig. 12) as well
as the systematic uncertainty of our model, we see the tension between our
results and the neutron star data as not significant. However, it may become
meaningful if the mentioned uncertainties decrease.

Our statistical uncertainties can be improved by certain future measurements,
assuming that the model will be able to describe them simultaneously. In
particular, flow systematics measured at several energies within
$\sqrt{s_{\textrm{NN}}}=2.5$–$5~{}\textrm{GeV}$, that is, the centrality, $y$,
and $p_{T}$ dependences of flow measured for different hadron species (proton,
lambda, deuteron, pion, etc.), will decrease statistical errors of the fit,
and also allow one to make the EOS parametrization more differential. We also
expect that there is a potential for constraining the EOS using measurements
of Hanbury-Brown-Twiss (HBT) correlations and dilepton spectra. Describing
these measurements will also require improvements of the model, leading to
reduced systematic errors.

The systematic errors of our simulations can only be alleviated by
improvements of the approach itself. In Sec. II.3 and Appendix D, we discuss
in detail the known weak points of our model and the possible influence of its
technical and theoretical features on the results. In particular, while we
found that varying the EOS has a large effect on the simulation results in the
studied energy region, other details of the simulation, not constrained within
this study, also lead to non-negligible effects, implying substantial
systematic errors. Most importantly, momentum-dependent potentials are missing
from our approach. While their inclusion is crucial for a differential
analysis of flow measurements (e.g., the $p_{T}$ dependence of $v_{2}$), their
absence in our model introduces a substantial systematic error for integrated
flow as well. The next important sources of the systematic error are the
deuteron production model as well as not including light nuclei other than
deuterons. Furthermore, Coulomb interactions and isospin-dependent potentials
were not present in the used model. While their effects seem to be rather
small at the studied energies, inclusion of these potentials is required for
an accurate assessment of the magnitude of the effect. They would also be
required for simulations at lower energies constraining the EOS at densities
$n_{B}<2n_{0}$, including meaningfully constraining $K_{0}$. Further, some of
the incurred systematic errors could be estimated by performing comparisons
between different transport simulations, where we note that for the error
estimate to be meaningful, it would be important to establish a set of minimum
requirements that each such simulation should satisfy, such as sufficiently
accurate energy and momentum conservation, reasonable relativistic properties,
and reproducing measured particle yields or spectra. Many such comparisons are
done within the Transport Model Evaluation Project; see Wolter  _et al._
(2022) for a review and further references. An overview of multiple potential
improvements and research directions in transport modeling was recently
provided in Sorensen  _et al._ (2023).

In summary, while future improvements on this work are possible both from the
theoretical and the experimental side, our current results prove that flow
measurements are very sensitive to the underlying EOS at densities
$n_{B}\in(2,4)n_{0}$ and that, with a better control of the simulations and of
the EOS as well as a more expansive Bayesian analysis, they can potentially
put a very tight constraint on the EOS of nearly symmetric nuclear matter.
Supported by model developments and together with studies of the EOS in
neutron stars as well as future tighter constraints on the symmetry energy,
heavy-ion data have the power to constrain the EOS of nuclear matter within a
sizable region of the QCD phase diagram.

###### Acknowledgements.

Some of results described in this work were presented at the 15th Workshop on
Particle Correlations and Femtoscopy at Michigan State University, and we wish
to thank the organizers for creating a great atmosphere for scientific
discussions. D.O. and A.S. thank Paweł Danielewicz, Roy Lacey, Bill Lynch,
Betty Tsang, Scott Pratt, Oleh Savchuk, Chun Shen, Sangwook Ryu, Volodymyr
Vovchenko, Jan Steinheimer, Nanxi Yao, Jaki Noronha-Hostler, and Veronica
Dexheimer for discussions of (selected) results and insightful comments. We
also thank the HEPData project, which allowed us to obtain experimental data
in a convenient way. D.O. thanks Rongrong Ma, Shaowei Lan, Hanna Zbroszczyk,
and Declan Keane for comments on experimental data. D.O., A.S., and L.M.
acknowledge support by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, under Grant No. DE-FG02-00ER41132. V.K. was
supported by the U.S. Department of Energy, Office of Science, Office of
Nuclear Physics, under Contract No. DE-AC02-05CH11231.

##  Appendix A The testparticle ansatz

In this appendix we obtain the test-particle equations of motion, Eqs. 7-8, by
employing the standard test particle ansatz. We base this derivation on a more
general case for off-shell particles which can be found in Leupold (1999).

As already seen in Eq. (5), the eight-dimensional distribution function
$f(x,p)$ can be written as

$\displaystyle
f(x^{\mu},p^{\mu})=\kappa(\Pi^{\mu})\tilde{f}(x^{\mu},\bm{p})~{},$ (22)

where $\kappa(\Pi^{\mu})$ is a factor ensuring that the momenta in the system
are on shell,

$\displaystyle\kappa(\Pi^{\mu})=2(2\pi)\Theta(\Pi^{0})\delta\big{(}\Pi^{\mu}\Pi_{\mu}-m^{2}\big{)}~{},$
(23)

and the continuous distribution $\tilde{f}$ is approximated by the test-
particle ansatz, Eq. (6), repeated here for convenience with the suppressed
index $i$,

$\displaystyle\tilde{f}(x,\bm{p})=\frac{1}{N_{T}}\sum_{j=1}^{N_{T}N}\delta^{(3)}\big{(}\bm{x}-\bm{x}_{j}(t)\big{)}\delta^{(3)}\big{(}\bm{p}-\bm{p}_{j}(t)\big{)}~{}.$
(24)

Substituting Eq. (22) into the Vlasov equation (that is, the Boltzmann
equation, Eq. (3), with the vanishing collision term, $I_{\textrm{coll}}=0$),
results in

$\displaystyle\kappa\Big{[}\Pi_{\mu}\partial^{\mu}_{x}\tilde{f}+\Pi_{\nu}\big{(}\partial_{\mu}^{x}A^{\nu}\big{)}\partial_{p}^{\mu}\tilde{f}\Big{]}$
$\displaystyle\hskip
28.45274pt+~{}\tilde{f}\Big{[}\Pi_{\mu}\partial^{\mu}_{x}\kappa+\Pi_{\nu}\big{(}\partial_{\mu}^{x}A^{\nu}\big{)}\partial_{p}^{\mu}\kappa\Big{]}=0~{}.$
(25)

We rewrite the expression in the first square bracket by explicitly separating
the zeroth components of the four-vectors,

$\displaystyle\Pi_{0}\partial^{0}_{x}\tilde{f}+\Pi_{\nu}\big{(}\partial_{0}^{x}A^{\nu}\big{)}\partial_{p}^{0}\tilde{f}$
$\displaystyle\hskip
28.45274pt+~{}\Pi_{i}\partial^{i}_{x}\tilde{f}+\Pi_{\nu}\big{(}\partial_{i}^{x}A^{\nu}\big{)}\partial_{p}^{i}\tilde{f}~{}.$
(26)

From Eq. (24) it is clear that $\partial{\tilde{f}}/\partial{p^{0}}=0$ and
$d{\tilde{f}}/d{t}=0$, with the latter allowing one to write
$\partial{\tilde{f}}/\partial{t}$ as

$\displaystyle\frac{\partial{\tilde{f}}}{\partial{t}}$
$\displaystyle=-\left(\frac{d{x^{k}}}{d{t}}\frac{\partial{}}{\partial{x^{k}}}+\frac{d{p^{k}}}{d{t}}\frac{\partial{}}{\partial{p^{k}}}\right)\tilde{f}$
$\displaystyle=\frac{1}{N_{T}}\sum_{j=1}^{N_{T}N}\left(-\frac{d{x^{k}}}{d{t}}\partial_{k}^{x}-\frac{d{p^{k}}}{d{t}}\partial_{k}^{p}\right)$
$\displaystyle\hskip
28.45274pt\times~{}\delta^{(3)}(\bm{x}-\bm{x}_{j})\delta^{(3)}(\bm{p}-\bm{p}_{j})~{}.$
(27)

Altogether, the Vlasov equation becomes

$\displaystyle\kappa\frac{1}{N_{T}}\sum_{j=1}^{N_{T}N}\bigg{[}\bigg{(}\Pi^{i}-\Pi_{0}\frac{d{x^{i}}}{d{t}}\bigg{)}\partial_{i}^{x}$
$\displaystyle\hskip
28.45274pt+~{}\bigg{(}-\Pi_{0}\frac{d{p^{i}}}{d{t}}+\Pi_{\nu}\big{(}\partial^{i}_{x}A^{\nu}\big{)}\bigg{)}\partial^{p}_{i}\bigg{]}$
$\displaystyle\hskip
56.9055pt\times~{}\delta^{(3)}(\bm{x}-\bm{x}_{j})\delta^{(3)}(\bm{p}-\bm{p}_{j})$
$\displaystyle\hskip
28.45274pt+~{}\tilde{f}\Big{[}\Pi_{\mu}\partial^{\mu}_{x}\kappa+\Pi_{\nu}\big{(}\partial_{\mu}^{x}A^{\nu}\big{)}\partial_{p}^{\mu}\kappa\Big{]}=0~{}.$
(28)

The above equation will be always satisfied if, for every test particle $j$,
the following sufficient conditions are met:

$\displaystyle\frac{d{(x_{j})^{i}}}{d{t}}=\frac{(\Pi_{j})^{i}}{(\Pi_{j})^{0}}~{},$
(29)
$\displaystyle\frac{d{(p_{j})^{i}}}{d{t}}=\frac{(\Pi_{j})^{\nu}}{(\Pi_{j})^{0}}\big{(}\partial^{i}_{x}A_{\nu}\big{)}~{},$
(30) $\displaystyle(\Pi_{j})^{\mu}(\Pi_{j})_{\mu}=m^{2}~{},$ (31)

where the last condition ensures that $\kappa=\textrm{const}$ and, therefore,
$\tilde{f}\big{[}\Pi_{\mu}\partial^{\mu}_{x}\kappa+\Pi_{\nu}\big{(}\partial_{\mu}^{x}A^{\nu}\big{)}\partial_{p}^{\mu}\kappa\big{]}=0$.

To bring Eqs. (29) and (30) into a fully relativistic form, we first note that
we can always write (here and in the following we have dropped the index $j$
for clarity)

$\displaystyle\frac{d{x^{0}}}{d{t}}=1=\frac{\Pi_{0}}{\Pi_{0}}~{},$ (32)

so that we can immediately write

$\displaystyle\frac{d{x^{\mu}}}{d{t}}=\frac{\Pi^{\mu}}{\Pi^{0}}~{},$ (33)

obtaining the same form as given in Eq. (7). Similarly, using
$p^{0}=\Pi^{0}+A^{0}=\sqrt{\big{(}\bm{p}-\bm{A}\big{)}^{2}+m^{2}}+A^{0}$, we
can write

$\displaystyle\frac{d{p^{0}}}{d{x_{0}}}$
$\displaystyle=\frac{\Pi^{i}}{\Pi_{0}}\frac{\partial{A_{i}}}{\partial{x_{0}}}+\frac{\partial{A^{0}}}{\partial{x_{0}}}=\frac{\Pi^{\nu}}{\Pi_{0}}\frac{\partial{A_{\nu}}}{\partial{x_{0}}}~{},$
(34)

which can be combined with Eq. (30) to yield

$\frac{d{p^{\mu}}}{d{t}}=\frac{\Pi_{\nu}}{\Pi^{0}}\big{(}\partial^{\mu}_{x}A^{\nu}\big{)}~{}.$
(35)

By substituting $p^{\mu}=\Pi^{\mu}+A^{\mu}$ and using

$\frac{d{A^{\mu}}}{d{t}}=\frac{d{x_{\nu}}}{d{t}}\frac{\partial{A^{\mu}}}{\partial{x_{\nu}}}=\frac{\Pi_{\nu}}{\Pi_{0}}\big{(}\partial^{\nu}_{x}A^{\mu}\big{)}~{},$
(36)

we obtain

$\frac{d{\Pi^{\mu}}}{d{t}}=\frac{\Pi_{\nu}}{\Pi^{0}}\Big{(}\partial^{\mu}_{x}A^{\nu}-\partial^{\nu}_{x}A^{\mu}\Big{)}~{}.$
(37)

Finally, defining

$F^{\mu\nu}\equiv\partial^{\mu}_{x}A^{\nu}-\partial^{\nu}_{x}A^{\mu}~{},$ (38)

we obtain

$\frac{d{\Pi^{\mu}}}{d{t}}=\frac{\Pi_{\nu}}{\Pi^{0}}F^{\mu\nu}~{},$ (39)

which agrees with Eq. (8).

It is worth noting that Eqs. (7) and (8) are analogous to the relativistic
equations of motion of a charge in an electromagnetic field, where the crucial
difference in our case is that the role of the field is played by the self-
consistently calculated baryon four-current instead of an external field.
Indeed, it is possible to rewrite Eq. (8) in the form

$\displaystyle\frac{d\bm{\Pi}}{dt}$
$\displaystyle=\Big{[}-\bm{\nabla}A^{0}-\partial^{0}\bm{A}\Big{]}+\frac{\bm{\Pi}}{~{}\Pi_{0}}\times\big{(}\bm{\nabla}\times\bm{A}\big{)}~{},$
(40)

further underscoring the analogy (a detailed derivation of the above standard
result can be found in Appendix I.4 of Sorensen (2021)).

##  Appendix B Derivation of the four-current and the energy-momentum tensor

This derivation largely follows Weber  _et al._ (1992), with two
modifications: in our case the vector field $A^{\mu}$ has no explicit momentum
dependence, and instead we are dealing with an arbitrary dependence of
$A^{\mu}$ on density.

We start with the BUU equation, Eq. (3), repeated here in a simplified
notation without the index $i$ enumerating the particle species,

$\Pi_{\mu}\partial^{\mu}_{x}f+\Pi_{\nu}\big{(}\partial_{\mu}^{x}A^{\nu}\big{)}\partial_{p}^{\mu}f=I_{\textrm{coll}}\,,$
(41)

where $f=f(x,p)$. We define $V=\frac{1}{2}\big{(}\Pi^{\mu}\Pi_{\mu}\big{)}$;
using $\Pi^{\mu}=p^{\mu}-A^{\mu}(x)$, we notice that

$\partial_{\mu}^{p}V=\Pi_{\mu}\hskip 14.22636pt\textrm{and}\hskip
14.22636pt\partial_{\mu}^{x}V=-\Pi_{\nu}\big{(}\partial_{\mu}^{x}A^{\nu}\big{)}~{},$
(42)

by means of which we rewrite the BUU equation as

$\big{(}\partial_{\mu}^{p}{V}\big{)}\partial_{x}^{\mu}f-\big{(}\partial_{\mu}^{x}V\big{)}\partial_{p}^{\mu}f=I_{\textrm{coll}}\,.$
(43)

To obtain the four-current, we integrate Eq. (43) over all possible four-
momenta, $g\int d^{4}p/(2\pi)^{4}$, where $g$ is the degeneracy factor. The
collision integral on the right-hand side of Eq. (43) vanishes due to particle
number conservation (an easy demonstration of this fact is shown in Appendix H
of Sorensen (2021)), and we are left with

$\displaystyle 0$
$\displaystyle=g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\bigg{(}\big{(}\partial_{\mu}^{p}{V}\big{)}\partial_{x}^{\mu}f-\big{(}\partial_{\mu}^{x}V\big{)}\partial_{p}^{\mu}f\bigg{)}$
$\displaystyle=g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\bigg{(}\big{(}\partial_{\mu}^{p}{V}\big{)}\partial_{x}^{\mu}f+\Big{[}\partial_{p}^{\mu}\big{(}\partial_{\mu}^{x}V\big{)}\Big{]}f\bigg{)}$
$\displaystyle=g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\partial_{x}^{\mu}\Big{[}\big{(}\partial_{\mu}^{p}{V}\big{)}~{}f\Big{]}$
$\displaystyle=\partial_{x}^{\mu}\bigg{[}g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\big{(}\partial_{\mu}^{p}{V}\big{)}~{}f\bigg{]}~{},$
(44)

where in the second and third equalities we used integration by parts and
$ab^{\prime}=(ab)^{\prime}-a^{\prime}b$, respectively. Recognizing the above
equality as the charge conservation equation $\partial^{\mu}j_{\mu}=0$, we can
identify the four-current as

$\displaystyle j_{\mu}$
$\displaystyle=g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\big{(}\partial_{\mu}^{p}{V}\big{)}~{}f~{}.$
(45)

At this point we turn our attention to the fact that the distribution function
$f$ is a function of eight independent variables $f=f(x^{\mu},p^{\mu})$. One
integrates over the energy $p^{0}$ by explicitly including the mass-shell
condition, Eq. (5), in the expression for $f(x,p)$:

$\displaystyle
f(x^{\mu},p^{\mu})=2(2\pi)\Theta(\Pi^{0})\delta\big{(}\Pi^{\mu}\Pi_{\mu}-m^{2}\big{)}~{}\tilde{f}(x^{\mu},\bm{p})~{}.$
(46)

With the identity
$\delta\big{[}g(x)\big{]}=\sum_{j}\frac{\delta(x-x_{j})}{|g^{\prime}(x_{j})|}$,
where $x_{j}$ are roots of $g(x)$, the mass-shell condition allows one to
write

$\displaystyle
g\int\frac{d^{4}\Pi}{(2\pi)^{4}}~{}f(x^{\mu},p^{\mu})=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{1}{\Pi^{0}}~{}\tilde{f}(x^{\mu},\bm{p})~{},$
(47)

where $\Pi^{0}=\sqrt{\bm{\Pi}^{2}+m^{2}}$. Consequently, additionally using
Eq. (42), we can also write the conserved current in the more well-known form

$\displaystyle j^{\mu}$
$\displaystyle=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{\Pi^{\mu}}{\Pi^{0}}~{}\tilde{f}~{}.$
(48)

The energy-momentum tensor is obtained in a similarly standard way, that is
via multiplying Eq. (43) by $p^{\mu}$ and integrating over $g\int
d^{4}p/(2\pi)^{4}$. The integral over $I_{\textrm{coll}}$ again vanishes, this
time due to the conservation of four-momentum in collisions (see Appendix H of
Sorensen (2021) for a simple demonstration), and one obtains

$\displaystyle 0$
$\displaystyle=g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\bigg{(}p^{\nu}\big{(}\partial_{\mu}^{p}{V}\big{)}\partial_{x}^{\mu}f-p^{\nu}\big{(}\partial_{\mu}^{x}V\big{)}\partial_{p}^{\mu}f\bigg{)}$
$\displaystyle=g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\bigg{(}\partial_{x}^{\mu}\Big{[}p^{\nu}\big{(}\partial_{\mu}^{p}{V}\big{)}f\Big{]}-\Big{[}\partial_{x}^{\mu}p^{\nu}\big{(}\partial_{\mu}^{p}{V}\big{)}\Big{]}f$
$\displaystyle\hskip
42.67912pt+~{}g^{\mu\nu}\big{(}\partial_{\mu}^{x}V\big{)}f+p^{\nu}\Big{[}\partial_{p}^{\mu}\big{(}\partial_{\mu}^{x}V\big{)}\Big{]}f\bigg{)}$
$\displaystyle=\partial_{x}^{\mu}~{}g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\big{(}\Pi^{\nu}+A^{\nu}\big{)}\big{(}\partial_{\mu}^{p}{V}\big{)}~{}f$
$\displaystyle\hskip
42.67912pt+~{}g^{\nu\mu}~{}g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\big{(}\partial_{\mu}^{x}V\big{)}~{}f~{},$
(49)

where in the first equality we again used
$ab^{\prime}=(ab)^{\prime}-a^{\prime}b$ and integration by parts, and in the
second equality we substituted $p^{\nu}=\Pi^{\nu}+A^{\nu}$. Because $A^{\nu}$
does not depend on momentum, we can pull it out of the integrals, which
together with Eq. (42) and the expression for the conserved four-current, Eqs.
(45) and (48), yield

$\displaystyle\partial_{x}^{\mu}\Bigg{[}g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\Pi^{\nu}\Pi_{\mu}~{}f+A^{\nu}j_{\mu}\Bigg{]}$
$\displaystyle\hskip
42.67912pt-~{}g^{\nu\mu}~{}j_{\lambda}\big{(}\partial_{\mu}^{x}A^{\lambda}\big{)}=0~{}.$
(50)

To extract the energy-momentum tensor from the above expression, one needs to
write the last term as a total derivative. Recalling that
$A^{\lambda}=\alpha(n)j^{\lambda}$ (where we suppress the baryon charge index
$B$) and $j_{\mu}j^{\mu}=n^{2}$, we can write

$\displaystyle j_{\lambda}\big{(}\partial_{\mu}^{x}A^{\lambda}\big{)}$
$\displaystyle=j_{\lambda}\Big{[}\partial_{\mu}^{x}\alpha(n)\Big{]}j^{\lambda}+j_{\lambda}\alpha(n)\Big{[}\partial_{\mu}^{x}j^{\lambda}\Big{]}$
$\displaystyle=n^{2}\Big{[}\partial_{\mu}^{x}\alpha(n)\Big{]}+\alpha(n)\bigg{(}\frac{1}{2}\partial_{\mu}^{x}\big{[}n^{2}\big{]}\bigg{)}$
$\displaystyle=\partial_{\mu}^{x}\Big{[}\alpha(n)n^{2}\Big{]}-\frac{1}{2}\alpha(n)\Big{(}\partial_{\mu}^{x}\big{[}n^{2}\big{]}\Big{)}~{}.$
(51)

The second term can be further rewritten as

$\displaystyle\alpha(n)n\Big{[}\partial_{\mu}^{x}n\Big{]}$
$\displaystyle=\bigg{[}\frac{\partial{}}{\partial{n}}\int_{0}^{n}dn^{\prime}~{}\alpha(n^{\prime})n^{\prime}\bigg{]}\bigg{[}\frac{\partial{n}}{\partial{x^{\mu}}}\bigg{]}$
$\displaystyle=\frac{\partial{}}{\partial{x^{\mu}}}\int_{0}^{n}dn^{\prime}~{}n^{\prime}\alpha(n^{\prime})~{}.$
(52)

Together with Eq. (51), the above allows one to write Eq. (50) as

$\displaystyle\partial^{x}_{\mu}\Bigg{[}g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\Pi^{\nu}\Pi^{\mu}~{}f+A^{\nu}j^{\mu}$
$\displaystyle\hskip
14.22636pt-~{}g^{\nu\mu}~{}\bigg{(}\alpha(n)n^{2}-\int_{0}^{n}dn^{\prime}~{}n^{\prime}\alpha(n^{\prime})\bigg{)}\Bigg{]}=0~{}.$
(53)

Because the energy-momentum tensor $T^{\mu\nu}$ satisfies
$\partial^{x}_{\nu}T^{\mu\nu}=0$, we can immediately identify

$\displaystyle T^{\mu\nu}$
$\displaystyle=g\int\frac{d^{4}p}{(2\pi)^{4}}~{}\Pi^{\mu}\Pi^{\nu}~{}f+A^{\mu}j^{\nu}$
$\displaystyle\hskip
14.22636pt-~{}g^{\mu\nu}~{}\bigg{(}nU(n)-\int_{0}^{n}dn^{\prime}~{}U(n^{\prime})\bigg{)}~{},$
(54)

where we substituted the expression for the single-particle rest-frame
potential as defined in Eq. (19). By using the mass-shell relation, Eq. (47),
substituting $d^{3}p\to d^{3}\Pi$ (which is allowed because $A^{\mu}$ does not
depend on momentum), and further using the fact that
$\Pi^{\mu}=p^{\mu}-A^{\mu}$ and the conserved current $j^{\mu}$ is given by
Eq. (48), one finally arrives at

$\displaystyle T^{\mu\nu}$
$\displaystyle=g\int\frac{d^{3}\Pi}{(2\pi)^{3}}~{}\frac{\Pi^{\mu}\Pi^{\nu}}{\Pi^{0}}~{}\tilde{f}+A^{\mu}j^{\nu}$
$\displaystyle\hskip
14.22636pt-~{}g^{\mu\nu}~{}\bigg{(}nU(n)-\int_{0}^{n}dn^{\prime}~{}U(n^{\prime})\bigg{)}~{}.$
(55)

The derived $T^{\mu\nu}$ agrees term by term with expressions in Sorensen and
Koch (2021), where they were derived in the special case of a polynomial form
of $\alpha(n)$.

##  Appendix C Thermodynamics and thermodynamic consistency

One can observe that the energy-stress tensor $T^{\mu\nu}$, derived in Eq.
(55), is symmetric in any frame, and in the Eckart frame, defined as the frame
in which $j^{\mu}=(n,0,0,0)$, it becomes diagonal,
$T^{\mu\nu}=\textrm{diag}(\mathcal{E},P,P,P)$. In that frame, Eq. (48) becomes
the well-known expression for the rest frame density,

$n=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\tilde{f}~{}.$ (56)

Up to this point, the form of the distribution function $\tilde{f}$ was not in
any way constrained (and, in particular, it could have been assumed to be
given by the test-particle ansatz, Eq. 6, making our formalism perfect for use
in hadronic transport). However, because now we are describing a system of
fermions in equilibrium, the distribution function $\tilde{f}$ can be assumed
to take the well-known form of the Fermi-Dirac distribution, where the single
particle energy includes the contribution from the mean-field potential $U(n)$
(see Appendix 1.C in Sorensen and Koch (2021) for a quick derivation, and Baym
and Pethick (2004) for a more complete discussion),

$\displaystyle\tilde{f}=\left(e^{\beta\big{(}\Pi^{0}+U(n)-\mu\big{)}}+1\right)^{-1}~{}.$
(57)

Furthermore, in the Eckart frame the field assumes the simple form
$A^{\mu}=\delta^{\mu}_{0}\alpha(n)n=\delta^{\mu}_{0}U(n)$, and the energy
density $\mathcal{E}$ and pressure $P$ are given by

$\displaystyle\mathcal{E}$
$\displaystyle=g\int\frac{d^{3}\Pi}{(2\pi)^{3}}~{}\Pi^{0}~{}\tilde{f}+\int_{0}^{n}dn^{\prime}~{}U(n^{\prime})~{},$
(58) $\displaystyle P$
$\displaystyle=g\int\frac{d^{3}\Pi}{(2\pi)^{3}}~{}\frac{\bm{\Pi}^{2}}{3\Pi^{0}}~{}\tilde{f}+nU(n)-\int_{0}^{n}dn^{\prime}~{}U(n^{\prime})~{}.$
(59)

Equations (56)-(59) fully describe the thermodynamics of our model, where the
parametrization of $U(n)$ can be chosen arbitrarily. The model is
thermodynamically consistent, and the system will evolve corresponding to the
test particle equations of motion, Eqs. (7)-(8).

We note that one can introduce a “shifted” chemical potential,

$\mu^{*}\equiv\mu-U(n)~{},$ (60)

known as the effective chemical potential. With $\mu^{*}$ at hand, it is
straightforward to see that the density, energy density, and the pressure can
be rewritten in the following way:

$\displaystyle n$ $\displaystyle=n_{\textrm{id}}(T,\mu^{*})~{},$ (61)
$\displaystyle\mathcal{E}$
$\displaystyle=\mathcal{E}_{\textrm{id}}(T,\mu^{*})+\int_{0}^{n}dn^{\prime}~{}U(n^{\prime})~{},$
(62) $\displaystyle P$
$\displaystyle=P_{\textrm{id}}(T,\mu^{*})+nU(n)-\int_{0}^{n}dn^{\prime}~{}U(n^{\prime})~{},$
(63)

where $n_{\textrm{id}}$, $\mathcal{E}_{\textrm{id}}$, and $P_{\textrm{id}}$
are the density, energy density, and pressure of an ideal Fermi gas with
chemical potential $\mu^{*}$, and the total energy density and pressure are
obtained by adding the field contributions. In particular, we can also see
that the entropy density is given by

$\displaystyle s(T,\mu)$
$\displaystyle\equiv\frac{1}{T}\Big{(}\mathcal{E}+P-\mu n\Big{)}$
$\displaystyle=\frac{1}{T}\Big{(}\mathcal{E}_{\textrm{id}}(T,\mu^{*})+P_{\textrm{id}}(T,\mu^{*})-\mu^{*}n\Big{)}~{},$
(64)

which is the entropy density of an ideal Fermi gas at the effective chemical
potential $\mu^{*}$.

As a first test of the thermodynamic consistency of Eqs. (61)-(63), we check
whether the thermodynamic relation $(\partial{P}/\partial{\mu})_{T}=n$ is
fulfilled. Integration by parts allows one to rewrite the ideal-Fermi-gas–like
contribution to the total pressure,

$\displaystyle P_{\textrm{id}}$
$\displaystyle=g\int\frac{d^{3}\Pi}{(2\pi)^{3}}~{}T~{}\ln\Big{[}1+e^{-\beta\big{(}\Pi^{0}+U(n)-\mu\big{)}}\Big{]}~{}.$
(65)

Using this form of $P$, it is straightforward to see that

$\displaystyle\left(\frac{\partial{P_{\textrm{id}}}}{\partial{\mu}}\right)_{T}$
$\displaystyle=g\int\frac{d^{3}\Pi}{(2\pi)^{3}}~{}T(-\beta)\bigg{[}\left(\frac{\partial{U(n)}}{\partial{\mu}}\right)_{T}-1\bigg{]}~{}\tilde{f}$
$\displaystyle=-n\left(\frac{\partial{U(n)}}{\partial{\mu}}\right)_{T}+n~{}.$
(66)

Differentiation of the remaining term in Eq. 59 yields

$\displaystyle\frac{\partial{}}{\partial{\mu}}\bigg{|}_{T}\left(nU(n)-\int_{0}^{n}dn^{\prime}~{}U(n^{\prime})\right)$
$\displaystyle=n\left(\frac{\partial{U(n)}}{\partial{\mu}}\right)_{T}~{},$
(67)

and, putting Eqs. (66) and (67) together, we confirm that
$\left(\frac{\partial{P}}{\partial{\mu}}\right)_{T}=n$ is satisfied.

For a second test, we similarly check that
$(\partial{P}/\partial{T})_{\mu}=s$. First, we calculate

$\displaystyle\left(\frac{\partial{P_{\textrm{id}}}}{\partial{T}}\right)_{\mu}$
$\displaystyle=\frac{P_{\textrm{id}}}{T}+g\int\frac{d^{3}\Pi}{(2\pi)^{3}}~{}\frac{\Pi^{0}+U(n)-\mu}{T}~{}\tilde{f}$
$\displaystyle\hskip
14.22636pt-~{}g\int\frac{d^{3}\Pi}{(2\pi)^{3}}~{}\left(\frac{\partial{U(n)}}{\partial{T}}\right)_{\mu}~{}\tilde{f}$
$\displaystyle=\frac{P_{\textrm{id}}}{T}+\frac{\mathcal{E}_{\textrm{id}}+nU(n)-\mu
n}{T}$ $\displaystyle\hskip
14.22636pt-~{}n\left(\frac{\partial{U(n)}}{\partial{T}}\right)_{\mu}~{}.$ (68)

A differentiation of the second term in the expression for the total pressure
yields

$\displaystyle\frac{\partial{}}{\partial{T}}\bigg{|}_{\mu}\left(nU(n)-\int_{0}^{n}dn^{\prime}~{}U(n^{\prime})\right)$
$\displaystyle=n\left(\frac{\partial{U(n)}}{\partial{T}}\right)_{\mu}~{},$
(69)

By putting the two above equations together and using Eq. (60) we immediately
obtain Eq. (64) and complete the proof of the thermodynamic consistency of the
model.

##  Appendix D Technical details of the simulations and tests of their
influence on flow observables

###  D.1 Density calculation

The rest frame baryon density is obtained from
$n_{B}^{2}=(j_{B})^{\mu}(j_{B})_{\mu}$, with the baryon current $j_{B}^{\mu}$
computed according to

$\displaystyle
j_{B}^{\mu}(\bm{r})=\frac{1}{N_{\textrm{ens}}}\frac{1}{N_{T}}\sum_{\textrm{ensembles}}\sum_{i=1}^{N_{T}}B_{i}\frac{p^{\mu}_{i}}{p^{0}_{i}}K_{G}(\bm{r}-\bm{r}_{i},u_{i})~{},$
(70)

where $B_{i}$ is the baryon number of the $i$th species,

$\displaystyle
K_{G}(\bm{r}-\bm{r}_{i},u,\sigma)=\frac{u_{0}}{(2\pi\sigma^{2})^{3/2}}\exp^{-\frac{(\bm{r}-\bm{r}_{i})^{2}+\left(\bm{u}\cdot(\bm{r}-\bm{r}_{i})\right)^{2}}{2\sigma^{2}}}$
(71)

is a Lorentz-contracted Gaussian smearing kernel as described in Oliinychenko
and Petersen (2016), with the width of the distribution set as $\sigma=1$ fm,
and

$\displaystyle u_{i}^{\mu}=(u_{i}^{0},\bm{u}_{i})=\frac{p_{i}^{\mu}}{m_{i}}$
(72)

is a test-particle’s four-velocity in the computational frame. The Gaussian
smearing kernel is cut off at $2\sigma$ and renormalized to give $\int
d^{3}r~{}K_{G}(\bm{r})=1$. Notice, however, that this does not guarantee that
the sum over lattice nodes, $\sum_{r\in\text{lattice}}K_{G}(r)$, is equal to
unity (this is because a sum over a finite number of values of $K_{G}(r)$ is
in practice a numerical integration of $K_{G}(r)$ using a rectangle rule,
which in the case of a Gaussian function introduces a numerical error); in
consequence, baryon number is not conserved on the lattice exactly. We have
also tested an alternative method of density calculation in which the smearing
kernel in each Cartesian direction is a triangular function with the base
proportional to the lattice spacing in that direction, often called
“triangular smearing”:

$\displaystyle
K_{t}(\bm{r}-\bm{r}_{i})=\frac{1}{(n^{3}l_{x}l_{y}l_{z})^{2}}g(\Delta
x)g(\Delta y)g(\Delta z)~{},$ (73)

with

$\displaystyle g(\Delta q)\equiv\big{(}nl_{q}-|\Delta
q|\big{)}\theta\Big{(}nl_{q}-|\Delta q|\Big{)}~{},$ (74)

where $l_{x}$, $l_{y}$, and $l_{z}$ are lattice spacings, and $n$ determines
the range of smearing in units of lattice spacings; we used $n=2$. The
advantage of the triangular smearing is that
$\sum_{r\in\text{lattice}}K_{t}(r)$ is equal to unity by construction (as
integration by the rectangle rule is exact for a linear function), while the
disadvantage is that it is not Lorentz covariant. Therefore, when triangular
smearing is used, one needs to Lorentz-contract the lattice itself in the beam
direction with the $\gamma$ factor corresponding to the beam energy (for
simulations performed in the center-of-mass frame of the colliding nuclei). In
the following, we will show that both types of smearing lead to comparable
results provided that $N_{\textrm{ens}}$ and $N_{T}$ are large enough to
ensure a reliable density calculation; in this way, we show that our results
do not depend on the technical implementation of the density calculation.

We note here that both of our smearing paradigms (and, in fact, any smearing
paradigm that we know of) lead to formal violations of causality. This is
because the smearing kernels for calculating the density at a space-time point
$(t,\bm{r})$ depend on the positions of the test particles in the vicinity of
the point $\bm{r}$ at the same time $t$, implying instantaneous propagation of
information about density changes within the smearing range.

For both types of smearing, the field derivatives $\partial^{\mu}A_{\nu}$ are
computed numerically on the lattice, including the time derivative terms. The
smearing kernels, mixed ensembles mode, and the equations of motion with the
relativistic vector mean field $A^{\mu}(n_{B})$ based on a polynomial form of
$\alpha(n_{B})$ (the VDF model) are included in the publicly available SMASH
2.1 software sma ; however, the parametrization of the vector mean field
$A^{\mu}(n_{B})$ with an arbitrary $c_{s}^{2}(n_{B},T=0)$ is not.

###  D.2 Effects related to $N_{T}$ and $N_{\textrm{ens}}$

Figure 15: Effect of the number of test particles per particle $N_{T}$ on the
proton elliptic flow (a) and the slope of the proton directed flow (b) in
Au+Au collisions at $E_{\textrm{lab}}=1$ GeV/nucleon. The product
$N_{T}N_{\textrm{ens}}=1000$ is kept constant, therefore the observed effect
is primarily due to the dependence of collisions on $N_{T}$. Based on this
result, for further simulations we choose $N_{T}=20$. Figure 16: Effect of the
number of ensembles $N_{\textrm{ens}}$ on the proton elliptic flow (a) and the
slope of the proton directed flow (b) in Au+Au collisions at
$E_{\textrm{lab}}=1$ GeV/nucleon. The number of test particles per particle is
set to $N_{T}=20$, which, together with the Gaussian smearing, seems to
provide a rather accurate result even at $N_{\textrm{ens}}=1$ (see text for
more details on the dependence of this result on the smearing paradigm).
However, for better accuracy, we choose $N_{\textrm{ens}}=50$ for further
simulations.

We performed a comparison of results from simulations utilizing different
numbers of test particles per particle $N_{T}$. To single out the effect of
$N_{T}$ on collisions only and not on the accuracy of the density calculation
(which we computed using the Gaussian smearing kernel), we kept the product
$N_{T}N_{\textrm{ens}}$ constant. Because the density and potentials are
computed based on all test particles in all ensembles, by keeping
$N_{T}N_{\textrm{ens}}=\textrm{const}$ we do not change the magnitude of the
numerical spatial fluctuations in density, and thus we obtain the same
accuracy of potential gradient calculations. In Fig. 15, one can see that the
effect of $N_{T}$ on proton flow observables is very substantial: as we show
in the main text, the change of $v_{2}(y^{\prime}=0)$ from -0.07 to -0.05 at
$E_{\textrm{lab}}=1$ GeV/nucleon due to varying $N_{T}$ from 1 to 10 is
comparable to a difference in $v_{2}(y^{\prime}=0)$ obtained due to using a
soft and a very stiff EOS. As the product $N_{T}N_{\textrm{ens}}$ is kept
constant in Fig. 15, the observed effect of $N_{T}$ is entirely due to the
dependence of collisions on $N_{T}$. The dependence of flow on $N_{T}$ is
attributed to spurious effects of collisions at nonzero range Cheng  _et al._
(2002), which become smaller at larger $N_{T}$ because cross sections scale
with $N_{T}$ according to $\sigma\to\sigma/N_{T}$ and, therefore, the
geometrical distance at which collisions occur decreases for larger $N_{T}$.
In agreement with this interpretation, we observe a saturation of
$v_{2}(y^{\prime}=0)$ as a function of $N_{T}$ which occurs at $N_{T}\approx
15$ (while keeping $N_{T}N_{\textrm{ens}}=1000$). Our results show that any
$N_{T}>15$ is a reasonable choice, and we take $N_{T}=20$ for our further
simulations.

When $N_{T}$ is fixed, increasing the number of ensembles $N_{\textrm{ens}}$
changes our results to a lesser extent, especially in the case of the elliptic
flow, as shown in Fig. 16. This is likely because the Gaussian smearing
together with the oversampling factor $N_{T}=20$ provide an acceptably smooth
spatial density profile even for $N_{\textrm{ens}}=1$, and the results
saturate already at $N_{\textrm{ens}}=5$. We found that the computational cost
increases almost proportionally to $N_{\textrm{ens}}$ (specifically, it is
slightly supralinear), while the simulation statistics increases exactly
proportionally to $N_{\textrm{ens}}$. Given the minimal additional computing
cost, we chose to be on the safe side and we took $N_{\textrm{ens}}=50$.

We note here that while similar tests of the influence of $N_{T}$ on flow
observables have been performed within SMASH in the past (see, e.g., Mohs  _et
al._ (2022)), the analysis presented here is qualitatively different as it
separates the influence of scaled-down collision cross sections from the
influence of the statistics used for calculating the mean field. In
particular, as can be seen in Figs. 15 and 16, the two effects have an
opposite influence on the magnitude of the elliptic flow in the explored
collision energy region. In contrast, Mohs  _et al._ (2022) shows the effect
of varying both the collision cross sections and the mean-field calculation
statistics at the same time. Previous works exploring the influence of
spurious fluctuations in the density calculation on the collision integral and
the mean-field calculation include Zhang  _et al._ (2018b); Colonna  _et al._
(2021); however, they do not address the effects on flow observables
explicitly.

We performed a similar test for simulations using the triangular smearing for
the density calculation, and in this case the results suggested that this type
of smearing demands a larger value of $N_{\textrm{ens}}$; this is because the
triangular smearing kernel is, by construction, less “smooth” than the
Gaussian smearing kernel, and so it requires more statistics to yield a smooth
density calculation. Nevertheless, in Fig. 17 we show that we obtain similar
results for the proton flow using triangular and Gaussian smearing at
$N_{T}=20$ and $N_{\textrm{ens}}=50$. This is a non-trivial observation,
because technically the two smearing paradigms are very different: in the
Gaussian smearing the kernel is contracted along the direction of motion of a
test particle, while in the triangular smearing the lattice is contracted
along the beam direction. Naturally, since physical results should not be
sensitive to the particular implementation of the density calculation,
comparing the results obtained using two different smearing kernels is a
quality test of our simulations. Finally, we note that the computing time is
almost the same for simulations using the two smearing kernels, with a 4%
increase in the case of the triangular smearing; although the triangular
smearing uses a kernel that is less computationally costly, in that case the
lattice, which is more dense in the $z$ direction due to being contracted,
needs a larger number of nodes in order to cover the same volume.

Figure 17: Comparison of the proton elliptic flow (a) and the slope of the
proton directed flow (b) in Au+Au collisions as a function of the collision
energy, obtained with simulations using different types of smearing kernels
for computing baryon density on the Cartesian lattice.

###  D.3 Nucleus initialization

Another important factor influencing the proton flow as obtained from the
hadronic transport simulations is the initialization of the nuclei. We employ
a default SMASH initialization as described in Weil  _et al._ (2016): in the
coordinate space each nucleon is sampled independently from the Woods-Saxon
distribution, while in the momentum space the nucleon momentum is sampled from
a uniform distribution in the Fermi-sphere, where the radius of the Fermi
sphere is computed from the local coordinate space density obtained
analytically from the Woods-Saxon distribution. This simple model of nucleus
initialization does not account for two physical phenomena: First, by
construction, the mean field does not influence the initialization of a
nucleus, and therefore the nucleus is not initialized exactly in a ground
state. In freely moving nuclei, this results in “breathing mode” oscillations
with the oscillation period on the order of 70–80 fm/$c$, as shown in Fig. 7
of Weil  _et al._ (2016). In simulations of heavy-ion collisions, the maximum
compression is reached at the latest at 13 fm/$c$ (for
$E_{\textrm{lab}}=1~{}\textrm{GeV/nucleon}$), as one can see in Fig. 1, which
is too short for the spurious oscillations to fully develop and thus does not
lead to sizable contributions (especially given the fact that the differences
in density due to the spurious oscillations are dwarfed by the changes in
density due to the compression of the nuclei). Second, correlated nucleon
pairs (as well as any other correlations) are not accounted for in the nucleus
initialization, causing the high-momentum tail of the momentum distribution to
be absent in our model, where the maximal momentum is limited by the Fermi
momentum at the center of the nucleus. While we do not expect these phenomena
to be important for the proton flow, we note that they may contribute to the
overall systematic uncertainty of our simulations.

###  D.4 Pauli blocking

In simulations of a dense baryonic medium one might expect Pauli blocking to
influence the observables. However, even at
$E_{\textrm{lab}}=1~{}\textrm{GeV/nucleon}$ the related effects turn out to be
virtually absent. We find no significant difference in proton spectra or flow
in simulations with and without Pauli blocking, which is in agreement with
Weil  _et al._ (2016). The implementation of Pauli blocking in SMASH
(described in Weil  _et al._ (2016)) is based on estimating of the phase space
density at the point of an action (where an action means a collision or a
decay) and rejecting the action with probability $\prod_{j}(1-f_{j})$, where
$f_{j}$ is the phase space density of particles $j$ exiting the action (a
collision or a decay). We find that in Au+Au collisions at
$E_{\textrm{lab}}=1~{}\textrm{GeV/nucleon}$ with the impact parameter $b=6$
fm, around 8% of collisions and decays are blocked; at
$E_{\textrm{lab}}=8~{}\textrm{GeV/nucleon}$ this number is reduced to around
3%.

Figure 18: Sensitivity of the proton flow to scaling of all cross sections in
Au+Au collisions at $b=6$ fm, simulated using SMASH. At $n_{B}\in(0,2)n_{0}$,
the mean-field is parametrized to reproduce the default Skyrme EOS, while at
higher densities the mean field is parametrized to yield
$c_{s}^{2}\big{[}n_{B}\in(2,3)n_{0}\big{]}=0.1$ and
$c_{s}^{2}\big{[}n_{B}>3n_{0}\big{]}=0.3$.

###  D.5 Absence of isospin, Coulomb, and momentum-dependent interactions

It is important to mention that the vector mean field $A^{\mu}$ is the only
field implemented in our simulations. Coulomb interactions are not included,
which may potentially influence our pion flow results Adamczewski-Musch  _et
al._ (2022). Isospin-dependent potentials likewise are not included, i.e.,
protons and neutrons feel the same mean field. At lower energies, in
particular below $E_{\textrm{lab}}=0.5~{}\textrm{GeV/nucleon}$, this cannot be
justified, but at $E_{\textrm{lab}}=1~{}\textrm{GeV/nucleon}$ and above it is
a reasonable approximation due to the fact that matter at midrapidity is more
isospin symmetric in that case and also because at high baryon density, the