# Phase transition amplification of proton number fluctuations in nuclear
collisions
from a transport model approach

Oleh Savchuk  Facility for Rare Isotope Beams, Michigan State University, East
Lansing, MI 48824 USA  Bogolyubov Institute for Theoretical Physics, Kyiv,
Ukraine  GSI Helmholtzzentrum für Schwerionenforschung GmbH, Planckstr. 1,
D-64291 Darmstadt, Germany  Roman V. Poberezhnyuk  Bogolyubov Institute for
Theoretical Physics, Kyiv, Ukraine  Frankfurt Institute for Advanced Studies,
Giersch Science Center, Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main,
Germany  Anton Motornenko  Frankfurt Institute for Advanced Studies, Giersch
Science Center, Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main, Germany  Jan
Steinheimer  Frankfurt Institute for Advanced Studies, Giersch Science Center,
Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main, Germany  Mark I. Gorenstein
Bogolyubov Institute for Theoretical Physics, Kyiv, Ukraine  Frankfurt
Institute for Advanced Studies, Giersch Science Center, Ruth-Moufang-Str. 1,
D-60438 Frankfurt am Main, Germany  Volodymyr Vovchenko  Institute for Nuclear
Theory, University of Washington, Box 351550, Seattle, WA 98195, USA
Frankfurt Institute for Advanced Studies, Giersch Science Center, Ruth-
Moufang-Str. 1, D-60438 Frankfurt am Main, Germany

(August 27, 2024)

###### Abstract

The time evolution of particle number fluctuations in nuclear collisions at
intermediate energies ($E_{\rm lab}=1.23-10A$ GeV) is studied by means of the
UrQMD-3.5 transport model. The transport description incorporates an
implementation of baryonic interactions through a density-dependent potential.
This allows for a consistent implementation of a first order phase transition
including a mechanically unstable region at large baryon density. The scaled
variance of the baryon and proton number distributions is calculated in the
central cubic spatial volume of the collisions at different times. A
significant enhancement of fluctuations associated with the unstable region is
observed. This enhancement persists to late times reflecting a memory effect
for the fluctuations. The presence of the phase transition has a much smaller
influence on the observable event-by-event fluctuations of protons in momentum
space.

##  I Introduction

The structure of the QCD phase diagram is one of the major open questions in
high energy physics with great theoretical and experimental efforts on the
line Busza  _et al._ (2018); Bzdak  _et al._ (2020). Most experimental
information for these studies comes from observables in nucleus-nucleus (A+A)
collisions. The hypothetical first order QCD phase transition at finite baryon
density ending in a critical point can be probed in A+A collisions at
intermediate energies. The dynamical description of A+A collisions is based
usually on hydrodynamics and phenomenological transport models Lin  _et al._
(2005); Petersen  _et al._ (2019); Bass  _et al._ (1998); Bleicher  _et al._
(1999); Bleicher and Bratkovskaya (2022); Li and Ko (2017); Steinheimer  _et
al._ (2018); Ohnishi and Randrup (1995). Phase transition should show itself
in the collapse of the directed flow Rischke  _et al._ (1995); Stoecker
(2005); Brachmann  _et al._ (2000a, b); Oliinychenko  _et al._ (2022),
fluctuations and correlations Stephanov  _et al._ (1998, 1999); Jeon and Koch
(2000); Asakawa  _et al._ (2000); Bzdak  _et al._ (2020), as well as light
nuclei enhancement Sun  _et al._ (2017). Another probe is electromagnetic
radiation, in particular measurements of dilepton emission, is a promising
tool for studying the properties of dense and hot matter Shuryak (1978);
McLerran and Toimela (1985); Bratkovskaya and Cassing (1997); Rapp and Wambach
(2000); Savchuk  _et al._ (2022a).

During the freeze-out stage of A+A collisions the interactions between
particles are expected to be weak. At this time, the multiplicities $\langle
N_{i}\rangle$ of the $i^{th}$ different measurable hadron species in A+A
collisions appear to be reasonably well described within the framework of an
ideal hadron resonance gas model Adamczyk  _et al._ (2017); Harabasz  _et al._
(2020); Motornenko  _et al._ (2021). Thermal fits are not sensitive to a
possible phase transition during the early stages of evolution in A+A
reaction, the extracted S/A could be in indirect ways Bumnedpan  _et al._
(2022). Therefore one should try to study event-by-event particle number
fluctuations Luo and Xu (2017); Bzdak  _et al._ (2020). Here, the hope is to
understand the equation of state (EoS) properties from measurements of
intensive combinations of the central moments of the event-by-event
multiplicity distribution of the proton number $N_{i}$, $\langle(\Delta
N_{i})^{2}\rangle\equiv\sigma^{2}$, $\langle(\Delta N_{i})^{3}\rangle$,
$\langle(\Delta N_{i})^{4}\rangle$, etc, where $\Delta N_{i}\equiv
N_{i}-\langle N_{i}\rangle$. In particular, the scaled variance $\omega$,
(normalized) skewness $S\sigma$, and kurtosis $\kappa\sigma^{2}$ of the
particle number distribution are defined as follows:

$\displaystyle\omega[N_{i}]$ $\displaystyle=\frac{\sigma^{2}}{\langle
N_{i}\rangle}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{\kappa_{2}}{\kappa_{1}}\leavevmode\nobreak\ ,$ (1) $\displaystyle
S\sigma[N_{i}]$ $\displaystyle=\frac{\langle(\Delta
N_{i})^{3}\rangle}{\sigma^{2}}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{\kappa_{3}}{\kappa_{2}}\leavevmode\nobreak\ ,$ (2)
$\displaystyle\kappa\sigma^{2}[N_{i}]$ $\displaystyle=\frac{\langle(\Delta
N_{i})^{4}\rangle-3\langle(\Delta
N_{i})^{2}\rangle^{2}}{\sigma^{2}}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{\kappa_{4}}{\kappa_{2}}\leavevmode\nobreak\ ,$ (3)

where $\kappa_{n}$ are the $n$-th cumulants of the $N_{i}$-distribution. Such
volume-independent (intensive) measures of number fluctuations (1-3) can also
be applied to conserved charges such as net baryon number $B$ and electric
charge $Q$.

The amplification of baryon number fluctuations due to the first order phase
transition has been predicted in fluid-dynamical simulations of heavy-ion
collisions Steinheimer and Randrup (2012, 2013). The equilibration time $\tau$
for fluctuations is typically larger than for mean quantities Asakawa  _et
al._ (2000); Jeon and Koch (2000). This is especially true for high order
fluctuation measures, such as the kurtosis in the vicinity of the critical
point (critical slowing down) Berdnikov and Rajagopal (2000); Bluhm  _et al._
(2020). Thus, one can expect a memory effect for particle number fluctuations
during the fast expansion processes: i.e. large fluctuations generated in the
intermediate stages of A+A collisions partially survive to the freeze-out
stage Kitazawa and Asakawa (2012); Mukherjee  _et al._ (2015).

Previously, the memory effect was addressed within hydrodynamics with explicit
evolution of fluctuations (see, e.g., Ref. Pradeep  _et al._ (2022) and
references therein) or evolution equations for the cumulants Mukherjee  _et
al._ (2015). Significant memory effects on the net proton cumulants were
predicted if the system trajectory passes the vicinity of the critical point
at supercritical temperatures. The goal of the present paper is to explore the
enhancement of event-by-event fluctuations within a dynamical description of
A+A collisions and analyze whether it may survive to the freeze-out stage. The
advantage of this approach is that the dynamical evolution of particle number
fluctuations, in phase space, is treated consistently throughout the system’s
evolution.

The paper is organized as follows. In Sec. II the UrQMD model with the two
types of EoS is formulated. In Sec. III the results for baryon and proton
number fluctuations are presented, both in coordinate and momentum space. The
summary in Sec. IV closes the paper.

##  II  UrQMD simulations with a phase transition

Figure 1:  Pressure as a function of density at $T=0$ for the two scenarios
under investigation. The CMF EoS which has no phase transition is shown as the
orange line and the blue line depicts the construction with a phase
transition. The range in density where the pressure gradient is negative is
called the spinodal region (indicated by the red shaded area).

The UrQMD transport model Bass  _et al._ (1998); Bleicher  _et al._ (1999);
Bleicher and Bratkovskaya (2022) describes A+A collisions in terms of the
explicit propagation of hadrons in phase space, their elastic and inelastic
two-body reactions, and decay of unstable particles. In its cascade version
the effective equation of state of UrQMD resembles a non-interacting hadron
resonance gas Bravina  _et al._ (1998); Zabrodin  _et al._ (2009). To include
a more realistic EoS and even a phase transition, a density dependent
potential is introduced in the QMD part of the simulation. The interactions
among baryons are implemented via a density dependent potential energy per
baryon $V(n_{B})$. This allows to incorporate any density-dependent EoS in the
non-relativistic Hamilton equations of motion (see Refs. Omana Kuttan  _et
al._ (2022); Steinheimer  _et al._ (2022) for more details). For the present
work we use an EoS which was derived from the Chiral SU(3)-flavor parity-
doublet Polyakov-loop quark-hadron mean-field model (CMF) Steinheimer  _et
al._ (2011a, b); Mukherjee  _et al._ (2017); Motornenko  _et al._ (2019), both
in its most recent version Motornenko  _et al._ (2020), and in a modified
version that contains an additional first order phase transition (PT), denoted
further as PT. The CMF model incorporates a realistic description of nuclear
matter with nuclear incompressibility of $K_{0}=267$ MeV, chiral symmetry
breaking in the hadronic and quark sectors, as well as an effective
deconfinement transition.

In the case of the PT model a simple augmentation is used in order to
implement an additional phase transition. To provide a metastable state at
large baryon densities the original mean-field potential of the CMF model is
truncated at density $n_{B}^{cut}=2.6n_{0}$, where $n_{0}=0.16$ fm-3 is the
nuclear matter saturation density. For $n_{B}>n_{B}^{cut}$ the potential is
then shifted by $\Delta n_{B}=2.6n_{0}$. The mean-field energy between
$n_{B}^{cut}<n_{B}<n_{B}^{cut}+\Delta n_{B}$ is interpolated by a third order
polynomial in order to create a second minimum in the energy per baryon
$V(n_{B})$ and to ensure that its derivative is a continuous function. Note,
that this procedure modifies the CMF EoS only at high baryon densities,
leaving low-density description consistent with nuclear matter properties and
lattice QCD constraints. For details see Refs. Savchuk  _et al._ (2022a); Li
_et al._ (2022); Steinheimer  _et al._ (2022) where this procedure was
introduced. Figure 1 compares the effective pressure of these two equations of
state at $T=0$ as a function of density. In this figure also the region of
mechanical instabilities is highlighted. If the system created in a collision
enters this area, the interactions will drive the system to rapidly separate
into the two coexisting phases, a mechanism well known as spinodal
decomposition. This phase separation will then lead to the formation of
clusters and consequently to an enhancement of the baryon fluctuations in
coordinate space. At this point it is important to note that the process of
phase separation and the properties (like the surface energy) of the clusters
which are created depends strongly on the finite range interactions present in
the system Randrup (2004). In the QMD approach used in this work, the
effective range of the interactions is governed by the width of the Gaussian
wave packages which are used to calculate the interaction density. While this
parameter was fixed to give a proper description of nuclear matter and nuclei
it can be argued that for another phase transition (quark-hadron for example)
this range may be different. We therefore expect that our results on the
properties of the baryon clumps created may not be very precise and only can
give a qualitative picture of the dynamics that are expected. Fortunately it
has been shown in a previous work, where the role of the finite range
interactions on the spinodal clumping in nuclear collisions was studied, it
was found that the quantitative effect on the density fluctuation is indeed
rather small Steinheimer and Randrup (2013). Thus, the following results may
be parameter dependent on a quantitative level, the general results of our
work are likely quite robust as observed in previous studies.

Figure 2:  The trajectory of central Au+Au collisions on the phase diagram in
the plane of baryon density $n_{B}$ and temperature $T$ at collision energies
$E_{\rm lab}=1.23,\,2,\,4,\,6,$ and $10A$ GeV (from bottom to top). The
trajectories correspond to the central cubic box of volume $V=27$ fm3. The
blue and orange lines represent the PT and CMF equations of state,
respectively. The diamonds and circles correspond to densities
$n_{B}=3.2\leavevmode\nobreak\ n_{0}$ and $n_{B}=1.2\leavevmode\nobreak\
n_{0}$, respectively. Red symbols correspond to times of $17fm/c$ and $20fm/c$
and black ones to $10fm/c$ and $16fm/c$. The red shaded area shows the
mechanically unstable region bounded by the spinodals for the PT EoS.

To study the time evolution of the bulk matter as well as the fluctuations
caused by the phase transition, we simulate $50000$ collision events for very
central ($b<2$ fm or $2\%$ most central) Au-Au collisions in the SIS18/SIS100
and RHIC-BES energy range of $E_{\rm lab}=1.23\div 10A$ GeV.

In particular, we study the time dependence of the net baryon, $B\equiv
N_{B}-N_{\bar{B}}$, and net proton, ${\rm p}\equiv N_{p}-N_{\bar{p}}$, numbers
inside a cubic volume of size $V=27$ fm3 located in the geometrical center in
the center of mass frame of Au+Au reaction. Note, that in the considered
region of collision energies $E_{\rm lab}$ below 10 AGeV one has
$N_{\bar{B}}\ll N_{B}$, thus the antibaryons and antiprotons play essentially
no role.

Figure 2 presents the UrQMD results for the event averaged expansion
trajectories of the central cubic volume in the $(n_{B},T)$-plane for
different beam energies. These trajectories are started at the maximal $n_{B}$
values achieved at given $E_{\rm lab}$. At this stage the momentum spectra of
baryons inside the central volume $V$ are consistent with the thermal Maxwell-
Boltzmann distribution. The temperature of the central volume can be extracted
from particle and energy densities in the cell. The trajectories end at
$n_{B}=n_{0}$. At this stage of the expansion one still has $N_{B}=n_{0}V\cong
4$ baryons inside the central volume.

In case of the PT model, the effects of the softening of the EoS occur in the
spinodal region. This results in a significant increase of the compression and
large baryon densities.

Also, the time evolution of the density is strongly affected by the EoS and
the times at which the system, created at $E_{\rm lab}=2A$ GeV, reaches $1.2$
and $3.2$ nuclear saturation densities is pointed out with the red and black
symbols. During its evolution the created system with a PT evolves through
mixed phase region which leads to a significantly longer expansion time and
possibly to an increase of the fluctuations in coordinate space. The purpose
of this work is to test to which extent these large fluctuations in the
coordinate space may survive to the late stages of the collisions and whether
they can be observed in a momentum space acceptance similar to experimental
detection.

Figure 3: The UrQMD results in central Au+Au collisions for the scaled
variance of net baryon (upper panel) and net proton (lower panel)
distributions inside a central cubic box of volume $V=27fm^{3}$ as a function
of the collision time in the center of mass system. The blue and orange lines
correspond to PT and CMF EoS, respectively. Dotted horizontal lines show the
Poisson baseline of $\omega=1$. The time $t=0$ is taken as the moment of
maximal baryon density in the central box. The lines are stopped at $t$ values
that correspond to $n_{B}=n_{0}$ in the central cell. The symbols for $E_{\rm
lab}=2A$ GeV have the same meaning as in Fig. 2. For each collision energy the
unstable spinodal region is shown by the shaded area.

##  III Event-by-Event Fluctuations

Baryon number fluctuations are sensitive probes of the interactions Li and Ko
(2017); Steinheimer  _et al._ (2018). However, several additional factors
affect these fluctuations, including baryon and electric charge conservation,
non-equilibrium dynamics, resonance formation and decay, as well as finite
size effects. The detailed microscopical UrQMD simulations naturally include
these effects.

Due to these effects, the optimal choice of the phase-space region, suitable
for observing the event-by-event fluctuations, is an important aspect. The
acceptance should be large enough compared to the correlation length relevant
for the important physical processes. However, if the acceptance is too large
the global conservation laws dominate Vovchenko  _et al._ (2020a);
Poberezhnyuk  _et al._ (2020); Vovchenko  _et al._ (2020b), i.e., by detecting
all final particles one does not permit any event-by-event fluctuations of the
conserved charges

###  III.1 Coordinate space

Figure 3 (upper panel) presents the UrQMD results for the scaled variance of
baryon number fluctuations, $\omega[B]$, in the central volume $V=3^{3}$ fm3
as functions of time for different collision energies. The results are
compared to the Poisson distribution baseline, $\omega=1$, which corresponds
to the production of uncorrelated particles and usually is approached late
when essentially no particles are present in the central volume. Deviations
from the Poisson baseline are evident at intermediate times. In the CMF case
one has $\omega<1$ which is driven by the global charge conservation effects.
However, in the case of PT, one observes $\omega>1$, and the maximum values
are reached when the system crosses the PT region. This behavior of
fluctuations reflects the enhancement of fluctuations due to the instabilities
associated with the first order PT.

The two EoS are indistinguishable at lower baryon densities, in particular at
$n=n_{0}$, reached at later stages of the collision. The scaled variance of
baryon number distribution at time moments corresponding to $n=n_{0}$ inside
the box is still different in the two scenarios despite the fact that the two
EoS become identical at this stage. We find that the scaled variance is
enhanced by a factor of 2 at $E_{\rm lab}=2A\,$GeV at $2n_{0}$ in the PT
scenario. Thus, in coordinate space one observes a memory effect for the
fluctuations in the fast expanding system.

Figure 4:  Baryon number densities in five different events in the transverse
coordinate plane $z=0$ at $t=17\leavevmode\nobreak\ $fm/$c$ in Au+Au central
collisions at $E_{\rm lab}=2A$ GeV for the PT EoS. The time moment of $t=17$
fm/$c$ corresponds to average $n_{B}=3.2\,n_{0}$ inside the central box when
the matter is in the mixed liquid-gas state. The projections of the central
box on the $z=0$ plane are shown by the white squares.

The enhancement of $\omega[B]$ due to the presence of the PT is most prominent
at $E_{\rm lab}=2A\,$GeV. This enhancement decreases rapidly with further
increase of collision energy. This is because the system spends less time
inside the PT region during the expansion at the higher collision energies.
The shorter time intervals become insufficient to generate large fluctuations.

In the experiment, it is problematic to measure all baryons due to
difficulties with the detection of neutral particles such as neutrons. For
this reason, fluctuations of the proton number are commonly used as a proxy
for the baryon number. The lower panel of Fig. 3 depicts the UrQMD results for
the scaled variance of the proton number distribution, $\omega[{\rm p}]$. This
quantity retains the qualitative features observed for $\omega[B]$ but shows
smaller deviations from the Poisson baseline of unity. This can be understood
as an additional acceptance effect modeled by a binomial distribution Kitazawa
and Asakawa (2012); Bzdak and Koch (2012); Savchuk  _et al._ (2020). We find
that $\omega[{\rm p}]$ is a good qualitative proxy of $\omega[B]$, while
quantitatively the fluctuation signals are suppressed by about 50%. Note, that
at early times a significant part of baryons is present as resonance states
which are not included in the proton number, reducing the fraction of protons
among all baryons.

Figure 5:  Baryon number distributions $P(B)$ inside the central box of volume
$V=27\leavevmode\nobreak\ $fm3 (upper panels) and inside $|y|<0.05$ rapidity
interval (lower) panels. Left panels correspond to the dense $n_{B}=3.2n_{0}$
system that for PT EoS is inside the unstable region bounded by spinodals.
Right panels correspond to the systems at $n_{B}=1.2n_{0}$ where both CMF and
PT EoS have the same properties. The fits of $P(B)$ for the PT EoS with two
Gaussian distributions – dotted and dashed lines – are shown in the upper
panels by the red solid lines.  Figure 6:  The UrQMD results for the scaled
variance $\omega[{\rm p}]$, skewness $S\sigma[{\rm p}]$, and kurtosis
$\kappa\sigma^{2}[{\rm p}]$ of proton number distribution as functions of the
rapidity acceptance interval $\Delta y$ for Au-Au collisions at
$E_{lab}=2A\,$GeV.

To understand the enhancement of fluctuations due to the PT let us consider
its possible dynamics during an expansion. Two different scenarios are
possible. In case of a slow change in thermodynamic properties, i.e., if the
system is constantly in equilibrium, one expects to see nucleation and growth
of the new phase. However, nucleation is usually suppressed by a potential
barrier as a function of spatial coordinate that needs to be crossed in order
to achieve a stable droplet of the new phase. This process generally requires
a long time which exceeds the time scales reached in heavy-ion collisions.
Another possibility is spinodal decomposition. When a system is inside the
spinodal region on the phase diagram, the associated mechanical instability
makes fast separation into phases possible. This separation creates regions of
high and low density that correspond to two stable states of matter and, on
even-by-event basis, can create high fluctuations in the coordinate space. The
central box of volume $V=27\,$fm3 in the coordinate space can, in each event,
contain either the gas or the liquid phase, leading to two peaks in the
particle number distribution as separate contributions from the two phases
Sorensen and Koch (2021).

To point out this effect, figure 4 shows five different, single event, baryon
density distributions obtained in UrQMD simulations with the PT EoS at $E_{\rm
lab}=2A\,$GeV in the transverse $(x,y)$ coordinate plane at $z=0\,$fm. The
time was chosen as $t=17$ fm/$c$ when the central volume $V=27\,$fm3 is inside
the unstable region. The projection of the volume $V$ on the $(x,y)$-plane is
shown in Fig. 4 by the white squares. One observes that the volume $V$ with
average baryon density $n_{B}=3.2\,n_{0}$ contains either a high (liquid) or
low (gas) baryon density, or clump vs. no-clump scenario. These different
states inside the central box change from one event to another, leading to
large event-by-event fluctuations of the baryon density. In the case of the
CMF EoS no large changes of the baryon density inside the central box are
observed.

The UrQMD results for the baryon number distributions in the central box at
$E_{lab}=2A\,$GeV are shown in Fig. 5 (a) and (b), corresponding to
$n_{B}=3.2\,n_{0}$ and $n_{B}=1.2\,n_{0}$, respectively. The probability to
observe $B$ baryons $P(B)$ for the CMF equation of state can be described with
a narrow, $\omega[B]<1$, bell shaped curve. In the case of the PT EoS the
$P(B)$ distribution is much wider, $\omega[B]>1,$ and asymmetric around the
mean. In order to emphasize the bimodal nature of this distribution we fit it
by a sum of two Gaussian distributions. Note that $B\approx N_{B}$ and
$N_{\overline{B}}\approx 0$ at $E_{\rm lab}=2A\,$GeV. At the
$n_{B}=1.2\,n_{0}$ a noticeable difference in the baryon number distributions
is still seen due to the memory effect of the dynamical expansion through the
spinodal region for the PT EoS.

###  III.2 Momentum space

The experimental observations are limited to momenta of produced hadrons,
therefore, any information about spatial correlations has to be extracted from
measurements performed in momentum space. This is feasible in the presence of
strong space-momentum correlations due to collective flow, in particular, due
to Bjorken longitudinal flow at high collision energies, see the recently
developed subensemble acceptance method Vovchenko  _et al._ (2020a);
Poberezhnyuk  _et al._ (2020); Vovchenko  _et al._ (2020b). The method is not
applicable in direct proximity to the critical point where correlation length
reaches the size of the system and inside a mixed phase region where one can
expect coexistence of two phases and strong non-equilibrium effects. At lower
collision energies, such as those considered in the present work, the
collective flow is, however, rather weaker. Thus, new ways to extract signal
of fluctuation enhancement are required. Here, we compute the baryon number
distribution in the momentum space by imposing rapidity cuts around
midrapidity in the center of mass system of Au+Au. To understand the effect of
the existing coordinate space fluctuations on momentum space observations,
figures 5 (c) and (d) present the $p_{t}$-integrated baryon number
distributions inside a $|y|<0.05$ rapidity interval. The times shown
correspond to those already presented in Figs. 5 (a) and (b) in the coordinate
space. In contrast to the enhanced fluctuations in the coordinate space due to
the PT, no qualitative difference between the $P(B)$ distributions with or
without the PT is observed. In both the CMF and PT cases the distributions are
single-peaked and close to the Poisson distribution. The suppression effects
due to global charge conservation lead to $\omega[B]$ being slightly smaller
than unity.

The fluctuation measures (1-3) for the final protons calculated in the UrQMD
model at $E_{\rm lab}=2\,A$GeV for $3\leavevmode\nobreak\ million$ events are
shown in Fig. 6 as functions of the rapidity acceptance interval $\Delta y$.
The statistical errors for the moments are estimated using the Delta method
Kendall and Stuart (1963) using the open source package sample-moments
Vovchenko (2021).

The scaled variance $\omega[{\rm p}]$ does not show large fluctuation effects.
One finds $\omega[{\rm p}]\cong 1\div 1.04$ at small $\Delta y\cong 0.2$ and
it decreases monotonously at larger values of $\Delta y$. As a function of the
rapidity acceptance interval $\omega[{\rm p}]$ it is only weakly influenced by
the presence of the PT. In addition, for the CMF EoS the values of
$\omega[{\rm p}]$ are even slightly larger than those for the PT case. This
surprising result can be understood as follows: In the collisions with the
CMF-EoS a larger pressure gradient in the initial overlap zone is created
which starts to create a more violent sidewards push of the participants into
the spectator region. These participants then have a higher chance to
rescatter with the remaining spectators. This increases the average number of
participants in the CMF scenario as compared to the PT case. The enhancement
of fluctuations in the CMF case is a result of the increased volume
fluctuations which is a dynamical effect and is only indirectly related to the
EoS. Furthermore, it can be observed that the cumulant ratios in both
scenarios drop below the binomial baseline at large rapidities. This
interesting observation is due to the rapidity dependence of iso-spin
randomization (or pion production), where protons with large rapidity are more
likely to stem from elastic scatterings. Since elastic scatterings cannot
change the isospin of the proton, the effect of isospin randomization is
significantly suppressed at higher rapidities and thus the conservation effect
dominates.

The results for the skewness and kurtosis shown in Figs. 6 (b) and (c) are
qualitatively very similar to the results for the scaled variance. The large
negative values of the kurtosis again result from the large volume
fluctuations with a strong contribution of conservation laws. This indicates
that the enhancement of fluctuations seen in the coordinate space is almost
completely washed out when one turns to the momentum space. Momentum space
fluctuations are dominated by the interplay between volume fluctuations and
conservation laws which can act in a non-trivial way as a function of rapidity
at these low beam energies.

##  IV Conclusions

The baryon and proton number fluctuations in a presence of the first order
phase transition (PT) were studied. We used the UrQMD-3.5 model with a
density-dependent mean field interactions to simulate Au+Au collisions at
intermediate energies. Two different equations of state, the CMF and phase
transition augmented CMF, were used. By construction, the two EoS have the
same properties at $n_{B}<2n_{0}$. However, they differ at higher baryon
densities due to the presence of the PT. In both scenarios, heavy-ion
collisions follow very similar trajectories in the $(n_{B},T)$ plane, the main
difference being as expected a longer expansion time in the presence of the
PT.

Despite similar trajectories on the phase diagram in the two scenarios, we
observe notable differences in baryon number fluctuations in the central
coordinate space volume. The fluctuations exhibit an enhancement in the PT
case, and it persists for some later stages of the collision when the system
has already left the spinodal region. The largest effect is observed at the
collision energy of $E_{\rm lab}=2A\,$GeV which is purely related to the
location of the unstable phase in the equation of state used.

On the other hand, the enhancement of proton number fluctuations due to the PT
is not seen when these fluctuations are analyzed in momentum space. The scaled
variance, skewness, and kurtosis of final state protons as functions of the
rapidity acceptance interval $\Delta y$ are not sensitive to the PT.

At large RHIC and LHC collision energies the A+A data demonstrate a monotonous
decrease of proton number fluctuations with $\Delta y$. This behavior is
considered to be the consequence of baryon conservation and excluded volume
repulsion effects Vovchenko (2022). The qualitatively different behavior of
increase of $\omega[{\rm p}]$ with $\Delta y$ from unity up to large values of
$\omega[{\rm p}]\cong 2.3$ at $\Delta y=1$ was reported by the HADES
Collaboration in Au+Au collisions at $E_{\rm lab}=1.23A$ GeV Adamczewski-Musch
_et al._ (2020).

The present study shows that these large baryon number fluctuations are
unlikely to be the result of spinodal amplification due to a phase transition.
An alternative explanation of the HADES data was suggested in the recent paper
Savchuk  _et al._ (2022b). To search for interesting physics in the EoS, one
should first exclude the event-by-event fluctuations of nucleon participants.
At small collision energy there is a significant number of light nuclei and
intermediate nuclear fragments in the final state. This can influence
fluctuations of the number of bare protons. Our UrQMD simulations do not yet
include light nuclei and nuclear fragment production which could be relevant
at intermediate collision energies Bratkovskaya  _et al._ (2022); Kireyeu  _et
al._ (2022). We hope to address the question of the role of nuclear clusters
for proton number fluctuations at HADES energy region in future studies.

###### Acknowledgements.

The authors thank S. Pratt, M. Gazdzicki, V. Kuznietsov, A. Rustamov, and H.
Stoecker 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. A.M. acknowledges the Stern–Gerlach
Postdoctoral fellowship of the Stiftung Polytechnische Gesellschaft. V.V. was
supported through the U.S. Department of Energy, Office of Science, Office of
Nuclear Physics, under contract number DE-FG02-00ER41132.

## References

  * Busza  _et al._ (2018) W. Busza, K. Rajagopal, and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * 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] .
  * Lin  _et al._ (2005) Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, Phys. Rev. C 72, 064901 (2005).
  * Petersen  _et al._ (2019) H. Petersen, D. Oliinychenko, M. Mayer, J. Staudenmaier, and S. Ryu, Nuclear Physics A 982, 399 (2019).
  * Bass  _et al._ (1998) S. A. Bass _et al._ , Prog. Part. Nucl. Phys. 41, 255 (1998), arXiv:nucl-th/9803035 .
  * Bleicher  _et al._ (1999) M. Bleicher _et al._ , J. Phys. G 25, 1859 (1999), arXiv:hep-ph/9909407 .
  * Bleicher and Bratkovskaya (2022) M. Bleicher and E. Bratkovskaya, Prog. Part. Nucl. Phys. 122, 103920 (2022).
  * Li and Ko (2017) F. Li and C. M. Ko, Phys. Rev. C 95, 055203 (2017), arXiv:1606.05012 [nucl-th] .
  * Steinheimer  _et al._ (2018) J. Steinheimer, Y. Wang, A. Mukherjee, Y. Ye, C. Guo, Q. Li, and H. Stoecker, Phys. Lett. B 785, 40 (2018), arXiv:1804.08936 [nucl-th] .
  * Ohnishi and Randrup (1995) A. Ohnishi and J. Randrup, Phys. Rev. Lett. 75, 596 (1995).
  * Rischke  _et al._ (1995) D. H. Rischke, Y. Pürsün, J. A. Maruhn, H. Stoecker, and W. Greiner, Acta Phys. Hung. A 1, 309 (1995), arXiv:nucl-th/9505014 .
  * Stoecker (2005) H. Stoecker, Nucl. Phys. A 750, 121 (2005), arXiv:nucl-th/0406018 .
  * Brachmann  _et al._ (2000a) J. Brachmann, S. Soff, A. Dumitru, H. Stoecker, J. A. Maruhn, W. Greiner, L. V. Bravina, and D. H. Rischke, Phys. Rev. C 61, 024909 (2000a), arXiv:nucl-th/9908010 .
  * Brachmann  _et al._ (2000b) J. Brachmann, A. Dumitru, H. Stoecker, and W. Greiner, Eur. Phys. J. A 8, 549 (2000b), arXiv:nucl-th/9912014 .
  * Oliinychenko  _et al._ (2022) D. Oliinychenko, A. Sorensen, V. Koch, and L. McLerran, (2022), arXiv:2208.11996 [nucl-th] .
  * Stephanov  _et al._ (1998) M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. Lett. 81, 4816 (1998), arXiv:hep-ph/9806219 [hep-ph] .
  * Stephanov  _et al._ (1999) M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. D60, 114028 (1999), arXiv:hep-ph/9903292 [hep-ph] .
  * Jeon and Koch (2000) S. Jeon and V. Koch, Phys. Rev. Lett. 85, 2076 (2000), arXiv:hep-ph/0003168 [hep-ph] .
  * Asakawa  _et al._ (2000) M. Asakawa, U. W. Heinz, and B. Muller, Phys. Rev. Lett. 85, 2072 (2000), arXiv:hep-ph/0003169 [hep-ph] .
  * Sun  _et al._ (2017) K.-J. Sun, L.-W. Chen, C. M. Ko, and Z. Xu, Phys. Lett. B 774, 103 (2017), arXiv:1702.07620 [nucl-th] .
  * Shuryak (1978) E. V. Shuryak, Phys. Lett. B 78, 150 (1978).
  * McLerran and Toimela (1985) L. D. McLerran and T. Toimela, Phys. Rev. D 31, 545 (1985).
  * Bratkovskaya and Cassing (1997) E. L. Bratkovskaya and W. Cassing, Nucl. Phys. A 619, 413 (1997), arXiv:nucl-th/9611042 .
  * Rapp and Wambach (2000) R. Rapp and J. Wambach, Adv. Nucl. Phys. 25, 1 (2000), arXiv:hep-ph/9909229 .
  * Savchuk  _et al._ (2022a) O. Savchuk, A. Motornenko, J. Steinheimer, V. Vovchenko, M. Bleicher, M. Gorenstein, and T. Galatyuk, (2022a), arXiv:2209.05267 [nucl-th] .
  * Adamczyk  _et al._ (2017) L. Adamczyk _et al._ (STAR), Phys. Rev. C 96, 044904 (2017), arXiv:1701.07065 [nucl-ex] .
  * Harabasz  _et al._ (2020) S. Harabasz, W. Florkowski, T. Galatyuk, t. Ma Lgorzata Gumberidze, R. Ryblewski, P. Salabura, and J. Stroth, Phys. Rev. C 102, 054903 (2020), arXiv:2003.12992 [nucl-th] .
  * Motornenko  _et al._ (2021) A. Motornenko, J. Steinheimer, V. Vovchenko, R. Stock, and H. Stoecker, Phys. Lett. B 822, 136703 (2021), arXiv:2104.06036 [hep-ph] .
  * Bumnedpan  _et al._ (2022) T. Bumnedpan, J. Steinheimer, M. Bleicher, A. Limphirat, and C. Herold, Phys. Lett. B 835, 137537 (2022), arXiv:2209.04096 [nucl-th] .
  * Luo and Xu (2017) X. Luo and N. Xu, Nucl. Sci. Tech. 28, 112 (2017), arXiv:1701.02105 [nucl-ex] .
  * 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] .
  * Berdnikov and Rajagopal (2000) B. Berdnikov and K. Rajagopal, Physical Review D 61 (2000), 10.1103/physrevd.61.105017.
  * Bluhm  _et al._ (2020) M. Bluhm _et al._ , Nucl. Phys. A 1003, 122016 (2020), arXiv:2001.08831 [nucl-th] .
  * Kitazawa and Asakawa (2012) M. Kitazawa and M. Asakawa, Phys. Rev. C 86, 024904 (2012), [Erratum: Phys.Rev.C 86, 069902 (2012)], arXiv:1205.3292 [nucl-th] .
  * Mukherjee  _et al._ (2015) S. Mukherjee, R. Venugopalan, and Y. Yin, Phys. Rev. C 92, 034912 (2015), arXiv:1506.00645 [hep-ph] .
  * Pradeep  _et al._ (2022) M. Pradeep, K. Rajagopal, M. Stephanov, and Y. Yin, Phys. Rev. D 106, 036017 (2022), arXiv:2204.00639 [hep-ph] .
  * Bravina  _et al._ (1998) L. V. Bravina _et al._ , Phys. Lett. B 434, 379 (1998), arXiv:nucl-th/9804008 .
  * Zabrodin  _et al._ (2009) E. E. Zabrodin _et al._ , J. Phys. G 36, 064065 (2009), arXiv:0902.4601 [hep-ph] .
  * Omana Kuttan  _et al._ (2022) M. Omana Kuttan, A. Motornenko, J. Steinheimer, H. Stoecker, Y. Nara, and M. Bleicher, Eur. Phys. J. C 82, 427 (2022), arXiv:2201.01622 [nucl-th] .
  * Steinheimer  _et al._ (2022) J. Steinheimer, A. Motornenko, A. Sorensen, Y. Nara, V. Koch, and M. Bleicher, (2022), arXiv:2208.12091 [nucl-th] .
  * Steinheimer  _et al._ (2011a) J. Steinheimer, S. Schramm, and H. Stocker, J. Phys. G 38, 035001 (2011a), arXiv:1009.5239 [hep-ph] .
  * Steinheimer  _et al._ (2011b) J. Steinheimer, S. Schramm, and H. Stocker, Phys. Rev. C 84, 045208 (2011b), arXiv:1108.2596 [hep-ph] .
  * Mukherjee  _et al._ (2017) A. Mukherjee, J. Steinheimer, and S. Schramm, Phys. Rev. C 96, 025205 (2017), arXiv:1611.10144 [nucl-th] .
  * Motornenko  _et al._ (2019) A. Motornenko, V. Vovchenko, J. Steinheimer, S. Schramm, and H. Stoecker, Nucl. Phys. A 982, 891 (2019), arXiv:1809.02000 [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] .
  * Li  _et al._ (2022) P. Li, J. Steinheimer, T. Reichert, A. Kittiratpattana, M. Bleicher, and Q. Li, (2022), arXiv:2209.01413 [nucl-th] .
  * Randrup (2004) J. Randrup, Phys. Rev. Lett. 92, 122301 (2004), arXiv:hep-ph/0308271 .
  * 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] .
  * 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] .
  * Vovchenko  _et al._ (2020b) V. Vovchenko, R. V. Poberezhnyuk, and V. Koch, JHEP 10, 089 (2020b), arXiv:2007.03850 [hep-ph] .
  * Bzdak and Koch (2012) A. Bzdak and V. Koch, Phys. Rev. C 86, 044904 (2012), arXiv:1206.4286 [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] .
  * Kendall and Stuart (1963) M. Kendall and A. Stuart,  _The Advanced Theory of Statistics: Distribution theory_, The Advanced Theory of Statistics (Hafner Publishing Company, 1963).
  * Vovchenko (2021) V. Vovchenko, “A header-only C++ library for mean and standard error estimation of moments and cumulants,” (2021), https://github.com/vlvovch/sample-moments [Online; accessed 18-October-2022].
  * Vovchenko (2022) V. Vovchenko, (2022), arXiv:2208.13693 [hep-ph] .
  * Adamczewski-Musch  _et al._ (2020) J. Adamczewski-Musch _et al._ (HADES), Phys. Rev. C 102, 024914 (2020), arXiv:2002.08701 [nucl-ex] .
  * Savchuk  _et al._ (2022b) O. Savchuk, R. V. Poberezhnyuk, and M. I. Gorenstein, Physics Letters B 835, 137540 (2022b).
  * Bratkovskaya  _et al._ (2022) E. Bratkovskaya, S. Glässel, V. Kireyeu, J. Aichelin, M. Bleicher, C. Blume, G. Coci, V. Kolesnikov, J. Steinheimer, and V. Voronyuk, “Midrapidity cluster formation in heavy-ion collisions,”  (2022).
  * Kireyeu  _et al._ (2022) V. Kireyeu, J. Steinheimer, J. Aichelin, M. Bleicher, and E. Bratkovskaya, “Deuteron production in ultra-relativistic heavy-ion collisions: A comparison of the coalescence and the minimum spanning tree procedure,”  (2022).