# Molecular dynamics analysis of particle number fluctuations in the mixed
phase of a first-order phase transition

Volodymyr A. Kuznietsov  Bogolyubov Institute for Theoretical Physics, 03680
Kyiv, Ukraine  Leibniz Institute For Solid State and Material Research, 01069
Dresden, Germany,  Frankfurt Institute for Advanced Studies, Giersch Science
Center, Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main, Germany  Oleh Savchuk
Facility for Rare Isotope Beams, Michigan State University, East Lansing, MI
48824 USA  Bogolyubov Institute for Theoretical Physics, 03680 Kyiv, Ukraine
GSI Helmholtzzentrum für Schwerionenforschung GmbH, Planckstr. 1, D-64291
Darmstadt, Germany  Frankfurt Institute for Advanced Studies, Giersch Science
Center, Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main, Germany  Roman V.
Poberezhnyuk  Bogolyubov Institute for Theoretical Physics, 03680 Kyiv,
Ukraine  Frankfurt Institute for Advanced Studies, Giersch Science Center,
Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main, Germany  Volodymyr Vovchenko
Physics Department, University of Houston, Box 351550, Houston, TX 77204, USA
Nuclear Science Division, Lawrence Berkeley National Laboratory, 1 Cyclotron
Road, Berkeley, California 94720, USA  Mark I. Gorenstein  Bogolyubov
Institute for Theoretical Physics, 03680 Kyiv, Ukraine  Frankfurt Institute
for Advanced Studies, Giersch Science Center, Ruth-Moufang-Str. 1, D-60438
Frankfurt am Main, Germany  Horst Stoecker  Frankfurt Institute for Advanced
Studies, Giersch Science Center, Ruth-Moufang-Str. 1, D-60438 Frankfurt am
Main, Germany  GSI Helmholtzzentrum für Schwerionenforschung GmbH, Planckstr.
1, D-64291 Darmstadt, Germany  Institute of Theoretical Physics, Goethe
Universität, Frankfurt, Germany

(August 28, 2024)

###### Abstract

Molecular dynamics simulations are performed for a finite non-relativistic
system of particles with Lennard-Jones potential. We study the effect of
liquid-gas mixed phase on particle number fluctuations in coordinate subspace.
A metastable region of the mixed phase, the so-called nucleation region, is
analyzed in terms of a non-interacting cluster model. Large fluctuations due
to spinodal decomposition are observed. They arise due to the interplay
between the size of the acceptance region and that of the liquid phase. These
effects are studied with a simple geometric model. The model results for the
scaled variance of particle number distribution are compared with those
obtained from the direct molecular dynamic simulations.

mixed phase, fluctuations, molecular dynamics

##  I Introduction

The endpoint of a first-order phase transition, noted as the critical point
(CP), occurs under different physical conditions, including most molecular and
ferromagnetic systems Landau and Lifshitz (1975); Greiner  _et al._ (2012),
nuclear matter Allen and Tildesley (2017), and potentially the hot QCD matter
at nonzero baryon density Stephanov  _et al._ (1999); Bzdak  _et al._ (2020).
In the thermodynamic limit, particle number fluctuations exhibit singular
behavior at the CP. These singularities are smeared out in finite-size
systems. Nevertheless, small systems also demonstrate specific features of
critical behavior such as enhancement of fluctuations Bernhardt  _et al._
(2022); Kuznietsov  _et al._ (2022).

Event-by-event fluctuations in nucleus-nucleus collisions are used as an
experimental tool to search for the QCD CP at finite baryon density Stephanov
_et al._ (1999); Bzdak  _et al._ (2020). The presence of the QCD CP should
manifest itself in the enhanced fluctuations of proton number Hatta and
Stephanov (2003) and possibly non-monotonic collision energy dependence of
non-Gaussian fluctuation measures Stephanov (2009, 2011). Measurements of
proton number fluctuations in nucleus-nucleus collisions have been performed
by different experiments such as STAR Adam  _et al._ (2021); Abdallah  _et
al._ (2021), HADES Adamczewski-Musch  _et al._ (2020), and ALICE Acharya  _et
al._ (2020). The measurements indicate a possible non-monotonic collision
energy dependence of the kurtosis of proton number Adam  _et al._ (2021) as
well as a possible enhancement of two-proton correlations over non-critical
baselines Vovchenko  _et al._ (2022) but conclusive evidence for the presence
of QCD CP is still lacking.

Figure 1: The liquid-gas region of the Lennard-Jones fluid phase diagram.
Horizontal dashed lines show the subcritical isotherm $\tilde{T}=0.76$ studied
in this work and the supercritical isotherm $\tilde{T}=1.06$ explored in Ref.
Kuznietsov  _et al._ (2022). Solid and dashed lines show the binodal and
spinodal lines, respectively. The blue and green regions correspond to the
nucleation and cavitation metastable parts of the mixed phase, respectively.
The grey area denotes the spinodal decomposition region. The black star
represents the CP. The squares denote the $(\tilde{n},\tilde{T})$ points where
the MD simulations in the mixed phase have been performed.

The grand canonical ensemble (GCE) of statistical mechanics is the most
suitable framework to study statistical fluctuations. Within this formulation,
the cumulants of particle number distribution are straightforwardly connected
to the chemical potential derivatives of thermodynamic potential. However, the
GCE can not be directly used for the conditions realized in the experiment
Koch (2010); Vovchenko (2021). Several essential restrictions should be taken
into account: (i) finite size of systems created in the experiment Berdnikov
and Rajagopal (2000); Poberezhnyuk  _et al._ (2020), (ii) influence of the
global conservation laws, for instance, baryon number conservation Bzdak  _et
al._ (2013); Vovchenko  _et al._ (2020a), and (iii) differences between
coordinate and momentum space acceptances. Recently the subensemble acceptance
method (SAM) to correct the fluctuation measurements for global conservation
laws has been developed Vovchenko  _et al._ (2020a, b); Poberezhnyuk  _et al._
(2020); Vovchenko (2022); Barej and Bzdak (2022a, b). This method is
applicable for statistical systems in the presence of interactions. In the
limit of ideal Maxwell-Boltzmann gas, it reduces to the binomial acceptance
correction procedure Bzdak  _et al._ (2013); Braun-Munzinger  _et al._ (2017);
Savchuk  _et al._ (2020).

In the present work we continue our studies Kuznietsov  _et al._ (2022) of
particle number fluctuations within molecular dynamics (MD) simulations of the
Lennard-Jones (LJ) fluid. The model considered here corresponds to an
interacting system of non-relativistic particles. The presence of both
attractive and repulsive interactions leads to a first-order liquid-gas phase
transition (LGPT). The MD simulations of the LJ fluid provide a microscopic
approach to fluctuations in a system with a phase transition. They also allow
one to study deviations from the baselines based on the GCE. This study thus
complements earlier analyses of correlations and fluctuations in the first-
order phase transition region performed using hadronic transport with mean
fields Sorensen and Koch (2021); Savchuk  _et al._ (2023) or fluid dynamics
with a finite-range term Steinheimer and Randrup (2012, 2013). With regard to
mean quantities the molecular dynamics of non-equilibrium finite systems was
studied previously in the context of heavy ion collisions in Refs. Molitoris
_et al._ (1984); Peilert  _et al._ (1988); Aichelin  _et al._ (1988); Peilert
_et al._ (1989); Hartnack  _et al._ (1998).

Our study is motivated by the measurements of baryon number fluctuations in
heavy-ion collisions to probe the QCD phase structure. In particular, the LJ
fluid can naturally model the nuclear liquid-gas transition between a dilute
gas of nucleons and clusters and the dense nuclear liquid, if one regards the
LJ particles as nucleon degrees of freedom. This nuclear LGPT is probed in
nuclear collisions at low energies Pochodzalla  _et al._ (1995); Natowitz  _et
al._ (2002); Karnaukhov  _et al._ (2003). Experiments at higher collision
energies, on the other hand, study the confinement-deconfinement transition,
which may contain a critical point and a line of first-order phase transition
at finite baryon density Stephanov  _et al._ (1999); Bzdak  _et al._ (2020).
The relevance of the LJ fluid to model the confinement-deconfinement
transition may seem less evident, given that it does not describe the expected
change of degrees of freedom from hadrons to quarks. Nevertheless, simulations
of the LJ fluid do provide useful guidance to understand the behavior of
baryon number fluctuations near the QCD CP, for two reasons: (i) the behavior
of baryon number fluctuations is universal near the QCD CP and governed by the
3D-Ising universality class Bzdak  _et al._ (2020) – the same universality
class that characterizes critical behavior in the LJ fluid Caillol (1998);
(ii) the LJ fluid simulations can test the validity of the model-independent
SAM procedure for subtracting the canonical ensemble effects on baryon number
cumulants, this is particularly relevant given that the finite-size effects,
that hinder the accuracy of the SAM, can be significant in the mixed phase
region.

This work focuses on fluctuations in the mixed-phase region of a first-order
phase transition. While significant attention has been given to higher-order
measures of fluctuations of conserved charges at supercritical temperatures
and in pure phases (see e.g. Refs. Vovchenko  _et al._ (2015a); Stephanov  _et
al._ (1999); Hatta and Ikeda (2003); Stephanov (2009, 2011); Mukherjee  _et
al._ (2017); Poberezhnyuk  _et al._ (2019); Motornenko  _et al._ (2020)), less
attention has been paid to the mixed phase. However, it is possible for a
system created in relativistic nucleus-nucleus collisions to enter the mixed
phase of a first-order phase transition under certain conditions. This is
particularly relevant because of the ongoing program of the HADES
collaboration at the GSI Helmholtzzentrum für Schwerionenforschung mbH to
measure higher-order net-proton and net-charge fluctuations in central Au+Au
collisions at collision energies of $0.2A-1.0A$ GeV. The system created in
these collisions may undergo freeze-out in the mixed phase of the nuclear
LGPT.

In our previous work Kuznietsov  _et al._ (2022), we studied a supercritical
isotherm, $T=1.06\,T_{c}$, observing a sizable increase of particle number
fluctuations near the critical particle number density $n\approx n_{c}$. In
the present work, we study particle number fluctuations along a subcritical
temperature $T=0.76\,T_{c}$ inside the liquid-gas mixed phase. First, we look
at the metastable part of the mixed phase – the so-called nucleation region.
The simulation results are compared to a simple model of non-interacting
particle clusters. Another part of the liquid-gas mixed phase – the spinodal
decomposition region – demonstrates anomalous large particle number
fluctuations. This happens at temperature $T$ and particle number density $n$
also far away from the CP. A simple analytical toy model is constructed to
clarify these effects.

Figure 2: The particle number distributions in a subvolume $V=\alpha V_{0}$
for the system with $N_{0}=400$ particles. The distributions $P(N)$ obtained
from the MD simulations at $\tilde{T}=0.76$ and different values of
$\tilde{n}$ inside the mixed phase.

The paper is organised as follows. The details of MD with LJ potential and the
results of the simulations for particle number fluctuations are presented in
Sec. II. A brief description of the mixed phase structure is described in Sec.
III. A simple model of non-interacting clusters in Sec. IV and a geometrical
toy model in Sec. V are developed to interpret the MD results in the
nucleation and spinodal decomposition regions, respectively. Summary in Sec.
VI closes the article

##  II Molecular dynamics with Lennard-Jones potential

We use molecular dynamics simulations of the classical non-relativistic system
of particles interacting via the Lenard-Jones (LJ) potential,

$\displaystyle V_{\rm
LJ}(r)=4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]\leavevmode\nobreak\
.$ (1)

The first term in Eq. (1) corresponds to the repulsive forces at short
distances whereas the second one describes the attractive interactions. The
parameter $\epsilon$ describes the depth of the attractive well, and $\sigma$
corresponds to the size of the particle, which also defines the distance
scale. It is convenient to introduce dimensionless reduced variables,

$\displaystyle V^{*}_{\rm LJ}(r^{*})=V_{\rm
LJ}(r)/\varepsilon=4\left((r^{*})^{-12}-(r^{*})^{-6}\right)\leavevmode\nobreak\
,$ (2)

with $r^{*}=r/\sigma$ being the reduced distance. The reduced thermodynamic
variables are the temperature $T^{*}=T/(\varepsilon)$, particle number density
$n^{*}=n\sigma^{3}$, and pressure $p^{*}=p\sigma^{3}/\varepsilon$. The
particle’s mass can be utilized to define the dimensionless time variable,
$t^{*}=t\sqrt{\varepsilon/(m\sigma^{2})}$.

Figure 3: The scaled variance $\tilde{\omega}$ (6) as a function of the
particle number density $\tilde{n}$ for $N_{0}=400$, $\alpha=0.2$ at
$\tilde{T}=0.76$ (a) and $\tilde{T}=1.06$ (b).

The LJ system possesses a rich phase diagram (see e.g. Ref. Stephan  _et al._
(2020) for an overview). At present, there are no direct analytical tools to
compute the phase diagram in the LJ system. Nevertheless, numerical methods
(see, e.g. Ref. Evans and da Gama (1979)) allow one to compute the approximate
locations of the LGPT binodal and spinodal lines, as well as the CP location.
This part of the phase diagram is of primary interest in the present work, and
it is shown in Fig. 1 in terms of the reduced temperature and density. The CP
location has been estimated from numerous MD simulations Stephan  _et al._
(2019)

$\displaystyle\begin{split}T^{*}_{c}&=1.321\pm 0.007\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \\\ n^{*}_{c}&=0.316\pm
0.005\leavevmode\nobreak\ ,\leavevmode\nobreak\ \leavevmode\nobreak\ \\\
p^{*}_{c}&=0.129\pm 0.005\leavevmode\nobreak\ .\end{split}$ (3)

In what follows, we use a set of dimensionless variables scaled by the
critical values

$\displaystyle\tilde{T}\leavevmode\nobreak\
\equiv\frac{T}{T_{c}}=\frac{T^{*}}{T^{*}_{\rm c}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\tilde{n}\equiv\frac{n}{n_{c}}=\frac{n^{*}}{n^{*}_{\rm c}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\tilde{p}\equiv\frac{p}{p_{c}}=\frac{p^{*}}{p^{*}_{\rm c}}\leavevmode\nobreak\
.$ (4)

The quantities (3) correspond to the thermodynamic limit when the system’s
volume $V\rightarrow\infty$. For finite systems, the physical meaning of the
LGPT and its CP should be treated with caution, as they are only rigorously
defined for infinite systems.

The MD simulations are performed by numerically integrating Newton’s equations
of motion using the Velocity Verlet integration method. The simulations are
done for a system of $N_{0}=400$ interacting particles in a cubic box of
volume $V_{0}$ with periodic boundary conditions with minimum image
convention.

In the mixed phase the time of reaching the thermal equilibrium can be rather
large (see Refs. Sanz  _et al._ (2013)). After the equilibration time,
$\tilde{t}_{\rm eq}=50$, the LJ system reaches a state with a stable
temperature111During all system evolution some temperature fluctuations can be
seen, but they are relatively small, so the mean value of temperature differs
from the desired by no more than $0.4\%$. (see Ref. Kuznietsov  _et al._
(2022)). The time of all simulations is $\tau=10^{6}$. This large time
interval guarantees small deviations (less than 1%) of the scaled variance in
independent simulations.

The total particle number $N_{0}$ in the entire volume is fixed. To study the
fluctuations of particle number one thus needs to choose a subvolume $V=\alpha
V_{0}$ ($0<\alpha<1$) of the whole volume. We choose a cubic subvolume placed
in the geometrical center of the system. From the MD simulations, we obtain
the normalized probability distribution $P(N)$ to observe $N$ particles in the
subvolume $V$.

Figure 4:  Different regions along the supercritical isotherm of the liquid-
gas transition: (a) gaseous phase, (b) nucleation, (c) spinodal decomposition,
and (d) cavitation.  Figure 5:  Possible position of the liquid phase in the
spinodal decomposition region relative to the acceptance subvolume (red dashed
square).

A useful measure of particle number fluctuations is the scaled variance:

$\displaystyle\omega=\frac{\left<N^{2}\right>-\left<N\right>^{2}}{\left<N\right>}.$
(5)

In MD simulations, the values $\left<N\right>$ and $\left<N^{2}\right>$ can be
calculated as time averages. In Fig. 2 we present the $P(N)$ distribution at
the subcritical temperature $\tilde{T}=0.76$ for several different particle
number densities $\tilde{n}$ inside the mixed phase. The total number of
particles is fixed as $N_{0}=400$ and the subvolume fraction is taken as
$\alpha=0.2$. From Fig. 2, one observes substantial deviations of the
resulting distributions from the Poisson distribution baseline. For
$\tilde{n}\approx 1$, a double-hump distribution is clearly observed.

Note that for any finite $\alpha$, fluctuations of $N$ in the subvolume $V$
will be influenced by the exact global conservation of the total particle
number $N_{0}$ in the full volume $V_{0}$. In the large volume limit, these
effects can be taken into account analytically Vovchenko  _et al._ (2020a).
One can defined a scaled variance $\tilde{\omega}$ corrected for exact $N_{0}$
conservation as

$\displaystyle\tilde{\omega}=\frac{\omega}{1-\alpha}\leavevmode\nobreak\ .$
(6)

The results for the corrected scaled variance $\tilde{\omega}$ as a function
of $\tilde{n}$ are presented in Fig. 3 for both (a) the subcritical and (b)
the supercritical isotherms $\tilde{T}=0.76$ and $\tilde{T}=1.06$,
respectively. All results are obtained for $N_{0}=400$ and $\alpha=0.2$, as in
Fig. 2.

One can immediately observe that fluctuations are much larger in the mixed
phase at $\tilde{T}=0.76$ compared to those along the temperature
$\tilde{T}=1.06$, slightly above the critical point. This indicates that,
although the fluctuations exhibit singular behavior in the vicinity of the CP,
they can be even larger in the mixed phase region away from the critical
point.

In the following sections, we provide a brief overlook of the structure of the
liquid-gas mixed phase and analyze the observed large values of
$\tilde{\omega}$ in the mixed phase in terms of simple analytical models.

##  III Mixed phase structure

One can specify three different regions inside the mixed phase: nucleation,
spinodal decomposition, and cavitation (see, e.g., Refs. Frenkel (1946);
Fisher (1967)). They are shown in Fig. 1 by blue, grey, and green colors,
respectively. Their microscopic structures are symbolically illustrated in
Fig. 4. The nucleation region includes a mixture of particles and small
clusters (liquid droplets), whereas the cavitation region is represented by
the liquid with small bubbles of the gaseous phase. In the context of heavy
ion collision clusters correspond to nuclear fragments whose distributions
were previously studied using MD in the case of expanding system in Refs.
Molitoris  _et al._ (1984); Stöcker and Greiner (1986); Peilert  _et al._
(1988); Aichelin  _et al._ (1988); Peilert  _et al._ (1989). Experimental
measurements of nuclear fragment mass distributions were used to probe the
nuclear LGPT and the CP (see, e.g., Refs. Lindenstruth (1993); Lindenstruth
_et al._ (1993); Pochodzalla  _et al._ (1995); Pochodzalla (1997)). The
nucleation and cavitation regions of the mixed phase correspond to the
metastable states. In the MD simulations one expects to achieve an
equilibrated steady state in these regions after a sufficiently long time. In
most cases, however, the time to reach complete equilibrium appears very long.
Note also that a strict physical meaning and location of the bounds of
different regions are dependent on the size of the system (see, e.g., Refs.
Klein and Unger (1983); Csernai and Neda (1994); Wedekind  _et al._ (2009))
and can be sensitive to the collective motion Pratt (1990); Kunde  _et al._
(1995); Pratt  _et al._ (1995).

The spinodal decomposition region is fundamentally different from the
metastable nucleation and cavitation ones (see, e.g., Refs. Lopez and Luebeck
(1989); Favvas and Mitropoulos (2008)). The LGPT manifests itself here as a
fast system separation into the gaseous and liquid phases. The equilibrium
states in this region (see, e.g., Ref. Elliott (1989)) are achievable in the
MD simulations. The heterogeneous structure of the spinodal decomposition
phase is illustrated in Fig. 5, showing a strong influence on the particle
number fluctuations obtained in the MD simulations. This is discussed in more
detail in Sec. V. One can note a principal difference between the
heterogeneous two-phase states in the spinodal decomposition region and the
homogeneous mixtures of particles plus clusters in the nucleation region and
liquid with gaseous bubbles in the cavitation region.

##  IV Mixture of particles and Clusters

Figure 6: Cluster probability distributions $P_{k}$ extracted from the MD of
Lennard-Jones fluid for $N_{0}=400$ and $\tilde{T}=0.76$ at gaseous binodal
$\tilde{n}=0.16$ (a) and gaseous spinodal $\tilde{n}=0.35$ (b). For
comparison, $P_{k}$ distributions are also presented for supercritical
temperature $\tilde{T}=1.06$.

To clarify some general features of the nucleation region, let us consider a
non-interacting multi-component gas of $k$-particle clusters ($k=1,2,\ldots$).
The GCE partition function reads

$\displaystyle\begin{split}&Z_{\rm GCE}=\prod_{k\geq
1}\sum_{N_{k}=0}^{\infty}\frac{\left(Vg(k)e^{\mu
k/T}\right)^{N_{k}}}{N_{k}!}(2\pi kmT)^{3N_{k}/2}\\\ &=\prod_{k\geq
1}\exp\left[V(2\pi km\,T)^{3/2}\,g(k)\,\exp\left(\frac{\mu
k}{T}\right)\right],\end{split}$ (7)

where $V$, $T$, and $\mu$ are, respectively, the system volume, temperature,
and chemical potential that corresponds to the total conserved number $N$ of
particles over all clusters; $g(k)$ is the “degeneracy” factor (number of
internal states of the $k$-th cluster), and $m$ the mass of a single particle,
such that the mass of a $k$-particle cluster equals $M_{k}=km$). The system is
considered to be in chemical equilibrium, thus $\mu_{k}=k\mu$. The CE
partition function $Z_{\rm CE}(V,T,N)$ of the cluster model (7) is considered
in Appendix A, where it is shown that the moments $\langle k^{l}\rangle$ of
the cluster distribution are identical between the CE and the GCE in the
thermodynamic limit.

The cluster distribution (i.e., the normalized probability to find the $k$th
cluster in the cluster system) can be written in a form

$\displaystyle P_{k}(T,\mu)\equiv\frac{\left<N_{k}\right>}{\sum\limits_{l\geq
1}\left<N_{l}\right>}=\frac{k^{3/2}g(k)\exp\left(\frac{\mu
k}{T}\right)}{\sum\limits_{l\geq 1}l^{3/2}g(l)\exp\left(\frac{\mu
l}{T}\right)}\leavevmode\nobreak\ ,$ (8)

where

$\displaystyle k\langle N_{k}\rangle=\frac{\partial\ln\left[Z_{\rm
GCE}(k)\right]}{\partial\mu}$ (9)

is the GCE average number of the $k$th clusters. The clusters pressure $p$ and
particle number density $n$ can be found as

$\displaystyle p$ $\displaystyle=(2\pi m)^{3/2}T^{5/2}\sum_{k\geq
1}k^{3/2}g(k)\exp\left(\frac{\mu k}{T}\right),$ (10) $\displaystyle n$
$\displaystyle=(2\pi mT)^{3/2}\sum_{k\geq 1}k^{5/2}g(k)\exp\left(\frac{\mu
k}{T}\right).$ (11)

Using Eqs. (8) and (10) one can rewrite the pressure as

$\displaystyle p=\frac{nT}{\langle k\rangle}\leavevmode\nobreak\ ,$ (12)

and the scaled variance $\omega_{\rm gce}$

$\displaystyle\omega_{\rm gce}$
$\displaystyle=T\left[\frac{dp}{dn}\right]^{-1}=\frac{T}{n}\left(\frac{\partial
n}{\partial\mu}\right)_{T}=\frac{\langle k^{2}\rangle}{\langle
k\rangle}\leavevmode\nobreak\ ,$ (13)

where we defined $\langle k^{l}\rangle\equiv\sum_{k\geq 1}k^{l}P_{k}$.
Therefore, the first two moments of the cluster probability distribution
$P_{k}$ define both the system pressure (12) and scaled variance (13). Due to
the evident inequalities, $\langle k\rangle\geq 1$ and $\langle
k^{2}\rangle\geq\langle k\rangle$, the results (12) and (13) demonstrate that
in the mixture of noninteracting $k$-th clusters ($k=1,2,\ldots$) the system
pressure becomes smaller and the scaled variance larger than the corresponding
ideal gas values $p_{\rm id}=nT$ and $\omega_{\rm id}=1$ with no cluster
formation, i.e., when $g(k=1)=1,\leavevmode\nobreak\ g(k>1)=0$. General
expression for cumulants $\kappa_{n}[N]$ of any order $n$ can be obtained:

$\displaystyle\kappa_{n}[N]=\frac{\partial^{n}\ln\left[Z_{\rm
GCE}\right]}{\partial\left(\mu\over T\right)^{n}}=\langle
k^{n}\rangle\sum\limits_{k\geq 1}\langle N_{k}\rangle.$ (14)

The model of noninteracting clusters discussed above can be considered as an
approximation for the LJ fluid in the nucleation region. The attractive part
of the LJ potential is responsible for the $k$th cluster formation. On the
other hand, the particle number density is still sufficiently small to justify
the absence of the repulsive interaction effects between clusters.

By definition, a cluster is a bound system of particles. There are several
ways to define clusters in molecular dynamics simulations Sator (2003);
Moretto  _et al._ (2011). In the following, we will use the Hill algorithm
Hill (1955). A pair of particles $i$ and $j$ is assumed to be bound if their
rest frame energy is negative,

$\displaystyle(\tilde{v}_{i}-\tilde{v}_{j})^{2}+\tilde{V}_{LJ}(|\tilde{r}_{i}-\tilde{r}_{j}|)<0.$
(15)

A given particle is assumed to belong to a cluster if it is bound to at least
one other particle in that cluster. Finding clusters is thus equivalent to
finding connected components in an undirected graph, whose vertices correspond
to particles and where all bound pairs of particles [i.e. the condition (15)
is satisfied] are connected by edges. We use depth-first search (DFS) to find
the connected components of the graph and thus identify all the clusters.

Utilizing the above procedure, one obtains the probability distribution
$P_{k}$ in a Lennard-Jones fluid for given $\tilde{n}$ and $\tilde{T}$ from MD
simulations. Examples of the extracted $P_{k}$ distributions for
$\tilde{T}=0.76$ and $\tilde{T}=1.06$ are shown in Figs. 6 (a) and (b) for
$\tilde{n}=0.16$ and $\tilde{n}=0.35$, respectively. The results indicate that
cluster formation becomes more significant when either temperature $\tilde{T}$
is decreased or particle number density $\tilde{n}$ is increased.

We then use the extracted $P_{k}$ distributions to evaluate $\langle k\rangle$
and $\langle k^{2}\rangle$ which we then plug into (13) to estimate the GCE
scaled variance in the cluster model. These results are compared with
$\tilde{\omega}$ calculated in a subvolume $V=\alpha V_{0}$ directly from MD
simulations. The cluster model results for $\tilde{T}=0.76$ are shown in Fig.
7 by the orange line. These results agree qualitatively with direct MD
simulations data (black line) in the range of densities
$0.16\lesssim\tilde{n}\lesssim 0.35$ corresponding to the nucleation region.
In particular, cluster formation explains the strong rise (approximately by a
factor of 20) of the scaled variance with $\tilde{n}$ in the nucleation
region. The cluster model, however, fails to describe the peak in
$\tilde{\omega}$ seen in MD simulations at higher densities, indicating its
breakdown in the spinodal region.

Figure 7:  The points connected by the solid line correspond to the MD results
for $N_{0}=400$ and $\alpha=0.2$ at $\tilde{T}=0.76$. The orange line
demonstrates the cluster model results in the nucleation region
$0.16\leq\tilde{n}\leq 0.35$. The dashed line shows the results of the
Minecraft model in the spinodal region $0.35\leq\tilde{n}\leq 1.75$, and the
dash-dotted line is its extension to the nucleation region.

##  V Fluctuations in the spinodal region

In Ref. Poberezhnyuk  _et al._ (2021), the GCE particle number fluctuations
were calculated in the mixed phase region. It was assumed that both the liquid
and gas phases are entirely inside the system volume $V_{0}$ that tends to
infinity. In MD simulations here, we instead study fluctuations in a subvolume
$V=\alpha V_{0}$, which corresponds to a different scenario. We thus develop
new models to understand qualitative features of the behavior observed in MD
simulations.

In the spinodal region, one assumes that the volume $V$ is partitioned into
volumes $V_{l}=xV$ and $V_{g}=yV$ occupied by the liquid and gaseous phases,
respectively ($0<x<1$, $y\equiv 1-x$). The corresponding particle number
densities in the liquid and gaseous phases are $\rho_{l}\equiv N_{l}/V_{l}$
and $\rho_{g}=N_{g}/V_{g}$. The $r$th moment of the particle number
distribution in the subvolume $V=\alpha V_{0}$ can then be presented as the
following:

$\displaystyle\langle
N^{r}\rangle=\langle(N_{l}+N_{g})^{r}\rangle=V^{r}\,\langle(x\rho_{l}+y\rho_{g})^{r}\rangle\leavevmode\nobreak\
.$ (16)

The fluctuating quantities are the densities $\rho_{l}$, $\rho_{g}$, and the
volume fraction $x$, whereas the volume $V$ is fixed. Following Refs.
Vovchenko  _et al._ (2015b); Satarov  _et al._ (2021) we assume that the
fluctuations of all these quantities are independent in the thermodynamic
limit, i.e.,
$\langle\rho_{l}^{k}\,\rho_{g}^{m}\,x^{n}\rangle=\langle\rho_{l}^{k}\rangle\,\langle\rho_{g}^{m}\rangle\,\langle
x^{n}\rangle$ for any non-negative integers $k$, $m$, and $n$.

The first moment ($r=1$), reduces via Eq. (16) to

$\displaystyle\langle N\rangle=x_{0}Vn_{l}+y_{0}Vn_{g}=Vn\leavevmode\nobreak\
,$ (17)

where $x_{0}=\langle x\rangle$ is the mean volume fraction occupied by the
liquid phase, $y_{0}\equiv 1-x_{0}$, and $n_{l}=\langle\rho_{l}\rangle$ and
$n_{g}=\langle\rho_{g}\rangle$ are the mean densities in the liquid and
gaseous phases, respectively. The particle number density is equal to
$n\equiv\langle N\rangle/V=N_{0}/V_{0}$. Equation (17) defines $x_{0}$ in
terms of the mean densities:

$\displaystyle x_{0}\equiv\langle
x\rangle=\frac{n-n_{g}}{n_{l}-n_{g}}\leavevmode\nobreak\ .$ (18)

At fixed temperature $T<T_{c}$ the mean densities of the liquid $n_{l}$ and
gaseous $n_{g}$ phases are assumed to remain constant with respect to system’s
particle number density $n$ in the spinodal region in the thermodynamic limit.
These quantities coincide with the corresponding values on the liquid (right)
and gaseous (left) binodals.

Using Eq. (16) one finds the variance of particle number distribution (see
Ref. Poberezhnyuk  _et al._ (2021) for details):

$\displaystyle{\rm Var}[N]\equiv\langle N^{2}\rangle-\langle N^{2}\rangle$
$\displaystyle={\rm Var}_{x}[N_{l}]\left(1+\frac{{\rm
Var}[x]}{x_{0}^{2}}\right)+{\rm Var}_{x}[N_{g}]\left(1+\frac{{\rm
Var}[x]}{y_{0}^{2}}\right)$ $\displaystyle\quad+V^{2}(n_{l}-n_{g})^{2}{\rm
Var}[x]\leavevmode\nobreak\ .$ (19)

Here ${\rm Var}_{x}[N_{l,g}]$ is the variance of $N_{l,g}$ at fixed volume
fraction $x$ and ${\rm Var}[x]$ is the variance of the $x$ distribution.

Suppose that there are several blobs of the liquid and gaseous phases, and all
of them are much smaller than the subvolume $V$. This would correspond to a
spatially homogeneous mixed phase. In this case, ${\rm Var}[x]$ is expressed
in terms of cumulants of $V_{l}$ distribution as ${\rm Var}[x]\equiv
V^{-2}{\rm Var}[V_{l}]$. In the thermodynamic limit, $V\rightarrow\infty$, all
cumulants of extensive quantities are proportional to the system volume, ${\rm
Var}_{x}[N_{\rm{l,g}}]\sim V$ and ${\rm Var}_{x}[V_{l,g}]\sim V$. Eq. (19)
reduces to

$\displaystyle{\rm Var}[N]$ $\displaystyle={\rm Var}_{x}[N_{l}]+{\rm
Var}_{x}[N_{g}]$ $\displaystyle+V^{2}(n_{l}-n_{g})^{2}{\rm
Var}[x]\leavevmode\nobreak\ ,$ (20)

where all terms are linear in $V$. The result (21) coincides with that
obtained for the GCE in Ref. Poberezhnyuk  _et al._ (2021), and it corresponds
to the finite values of the scaled variance at $T<T_{c}$ in the thermodynamic
limit.

Note that the above derivation is based on the assumption of homogeneity. This
assumption is valid for pure phases. In the mixed-phase region, however, this
assumption may only be reasonable when applied to long-lived metastable
phases. Such a configuration of the system, however, can not be viewed as an
equilibrium configuration in a region of spinodal decomposition. There, the
sizes of the liquid and gaseous blobs are both of the order of the total
volume $V_{0}$, and their volumes are comparable to the subvolume $V$. Thus,
the whole picture is strongly heterogeneous (see Fig. 4 (c) and Fig. 5). As a
consequence, ${\rm Var}[x]$ becomes volume independent, thus, the last term in
Eq. (19) is quadratic in $V$ and makes the dominant contribution to
fluctuations. Leaving only this last term, one obtains:

$\displaystyle{\rm Var}[N]=V^{2}(n_{l}-n_{g})^{2}{\rm
Var}[x]\leavevmode\nobreak\ ,$ (21)

and

$\displaystyle\begin{split}\tilde{\omega}[N]&=\frac{{\rm
Var}[N]}{(1-\alpha)\langle N\rangle}\\\
&=\alpha(1-\alpha)N_{0}\frac{(n_{l}-n_{g})^{2}}{n^{2}}{\rm
Var}[x]\leavevmode\nobreak\ .\end{split}$ (22)

This result indicates that $\tilde{\omega}[N]$ scales with $N_{0}$, i.e. the
scaled variance diverges in the thermodynamic limit. We checked that for
$N_{0}\gg 400$ the substantial increase of $\tilde{\omega}$ is observed within
MD simulations, however the scaling behaviour for fluctuations is out of the
scope of the present paper. In the following, we present estimates for ${\rm
Var}[x]$.

Small $\alpha$ limit. At $\alpha\ll 1$ one has $V_{\rm l}\gg V$ and $V_{\rm
g}\gg V$. This means that one can neglect the events when both phases are
simultaneously present inside the subvolume $V$, and the whole subvolume is
entirely inside either the gaseous or liquid phase. The probability
distribution $P[x]$ thus reads

$\displaystyle P[x]=x_{0}\,\delta(1-x)+y_{0}\,\delta(x)\leavevmode\nobreak\ .$
(23)

This means that one can neglect the events when both phases are simultaneously
present inside the subvolume $V$. From Eq. (23) one finds

$\displaystyle{\rm Var}[x]=x_{0}y_{0}.$ (24)

The maximal value of ${\rm Var}[x]=0.25$ is reached at $x_{0}=\sqrt[3]{0.5}$.
Using Eqs. (18) and (21) one obtains:

$\displaystyle{\rm Var}[N]=V^{2}(n-n_{g})(n_{l}-n)\leavevmode\nobreak\ .$ (25)

One sees that the scaled variance of particle number distribution is indeed
divergent inside the mixed phase in the thermodynamic limit, scaling with the
subvolume $\tilde{\omega}\sim V$.

Minecraft model.222This name is inspired by the popular 3D video game. Now let
us calculate ${\rm Var}[x]$ when the sizes of the volume, subvolume, and blobs
are all comparable. For that we consider a simple ”geometric” toy model of the
cubic system with unit volume which contains both liquid and gaseous phases
(see Fig. 8). The cubic subvolume $V=\alpha$ is located in the center of the
system with coordinates $(a_{\rm x},a_{\rm y},a_{\rm z})=(0,0,0)$. The edge
length of the subvolume is $a=\sqrt[3]{\alpha}$. All liquid is condensed into
a single blob which freely moves within the system. Here we neglect the
effects of a geometric form and assume that this blob has a shape of a perfect
cube. The volume of the cube of liquid is $V_{l}=x_{0}$. Correspondingly its
edge length is $b=\sqrt[3]{x_{0}}$. The system has periodic boundary
conditions, therefore, the coordinates $(b_{x},b_{y},b_{z})$ of the center of
the cube of liquid are limited by $-\frac{1}{2}<b_{\rm x},b_{\rm y},b_{\rm
z}<\frac{1}{2}$. The fraction $x$ of the subvolume occupied by the liquid
phase is the overlap volume between the cubic subvolume and the cubic liquid
divided by the subvolume $V=\alpha$.

Figure 8: The illustration of the Minecraft toy model of an equilibrium system
in the unstable region of the mixed phase. The subsystem is shown by the grey
cube in the center while the green cube represents the liquid. The remaining
space of the system is occupied by gas.

The system has three degrees of freedom – the coordinates of the liquid cube
$(b_{x},b_{y},b_{z})$. Since the cube center is uniformly distributed over
$-\frac{1}{2}<b_{\rm x},b_{\rm y},b_{\rm z}<\frac{1}{2}$, the three
coordinates are independent. The fraction $x$ as a function of these three
coordinates and can be written as

$\displaystyle x=\frac{f(b_{\rm x})f(b_{\rm y})f(b_{\rm z})}{\alpha}.$ (26)

Here $f(b_{\rm i})$ is the overlap of liquid blob with subvolume along the
coordinate $i$ as a function of $b_{\rm i}$.

The mean value $\langle x\rangle$ can be found as

$\displaystyle\langle x\rangle=\frac{1}{\alpha
v}\left(\int_{-1/2}^{1/2}f(b_{\rm i}){\rm d}b_{\rm i}\right)^{3}=x_{0}$ (27)

where $v=1$ is the volume of the system. Similarly, one can calculate the
variance of $x$:

$\displaystyle{\rm Var}[x]=-b^{6}$ (28)
$\displaystyle+\left[\frac{b^{2}(3a-b)+\Theta_{a+b-1}(a+b-1)^{3}+\Theta_{b-a}(b-a)^{3}}{3a^{2}}\right]^{3}$

where $\Theta_{...}\equiv\Theta[...]$ is the step function and, as before,
$a=\sqrt[3]{\alpha}$ and $b=\sqrt[3]{x_{0}}$. One sees that Eq. (28) reduces
to Eq. (24) when $\alpha\rightarrow 0$. In other limiting cases ${\rm
Var}[x]\rightarrow 0$ when $\alpha\rightarrow 1$ or $x_{0}\rightarrow 1$ or
$x_{0}\rightarrow 0$. ${\rm Var}[x]$ as a function of $x_{0}$ and $\alpha$ is
shown in Fig. 9.

The scaled variance $\tilde{\omega}[N]$ given by Eq. (22), with ${\rm Var}[x]$
estimated using the Minecraft model, Eq. (28), is shown in Fig. 7 in spinodal
and nucleation regions by dashed and dotted lines, respectively.

Figure 9: The variance of the volume fraction occupied by the liquid phase,
${\rm Var}[x]$, as a function of $\langle x\rangle\equiv x_{0}$ and $\alpha$
calculated in the Minecraft model [Eq. (28)]. The dashed line corresponds to
the maximum value of ${\rm Var}[x]$ at fixed $\alpha$.

##  VI Summary

We studied particle number fluctuations inside the mixed phase of a liquid-gas
phase transition by utilizing molecular dynamics simulations of the Lennard-
Jones fluid. The simulations were performed for $N_{0}=400$ particles in a
cubic box with periodic boundary conditions. The fluctuations are studied
inside a cubic subvolume $V=0.2\,V_{0}$ located in the geometrical center of
the system.

First, we briefly explore the supercritical temperature, where one observes
the approximate Gaussian shape of the $P(N)$ distribution. The scaled variance
$\tilde{\omega}$ characterizes the width of the $P(N)$ distribution. It first
increases with density from $\tilde{\omega}\approx 1$ at small $\tilde{n}$ to
its maximum above unity around the critical density $\tilde{n}=1$, and then it
decreases with $\tilde{n}$ to small values $\tilde{\omega}<1$. This is
illustrated in Fig. 3 (b).

The situation differs in the mixed phase, $\tilde{T}<1$. The structure of the
$P(N)$ distribution is significantly more intricate. For $\tilde{n}\approx 1$,
the distribution is bi-modal, as shown in Fig. 2. The corresponding variance
of particle number is much more significant compared to pure phases [Fig. 3
(a)].

To understand the qualitative features of the observed behavior, we formulate
two phenomenological toy models.

The first model describes the system as non-interacting multi-component gas of
$k$-particle clusters, taking the cluster probability distribution $P_{k}$
directly from the MD simulations as input. The model describes semi-
quantitatively the rapid increase of $\tilde{\omega}$ with density in the
nucleation region of the mixed phase, i.e., the region between the gaseous
binodal and spinodal [Fig. 7].

The second model – the Minecraft model – is formulated for the spinodal region
of the mixed phase. The particles are separated into two phases, namely, the
liquid blob with volume $V_{\rm{l}}$ surrounded by gas. The size of the blob
$V_{\rm{l}}$ can be expressed through the total density of the system
$\tilde{n}$ and densities on the binodals. The Minecraft model considers the
geometrical effects that become important when the volumes $V_{\rm{l,g}}$ and
$V$ are of comparable size. With this consideration, the model indicates that
$\tilde{\omega}\sim N_{0}\rightarrow\infty$, thus the variance is divergent in
the thermodynamic limit inside the spinodal region.

The present work is motivated by the study of event-by-event fluctuations in
nucleus-nucleus collisions to probe the phase structure of QCD. Our MD
simulations inside the mixed phase were performed for $400$ particles, while
the fluctuations were studied in the subvolume $V=0.2\,V_{0}$. These two
parameters correspond to typical total numbers of nucleons and the percentage
of accepted final particles in heavy-ion collisions. The results indicate that
large fluctuations of particle number in coordinate space can be interpreted
as a signal of the spinodal region of the first-order phase transition.
However, there are significant differences between our calculations and heavy-
ion collisions. One difference is that in heavy-ion collisions, particles are
not detected during the equilibrium phase of the collision but only after they
fly away to the detector. Another difference is that particle momenta, not the
coordinates, are measured in the experiment. We plan to address these issues
by performing MD simulations of expanding systems.

Our simulations provide a first microscopic model test of the subensemble
acceptance method (SAM) Vovchenko  _et al._ (2020a, b) in the mixed phase
region of a first-order phase transition. The SAM is a method for correcting
the baryon number cumulants in heavy-ion collisions, which is model-
independent in the thermodynamic limit, and it was previously tested in the
crossover region at supercritical temperatures in Ref. Kuznietsov  _et al._
(2022). Our simulations reveal that the SAM remains accurate in metastable
regions of the phase diagram but breaks down in the spinodal decomposition
region. The reason is that the finite-size effects remain sizable even in
large systems in this region of the phase diagram. The treatment of the
canonical effects in the spinodal region is thus more complex. It will require
appropriate generalizations of the SAM, such as including macroscopic
geometrical effects encompassed by the Minecraft model introduced here.

Another future avenue is generalizing the presented analysis to higher-order
moments of particle number distributions, such as skewness and kurtosis.

## Acknowledgements

The authors are thankful to Jeroen van den Brink, Volker Koch, Flavio
Nogueira, Scott Pratt and Jan Steinheimer for fruitful comments and
discussions. O.S. acknowledges the scholarship grant from the GET$\\_$INvolved
Programme of FAIR/GSI and support by the Department of Energy Office of
Science through grant no. DE-FG02-03ER41259. This work is supported by the
National Academy of Sciences of Ukraine, Grant No. 0122U200259. M.I.G. and
R.V.P. acknowledge the support from the Alexander von Humboldt Foundation.
This work was supported by a grant from the Simons Foundation (Grant Number
1039151). H.St. appreciates the Judah M. Eisenberg Professur Laureatus of the
Walter Greiner Gesellschaft/Förderverein für physikalische Grundlagenforschung
Frankfurt, and the Fachbereich Physik at Goethe Universität.

##  Appendix A Cluster partition function in the CE

For a system of non-interacting multi-component gas of $k$th particle
clusters, the canonical ensemble (CE) partition function reads

$\displaystyle Z_{\rm CE}(V,T,N)$
$\displaystyle=\prod_{k=1}^{N}\sum_{N_{k}\geq
0}\frac{\left(Vg(k)\right)^{N_{k}}}{N_{k}!}(2\pi kmT)^{3N_{k}/2}$
$\displaystyle\times\delta\left[N-\sum_{k=1}^{N}kN_{k}\right]\leavevmode\nobreak\
.$ (29)

Applying the integral form of the Kronecker symbol,

$\displaystyle\delta\left[N-\sum_{k=1}^{N}kN_{k}\right]=\int\limits_{0}^{2\pi}\frac{d\varphi}{2\pi}\exp\left[i\varphi\left(N-\sum_{k=1}^{N}kN_{k}\right)\right]\leavevmode\nobreak\
,$ (30)

to Eq. (29), one obtains

$\displaystyle Z_{\rm
CE}(V,T,N)=\int\limits^{2\pi}_{0}\frac{d\varphi}{2\pi}\leavevmode\nobreak\
e^{-i\varphi N}\exp\left[\sum_{k\geq 1}r(k)e^{i\varphi k}\right].$ (31)

Here $r(k)\equiv Vg(k)(2\pi kmT)^{3/2}$.

Using the Maclaurin expansion, one has

$\displaystyle\exp\left[\sum_{k\geq 1}r(k)e^{i\varphi
k}\right]=\sum_{l=0}^{\infty}\frac{B_{l}(r(1),\dots,l!r(l))}{l!}e^{i\varphi
l},$ (32)

where $B_{l}$ are Bell polinomials Bell (1927).

Substituting (32) into (31) gives

$\displaystyle Z_{\rm
CE}(V,T,N)=\frac{B_{N}(r(1),\dots,N!r(N))}{N!}\leavevmode\nobreak\ .$ (33)

From the above equations one finds the GCE partition function

$\displaystyle Z_{\rm GCE}=\sum_{N=0}^{\infty}Z_{\rm CE}\exp\left(\frac{\mu
N}{T}\right)=\prod_{k=1}^{N}\exp\left(r(k)e^{\mu
k/T}\right)\leavevmode\nobreak\ ,$ (34)

which coincides with Eq. (7). $Z_{\rm CE}$ can be expressed in terms of
$Z_{\rm GCE}$ through the Mellin transformation

$\displaystyle Z_{CE}=\int\limits^{c+i\infty}_{c-i\infty}Z_{GCE}e^{-\mu
N/T}d\mu\leavevmode\nobreak\ .$ (35)

The integral (35) can be evaluated in the large $N$ limit using the steepest
descent method Debye (1909). Therefore,

$\displaystyle Z_{\rm CE}(V,T,N)\leavevmode\nobreak\ $
$\displaystyle\approx\leavevmode\nobreak\ \sqrt{\frac{2\pi
T^{2}}{\sum_{k=1}^{N}k^{2}r(k)e^{\mu_{0}k/T}}}$
$\displaystyle\times\leavevmode\nobreak\
\exp\left(\sum_{k=1}^{N}r(k)e^{\mu_{0}k/T}-\frac{\mu_{0}}{T}N\right)\leavevmode\nobreak\
,$ (36)

where $\mu_{0}(T,N)$ can be found from the saddle point equation

$\displaystyle\sum_{k=1}^{N}kr(k)e^{\mu_{0}k/T}-N=0\leavevmode\nobreak\ .$
(37)

Equation (36) indicated that the $j$-th moment of $k$th cluster distribution
in the large $N$ limit reads

$\displaystyle\left<k^{j}\right>_{\rm CE}\leavevmode\nobreak\
=\leavevmode\nobreak\ \left<k^{j}\right>_{\rm GCE}\leavevmode\nobreak\
+\leavevmode\nobreak\ O(N^{-1})\leavevmode\nobreak\ .$ (38)

This result shows that all moments $j=1,2,\ldots$ of the $k$-th cluster
distribution $(k=1,\ldots,N)$ are the same in the CE and GCE in the
thermodynamic limit $N\rightarrow\infty$.

If $N$ is large, this justifies the use of $P_{k}$ probabilities from MD
simulations as input into the calculations of fluctuations in the GCE using
formulas (12) and (13).

## References

  * Landau and Lifshitz (1975) L. D. Landau and E. M. Lifshitz,  _Statistical Physics_ (Pergamon, Oxford, 1975).
  * Greiner  _et al._ (2012) W. Greiner, D. Rischke, L. Neise, and H. Stöcker,  _Thermodynamics and Statistical Mechanics_, Classical Theoretical Physics (Springer New York, 2012).
  * Allen and Tildesley (2017) M. P. Allen and D. J. Tildesley,  _Computer simulation of liquids_ (Oxford university press, 2017).
  * Stephanov  _et al._ (1999) M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. D 60, 114028 (1999), arXiv:hep-ph/9903292 .
  * Bzdak  _et al._ (2020) A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, Phys. Rept. 853, 1 (2020), arXiv:1906.00936 [nucl-th] .
  * Bernhardt  _et al._ (2022) J. Bernhardt, C. S. Fischer, and P. Isserstedt, (2022), arXiv:2208.01981 [hep-ph] .
  * Kuznietsov  _et al._ (2022) V. A. Kuznietsov, O. Savchuk, M. I. Gorenstein, V. Koch, and V. Vovchenko, Phys. Rev. C 105, 044903 (2022), arXiv:2201.08486 [hep-ph] .
  * Hatta and Stephanov (2003) Y. Hatta and M. A. Stephanov, Phys. Rev. Lett. 91, 102003 (2003), [Erratum: Phys.Rev.Lett. 91, 129901 (2003)], arXiv:hep-ph/0302002 .
  * Stephanov (2009) M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009), arXiv:0809.3450 [hep-ph] .
  * Stephanov (2011) M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011), arXiv:1104.1627 [hep-ph] .
  * Adam  _et al._ (2021) J. Adam _et al._ (STAR), Phys. Rev. Lett. 126, 092301 (2021), arXiv:2001.02852 [nucl-ex] .
  * Abdallah  _et al._ (2021) M. Abdallah _et al._ (STAR), Phys. Rev. C 104, 024902 (2021), arXiv:2101.12413 [nucl-ex] .
  * Adamczewski-Musch  _et al._ (2020) J. Adamczewski-Musch _et al._ (HADES), Phys. Rev. C 102, 024914 (2020), arXiv:2002.08701 [nucl-ex] .
  * Acharya  _et al._ (2020) S. Acharya _et al._ (ALICE), Phys. Lett. B 807, 135564 (2020), arXiv:1910.14396 [nucl-ex] .
  * Vovchenko  _et al._ (2022) V. Vovchenko, V. Koch, and C. Shen, Phys. Rev. C 105, 014904 (2022), arXiv:2107.00163 [hep-ph] .
  * Koch (2010) V. Koch, “Hadronic Fluctuations and Correlations,” in  _Relativistic Heavy Ion Physics_, edited by R. Stock (2010) pp. 626–652, arXiv:0810.2520 [nucl-th] .
  * Vovchenko (2021) V. Vovchenko, in  _International Conference on Critical Point and Onset of Deconfinement_ (2021) arXiv:2110.02446 [nucl-th] .
  * Berdnikov and Rajagopal (2000) B. Berdnikov and K. Rajagopal, Phys. Rev. D 61, 105017 (2000), arXiv:hep-ph/9912274 .
  * Poberezhnyuk  _et al._ (2020) R. V. Poberezhnyuk, O. Savchuk, M. I. Gorenstein, V. Vovchenko, K. Taradiy, V. V. Begun, L. Satarov, J. Steinheimer, and H. Stoecker, Phys. Rev. C 102, 024908 (2020), arXiv:2004.14358 [hep-ph] .
  * Bzdak  _et al._ (2013) A. Bzdak, V. Koch, and V. Skokov, Phys. Rev. C 87, 014901 (2013), arXiv:1203.4529 [hep-ph] .
  * Vovchenko  _et al._ (2020a) V. Vovchenko, O. Savchuk, R. V. Poberezhnyuk, M. I. Gorenstein, and V. Koch, Phys. Lett. B 811, 135868 (2020a), arXiv:2003.13905 [hep-ph] .
  * Vovchenko  _et al._ (2020b) V. Vovchenko, R. V. Poberezhnyuk, and V. Koch, JHEP 10, 089 (2020b), arXiv:2007.03850 [hep-ph] .
  * Vovchenko (2022) V. Vovchenko, Phys. Rev. C 105, 014903 (2022), arXiv:2106.13775 [hep-ph] .
  * Barej and Bzdak (2022a) M. Barej and A. Bzdak, Phys. Rev. C 106, 024904 (2022a), arXiv:2205.05497 [hep-ph] .
  * Barej and Bzdak (2022b) M. Barej and A. Bzdak, (2022b), arXiv:2210.15394 [hep-ph] .
  * Braun-Munzinger  _et al._ (2017) P. Braun-Munzinger, A. Rustamov, and J. Stachel, Nucl. Phys. A 960, 114 (2017), arXiv:1612.00702 [nucl-th] .
  * Savchuk  _et al._ (2020) O. Savchuk, R. V. Poberezhnyuk, V. Vovchenko, and M. I. Gorenstein, Phys. Rev. C 101, 024917 (2020), arXiv:1911.03426 [hep-ph] .
  * Sorensen and Koch (2021) A. Sorensen and V. Koch, Phys. Rev. C 104, 034904 (2021), arXiv:2011.06635 [nucl-th] .
  * Savchuk  _et al._ (2023) O. Savchuk, R. V. Poberezhnyuk, A. Motornenko, J. Steinheimer, M. I. Gorenstein, and V. Vovchenko, Phys. Rev. C 107, 024913 (2023), arXiv:2211.13200 [hep-ph] .
  * Steinheimer and Randrup (2012) J. Steinheimer and J. Randrup, Phys. Rev. Lett. 109, 212301 (2012), arXiv:1209.2462 [nucl-th] .
  * Steinheimer and Randrup (2013) J. Steinheimer and J. Randrup, Phys. Rev. C 87, 054903 (2013), arXiv:1302.2956 [nucl-th] .
  * Molitoris  _et al._ (1984) J. J. Molitoris, J. B. Hoffer, H. Kruse, and H. Stoecker, Phys. Rev. Lett. 53, 899 (1984).
  * Peilert  _et al._ (1988) G. Peilert, A. Rosenhauer, H. Stoecker, W. Greiner, and J. Aichelin, Mod. Phys. Lett. A 3, 459 (1988).
  * Aichelin  _et al._ (1988) J. Aichelin, G. Peilert, A. Bohnet, A. Rosenhauer, H. Stöcker, and W. Greiner, Phys. Rev. C 37, 2451 (1988).
  * Peilert  _et al._ (1989) G. Peilert, H. Stöcker, W. Greiner, A. Rosenhauer, A. Bohnet, and J. Aichelin, Phys. Rev. C 39, 1402 (1989).
  * Hartnack  _et al._ (1998) C. Hartnack, R. K. Puri, J. Aichelin, J. Konopka, S. A. Bass, H. Stoecker, and W. Greiner, Eur. Phys. J. A 1, 151 (1998), arXiv:nucl-th/9811015 .
  * Pochodzalla  _et al._ (1995) J. Pochodzalla _et al._ , Phys. Rev. Lett. 75, 1040 (1995).
  * Natowitz  _et al._ (2002) J. B. Natowitz, K. Hagel, Y. Ma, M. Murray, L. Qin, R. Wada, and J. Wang, Phys. Rev. Lett. 89, 212701 (2002), arXiv:nucl-ex/0204015 .
  * Karnaukhov  _et al._ (2003) V. A. Karnaukhov _et al._ , Phys. Rev. C 67, 011601 (2003), arXiv:nucl-ex/0302006 .
  * Caillol (1998) J. Caillol, The Journal of chemical physics 109, 4885 (1998).
  * Vovchenko  _et al._ (2015a) V. Vovchenko, D. V. Anchishkin, M. I. Gorenstein, and R. V. Poberezhnyuk, Phys. Rev. C 92, 054901 (2015a), arXiv:1506.05763 [nucl-th] .
  * Hatta and Ikeda (2003) Y. Hatta and T. Ikeda, Phys. Rev. D 67, 014028 (2003), arXiv:hep-ph/0210284 .
  * Mukherjee  _et al._ (2017) A. Mukherjee, J. Steinheimer, and S. Schramm, Phys. Rev. C 96, 025205 (2017), arXiv:1611.10144 [nucl-th] .
  * Poberezhnyuk  _et al._ (2019) R. Poberezhnyuk, V. Vovchenko, A. Motornenko, M. I. Gorenstein, and H. Stoecker, Phys. Rev. C 100, 054904 (2019), arXiv:1906.01954 [hep-ph] .
  * Motornenko  _et al._ (2020) A. Motornenko, J. Steinheimer, V. Vovchenko, S. Schramm, and H. Stoecker, Phys. Rev. C 101, 034904 (2020), arXiv:1905.00866 [hep-ph] .
  * Stephan  _et al._ (2020) S. Stephan, J. Staubach, and H. Hasse, Fluid Phase Equilibria 523, 112772 (2020).
  * Evans and da Gama (1979) R. Evans and M. T. da Gama, Molecular Physics 38, 687 (1979), https://doi.org/10.1080/00268977900101971 .
  * Stephan  _et al._ (2019) S. Stephan, M. Thol, J. Vrabec, and H. Hasse, Journal of Chemical Information and Modeling 59, 4248 (2019).
  * Sanz  _et al._ (2013) E. Sanz, C. Vega, J. R. Espinosa, R. Caballero-Bernal, J. L. F. Abascal, and C. Valeriani, Journal of the American Chemical Society 135, 15008 (2013).
  * Frenkel (1946) J. Frenkel, “Kinetic theory of liquids,” Oxford: at the Clarendon Press, 488 p. (1946).
  * Fisher (1967) M. E. Fisher, Physics Physique Fizika 3, 255 (1967).
  * Stöcker and Greiner (1986) H. Stöcker and W. Greiner, Physics Reports 137, 277 (1986).
  * Lindenstruth (1993) V. Lindenstruth, GSI report 93-18 (1993).
  * Lindenstruth  _et al._ (1993) Lindenstruth _et al._ , Prepared for the US Department of Energy under Contract No. DE-AC03-76SF00098 , 53 (1993).
  * Pochodzalla (1997) J. Pochodzalla, Prog. Part. Nucl. Phys. 39, 443 (1997).
  * Klein and Unger (1983) W. Klein and C. Unger, Phys. Rev. B 28, 445 (1983).
  * Csernai and Neda (1994) L. P. Csernai and Z. Neda, Phys. Lett. B 337, 25 (1994).
  * Wedekind  _et al._ (2009) J. Wedekind, G. Chkonia, J. Wölk, R. Strey, and D. Reguera, The Journal of Chemical Physics 131, 114506 (2009), https://doi.org/10.1063/1.3204448 .
  * Pratt (1990) S. Pratt, Phys. Rev. A 42, 7447 (1990).
  * Kunde  _et al._ (1995) G. J. Kunde _et al._ , Phys. Rev. Lett. 74, 38 (1995).
  * Pratt  _et al._ (1995) S. Pratt, C. Montoya, and F. Ronning, Phys. Lett. B 349, 261 (1995).
  * Lopez and Luebeck (1989) J. Lopez and G. Luebeck, Physics Letters B 219, 215–221 (1989).
  * Favvas and Mitropoulos (2008) E. Favvas and A. Mitropoulos, Journal of Engineering Science and Technology Review 1 (2008), 10.25103/jestr.011.05.
  * Elliott (1989) C. M. Elliott, “The cahn-hilliard model for the kinetics of phase separation,” in  _Mathematical Models for Phase Change Problems_, edited by J. F. Rodrigues (Birkhäuser Basel, Basel, 1989) pp. 35–73.
  * Sator (2003) N. Sator, Physics Reports 376, 1 (2003).
  * Moretto  _et al._ (2011) L. G. Moretto, J. B. Elliott, L. Phair, and P. T. Lake, Journal of Physics G: Nuclear and Particle Physics 38, 113101 (2011).
  * Hill (1955) T. L. Hill, The Journal of Chemical Physics 23, 617 (1955), https://doi.org/10.1063/1.1742067 .
  * Poberezhnyuk  _et al._ (2021) R. V. Poberezhnyuk, O. Savchuk, M. I. Gorenstein, V. Vovchenko, and H. Stoecker, Phys. Rev. C 103, 024912 (2021), arXiv:2011.06420 [hep-ph] .
  * Vovchenko  _et al._ (2015b) V. Vovchenko, D. V. Anchishkin, and M. I. Gorenstein, J. Phys. A 48, 305001 (2015b), arXiv:1501.03785 [nucl-th] .
  * Satarov  _et al._ (2021) L. M. Satarov, R. V. Poberezhnyuk, I. N. Mishustin, and H. Stoecker, Phys. Rev. C 103, 024301 (2021), arXiv:2009.13487 [nucl-th] .
  * Bell (1927) E. T. Bell, Annals of Mathematics 29, 38 (1927).
  * Debye (1909) P. Debye, “Näherungsformeln für die Zylinderfunktionen für grosse Werte des Arguments und unbeschränkt veranderliche Werte des Index,”  (1909).