to a non-zero expectation values of the meson fields.

To calculate properties of the system described by the Walecka Lagrangian, Eq.
(1), the wave-function $\psi(x)$ must be found. Since the fields are strongly-
interacting, the couplings $g_{\omega}$ and $g_{\sigma}$ are large, which
precludes obtaining solutions using trivial perturbation theory. Several
simplifying assumptions can be made to proceed. First, one can take the local
density approximation, in which one neglects the derivative terms present in
Eqs. (5) and (6), leading to

$\displaystyle m_{\sigma}^{2}\sigma=g_{\sigma}\bar{\psi}\psi~{},$ (7)
$\displaystyle
m_{\omega}^{2}\omega^{\mu}=g_{\omega}\bar{\psi}\gamma^{\mu}\psi~{}.$ (8)

Formally, the local density approximation is equivalent to the limit of
$m_{\sigma}\to\infty,~{}m_{\omega}\to\infty$, and physically it implies that
the meson masses are significantly larger than any gradients present in the
system; importantly, while this implies that the system is somewhat smooth,
the density is still allowed to vary. Next, the non-zero expectation values of
the meson fields, already mentioned above, suggest that the fields can be
rewritten as sums of expectation values and fluctuations around those values,

$\displaystyle\sigma=\bar{\sigma}+\sigma^{\prime}~{},$ (9)
$\displaystyle\omega_{\mu}=\bar{\omega}_{\mu}+\omega^{\prime}_{\mu}~{},$ (10)

where $\bar{\sigma}$ and $\bar{\omega}_{\mu}$ are field averages. Further, if
the rest frame baryon density $n_{B}$ is large, then the source terms on the
right-hand side of Eqs. (7) and (8) are also large. In this case, it is a
valid approximation to consider the meson fields as well-described by
classical fields and replace their source terms by expectation values of the
scalar baryon density and the baryon 4-current,
$\langle\bar{\psi}\psi\rangle\equiv n_{s}$ and
$\langle\bar{\psi}\gamma^{\mu}\psi\rangle\equiv j^{\mu}$ (note that this is
exactly the way one treats classical electromagnetic fields). This approach,
known as the mean-field approximation, is equivalent to neglecting the meson
field fluctuations in Eqs. (9) and (10), and leads to replacing Eqs. (7) and
(8) by their expectation values,

$\displaystyle m_{\sigma}^{2}\bar{\sigma}=g_{\sigma}n_{s}~{},$ (11)
$\displaystyle m_{\omega}^{2}\bar{\omega}^{\mu}=g_{\omega}j^{\mu}~{}.$ (12)

Physically, the mean-field approximation corresponds to a situation in which
the propagation of the nucleon field, rather than being affected by the
particular values of the meson fields at a given point in space, can be
treated as a practically independent propagation in an average constant
background of the meson fields. Computationally, this means that the meson-
field operators are reduced to complex-valued fields.

Eqs. (11) and (12) can be used to replace the meson fields in the nucleon
field equation, Eq. (4), by the baryon vector and scalar currents,

$\displaystyle\Bigg{[}\gamma^{\mu}\bigg{(}i\partial_{\mu}-\frac{g_{\omega}^{2}}{m_{\omega}^{2}}j_{\mu}(x)\bigg{)}-\bigg{(}m_{N}-\frac{g_{\sigma}^{2}}{m_{\sigma}^{2}}n_{s}(x)\bigg{)}\Bigg{]}\psi(x)=0~{};$
(13)

here, the terms have been rearranged in a way which makes it apparent that the
coupling to the vector field is a “minimal coupling”, appearing through a
covariant derivative term (similarly as in electrodynamics), while the
coupling to the scalar field results in an effective mass,

$\displaystyle m^{*}\equiv
m_{N}-\frac{g_{\sigma}^{2}}{m_{\sigma}^{2}}n_{s}~{}.$ (14)

In the case of equilibrated, infinite nuclear matter the system is uniform
(translationally and rotationally invariant), so that
$\bar{\omega}_{\mu}=\delta^{0}_{\mu}\bar{\omega}_{0}$ (as the rotational
invariance of the system warrants that $\bar{\omega}_{i}=0$), which means that
all currents in the system vanish, $\bm{j}=0$. In that case, Eq. (13) reduces
to

$\displaystyle\Bigg{[}\gamma^{\mu}\bigg{(}i\partial_{\mu}-\delta_{\mu}^{0}\frac{g_{\omega}^{2}}{m_{\omega}^{2}}j_{0}\bigg{)}-m^{*}\Bigg{]}\psi(x)=0~{}.$
(15)

Note that now the nucleon field equation, Eq. (15), can be rewritten as
$\gamma^{\mu}\partial_{\mu}\psi(x)=C\psi(x)$, from which it follows that
$\psi(x)$ is a plane wave; as mentioned at the beginning of this section, one
can expect this to be a reasonable result based on the wide applicability of
the nuclear shell model. Looking for normal-mode solutions,
$\psi(x)=U(p)e^{-ip^{\nu}x_{\nu}}$, further leads to

$\displaystyle\Big{[}\bm{\gamma}\cdot\bm{p}+m^{*}\Big{]}U(p)=\gamma^{0}\bigg{(}p_{0}-\frac{g_{\omega}^{2}}{m_{\omega}^{2}}j_{0}(x)\bigg{)}U(p)~{},$
(16)

and multiplying both sides by $\gamma^{0}$ and squaring yields the eigenvalue
relation

$\displaystyle
p^{0}=\pm\sqrt{\bm{p}^{2}+{m^{*}}^{2}}+\frac{g_{\omega}^{2}}{m_{\omega}^{2}}j^{0}~{},$
(17)

in which the influence of the vector and scalar interactions on the nucleon
energy are clear. Finally, the general plane wave solution for the nucleon
field operator $\psi(x)$ can always be expanded in the basis of solutions to
the Dirac equation,

$\displaystyle\psi(x)=\frac{1}{\sqrt{V}}\sum_{\bm{p}}\sum_{s}\bigg{[}U(\bm{p},s)a_{\bm{p},s}e^{i\bm{p}\cdot\bm{x}}+V(-\bm{p},s)b^{\dagger}_{\bm{p},s}e^{-i\bm{p}\cdot\bm{x}}\bigg{]}~{},$
(18)

where the amplitudes for the normal modes, $a_{\bm{p},s}$ and
$b^{\dagger}_{\bm{p},s}$, are destruction and creation operators satisfying
the anticommutation relations

$\displaystyle\big{\\{}c_{\bm{p},s},c^{\dagger}_{\bm{p}^{\prime},s^{\prime}}\big{\\}}=\delta_{\bm{p}\bm{p}^{\prime}}\delta_{ss^{\prime}}~{},\hskip
28.45274ptc=\\{a,b\\}~{}.$ (19)

Importantly, it is possible to solve Eq. (13) without making an assumption of
uniform nuclear matter, that is for $\bm{j}\neq 0$. This is done using the
relativistic Wigner function formalism [146], and leads to the eigenvalues

$\displaystyle
p^{0}=\pm\sqrt{\bigg{(}\bm{p}-\frac{g_{\omega^{2}}}{m_{\omega}^{2}}\bm{j}(x)\bigg{)}^{2}+{m^{*}}^{2}}+\frac{g_{\omega}^{2}}{m_{\omega}^{2}}j^{0}~{}.$
(20)

Here, it is apparent that the mean-field vector interaction results in the
emergence of an effective nucleon momentum, $\Pi^{\mu}\equiv
p^{\mu}-\frac{g_{\omega}^{2}}{m_{\omega}^{2}}j^{\mu}(x)$, often called a
kinetic momentum. (It’s worth noting that while the canonical momentum
$p^{\mu}$ is the momentum subject to momentum conservation, the kinetic
momentum determines the motion of the nucleon.) The effects of vector and
meson fields can be then summarized in terms of “shifts” of the energy
(including the effective mass) and momentum of the nucleon.

From Eq. (20) it is also evident that the effects due to the mean-field vector
potential “mimic” the influence of the electromagnetic four-potential
$(\phi,\bm{A})$. This is to be expected based on the “minimal coupling” of the
mean-field vector field to the nucleon field, and we further elucidate this
fact in Appendix 12, where in particular in Section 12.C we show that the
mean-field Dirac equation, Eq. (13), leads to a Pauli-like equation in the
non-relativistic limit.

To obtain the energy-momentum tensor of the theory, one uses the Noether
theorem [164],

$\displaystyle
T^{\mu\nu}\equiv\sum_{a}\frac{\partial{\mathcal{L}}}{\partial{\left(\frac{\partial{\phi_{a}}}{\partial{x^{\mu}}}\right)}}\frac{\partial{\phi_{a}}}{\partial{x_{\nu}}}-g^{\mu\nu}\mathcal{L}~{},$
(21)

where the sum goes over the fields present in the Lagrangian; in particular,
$T_{W}^{00}\equiv\mathcal{E}_{W}$, the energy density. Using the nucleon field
expanded in terms of the mean-field Dirac equation solutions, Eq. (18), leads
to summations over contributions from given momentum states and the
corresponding probabilities of these states being occupied. In the language of
thermodynamics, this is nothing else than sums weighted with a distribution
function of the system, and thus one arrives at

$\displaystyle\mathcal{E}_{W}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\sqrt{\Pi^{2}+{m^{*}}^{2}}~{}f_{\bm{p}}+\frac{g_{\omega^{2}}}{m_{\omega}^{2}}(j^{0})^{2}+\frac{1}{2}\left(\frac{g_{\sigma}^{2}}{m_{\sigma}^{2}}\right)n_{s}^{2}-\frac{1}{2}\left(\frac{g_{\omega}^{2}}{m_{\omega}^{2}}\right)j_{\mu}j^{\mu}~{},$
(22)

where $g$ is the degeneracy of the system (for isospin-symmetric nuclear
matter, $g=g_{\textrm{spin}}g_{\textrm{iso}}=4$) and $f_{\bm{p}}$ is the
distribution function, and where the vector and scalar currents are given by

$\displaystyle
j^{\mu}\equiv\langle\bar{\psi}\gamma^{\mu}\psi\rangle=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{\Pi^{\mu}}{\Pi^{0}}~{}f_{\bm{p}}$
(23)

and

$\displaystyle
n_{s}\equiv\langle\bar{\psi}\psi\rangle=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{m^{*}}{\sqrt{\bm{\Pi}^{2}+{m^{*}}^{2}}}~{}f_{\bm{p}}~{}.$
(24)

Note that in view of Eq. (14), Eq. (24) is in fact a self-consistent equation
for the scalar density, which has to be solved numerically for a given
temperature $T$ and baryon density $n_{B}$ (or equivalently baryon chemical
potential $\mu_{B}$).

In the case of uniform nuclear matter, the Walecka energy density becomes

$\displaystyle\mathcal{E}_{W}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\sqrt{p^{2}+{m^{*}}^{2}}~{}f_{\bm{p}}+\frac{1}{2}\left(\frac{g_{\omega}^{2}}{m_{\omega}^{2}}\right)n_{B}^{2}+\frac{1}{2}\left(\frac{g_{\sigma}^{2}}{m_{\sigma}^{2}}\right)n_{s}^{2}~{}.$
(25)

Importantly, both ratios of the coupling constants to the meson masses,
$C_{\omega}=\left(\frac{g_{\omega}^{2}}{m_{\omega}^{2}}\right)$ and
$C_{\sigma}=\left(\frac{g_{\sigma}^{2}}{m_{\sigma}^{2}}\right)$, are unknown
constants (the meson masses are in fact known relatively well from the
experiment, however, the coupling constants are not, making the ratios
unconstrained). These constants are established by demanding that the binding
energy per particle,
$B=\frac{E_{W}}{A}-m_{N}=\frac{\mathcal{E}_{W}}{n_{B}}-m_{N}$, has a minimum
value of $B_{0}$ at $n_{B}=n_{0}$. In this way, the mean-field Yukawa
couplings are fitted to be used as effective couplings that capture, at least
up to the first order, many-body effects not included in the original
Lagrangian.

Having fitted the coupling constants $C_{\omega}$ and $C_{\sigma}$, one can
proceed to calculate other properties of nuclear matter. In particular, one of
the early successes of the Walecka model was not only that it allowed for a
correct description of the saturation properties of ordinary nuclear matter (a
feature that was not attained in nonrelativistic models of nucleon-nucleon
interactions), but also that it accurately identified the critical point of
the nuclear liquid-gas phase transition (see Section 1 for more details).

####  2 Quasiparticles

Within the approach used in the Walecka model, outlined in the previous
section, the starting point is a microscopic model. Employing the mean-field
approximation is equivalent to assuming that the states of the system are
single-particle wave functions, with dispersion relations “shifted” due to
uniform background meson fields, Eq. (20). Based on this picture, one can
obtain macroscopic observables, such as the energy density, which are then
used to fit the coupling constants such that they reproduce macroscopic
properties of nuclear matter. The use of plane waves (corresponding to freely
propagating particles) with some of the free-particle properties “shifted” due
to the interactions is equivalent to using “quasiparticles”, a concept that we
introduce in this section.

The notion of quasiparticles originates from the many-body theory, which aims
at describing phenomena occurring in systems with many degrees of freedom.
This includes, for example, the behavior of metals or the behavior of liquid
helium. In such systems there are macroscopic numbers of elementary
constituents, such as ions and electrons of a crystal lattice, or liquid 3He
atoms. These constituents interact with each other in some very complex ways,
which are not only challenging to capture on the level of particle-particle
interactions, but would be absolutely impossible to handle when considering a
macroscopic sample of the system. However, it can be argued that such detailed
description of all interactions occurring in a given system is simply
unnecessary for capturing the macroscopic properties of that system, which
largely reflect average properties of the interactions in the system as a
whole.

This average behavior of the system can be described by means of collective
excitations. One of the most illustrative examples of collective excitations
is provided by vibrations of a lattice of $N$ atoms. If the oscillations are
assumed to be small, such lattice can be represented by a set of coupled
harmonic oscillators. Through introducing normal coordinates one can arrive at
a description in which the system is represented by $3N$ linear oscillators,
each described by a given eigenfrequency $\omega_{i}$. From quantum mechanics,
the energy of such a system is then given by the formula

$\displaystyle E=\sum_{i=1}^{3N}\Big{(}n_{i}+\frac{1}{2}\Big{)}\omega_{i}~{},$
(26)

and taking different sets of values of $n_{i}$ leads to constructing the
energy spectrum of the system. In the case where the potential energy of the
lattice takes into the account anharmonic terms (which become more important
with the amplitude of the vibrations, that is with increasing temperature),
the system will start to exhibit non-zero probability of transitioning between
states with different values of the $n_{i}.$

On the other hand, lattice vibrations can be always decomposed into a sum of
monochromatic waves propagating in the crystal. Based on the correspondence
principle of quantum mechanics, a plane wave can be associated with a moving
“particle”: the wave vector $\bm{k}$ is directly related to the momentum of
such a particle, and the frequency $\omega(\bm{k})$ is related to its energy.
Any excited state of the lattice can be described as a set of such “particles”
moving throughout the volume of the considered solid, and so the energy levels
of the system will be very similar to energy states of an ideal gas of these
“particles”, called “phonons” (a name first introduced by Igor Tamm in 1932,
in an analogy to “photon” and based on a Greek word meaning “sound”). In
consequence, one can reinterpret the number $n_{i}$ as the number of phonons
in the state $i$. Because $n_{i}$ can take any values, evidently phonons are
subject to Bose-Einstein statistics (this holds also in the case when the
lattice atoms have half-odd-integer spin). The number of phonons increases
with temperature, which leads to an increased role of phonon-phonon
interactions (in analogy to an increasing importance of anharmonic terms). In
particular, one can reinterpret the finite probability of the system going
from one energy state to another, characterized by different values of
$\\{n_{i}\\}$, in terms of different interactions between the phonons, such as
scattering, creation, or decay.

Based on the above, one can see that the excitations of a lattice, energy
levels of which are a result of an incomprehensible number of interactions
between a macroscopic number of particles, can to a good approximation be
constructed out of energy levels of free “particles”, which are entirely
independent for sufficiently low temperatures and are weakly-interacting for
higher temperatures. These “particles” are often called “collective
excitations” or “quasiparticles” (the former is often used when the
“particles” are governed by the Bose-Einstein statistics, while the latter
when they are governed by the Fermi-Dirac statistics, although this is by no
means a strict rule, and we will not always follow it). Intuitively,
quasiparticles can be understood as emergent phenomena occurring when a
microscopically complex system of “real” particles can be described as if it
was made of different, weakly interacting “quasiparticles” in free space,
which are characterized by a particular value of the momentum $\bm{p}$ and the
corresponding dispersion relation $\varepsilon(\bm{p})$. A widely known
example of a system that can be described using this concept is a
semiconductor, in which the behavior of an electron interacting with the ion
lattice can be described as a motion of a free electron with a different,
“effective” mass.

One needs to remember that the collective excitations are born out of the
collective behavior of the studied system, and as such do not exist “on their
own”; rather, their properties are determined by the state of the system.
Importantly, the number of quasiparticles can, but need not be equal to the
total number of “real” particles considered. Furthermore, because of a finite
probability of a transition from one state of the system to another, each
described by a different set of quasiparticles, the quasiparticles are
characterized by a finite decay width. Finally, it should be stressed that the
concept of collective excitations only makes sense for length scales larger
then some fundamental scale of the problem such as the lattice spacing; for
example, if one resolves the individual microscopic motions of the lattice
ions, the description in terms of phonons becomes meaningless.

In view of the above, it is clear that the nucleon fields described in the
mean-field Walecka model are in fact quasiparticles. We can think of them as
nucleons “dressed” by their interactions with mesons, which results in a
“quasi-nucleon” field characterized by an effective, density- and temperature-
dependent mass, and an effective (kinetic) momentum, both reflecting the fact
that each nucleon is interacting with baryon currents created by all other
nucleons in the system. The next section is devoted to introducing a
quasiparticle theory for spin-$\frac{1}{2}$ particles which indeed can be used
to reproduce Walecka model results.

####  3 Landau Fermi-liquid theory

A system of non-interacting particles of spin-$\frac{1}{2}$ (which are
therefore governed by the Fermi-Dirac statistics) is called a Fermi gas. A
system of interacting spin-$\frac{1}{2}$ particles at low temperatures is
called a Fermi liquid. Examples of such systems include liquid 3He, electrons
in metals, or nuclear matter. One of the most fruitful phenomenological
theories of Fermi liquids was developed by Lev Landau between 1956 and 1959
[165, 166, 167], and is nowadays known as the Landau Fermi-liquid theory.

Landau’s theory is based on an assumption that the excitation spectrum of a
Fermi liquid has a structure similar to that of an ideal Fermi gas. This is a
strong assumption on the dominating type of interactions occurring between the
particles of the liquid which must preserve their quantum statistics (an
example of a system in which the opposite is the case is deuterium atoms, in
which nucleons, governed by the Fermi-Dirac statistics, interact in a way
leading to the formation of molecules governed by the Bose-Einstein
statistics). Furthermore, the theory assumes that the number of particles in
the liquid and the number of quasiparticles are equal, and the Fermi momentum
of the system is related to the particle density of the liquid, $n$, in a way
completely mirroring the ideal gas case,

$\displaystyle p_{F}=\left(\frac{6\pi^{2}n}{g}\right)^{1/3}~{},$ (27)

where $g$ is the degeneracy (the above equation follows from the fact that the
particle density $n$ is given by the momentum integral of the distribution
function $f_{\bm{p}}$, $n=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}f_{\bm{p}}$, which
for the ideal Fermi gas at $T=0$, where the system occupies all momentum
states from 0 to $p_{F}$ and no momentum states above $p_{F}$, reduces to
$n_{B}=\frac{g}{2\pi^{2}}\int_{0}^{p_{F}}dp~{}p^{2}$).

The bulk behavior of a system of quasiparticles is assumed to be described by
a quasiparticle distribution function, constructed based on the one-to-one
correspondence between quasiparticles and “real” particles. This
correspondence can be introduced formally as follows: One begins by
considering an ideal (nonrelativistic) Fermi gas, in which the dispersion
relation (neglecting spin) is

$\displaystyle\varepsilon_{\bm{p}}^{\textrm{free}}=\frac{\bm{p}^{2}}{2m}~{}.$
(28)

The state of the system as whole can be specified by giving the number of
particles $N_{\bm{p}}=\\{0,1\\}$ in each of the single-particle states defined
by a specific value of the momentum $\bm{p}$. Thus for example in the ground
state, each of the states with momenta less than the Fermi momentum, $p_{F}$,
is occupied ($N_{\bm{p}}=1$), and all other states are empty ($N_{\bm{p}}=0$).
One can then imagine that interactions in the system are slowly turned on in
such a way that the process is adiabatic. Quantum mechanics shows that while
such an adiabatic change will lead to a distortion of the energy levels, it
will preserve their number. This means that the distribution function
$N_{\bm{p}}$, while also smoothly distorted, preserves its functional form.
Now, however, the dispersion relation $\varepsilon_{\bm{p}}^{\textrm{int}}$
takes interactions into the account, and so it differs from that of a free
particle, Eq. (28). (Note that the dependence of the construction of the
quasiparticle distribution function on the one-to-one correspondence between
quasiparticles and “real” particles means that the formalism is not
appropriate for describing phenomena in which the number of particles in the
system changes as a result of the interactions, such as formation or
dissolution of bound states.) For describing macroscopic properties of a Fermi
liquid, it is sufficient to use a mean or smoothed quasiparticle distribution
function, often denoted by $f_{\bm{p}}$, which is an average of $N_{\bm{p}}$
over a group of neighboring single-particle states. While $N_{\bm{p}}$ is a
discrete function of $\bm{p}$, $f_{\bm{p}}$ and is a continuous function of
$\bm{p}$.

It follows from the above that the quasiparticle distribution function
$f_{\bm{p}}$ has the same functional form as the Fermi-Dirac distribution
function; a more formal derivation of this fact can be found in Appendix 13.
Importantly, the Fermi distribution function satisfies the constraint that the
quasiparticles are subject to the Pauli exclusion principle, which is an
essential ingredient of the theory. At low temperatures ($T\ll p_{F}$), all
states deep in the Fermi sphere are occupied, so that available states only
appear in the vicinity of the Fermi surface. Excitation from one of the states
from deep within the Fermi sphere to one of the available states is then
strongly suppressed by a factor of $1/\Delta E$, where $\Delta E$ is the
energy difference between the two considered states. In consequence, only the
states in the vicinity of the Fermi surface can actively contribute to the
dynamics of the system, while the rest of the states constitute a necessary
(from the point of view of Pauli blocking) but kinematically suppressed
“background”.

As the next step in developing the theory, Landau assumed that the interaction
between the quasiparticles takes the form of a self-consistent field, whose
effect on any given quasiparticle is considered to be produced by the
remaining quasiparticles. He then proceeded to construct a theory of
quasiparticles in which energy and momentum were conserved. To this end, he
utilized the Boltzmann equation for quasiparticles,

$\displaystyle\frac{df_{\bm{p}}}{dt}=\frac{\partial{f_{\bm{p}}}}{\partial{t}}+\frac{d\bm{x}}{dt}\frac{\partial{f_{\bm{p}}}}{\partial{\bm{x}}}+\frac{d\bm{p}}{dt}\frac{\partial{f_{\bm{p}}}}{\partial{\bm{p}}}=\left(\frac{d{f_{\bm{p}}}}{d{t}}\right)_{\textrm{coll}}~{},$
(29)

where $f_{\bm{p}}$ denotes the quasiparticle distribution function and
$\big{(}d{f_{\bm{p}}}/d{t}\big{)}_{\textrm{coll}}$ is the change in the
distribution function due to collisions, often referred to as the “collision
term” or the “collision integral” (a more thorough discussion of the Boltzmann
equation can be found in Chapter 3 and Appendix 23). From Hamilton’s equations
we know that $d\bm{x}/dt\equiv\partial H_{(1)}/\partial\bm{p}$ and
$d\bm{p}/dt=-\partial H_{(1)}/\partial\bm{x}$, where $H_{(1)}$ is the single-
particle Hamiltonian. Equating $H_{(1)}$ with the energy of a quasiparticle
$\varepsilon_{\bm{p}}$ yields

$\displaystyle\frac{\partial{f_{\bm{p}}}}{\partial{t}}-\frac{\partial{\varepsilon_{\bm{p}}}}{\partial{p_{i}}}\frac{\partial{f_{\bm{p}}}}{\partial{x^{i}}}+\frac{\partial{\varepsilon_{\bm{p}}}}{\partial{x_{i}}}\frac{\partial{f_{\bm{p}}}}{\partial{p^{i}}}=\left(\frac{d{f_{\bm{p}}}}{d{t}}\right)_{\textrm{coll}}~{},$
(30)

where we utilize the Einstein summation convention as well as co- and
contravariant vector notation (as summarized in Appendix 7), the latter of
which leads to a change in signs as compared with Eq. (29).

To ensure momentum conservation, Landau wanted to obtain an expression for the
continuity of the momentum density (leading to the conservation of the
quasiparticle momentum by assuming that the momentum density flux is zero on
the surface of the integrated volume), which formally can be written as

$\displaystyle\frac{\partial{}}{\partial{t}}~{}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{i}~{}f_{\bm{p}}+\frac{\partial{}}{\partial{x^{k}}}~{}\Pi^{ik}=0~{},$
(31)

where $\Pi^{ik}$ is the momentum density flux tensor of an as of yet
unspecified form. To arrive at the momentum conservation law through the
Boltzmann equation, Eq. (30), one multiplies both sides of the equation by
$p^{k}$ and integrates over all possible momenta,

$\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{k}~{}\frac{\partial{f_{\bm{p}}}}{\partial{t}}-g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{k}\frac{\partial{\varepsilon_{\bm{p}}}}{\partial{p_{i}}}~{}\frac{\partial{f_{\bm{p}}}}{\partial{x^{i}}}+g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{k}\frac{\partial{\varepsilon_{\bm{p}}}}{\partial{x_{i}}}~{}\frac{\partial{f_{\bm{p}}}}{\partial{p^{i}}}=0~{};$
(32)

here, the right hand side vanishes as the integral
$g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{k}~{}\big{(}\frac{d{f_{\bm{p}}}}{d{t}}\big{)}_{\textrm{coll}}$
is equal zero due to the conservation of momentum in collisions (see Appendix
14 for more details). Using the product rule for derivatives and the fact that
the distribution function $f_{\bm{p}}$ vanishes at the boundary,
$f_{\bm{p}}\big{(}|\bm{p}|=+\infty\big{)}=0$, it is possible to rewrite Eq.
(32) as (a similar calculation is performed in detail in Appendix 15.B)

$\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{k}~{}\frac{\partial{f_{\bm{p}}}}{\partial{t}}-\frac{\partial{}}{\partial{x^{i}}}~{}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{k}\frac{\partial{\varepsilon_{\bm{p}}}}{\partial{p_{i}}}~{}f_{\bm{p}}$
$\displaystyle\hskip
28.45274pt-~{}\frac{\partial{}}{\partial{x_{i}}}g^{k}_{~{}i}~{}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}~{}f_{\bm{p}}+g^{k}_{~{}i}~{}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}~{}\frac{\partial{f_{\bm{p}}}}{\partial{x_{i}}}=0~{},$
(33)

where we have used
$\frac{\partial{p^{k}}}{\partial{p^{i}}}=\delta^{k}_{i}=g^{kj}g_{ji}=g^{k}_{~{}i}$.
The first three terms in the above equation are already of the from required
by Eq. (31), however, the same cannot be achieved trivially for the last term.
The only way to write that term as a divergence of some quantity is to assume
that the integrand is a functional differential of some functional
$\mathcal{E}$, that is to assume that

$\displaystyle\delta\mathcal{E}[f_{\bm{p}}]=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}~{}\delta
f_{\bm{p}}~{}.$ (34)

If this is the case, then one achieves the goal of bringing Eq. (33) to the
form of Eq. (31) (which is equivalent to establishing that momentum density is
conserved), with

$\displaystyle\Pi^{ik}=-g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{i}\frac{d\varepsilon_{\bm{p}}}{dp_{k}}~{}f_{\bm{p}}+g^{ik}~{}\left[\mathcal{E}-g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}~{}f_{\bm{p}}\right]~{}.$
(35)

Looking at Eq. (34), it is quite natural to assume that $\mathcal{E}$ is the
energy density of the system. Then the quasiparticle energy
$\varepsilon_{\bm{p}}$, which based on the above derivation is a functional
derivative of the energy density,

$\displaystyle\varepsilon_{\bm{p}}\equiv\frac{\delta\mathcal{E}}{\delta
f_{\bm{p}}}~{},$ (36)

can be interpreted as the change in the energy density of the system due to an
addition of a quasiparticle with a momentum $\bm{p}$. So defined, the
quasiparticle energy naturally contains information about the interactions in
the liquid. Consequently, the energy of a quasiparticle includes potential
energy terms arising due to the interactions with all other quasiparticles,
and from this it follows that the energy of the system as a whole is not equal
to the sum of the energies of individual quasiparticles. (Note that this is
similar to the case of a classical electrostatic system, where in the
calculation of the total electrostatic potential energy $U$ one needs to avoid
summing over the same particle pair twice,
$U=\sum_{i=0}^{N}\sum_{j>i}U_{ij}=\frac{1}{2}\sum_{i=0}^{N}\sum_{j\neq
i}U_{ij}$, where $U_{ij}=(4\pi\varepsilon_{0})^{-1}q_{i}q_{j}/r_{ij}$ is the
potential energy of the particle $i$ due to the presence of the particle $j$,
with $q_{i}$ and $q_{j}$ denoting charges of the $i$-th and the $j$-th
particle and $r_{ij}$ denoting the distance between these particles. Defining
the overall potential energy of the $i$-th particle as $U_{i}=\sum_{j\neq
i}U_{ij}$, the potential energy of the system can be rewritten as
$U=\frac{1}{2}\sum_{i}U_{i}$. Then the total energy of the system can be
expressed as
$E=\sum_{i=0}^{N}\big{(}K_{i}+U_{i}\big{)}-\frac{1}{2}\sum_{i}U_{i}$, where
$K_{i}$ is the kinetic energy contribution from the $i$-th particle and
$K_{i}+U_{i}=\varepsilon_{i}$, the total energy of that particle. The term
$-\frac{1}{2}\sum_{i}U_{i}$ is known as the “counterterm”, and its purpose is
to cancel out the double-counting of potential energy terms that occurs in the
sum over the energies of all particles. In a complete analogy, the energy
density of a system of quasiparticles always has the form
$\mathcal{E}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}+\mathcal{E}_{\textrm{counter}}$,
where $\mathcal{E}_{\textrm{counter}}$ is the counterterm.)

With Eq. (36) established, it is straightforward to show that the developed
Landau Fermi-liquid theory conserves not only the quasiparticle momentum, but
also the number of quasiparticles and the quasiparticle energy. The
conservation of the number of quasiparticles is obtained by integrating both
sides of Eq. (30) over all momenta; here, the right-hand side of the obtained
expression vanishes due to the conservation of the number of particles in
collisions (see Appendix 14 for more details), and after some algebra we get

$\displaystyle\frac{\partial{}}{\partial{t}}\int\frac{d^{3}p}{(2\pi)^{3}}~{}f-\frac{\partial{}}{\partial{x^{i}}}\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{\partial{\varepsilon_{\bm{p}}}}{\partial{p_{i}}}~{}f=0~{}.$
(37)

Using
$-\partial{\varepsilon_{\bm{p}}}/\partial{p_{i}}=\partial{\varepsilon_{\bm{p}}}/\partial{p^{i}}\equiv
v^{i}$, where $v^{i}$ is the velocity of a quasiparticle, we immediately
obtain

$\displaystyle\frac{\partial{}}{\partial{t}}\int\frac{d^{3}p}{(2\pi)^{3}}~{}f+\frac{\partial{}}{\partial{x^{i}}}\int\frac{d^{3}p}{(2\pi)^{3}}~{}v^{i}~{}f=0~{},$
(38)

where identifying $\int\frac{d^{3}p}{(2\pi)^{3}}~{}f$ and
$\int\frac{d^{3}p}{(2\pi)^{3}}~{}v^{i}~{}f$ with the number density $n$ and
the number current $j^{i}$, respectively, reveals Eq. (38) as the continuity
equation for the quasiparticle number density. Similarly, the law of the
conservation of the quasiparticle energy is obtained by multiplying both sides
of Eq. (30) by $\varepsilon_{\bm{p}}$ and integrating over all momenta. As
before, the right-hand side of the obtained expression vanishes (see Appendix
14 for more details), leading to

$\displaystyle\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}~{}\frac{\partial{f}}{\partial{t}}-\frac{\partial{}}{\partial{x^{i}}}\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}\frac{\partial{\varepsilon_{\bm{p}}}}{\partial{p_{i}}}~{}f=0~{}.$
(39)

Using Eq. (34) to rewrite the first term and the definition of the
quasiparticle velocity $v^{i}$, the above equation becomes

$\displaystyle\frac{\partial{\mathcal{E}}}{\partial{t}}+\frac{\partial{}}{\partial{x^{i}}}\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}v^{i}~{}f=0~{},$
(40)

where identifying
$\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}v^{i}~{}f$ as the
quasiparticle energy density flux yields the continuity equation for the
quasiparticle energy density.

Let us further discuss two features of the Landau theory.

First, because the system is interacting, the concept of quasiparticles only
has a well-defined meaning for energies close to the Fermi energy
$\varepsilon_{F}$, or equivalently momenta close to the Fermi momentum
$p_{F}$. This is because for excitations far away from the Fermi surface,
processes such as scattering into other states become kinematically allowed.
In the “real” system, these correspond to a finite probability for the system
to transition from one state to another, leading to damping (finite decay
width) of the excitation that the quasiparticles are to represent. This means
that the description in terms of quasiparticles is strictly valid only for
temperatures which are low in comparison with the Fermi momentum. (Note,
however, a recent work which shows that the expectation of broadening and
disappearance of quasiparticle states at high energies, correct for weak
interactions, is invalid when strong interactions are considered [168].)

Second, one notes that the above derivation assumes a uniform distribution of
quasiparticles in space. Such an assumption can be considered valid as long as
any spatial non-uniformity occurs at distances not probed by quasiparticles,
that is distances that exceed the quasiparticle wavelength. At the Fermi
surface, the momenta are close to the Fermi momentum
$p_{F}=\big{[}\left({6\pi^{2}N}/{gV}\right)\big{]}^{1/3}$, and consequently
the quasiparticle wavelengths are of order $\left({V}/{N}\right)^{1/3}$, which
is the average inter-particle spacing. It follows that the system only needs
to be approximately uniform on a length scale of the average distance between
the particles, which in practice doesn’t impose any restrictions for systems
at low temperatures.

How does one use the Landau-Fermi liquid theory in practice? First, one
postulates the energy density of the system, and in particular its dependence
on the distribution function, $\mathcal{E}=\mathcal{E}[f_{\bm{p}}]$. (Note
that this is very much different than the strategy used in microscopic
approaches, which start from microscopic interactions and then endeavor to
calculate macroscopic observables such as the energy density or pressure.)
Next, one calculates the quasiparticle energy using Eq. (36). Then one can
also calculate the second functional derivative of the energy density,

$\displaystyle\phi\equiv\frac{\delta\mathcal{E}}{\delta
f_{\bm{p}^{\prime}}\delta f_{\bm{p}}}~{},$ (41)

which can be interpreted as a quasiparticle interaction energy and generally
depends on the two quasiparticles’ momenta, $\bm{p}$ and $\bm{p^{\prime}}$, as
well as their spins, $s$ and $s^{\prime}$. By neglecting the spin dependence
and assuming that both quasiparticles are close to the Fermi surface,
$|\bm{p}|\approx|\bm{p}^{\prime}|\approx p_{F}$, one has
$\phi=\phi(p_{F},\theta)$, where
$\theta=\measuredangle(\bm{p},\bm{p}^{\prime})$. It is then possible to expand
$\phi$ in spherical harmonics, and coefficients of this expansion are known as
the Landau coefficients; often, it is enough to know the two first
coefficients of the expansion, $F_{0}$ and $F_{1}$. Many properties of matter,
such as the quasiparticle velocity at the Fermi surface or the zero sound, can
then be calculated and expressed in terms of the Landau coefficients. One can
utilize experimental measurements of some of these quantities to obtain
phenomenological values of $F_{0}$ and $F_{1}$, and then use these to predict
the values of other experimental observables. In general, this has been a very
successful approach in describing phenomena occurring in Fermi liquids,
including the behavior of liquid 3He [169], electrons in metals [170], and
nuclear matter [171].

In a 1975 paper, Gordon Baym and Siu Chin [172] generalized the Landau Fermi-
liquid theory to include relativistic effects. Soon after, the formalism was
used by Tetsuo Matsui [173] as an alternative way of obtaining some of the
results predicted by the Walecka model, such as the incompressibility or the
speed of sound of nuclear matter. This further underscored the equivalency
between the Landau Fermi-liquid theory and the Walecka model field theory in
the mean-field approximation. In particular, both of these approaches are
appropriate in sufficiently dense regimes, where excitations are low-energetic
in comparison with the Fermi surface and the collective effects due to a self-
consistent background are more important than contributions from fluctuations
in the fields. For nuclear matter, the Fermi momentum at saturation density is
$p_{F}(n_{0})\approx 260\ \textrm{MeV}$, resulting in a relatively substantial
range of temperatures in which the description in terms of quasiparticles is
appropriate.

### 2\. Relativistic Landau Fermi-liquid theory for hadronic transport

In the previous section, we reviewed two phenomenological methods of
describing nuclear matter: First, we introduced the Walecka model as an
example of a Lagrangian-based, self-consistent approach at the mean-field
level. Then we introduced the Landau Fermi-liquid theory, in which the
relevant physics is entirely encoded in the postulated energy density and its
variations with respect to the distribution function governing the system.
Both methods utilize the concept of quasiparticles, which can be thought of as
“dressed” real particles whose properties, defined through self-consistent
equations, incorporate the effects of interactions. Based on that, we
concluded that the two approaches describe the same physics. The advantage of
the Landau Fermi-liquid theory is that one can freely postulate any energy
density functional (as long as it transforms like the $00$-component of the
energy-momentum tensor under a Lorentz boost and includes correct
counterterms) and rely on the formalism to yield a description of the studied
system in terms of quasiparticles whose dynamics conserve particle number,
energy, and momentum. At the same time, the advantage of the Walecka model is
that it provides a clear example of what such an energy density functional
could look like based on a Lagrangian incorporating nucleon-meson
interactions.

In the following, we will draw on both approaches to form a flexible equation
of state for nuclear matter. We will use the energy density stemming from the
Walecka model as an inspiration for postulating a generalized energy density
in a model with an arbitrary number of scalar- and vector-type interactions.
We will then obtain the single-particle energies of the quasiparticles using
the Landau Fermi-liquid theory, from which all other properties of the system
will be derived. Unlike in the original Landau Fermi-liquid approach, which
uses the properties of the interactions between the quasiparticles (as encoded
in the Landau parameters) to constrain the free parameters of the model, we
will instead fit these parameters to reproduce a set of chosen bulk properties
of nuclear matter (this will take place in Chapter 2). In this regard, the
formalism we develop is similar to a large class of effective approaches to
describing nuclear matter using self-consistent models based on the density
functional theory (DFT) [174, 175]. Such models are a starting point for
numerous Skyrme-like potentials of varying degrees of complexity which are
successfully applied in low-energy nuclear physics [176].

### 3\. Relativistic density functional equation of state for nuclear matter

The contents of this section are a generalization of the model presented in
Ref. [177].

####  1 Formalism

We adopt the relativistic Landau Fermi-liquid theory [172] with vector- and
scalar-density–dependent interactions as the basis for constructing a vector
and scalar density functional (VSDF) model of dense nuclear matter EOS. To
simplify the notation, we will introduce a VSDF model with a single vector-
current- and a single scalar-current–dependent interaction term, however, it
is straightforward to generalize to a model with multiple interaction terms of
both kinds, which we do in the next subsection. Some of the details of the
derivation will be referred to Appendix 15.

We introduce the energy density of a system composed of one species of
fermions of rest mass $m_{0}$, interacting through a single mean-field vector
interaction term and a single mean-field scalar interaction term,

$\displaystyle\mathcal{E}_{\textrm{v1s1}}$ $\displaystyle=$ $\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\epsilon_{\textrm{kin}}^{*}~{}f_{\bm{p}}+C_{1}\big{(}j_{\mu}j^{\mu}\big{)}^{\frac{b_{1}}{2}-1}\big{(}j^{0}\big{)}^{2}$
(42) $\displaystyle\hskip
14.22636pt-~{}g^{00}C_{1}\left(\frac{b_{1}-1}{b_{1}}\right)\big{(}j_{\mu}j^{\mu}\big{)}^{\frac{b_{1}}{2}}+G_{1}\left(\frac{d_{1}-1}{d_{1}}\right)n_{s}^{d_{1}}~{}.$

Here, the subscript “v1s1” underscores the fact that we have one vector and
one scalar interaction term, $g$ is the degeneracy,
$\epsilon_{\textrm{kin}}^{*}$ is the kinetic energy of a single quasiparticle,

$\displaystyle\epsilon_{\textrm{kin}}^{*}=\sqrt{\bigg{(}\bm{p}-C_{1}\big{(}j_{\mu}j^{\mu}\big{)}^{\frac{b_{1}}{2}-1}\bm{j}\bigg{)}^{2}+{m^{*}}^{2}}~{},$
(43)

$\bm{j}$ and $j^{0}$ are the spatial and temporal components of the number
current $j^{\mu}$, given by

$\displaystyle\bm{j}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{\bm{p}-C_{1}\big{(}j_{\mu}j^{\mu}\big{)}^{\frac{b_{1}}{2}-1}\bm{j}}{\epsilon_{\textrm{kin}}^{*}}~{}f_{\bm{p}}$
(44)

and

$\displaystyle j^{0}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}f_{\bm{p}}~{},$ (45)

respectively, $f_{\bm{p}}$ is the quasiparticle distribution function,
$m^{*}=m^{*}(x)$ is an effective mass given by a self-consistent equation,

$\displaystyle m^{*}=m_{0}-G_{1}n_{s}^{d_{1}-1}~{},$ (46)

where the number density $n_{s}$ is given by

$\displaystyle
n_{s}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{m^{*}}{\epsilon_{\textrm{kin}}^{*}}~{}f_{\bm{p}}~{},$
(47)

$g^{00}$ is the $00$-component of the metric tensor (see Appendix 7 for the
metric convention), and finally $C_{1}$, $G_{1}$, $b_{1}$, $d_{1}$ are
constants specifying the interaction, as of yet undetermined. The coefficients
and powers of the interaction terms have number indeces in anticipation of
adding more interaction terms. We note that Eq. (42) reduces to the Walecka
model energy density, Eq. (22), for $b_{1}=2$ and $d_{1}=2$ (and an
appropriate identification of the coefficients of the interaction terms).

The quasiparticle energy $\varepsilon_{\bm{p}}$ is obtained from the
definition (details of the calculation can be found in Appendix 15.A),

$\displaystyle\varepsilon_{\bm{p}}^{*}\equiv\frac{\delta\mathcal{E}_{\textrm{v1s1}}}{\delta
f_{\bm{p}}}=\epsilon_{\textrm{kinetic}}^{*}+C_{1}\big{(}j_{\mu}j^{\mu}\big{)}^{\frac{b_{1}}{2}-1}j_{0}~{}.$
(48)

Also here setting $b_{1}=2$ and $d_{1}=2$ reduces the quasiparticle energy,
Eq. (48), to the expression known from the Walecka model, Eq. (20).

To simplify the notation, we introduce a vector field,

$\displaystyle
A^{\lambda}_{(1)}=A^{\lambda}(x;C_{1},b_{1})=C_{1}\big{(}j_{\mu}j^{\mu}\big{)}^{\frac{b_{1}}{2}-1}j^{\lambda}~{},$
(49)

which allows us to concisely write

$\displaystyle\varepsilon_{\bm{p}}^{*}=\sqrt{\Big{(}\bm{p}-\bm{A}_{(1)}\Big{)}^{2}+{m^{*}}^{2}}+A_{0}^{(1)}$
(50)

and

$\displaystyle\mathcal{E}_{v1s1}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}^{*}~{}f_{\bm{p}}-g^{00}\left(\frac{b_{1}-1}{b_{1}}\right)A_{\lambda}^{(1)}j^{\lambda}+G_{1}\left(\frac{d_{1}-1}{d_{1}}\right)n_{s}^{d_{1}}~{}.$
(51)

For clarity, in the following derivation we will temporarily suppress the
interaction term index $(1)$ and refer to this variable simply as
$A^{\lambda}(x)$.

Quite generally, the quasiparticle energy is equivalent to the single-particle
Hamiltonian, $\varepsilon_{\bm{p}}^{*}=H_{(1)}$. Given the single-particle
Hamiltonian, the equations of motion governing the evolution of quasiparticles
follow immediately from Hamilton’s equations,

$\displaystyle\frac{dx^{i}}{dt}$ $\displaystyle\equiv$
$\displaystyle-\frac{\partial{H_{(1)}}}{\partial{p_{i}}}=-\frac{\partial{\varepsilon_{\bm{p}}^{*}}}{\partial{p_{i}}}~{},$
(52) $\displaystyle\frac{dp^{i}}{dt}$ $\displaystyle\equiv$
$\displaystyle\frac{\partial{H_{(1)}}}{\partial{x_{i}}}=\frac{\partial{\varepsilon_{\bm{p}}^{*}}}{\partial{x_{i}}}~{}.$
(53)

where we utilize the co- and contravariant vector notation as summarized in
Appendix 7. Inserting Eqs. (52) and (53) into the Boltzmann equation (where we
use the Einstein summation convention) gives

$\displaystyle\frac{\partial{f_{\bm{p}}}}{\partial{t}}-\frac{\partial{\varepsilon_{\bm{p}}^{*}}}{\partial{p_{i}}}\frac{\partial{f_{\bm{p}}}}{\partial{x^{i}}}+\frac{\partial{\varepsilon_{\bm{p}}^{*}}}{\partial{x_{i}}}\frac{\partial{f_{\bm{p}}}}{\partial{p^{i}}}=\left(\frac{d{f_{\bm{p}}}}{d{t}}\right)_{\textrm{coll}}~{},$
(54)

where $\big{(}d{f_{\bm{p}}}/d{t}\big{)}_{\textrm{coll}}$ is the change in the
distribution function due to collisions occurring in the system (a thorough
discussion of the Boltzmann equation can be found in Chapter 3 and Appendix
23). Multiplying both sides of Eq. (54) by
$X=\\{1,\varepsilon_{\bm{p}}^{*},p^{j}\\}$ and integrating over all momenta,
$g\int\frac{d^{3}p}{(2\pi)^{3}}$, yields the conservation laws for
quasiparticle number ($X=1$), energy ($X=\varepsilon_{\bm{p}}^{*}$), and
momentum ($X=p^{j}$). In particular, terms entering the energy and momentum
conservation laws can be identified with the components of the energy-momentum
tensor (see Appendix 15.B for the derivation), given by

$\displaystyle T^{00}=\mathcal{E}~{},$ (55) $\displaystyle
T^{0i}=-g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}^{*}\frac{\partial{\varepsilon_{\bm{p}}^{*}}}{\partial{p_{i}}}~{}f_{\bm{p}}~{},$
(56) $\displaystyle
T^{i0}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{i}~{}f_{\bm{p}}~{},$ (57)
$\displaystyle
T^{ij}=-g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{i}\frac{\partial{\varepsilon_{\bm{p}}^{*}}}{\partial{p_{j}}}~{}f_{\bm{p}}+g^{ij}\left(\mathcal{E}-g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}^{*}~{}f_{\bm{p}}\right)~{}.$
(58)

One can show [172] that such obtained $T^{\mu\nu}$ has the correct
transformation properties under a Lorentz boost. Additionally, the
conservation of energy and momentum, $\partial_{\nu}T^{\mu\nu}=0$, is ensured
by construction.

Inserting the expression for the quasiparticle energy, Eq. (50), into Eqs.
(52) and (53) yields

$\displaystyle\frac{dx^{i}}{dt}$ $\displaystyle=$
$\displaystyle-\frac{\partial{\varepsilon_{\bm{p}}^{*}}}{\partial{p_{i}}}=\frac{p^{i}-A^{i}}{\epsilon_{\textrm{kin}}^{*}}~{},$
(59) $\displaystyle\frac{dp^{i}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\partial{\varepsilon_{\bm{p}}^{*}}}{\partial{x_{i}}}=\frac{\big{(}p^{k}-A^{k}\big{)}\frac{\partial{A_{k}}}{\partial{x_{i}}}+m^{*}\frac{\partial{m^{*}}}{\partial{x_{i}}}}{\epsilon_{\textrm{kin}}^{*}}+\frac{\partial{A_{0}}}{\partial{x_{i}}}~{},$
(60)

while the energy-momentum tensor becomes

$\displaystyle T^{00}=\mathcal{E}_{v1s1}~{},$ (61) $\displaystyle
T^{0i}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}^{*}~{}\frac{p^{i}-A^{i}}{\epsilon_{\textrm{kin}}^{*}}~{}f_{\bm{p}}~{},$
(62) $\displaystyle
T^{i0}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{i}~{}f_{\bm{p}}~{},$ (63)
$\displaystyle
T^{ij}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{i}~{}\frac{p^{j}-A^{j}}{\epsilon_{\textrm{kin}}^{*}}~{}f_{\bm{p}}+g^{ij}\left(\mathcal{E}-g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\varepsilon_{\bm{p}}^{*}~{}f_{\bm{p}}\right)~{}.$
(64)

In particular, we can immediately check that the energy current $T^{0i}$ is
equal to the momentum density $T^{i0}$, as required by hydrodynamics,

$\displaystyle T^{0i}$ $\displaystyle=$ $\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\Big{[}\epsilon_{\textrm{kinetic}}^{*}+A_{0}\Big{]}\left[\frac{p^{i}-A^{i}}{\epsilon_{\textrm{kinetic}}^{*}}\right]~{}f_{\bm{p}}$
(65) $\displaystyle=$ $\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{i}~{}f_{\bm{p}}-A^{i}~{}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}~{}f_{\bm{p}}+A_{0}~{}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{p^{i}-A^{i}}{\epsilon_{\textrm{kinetic}}^{*}}~{}f_{\bm{p}}$
$\displaystyle=$ $\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{i}~{}f_{\bm{p}}-A^{i}j^{0}+A_{0}j^{i}$
$\displaystyle=$ $\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}p^{i}~{}f_{\bm{p}}=T^{i0}~{},$

where we have used the definition of vector current, Eq. (44), the definition
of vector density, Eq. (45), and the fact that from the definition of
$A^{\lambda}$, Eq. (49), we have $-A^{i}j^{0}+A_{0}j^{i}=0$.

For completeness, we note that the quasiparticle number conservation law as
obtained by taking the zeroth moment ($X=1$) of the Boltzmann equation, Eq.
(54), is

$\displaystyle\frac{\partial{}}{\partial{t}}~{}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}f_{\bm{p}}-\frac{\partial{}}{\partial{x^{i}}}~{}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{\partial{\varepsilon_{\bm{p}}^{*}}}{\partial{p_{i}}}~{}f_{\bm{p}}=0~{}.$
(66)

Using Eq. (59), the above equation becomes

$\displaystyle\frac{\partial{}}{\partial{t}}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}f_{\bm{p}}+\frac{\partial{}}{\partial{x^{i}}}g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{p^{i}-A^{i}}{\epsilon_{\textrm{kin}}^{*}}~{}f_{\bm{p}}=\partial_{t}j^{0}+\partial_{i}j^{i}=0~{},$
(67)

where we used the defining equations for the number current $j^{i}$ and
density $j^{0}$, Eqs. (44) and (45). The above result can be taken as a
confirmation that $j^{i}$ and $j^{0}$ have been correctly defined.

####  2 Multiple terms, kinetic momentum, and manifestly covariant equations
of motion

Having derived the properties of the VSDF model with one vector and one scalar
interaction term, we can easily extend the formalism to an arbitrary number of
interaction terms of each kind. Here, we are dealing with multiple vector
fields labeled by an index $k$ (where one should remember that $k$ is a
numbering subscript and not a space-time index),

$\displaystyle
A_{k}^{\lambda}(x;C_{k},b_{k})=C_{k}\big{(}j_{\nu}j^{\nu}\big{)}^{\frac{b_{k}}{2}-1}j^{\lambda}~{},$
(68)

in terms of which a general energy density with $K$ vector-density–dependent
terms and $M$ scalar-density–dependent terms can be written as

$\displaystyle\mathcal{E}_{(K,M)}$ $\displaystyle=$ $\displaystyle\int
d^{3}\tilde{p}~{}\epsilon^{*}_{\textrm{kin}}~{}f_{\bm{p}}+\sum_{k=1}^{K}A^{0}_{k}j_{0}-g^{00}\sum_{k=1}^{K}\left(\frac{b_{k}-1}{b_{k}}\right)A_{k}^{\mu}j_{\mu}+\sum_{m=1}^{M}G_{m}\left(\frac{d_{m}-1}{d_{m}}\right)n_{s}^{d_{m}}$
(69) $\displaystyle=$ $\displaystyle\int
d^{3}\tilde{p}~{}\varepsilon^{*}_{\bm{p}}~{}f_{\bm{p}}-g^{00}\sum_{k=1}^{K}\left(\frac{b_{k}-1}{b_{k}}\right)A_{k}^{\mu}j_{\mu}+\sum_{m=1}^{M}G_{m}\left(\frac{d_{m}-1}{d_{m}}\right)n_{s}^{d_{m}}~{}.$

We note that taking $K=2$, $M=0$ and evaluating $\mathcal{E}_{2,0}$ in the
rest frame results in the interaction of the same form as in commonly used
parametrizations of the Skyrme model (see, e.g., Ref. [178]); in particular,
using $b_{1}=2$ and $b_{2}=3$ ($b_{2}=\frac{13}{6}$) leads to the stiff (soft)
parametrization. In our approach, however, we leave the interaction
parameters, including the powers of vector and scalar number densities in the
interaction terms, unspecified until a later time, when we will fit them to
match chosen properties of nuclear matter.

The generalization of the remaining parts of the VSDF model is
straightforward, and in particular we arrive at the quasiparticle energy,

$\displaystyle\varepsilon_{\bm{p}}^{*}=\sqrt{\bigg{(}\bm{p}-\sum_{k=1}^{K}\bm{A}_{k}\bigg{)}^{2}+{m^{*}}^{2}}+\sum_{k=1}^{K}A^{0}_{k}~{},$
(70)

with the conserved current and the effective mass given by

$\displaystyle\bm{j}=\int
d^{3}\tilde{p}~{}\frac{\bm{p}-\sum_{k=1}^{K}\bm{A}_{k}}{\epsilon^{*}_{\textrm{kin}}}~{}f_{\bm{p}}$
(71)

and

$\displaystyle m^{*}=m_{0}-\sum_{m=1}^{M}G_{m}n_{s}^{d_{m}-1}~{},$ (72)

respectively, while the equations of motion become

$\displaystyle\frac{dx^{i}}{dt}$ $\displaystyle=$
$\displaystyle\frac{p^{i}-\sum_{k=1}^{K}(A_{k})^{i}}{\epsilon_{\textrm{kin}}^{*}}~{},$
(73) $\displaystyle\frac{dp^{i}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\bigg{(}p^{j}-\sum_{k=1}^{K}(A_{k})^{j}\bigg{)}\left(\sum_{k=1}^{K}\frac{\partial{(A_{k})_{j}}}{\partial{x_{i}}}\right)+m^{*}\frac{\partial{m^{*}}}{\partial{x_{i}}}}{\epsilon_{\textrm{kin}}^{*}}+\sum_{k=1}^{K}\frac{\partial{A_{k}^{0}}}{\partial{x_{i}}}~{}.$
(74)

We stress that the generalization to $(K,M)$ interaction terms preserves the
conservation laws and the relativistic covariance of the $T^{\mu\nu}$ tensor.

Finally, the equations of motion, Eqs. (73) and (74), can be rewritten in a
manifestly covariant way. First, we rewrite Eq. (70) as

$\displaystyle\varepsilon_{\bm{p}}^{*}-\sum_{k=1}^{K}A^{0}_{k}=\sqrt{\bigg{(}\bm{p}-\sum_{k=1}^{K}\bm{A}_{k}\bigg{)}^{2}+{m^{*}}^{2}}~{}.$
(75)

It is then natural to define a quantity known as the kinetic momentum
$\Pi^{\mu}$ [146],

$\displaystyle\Pi^{\mu}\equiv p^{\mu}-\sum_{k=1}^{K}(A_{k})^{\mu}~{},$ (76)

which by construction satisfies

$\displaystyle\Pi^{0}=\sqrt{\bm{\Pi}^{2}+m^{2}}~{}.$ (77)

Note that in terms of the kinetic momentum $\Pi^{\mu}$, the baryon 4-current
can be naturally defined as

$\displaystyle
j^{\mu}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{\Pi^{\mu}}{\Pi_{0}}~{}f_{\bm{p}}~{}.$
(78)

Using the kinetic momentum, one can rewrite the equations of motion as (see
Appendix 15.C for details)

$\displaystyle\frac{dx^{\mu}}{dt}=\frac{\Pi^{\mu}}{\Pi_{0}}~{},$ (79)
$\displaystyle\frac{d\Pi^{\mu}}{dt}=\sum_{\nu}\frac{\Pi_{\nu}}{\Pi_{0}}\sum_{k=1}^{K}\Big{[}\partial^{\mu}(A_{k})^{\nu}-\partial^{\nu}(A_{k})^{\mu}\Big{]}+\frac{m^{*}}{\Pi_{0}}\frac{\partial{m^{*}}}{\partial{x_{i}}}~{}.$
(80)

It is apparent that the first force term on the right-hand side of Eq. (80),
stemming from the vector interaction terms, is similar to the electromagnetic
force known from covariantly formulated electrodynamics [179], except that in
our case there is an arbitrary number $K$ of vector fields. In particular,
this means that the equation of motion in Eq. (80) naturally includes a
Lorentz force (see Appendix 15.D for an explicit demonstration).

####  3 Thermodynamics and thermodynamic consistency

The thermodynamic properties of the VSDF model are obtained in the following
way. In equilibrium, the distribution function must have the form of the
Fermi-Dirac distribution (for details, see Appendix 13),

$\displaystyle
f_{\bm{p}}=\frac{1}{e^{\beta(\varepsilon_{\bm{p}}^{*}-\mu)}+1}~{},$ (81)

where $\beta$ is the inverse temperature, $\beta=1/T$, and $\mu$ is the
chemical potential. We consider the energy-momentum tensor in the rest frame,
where it has the form

$\displaystyle T^{\mu\nu}=\textrm{diag}\big{(}\mathcal{E},P,P,P\big{)}$ (82)

and where the spatial components of the current vanish, $j^{i}=0$, while
$j_{\mu}j^{\mu}=n^{2}$, with $n$ being the rest frame density. Then the
pressure is given by

$\displaystyle P_{(K,M)}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\sum_{k}T^{kk}\bigg{|}_{\begin{subarray}{c}\text{rest}\\\
\text{frame}\end{subarray}}=$ (83) $\displaystyle=$ $\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}T~{}\ln\bigg{[}1+e^{-\beta\big{(}\varepsilon_{\bm{p}}^{*}-\mu_{B}\big{)}}\bigg{]}_{\begin{subarray}{c}\text{rest}\\\
\text{frame}\end{subarray}}$ $\displaystyle\hskip
14.22636pt+~{}\sum_{k=1}^{K}C_{k}\left(\frac{b_{k}-1}{b_{k}}\right)n^{b_{k}}-\sum_{m=1}^{M}G_{m}\left(\frac{d_{m}-1}{d_{m}}\right)n_{s}^{d_{m}}~{}.$

We note that in an equilibrated system, vector-density–dependent interactions
can be described in terms of a shift of the chemical potential $\mu$. Using
the expression for the quasiparticle energy, Eq. (70), in the exponent of the
distribution function, Eq. (81), we see that it is convenient to define an
effective chemical potential $\mu^{*}$ through

$\displaystyle\mu^{*}\equiv\mu-\sum_{k=1}^{K}A^{0}_{k}~{},$ (84)

so that

$\displaystyle\varepsilon_{\bm{p}}^{*}-\mu=\epsilon_{\textrm{kin}}^{*}-\mu^{*}=\Pi^{0}-\mu^{*}~{}.$
(85)

Consequently, the dependence of the thermal part of the pressure, Eq. (83), on
temperature $T$ and effective chemical potential $\mu^{*}$ is just like that
of an ideal Fermi gas, with the effects due to scalar interactions “folded
into” the effective mass $m^{*}$ and the effects due to vector interactions
“folded into” the effective chemical potential $\mu^{*}$.

The grand canonical potential is related to the pressure through
$\Omega(T,\mu,V)=-PV$, and we can immediately calculate the entropy density,

$\displaystyle s$ $\displaystyle\equiv$
$\displaystyle-\frac{1}{V}\left(\frac{d\Omega}{dT}\right)_{V,\mu}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\bigg{(}\ln\Big{[}1+e^{-\beta(\varepsilon_{\bm{p}}^{*}-\mu)}\Big{]}+\frac{\varepsilon_{\bm{p}}^{*}-\mu}{T}~{}f_{\bm{p}}\bigg{)}~{},$
(86)

and the number density,

$\displaystyle
n\equiv-\frac{1}{V}\left(\frac{d\Omega}{d\mu}\right)_{V,T}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}f_{\bm{p}}~{},$
(87)

where the latter equation proves the correct normalization of the distribution
function $f_{\bm{p}}$. We note that the entropy density, Eq. (86), can be
rewritten as

$\displaystyle
s=\frac{1}{T}\Big{[}\big{(}P_{FG}^{*}+\mathcal{E}_{FG}^{*}\big{)}-\mu^{*}n\Big{]}=s_{FG}^{*}~{},$
(88)

where we marked the dependence of the ideal Fermi gas pressure, energy
density, and entropy density on $m^{*}$ and $\mu^{*}$ by a star superscript.
Thus we see that the entropy density in the VSDF model is just like that of an
ideal Fermi gas with the effective mass $m^{*}$ and effective chemical
potential $\mu^{*}$.

Minimizing the thermodynamic potential (which in the grand canonical ensemble
is equivalent to finding the extremum of the pressure) with respect to the
auxiliary fields, $n$ and $n_{s}$, confirms the defining equations for the
effective mass and the baryon density, Eqs. (72) and (78) (see Appendix 15.E
for more details). Finally, calculating the energy density using the
fundamental equation of thermodynamics, $\mathcal{E}\equiv sT-P+\mu n$, yields
Eq. (69) evaluated in the rest frame, thus proving that the VSDF model is
thermodynamically consistent.

### 4\. Summary

In this chapter, we have given a short review of the mean-field approach to
describing nuclear matter and of the Landau Fermi-liquid theory, followed by a
derivation of the flexible vector and scalar density functional (VSDF) model.
The VSDF model, which obeys Lorentz covariance, preserves conservation laws,
and is shown to be thermodynamically consistent, allows one to construct a
parametrized dense nuclear matter EOS. In particular, the energy density and
pressure as obtained in the VSDF model can be used to fix the thermodynamic
properties of the system, while the corresponding quasiparticle energies and
quasiparticle equations of motion can be implemented in a hadronic transport
simulation.

For convenience, we recall the most important equations from the derivation
above. Starting from the energy density of the form

$\displaystyle\mathcal{E}_{(K,M)}$ $\displaystyle=$ $\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\epsilon^{*}_{\textrm{kin}}~{}f_{\bm{p}}+\sum_{k=1}^{K}A^{0}_{k}j_{0}-g^{00}\sum_{k=1}^{K}\left(\frac{b_{k}-1}{b_{k}}\right)A_{k}^{\mu}j_{\mu}$
(89) $\displaystyle\hskip
28.45274pt+~{}\sum_{m=1}^{M}G_{m}\left(\frac{d_{m}-1}{d_{m}}\right)n_{s}^{d_{m}}~{},$

where the vector $A_{k}^{\lambda}$ field is defined as

$\displaystyle
A_{k}^{\lambda}(x;C_{k},b_{k})=C_{k}\big{(}j_{\nu}j^{\nu}\big{)}^{\frac{b_{k}}{2}-1}j^{\lambda}~{},$
(90)

the kinetic energy is given by

$\displaystyle\epsilon_{\textrm{kin}}^{*}=\sqrt{\bigg{(}\bm{p}-\sum_{k=1}^{K}\bm{A}_{k}\bigg{)}^{2}+{m^{*}}^{2}}~{},$
(91)

and the conserved current and the scalar density are given by

$\displaystyle
j^{\mu}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{p^{\mu}-\sum_{k=1}^{K}A_{k}^{\mu}}{\epsilon^{*}_{\textrm{kin}}}~{}f_{\bm{p}}$
(92)

and

$\displaystyle
n_{s}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\frac{m^{*}}{\epsilon_{\textrm{kin}}^{*}}~{}f_{\bm{p}}~{}.$
(93)

respectively, with the effective mass satisfying

$\displaystyle m^{*}=m_{0}-\sum_{m=1}^{M}G_{m}n_{s}^{d_{m}-1}~{},$ (94)

we arrive at the quasiparticle energy

$\displaystyle\varepsilon_{\bm{p}}^{*}=\sqrt{\bigg{(}\bm{p}-\sum_{k=1}^{K}\bm{A}_{k}\bigg{)}^{2}+{m^{*}}^{2}}+\sum_{k=1}^{K}A^{0}_{k}~{},$
(95)

the quasiparticle equations of motion,

$\displaystyle\frac{dx^{i}}{dt}$ $\displaystyle=$
$\displaystyle\frac{p^{i}-\sum_{k=1}^{K}(A_{k})^{i}}{\epsilon_{\textrm{kin}}^{*}}~{},$
(96) $\displaystyle\frac{dp^{i}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\bigg{(}p^{j}-\sum_{k=1}^{K}(A_{k})^{j}\bigg{)}\left(\sum_{k=1}^{K}\frac{\partial{(A_{k})_{j}}}{\partial{x_{i}}}\right)+m^{*}\frac{\partial{m^{*}}}{\partial{x_{i}}}}{\epsilon_{\textrm{kin}}^{*}}+\sum_{k=1}^{K}\frac{\partial{A_{k}^{0}}}{\partial{x_{i}}}~{},$
(97)

and the pressure of the system,

$\displaystyle P_{(K,M)}$ $\displaystyle=$ $\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}T~{}\ln\bigg{[}1+e^{-\beta\big{(}\varepsilon_{\bm{p}}^{*}-\mu_{B}\big{)}}\bigg{]}_{\begin{subarray}{c}\text{rest}\\\
\text{frame}\end{subarray}}$ (98) $\displaystyle\hskip
14.22636pt+~{}\sum_{k=1}^{K}C_{k}\left(\frac{b_{k}-1}{b_{k}}\right)n^{b_{k}}-\sum_{m=1}^{M}G_{m}\left(\frac{d_{m}-1}{d_{m}}\right)n_{s}^{d_{m}}~{}.$

We again note that all of the presented expressions reduce to the Walecka
model for $K=M=1$ and $b_{1}=d_{1}=2$.

Let us once more stress that within the Landau Fermi-liquid theory, the form
of the energy density $\mathcal{E}$, Eq. (89), determines all other quantities
of interest. A natural question that arises here is: Does the particular
choice of $\mathcal{E}$ presented in Eq. (89) provide maximal flexibility? If
one considers a mean-field model with self-consistently defined baryon current
and effective mass, and one demands that the coefficients of the interaction
are constant, then it can be shown (see Appendix 16) that the interaction
terms must have the form of polynomials of vector and scalar number densities,
$n$ and $n_{s}$. While it is not entirely inconceivable to construct a model
within which the interaction coefficients are functions of $n$ and $n_{s}$,
such model would present significant computational challenges. Therefore for
the purposes of this thesis, which centers on applications in hadronic
transport models, an energy density of the form displayed in Eq. (89) provides
a maximal flexibility for a mean-field model of vector and scalar
interactions.

##  Chapter 2 Parametrization and model results

The contents of this chapter are largely based on Ref. [177]. We note,
however, that unlike in Ref. [177], we include antiparticles in our
description of nuclear matter.

The VSDF model (see Section 4 for the model summary) has been constructed with
applications to hadronic transport simulations in mind. Most importantly, in
opposition to theories in which the chemical potential $\mu$ is the
independent variable, we chose the number density $n$, naturally accessible in
hadronic transport, to be the dynamical variable of the theory (the scalar
number density is a function of $n$ and temperature $T$). The constructed EOS
is easily parametrizable through fitting the interaction parameters
$(C_{k},b_{k},G_{m},d_{m})$ to reproduce a chosen set of nuclear matter
properties. Here, the number of vector and scalar interaction terms, $K$ and
$M$ in Eq. (89), is entirely determined by the number of constraints that one
want to impose on the EOS. In particular, it is straightforward to construct a
family of EOSs that on the one hand reproduces the known properties of
ordinary nuclear matter, and on the other allows one to postulate and explore
critical behavior in dense nuclear matter over vast regions of the phase
diagram. The former will ensure that the model takes into the account the
behavior of nuclear matter known from the experiment (see Section 1 for more
details), while the latter will allow for meaningful comparisons of the
influence of different possible properties of the QGP phase transition on
observables. Such comparisons can be made, among others, through Bayesian
analysis [180, 181].

Apart from the flexibility of the EOS allowing for a systematic exploration of
possible phase diagrams, another factor to consider in simulations is
numerical efficiency. This is especially important for afterburner studies,
which often require simulating thousands or even millions of events, depending
on the studied observables, to reach statistical significance comparable with
the experiment. Here, vector-type interactions are more convenient for
numerical evaluation of mean-field potentials than scalar-type interactions,
which require solving a self-consistent equation at each point where mean-
fields are calculated (usually, this means at every node of a discrete
3-dimensional lattice representing the fields). Therefore, as an initial step
in studies of critical behavior in hadronic transport, in this thesis we
develop the parametrization of the flexible VSDF EOS that includes only
interactions of the vector type, that is we set $M=0$ in Eqs. (89-98). In this
case, the energy density reduces to a vector density functional (VDF) model of
the EOS. While using an energy density that is a VDF instead of a VSDF leads
to a smaller flexibility of the obtained family of EOSs (which is connected
with the fact that vector interactions depend only on the baryon density,
which in our model, along with the temperature $T$, is one of the two
independent variables), it nevertheless allows for creating a large family of
EOSs that can be easily utilized in simulations and whose behavior can be
easily understood. (We note that in spite of the significantly increased
numerical cost of simulations utilizing scalar mean-field interactions, this
cost is not entirely prohibitive, especially for observables requiring lower
statistics such as particle yields or elliptic flow, and future studies will
be devoted to this area; among the existing approaches, a prominent example of
simulations with scalar-type interactions are those using equations of motion
stemming from the Walecka model, see Refs. [182, 146, 183].)

To apply our approach to studies of heavy-ion collisions and the search for
the location of the conjectured QCD critical point, we want to parametrize the
VDF model so that it describes hadronic matter whose phase diagram contains
two first-order phase transitions. The first of these is the experimentally
observed low-temperature, low-density phase transition in nuclear matter,
sometimes known as the nuclear liquid-gas transition. The second is a
postulated high-temperature, high-density phase transition that is intended to
correspond to the QCD phase transition.

We want to stress that while the latter may, in principle, coincide with the
location of the phase transition in the real QCD phase diagram, its nature is
fundamentally different. This is because within Landau Fermi-liquid theory,
unlike in QCD, the degrees of freedom do not change across the phase
transition. Interestingly, this is also the case in some other approaches to
the QCD EOS, for example in models based on quarkyonic matter [184], where the
active degrees of freedom at the Fermi surface remain hadronic even after
quark degrees of freedom appear (although the extent to which such dynamics
may be captured in the VDF or VSDF model remains to be seen). The nature of
the phase transition that we can simulate in the VDF (or VSDF) model is that
of going from a less organized to a more organized state. This is easily
visualized in the case of the transition from gas to liquid (nucleon gas to
nuclear drop). In the case of the high-temperature, high-density phase
transition, we may think of it as a transition from a fluid to an even more
dense, and more organized, fluid (nuclear matter to quark matter). This
interpretation is supported by the functional dependence of entropy per
particle on the order parameter (number density), which decreases across the
phase transition from a less dense to a more dense state (for an extended
discussion, see Ref. [185]). For brevity, in the following we will refer to
the high-temperature, high-density phase transition within the VDF model as
the “QGP-like” or “quark-hadron” phase transition, with the expectation that
it is understood as a useful moniker rather than a statement on the nature of
the described transformation.

We also emphasize that the degrees of freedom used in the VDF model agree with
those expected after hadronization by design. Since the VDF model is
ultimately intended to be used in the hadronic afterburner stage of a hybrid
heavy-ion collision simulation, the issue of hadronic degrees of freedom
present above the QGP-like phase transition will never arise in realistic
calculations. At the same time, in parts of the phase diagram close to the
critical region, the hadronic systems studied will display effects typical for
approaching a phase transition, such as for example those expected due to the
softening of the EOS.

### 1\. Parametrization

We present a rather simplified version of the VDF model, in which we choose
the degrees of freedom to be those of isospin symmetric nuclear matter, that
is nucleons with nucleon mass $m_{N}=938$ MeV and degeneracy factor $g_{N}=4$.
Additionally, we also consider the case including thermally induced $\Delta$
resonances, whose mass and degeneracy factor are taken to be $m_{\Delta}=1232$
MeV and $g_{\Delta}=16$, respectively. We note that the model can be easily
extended to include arbitrarily many baryon resonances (which is done through
a substitution $gf_{\bm{p}}\to\sum_{i}g_{i}f^{(i)}_{\bm{p}}$), however, a
thorough study of the corresponding effects is beyond the scope of this
thesis. Finally, we modify the defining expressions of our model, Eqs.
(89-98), to include contributions from antiparticles. For example, in this
case $j^{0}$ represents the net baryon density given by

$\displaystyle j^{0}\equiv
j^{0}_{\textrm{baryons}}-j^{0}_{\textrm{antibaryons}}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\Big{[}f_{\bm{p}}-\bar{f}_{\bm{p}}\Big{]}~{},$
(1)

where we introduce the distribution function for antibaryons,

$\displaystyle\bar{f}_{\bm{p}}(T,\mu_{B})\equiv\frac{1}{e^{\beta\big{(}\varepsilon_{\bm{p}}+\mu_{B}\big{)}}+1}=f_{\bm{p}}(T,-\mu_{B})~{},$
(2)

and the remaining VDF model equations are likewise modified, in particular the
equations for the energy density

$\displaystyle\mathcal{E}_{(K,M)}$ $\displaystyle=$ $\displaystyle g\int
d^{3}\tilde{p}~{}\epsilon^{*}_{\textrm{kin}}~{}\Big{[}f_{\bm{p}}+\bar{f}_{\bm{p}}\Big{]}+\sum_{k=1}^{K}A^{0}_{k}j_{0}-g^{00}\sum_{k=1}^{K}\left(\frac{b_{k}-1}{b_{k}}\right)A_{k}^{\mu}j_{\mu}$
(3)

and the pressure,

$\displaystyle P_{(K,M)}$ $\displaystyle=$ $\displaystyle
g\int\frac{d^{3}p}{(2\pi)^{3}}~{}T~{}\ln\bigg{[}1+e^{-\beta\big{(}\varepsilon_{\bm{p}}^{*}-\mu_{B}\big{)}}\bigg{]}+g\int\frac{d^{3}p}{(2\pi)^{3}}~{}T~{}\ln\bigg{[}1+e^{-\beta\big{(}\varepsilon_{\bm{p}}^{*}+\mu_{B}\big{)}}\bigg{]}$
(4) $\displaystyle\hskip
14.22636pt+~{}\sum_{k=1}^{K}C_{k}\left(\frac{b_{k}-1}{b_{k}}\right)n_{B}^{b_{k}}~{},$

where $n_{B}\equiv\sqrt{j_{\mu}j^{\mu}}$ is the rest frame net baryon number
density.

In a system that undergoes two first-order phase transitions, the pressure
exhibits two mechanically unstable regions, known as the spinodal regions (see
Appendix 8 for more details), defined by the condition that the first
derivative of the pressure with respect to the order parameter is negative
[186, 187]. In a minimal VDF model realizing such behavior, the pressure needs
to be a four-term polynomial in the order parameter, and thus we adopt a
version of the VDF model in which we utilize four interaction terms. (We note
that to describe only one of the phase transitions mentioned above, it is
enough to adopt a model with two interaction terms. In the case of the nuclear
liquid-gas phase transition, the resulting model will be similar to many
Skyrme-based parametrizations of the EOS [176]). The energy density, Eq. (3),
is easily adapted to include $K=4$ interaction terms, and in the rest frame we
have

$\displaystyle\mathcal{E}\big{|}_{\begin{subarray}{c}\text{rest}\\\
\text{frame}\end{subarray}}=g\int\frac{d^{3}p}{(2\pi)^{3}}~{}\epsilon_{\textrm{kin}}~{}\Big{[}f_{\bm{p}}+\bar{f}_{\bm{p}}\Big{]}+\sum_{i=1}^{4}\frac{C_{i}}{b_{i}}n_{B}^{b_{i}}~{}.$
(5)

Similarly, we take $K=4$ in the expression for the pressure, Eq. (4).

To construct an EOS with a general QGP-like phase transition properties while
ensuring that the known properties of ordinary nuclear matter are well
reproduced, we choose the following constraints to fix the eight free
parameters $\\{b_{1},b_{2},b_{3},b_{4},C_{1},C_{2},C_{3},C_{4}\\}$ in the VDF
model:
1) the position of the minimum of energy per particle of nuclear matter at the
saturation density $n_{B}=n_{0}$,

$\displaystyle\frac{d\left(\frac{\mathcal{E}_{(K=4)}}{n_{B}}\right)}{dn_{B}}\bigg{|}_{\begin{subarray}{c}T=0\\\
n_{B}=n_{0}\end{subarray}}=0~{},$ (6)

2) the value of the binding energy at the minimum,

$\displaystyle\frac{\mathcal{E}_{(K=4)}}{n_{B}}\bigg{|}_{\begin{subarray}{c}T=0\\\
n_{B}=n_{0}\end{subarray}}-m_{N}=B_{0}~{},$ (7)

3, 4) the position of the critical point
$\big{(}T^{(N)}_{c},n^{(N)}_{c}\big{)}$ for the nuclear liquid-gas phase
transition,

$\displaystyle\frac{dP}{dn_{B}}\Big{(}T=T^{(N)}_{c},n_{B}=n^{(N)}_{c}\Big{)}=0~{},$
(8)
$\displaystyle\frac{d^{2}P}{dn_{B}^{2}}\Big{(}T=T^{(N)}_{c},n_{B}=n^{(N)}_{c}\Big{)}=0~{},$
(9)

5, 6) the position of the critical point
$\big{(}T^{(Q)}_{c},n^{(Q)}_{c}\big{)}$ for the quark-hadron phase transition,

$\displaystyle\frac{dP}{dn_{B}}\Big{(}T=T^{(Q)}_{c},n_{B}=n^{(Q)}_{c}\Big{)}=0~{},$
(10)
$\displaystyle\frac{d^{2}P}{dn_{B}^{2}}\Big{(}T=T^{(Q)}_{c},n_{B}=n^{(Q)}_{c}\Big{)}=0~{},$
(11)

7, 8) the position of the lower (left) and upper (right) boundaries of the
spinodal region, $\eta_{L}$ and $\eta_{R}$, for the quark-hadron phase
transition at $T=0$,

$\displaystyle\frac{dP}{dn_{B}}\Big{(}T=0,n_{B}=\eta_{L}\Big{)}=0~{},$ (12)
$\displaystyle\frac{dP}{dn_{B}}\Big{(}T=0,n_{B}=\eta_{R}\Big{)}=0~{}.$ (13)

(Explicit forms of the above equations are shown in Appendix 17.) The set of
quantities
$(n_{0},B_{0},T_{c}^{(N)},n_{c}^{(N)},T_{c}^{(Q)},n_{c}^{(Q)},\eta_{L},\eta_{R})$
is referred to as the characteristics of an EOS.

Table 1: Example characteristics $\big{(}T_{c}^{(Q)},n_{c}^{(Q)},\eta_{L},\eta_{R}\big{)}$ of the QGP-like phase transition: critical temperature $T_{c}^{(Q)}$, critical baryon number density $n_{c}^{(Q)}$, and the boundaries of the spinodal region at $T=0$, $\eta_{L}$ and $\eta_{R}$. The corresponding parameter sets can be found in Appendix 18. Characteristics in sets I-V are obtained based on systems composed only of nucleons, while in set VI we consider a system composed of nucleons and thermally produced $\Delta$-resonances. We also show the incompressibility at saturation density and zero temperature $K_{0}$ calculated for the parametrized EOSs. set |  $T_{c}^{(Q)}[\textrm{MeV}]$ |  $n_{c}^{(Q)}[n_{0}]$ |  $\eta_{L}[n_{0}]$ |  $\eta_{R}[n_{0}]\hskip 14.22636pt$ | particle species |  $K_{0}[\textrm{MeV}]$
---|---|---|---|---|---|---
I | 50 | 3.0 | 2.70 | 3.22 | N | 260
II | 50 | 3.0 | 2.85 | 3.12 | N | 279
III | 50 | 4.0 | 3.90 | 4.08 | N | 280
IV | 100 | 3.0 | 2.50 | 3.32 | N | 261
V | 125 | 4.0 | 3.60 | 4.28 | N | 278
VI | 125 | 4.0 | 3.60 | 4.28 | N + $\Delta$ | 277

We choose the properties of the ordinary nuclear matter, encoded in conditions
(6-9), based on experimentally determined values [12, 153]:

$\displaystyle n_{0}=0.160~{}\textrm{fm}^{-3}~{},\hskip
14.22636ptB_{0}=-16.3~{}\textrm{MeV}~{},$ (14) $\displaystyle
T^{(N)}_{c}=18~{}\textrm{MeV}~{},\hskip
14.22636ptn^{(N)}_{c}=0.06~{}\textrm{fm}^{-3}~{}.$ (15)

On the other hand, the properties of dense nuclear matter, $n_{B}\gg n_{0}$,
are only weakly constrained by experiment at this time. We are then in a
position to create a family of possible EOSs based on a number of different
postulated characteristics (10-13), while ensuring that nuclear matter
properties are preserved. The resulting family of EOSs encompasses QGP-like
phase transition characteristics spanning vast regions of the dense nuclear
matter phase diagram. This allows for a systematic comparison with
experimental data, with the goal of constraining the properties of allowed
EOSs to a small subfamily with qualitatively similar properties.

In the remainder of this chapter, we illustrate properties of the VDF model by
discussing key results for a few representative EOSs which exhibit chosen sets
of the QGP-like phase transition characteristics
$\big{(}T_{c}^{(Q)},n_{c}^{(Q)},\eta_{L},\eta_{R}\big{)}$ listed in Table 1.
The corresponding parameter sets can be found in Appendix 18.

### 2\. VDF model results: Pressure, the speed of sound, and energy per
particle

Figure 1: An example of the fitting procedure used for the VDF model. Pressure
is plotted as a function of baryon number density at three significant
temperatures ($T=0$, nuclear critical temperature $T_{c}^{(N)}$, and quark-
hadron critical temperature $T_{c}^{(Q)}$) for an EOS with characteristics
from set I, see Table 1 or the legend (where $T_{c}^{(Q)}$ is given in MeV,
while the critical density $n_{c}^{(Q)}$ and the boundaries of the spinodal
region at $T=0$, $\eta_{L}$ and $\eta_{R}$, are given in units of saturation
density, $n_{0}=0.160\ \textrm{fm}^{-3}$). Specific points at which the
parameters of the EOS are fixed are indicated on the plot; see text for more
details. Figure from [177].

We illustrate the fitting procedure in Fig. 1, where we show pressure as a
function of baryon number density at three significant temperatures ($T=0$,
nuclear critical temperature $T_{c}^{(N)}$, and quark-hadron critical
temperature $T_{c}^{(Q)}$) for an EOS with characteristics from set I (see
Table 1) and where we also indicate the location of key features that
determine the fit parameters. At temperature $T=0$, conditions (6) and (7) are
applied at the saturation density of nuclear matter, marked with a blue
circle; we note here that because at $T=0$ the pressure is given by

$\displaystyle P\equiv
n_{B}^{2}\frac{d}{dn_{B}}\left(\frac{\mathcal{E}}{n_{B}}\right)~{},$ (16)

condition (6) is equivalent to demanding that $P=0$. Also at $T=0$, conditions
(12) and (13) fix the positions of the lower (left) and upper (right) boundary
of the high density spinodal region, $\eta_{L}$ and $\eta_{R}$; these are
denoted with blue diamonds. At the critical point of nuclear matter,
$T=T_{c}^{(N)}$ and $n_{B}=n_{c}^{(N)}$, denoted with a green square,
conditions (8) and (9) are enforced. Finally, conditions (10) and (11) are
applied to set the position of the QGP-like critical point
$\big{(}T_{c}^{(Q)},n_{c}^{(Q)}\big{)}$, denoted with a red star. Note that in
general, the values of pressure at baryon densities corresponding to ordinary
nuclear matter are significantly smaller than the values at large baryon
densities, including the region where the QGP-like phase transition and the
corresponding softening of the EOS occur.

Fig. 2 shows pressure as a function of baryon number density at zero
temperature, where the curves correspond to all EOSs defined by the sets of
characteristics listed in Table 1. While most of the results are calculated in
the presence of nucleons only, the thin dotted red line shows pressure for a
system with both nucleons (protons and neutrons) and thermally excited
$\Delta$ resonances. By construction, all of the EOSs display the same
behavior for baryon number densities corresponding to ordinary nuclear matter
and only start differing from each other in regions currently not constrained
by experimental data, $n_{B}\gtrsim 1.5n_{0}$.

Figure 2: Pressure is plotted as a function of baryon number density at
temperature $T=0$ for all sets of characteristics listed in Table 1. In the
legend, the critical temperature of the QGP-like phase transition
$T_{c}^{(Q)}$ is given in MeV, while the critical density $n_{c}^{(Q)}$ and
the boundaries of the spinodal region at $T=0$, $\eta_{L}$ and $\eta_{R}$, are
given in units of saturation density, $n_{0}=0.160\ \textrm{fm}^{-3}$. The
hardness of the EOSs is noticeable for densities above the quark-hadron
transition regions, and is a consequence of employing interaction terms with
high powers ($b_{i}>2$) of baryon number density $n_{B}$ (see text for
details).

A few regularities are apparent in the behavior of the pressure curves at zero
temperature in regions corresponding to the QGP-like phase transition. Let us
focus on the value of the pressure at the lower boundary of the spinonal
region $P(\eta_{L})$ (which is directly related to the average value of the
pressure across the transition region), and compare its values for sets of
characteristics between which only one property of the QGP-like phase
transition changes substantially. First, $P(\eta_{L})$ increases with the
critical baryon number density $n_{c}^{(Q)}$, which can be seen by comparing
the pressure curves for the second and third sets of characteristics
(delineated with medium dashed green and thin dashed magenta lines,
respectively). Second, $P(\eta_{L})$ decreases with the critical temperature
$T_{c}^{(Q)}$, as evidenced by pressure curves for the first and fourth sets
of characteristics (delineated with thick dashed orange and solid purple
lines, respectively). Third, $P(\eta_{L})$ decreases with the width of the
spinodal region, $\Delta\eta=\eta_{R}-\eta_{L}$, which can be seen by
comparing pressure curves for the first and second sets of characteristics
(thick dashed orange and medium dashed green lines, respectively).
Furthermore, the magnitude of the drop in the pressure across the spinodal
region $\Delta P=P(\eta_{R})-P(\eta_{L})$, increases with the critical
temperature, as seen by comparing curves for the first and fourth sets of
characteristics (thick dashed orange and solid purple lines, respectively).
Importantly, these features create a physical bound on which QGP-like
transitions are allowed in the VDF model. A transition with a wide spinodal
region, with a critical point at a relatively low baryon number density and at
the same time a relatively high critical temperature can often be excluded, as
it leads to such a significant drop in the pressure across the spinodal region
that the pressure becomes negative in some parts of the quark-hadron
coexistence region, which would correspond to an unphysical “QGP bound state”.
This is because at $T=0$ the pressure is given by Eq. (16), and locally
negative pressure implies that there exists a baryon density for which
$d\left({\mathcal{E}}/{n_{B}}\right)/dn_{B}=0$ and
$d^{2}\left({\mathcal{E}}/{n_{B}}\right)/dn_{B}^{2}>0$, corresponding to a
local minimum in energy per particle, ${\mathcal{E}}/{n_{B}}$. While such a
minimum is in fact expected in the region of the phase diagram corresponding
to ordinary nuclear matter, where
$d\left({\mathcal{E}}/{n_{B}}\right)/dn_{B}=0$ at the nuclear saturation
density $n_{0}$, it is forbidden for large baryon number densities, where it
would correspond to a metastable or even stable state of QGP. As an example,
most obtained phase transitions with $n_{c}^{(Q)}=2.5n_{0}$ and
$T_{c}^{(Q)}\geq 125\ \textrm{MeV}$ are rejected based on this argument.

Figure 3: The isentropic speed of sound $c_{\sigma}^{2}$ at $T=0$ as a
function of baryon number density, plotted for all sets of characteristics
listed in Table 1. In the legend, the critical temperature of the QGP-like
phase transition $T_{c}^{(Q)}$ is given in MeV, while the critical density
$n_{c}^{(Q)}$ and the boundaries of the spinodal region at $T=0$, $\eta_{L}$
and $\eta_{R}$, are given in units of saturation density, $n_{0}=0.160\
\textrm{fm}^{-3}$. The speed of sound becomes acausal for relatively large
baryon number densities above the quark-hadron transition region, which is a
consequence of the hardness of the equation of state in the same region (see
Fig. 2). See text for details.

Next, it is easy to notice that the pressure rises rapidly after leaving the
quark-hadron transition region. This hardness of the EOS is a general feature
of models based on high powers of baryon number density (specifically, with
exponents higher than 2), and is ubiquitous among various Skyrme-type models
(see, e.g., Ref. [188]). In fact, it can be shown that any relativistic
Lagrangian with vector-type interactions leading, in the mean-field
approximation, to terms of the form $n_{B}^{\alpha}$, where $\alpha>2$,
results in acausal phenomena at high baryon number densities [189]. Indeed,
Fig. 3 shows the isentropic speed of sound squared
$\left(c_{\sigma}/c\right)^{2}$ at $T=0$ for the chosen sets of phase
transition characteristics listed in Table 1 (we note that at $T=0$, the
isothermal and isentropic speeds of sound are identical; for the derivation of
the formula used to calculate the speed of sound, see Appendix 20). The speed
of sound squared is negative within the spinodal region, as expected for a
first-order phase transition [187], while for large baryon number densities
above the quark-hadron phase transition it eventually becomes acausal.
Although the latter feature of the speed of sound is not ideal, such
pathological behavior of the EOS can be, in fact, expected outside of the
region in which its parameters are fitted. Moreover, it is not going to pose
challenges to using the VDF model in a hadronic afterburner, where nuclear
matter is simulated at densities below the coexistence region which in general
are not affected by this problem. Out of the abundance of caution, in creating
parameter sets we make sure that the speed of sound preserves causality for
all baryon number densities below the upper boundary of the quark-hadron
coexistence region. We note that in some of the studied phase transitions this
still allows for a speed of sound above the conformal bound of
$\left({c_{\sigma}}/{c}\right)^{2}\leq{1}/{3}$; it is presently unclear
whether this bound is satisfied in dense nuclear matter, see, e.g., Refs. [15,
190, 184, 191].

Figure 4: The binding energy at $T=0$ as a function of baryon number density,
plotted for all sets of characteristics listed in Table 1. In the legend, the
critical temperature of the QGP-like phase transition $T_{c}^{(Q)}$ is given
in MeV, while the critical density $n_{c}^{(Q)}$ and the boundaries of the
spinodal region at $T=0$, $\eta_{L}$ and $\eta_{R}$, are given in units of
saturation density, $n_{0}=0.160\ \textrm{fm}^{-3}$. The degree of the
softening in energy per particle at high baryon number density is directly
related to the width of the spinodal region of a given EOS (see text for more
details).

Finally, in Fig. 4 we show the binding energy $\mathcal{E}_{(4)}/n_{B}-m_{N}$
at $T=0$ as a function of baryon number density for EOSs corresponding to all
sets of characteristics listed in Table 1. All obtained EOSs describe the same
physics in the region $n_{B}\lesssim 1.5n_{0}$, where the behavior of nuclear
matter is relatively well known; in particular, all curves reproduce the value
of the chosen binding energy at nuclear matter saturation as well as the
location of the saturation density, Eq. (14). On the other hand, at high
densities the binding energy displays a softening related to the postulated
QGP-like phase transition, which is different for each of the considered EOS.
We note that the extent of this softening is directly related to the width of
the spinodal region of a given EOS. This can again be seen from the fact that
at zero temperature the pressure is given by Eq. (16), from which it
immediately follows that the curvature of the energy density,
$\left({d^{2}\mathcal{E}}/{dn_{B}^{2}}\right)$, must be negative in the
spinodal region; consequently, the region over which
$\left({d^{2}\mathcal{E}}/{dn_{B}^{2}}\right)<0$ is related to
$(\eta_{L},\eta_{R})$.

Figure 5: Pressure as a function of baryon number density at $T=0$ for 283 VDF
EOSs. Green, magenta, and blue curves mark EOSs with the critical density of
$n_{c}=3.0n_{0}$, $4.0n_{0}$, and $5.0n_{0}$, respectively, while light,
medium dark, and dark curves mark EOSs with the critical temperature
$T_{c}=50$, $100$, and $150\ \textrm{MeV}$. For each $(n_{c},T_{c})$, the left
boundary of the spinodal region at $T=0$ was set a series of values given by
$\eta_{L}=n_{c}-0.1-i\times 0.05n_{0}$, with $i=\\{0,1,2,\dots\\}$. Small
variations in the right boundary of the spinodal region at zero temperature
$\eta_{R}$, explained further in Section 3, result in “bundles” of EOSs that
correspond to slightly differing $\eta_{R}$’s for a given set of the EOS
characteristics $(n_{c},T_{c},\eta_{L})$.

Although we have only shown results corresponding to a few possible QGP-like
phase transitions, it should be understood that arbitrarily many versions of
the dense nuclear matter EOS can be obtained in the VDF model, as can be seen
in Fig. (5), which shows pressure curves at $T=0$ for all EOSs obtained by
demanding the critical density of $n_{c}^{(Q)}\in\\{3.0,4.0,5.0\\}\ [n_{0}]$
(green, magenta, and blue curves, respectively), the critical temperature of
$T_{c}^{(Q)}\in\\{50,100,150\\}\ [\textrm{MeV}]$ (light, medium, and dark
curves, respectively), and the lower boundary of the spinodal region varied
according to $\eta_{L}=n_{c}-0.1-i\times 0.05n_{0}$, where
$i=\\{0,1,2,\dots\\}$. While the obtained EOSs vary widely in the high baryon
density region, by construction they all reproduce the same physics in the
range of baryon number densities corresponding to ordinary nuclear matter,
$n_{B}\lesssim 1.5n_{0}$. In fact, fitting the VDF model to reproduce the
experimental values of the saturation density, the binding energy, and the
nuclear critical point gives a remarkably good prediction for the value of
pressure at the nuclear critical point $P_{c}$ and the value of
incompressibility at the saturation point $K_{0}$ as compared with experiment
and against other models (summarized in Table 2). This is partially expected,
as the value of the incompressibility $K_{0}$ depends strongly on the critical
temperature [192], and the latter is used in the VDF model as an input
parameter fixed based on the known experimental results. Nevertheless, it is
noteworthy that the minimal VDF model, based on a few characteristics taken at
their experimentally established values (here $n_{0}$, $B_{0}$, $T_{c}^{(N)}$,
$n_{c}^{(N)}$), leads to predictions for other properties of nuclear matter
agreeing remarkably well with experimental data. Apparently, constraining four
properties of the EOS is enough to reproduce the thermodynamic behavior of
nuclear matter in the fitted region. The same could be true in the case of
nuclear matter at high baryon number density. We may be hopeful that
postulating QGP-like phase transition characteristics that happen to lay close
to their true QCD values will lead to a VDF model parametrization correctly
describing other properties of dense nuclear matter in the transition region.
We expect that this correct description would manifest itself through
agreement of simulation results with experimental data.

Table 2: Comparison of the nuclear phase transition critical temperature $T_{c}^{(N)}\ [\textrm{MeV}]$, the critical baryon number density $n_{c}^{(N)}\ [\textrm{fm}^{-3}]$, pressure at the critical point $P_{c}\ [\textrm{MeV fm}^{-3}]$, and incompressibility $K_{0}\ [\textrm{MeV}]$ obtained in experiment [12] and in various models: the Walecka model [154, 155], the quantum Van der Waals model [193] (denoted with “Quantum VdW”), the VDF model with nuclear phase transition only (using two interaction terms and denoted by “V2DF (N)”), and the VDF model with both nuclear and quark-hadron phase transitions (using four interaction terms and denoted by “V4DF (N+Q)”). For the latter, values of $P_{c}$ and $K_{0}$ are given as averages calculated across all obtained EOSs for quark-hadron critical temperatures $T_{c}^{(Q)}\in\\{50,100,150\\}\ [\textrm{MeV}]$ and critical baryon number densities $n_{c}^{(Q)}\in\\{3.0,4.0,5.0\\}\ [n_{0}]$, see Fig. 5. Asterisks mark input values of the models. |  Experiment |  Walecka | Quantum VdW |  V2DF (N) |  V4DF (N+Q)
---|---|---|---|---|---
$T_{c}^{(N)}$ | $17.9\pm 0.4$ | 18.9 | 19.7 | 18* | 18*
$n_{c}^{(N)}$ | $0.06\pm 0.01$ | 0.070 | 0.072 | 0.06* | 0.06*
$P_{c}$ | $0.31\pm 0.07$ | 0.48 | 0.52 | 0.311 | $\langle P_{c}\rangle=0.3067\pm 0.0014$
$K_{0}$ | 230-315 | 553 | 763 | 282 | $\langle K_{0}\rangle=273.1\pm 5.1$

### 3\. VDF model results: Phase diagrams

The phase diagrams for the EOSs corresponding to the characteristics listed in
Table 1 are shown in Figs. 6 and 7. Solid and dashed lines represent the
boundaries of the coexistence and spinodal regions, respectively. The
coexistence and spinodal regions of the nuclear phase transition, depicted
with black lines, are common for all used EOSs by construction. For sets V and
VI, which differ only by the number of species contributing to the baryon
number density $n_{B}$, the phase diagrams in Fig. 6 coincide, which is
expected.

Figure 6: Phase diagrams in the $(T,n_{B})$ plane for sets of characteristics
listed in Table 1. Solid (dashed) lines represent the boundaries of the
coexistence (spinodal) regions. In the legend, the critical temperature of the
QGP-like phase transition $T_{c}^{(Q)}$ is given in MeV, while the critical
baryon number density $n_{c}^{(Q)}$ and the boundaries of the spinodal region,
$\eta_{L}$ and $\eta_{R}$, are given in units of saturation density,
$n_{0}=0.160\ \textrm{fm}^{-3}$. The coexistence and spinodal regions of the
nuclear phase transition, depicted with black lines, are common to all shown
EOSs. The phase diagrams of sets V and VI coincide, which is expected.

It is immediately apparent that the QGP-like coexistence curves in the phase
diagrams all look alike. This is a consequence of our choice to employ only
interactions depending on vector baryon number density, as in this case the
dependence of the thermal part of the pressure on temperature $T$ and
effective chemical potential $\mu^{*}$ is just like that of an ideal Fermi
gas, as can be seen from Eq. (98). Consequently, all VDF EOSs display similar
behavior with increasing temperature $T$. This can be especially easily seen
in the $(T,\mu_{B})$ phase diagram (Fig. 7), where the coexistence lines
exhibit the exact same curvature. An exception from this behavior as shown on
this phase diagram is the curve calculated for a system with both nucleons and
thermally produced $\Delta$ resonances (denoted with a red line), which bends
more forcefully towards the $\mu_{B}=0$ axis as the temperature increases.
This is to be expected, as including an additional baryon species lowers the
value of the baryon chemical potential for a given baryon number density.
Including even more baryon species would strengthen this effect.

Figure 7: Phase diagrams in the $(T,\mu_{B})$ plane for sets of
characteristics listed in Table 1. The legend is the same as in Fig. 6. Solid
lines represent the coexistence lines. The coexistence line of the nuclear
phase transition, depicted with a solid black line, is common to all sets of
characteristics. Also shown are chemical freeze-out points obtained in
experiment and a parametrization of the freeze-out line from [51]. The
noticeably different curvature of the coexistence line at high temperatures
for set IV (solid red line) is due to employing both nucleons and thermally
induced Delta resonances.

Another feature, easily discerned in the $(T,n_{B})$ phase diagram (Fig. 6),
is that the spinodal regions $[\eta_{L},\eta_{R}]$ (and likewise the
coexistence regions $[n_{L},n_{R}]$) are always approximately centered around
the critical baryon number density $n_{c}^{(Q)}$. This is again a consequence
of the fact that vector-like interactions do not have a temperature
dependence, so that the thermal contribution to the pressure is just like that
of an ideal gas of fermions with mass $m_{N}$ (for a detailed explanation see
Appendix 19). As a result, the critical baryon number density $n_{c}^{(Q)}$
and the boundaries of the spinodal region, $\eta_{L}$ and $\eta_{R}$, are not
independent. In consequence, we have effectively one less free parameter. For
example, once we set the ordinary nuclear matter properties, the critical
point of the quark-hadron phase transition, and the lower spinodal boundary at
$T=0$, $\eta_{L}$, the right spinodal boundary at $T=0$, $\eta_{R}$, is
practically fixed. This could be already observed in Fig. 5, where the minimal
allowed changes in $\eta_{R}$, given $(n_{c},T_{c},\eta_{L})$, result in
easily discernible “bundles” of EOSs.

We expect that all these regularities in the behavior of the spinodal and
coexistence lines would not be as prominent if scalar-type interactions were
included, rendering the thermal part of the pressure non-trivial. In
particular, we expect that this would allow us to obtain coexistence regions
bending towards the $n_{B}=0$ axis in the $(T,n_{B})$ plane, which would
correspond to an even stronger tendency to bend towards the $\mu_{B}=0$ axis
in the $(T,\mu_{B})$ plane. This expectation is based on the fact that,
typically, scalar interactions result in a small effective mass, which in
addition decreases with temperature, and that in turn produces a relatively
larger thermal contribution to the pressure for a given $n_{B}$ and $T$. As a
result, such phase transitions would more significantly affect the region of
the phase diagram covered by the Beam Energy Scan program at RHIC.
Parametrizations of the VSDF model leading to such effects are planned for the
near future.

### 4\. VDF model results: Cumulants of baryon number

In analyses of heavy-ion collision experiments, considerable attention has
been paid to cumulants of the baryon number distribution. In the grand
canonical ensemble, the $j$th cumulant of the baryon number $\kappa_{j}$ can
be calculated from

$\displaystyle\kappa_{j}=VT^{j-1}\frac{d^{j}P}{d\mu_{B}^{j}}~{},$ (17)

and the first four cumulants in terms of derivatives with respect to the
baryon number density were already given in Eqs. (28-31), repeated here for
convenience,

$\displaystyle\kappa_{1}=Vn_{B}=N_{B}~{},$ (18)
$\displaystyle\kappa_{2}=\frac{VTn_{B}}{\left(\frac{dP}{dn_{B}}\right)_{T}}~{},$
(19)
$\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]~{},$
(20)
$\displaystyle\kappa_{4}=\frac{VT^{3}n_{B}}{\left(\frac{dP}{dn_{B}}\right)_{T}^{3}}\left[1-\frac{4n_{B}}{\left(\frac{dP}{dn_{B}}\right)_{T}}\left(\frac{d^{2}P}{dn_{B}^{2}}\right)_{T}+\frac{3n_{B}^{2}}{\left(\frac{dP}{dn_{B}}\right)_{T}^{2}}\left(\frac{d^{2}P}{dn_{B}^{2}}\right)_{T}^{2}-\frac{n_{B}^{2}}{\left(\frac{dP}{dn_{B}}\right)_{T}}\left(\frac{d^{3}P}{dn_{B}^{3}}\right)_{T}\right]~{}.$
(21)

In the VDF model, analytic expressions for derivatives of the pressure
$d{{}^{j}P}/d{n_{B}^{j}}$ are available, and in particular we provide these
formulas up to $j=3$ in Appendix 21. Cumulants may also be expressed in terms
of derivatives of the baryon chemical potential
$d{{}^{j}\mu_{B}}/d{n_{B}^{j}}$; we provide these formulas up to the sixth-
order cumulant in Appendix 22. Which expressions are used is purely a question
of preference.

Figure 8: Contour plots of the cumulant ratio $\kappa_{2}/\kappa_{1}$ for the
fourth (IV) EOS listed in Table 1. Coexistence (spinodal) regions are denoted
with black (yellow) lines, and critical points are marked with yellow dots.
White regions correspond to $\kappa_{2}/\kappa_{1}=1\pm 0.03$. Legend entries
denote upper limits of values of the cumulant ratios.

While cumulants measure the derivatives of the EOS, Eq. (17), at the same time
they can be related to the moments of the baryon number distribution, see Eqs.
(23-26). Because cumulants of the proton number distribution (serving as a
proxy for the cumulants of the baryon number distribution) can be measured in
experiment, cumulants provide one of the strongest links between theoretical
predictions and experimental data. The experimental values of cumulants are
expected to be influenced by enhanced fluctuations of conserved charges in the
vicinity of the critical point, rendering them a signal for the existence of
the critical point and a first-order phase transition in QCD [94, 95, 96]. In
particular, it is argued that, for systems crossing the phase diagram close to
and above the critical point, the sign of the third-order cumulant
$\kappa_{3}$ will change once [97], while the sign of the fourth-order
cumulant $\kappa_{4}$ will change twice [98]. This is explained at length in
Section 6.

Figure 9: Contour plots of the cumulant ratio $\kappa_{3}/\kappa_{2}$ for the
fourth (IV) EOS listed in Table 1. Coexistence (spinodal) regions are denoted
with black (yellow) lines, and critical points are marked with yellow dots.
White regions correspond to $\kappa_{3}/\kappa_{2}=1\pm 0.03$. Legend entries
denote upper (lower) limits of positive (negative) values of the cumulant
ratios.

The explicit volume dependence of the cumulants, which is typically divided
out in theoretical calculations, is difficult to control in experiment.
Therefore, it is customary to consider ratios of cumulants, most commonly

$\displaystyle\frac{\sigma^{2}}{\mu}=\frac{\kappa_{2}}{\kappa_{1}}~{},\hskip
28.45274ptS\sigma=\frac{\kappa_{3}}{\kappa_{2}}~{},\hskip
28.45274pt\kappa\sigma^{2}=\frac{\kappa_{4}}{\kappa_{2}}~{},$ (22)

where $\mu$ denotes the mean, $\sigma^{2}$ denotes variance, $S$ denotes
skewness, and $\kappa$ denotes excess kurtosis. Preliminary results from the
Beam Energy Scan indeed suggest that the fourth-order cumulant ratio
${\kappa_{4}}/{\kappa_{2}}$ exhibits non-monotonic behavior with the collision
energy [104].

Figure 10: Contour plots of the cumulant ratio $\kappa_{4}/\kappa_{2}$ for the
fourth (IV) EOS listed in Table 1. Coexistence (spinodal) regions are denoted
with black (yellow) lines, and critical points are marked with yellow dots.
White regions correspond to $\kappa_{4}/\kappa_{2}=1\pm 0.03$. Legend entries
denote upper (lower) limits of positive (negative) values of the cumulant
ratios.

In this section, we will focus on results for the fourth (IV) set of
characteristics listed in Table 1. The choice of this set is arbitrary and
does not reflect any preference for the location of the QCD critical point,
but simply serves as an illustration of the properties of the VDF model which
are qualitatively comparable for all obtained EOSs. In Figs. 8, 9, and 10 we
plot the cumulant ratios (22) in the $(T,n_{B})$ and $(T,\mu_{B})$ planes; we
note here that gray areas on these diagrams signify regions in which either
the cumulant calculation is invalid (inside the spinodal regions, visible in
the $(T,n_{B})$ diagrams) or where data has not been produced (at very small
values of $n_{B}$, which affects the $(T,\mu_{B})$ diagrams). Dramatic
increase in magnitudes of cumulant ratios as well as sudden changes in sign,
observed in regions close to and above the critical point, agree with the
expectations for the behavior of the cumulants described above and in Section
6. Interestingly, the effects of the nuclear phase transition are clearly
present even at very high temperatures (this has been also observed in Ref.
[194]). This raises the question to what extent the presence of the nuclear
phase transition affects the interpretation of experimental data, either by
damping the signal originating at the QGP phase transition, or by acting as an
imposter. Such questions could be answered by comparing outcomes of
simulations utilizing a VDF EOS with either nuclear phase transition only, or
both nuclear and quark-hadron phase transitions. Studies of this type are
planned for future research and are beyond the scope of this thesis.

### 5\. Summary

In this chapter, we parametrized the vector density functional (VDF) model to
describe two first-order phase transitions: the experimentally observed
nuclear liquid-gas phase transition, and a postulated high-temperature, high-
density phase transition intended to model the QGP phase transition. We then
described gross properties of the obtained family of EOSs based on several
EOSs which are representative of the range of EOSs available in the
parametrized family. In particular, we focused on the behavior of the
pressure, the speed of sound, and the energy per particle as functions of the
baryon density, as well as on the features of the phase diagrams in the
$(T,n_{B})$ and $(T,\mu_{B})$ planes and the cumulants of baryon number.

##  Chapter 3 Implementation in a hadronic transport code

A significant part of work underlying this thesis was devoted to implementing
and testing the VDF equations of motion, Eqs. (79) and (80) with $M=0$, in the
hadronic transport code SMASH [38]. This implementation is available in SMASH
starting with the SMASH-2.1 release. As a hadronic transport code, SMASH
simulates hadronic non-equilibrium dynamics through numerically solving the
Boltzmann equation. Below, we briefly sketch the theoretical foundations of
hadronic transport simulations and discuss some of the technical details of
the implementation. A brief introduction to the kinetic theory is provided in
Appendix 23.

### 1\. Transport simulations

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

In statistical physics, a transport equation describes the time evolution of a
distribution function of a system of $N$ particles $f_{N}(\bm{x},\bm{p})$,
where we use a short-hand notation $\bm{x}=(\bm{x}_{1},\dots,\bm{x}_{N})$ and
$\bm{p}=(\bm{p}_{1},\dots,\bm{p}_{N})$. The distribution
$f_{N}(\bm{x},\bm{p})$ can be understood as the density of particles in the
phase space, or equivalently as the probability that the $i$-th particle has a
given position $\bm{x}_{i}$ and a given momentum $\bm{p}_{i}$, $i=1,\dots,N$.
From the Liouville theorem (see Appendix 23.C) we know that the time evolution
of $f_{N}$ satisfies

$\displaystyle\frac{df_{N}}{dt}=\left(\frac{\partial{f_{N}}}{\partial{t}}\right)_{\textrm{drift}}+\frac{d\bm{x}}{dt}\frac{\partial{f_{N}}}{\partial{\bm{x}}}+\frac{d\bm{p}}{dt}\frac{\partial{f_{N}}}{\partial{\bm{p}}}=0~{}.$
(1)

In practice, many of the most relevant features of the system are described by
the 1-body distribution function $f_{1}(\bm{x}_{1},\bm{p}_{1})$, which gives
the probability that any of the $N$ particles in the system has a given
position $\bm{x}_{1}$ and a given momentum $\bm{p}_{1}$. The time evolution of
$f_{1}(\bm{x}_{1},\bm{p}_{1})$ is given by

$\displaystyle\frac{\partial{f_{1}}}{\partial{t}}+\frac{d{\bm{x}_{1}}}{d{t}}\frac{\partial{f_{1}}}{\partial{\bm{x}_{1}}}+\frac{d{\bm{p}_{1}}}{d{t}}\frac{\partial{f_{1}}}{\partial{\bm{p}_{1}}}=I_{\textrm{coll}}~{},$
(2)

where on the right-hand side we introduced the collision integral
$I_{\textrm{coll}}$, which is an often-used common notation for all terms
contributing to changes in $f_{1}$ due to particle-particle collisions
(particles colliding into or out of the phase space volume element) and
particle transformations such as decays or resonance formation. If the
collision term is constructed according to the assumption of molecular chaos
(meaning that particles are not correlated neither before nor after a
collision takes place), also known as the Stosszahlansatz, then Eq. (2) is
known as the Boltzmann equation. A formal derivation of the Boltzmann equation
from the Liouville theorem, Eq. (1), can be obtained using the BBGKY hierarchy
of equations (see Appendix 23.E for a simplified derivation of the Boltzmann
equation from the BBGKY hierarchy, which is itself derived in Appendix 23.D).
In the following, we will suppress the subscript “1” and assume that
$f(\bm{x},\bm{p})$ always refers to a 1-body distribution function.

In the context of heavy-ion collision simulations, the Boltzmann equation is
often also called the Boltzmann-Uehling-Uhlenbeck (BUU) equation, the Vlasov-
Uehling-Uhlenbeck (VUU) equation, or the Boltzmann-Nordheim-Vlasov (BNV)
equation. The specification comes from considering the phase space
distribution in the presence of a self-consistent mean field
$U(\bm{x},\bm{p})$ entering the dynamics of the system through the single-
particle Hamiltonian $H_{(1)}=\sqrt{\bm{p}^{2}+m^{2}}+U(\bm{x},\bm{p})$,

$\displaystyle\bigg{[}\frac{\partial{}}{\partial{t}}+\frac{\partial{H_{(1)}}}{\partial{\bm{p}}}\bm{\nabla}_{\bm{x}}-\frac{\partial{H_{(1)}}}{\partial{\bm{x}}}\bm{\nabla}_{\bm{p}}\bigg{]}f(t,\bm{x},\bm{p})=I_{\textrm{coll}}~{},$
(3)

where we inserted Hamilton’s equations of motion into the Boltzmann equation
(2). In cases where there are no collisions, $I_{\textrm{coll}}=0$, the term
“Vlasov equation” is also often used.

The function $f(t,\bm{x},\bm{p})$ is a continuous distribution function for a
given total number of nucleons, denoted by $A$, and solving the Boltzmann
equation is equivalent to obtaining the time evolution of
$f(t,\bm{x},\bm{p})$. (We note here that in general, heavy-ion collision
simulations evolve hundreds of hadron species, described by a system of
hundreds of coupled Boltzmann equations; the following discussion, treating
explicitly the Boltzmann equation in the presence of nucleons only, can be
straightforwardly generalized to include all necessary particle species.)
Numerically, given the initial condition in form of the distribution function
at some time $t_{0}$, $f(t_{0},\bm{x}_{0},\bm{p}_{0})$, we solve for the
distribution at a slightly later time $t=t_{0}+\delta t$, $f(t_{0}+\delta
t,\bm{x},\bm{p})$, and repeat the process until a final time
$t=t_{\textrm{end}}$ is reached. In more detail, the numerical solution of the
VUU equation is achieved within a numerical approach known as the method of
test particles [195], which is based on the assumption that the continuous
phase space distribution $f(t,\bm{x},\bm{p})$ of a system of $A$ nucleons can
be approximated by the distribution of a large number $N$ of discrete test
particles with phase space coordinates
$\big{(}\bm{x}_{i}(t),\bm{p}_{i}(t)\big{)}$,

$\displaystyle
f(t,\bm{x},\bm{t})\approx\frac{1}{N_{T}}\sum_{i=1}^{N}\delta\Big{(}\bm{x}-\bm{x}_{i}(t)\Big{)}\delta\Big{(}\bm{p}-\bm{p}_{i}(t)\Big{)}~{}.$
(4)

Here, $N_{T}$ is the number of test particles per nucleon and $N=N_{T}A$. Each
test particle contributes to the total charge in the system with a charge of
the corresponding real particle divided by $N_{T}$ (for example, a “nucleon
test particle” contributes a baryon charge of ${1}/{N_{T}}$), so that the
total charge in the simulation equals that of a system of $A$ particles. The
cornerstone of the method of test particles is the realization that if we
demand that these test particles are propagated according to

$\displaystyle\frac{d\bm{x}}{dt}=\frac{\partial{H_{(1)}}}{\partial{\bm{p}}}~{},\hskip
42.67912pt\frac{d\bm{p}}{dt}=-\frac{\partial{H_{(1)}}}{\partial{\bm{x}}}~{},$
(5)

then the Vlasov equation, which is the left-hand side of Eq. (3), immediately
follows from the Liouville theorem. In other words, given some configuration
of the positions and momenta of the test particles at a time $t_{0}$, which
approximates the phase space distribution $f(t_{0},\bm{x},\bm{p})$ according
to Eq. (4), we can propagate the test particles over a time $dt$ using the
Hamilton equations, Eq. (5), and use the configuration of positions and
momenta of the propagated test particles to approximate
$f(t_{0}+dt,\bm{x},\bm{p})$. We note here that in SMASH, the equations of
motion propagate the kinetic momentum of particles, see Eq. (76), which is
numerically more straightforward. An alternative approach, in which the
canonical momenta are propagated, is possible [196].

In practice, there exist two ways of realizing the method of test particles in
hadronic transport. Within the first approach, one initializes a system with
$N_{T}A$ test particles, which are then propagated according to the equations
of motion. Scatterings are performed according to cross sections that are
scaled as $\sigma/N_{T}$, where $\sigma$ is the physical cross section, which
ensures that the average number of scatterings per test particle is the same
as in a system of $A$ particles. Each test particle carries the charge of the
corresponding real particle, but contributes to the density with a scaling
factor of $1/N_{T}$, which preserves the total baryon number evolved in the
simulation and, for $N_{T}\gg 1$, results in a mean field that is a smoothed
out version of the mean field corresponding to $A$ particles. This approach is
sometimes referred to as the “full ensemble”.

An alternative approach is known as “parallel ensembles” [197]. In this
method, $N_{T}$ instances of a system of $A$ particles are created. Particles
in each instance are propagated according to the equations of motion, and
scatterings are performed using the physical cross section $\sigma$. Each test
particle carries the charge of the corresponding real particle, but
contributes to the density and the corresponding mean field, which are
calculated by summing contributions from all $N_{T}$ instances of the system,
with a scaling factor of $1/N_{T}$. Evolving the $N_{T}$ systems with mean
fields calculated in this fashion means that the systems are not in fact
independent, and their evolution due to the mean fields is shared. At the same
time, this approach is computationally much more efficient, as collision
searches are performed only within individual instances of the system, thus
reducing the numerical cost of performing the collisions within one event by a
factor on the order of $\mathcal{O}(N_{T}^{2})$, and the overall cost of the
simulation by a factor on the order of $\mathcal{O}(N_{T})$ (this is because
the parallel ensembles mode evolves $N_{T}$ events at the same time, while the
full ensemble mode evolves only one event).

It can be checked that these two simulation paradigms lead to the same results
in typical cases [198]. In this thesis we utilize the full ensemble approach
to the test particle method, as this is the primary method supported in SMASH.
Recently, the option to run the simulations in the parallel ensembles mode has
been added to SMASH [199], and it is currently being tested.

Let us discuss in some detail the numerical method used to evolve the test
particles in time. In SMASH, the positions and momenta of particles are
advanced using the leapfrog algorithm. In this method, the evolutions of
positions and momenta are “interlaced” in the following way: We start from
some initial positions $\bm{x}_{0}$ and momenta $\bm{p}_{0}$. The momenta at
time $t_{0}$ are advanced by half a time step to time $t_{0}+\frac{1}{2}dt$
using forces calculated based on positions at the initial time
$F(\bm{x}_{0})$, $\bm{p}_{0}\to\bm{p}_{1/2}$. Then, the positions of particles
are advanced through the time step $\Delta t$ using the values of the momenta
$\bm{p}_{1/2}$, $\bm{x}_{0}\to\bm{x}_{1}$. These positions are used to
calculate the forces on test particles $F(\bm{x}_{1})$, which are then used to
update the values of momenta from $\bm{p}_{1/2}$ to $\bm{p}_{3/2}$. Then
positions $\bm{x}_{1}$ are advanced through the second time step $\Delta t$
using the momenta $\bm{p}_{3/2}$, and so on. Thus within the leapfrog method,
the positions and momenta are updated at interleaved time points, so that the
evolutions of $\bm{x}$ and $\bm{p}$ “leapfrog” over each other. This is
markedly different from Euler integration, in which we would advance the
positions from $\bm{x}_{0}$ to $\bm{x}_{1}$ using momenta $\bm{p}_{0}$, and
advance the momenta from $\bm{p}_{0}$ to $\bm{p}_{1}$ using forces calculated
based on the particle positions $\bm{x}_{0}$ (see Fig. 1).

Initial state: $x(t_{0})$, $p(t_{0})$ Take $x_{0}$, $p_{0}$ Calculate the
force: $F(x_{0})$ Propagate by $\Delta t$: $p_{1}=p_{0}+F(x_{1})\Delta t$,
$x_{1}=x_{0}+p_{0}\Delta t$, $t_{1}=t_{0}+\Delta t$ Does
$t_{1}=t_{\textrm{end}}$? Substitute new initial state: $x_{0}=x_{1}$,
$p_{0}=p_{1}$, $t_{0}=t_{1}$ no Final state: $x_{1}$, $p_{1}$ yes Initial
state: $x(t_{0})$, $p(t_{0})$ Calculate the force: $F\big{(}x(t_{0})\big{)}$
Propagate $p(t_{0})$ by $\frac{1}{2}\Delta t$:
$p_{1/2}=p_{0}+F(x_{1})\frac{\Delta t}{2}$ Take $x_{0}$, $p_{1/2}$ Propagate
$x_{0}$ by $\Delta t$: $x_{1}=x_{0}+p_{1/2}\Delta t$, calculate the force:
$F(x_{1})$ Propagate $p_{1/2}$ by $\Delta t$: $p_{3/2}=p_{1/2}+F(x_{1})\Delta
t$, update time: $t_{1}=t_{0}+\Delta t$ Does $t_{1}=t_{\textrm{end}}$?
Substitute new initial state: $x_{0}=x_{1}$, $p_{1/2}=p_{3/2}$, $t_{0}=t_{1}$
no Final state: $x_{1}$, $p_{3/2}$ yes Figure 1: A schematic representation of
the Euler (top diagram) and leapfrog (bottom diagram) integration algorithms.

In fact, it can be shown that the leapfrog algorithm, although it requires the
same number of steps as the Euler integration (which is a first-order method),
is a second-order method, and consequently it leads to more stable solutions.
Indeed, we have tested the quality of solving the equations of motion using
the SMASH leapfrog algorithm against the 4th-order Runge-Kutta integrator by
tracking the energy conservation throughout the simulations (although energy
should be conserved exactly by the Boltzmann equation, various effects
connected with the method of integration, density calculation, and other
elements of the simulation result in deviations from this theoretical
expectation; studies of this nature are further discussed in Section 4). We
confirmed that the two algorithms behave comparably for the standard
simulation time step of $\Delta t=0.1\ \textrm{fm}/c$. At the same time, the
leapfrog algorithm is considerably more numerically efficient. We note here
for completeness that the leapfrog algorithm as implemented in SMASH misses
the first step of advancing the initial momenta by half a time step. Such an
omission could be disastrous for applications of the leapfrog integrator to
problems with precise initial conditions, however, in SMASH the initial
conditions are always created by random sampling of the desired phase space
distribution with a large number of test particles. Advancing the initial
momenta by half a time step is in this case equivalent to another random
sample of the phase space distribution, and in result has no tangible
consequences; this has been confirmed in tests.

Finally, we note that while the momentum changes due to the presence of
potentials are addressed in SMASH at evolution times separated by a constant
time interval $\Delta t$, particle-particle collisions and decays can take
place at any time during advancing the particle positions. As a result, in
SMASH the propagation within one time step is performed in the following way:

  1. 1.

At $t=t_{0}$: Take particles’ momenta $\bm{p}$; find all collisions (and
decays) that will occur based on particle trajectories given by
$\bm{x}(t_{0})+\bm{p}\Delta t$, and order them according to the time at which
they occur; identify the time of the first occurring collision (or decay)
$t^{\prime}$.

  2. 2.

Propagate the particles in straight lines until $t=t^{\prime}$.

  3. 3.

At $t=t^{\prime}$: perform the collision (decay), which in particular will
produce new momenta for the interacting particles; find collisions (decays)
for the particles that have just interacted, include them in the list of all
collisions (decays); identify the time of the next occurring collision (decay)
$t^{\prime\prime}$.

  4. 4.

Propagate the particles until $t=t^{\prime\prime}$, and repeat the above.
Repeat steps 2-3 until $t=t_{0}+\Delta t$ is reached.

  5. 5.

At $t=t_{0}+\Delta t$: Update particles’ momenta according to the equations of
motion, with forces calculated at the propagated positions,
$F\big{(}\bm{x}(t+\Delta t)\big{)}$.

Steps 1-5 are then repeated for each time step until the end of the
simulation.

In the next section, we discuss some details of the density and mean field
calculation in SMASH. There are, however, many other features and properties
of SMASH not presented in this brief overview. A complete SMASH User Guide can
be found online [200], and similarly for the SMASH development documentation
[201].

### 2\. Calculation of mean fields

Parts of this section are based on Ref. [177].

For an evolution without mean fields (that is for an evolution which does not
require density calculations) it is most natural to take the number of test
particles exactly corresponding to the actual number of nucleons present in
the system, $N=A$, $N_{T}=1$. However, employing mean-fields dependent on
local density and its gradients requires adopting an approach in which the
statistical noise due to a finite number of test particles is suppressed. This
is especially important in the case of models with competing repulsive and
attractive potentials of large magnitudes, where relatively small numerical
fluctuations can produce significant errors in the mean field calculation.
Thus, for example, for studies of nuclear matter with average density around
the saturation density $n_{0}$, a number of test particles per nucleon
$N_{T}=100$ is often used.

The local baryon current can be defined on a discrete lattice, and the current
at a given lattice node can be calculated as a sum of contributions from all
test particles which are in the volume element $V_{i}$ corresponding to that
lattice node,

$\displaystyle j^{\mu}(\bm{r}_{i})=\frac{1}{N_{T}}\frac{1}{V_{i}}\sum_{k\in
V_{i}}B_{k}\frac{\Pi^{\mu}_{k}}{\Pi^{0}_{k}}~{},$ (6)

where $B_{k}$ is the baryon charge of the $k$-th test particle. This
prescription, known as the particle-in-cell method, naturally reproduces the
baryon number in a given volume element,

$\displaystyle B(i)=j^{0}(i)V_{i}=\frac{N(i)}{N_{T}}~{},$ (7)

where $N(i)$ is the (net) number of test particles in $V_{i}$. In practice, in
order for the local densities and currents to be smooth enough, a prescription
is used in which currents at a given lattice node $i$ are weighted sums of
contributions from all test particles in some chosen volume $V_{s}$ around
that lattice node $i$, with $V_{s}$ taken to be larger than the volume element
$V_{i}$, $V_{s}>V_{i}$. In this case the baryon current is given by

$\displaystyle j^{\mu}(\bm{r}_{i})=\sum_{k\in
V_{s}}B_{k}\frac{\Pi^{\mu}_{k}}{\Pi^{0}_{k}}S(\bm{r}_{i}-\bm{r}_{k})~{},$ (8)

where the weight $S(\bm{r}_{i}-\bm{r}_{k})$ is known as the smearing function
(often also referred to as the smearing factor or the smearing kernel),
normalized such that

$\displaystyle\sum_{i}V_{i}~{}S(\bm{r}_{i}-\bm{r}_{k})=\frac{1}{N_{T}}~{}.$
(9)

Note that the volume $V_{s}$ from which test particles contribute to the
density at a given node effectively establishes the range of the interaction
and therefore it is a parameter of the algorithm that has a physical
relevance; while this smearing range can be adjusted to reproduce certain
physical behaviors of studied systems, in most applications it is often set at
around $r\approx 2\ \textrm{fm}$ [202]. Having calculated the current
$j^{\mu}$, current gradients are usually calculated by taking a finite
difference of baryon current values at adjacent lattice nodes, for example

$\displaystyle\frac{dj_{0}(x_{i},y_{j},z_{k})}{dx}=\frac{j_{0}(x_{i}+1,y_{j},z_{k})-j_{0}(x_{i}-1,y_{j},z_{k})}{2L_{x}}~{},$
(10)

where $L_{x}$ is the lattice spacing in the $x$-direction. Consequently,
derivatives calculated in hadronic transport can only resolve structures of
size bigger than the lattice spacing.

To achieve a higher resolution, one can in principle use a smaller lattice
spacing, however, the computational cost of using finer lattices is often
prohibitive. This is because finer lattices not only mean that there are more
lattice nodes at which the baryon current is computed, but more importantly,
the number of test particles per particle must be increased so that the number
of test particles per lattice cell is high enough to ensure a smooth density
calculation. This significantly affects the time cost of the simulation, and
the effect is especially pronounced for simulations run in the full ensemble
mode, where the number of performed particle-particle collisions scales as
$N_{T}^{2}$. As a result, one needs to find a compromise between the
resolution and a viable computation time. In practice, the lattice spacing in
SMASH is often chosen at a typical value of $1\ \textrm{fm}$ (which is also
the case in this work), and so the constraint on the resolution coming from
the range of the smearing factor is typically more important.

The time derivatives of densities and currents are most often also calculated
using the finite difference method, based on the difference between the values
of the baryon current at times $t_{0}+\delta t$ and $t_{0}$ at each node. The
resolution here naturally depends on the time step $\Delta t$ used to advance
the evolution.

Finally, because the positions of the test particles take continuous values,
while the lattice nodes are located at discrete positions $\bm{x}_{i}$, a
question arises what density should be used for the mean-field evolution of a
test particle located at a position $\bm{x}\neq\bm{x}_{i}~{}\forall i$.
Currently in SMASH this is treated within a zeroth-order approximation, where
the density at $\bm{x}$ is assumed to be equal to the density at the nearest
lattice node. Naturally, this approximation is good only as long as the
density doesn’t change appreciably over distances on the order of
$L_{i}/\sqrt{2}$.

The existing transport codes employ various numerical approaches for
calculating the effects of the mean-field potentials (for an exhaustive
comparison of mean-field calculations in nuclear matter simulated in several
different transport codes, see Ref. [202]). In the case of results presented
in Chapter 4, we employ a triangular smearing function, originating from the
lattice Hamiltonian method of solving nuclear dynamics [203], and derivatives
based on the finite difference method. Below, we briefly discuss and compare
three smearing schemes available in SMASH together with some additional
details on the calculation of derivatives. The Gaussian smearing scheme is the
original smearing available in SMASH, while the discrete and triangular
smearing have been added as part of this work. Likewise, we introduced the
option to calculate the gradients of currents using the chain rule or direct
finite difference derivatives scheme. All of the options discussed below are
available starting with the SMASH-2.1 release. We note here for clarity that
although we always refer to calculating the baryon current, in SMASH it is
also possible to calculate currents of different particle charges such as,
e.g., isospin.

####  1 Discrete smearing

Within a smearing scheme which in SMASH is referred to as “discrete smearing”,
each test particle contributes a fraction $F_{\textrm{c}}=1/a$ of its charge
to the lattice cell in which it is located (center cell) and a fraction
$F_{\textrm{n}}=(1-F_{\textrm{c}})/6=(a-1)/6a$ of its charge to each of its
six nearest neighbor lattice cells, that is nearest lattice cells in the $-x$,
$+x$, $-y$, $+y$, $-z$, and $+z$ directions. In order for
$F_{\textrm{c}}>F_{\textrm{n}}$ to hold, the discrete smearing parameter must
satisfy $F_{\textrm{c}}\geq 1/7$; an often used choice is