$F_{\textrm{c}}=1/3$. Such weighted contributions are then further scaled by
the number of test particles per particle $N_{T}$ (if it hasn’t been already
taken into the account in the charge of the test particle) and the volume of
the cells $V_{c}$ to obtain the baryon current at each node. This smearing
scheme has the advantage that the integration of the baryon density calculated
on the lattice equals the number of baryon charges present within the bounds
of the lattice exactly, meaning that the smearing conserves the baryon number.
Additionally, it is very computationally efficient. However, at the same time
it often requires using a very large number of test particles per particle
$N_{T}$, on the order of $N_{T}\approx 1000$, to provide reliable results.
This considerably lowers the efficiency of this scheme in the full ensemble
mode.

####  2 Triangular smearing

In the smearing scheme referred to in SMASH as “triangular smearing”,
originating from the lattice Hamiltonian method of solving nuclear dynamics
[203], each test particle contributes to the baryon current at a given node
with a smearing factor

$\displaystyle
S(\bm{r}_{i}-\bm{r})=\frac{B}{N_{T}\left(nL_{x}\right)^{2}\left(nL_{y}\right)^{2}\left(nL_{z}\right)^{2}}g(\Delta
x)g(\Delta y)g(\Delta z)~{},$ (11)

where

$\displaystyle g(\Delta q)=\Big{(}nL_{q}-|\Delta
q|\Big{)}\theta\Big{(}nL_{q}-|\Delta q|\Big{)}~{}.$ (12)

Here, $\bm{r}_{i}$ is the position of the node, $\bm{r}$ is the position of
the test particle, $B$ is the baryon charge of the contributing test-particle
(not scaled by $N_{T}$), $L_{x}$, $L_{y}$, $L_{z}$ are the lattice spacings in
the $x$-, $y$-, and $z$-direction, the smearing parameter $n$ is an integer
which determines the range of $S$ in terms of the lattice spacings, and
$\Delta x$, $\Delta y$, $\Delta z$ are the $x$, $y$, and $z$ components of
$\bm{r}_{i}-\bm{r}$. Because the triangular smearing is linear, the summation
of the resulting baryon density over volume integrates to the number of
baryons present within the bounds of the lattice exactly (this is because
Euler integration is exact for a linear function). Additionally, the
triangular smearing is very efficient numerically, and a smooth density
calculation is obtained for a standard number of test particles per particle
on the order of $N_{T}\approx 100$.

####  3 Gaussian smearing

Finally, the type of smearing referred to in SMASH as “Gaussian smearing” and
present in the original release of SMASH uses a Lorentz covariant Gaussian
kernel as its basis [204]. Starting from a smearing kernel in the rest frame
of a test particle defined as

$\displaystyle
K(\bm{r}-\bm{r}_{i})=\frac{1}{\big{(}2\pi\sigma^{2}\big{)}^{3/2}}\exp\left(-\frac{\big{(}\bm{r}-\bm{r_{i}}\big{)}^{2}}{2\sigma^{2}}\right)~{},$
(13)

one can perform a boost to the computational frame in the usual manner.
Additionally adopting a normalization ensuring that volume integrations over
$K(\bm{r}-\bm{r}_{i})$ are also Lorentz covariant leads to the expression for
the baryon current in the computation frame of the form

$\displaystyle
j^{\mu}_{\textrm{comp}}(\bm{r})=\frac{\gamma}{N_{T}\big{(}2\pi\sigma^{2}\big{)}^{3/2}}\sum_{i=1}^{N}B_{i}\frac{\Pi^{\mu}_{i}}{\Pi_{0,i}}\exp\left(-\frac{\big{(}\Lambda_{i}(\bm{r}-\bm{r}_{i})\big{)}^{2}}{2\sigma^{2}}\right)~{},$
(14)

where $\gamma=\big{(}1-v^{2}\big{)}^{-1/2}$ is the Lorentz factor, the sum is
performed over all test particles, $B_{i}$ is the baryon charge of the $i$-th
test particle (not scaled by $N_{T}$), and $\Lambda_{i}$ is the Lorentz
transformation matrix from the test particle’s rest frame to the computational
frame. The Gaussian smearing is characterized by two parameters: the first of
them is the Gaussian width, usually taken to be $\sigma=1$ fm, and the second
is a finite cut-off radius expressed in units of the Gaussian width, often
taken to be $r_{\textrm{c}}=4\sigma$, which is the radius beyond which the
test particle’s contribution to the density is considered negligible. This
cut-off is introduced to save the computation time and not evaluate near-zero
contributions to density at distances far from a given test particle; for a
given chosen $r_{\textrm{c}}$, the Gaussian smearing kernel is again
normalized such that it integrates to $1/N_{T}$, as it should.

While the use of exponentials is numerically very costly, the Gaussian
smearing has the advantage in that it is Lorentz covariant and therefore
properly treats relativistic effects through the contraction of the smearing
kernel in the direction of a test particle’s motion (in other smearing
schemes, these effects are often dealt with by using ad-hoc methods such as,
for example, increasing the fineness of the lattice in a direction in which
the system is expected to be contracted). Additionally, using the fact that
the Gaussian smearing kernel is everywhere differentiable, in SMASH the
derivatives of the baryon current can be calculated analytically (naturally,
only in the case of the Gaussian smearing). For example, the gradient of the
baryon current, Eq. (14), is given by

$\displaystyle\frac{dj^{\mu}}{d\bm{r}}=\frac{\gamma}{N_{T}\big{(}2\pi\sigma^{2}\big{)}^{3/2}}\sum_{i=1}^{N}B_{i}\frac{\Pi^{\mu}_{i}}{\Pi_{0,i}}\exp\left(-\frac{(\Lambda_{i}(\bm{r}-\bm{r}_{i}))^{2}}{2\sigma^{2}}\right)\frac{\Lambda_{i}^{2}(\bm{r}-\bm{r}_{i})}{\sigma^{2}}~{},$
(15)

where we can see that the derivative has been applied directly to the smearing
function. Consequently, these “Gaussian gradients” do not depend on the size
of the lattice spacing, but rather on the Gaussian smearing parameters
$\sigma$ and $r_{\textrm{c}}$. In certain situations this may prove to be
problematic: tests have shown that because of this feature, rapidly varying
gradients are not reproduced well due to the fact that the analytic
expressions for derivatives are “smoothed out” over large distances.
Additionally, this derivative scheme misses a possibly important term when
mean fields are present, as in this case the kinetic momentum has itself a
spacetime-dependence through
$\Pi_{i}^{\mu}=p_{i}^{\mu}-\sum_{k}A_{k}^{\mu}(x)$, see Eq. (76); we show this
explicitly below.

Finally, the complex functional form of the Gaussian smearing factor means
that the sum of its contributions to discrete lattice nodes does not add up to
$1/N_{T}$ (this is because the Euler integration is not exact in the case of a
smearing kernel of the form Eq. (13)), which means that the baryon number is
not completely preserved on the lattice.

####  4 Gaussian smearing and the time derivative of vector baryon current

Let us consider the time derivative of the vector baryon current calculated
within the Gaussian smearing scheme, Eq. (14), which we rewrite here with
explicit time dependence,

$\displaystyle\bm{j}\big{(}\bm{x},t\big{)}=\sum_{i}\frac{\bm{\Pi}_{i}(t)}{\Pi_{0,i}(t)}K\Big{(}\bm{x}-\bm{x}_{i}(t)\Big{)}~{},$
(16)

where $\bm{\Pi}_{i}$ and $\Pi_{0,i}$ denote the spatial and temporal
components of the kinetic momentum $\Pi^{\mu}$, respectively, and $K$ is the
smearing factor which includes the baryon number coefficient $B_{i}$ as well
as the normalization for clarity of calculations. The current at the time
$t+\Delta t$ is given by

$\displaystyle\bm{j}\big{(}\bm{x},t+\Delta
t\big{)}=\sum_{i}\frac{\bm{\Pi}_{i}(t+\Delta t)}{\Pi_{0,i}(t+\Delta
t)}K\Big{(}\bm{x}-\bm{x}_{i}(t+\Delta t)\Big{)}~{},$ (17)

which we can expand around $t$,

$\displaystyle\bm{j}\big{(}\bm{x},t+\Delta t\big{)}$ $\displaystyle=$
$\displaystyle\sum_{i}\frac{\bm{\Pi}_{i}(t)+\Delta
t\frac{\partial{\bm{\Pi}_{i}}}{\partial{t}}}{\Pi_{0,i}(t)+\Delta
t\frac{\partial{\Pi_{0,i}(t)}}{\partial{t}}}\left[K\Big{(}\bm{x}-\bm{x}_{i}(t)\Big{)}+\Delta
t\frac{\partial{\bm{x}_{i}}}{\partial{t}}\frac{\partial{K}}{\partial{\bm{x}_{i}}}\right]+\mathcal{O}\Big{(}(\Delta
t)^{2}\Big{)}$ (18) $\displaystyle\approx\sum_{i}\frac{\bm{\Pi}_{i}(t)+\Delta
t\frac{\partial{\bm{\Pi}_{i}}}{\partial{t}}}{\Pi_{0,i}(t)}\left(1-\frac{\Delta
t}{\Pi_{0,i}}\frac{\bm{\Pi}_{i}}{\Pi_{0,i}}\frac{\partial{\bm{\Pi}_{i}(t)}}{\partial{t}}\right)\left[K\Big{(}\bm{x}-\bm{x}_{i}(t)\Big{)}+\Delta
t\frac{\partial{\bm{x}_{i}}}{\partial{t}}\frac{\partial{K}}{\partial{\bm{x}_{i}}}\right]~{}.$

Keeping only terms up to linear order in $\Delta t$ results in

$\displaystyle\bm{j}\big{(}\bm{x},t+\Delta t\big{)}$ $\displaystyle\approx$
$\displaystyle\sum_{i}\frac{\bm{\Pi}_{i}(t)}{\Pi_{0,i}(t)}K\Big{(}\bm{x}-\bm{x}_{i}(t)\Big{)}$
(19) $\displaystyle+~{}\Delta
t\sum_{i}\frac{1}{\Pi_{0,i}(t)}\left[\frac{m_{i}^{2}}{\Pi_{0,i}^{2}}\right]\frac{\partial{\bm{\Pi}_{i}(t)}}{\partial{t}}K\Big{(}\bm{x}-\bm{x}_{i}(t)\Big{)}+\Delta
t\sum_{i}\frac{\bm{\Pi}_{i}(t)}{\Pi_{0,i}(t)}\left[\frac{\partial{\bm{x}_{i}}}{\partial{t}}\frac{\partial{K}}{\partial{\bm{x}_{i}}}\right]~{},$

where we have used the fact that $\Pi_{0,i}^{2}=\bm{\Pi}_{i}^{2}+m^{2}$. We
can then calculate

$\displaystyle\frac{\partial{}}{\partial{t}}\bm{j}\big{(}\bm{x},t\big{)}$
$\displaystyle\equiv$ $\displaystyle\lim_{\Delta t\to
0}\frac{\bm{j}\big{(}\bm{x},t+\Delta
t\big{)}-\bm{j}\big{(}\bm{x},t\big{)}}{\Delta t}$ (20) $\displaystyle=$
$\displaystyle\sum_{i}\frac{1}{\Pi_{0,i}(t)}\left[\frac{m_{i}^{2}}{\Pi_{0,i}^{2}}\right]\frac{\partial{\bm{\Pi}_{i}(t)}}{\partial{t}}K\Big{(}\bm{x}-\bm{x}_{i}(t)\Big{)}+\sum_{i}\frac{\bm{\Pi}_{i}(t)}{\Pi_{0,i}(t)}\left[\frac{\partial{\bm{x}_{i}}}{\partial{t}}\frac{\partial{K}}{\partial{\bm{x}_{i}}}\right]~{}.$

By comparing with Eq. (15), adapted to the case of the derivative with respect
to time, one can see that the analytical derivative used in the Gaussian
smearing scheme includes only the second term in the above expression. This,
in particular, misses contributions coming directly from changes in the
momentum due to the mean-field potential.

####  5 Chain rule and direct finite difference derivatives

A convenient form of the equation of motion for the kinetic momentum in the
VSDF model, Eq. (80), has been derived in Appendix 15.D, Eq. (58), which we
rewrite here for the VDF case,

$\displaystyle\frac{d\bm{\Pi}}{dt}=\sum_{k=1}^{K}\left\\{-\Big{(}\bm{\nabla}A_{k}^{0}+\partial^{0}\bm{A}_{k}\Big{)}+\frac{\bm{\Pi}}{\Pi_{0}}\times\big{(}\bm{\nabla}\times\bm{A}_{k}\big{)}\right\\}~{}.$
(21)

The finite difference derivatives of the vector field
$\sum_{k=1}^{K}A_{k}^{\mu}$, where $A_{k}^{\mu}$ is given by Eq. (90),
repeated here for convenience

$\displaystyle A_{k}^{\mu}(x;C_{k},b_{k})=C_{k}n_{B}^{b_{k}-2}j^{\mu}~{},$
(22)

can be performed in two ways. We can apply the chain rule to the derivative
expressions in Eq. (21), obtaining

$\displaystyle\sum_{k=1}^{K}\bm{\nabla}A_{k}^{~{}0}=\sum_{k=1}^{K}C_{k}n_{B}^{b_{k}-3}\bigg{(}(b_{k}-2)\Big{(}\bm{\nabla}n_{B}\Big{)}j^{0}+n_{B}\Big{(}\bm{\nabla}j^{0}\Big{)}\bigg{)}~{},$
(23)
$\displaystyle\sum_{k=1}^{K}\frac{\partial{}}{\partial{t}}\bm{A}_{k}=\sum_{k=1}^{K}C_{k}n_{B}^{b_{k}-3}\bigg{(}(b_{k}-2)\left(\frac{\partial{n_{B}}}{\partial{t}}\right)\bm{j}+n_{B}\left(\frac{\partial{\bm{j}}}{\partial{t}}\right)\bigg{)}~{},$
(24)
$\displaystyle\sum_{k=1}^{K}\bm{\nabla}\times\bm{A}_{k}=\sum_{k=1}^{K}C_{k}n_{B}^{b_{k}-3}\bigg{(}(b_{k}-2)\Big{(}\bm{\nabla}n_{B}\Big{)}\times\bm{j}+n_{B}\Big{(}\bm{\nabla}\times\bm{j}\Big{)}\bigg{)}~{},$
(25)

and then use the finite difference derivatives of the rest frame density
$n_{B}$ and current $j^{\mu}$ at a given node to obtain
$\sum_{k=1}^{K}\bm{\nabla}A_{k}^{~{}0}$,
$\sum_{k=1}^{K}\partial{\bm{A}_{k}}/\partial{t}$, and
$\sum_{k=1}^{K}\bm{\nabla}\times\bm{A}_{k}$ at that node. Alternatively, we
can use the values of $n_{B}$ and $j^{\mu}$ at a given node to compute
$\sum_{k=1}^{K}A_{k}^{\mu}$ at that node, and then compute the finite
difference derivatives of the relevant terms. As an example, let us consider
the gradient of the zeroth component of the field
$\sum_{k=1}^{K}\bm{\nabla}A_{k}^{0}$. The $z$-direction component of the
gradient is then simply given by

$\displaystyle\sum_{k=1}^{K}\nabla_{z}A_{k}^{0}(z_{0})=\frac{\sum_{k=1}^{K}A_{k}^{0}(z_{0}+\Delta
z)-\sum_{k=1}^{K}A_{k}^{0}(z_{0}-\Delta z)}{2\Delta z}~{},$ (26)

where $\Delta z$ is the lattice spacing the $z$-direction.

We refer to the former method of calculating derivatives as the chain rule
finite difference derivative and to the latter as the direct finite difference
derivative. In our studies, the direct finite difference derivative has proven
to be numerically more stable, as evaluated by the degree to which energy is
conserved in the simulation.

### 3\. Potentials in SMASH

SMASH was introduced with a density-based mean-field potential of the Skyrme
parametrization type, which we will refer to simply as the “Skyrme” potential.
Based on a semi-relativistic generalization of the rest-frame Skyrme
potential, the equation of motion for the kinetic momentum as used in the
SMASH Skyrme mean-field mode is given by

$\frac{d\bm{\Pi}}{dt}=\frac{\partial U}{\partial
n_{B}}\bigg{[}-\bigg{(}\bm{\nabla}j^{0}+\frac{\partial{\bm{j}}}{\partial{t}}\bigg{)}+\frac{d\bm{x}}{dt}\times\big{(}\bm{\nabla}\times\bm{j}\big{)}\bigg{]},$
(27)

where $n_{B}\equiv\sqrt{j_{\mu}j^{\mu}}$ is the rest frame baryon density and
$U$ is the rest frame single-particle potential given by

$\displaystyle
U(n_{B})=a\left(\frac{n_{B}}{n_{0}}\right)+b\left(\frac{n_{B}}{n_{0}}\right)^{\tau}~{},$
(28)

where $n_{0}$ is the saturation density, so that

$\displaystyle\frac{\partial U}{\partial
n_{B}}=\frac{a}{n_{0}}+b\tau\frac{n_{B}^{\tau-1}}{n_{0}^{\tau}}~{}.$ (29)

(We note here that SMASH allows the user to vary the values of the Skyrme
parameters $a$, $b$, and $\tau$, but uses a fixed value of saturation density
$n_{0}=0.168\ \textrm{fm}^{-3}$.) Importantly, while Eq. (27) has a form
similar to that of the equations of motion in relativistic electrodynamics, it
is not fully relativistic, as we show below.

Within the work presented in this thesis, we implemented the VDF potential and
the corresponding fully relativistic equations of motion in SMASH. Using the
VDF model, it is straightforward to reproduce the Skyrme potential of the form
given in Eq. (28) by using the vector field
$\sum_{k=1}^{K}\bm{\nabla}A_{k}^{\mu}=\sum_{k=1}^{K}\tilde{C}_{k}(n_{B}/n_{0})^{b_{k}-2}(j^{\mu}/n_{0})$
(see Appendix 18 for more details on the conversion of the VDF interaction
coefficients to units of energy) with

$\displaystyle K=2~{},\hskip 14.22636pt\tilde{C}_{1}=a~{},\hskip
14.22636ptb_{1}=2~{},\hskip 14.22636pt\tilde{C}_{2}=b~{},\hskip
14.22636ptb_{2}=\tau+1~{},\hskip 14.22636ptn_{0}=0.168\ \textrm{fm}^{-3}~{}.$
(30)

This possibility is very useful for comparisons of simulation results obtained
using the two potentials. On the theoretical level, comparing Eq. (27) with
Eq. (21), supported by Eqs. (23-25), immediately reveals that the SMASH Skyrme
equation of motion, Eq. (27), misses contributions proportional to
$(\bm{\nabla}n_{B})j^{0}$, $(\partial{n_{B}}/\partial{t})\bm{j}$, and
$\bm{\nabla}n_{B}\times\bm{j}$. Notably, the VDF equation of motion, Eq. (21),
becomes the SMASH equation of motion, Eq. (27), for $j^{0}\approx n_{B}$,
which is the nonrelativistic limit.

In addition to the Skyrme potential, SMASH also includes the option to use a
symmetry potential, which is of relevance for systems characterized by a non-
zero isospin asymmetry. This possibility is at this time only available for
simulations using the Skyrme potential and therefore will not be discussed
here.

### 4\. Comparison of different numerical procedures

As part of the implementation of the VDF equations of motion and different
baryon current calculation schemes in SMASH, we thoroughly compared their
behavior against the original SMASH Skyrme potentials and the Gaussian
smearing scheme. In particular, we extensively studied deviations from perfect
energy conservation in systems initialized at different points of the phase
diagram, with different EOSs, characterized by various geometries, using
different smearing and derivative calculation schemes, and employing different
mean fields.

We note that problems related to energy conservation are present in many
hadronic transport codes and can be directly linked not only to the
integration algorithm for the equations of motion, but also to the particular
details of density and density gradient calculations (including the influence
of statistical noise fluctuations on the magnitude of local density gradients)
as well as to the readout of these quantities from the discrete lattice.
Although methods ensuring energy conservation are known [205] for transport
simulations dealing with nonrelativistic systems, an application of such
approaches in relativistic transport codes has not been attempted.

One of the best ways to test the robustness of the integration of equations of
motion is to study systems evolving under extreme conditions. An example of
such a situation is uniform nuclear matter initialized at a low temperature
inside the spinodal region of the nuclear phase transition. This system is
mechanically unstable, and because the initialized matter is not perfectly
uniform due to the finite number of test particles used to sample the uniform
density profile, there are initial density fluctuations which over time will
be magnified by the mean field, leading to a spontaneous separation into two
coexisting phases: a nuclear liquid and a gas of nucleons (more details on
this and other simulations mentioned in this section will be provided in
Chapter 4). This separation causes the system to exhibit steep density
gradients, and the degree to which the density calculation and the integration
of the equations of motion deals well with these gradients is reflected in the
deviation of the energy conservation from the ideal case, where the energy of
the whole system is conserved over time, $dE_{\textrm{system}}/dt=0$. This
deviation is best measured as the energy difference per particle, denoted here
as $\Delta E={(E_{\textrm{final}}-E_{\textrm{initial}})}/{N}$.

While the total energy of the system should stay constant, the kinetic energy
and the mean-field energy are allowed to change their values, and in fact it
is expected that they do so, since the non-optimal initial value of the mean-
field energy is what drives the system to undergo a phase separation in the
first place. Here, the evolution of the system can be, to a zeroth degree,
captured within the ratio of the final mean-field energy of the system to the
initial mean-field energy,
$R_{\textrm{MF}}=E_{\textrm{MF,final}}/E_{\textrm{MF,initial}}$.

Finally, an important part of testing a given scheme for integrating the
equations of motion is keeping track of the amount of time it takes to
simulate a single event $t_{\textrm{sim}}$, as typical simulations of heavy-
ion collisions require on the order of 1,000 or even 10,000 events. For the
ballpark calculations presented here we used a personal computer with a 2.6
GHz processor and 16 GB of memory.

In Tables 1, 2, and 3 we show the values of $\Delta E$ (in
$\textrm{MeV}/\textrm{nucleon}$), $R_{\textrm{MF}}$, and $t_{\textrm{sim}}$
(in minutes), obtained from simulations initializing uniform nuclear matter in
a box of side $10\ \textrm{fm}$ with periodic boundary conditions (SMASH Box
modus) at $T=1\ \textrm{MeV}$ and $n_{B}=0.25\ n_{0}$, and averaged over 10
events. The lattice used for the density calculation had 10 nodes per side, so
that $L_{x}=L_{y}=L_{z}=1\ \textrm{fm}$, and we turned off all collisions and
decays between the test particles to focus on the effects of the mean-field
calculation. We tested different smearing schemes: Gaussian, triangular, and
discrete, at different values of their respective parameters (with the
Gaussian width always equal $\sigma=1\ \textrm{fm}$), as well as different
derivative calculation schemes: Gaussian, chain rule finite difference, and
direct finite difference, for both the original SMASH Skyrme potential, Eq.
(28), and the implemented VDF potential reproducing the SMASH Skyrme
potential, Eq. (30), with Skyrme parameters taken at their usual values of
$a=-209.2\ \textrm{MeV}$, $b=156.4\ \textrm{MeV}$, and $\tau=1.35$. The
difference between the two models, as discussed in Section 3, is that the
Skyrme potential is semi-relativistic, while the VDF potential is fully
relativistic.

Table 1: Values of $\Delta E$ (in MeV/nucleon), $R_{\textrm{MF}}$, and
$t_{\textrm{sim}}$ (in minutes) for uniform nuclear matter simulated in a box
with periodic boundary conditions, initialized in the spinodal region of the
nuclear phase transition, as calculated using the Skyrme potential for
different smearing schemes (“Sm.”, with abbreviations for Gaussian (“G.”),
triangular (“tr.”), and discrete (“dis.”) smearing), values of smearing
parameters (“par.”), numbers of test particles (“$N_{T}$”), and derivative
modes (“$\bm{\nabla}$”). Gaussian smearing always used a width of $\sigma=1\
\textrm{fm}/c$. Uniform matter in a box, Skyrme potential: $T=1\
\textrm{MeV}$, $n_{B}=0.25\ n_{0}$, $t_{\textrm{end}}=200\ \textrm{fm}$.
---
Sm. | par. | $N_{T}$ | Gaussian $\bm{\nabla}$ | finite $\bm{\nabla}$ (chain rule)
$\Delta E$ | $R_{\textrm{MF}}$ | $t_{\textrm{sim}}$ | $\Delta E$ | $R_{\textrm{MF}}$ | $t_{\textrm{sim}}$
G. | $r_{\textrm{c}}=4$ | 20 | $2.40\pm 0.21$ | $1.36\pm 0.04$ | 0.74 | $1.96\pm 0.21$ | $1.34\pm 0.02$ | 0.53
G. | $r_{\textrm{c}}=2$ | 20 | $1.30\pm 0.24$ | $1.15\pm 0.02$ | 0.10 | $2.39\pm 0.23$ | $1.39\pm 0.03$ | 0.08
G. | $r_{\textrm{c}}=4$ | 200 | $1.85\pm 0.11$ | $1.39\pm 0.01$ | 7.22 | $1.25\pm 0.09$ | $1.37\pm 0.02$ | 4.68
G. | $r_{\textrm{c}}=2$ | 200 | $-0.04\pm 0.07$ | $1.08\pm 0.02$ | 1.11 | $1.41\pm 0.12$ | $1.45\pm 0.02$ | 0.78
tr. | $n=2$ | 200 | — | — | — | $1.40\pm 0.08$ | $1.43\pm 0.02$ | 0.54
tr. | $n=4$ | 200 | — | — | — | $0.07\pm 0.04$ | $1.01\pm 0.00$ | 2.60
dis. | $F_{\textrm{c}}=\frac{1}{3}$ | 200 | — | — | — | $3.66\pm 0.20$ | $1.43\pm 0.02$ | 0.21
dis. | $F_{\textrm{c}}=\frac{1}{3}$ | 1000 | — | — | — | $2.02\pm 0.19$ | $1.51\pm 0.02$ | 2.96
Table 2: Values of $\Delta E$ (in MeV/nucleon), $R_{\textrm{MF}}$, and
$t_{\textrm{sim}}$ (in minutes) for uniform nuclear matter simulated in a box
with periodic boundary conditions, initialized in the spinodal region of the
nuclear phase transition, as calculated using the VDF potential for different
smearing schemes (“Sm.”, with abbreviations for Gaussian (“G.”), triangular
(“tr.”), and discrete (“dis.”) smearing), values of smearing parameters
(“par.”), numbers of test particles (“$N_{T}$”), and derivative modes
(“$\bm{\nabla}$”). Gaussian smearing always used a width of $\sigma=1\
\textrm{fm}/c$. Uniform matter in a box, VDF potential: $T=1\ \textrm{MeV}$,
$n_{B}=0.25\ n_{0}$, $t_{\textrm{end}}=200\ \textrm{fm}$.
---
Sm. | par. | $N_{T}$ | Gaussian $\bm{\nabla}$ | finite $\bm{\nabla}$ (chain rule)
$\Delta E$ | $R_{\textrm{MF}}$ | $t_{\textrm{sim}}$ | $\Delta E$ | $R_{\textrm{MF}}$ | $t_{\textrm{sim}}$
G. | $r_{\textrm{c}}=4$ | 20 | $2.26\pm 0.17$ | $1.35\pm 0.02$ | 0.72 | $1.88\pm 0.22$ | $1.35\pm 0.03$ | 0.46
G. | $r_{\textrm{c}}=2$ | 20 | $1.29\pm 0.34$ | $1.14\pm 0.03$ | 0.11 | $2.47\pm 0.21$ | $1.39\pm 0.03$ | 0.08
G. | $r_{\textrm{c}}=4$ | 200 | $1.73\pm 0.12$ | $1.39\pm 0.01$ | 7.14 | $1.24\pm 0.09$ | $1.36\pm 0.02$ | 4.48
G. | $r_{\textrm{c}}=2$ | 200 | $-0.02\pm 0.09$ | $1.08\pm 0.02$ | 1.11 | $1.38\pm 0.13$ | $1.45\pm 0.02$ | 0.79
tr. | $n=2$ | 200 | — | — | — | $1.38\pm 0.09$ | $1.43\pm 0.02$ | 0.50
tr. | $n=4$ | 200 | — | — | — | $0.07\pm 0.04$ | $1.01\pm 0.00$ | 2.46
dis. | $F_{\textrm{c}}=\frac{1}{3}$ | 200 | — | — | — | $4.25\pm 0.18$ | $1.40\pm 0.01$ | 0.23
dis. | $F_{\textrm{c}}=\frac{1}{3}$ | 1000 | — | — | — | $2.11\pm 0.18$ | $1.51\pm 0.02$ | 2.92
Table 3: Values of $\Delta E$ (in MeV/nucleon), $R_{\textrm{MF}}$, and
$t_{\textrm{sim}}$ (in minutes) for uniform nuclear matter simulated in a box
with periodic boundary conditions, initialized in the spinodal region of the
nuclear phase transition, as calculated using the VDF potential for different
smearing schemes (“Sm.”, with abbreviations for Gaussian (“G.”), triangular
(“tr.”), and discrete (“dis.”) smearing), values of smearing parameters
(“par.”), numbers of test particles (“$N_{T}$”), and derivative modes
(“$\bm{\nabla}$”). Gaussian smearing always used a width of $\sigma=1\
\textrm{fm}/c$. Uniform matter in a box, VDF potential: $T=1\ \textrm{MeV}$,
$n_{B}=0.25\ n_{0}$, $t_{\textrm{end}}=200\ \textrm{fm}$.
---
Sm. | par. | $N_{T}$ | finite $\bm{\nabla}$ (chain rule) | finite $\bm{\nabla}$ (direct)
$\Delta E$ | $R_{\textrm{MF}}$ | $t_{\textrm{sim}}$ | $\Delta E$ | $R_{\textrm{MF}}$ | $t_{\textrm{sim}}$
G. | $r_{\textrm{c}}=4$ | 20 | $1.88\pm 0.22$ | $1.35\pm 0.03$ | 0.46 | $1.85\pm 0.25$ | $1.36\pm 0.02$ | 0.46
G. | $r_{\textrm{c}}=2$ | 20 | $2.47\pm 0.21$ | $1.39\pm 0.03$ | 0.08 | $2.24\pm 0.21$ | $1.41\pm 0.04$ | 0.09
G. | $r_{\textrm{c}}=4$ | 200 | $1.24\pm 0.09$ | $1.36\pm 0.02$ | 4.48 | $1.19\pm 0.07$ | $1.38\pm 0.02$ | 4.59
G. | $r_{\textrm{c}}=2$ | 200 | $1.38\pm 0.13$ | $1.45\pm 0.02$ | 0.79 | $1.28\pm 0.08$ | $1.47\pm 0.02$ | 0.78
tr. | $n=2$ | 200 | $1.38\pm 0.09$ | $1.43\pm 0.02$ | 0.50 | $1.29\pm 0.09$ | $1.45\pm 0.01$ | 0.48
tr. | $n=4$ | 200 | $0.07\pm 0.04$ | $1.01\pm 0.00$ | 2.46 | $0.07\pm 0.03$ | $1.01\pm 0.00$ | 2.55
dis. | $F_{\textrm{c}}=\frac{1}{3}$ | 200 | $4.25\pm 0.18$ | $1.40\pm 0.01$ | 0.23 | $4.13\pm 0.24$ | $1.41\pm 0.02$ | 0.23
dis. | $F_{\textrm{c}}=\frac{1}{3}$ | 1000 | $2.11\pm 0.18$ | $1.51\pm 0.02$ | 2.92 | $2.04\pm 0.15$ | $1.52\pm 0.02$ | 2.90

The direct finite difference derivatives are only available for the VDF model,
and Gaussian derivatives only make sense when used together with the Gaussian
smearing. Additionally, finite difference derivatives properly account for the
time derivative of the baryon current, as already discussed in Section 3.
Table 1 focuses on results using the Skyrme potential, while Tables 2 and 3
show results using the VDF potential (the first of these tables compares the
performance of Gaussian versus finite difference chain rule derivatives, and
the second compares the performance of chain rule versus direct finite
difference derivatives).

For an easier analysis of the simulations results shown in Tables 1, 2, and 3,
we color coded the most important cases. A red font marks simulation schemes
which result in a completely incorrect evolution. For all of these cases, in
each event the system fails to separate at all or the degree of separation is
very low, which can be seen from the low values of $R_{\textrm{MF}}$. This
failure to undergo the spinodal decomposition is connected to a faulty
combination of density and gradient calculation schemes; for example, in the
case of the triangular smearing with the range parameter $n=4$ (corresponding
to the smearing range equal $4L_{i}=4\ \textrm{fm}$), the intrinsic density
variations due to the finite number of test particles are spread over a
considerable volume, so that fluctuations at any given part of the box never
become large enough to sustain a spinodal breakup. With orange font we denoted
simulation schemes which result in the deviation from energy conservation
larger than $\Delta E=2\ \textrm{MeV}/\textrm{nucleon}$, while green font
denotes simulation schemes where $\Delta E\leq
1.5\textrm{MeV}/\textrm{nucleon}$. The best performing simulations are those
utilizing finite difference derivatives and either the Gaussian smearing
scheme at $N_{T}=200$ and $r_{\textrm{cut}}=\\{2\sigma,4\sigma\\}$ or the
triangular smearing with $N_{T}=200$ and the triangular smearing range
parameter $n=2$. We see that using Gaussian derivatives is strongly disfavored
in most cases, which is understood to be a consequence of not properly
resolving the density gradients in this derivative scheme; indeed, by
construction, Gaussian derivatives are affected by smearing in the same manner
as the density calculation, making the effects of density fluctuations highly
non-local. We also found that among the best performing simulation schemes
mentioned above, characterized by $\Delta E\leq
1.5\textrm{MeV}/\textrm{nucleon}$, the simulation time for the Gaussian
smearing with $r_{\textrm{cut}}=4\sigma$ is significantly longer than in the
other two cases. Additionally, simulations utilizing the VDF model show that
using direct finite difference derivatives leads to a further improvement in
energy conservation.

To conclude our study of violations of the energy conservation in hadronic
transport, we simulated systems initialized at different points of the phase
diagram, where we also varied the simulation time $t_{\textrm{end}}$ and the
time step $\Delta t$, both measured in $\textrm{fm}/c$. In these simulations,
we used a parametrization of the VDF model reproducing the fourth (IV) set of
characteristics listed in Table 1 (the parameters of this potential can be
found in Appendix 18), which describes nuclear matter with two phase
transitions: one corresponding to the ordinary nuclear matter liquid-gas phase
transition, and one corresponding to a possible QGP-like phase transition with
the critical temperature $T_{c}^{(Q)}=100\ \textrm{MeV}$, the critical density
$n_{c}^{(Q)}=3.0n_{0}$, and the boundaries of the spinodal region at $T=0$
given by $\eta_{L}=2.50n_{0}$ and $\eta_{R}=3.32n_{0}$, where $n_{0}$ is the
saturation density of nuclear matter, $n_{0}=0.160\ \textrm{fm}^{-3}$ (this
and similar EOSs are discussed in Chapter 2, while a detailed discussion of
the results of simulations using the chosen EOS is given in Chapter 4). For
the density and gradients calculation, we used triangular smearing with the
triangular smearing parameter $n=2$ and finite difference direct derivatives.
As previously, the systems were initialized in a box of size $10\ \textrm{fm}$
with periodic boundary conditions, the lattice used had 10 nodes per side, and
we averaged over 10 events.

Table 4: Values of $\Delta E$ (in MeV/nucleon), $R_{\textrm{MF}}$, and $t_{\textrm{sim}}$ (in minutes) for uniform nuclear matter simulated in a box with periodic boundary conditions, initialized for a series of baryon densities (“$n_{B}$”, in units of the saturation density $n_{0}=0.160\ \textrm{fm}^{-3}$), temperatures (“$T$”, in MeV), and evolution times (“$t_{\textrm{end}}$”, in $[\textrm{fm}/c]$) while varying the number of test particles (“$N_{T}$”) and the time step (“$\Delta t$”, in $[\textrm{fm}/c]$), as calculated using the VDF potential with two phase transitions, triangular smearing with $n=2$, and direct finite difference derivatives. $n_{B}$ |  $T$ |  $t_{\textrm{end}}$ |  $N_{T}$ |  $\Delta t$ |  $\Delta E$ |  $R_{\textrm{MF}}$ |  $t_{\textrm{sim}}$
---|---|---|---|---|---|---|---
0.25 | 1 | 200 | 200 | 1.0 | $1.496\pm 0.121$ | $1.455\pm 0.012$ | 0.09
0.25 | 1 | 200 | 200 | 0.1 | $1.435\pm 0.121$ | $1.463\pm 0.012$ | 0.53
0.25 | 1 | 200 | 200 | 0.01 | $1.424\pm 0.117$ | $1.458\pm 0.011$ | 4.71
0.25 | 1 | 400 | 200 | 0.1 | $1.784\pm 0.114$ | $1.465\pm 0.014$ | 1.04
0.25 | 1 | 200 | 500 | 0.1 | $1.295\pm 0.104$ | $1.478\pm 0.014$ | 1.64
0.25 | 25 | 200 | 200 | 0.1 | $0.138\pm 0.023$ | $1.002\pm 0.001$ | 0.54
1.00 | 1 | 200 | 200 | 0.1 | $0.012\pm 0.001$ | $1.0\pm 1.2\times 10^{-5}$ | 3.27
3.00 | 1 | 50 | 50 | 0.1 | $1.903\pm 0.087$ | $1.054\pm 0.004$ | 0.59
3.00 | 50 | 50 | 50 | 0.1 | $0.714\pm 0.040$ | $1.031\pm 0.004$ | 1.06
3.00 | 125 | 50 | 50 | 0.1 | $0.324\pm 0.042$ | $1.004\pm 0.001$ | 0.82

A representative collection of results of these simulations is gathered in
Table 4. First, in the case of nuclear matter initialized in the spinodal
region of the nuclear phase transition, at $n_{B}=0.25n_{0}$ and $T=1\
\textrm{MeV}$, we see that the degree to which the energy is not conserved
depends on the time step $\Delta t$, with smaller unphysical contributions to
the energy for smaller values of $\Delta t$. Consequently, the violation of
energy conservation is seen to occur partially due to errors stemming from the
numerical integration of the equations of motion. However, very small time
steps $\Delta t$ lead to significant increases in the simulation time,
therefore obtaining greater accuracy must be balanced with practical
considerations.

We also find that doubling the duration of a simulation $t_{\textrm{end}}$
does not lead to doubling of the unphysical energy gain $\Delta E$. This
underscores the fact that the rapid non-equilibrium evolution of the system
undergoing a spontaneous spinodal decomposition leads to the most significant
contributions to the violation of energy conservation. In other words, once
the spinodal decomposition has taken place (which for the case of the nuclear
spinodal decomposition at $T=1\ \textrm{MeV}$, considered here, occurs around
$t\approx 100\ \textrm{fm}/c$), the evolution of the system is shown to be
more smooth.

Next, we consider the case of nuclear matter initialized deep inside the
spinodal region of the QGP-like phase transition, at $n_{B}=3.0n_{0}$ and
$T=1\ \textrm{MeV}$. Similarly to the spontaneous spinodal decomposition in
ordinary nuclear matter, we can see that the unphysical energy gain $\Delta E$
is substantial. (We note that the relatively small values of $R_{\textrm{MF}}$
obtained in this case do not indicate an incorrect evolution, unlike in the
example of the spinodal decomposition in ordinary nuclear matter. This can be
understood as follows: After the nuclear spinodal fragmentation occurs, nearly
all test particles can be found in the region where the nuclear liquid drop
forms, and consequently the change in the mean-field energy of the system
roughly corresponds to the difference between the mean-field energy at
$n_{B}=0.25n_{0}$ and the mean-field energy at $n_{B}=n_{0}$. In contrast,
during the spinodal decomposition due to the QGP-like phase transition, the
system separates into a “less dense” ($n_{B}\approx 2.13n_{0}$) and a “more
dense” ($n_{B}\approx 3.57n_{0}$) phase, and the corresponding changes in the
mean-field energies of these parts of the system are opposite in sign.
Consequently, the total change in the mean-field energy is relatively smaller
in the latter case.) At the same time, $\Delta E$ decreases appreciably for an
initialization at $n_{B}=3.0n_{0}$ and $T=50\ \textrm{MeV}$; based on the
phase diagram of the chosen VDF EOS (Fig. 6), we can see that at this higher
temperature the widths of the spinodal and coexistence regions are smaller
than at $T=1\ \textrm{MeV}$, which likely leads to less steep density
gradients developing in the system and through that prevents large unphysical
energy gains. Furthermore, the results obtained for systems initialized above
the critical point of either of the transitions ($T=25\ \textrm{MeV}$ and
$T=125\ \textrm{MeV}$ for the nuclear and QGP-like phase transition,
respectively) show values of $\Delta E$ that are fractions of the values
obtained at $T=1\ \textrm{MeV}$. This is similarly understood to be a
consequence of the fact that at high temperatures, the thermal noise prevents
the system from forming any sustained density gradients. Note that by
construction, the VDF potential does not depend on the temperature; therefore,
the differences observed in the deviations from energy conservation at
different temperatures can be entirely assigned to the presence or lack of
non-uniform density regions in the simulated systems.

Finally, we can easily see that $\Delta E$ is negligible for the
initialization at the equilibrium point of nuclear matter, that is at
$n_{B}=n_{0}$ and $T=1\ \textrm{MeV}$. This is, naturally, connected with the
fact that at the equilibrium point of nuclear matter the forces acting on the
system are vanishing and, at the same time, any fluctuations away from
equilibrium decay fast, leading to only minuscule violation of the
conservation of energy.

In general, we conclude that contributions to $\Delta E$ are large for systems
initialized in regions of the phase diagram where the density gradients, and
therefore forces acting on the test particles, are also large, e.g., inside
the spinodal region of a phase transition. Conversely, energy conservation is
very satisfactory when forces acting on the test particles are small, e.g., in
regions of the phase diagram where nuclear matter is thermodynamically stable.

### 5\. Summary

In this chapter, we have given a brief overview of hadronic transport
simulations, including the description of the method of test particles, and we
have also described some of the computational details specific to simulations
performed with the hadronic transport code SMASH. In particular, we focused on
the aspects of the simulation connected with the calculation of the baryon
density, and we described the density calculation schemes available in SMASH,
many of which were introduced as part of this work. Next, we described the
types of mean-field potentials available in SMASH, including the
implementation of the fully relativistic VDF model potentials, which was part
of the work underlying this thesis. Finally, we provided a comparison of
different numerical schemes for calculating mean fields and mean-field
evolution in SMASH, based on the degree to which the integration of the
equations of motion satisfies energy conservation. We demonstrated that both
the implemented VDF potentials as well as the implemented density and mean-
field calculation schemes perform better than the original options available
in SMASH; in addition, these schemes are shown to be computationally more
efficient. We concluded by discussing the numerical performance of simulations
initialized at different points of the phase diagram.

##  Chapter 4 Results from hadronic transport simulations

The contents of this chapter are largely based on Ref. [177].

We investigate simulations of nuclear matter in SMASH [38] (version 1.8 [206])
realized in a box with periodic boundary conditions. Such studies are best
suited for testing the thermodynamic behavior following from equations of
motion with mean-field interactions, as well as for exploring observables
sensitive to critical phenomena in a scenario in which matter is allowed to
equilibrate. While admittedly systems considered here cannot be reproduced in
the laboratory, insights gained in this study will provide a useful stepping
stone to understanding results of simulations of heavy-ion collisions
utilizing the VDF EOS, planned for future work.

### 1\. Analysis

In contrast to heavy-ion collision experiments, semi-classical hadronic
transport simulations have an access to the positions of individual particles.
Consequently, observables that can be used as a measure of the collective
behavior of the system include the spatial pair correlation function and the
distribution of particles in coordinate space. Below, we describe the details
of extracting these observables.

####  1 Pair distribution function

The radial distribution function $g(r)$ gives the probability of finding a
particle at a distance $r$ from a reference particle. While in select simple
cases it can be calculated analytically, in practice, given a distribution of
particles, $g(r)$ is obtained by determining the distance between the
reference particle and all other particles and constructing a corresponding
histogram. Thus for finding the radial distribution about the $i$th
(reference) particle at a given distance $r$, we count all particles within an
interval $\Delta r$ around $r$, which can be written as

$\displaystyle g_{i}(r,\Delta r)=\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}\theta\Big{(}r+\Delta
r-\mathcal{R}_{ij}\Big{)}\theta\Big{(}\mathcal{R}_{ij}-(r-\Delta
r)\Big{)}~{}.$ (1)

Here, the sum is performed over all particles (with the exception of the $i$th
particle) which we index by $j$, $N$ is the total number of particles,
$\theta$ is the Heaviside theta function, and
$\mathcal{R}_{ij}=|\bm{r}_{i}-\bm{r}_{j}|$, where $\bm{r}_{i}$ is the position
of the reference particle and $\bm{r}_{j}$ is the position of the $j$th
particle. The role of the Heaviside theta functions is to only allow
contributions from particles whose positions are within a distance
$\mathcal{R}_{ij}\in(r-\Delta r,r+\Delta r)$ from the reference particle. The
obtained histogram is then normalized with respect to an ideal gas, whose
radial distribution histogram is that of completely uncorrelated particles,
$g_{0}(r)\propto n~{}4\pi r^{2}~{}dr$, where $n$ denotes density.

Figure 1: Three particles positioned at the vertices of a unilateral triangle
of side $r_{0}$. The arrows point to particles at a distance $r_{0}$ from a
given reference particle (color coded). The figure helps to interpret the pair
distribution function $\widetilde{g}(r_{0},\Delta r)$ as the number of
distinct pairs of particles separated by a distance $r_{0}\pm\Delta r$ (see
text for details).

We can further define the radial distribution function of all distinct pairs
in the system, or the pair distribution function,

$\displaystyle\widetilde{g}(r,\Delta r)$ $\displaystyle=$
$\displaystyle\mathcal{N}~{}\sum_{i=1}^{N}g_{i}(r,\Delta
r)=\frac{\mathcal{N}}{2}\sum_{i=1}^{N}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}\theta\Big{(}r+\Delta
r-\mathcal{R}_{ij}\Big{)}\theta\Big{(}\mathcal{R}_{ij}-(r-\Delta
r)\Big{)}~{},$ (2)

where the factor of $1/2$ appears to avoid counting any of the particle pairs
twice, and where $\mathcal{N}$ is a normalization factor, so far unspecified
(as already mentioned above, in practice the radial distribution function is
compared to that of an ideal gas, in which case the normalization factors
cancel out). It can be easily seen that $\widetilde{g}(r,\Delta r)$ is equal
to the number of distinct (hence the factor of $1/2$) particle pairs separated
by a distance $r\pm\Delta r$. As a simple example, consider a system composed
of 3 particles whose relative distances are equal,
$\mathcal{R}_{ij}=r_{0}~{}\forall~{}i,j$; naturally, these particles must sit
at the vertices of a unilateral triangle of side $r_{0}$, see Fig. 1. The pair
distribution function $\widetilde{g}(r_{0},\Delta r)$ is equal to half the sum
over the number of particles that each reference particle sees at a distance
$r_{0}\pm\Delta r$ away. Clearly, the first reference particle sees 2
particles at a distance $r_{0}$, and likewise for the second and third
reference particles, so that in this case $\widetilde{g}(r_{0},\Delta
r)=(3\times 2)/2$, that is the number of distinct particle pairs separated by
a distance $r_{0}\pm\Delta r$.

The pair distribution function of an ideal gas, denoted by
$\widetilde{g}_{0}(r)$, is related to $g_{0}(r)$ through
$\widetilde{g}_{0}(r)\simeq(N/2)g_{0}(r)$, where $N$ is the total number of
particles in the system, which stems directly from the fact that the total
number of distinct pairs in the system is equal $N(N-1)/2$. For simulations in
a box with periodic boundary conditions, however, this relationship becomes
more complicated for distances $r>L/2$, where $L$ is the side length of the
box, due to geometry effects (see the discussion below). For this reason and
because in simulations presented in this work we initialize the systems
uniformly, in our analysis we use the $t=0$ histogram as the reference pair
distribution function, $\widetilde{g}_{0}=\widetilde{g}(t=0)$.

We stress that taking the pair distribution function of a uniform system as
the reference ensures that the normalized pair distribution function
$\widetilde{g}/\widetilde{g}_{0}$ is sensitive to density fluctuations in the
system. A prominent example here is the spinodal breakup, where a spontaneous
separation into two coexistent phases with different densities occurs. If the
system is confined to some constant volume $V$, then the average density of
the system is the same before and after the spinodal decomposition takes
place. However, local fluctuations in the number of particles will be visible
in the pair distribution function, as more particle pairs reside inside a high
density region as compared to a low density region.

While the spinodal decomposition is the most obvious example of a situation
where $\widetilde{g}/\widetilde{g}_{0}\neq 1$, the normalized pair
distribution function deviates from unity for any system in which the
interactions between the particles affect their collective behavior. In
particular, at small $r$, the normalized pair distribution function satisfies
$\widetilde{g}/\widetilde{g}_{0}>1$ for correlated particles and
$\widetilde{g}/\widetilde{g}_{0}<1$ for anti-correlated particles (see
Appendix 24 for a detailed derivation), which corresponds to attractive and
repulsive interactions between the particles, respectively. Since the number
of particles and thus the number of pairs is conserved, one sees an opposite
trend at intermediate to large distances.

We note that in our simulations the range of $r$ over which
$\widetilde{g}(r)/\widetilde{g}_{0}(r)$ significantly deviates from 1 is
related to the range of the interaction, which is determined by the smearing
range in the density calculation (for more details see Section 2).

##### Pair distribution function in a system with periodic boundary conditions

Importantly, for a system with periodic boundary conditions the radial
distance between two particles $\mathcal{R}$ is not uniquely defined. This is
because for any reference particle the distance to any other particle can be
calculated using the position of that other particle in the original box or in
any of its 26 equivalent images. We adopt a prescription in which the smallest
distance between particles is used in calculating the pair distribution
function $\widetilde{g}$ (known as the minimum image criterion). This smallest
distance can range from $\mathcal{R}_{\textrm{min}}=0$ to
$\mathcal{R}_{\textrm{max}}={\sqrt{3}L}/{2}$, where $L$ is the side length of
the box.

Figure 2: A 2D diagram of geometry effects affecting the pair distribution
function of a uniform system in a box with periodic boundary conditions. We
draw the simulated box, with the reference particle denoted with a blue
triangle and three other particles denoted with a green square, a red dot, and
a red square, as well as its 8 equivalent images (faded colors). A green
dashed circle of radius $L/2$ and a red dashed circle of radius $L/\sqrt{2}$
are drawn around the reference particle. Distances from the reference particle
to other particles which contribute to the pair distribution function are
implied with gray arrows.

That said, even for a uniform and uncorrelated system the geometry of the
problem affects the number of particles that can be encountered at the maximal
distance $\mathcal{R}_{\textrm{max}}$. Let us show this based a 2D
illustration of this issue, see Fig. 2. To calculate the minimal distance
between any two particles, we consider the original system (in the middle) and
its 8 equivalent images (the latter are indicated by faded colors). We
consider a particle at a position denoted with a blue triangle, and to aid the
analysis we draw a green dashed circle of radius $L/2$ and a red dashed circle
of radius $L/\sqrt{2}$ around it. Consider two other particles, whose
positions are denoted with a red square and a red dot, separated from the blue
triangle by an equal distance $\mathcal{R}$ such that
$L/2<\mathcal{R}<L/\sqrt{2}$. For the particle pair composed of the blue
triangle and the red dot, the line separating the particles lies along the
diagonal of the box; since $\mathcal{R}$ is already the smallest distance
between the blue triangle and any of the equivalent images of the red dot,
this particle pair contributes to the pair distribution function at
$\mathcal{R}$. For the particle pair composed of the blue triangle and the red
square, the line separating the particles lies along one of the principle
directions of the box; in contrast to the previous case, there is an
equivalent image of the red square characterized by a smaller distance to the
blue triangle $\mathcal{R}^{\prime}<\mathcal{R}$; consequently, this pair
contributes to the pair distribution function at $\mathcal{R}^{\prime}$. In
fact, particles pairs lying along the diagonal and along one of the principal
box directions are special cases. Notice, in particular, that the only points
for which it is possible to have $\mathcal{R}=\mathcal{R}_{\textrm{max}}$ are
points on the diagonal of the box; for any points separated by
$\mathcal{R}_{\textrm{max}}$ that are not on the diagonal, there exists a
smaller $\mathcal{R}$ obtained by using the position of the second particle
from one of the equivalent box images. This problem also affects, to a
proportionally lesser extent, interparticle distances $\mathcal{R}$ in the
range ${L}/{2}<\mathcal{R}<\mathcal{R}_{\textrm{max}}$. Only in the case of
particles which are ${L}/{2}$ or less apart the geometry of the box never
affects the pair distribution function, as can be seen from the example of the
particle depicted with a green square.

Figure 3: Pair distribution function of spatially uniform matter in a cubic
box with periodic boundary conditions (red histogram). To obtain the
histogram, a cubic box of side length $L=10\ \textrm{fm}$ was initialized with
proton and neutron numbers $N_{p}=N_{n}=20$ (corresponding to a baryon density
of $0.25n_{0}$), and the number of test particles per particle was set at
$N_{T}=200$, resulting in the test particle density of $n=8\
\textrm{fm}^{-3}$. For $r<L/2$, the pair distribution function is given by the
pair distribution function of a uniform ideal gas,
$\widetilde{g}(r)\simeq(N/2)~{}n~{}4\pi~{}r^{2}~{}\Delta r$ (solid blue line),
where $\Delta r=0.433$ is the histogram bin width. For $r>L/2$, geometric
effects related to the periodic boundary conditions start playing a role; see
text for more details.

This influence of finite size effects on the pair distribution function can be
clearly seen in Fig. 3, which shows the pair distribution function for a
uniform matter in a cubic box of side length $L=10\ \textrm{fm}$ and subject
to periodic boundary conditions. In infinite matter, the pair distribution
function of uncorrelated particles grows like $r^{2}$. However, finite
geometry effects described above introduce an effective cut on the
distribution starting at ${L}/{2}=5\ \textrm{fm}$, explaining the shape of the
presented histogram.

Figure 4: A 2D diagram of the geometry effects affecting the normalized pair
distribution function of a non-uniform system. A clump of particles (blue
circle) surrounded by vacuum is shown in a box of side $L$ (middle box), and
positions of three hypothetical particles are indicated with a green square, a
blue circle, and a red triangle. We indicate a range over which the pair
distribution function is calculated for the particle in the middle of the
clump (blue dot) with a circle of radius $L/\sqrt{2}$ centered around that
particle, and we also show the 8 equivalent images of the box (faded colors).
The smallest equivalent distance between the particles represented by the red
triangle and the green box is realized “across the vacuum”, either from the
red triangle to the faded green square, or from the green square to the faded
red triangle. Such a situation will arise for many particles on the edge of
the clump, leading to a peculiar behavior of the pair distribution function at
distances comparable to $L/\sqrt{2}$. See text for more details.

Naturally, geometry and periodic boundary conditions also play a role in the
shape of the normalized pair distribution function in non-uniform systems. As
an easy illustration, consider a 2D system composed of particles forming a
spherical “clump” of matter surrounded by vacuum, see Fig. 4. It is easy to
see that if the diameter of the clump satisfies $D>{L}/{2}$, then for some of
the particles forming the clump the smallest distance to some of the other
particles in that same clump will be “across the vacuum”, to one of the
equivalent mirror images of the latter particles; in the figure this is
showcased by the example of particle positions represented by the green square
and the red triangle. In result, on average, the normalized pair distribution
function for this system will display a strong correlation for small values of
$r$ (reflecting the existence of the clump), then an anticorrelation for
intermediate values of $r$ (reflecting the vacuum), and finally a significant
correlation for values of $r$ approaching $\mathcal{R}_{\textrm{max}}$
(reflecting the smallest distances “across the vacuum”).

In our simulations, which we will discuss at length in the following Sections,
we clearly observe the above described effect in the case of nuclear spinodal
decomposition, which yields a nuclear drop surrounded by a nearly perfect
vacuum. The diameter of the nuclear drop turns out to satisfy $D>{L}/{2}$,
which means that for some of the test particles belonging to that drop, the
smallest distance to some of the other test particles in that same drop is
obtained “across the vacuum”. This explains the rise in the normalized
distribution function for $r>{L}/{2}$, which can be seen, e.g., on the right-
hand panel in Fig. 8. The magnitude of this effect depends on the drop
diameter $D$.

Finally, we note that the artifacts produced by the geometry of the simulation
and periodic boundary conditions do not present a significant complication in
analyzing critical behavior if we resolve to only probe the system at length
scales on the order of ${L}/{2}$ or smaller.

##### Pair distribution function and Boltzmann transport

One may ask whether calculating a pair distribution function in hadronic
transport is justified in view of the fact that the BUU equation explicitly
evolves a one-body distribution function which does not carry any information
about the two-body distribution, usually employed in the description of two-
particle correlations (see Appendix 23.E for a derivation of the Boltzmann
equation from the BBGKY hierarchy, where the relevant steps leading to
neglecting the two-body contributions are shown). While this may appear to be
problematic, a closer look reveals that our analysis is correct. First, one
needs to note that hadronic transport simulations only solve the Boltzmann
equation exactly in the limit of an infinite number of test particles per
particle $N_{T}$. The finite number of test particles employed in simulations
leads to intrinsic numerical fluctuations. These numerical fluctuations are of
statistical nature, similarly to variances of microscopic observables, and
likewise, through both scattering and mean fields, they can become a seed for
collective behavior such as spontaneous spinodal decomposition. Such effects
have been described, e.g., in Ref. [207] (see also Refs. [208, 209]), where
fluctuation observables calculated using hadronic transport with the method of
test particles agree with both theoretical predictions and experimental
results. Additionally, it was established that for large enough $N_{T}$ (which
the authors of that particular study found to be $N_{T}\gtrsim 40$) the
numerical noise intrinsic to the method of test particles is negligible, while
the correct statistical fluctuations are preserved.

It is possible to construct a Boltzmann-Langevin extension of the standard BUU
equation, which ensures that the simulated fluctuations are physically correct
(see, e.g., [210]). However, it has been found that, for example, in the case
of the nuclear spinodal fragmentation the source of the noise seeding the
spinodal decomposition is not essential, and it is possible to develop good
approximations to the Boltzmann-Langevin equation that are numerically
favorable, including the method of test particles [187].

We note here that a particular problem that arises in the method of test
particles is that the fluctuations in the events, simulating the evolution of
$N_{T}N_{B}$ test particles, are suppressed by a factor of $N_{T}$. The
authors of [207] dealt with this issue by employing the method of parallel
ensembles at final simulation times, that is a posteriori, which allows one to
obtain events with the number of test particles corresponding to the physical
baryon number $N_{B}$ (we briefly describe this method in Section 1, while
Appendix 25 explains the a posteriori application of the method).

Based on the above it is apparent that the distribution function obtained
through hadronic transport simulations, and in particular through the method
of test particles, contains information not only about the mean of the
distribution function $\langle f_{\bm{p}}\rangle$, but also about its
fluctuations. Consequently, calculating fluctuation observables such as the
pair distribution function is well-defined in hadronic transport. Some
questions regarding the quantitative behavior of fluctuation observables
obtained in simulations using the number of test particles $N_{T}>1$ remain,
in particular regarding the specific methods used to connect fluctuations in
systems evolving $N_{T}N_{B}$ particles as compared to systems evolving
$N_{B}$ particles. For that reason we refrain from making quantitative
statements at this time, and focus on the qualitative behavior of the pair
distribution functions. Future work will be devoted to a quantitative analysis
of this problem, and in Section 3 we give a short overview of the effects due
to this issue.

####  2 Number distribution functions

A complementary method of analyzing the collective behavior in a simulation
utilizes coordinate space number distribution functions. To calculate number
distribution functions, we divide the simulation box into $C$ cells of side
length $\Delta l$ (also referred to as cell width), and construct a histogram
of the number of cells in which the number of particles lies in a given
interval $N_{i}\pm\Delta N$, where $N_{i}$ is the central value of the $i$th
bin. We note that we scale the entries by the total number of cells $C$ so
that the resulting histogram is a properly normalized representation of the
corresponding probability distribution. We also note that in the subsequent
parts of this work we scale the histogram entries by the volume of the cells
$(\Delta l)^{3}$ in order to obtain the histogram as a function of number
density.

The test-particle evolution in SMASH is governed by the mean field, which
depends on the underlying continuous baryon number density for a given baryon
number $N_{B}$,

$\displaystyle
n_{B}(\bm{x};N_{B})=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}f(\bm{x},\bm{p})~{}.$ (3)

Formally, hadronic transport can give access to $n_{B}(\bm{x};N_{B})$ through
solving the Boltzmann equation, Eq. (3), in the limit of infinitely many test
particles per particle, and substituting the obtained distribution function
$f(\bm{x},\bm{p})$ in Eq. (3). Below, we present three number distribution
functions accessible in practice given the finite number of test particles
used.

##### Test particle number distribution function

Hadronic transport simulations of nuclear matter are realized through evolving
$N=N_{B}N_{T}$ test particles in space and time (where $N_{B}$ is the baryon
number in the simulation and $N_{T}$ is the number of test particles per
particle), giving a direct access to a discrete test particle number
distribution function. This distribution can be written as a probability of
obtaining a cell contributing to the $i$th bin of the histogram with a center
value $N_{i}$ (that is, a cell with $N\in(N_{i}-\Delta N,N_{i}+\Delta N)$ test
particles),

$\displaystyle P_{N}(N_{i})$ $\displaystyle=$ $\displaystyle
P\Big{(}N_{i},N(N_{B},N_{T}),\Delta l\Big{)}=$ (4) $\displaystyle=$
$\displaystyle\frac{N^{(i)}_{c}\Big{(}N(N_{B},N_{T}),\Delta l\Big{)}}{C}~{}.$
(5)

Here, $C$ is the total number of cells used and $N^{(i)}_{c}$ is the number of
cells containing a number of test particles $N$ within the range
$N_{i}\pm\Delta N$. We note that the number of test particles in any given
cell depends both on the baryon number evolved in the simulation $N_{B}$ and
the number of test particles per particle $N_{T}$. We also stress that the
distribution $P_{N}$ depends on the scale (chosen cell width $\Delta l$) at
which the system is analyzed.

##### Continuous baryon number distribution function

The discrete test particle distribution function, Eq. (4), can be thought of
as having been obtained through sampling from the underlying continuous baryon
number distribution function with a finite number $N_{T}N_{B}$ of test
particles. Given access to the underlying baryon number distribution, one
could use it directly to create a corresponding histogram. Indeed, the number
of baryons at a cell at position $\bm{x}_{k}$ is given by the integral of the
continuous baryon number density, Eq. (3),

$\displaystyle B(\bm{x}_{k})=\int_{V_{k}=(\Delta
l)^{3}}dV~{}n_{B}(\bm{x},N_{B})~{},$ (6)

where $k$ indexes the histogram cells. Adding contributions from all cells
yields the total baryon number in the system, $B$. We can then construct a
probability distribution function for encountering a cell with a given number
of baryons $N_{i}$,

$\displaystyle P_{B}(N_{i})=\frac{N_{c}^{(i)}\Big{(}N_{i},B,\Delta
l\Big{)}}{C}~{},$ (7)

where $N_{c}^{(i)}$ is the number of cells containing a number of baryons $N$
within the range $N_{i}\pm\Delta N$.

For a large number of test particles per particle $N_{T}$, statistical
observables calculated using the test particle number distribution, with the
number of test particles in a given sample scaled by ${1}/{N_{T}}$, are a very
good approximation to the underlying continuous baryon number distribution
[211]. That is, it can be shown that

$\displaystyle
P_{B}(N_{i})=\lim_{N_{T}\to\infty}P\bigg{(}N_{i},\frac{N(N_{B},N_{T})}{N_{T}},\Delta
l\bigg{)}~{}.$ (8)

Given that in our simulations we use sufficiently large numbers of test
particles per particle $N_{T}$, we will refer to histograms constructed
through the prescription on the right-hand side of Eq. (8) as the continuous
baryon number distribution function (or just baryon number distribution
function) $P_{B}(N_{i})$, with the understanding that it is only exact in the
limit $N_{T}\to\infty$.

##### Physical baryon number distribution function

Both the test particle and the continuous baryon number distribution
functions, Eqs. (4) and (8), are markedly different from the physical baryon
number distribution function corresponding to a discrete baryon number
$N_{B}$. Here we can intuitively think of the physical baryon number
distribution function as obtained through sampling from the underlying
continuous baryon number distribution with $N_{B}$ test particles,

$\displaystyle P_{N_{B}}(N_{i})=P\Big{(}N_{i},N(N_{B},N_{T}=1),\Delta
l\Big{)}~{}.$ (9)

Strictly speaking, the physical baryon number distribution function could be
obtained in transport by solving the Boltzmann equation in the limit of
infinitely many test particles per particle, thus obtaining the underlying
continuous baryon number distribution function, Eq. (3), and sampling
$n_{B}(\bm{x},N_{B})$ with $N_{B}$ particles. While such an approach is in
principle valid, in practice it is impractical due to the enormous numerical
coast necessary to approximate the limit of infinitely many test particles
well. Alternatively, one can turn to the concept of parallel ensembles
(introduced in Section 1). It can be shown that the test particle distribution
obtained within a parallel ensembles mode can serve as a proxy for the
physical baryon number distribution. To reiterate, within the concept of
parallel ensembles, a simulation corresponding to $N_{B}$ baryons with $N_{T}$
test particles per baryon is divided into $N_{T}$ events with $N_{B}$ test
particles each. These $N_{T}$ events are not independent, as they share a
common mean field. Nevertheless, at the end of the simulation we have access
to $N_{T}$ events with the test particle number exactly corresponding to the
baryon number in the “real” system. That is, each of the $N_{T}$ events is
described by the probability distribution function
$P_{N_{B}}(N_{i})=P\Big{(}N_{i};N(N_{B},N_{T}=1);\Delta l\Big{)}$. Observables
calculated using $P_{N_{B}}(N_{i})$ are probably the closest to those one
would find in an experiment if one could measure positions of the particles.

### 2\. Infinite matter simulation results

To simulate isospin-symmetric infinite nuclear matter, we initialize equal
numbers of proton and neutron test particles in a box with periodic boundary
conditions. The side length of the box is taken to be $L=10\ \textrm{fm}$;
this is informed by the fact that with periodic boundary conditions, the box
can be kept relatively small with no significant finite-size effects. The time
step used in the simulation needs to be small enough to resolve all gradients
occurring during the evolution (intuitively speaking, a test particle should
not “jump over” a potential gradient within a single time step). We found that
a time step of $\Delta t=0.1\ \textrm{fm}/c$ is small enough to satisfy this
condition, and it correctly solves the equations of motion, Eqs. (79) and
(80), using the leapfrog algorithm. The mean field is calculated using the
triangular smearing with triangular smearing parameter $n=2$ and finite
difference direct derivatives (for more details see Chapter 3) on a lattice
with lattice spacing $a=1\ \textrm{fm}$, which has been tested to be
sufficiently fine for accurately resolving mean-field gradients. To ensure
smooth density and mean-field gradient calculations, we utilize a large number
of test particles per particle, specifically, we use $N_{T}=200$ for ordinary
nuclear matter (Section 1) and $N_{T}=50$ for dense nuclear matter (Sections 2
and 3). Using different numbers of test particles in these two cases is
justified by the fact that smooth density and density gradient calculations
are ensured when the average number of test particles encountered in a cell of
the lattice $N_{\textrm{avg}}$ is large enough. As an example, within the
described setup, this number will be equal to $N_{\textrm{avg}}=8$ for
ordinary nuclear matter at $n_{B}=0.25n_{0}$, and equal to
$N_{\textrm{avg}}=24$ for dense nuclear matter at $n_{B}=3n_{0}$. We choose
$N_{\textrm{avg}}$ to be bigger in the case of dense nuclear matter as mean
fields encountered in that region of the phase diagram are significantly
larger and require an even more smooth gradient computation.

For studying the thermodynamic behavior of nuclear matter, we are simulating
systems in which all collision and decay channels are turned off. We have
checked that the thermodynamic effects described here persist when collisions
are allowed, and in this work we choose to omit them because our goal is to
study mean-field dynamics. As in Section 4, we are considering only one of the
many EOSs accessible within the VDF model, namely, the one corresponding to
the fourth (IV) set of characteristics listed in Table 1. The choice of this
set is arbitrary and serves as an illustration of the properties of the VDF
model which are qualitatively comparable for all obtained EOSs.

####  1 Ordinary nuclear matter

We investigate the behavior of systems initialized at temperatures and baryon
number densities specific to ordinary nuclear matter to validate the
implementation of the VDF model in SMASH [38]. For illustrative purposes, we
discuss results for a single simulation run, that is one event. Remarkably,
the thermodynamic behavior of the system is apparent already for this minimal
statistics. This is a consequence of the large number of test particles per
particle used ($N_{T}=200$), as well as the fact that the investigated effects
are characterized by large fluctuations, which result in clear signals.

Figure 5: Baryon number distribution, scaled by the volume of the cell and
shown in units of the saturation density, $n_{0}=0.160\ \textrm{fm}^{-3}$. The
cell width is $\Delta l=2\ \textrm{fm}$. Red histograms show distributions at
initialization ($t=0$), while blue shaded histograms show distributions at the
end of the evolution, $t=t_{\textrm{end}}$. Left panel: Nuclear matter
initialized at the saturation density $n_{0}$ and temperature $T=1\
\textrm{MeV}$, evolved until $t_{\textrm{end}}=200\ \textrm{fm}/c$. The
system, initialized in a stable configuration, remains in the same state at
$t_{\textrm{end}}$. Right panel: Nuclear matter initialized in the spinodal
region of the nuclear phase transition, at baryon density $n_{B}=0.25n_{0}$
and temperature $T=1\ \textrm{MeV}$, evolved until $t_{\textrm{end}}=100\
\textrm{fm}/c$. The system, initialized in a mechanically unstable region of
the phase diagram, spontaneously separates into a nucleon gas and a nuclear
liquid drop with central density $n_{B}\approx n_{0}$. Figure from [177].

To start, we initialize symmetric nuclear matter at saturation density
$n_{B}=n_{0}$, which for the box setup described above corresponds to the
number of protons and neutrons $N_{p}=N_{n}=80$, and at temperature $T=1\
\textrm{MeV}$. Except for a slight increase in temperature from the degenerate
limit, which is not significant enough to introduce any relevant changes, this
is the saturation point of nuclear matter, where the configuration is stable.
We let the simulation evolve until $t_{\textrm{end}}=200\ \textrm{fm}/c$ and
investigate whether the equilibrium is preserved by hadronic transport.

Figure 6: Time evolution of kinetic energy per particle (left panel), mean-
field energy per particle (middle panel), and total (binding) energy per
particle (right panel) for a system initialized at nuclear saturation density
$n_{0}=0.160\ \textrm{fm}^{-3}$ and temperature $T=1\ \textrm{MeV}$. The
binding energy per particle at initialization, $E_{B}(t=0)\approx-16.218\
\textrm{MeV}$, is within 0.1% from the value expected from model calculations.
The mean-field energy oscillates slightly throughout the evolution, reflecting
local fluctuations in density, but its average value remains constant. The
increase in the kinetic energy per particle in time, which also causes the
increase in the total energy per particle in time, is an unwanted feature of
the simulation. A slight violation of the conservation of energy is a common
feature of many hadronic transport codes, and is connected to the choice of
the integration method for the equations of motion as well as to details of
density and density gradient calculations (see Section 4 for more details).
Figure from [177].

To address this question, we examine the continuous baryon number distribution
function (for details, see Section 2), which we calculate using the cell width
$\Delta l=2\ \textrm{fm}$; we scale the histogram entries by the volume of the
cell to obtain the distribution in units of the baryon number density, and
further scale the results to express them in units of the saturation density,
$n_{0}=0.160\ \textrm{fm}^{-3}$. As expected for matter in equilibrium, the
baryon number distribution remains unchanged throughout the evolution, as can
be seen in the left panel of Fig. 5. Note that the finite width of the shown
distribution is a direct consequence of using a finite number of test
particles; we will address this issue and related effects in Section 3. We
find that throughout the simulation, the binding energy per particle agrees
with the theoretically obtained value within $0.1\%$, see Fig. 6 (see Section
4 for more details on energy conservation in hadronic transport). An in-depth
discussion of the mean-field response to fluctuations around nuclear
saturation density, comparing the results from several transport codes
including SMASH utilizing the VDF model, can be found in Ref. [202].

Figure 7: Visualization of spinodal decomposition in a system initialized
inside the nuclear spinodal region (at temperature $T=1\ \textrm{MeV}$ and
baryon number density $n_{B}=0.25n_{0}$, where $n_{0}$ is the saturation
density, $n_{0}=0.160\ \textrm{fm}^{-3}$). In the left panel one sees uniform
nuclear matter at initialization, in the middle panel the mechanically
unstable system is “collapsing” around a local random density fluctuation, and
in the right panel the system has reached equilibrium where matter consists of
a dense “nuclear drop” surrounded by a very dilute gas of nucleons. The
effects due to periodic boundary conditions are clearly visible. We note that
the black dots do not represent individual test particles, but rather are
indicative of relative test particle density in a given region of the box.

Next, we study nuclear matter inside the spinodal region of the nuclear phase
transition. Specifically, we initialize the system with the number of protons
and neutrons $N_{p}=N_{n}=20$, corresponding to a baryon number density
$n_{B}=0.25n_{0}$, at temperature $T=1\ \textrm{MeV}$. We let the system
evolve until $t_{\textrm{end}}=100\ \textrm{fm}/c$. The spinodal region is
both thermodynamically and mechanically unstable, and so we expect that local
density fluctuations will drive the matter to separate into two coexisting
phases: a dense phase, also known as a nuclear drop, and a dilute phase which
is a nucleon gas. That this indeed happens can be seen on the right panel in
Fig. 5, which shows the change in the baryon number distribution function due
to the system’s separation into two coexisting phases. The distribution,
initially centered at $n_{B}=0.25n_{0}$, at the end of the evolution has a
large contribution at $n_{B}\approx 0$, corresponding to the majority of box
cells that after the separation are empty or nearly empty, and a long tail
reaching out to $n_{B}\approx n_{0}$, which corresponds to the center of the
nuclear drop. Even more directly, this separation can be seen in a
visualization of the decomposition process, shown at three chosen evolution
times, $t=0,50,100\ \textrm{fm}/c$, in Fig. 7.

Figure 8: Time evolution of the pair distribution function for a system
initialized inside the nuclear spinodal region (at temperature $T=1\
\textrm{MeV}$ and baryon number density $n_{B}=0.25n_{0}$, where $n_{0}$ is
the saturation density, $n_{0}=0.160\ \textrm{fm}^{-3}$). The $t=0$ plot
(left) shows the distribution at initialization, while plots at $t=50$ and
$100\ \textrm{fm}/c$ (middle and right, respectively) show normalized
distributions. Spontaneous spinodal decomposition occurs at $t>0$ and leads to
a formation of a nuclear drop surrounded by a near-perfect vacuum, resulting
in a strong correlation between particles clustered within the drop. See
Section 1 for a discussion of the influence of finite-size effects on the
shape and large-distance behavior of the pair distribution function. Figure
from [177].

We then proceed to calculate the pair distribution function (for details, see
Section 1) for the system initialized in the spinodal region of nuclear
matter. The results are shown in Fig. 8. Here, the three panels correspond to
three time slices of the evolution: $t=0,50,100\ \textrm{fm}/c$. The $t=0$
plot (left) shows the pair distribution function, Eq. (2), at initialization
$\widetilde{g}_{0}(r,\Delta r)$, while plots at $t=50$ and $100\
\textrm{fm}/c$ (middle and right, respectively) show normalized pair
distribution functions $\widetilde{g}(r,\Delta r)/\widetilde{g}_{0}(r,\Delta
r)$. The time evolution of the pair distribution function shows that during
the spinodal decomposition the test particles cluster into the nuclear drop.
The half width at half maximum of the pair distribution function is about $2\
\textrm{fm}$, which corresponds to the density smearing range used, see
Section 1 as well as Section 2 for more details. The influence of the periodic
boundary conditions on the shape and behavior of the pair distribution
function at large inter-particle distances is also discussed in Section 1.
Additionally, the number of drops that form during spinodal decomposition
depends both on the size of the box and the size of a drop, where the latter
again depends on the smearing range used in density calculation.

Altogether, the results presented above demonstrate that the VDF equations of
motion implemented in SMASH reproduce the expected bulk behavior of ordinary
nuclear matter.

####  2 Dense nuclear matter and the QGP-like phase transition

For simulations of critical behavior in dense symmetric nuclear matter, we run
$N_{\textrm{ev}}=500$ events and average the results, calculated event-by-
event. We first initialize the system at $n_{B}=3n_{0}$, which corresponds to
the number of protons and neutrons $N_{p}=N_{n}=240$, and at temperature $T=1\
\textrm{MeV}$. It can be seen in Figs. 6 and 8-10 that this corresponds to
initializing dense nuclear matter inside the spinodal region of the QGP-like
phase transition described by the EOS employed (the fourth (IV) set of
characteristics listed in Table 1). We evolve the system until
$t_{\textrm{end}}=50\ \textrm{fm}/c$, which is sufficient for reaching
equilibrium after a spinodal decomposition at high baryon number densities,
since due to considerably larger values of the mean-field forces on test
particles the density instabilities develop more rapidly.

Figure 9: Time evolution of the baryon number distribution, scaled by the
volume of the cell and shown in units of the saturation density of nuclear
matter, $n_{0}=0.160\ \textrm{fm}^{-3}$, for a system initialized inside the
quark-hadron spinodal region (at baryon number density $n_{B}=3n_{0}$ and
temperature $T=1\ \textrm{MeV}$), averaged over $N_{\textrm{ev}}=500$ events.
The cell width is chosen at $\Delta l=2\ \textrm{fm}$. Histograms delineated
with red curves correspond to the baryon distribution at initialization
($t=0$), while histograms delineated and shaded with blue curves correspond to
baryon distributions at a chosen time during the evolution ($t=\\{25,50\\}\
\textrm{fm}/c$). The system, initialized in a mechanically unstable region of
the phase diagram, undergoes a spontaneous separation into a less dense and a
more dense nuclear liquid (see Section 1 for more discussion), resulting in a
double-peaked baryon number distribution. The green arrows point to values of
baryon number densities corresponding to the boundaries of the coexistence
region at $T=1\ \textrm{MeV}$, $n_{L}=2.13n_{0}$ and $n_{R}=3.57n_{0}$. Figure
from [177].

Figure 10: Time evolution of the pair distribution function for a system
initialized inside the QGP-like spinodal region (at baryon number density
$n_{B}=3.0n_{0}$ and temperature $T=1\ \textrm{MeV}$), averaged over
$N_{\textrm{ev}}=500$ events. The $t=0$ plot (first panel) shows the pair
distribution at initialization, while plots at $t=10,\ 20,\ 30\ \textrm{fm}/c$
(second, third, and fourth panels) show normalized pair distributions.
Spontaneous spinodal decomposition occurs at $t>0$ and leads to a formation of
two coexisting phases: a less dense and a more dense nuclear liquid. The
increased relative concentration of particles in the more dense phase results
in an elevated normalized pair distribution at small distances. Figure from
[177].

In Fig. 9, we show the evolution of the baryon number distribution (see
Section 2). The cell width is chosen at $\Delta l=2\ \textrm{fm}$, and the
histogram entries are scaled by the volume of the cell in order to be given in
units of the baryon number density; we then further scale the results to
express them in units of the saturation density, $n_{0}=0.160\
\textrm{fm}^{-3}$. In the figure, the red curve corresponds to the
distribution at time $t=0$, while the blue curves delineate the distribution
at times $t>0$. At $t=0$, the distribution is peaked at the initialization
density $n_{B}=3n_{0}$, with its width reflecting the finite number
statistics. In the course of the evolution the system separates into two
coexisting phases, a “less dense” and a “more dense” nuclear liquid (see
Section 1 for more discussion). As a result, the baryon distribution displays
two peaks largely coinciding with the theoretical values of the coexistence
region boundaries, $n_{L}=2.13n_{0}$ and $n_{R}=3.57n_{0}$, indicated by the
green arrows. We find that the prominence of the peaks depends slightly on the
choice of the EOS. For example, an equation of state with the same value of
critical density $n_{c}^{(Q)}$ and the same spinodal region
$(\eta_{L},\eta_{R})$, but a higher critical temperature $T_{c}^{(Q)}$ will
correspond to a more negative slope of the pressure in the spinodal region
and, correspondingly, to stronger mean-field forces inside the spinodal
region, leading to more prominent peaks.

Next, in Fig. 10 we show the evolution of the pair distribution function.
Similarly as in the case of nuclear spinodal decomposition, the “hadron-quark”
spinodal decomposition leads to a pair distribution function indicating the
formation of two phases of different densities. Unlike in nuclear spinodal
decomposition, where drops of a “nuclear liquid” form in vacuum, in this case
we have drops of a “more dense liquid” submerged in a “less dense liquid” (for
a detailed discussion, see Section 1). Consequently, the absolute values of
the normalized pair distribution function
$\widetilde{g}(r)/\widetilde{g}_{0}(r)$ are much smaller for the case of the
“hadron-quark” spinodal decomposition, as the difference between the number of
test particle pairs occupying the dense and dilute regions is less pronounced
in this case. Nevertheless, the effect, although small, is clearly
distinguishable and statistically significant.

We note here that a phase separation is such a distinct behavior of the system
that the baryon distribution function and the pair distribution function as
shown in Figs. 9 and 10, respectively, can be largely recovered even in the
case of minimal statistics, that is for one event. However, effects at and
around the critical point, as discussed below, are much more subtle and
require a relatively large number of events.

To conclude our study of dense nuclear matter in SMASH, we want to investigate
the behavior of systems initialized at various points of the phase diagram
above the critical point, inspired by possible phase diagram trajectories of
heavy-ion collisions at different beam energies. Specifically, we initialize
the system at one chosen temperature and a series of baryon number densities

$\displaystyle T=125~{}\textrm{MeV},\hskip
5.0ptn_{B}\in\\{2.0,2.5,3.0,3.5,4.0\\}n_{0}~{}.$ (10)

In contrast with most of the previous examples, systems initialized in this
region of the phase diagram are thermodynamically stable, and there are
specific predictions for the behavior of thermodynamic observables such as
ratios of cumulants of baryon number (see Figs. 8-10). In the upper panel of
Fig. 11, we show values of the second-order cumulant ratio,
${\kappa_{2}}/{\kappa_{1}}$, as calculated from the VDF model, both in the
$(T,n_{B})$ and the $(T,\mu_{B})$ plane. The dots on the cumulant diagrams
mark the points at which we initialize the system, specified in Eq. (10), and
are intended to guide the eye toward the corresponding normalized pair
distribution plots at the end of the evolution, $t=t_{\textrm{end}}$,
displayed in the lower panel of the same figure. The deviation of values of
the normalized pair distributions at small distances from 1 (where 1
corresponds to a system of non-interacting particles) directly follows the
deviation of values of the second-order cumulant ratio
${\kappa_{2}}/{\kappa_{1}}$ from the Poissonian limit of 1,

$\displaystyle\frac{\widetilde{g}~{}(0,\Delta r)}{\widetilde{g}_{0}(0,\Delta
r)}>1~{}\Leftrightarrow~{}\frac{\kappa_{2}}{\kappa_{1}}>1~{},$ (11)
$\displaystyle\frac{\widetilde{g}~{}(0,\Delta r)}{\widetilde{g}_{0}(0,\Delta
r)}<1~{}\Leftrightarrow~{}\frac{\kappa_{2}}{\kappa_{1}}<1~{}.$ (12)

We show a detailed derivation of this fact in Appendix 24. It is clear that a
two-particle correlation corresponds to a value of the cumulant ratio
${\kappa_{2}}/{\kappa_{1}}>1$, and a two-particle anticorrelation corresponds
to a value of the cumulant ratio ${\kappa_{2}}/{\kappa_{1}}<1$. This behavior
is exactly reflected in Fig. 11.

Figure 11: Comparison of the cumulant ratio $\kappa_{2}/\kappa_{1}$ (upper
panel), calculated in the VDF model, and the normalized pair distribution
function $\widetilde{g}/\widetilde{g}_{0}$ at $t=30\ \textrm{fm}/c$ (lower
panel), for a series of initialization points, Eq. (10), marked with green
dots on the cumulant diagrams. The description of the cumulant diagrams is the
same as in Fig. 8. The deviation of $\widetilde{g}/\widetilde{g}_{0}$ from the
behavior of an uncorrelated system (red line) directly follows the deviation
of $\kappa_{2}/\kappa_{1}$ from the Poissonian limit of 1. Figure from [177].

We want to stress that the pair distributions shown in Fig. 11 develop
relatively fast. In Fig. 10, where we explored the behavior of a system
initialized at a temperature $T=1\ \textrm{MeV}$, one can see by comparing the
second and the fourth panels that already at $t=10\ \textrm{fm}/c$ a
significant part of the pair distribution has developed. This effect is
further magnified at higher temperatures, where relatively larger momenta of
the test particles result in a faster propagation of effects related to mean
fields, and consequently the system settles faster as well. For systems shown
in Fig. 11, we have verified that the majority of the pair distribution
function development occurs within $\Delta t=3\ \textrm{fm}/c$.

These results show not only that hadronic transport is sensitive to critical
behavior of systems evolving above the critical point, but also that this
behavior is exactly what is expected based on the underlying model. Moreover,
we note that the behavior of both the second-order cumulant and the pair
distribution function across the region of the phase diagram affected by the
critical point is remarkably distinct. It is evident that an equilibrated
system traversing the phase diagram through the series of chosen points, Eq.
(10), follows a clear pattern: first displaying anticorrelation, then
correlation, and then again anticorrelation. Thus already the second-order
cumulant ratio presents sufficient information to explore the phase diagram,
and, provided that correlations in the coordinate space are transformed into
correlations in the momentum space during the expansion of the fireball, this
pattern may be utilized to help locate the QCD critical point, in addition to
signals carried by the third- [97] and fourth-order [98] cumulant ratios. This
may prove to be especially important given that the quantity observed in
heavy-ion collision experiments is not the net baryon number, but the net
proton number. In calculations of the net baryon number cumulants based on the
net proton number cumulants, the higher order observables are increasingly
more affected by Poisson noise [102]. In view of this, the second-order
cumulant ratio (or equivalently the two-particle correlation) could be
considered among the key observables utilized in the search for the QCD
critical point, and it remains to be seen if this somewhat smaller signal (as
compared to higher order cumulant ratios) is nevertheless noteworthy due to
the much higher precision with which it can be measured in experiments.

####  3 Effects of finite number statistics

Qualitative and quantitative features of observables are influenced by the
finite number of particles in the analyzed samples. When analyzing observables
such as the baryon distribution, one has to keep in mind that fluctuations due
to finite number statistics may wash out the expected signals. This is not
only a numerical problem but, as we shall discuss below, it is also an issue
relevant for experiments.

Figure 12: Time evolution of the continuous baryon number distribution, scaled
by the volume of the cell and shown in units of the saturation density of
nuclear matter, $n_{0}=0.160\ \textrm{fm}^{-3}$, for the same system as
described in Fig. 9. Here, the cell width is $\Delta l=1\ \textrm{fm}$. The
red histogram shows the baryon distribution at initialization ($t=0$), while
the blue shaded histogram shows the distribution at the end of the evolution
($t_{\textrm{end}}=50\ \textrm{fm}/c)$. Nuclear matter, initialized in a
mechanically unstable region of the phase diagram, spontaneous separates into
a less dense and a more dense nuclear liquid, and the distribution function
becomes wider; however, due to the size of the binning cell, the average
number of test particles in a cell is small and consequently the double-peaked
structure, clearly seen on the right panel in Fig. 9, is washed out by
Poissonian fluctuations. Figure from [177].

Figure 13: Time evolution of the continuous baryon distribution, scaled by the
volume of the cell and shown in units of the saturation density of nuclear
matter, $n_{0}=0.160\ \textrm{fm}^{-3}$, for the same system as described in
Figs. 9 and 12, but calculated using the parallel ensembles method; the
results are averaged over $N_{\textrm{ev}}^{(\textrm{parallel})}=N_{T}\times
N_{\textrm{ev}}=25,000$ events. Cell width is $\Delta l=2\ \textrm{fm}$. The
red histogram shows the distribution at initialization ($t=0$), while the blue
shaded histogram shows the distribution at the end of the evolution
($t_{\textrm{end}}=50\ \textrm{fm}/c)$. Nuclear matter, initialized in a
mechanically unstable state, spontaneously separates into a less dense and a
more dense nuclear liquid, and the distribution widens; however, small numbers
of particles in cells used to construct the histogram and related finite
number effects wash out the structure clearly seen in the right panel of Fig.
9. Figure from [177].

First, we discuss this subject in the context of the choice of binning width.
In particular, the double-peak structure in the baryon number distribution
shown in the right panel of Fig. 9 depends on the size of the cell used to
construct the histogram, chosen to be $\Delta l=2\ \textrm{fm}$. In this case,
the Poissonian finite number statistics superimposed on the underlying baryon
distribution is characterized by a certain width $\sigma_{(\textrm{2~{}fm})}$.
If we reduce the cell width $\Delta l$ by a factor of 2, the average number of
particles in a cell is reduced by a factor of 8. Consequently, the width of
the Poissonian fluctuations will be
$\sigma_{(\textrm{1~{}fm})}=2\sqrt{2}\sigma_{(\textrm{2~{}fm})}$, which is
considerably larger than previously and which in fact washes out the double-
peak structure. This can be seen in Fig. 12, where we show the baryon number
distribution for a sampling cell width of $\Delta l=1\ \textrm{fm}$ for the
same events as used to create Fig. 9; the red and blue lines correspond to the
distribution at time $t=0$ and $t_{\textrm{end}}=50\ \textrm{fm}/c$,
respectively. For the system at hand, the Poissonian widths in the two cases,
in terms of baryon density, were $\sigma_{(\textrm{2~{}fm})}=0.22n_{0}$ and
$\sigma_{(\textrm{1~{}fm})}=0.62n_{0}$. If we then estimate the full width at
half maximum as approximately given by $2.355\sigma$ (the full width at half-
maximum of a normal distribution), it is clear that in the case of the cell
width $\Delta l=1\ \textrm{fm}$, the full width is comparable with the
separation of the peaks given by the width of the coexistence region,
$n_{R}-n_{L}=1.44n_{0}$. As a result, the two-peak structure cannot be
resolved for this sampling statistics. Let us note here that decreasing the
volume of the cells $(\Delta l)^{3}$ can be done without penalty if one
proportionally increases the number of test particles per particle, $N_{T}$.
Conversely, decreasing the number of test particles per particle $N_{T}$
exacerbates the effects of finite number statistics.

Figure 14: Pair distribution function at $t=30\ \textrm{fm}/c$ for a system
initialized inside the quark-hadron spinodal region (at baryon number density
$n_{B}=3.0n_{0}$ and temperature $T=1\ \textrm{MeV}$), calculated with the
parallel ensembles method; the results are averaged over
$N_{\textrm{ev}}^{(\textrm{parallel})}=N_{T}\times N_{\textrm{ev}}=25,000$
events. Spontaneous spinodal decomposition leads to a formation of two
coexisting phases: a “less dense” and a “more dense” nuclear liquid. The
increased relative concentration of particles in the “more dense” phase
results in an elevated normalized pair correlation at small distances. The
correlation is exactly the same as shown on the rightmost panel in Fig. 10.
Figure from [177].

While this discussion may appear to be of purely numerical nature,
experimental data are similarly affected by finite number statistics. In
experiments, one always deals with exactly $N_{B}$ particles per event, which
in our simulations corresponds to $N_{T}=1$. Naturally, it must lead to a
distribution in which any possible peaks are even more washed out. This can be
seen in Fig. 13, where we show results for the case of $N_{T}=1$ and $\Delta
l=2\ \textrm{fm}$; the red and blue lines correspond to the distribution at
time $t=0$ and $t_{\textrm{end}}=50\ \textrm{fm}/c$, respectively. Here, in
order to ensure that we are comparing systems with identical dynamics, we used
the same simulation data as in Figs. 9 and 12, but this time we accessed the
baryon number distribution corresponding to $N_{T}=1$ using the parallel
ensembles method (for details, see Section 3 and Appendix 25). Not
surprisingly, the signal is almost entirely washed out and only a slight
broadening of the distribution is discernible. We note that increasing the
number of events does not resolve this issue, as the resolution is determined
by Poissonian fluctuations in individual events. Consequently, one needs to
devise other methods to extract the information about the underlying baryon
distribution; this task is beyond the scope of this work.

Finally, we note that the pair distribution function is less affected by
finite number statistics. In Fig. 14, we show the pair distribution function
calculated within the parallel ensembles method, which is nearly identical to
the pair distribution function calculated in the full ensemble, Fig. 10.
Indeed, the normalized pair distribution function is not determined by the
total number of test particles in an event or in a given subvolume of the
system, but by relations between any two test particles. The only difference
between the pair distribution functions obtained within the two methods is in
the error bars, which are larger in the parallel ensembles case due to smaller
statistics: the number of pairs in the full ensemble is given by
$N_{\textrm{ev}}(N_{B}N_{T})^{2}$, while in the parallel ensembles it’s equal
$N_{T}N_{\textrm{ev}}(N_{B})^{2}$. Obtaining the same pair distribution
function demonstrates that the physics accessible in the full ensemble and the
parallel ensembles approach is the same.

### 3\. Summary

Results from simulations in SMASH, presented above, demonstrate that critical
behavior in dense nuclear matter can be studied within a hadronic transport
approach equipped with interactions corresponding to a chosen EOS. In
particular, we have shown that systems initialized in unstable regions of the
phase diagram undergo spontaneous spinodal decomposition, followed by an
evolution towards an equilibrated mixture of two coexisting phases with
compositions matching the predictions from the underlying EOS. The correct
description of both the thermodynamics and the non-equilibrium phenomena
implies that hadronic transport can be used as a tool with unique capabilities
to investigate the dynamic evolution of matter created in heavy-ion
collisions.

We have also shown that for systems initialized at various points of the phase
diagram, including in the vicinity of a critical point, the pair distribution
functions calculated from hadronic transport simulation data follow
theoretical expectations based on the second-order cumulant ratio
${\kappa_{2}}/{\kappa_{1}}$. In particular, as the baryon number density (and,
consequently, baryon chemical potential) is increased in the region of the
phase diagram affected by the critical point, the pair distribution function
follows a clear pattern: displaying first anticorrelation, then correlation,
and then again anticorrelation. The behavior of two-particle correlations
(and, on the theoretical side, of the second-order cumulant ratio
${\kappa_{2}}/{\kappa_{1}}$) is therefore a clear signature of crossing the
phase diagram above the critical point. This is especially important in view
of the experimental search for the QCD critical point, as lower order
statistical observables, such as ${\kappa_{2}}/{\kappa_{1}}$, are more likely
to be measured with accuracy sufficient for discerning signals of critical
behavior.

Further, we have shown that finite number statistics affects both the
qualitative and quantitative features of statistical observables. We have
discussed two complementary simulation paradigms (described in Section 1 and
in Appendix 25), within which hadronic transport gives access to both the
continuous baryon number distribution, employed in theoretical calculations,
and the physical baryon number distribution relevant to experimental results.
Though driven by the same physics, these distributions lead to starkly
different values for integrated statistical observables, which has
consequences for comparisons with experimental data.

##  Chapter 5 Beyond hadronic transport: the speed of sound in heavy-ion
collisions

The contents of this chapter are largely based on Ref. [212].

One of the ultimate goals of heavy-ion collision studies is connecting
experimental observables to the underlying EOS of dense nuclear matter. Here,
we present a method that may allow an estimate of the value of the speed of
sound as well as its logarithmic derivative with respect to the baryon number
density in matter created in heavy-ion collisions. To this end, we utilize
well-known observables: cumulants of the baryon number distribution. In
analyses aimed at uncovering the phase diagram of strongly interacting matter,
cumulants gather considerable attention as their qualitative behavior along
the explored range of collision energies is expected to aid in detecting the
QCD critical point (see Sections 6 and 4 for more details). We show that the
cumulants may also reveal the behavior of the speed of sound in the
temperature and baryon chemical potential plane. We demonstrate the
applicability of such estimates within two models of nuclear matter: the
Walecka model and the VDF model presented in this thesis (see Section 3 and
Chapter 2), and we explore what might be understood from known experimental
data.

### 1\. Background

The speed of sound $c_{s}$ is one of the fundamental properties of any
substance. In fluids, it is the velocity with which a longitudinal compression
wave propagating through the medium, and its square can be computed as the
ratio of a change in the pressure $P$ corresponding to a change in the energy
density $\mathcal{E}$. As such, it is directly related to a number of
thermodynamic properties of the system, including its equation of state (EOS).

In dense nuclear matter, $c_{s}$ is of particular interest to neutron star
research: its behavior as a function of baryon number density $n_{B}$ has a
direct influence on the mass-radius relationship and, consequently, on the
maximum possible mass of a neutron star [13]. Current neutron star data
suggest that $c_{s}$ rises significantly for an $n_{B}$ larger than the
nuclear saturation density $n_{0}$, and that it approaches and perhaps exceeds
$c_{s}\sim 1/\sqrt{3}$ at densities as low as a few times that of normal
nuclear matter. This possibility was first suggested in Ref. [190] and
followed by multiple studies on the subject, see, e.g., Refs. [213, 184, 15].

Presently, heavy-ion collisions are the only means of studying dense nuclear
matter in a laboratory. Many experiments probing nuclear matter at high
$n_{B}$, such as the Beam Energy Scan (BES) program at the Relativistic Heavy
Ion Collider (RHIC), put special significance on the search for the QCD
critical point (CP). Also in this case the behavior of $c_{s}$ conveys
relevant information: a crossover transition is characterized by a local
minimum of the speed of sound, whereas at the CP and on the associated
spinodal lines of the first-order phase transition, the speed of sound
vanishes. Indeed, lattice QCD (LQCD) results show that at vanishing baryon
chemical potential, $\mu_{B}=0$, a local minimum in $c_{s}$ occurs at
temperature $T_{0}=156.5\pm 1.5$ MeV [25] (see also Ref. [27]), corresponding
to a crossover transition between hadron gas and quark-gluon plasma (QGP).

To date, a few attempts have been made to evaluate $c_{s}$ from heavy-ion
collision data. In Ref. [214], $c_{s}$ is estimated in ultrarelativistic
collisions based on the proportionality of entropy density $s$ and temperature
$T$ to charged particle multiplicity and mean transverse momentum,
respectively, and the obtained value agrees with LQCD results. Unfortunately,
the method used relies on the fact that at highest collision energies
$\mu_{B}\approx n_{B}\approx 0$, making this approach inapplicable for most of
the energy range covered at RHIC. At values of the baryon chemical potential
relevant to the BES program, $c_{s}$ was estimated in Ref. [215], where both
the Landau model as well as hybrid hydrodynamics and hadronic transport
simulations were used to reproduce the widths of the negatively charged pion
rapidity distribution. That study purports to locate a robust minimum in
$c_{s}$ within the collision energy range $\sqrt{s_{NN}}=4\textrm{-}9\
\textrm{GeV}$.

Below, we suggest a novel approach to exploring the behavior of $c_{s}$ by
using cumulants of the baryon number distribution. The sensitivity of the
cumulants to the EOS near the CP [97, 98], which makes them central
observables pursued in the Beam Energy Scan, follows directly from their
sensitivity to derivatives of the pressure with respect to $\mu_{B}$. The key
observation we make here is that, besides the vicinity of the CP, cumulants of
the baryon number distribution provide rich information about the EOS at all
points of the phase diagram, and in particular they allow a measurement of
$c_{s}$ in matter created in heavy-ion collisions.

### 2\. Relation between cumulants and the speed of sound

Cumulants of net baryon number $\kappa_{j}$ can be obtained from
$\kappa_{j}=VT^{j-1}\left(d^{j}P/d\mu_{B}^{j}\right)_{T}$, where $V$ is the
volume. Expressed in terms of derivatives with respect to $n_{B}$, the first
three cumulants are given by Eqs. (28-30), repeated here for convenience,

$\displaystyle\kappa_{1}=Vn_{B}~{},$ (1)
$\displaystyle\kappa_{2}=\frac{VTn_{B}}{\left(\frac{dP}{dn_{B}}\right)_{T}}~{},$
(2)
$\displaystyle\kappa_{3}=\frac{VT^{2}n_{B}}{\left(\frac{dP}{dn_{B}}\right)_{T}^{2}}\left[1-\frac{n_{B}}{\left(\frac{dP}{dn_{B}}\right)_{T}}\left(\frac{d^{2}P}{dn_{B}^{2}}\right)_{T}\right]~{}.$
(3)

Importantly, cumulants can be also given in terms of moments of the net baryon
distribution, which can be directly measured in experiment; in particular, for
$j\leq 3$,
$\kappa_{j}\equiv\big{\langle}\big{(}N_{B}-\big{\langle}N_{B}\big{\rangle}\big{)}^{j}\big{\rangle}$
(for more details, see Sections 6 and 4).

The exact definition of $c_{s}$ depends on specifying which properties of the
system can be considered constant during the propagation of the compression
wave. One often uses the speed of sound at constant entropy $S$ per net baryon
number $N_{B}$, $c_{\sigma}^{2}\equiv\left(dP/d\mathcal{E}\right)_{\sigma}$,
where $\sigma=S/N_{B}$. Similarly, the speed of sound at constant temperature
is $c_{T}^{2}\equiv\left(dP/d\mathcal{E}\right)_{T}$. These two variants have
specific regions of applicability. For example, the propagation of sound in
air is governed by adiabatic compression, so that using $c_{\sigma}^{2}$ is
appropriate. On the other hand, when there is a temperature reservoir (e.g.,
in porous media) or when the cooling timescale is very fast compared with the
sound wave period (as is the case, e.g., for an interstellar medium subject to
radiative cooling), $c_{T}^{2}$ is applicable.

Explicitly, $c_{\sigma}^{2}$ and $c_{T}^{2}$ can be written as (see Appendix
20 for a detailed derivation)

$\displaystyle
c_{\sigma}^{2}=\frac{\Big{(}\frac{dP}{dn_{B}}\Big{)}_{T}\Big{(}\frac{ds}{dT}\Big{)}_{n_{B}}+\Big{(}\frac{dP}{dT}\Big{)}_{n_{B}}\bigg{[}\frac{s}{n_{B}}-\Big{(}\frac{ds}{dn_{B}}\Big{)}_{T}\bigg{]}}{\Big{(}\frac{sT}{n_{B}}+\mu_{B}\Big{)}\Big{(}\frac{ds}{dT}\Big{)}_{n_{B}}}$
(4)

and

$\displaystyle
c_{T}^{2}=\frac{\Big{(}\frac{dP}{dn_{B}}\Big{)}_{T}}{T\Big{(}\frac{ds}{dn_{B}}\Big{)}_{T}+\mu_{B}}~{},$
(5)

respectively. In the limit $T\to 0$, the above expressions both lead to

$\displaystyle
c^{2}\Big{|}_{T=0}=\frac{1}{\mu_{B}}\bigg{(}\frac{dP}{dn_{B}}\bigg{)}_{T}~{}.$
(6)

Consequently, for $(\mu_{B}/T)\gg 1$, the values of $c_{\sigma}^{2}$ and
$c_{T}^{2}$ should largely coincide. Moreover, Eq. (5) can be transformed to
express $c_{T}^{2}$ as a function of the cumulants, Eqs. (1-3) ,

$\displaystyle
c_{T}^{2}=\left[\bigg{(}\frac{\partial{\log\kappa_{1}}}{\partial{\log
T}}\bigg{)}_{\mu_{B}}+\frac{\mu_{B}}{T}\frac{\kappa_{2}}{\kappa_{1}}\right]^{-1}~{}.$
(7)

Therefore, cumulants can be used to compute $c_{T}^{2}$ if they are known as
functions of $T$ and $\mu_{B}$. The first term in Eq. (7) is challenging to
estimate from experimental data, however, it can be shown to be negligible for
a degenerate Fermi gas, $(\mu_{B}/T)\gg 1$, where it constitutes an order
$\left(T/\mu_{B}\right)^{2}$ correction; in that case we can write

$\displaystyle c_{T}^{2}\approx\frac{T\kappa_{1}}{\mu_{B}\kappa_{2}}~{}.$ (8)

We note that Eq. (8) provides an upper limit to the value of $c_{T}^{2}$ as
long as $(\partial\log\kappa_{1}/\partial\log T)_{\mu_{B}}>0$.

Using Eq. (5), one can also calculate the logarithmic derivative of
$c_{T}^{2}$,

$\displaystyle\bigg{(}\frac{d\ln c_{T}^{2}}{d\ln
n_{B}}\bigg{)}_{T}=\frac{n_{B}\Big{(}\frac{d^{2}P}{dn_{B}^{2}}\Big{)}_{T}}{\Big{(}\frac{dP}{dn_{B}}\Big{)}_{T}}-\frac{\Big{(}\frac{dP}{dn_{B}}\Big{)}_{T}+Tn_{B}\Big{(}\frac{d^{2}s}{dn_{B}^{2}}\Big{)}_{T}}{\mu_{B}+T\Big{(}\frac{ds}{dn_{B}}\big{)}_{T}}~{}.$
(9)

It is again possible to represent the above equation in terms of the
cumulants,

$\displaystyle\bigg{(}\frac{d\ln c_{T}^{2}}{d\ln
n_{B}}\bigg{)}_{T}+c_{T}^{2}=1-\frac{\kappa_{3}\kappa_{1}}{\kappa_{2}^{2}}-c_{T}^{2}\bigg{(}\frac{d\ln(\kappa_{2}/T)}{d\ln
T}\bigg{)}_{n_{B}}~{},$ (10)

and neglecting the last term on the right-hand side yields

$\displaystyle\left(\frac{d\ln c_{T}^{2}}{d\ln
n_{B}}\right)_{T}+c_{T}^{2}\approx
1-\frac{\kappa_{3}\kappa_{1}}{\kappa_{2}^{2}}~{}.$ (11)

This approximation is again valid for $(\mu_{B}/T)\gg 1$, and the correction
due to the neglected term is likewise of order $\left(T/\mu_{B}\right)^{2}$.

It is worth noting that in the opposite limit, $\mu_{B}\to 0$, Eq. (7) reveals
a similarly simple form, $c_{T}^{2}=\left(d\ln\kappa_{2}/d\ln
T\right)^{-1}_{\mu_{B}=0}$ (see Appendix 26 for the derivation), suggesting
that $c_{T}^{2}$ can be estimated in ultrarelativistic heavy-ion collisions,
provided measurements of $\kappa_{2}$ are available at different temperatures
with sufficient precision. It might be possible to achieve this with data from
a combination of centralities, energies, collision species, or rapidity
ranges. At this time, however, we are interested in utilizing Eqs. (8) and
(11) applied to collisions at medium and low energies.

Figure 1: Model study of regions of applicability of Eqs. (8) and (11). The
left (right) panels show results obtained in the VDF (Walecka) model. The
upper and lower panels show quantities entering Eq. (8) and Eq. (11),
respectively. Results at $T=50,100,150,200\ \textrm{MeV}$ are given by blue
and green solid lines, dark and light purple long-dashed lines, red and pink
short-dashed lines, and orange and brown dash-dotted lines. For each $T$, the
thickest lines show the exact results and the medium-thick lines show the
approximations, given by the right-hand sides of Eqs. (8) and (11).
Additionally, in upper panels the thinnest lines show $c_{\sigma}^{2}$. For
both models, Eq. (8) is valid for $T\lesssim 100\ \textrm{MeV}$ and
$\mu_{B}\gtrsim 600\ \textrm{MeV}$, while Eq. (11) is valid for all values of
$T$ and $\mu_{B}\gtrsim 200\ \textrm{MeV}$, with the exception of the Walecka
model at $T=200\ \textrm{MeV}$, where a phase transition to an almost massless
gas of nucleons decreases the applicability of both Eqs. (8) and (11). Figure
from [212].

### 3\. Validation

We are interested in finding the limitations of the low-temperature
approximation used to derive Eqs. (8) and (11), and for this we turn to
effective models. Because the regions of the phase diagram where we expect to
apply our formulas are described by hadronic degrees of freedom, we choose two
models of dense nuclear matter: the VDF model with two phase transitions (see
Chapter 2) and the Walecka model [154, 155]. The VDF model utilizes only
interactions of the vector type, while the Walecka model employs both vector-
and scalar-type interactions. Both models describe the nuclear liquid-gas
phase transition, while the VDF model additionally describes a conjectured
high-density, high-temperature phase transition modeling the QGP phase
transition. In this study, the QGP-like phase transition is chosen to exhibit
a CP at $T_{c}=100\ \textrm{MeV}$ and $n_{c}=3n_{0}$, with the $T=0$
boundaries of the spinodal region in the $(T,n_{B})$ plane given by
$n_{B,\textrm{left spinodal}}(T=0)\equiv\eta_{L}=2.5n_{0}$ and
$n_{B,\textrm{right spinodal}}(T=0)\equiv\eta_{R}=3.32n_{0}$, where
$n_{0}=0.160\ \textrm{fm}^{-3}$ (which altogether correspond to the fourth
(IV) set of characteristics listed in Table 1); this choice is arbitrary and
serves as a plausible example.

To test the derived expressions, we plot both sides of Eqs. (8) and (11) as
functions of $\mu_{B}$ at a series of temperatures in Fig. 1. We note that the
explored temperature range reaches beyond the region of the phase diagram
where hadronic models plausibly describe matter created in heavy-ion
collisions; nevertheless, it is still instructive to test the validity of our
approximations at these extreme conditions.

We use natural units in which the speed of light in vacuum is $c=1$. We note
that in the VDF model, $c_{s}$ quickly becomes acausal for $\mu_{B}$ exceeding
the values defining the coexistence line of the QGP-like phase transition. As
we already detailed in Section 2, it is an expected behavior in models using
interactions dependent on high powers of $n_{B}$, and while generally not
ideal, it does not affect the current analysis.

In all panels in Fig. 1, the exact model calculations show expected features
as functions of $\mu_{B}$. In the upper left panel, showing both $c_{T}^{2}$
and $c_{\sigma}^{2}$ in the VDF model, at small $\mu_{B}$ we see a softening
of the EOS due to the influence of the nuclear CP, followed by an increase at
densities of the order of $n_{0}$, then a deep dive in $c_{s}^{2}$ caused by
the QGP-like phase transition, and finally a steep rise for high values of
$\mu_{B}$. In the upper right panel, showing $c_{T}^{2}$ and $c_{\sigma}^{2}$
in the Walecka model, we similarly observe a soft EOS at small $\mu_{B}$,
while the value of $c_{s}^{2}$ goes asymptotically to 1 for large $\mu_{B}$.
Additionally, for $T=200\ \textrm{MeV}$, the Walecka model shows effects due
to a phase transition in the nucleon-antinucleon plasma, occurring around
$T\approx 190\ \textrm{MeV}$ and $n_{B}=0$; above this transition, the model
describes an almost noninteracting gas of nearly massless nucleons [216]. The
behavior of the curves in the lower panels, showing $c_{T}^{2}+\big{(}{d\ln
c_{T}^{2}}/{d\ln n_{B}}\big{)}_{T}$, can be directly traced to the behavior of
the curves in the upper panels. In particular, for the VDF model we observe
strong divergences due to the softening of the EOS in the QGP-like phase
transition region.

Comparing the exact results to the approximations, we see that while Eq. (8)
is valid for small to moderate temperatures, $T\lesssim 100\ \textrm{MeV}$,
and moderate to high baryon chemical potentials, $\mu_{B}\gtrsim 600\
\textrm{MeV}$, it behaves poorly, both qualitatively and quantitatively, for
$T$ and $\mu_{B}$ corresponding to regions of the phase diagram probed by
moderately to highly energetic heavy-ion collisions (upper panels). On the
other hand, the approximation introduced in Eq. (11) is qualitatively valid
for most of the probed $T$ and $\mu_{B}$, with the exception of regions
characterized by low chemical potentials, $\mu_{B}\lesssim 200\ \textrm{MeV}$
(lower panels).

### 4\. Experimental data and interpretation

We proceed to apply Eqs. (8) and (11) to heavy-ion collision data. We consider
cumulants of the net proton number and chemical freeze-out parameters,
$(T_{\textrm{fo}},\mu_{\textrm{fo}})$, in collisions at 0–5% centrality,
determined by the solenoidal tracker at RHIC (STAR) [104] and high acceptance
dielectron spectrometer (HADES) [217, 218] experiments, and we use them to
plot Eqs. (8) and (11) (red triangles, left and right panel in Fig. 2,
respectively) against $\mu_{B}$.

Figure 2: Comparison of the right-hand sides of Eq. (8) (left panel) and Eq.
(11) (right panel) for experimental data (red triangles), ideal gas at the
freeze-out (small gray circles), the VDF model at the freeze-out (light green
stars), and the VDF model at a set of points chosen to reproduce the data
(dark purple stars); exact results, that is, the left-hand sides of Eqs. (8)
and (11), are shown for the two cases considered in the VDF model (green and
purple circles). The data points for the matched VDF results are chosen to
reproduce experimental values of $1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$ (see
Fig. 3). We note that at $\sqrt{s}=2.4\ \textrm{GeV}$, matching the value of
$1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$ exactly is possible, but would place
the matched point close to the nuclear liquid-gas CP, which we find unlikely.
Figure from [212].

Based on the model validation study presented in Section 3, we trust the
results presented in the left panel of Fig. 2, approximating $c_{T}^{2}$, only
for the lowest collision energy, $\sqrt{s}=2.4\ \textrm{GeV}$ from the HADES
experiment. At this collision energy, the value of $c_{T}^{2}$ as obtained
from Eq. (8) is small: less than half of the ideal gas value (see insert). At
the same time, the value of $1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$ (shown in
the right panel in Fig. 2), which we assume is dominated by $\left({d\ln
c_{T}^{2}}/{d\ln n_{B}}\right)_{T}$, drops with decreasing collision energy to
reach a minimum at the lowest STAR point, $\sqrt{s}=7.7\ \textrm{GeV}$, and
then steeply rises for the HADES point. This could mean that in matter created
in $\sqrt{s}=7.7\ \textrm{GeV}$ collisions, $c_{T}^{2}$ is approximately
constant as a function of $n_{B}$, while in $\sqrt{s}=2.4\ \textrm{GeV}$
collisions the matter is characterized by a small $c^{2}_{T}$ which
nevertheless has a large positive slope as a function of $n_{B}$.

Figure 3: Contour plot of $1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$ in the VDF
model. Yellow (black) lines correspond to the spinodal (coexistence) lines,
while white contours mark regions where
$1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}=1\pm 0.03$. Light green stars denote
experimentally measured freeze-out parameters
$(T_{\textrm{fo}},\mu_{\textrm{fo}})$, while dark purple stars denote points
where $1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$, taken along lines informed by
average phase diagram trajectories for STAR collision energies [219], matches
the experimentally measured values for a given collision energy. Figure from
[212].

To further understand this behavior, we study the dependence of
$1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$ on $\mu_{B}$ and $T$ within the VDF
model, shown as a contour plot in Fig. 3. Remarkably, the softening of the
EOS, leading to negative values of $\left({d\ln c_{T}^{2}}/{d\ln
n_{B}}\right)_{T}$, occurs in two regions of the phase diagram, corresponding
to the ordinary nuclear matter phase transition and to the conjectured QGP-
like phase transition. In order to compare model results to experiment, we
need to choose at which $T$ and $\mu_{B}$ to take values of
$1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$. A natural choice is to use values at
the freeze-out points $(T_{\textrm{fo}},\mu_{\textrm{fo}})$ (marked in Fig. 3
with light green stars), however, as can be seen in the right panel of Fig. 2,
the VDF model values of $1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$ taken at these
points (also denoted with light green stars) do not lead to an agreement with
experimental data; in particular, the biggest discrepancy occurs for the HADES
point. (We note here that the values of the freeze-out parameters
$(T_{\textrm{fo}},\mu_{\textrm{fo}})$ are established with hadron interactions
neglected, and the degree to which this affects our results may vary
significantly across the phase diagram.) However, critical fluctuations are
shown to exhibit a large relaxation time [105, 139, 220], and as such their
measured values could be affected by stages of the collision preceding the
chemical freeze-out. Guided by this insight and aided by the average phase
diagram trajectories of hybrid simulations of heavy-ion collisions [219], we
consider values of $1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$ at points in the
phase diagram corresponding to slightly earlier stages of the collisions.
Additionally demanding that the values of
$1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$ reproduce the experimentally obtained
values for a given collision energy results in points denoted with dark purple
stars in both Fig. 2 and in Fig. 3 (where we choose to show points only for
collisions at low energies, for which using the model is justified); here the
exception from the exact matching occurs for the HADES point, for which we
prioritized choosing a point in a reasonable vicinity of the measured freeze-
out over obtaining a value equal to the experimental data. Comparing to the
exact model results for $c_{T}^{2}$ and $c_{T}^{2}+\big{(}{d\ln
c_{T}^{2}}/{d\ln n_{B}}\big{)}_{T}$, also displayed in Fig. 2, as well as to
the upper left panel of Fig. 1, we can confirm that at the point reproducing
the experimental value of $1-\kappa_{3}\kappa_{1}/\kappa_{2}^{2}$ for the
lowest STAR energy, $c_{T}^{2}$ is nearly constant as a function of $\mu_{B}$
(thick short-dashed line, $T=150\ \textrm{MeV}$, at $\mu_{B}\approx 410\
\textrm{MeV}$ in Fig. 1), while at the point reproducing the result for the
HADES energy, $c_{T}^{2}$ increases sharply with $\mu_{B}$ (thick solid line,
$T=50\ \textrm{MeV}$, at $\mu_{B}\approx 850\ \textrm{MeV}$ in Fig. 1).

Naturally, the choice of $T$ and $\mu_{B}$ at which we compare model
calculations with STAR and HADES cumulant data is driven by the wish to match
the experimental results, and as such it serves mainly to show that baryon
number cumulants measured in heavy-ion collisions can be connected to the
speed of sound in hot and dense nuclear matter. Whether values of higher order
cumulants are indeed significantly affected by stages of the evolution
preceding the freeze-out needs to be further investigated (for recent
developments, see Ref. [106, 221]). Moreover, while experiments measure proton
number cumulants, the VDF model provides baryon number cumulants, putting more
strain on a direct interpretation of our results. Effects due to baryon number
conservation should likewise be important [222, 103]. Finally, it is
understood that our model results may not be applicable in regions of the
phase diagram where quark and gluon degrees of freedom become increasingly
relevant.

Nonetheless, hadronic models are well-justified for describing low-energy
collisions whose evolution is dominated by the hadronic stage, and it is the
results in this region of the phase diagram that we want to stress. The
comparison between the experimental data and the VDF model suggests that
collisions at the lowest STAR and HADES energies may be probing regions of the
phase diagram where the cumulants of the baryon number tell us more about
hadronic physics than the QCD CP. In particular, the change in the sign of
$\kappa_{3}$, predicted to take place in the vicinity of a critical point [97]
and apparent in the HADES data (see the right panel of Fig. 2), may mark the
region of the phase diagram affected by the nuclear liquid-gas phase
transition. If this is true, it may be worthwhile to study the cumulants at
even lower collision energies, starting from $0.1\ \textrm{GeV}$ projectile
kinetic energy, and obtain the speed of sound around the nuclear liquid-gas
CP. Conversely, at higher energies it could be possible to use collisions at
different centralities and different rapidity windows to estimate the
neglected terms in Eqs. (8) and (11), and obtain a stronger estimate for the
speed of sound in the respective regions of the phase diagram.

### 5\. Summary

In this chapter, we have used cumulants of the baryon number distribution to
estimate the isothermal speed of sound squared and its logarithmic derivative
with respect to the baryon number density. The approximations and the model
comparisons we considered here apply to experiments at low energies, however,
the presented method can be used at any collision energy provided that
measurements of cumulants of the baryon number distribution as well as their
temperature dependence are available. Further studies of effects due to
dynamics, in particular using state-of-the-art simulations, will be essential
in determining the extent to which the proposed method provides a reliable
extraction of sound velocities and their derivatives. Nevertheless, this
result provides a new approach to obtaining information about fundamental
properties of nuclear matter studied in heavy-ion collisions, with
consequences for both the search for the QCD critical point and neutron star
studies.

##  Chapter 6 Thesis summary

The work described in this thesis is devoted to developing, testing, and
applying a parametrizable equation of state (EOS) of dense nuclear matter to
studies centered on the phase diagram of quantum chromodynamics (QCD). The
presented formalism has the potential to help interpret the measurements
obtained from heavy-ion collisions experiments and, in particular, to reveal
the thermodynamics of strongly-interacting QCD matter. Because the constructed
family of EOSs shows a substantial flexibility in postulating the properties
of the possible QCD phase transition at high temperatures and high baryon
densities, it is perfectly suited for use in systematic studies of effects of
different dense nuclear matter EOS on final state observables (e.g., using
Bayesian analysis) and can facilitate meaningful comparisons of simulation
results with experimental data.

Below, we provide a summary of the developments presented in this work,
followed by a brief outlook on future research directions.

In Chapter Density Functional Equation of State  and Its Application to the
Phenomenology of Heavy-Ion Collisions, we have given an overview of the search
for the conjectured first-order phase transition in dense QCD matter and the
most relevant results obtained within the Phase I of the Beam Energy Scan
program at the Relativistic Heavy-Ion Collider. We have also discussed the
state-of-the-art approaches to simulations of heavy-ion collisions and the
indispensable role they will play in the interpretation of experimental data
from the Phase II of the Beam Energy Scan.

In Chapter 1, after an introductory overview of the mean-field and Landau
Fermi-liquid theory approaches to describing strongly-interacting many-body
systems, we have developed an easily parametrizable EOS with an arbitrary
number of scalar- and vector-density-functional–based interactions (the VSDF
model). We have also shown that this model leads to Lorentz-covariant
equations of motion, preserves conservation laws, and is thermodynamically
consistent.

In Chapter 2, we first discussed the parametrization procedure used to obtain
a family of EOSs based on the version of the model introduced in Chapter 1
that utilizes only interactions of the vector type (the VDF model). With
application to studies of the phase diagram of QCD in mind, we chose
parametrizations of the VDF model with two phase transitions: one
corresponding to the well-known low-temperature, low-density phase transition
in ordinary nuclear matter, and one corresponding to a conjectured high-
temperature, high-density QGP-like phase transition, meant to model the
transition from matter described by hadronic degrees of freedom to matter
described by quark and gluon degrees of freedom. Then, based on several EOSs
chosen to represent the flexibility of the model, we discussed at length the
thermodynamic properties of described systems.

In Chapter 3, we have discussed the implementation of the equations of motion
stemming from the VDF model in the hadronic transport code SMASH, which we
started with an overview of modeling many-body mean-field dynamics within
hadronic transport approaches. In particular, we provided a detailed
description of several different density and mean-field calculation algorithms
available within SMASH, a few of which were introduced as part of the work
underlying this thesis, and we included a comparison of the numerical
performance of these options.

In Chapter 4, we have presented results of SMASH simulations of infinite
nuclear matter, realized as matter in a cubic box with periodic boundary
conditions. While the development of the V(S)DF model and the introduction of
its supporting framework in SMASH have been done to enable large-scale
comparisons between experimental data and simulations spanning a broad family
of EOSs, the investigation of the qualitative behavior of dense nuclear matter
presented in this thesis has been facilitated by using one specific
parametrization of the dense nuclear matter EOS. After detailing the analysis
methods used to study the simulation data, we discussed results pertaining to
two distinct regions of the QCD phase diagram. First, we discussed simulations
of nuclear matter in the density range corresponding to the ordinary nuclear
matter, and we confirmed that the implemented VDF mean-field potentials lead
to the correct thermodynamic behavior of nuclear matter as simulated in SMASH.
Then we discussed results of simulations run in the high baryon density
regime, affected by a conjectured QCD phase transition, stressing the
excellent agreement of the behavior of the simulated systems with the
thermodynamic predictions stemming from the underlying EOS. Furthermore, we
have shown that the calculated fluctuation observables reproduce the
thermodynamic behavior expected in the vicinity of the critical point based on
the model predictions, proving that hadronic transport can be sensitive to
effects due to phase transitions. In addition, we have discussed effects due
to finite number statistics, whose thorough understanding is crucial to
connecting theoretical calculations as well as simulation results to values of
fluctuations measured in experiment.

In Chapter 5, we have discussed a novel method of estimating the speed of
sound as well as its logarithmic derivative with respect to the baryon number
density in matter created in heavy-ion collisions. The method, based on a
connection between the isothermal speed of sound and the susceptibilities of
pressure characterizing the system, utilizes one of the most prominent
observables used in the search for the QCD critical point: the cumulants of
the baryon number distribution. The results of our analysis, based on
available data on the cumulants measured by the STAR and HADES experiments,