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

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 baryonic interactions through a density-dependent potential. This allows for an 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 [1,2]. 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 [3][4][5][6][7][8][9][10]. Phase transition should show itself in the collapse of the directed flow [11][12][13][14][15], fluctuations and correlations [2,[16][17][18][19], as well as light nuclei enhancement [20]. Another probe is electromagnetic radiation, in particular measurements of dilepton emission, is a promising tool for studying the properties of dense and hot matter [21][22][23][24][25].
During the freeze-out stage of A+A collisions the interactions between particles are expected to be weak. At this time, the multiplicities N i 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 [26][27][28]. 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 [29]. Therefore one should try to study event-by-event particle number fluctuations [2,30]. 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 , (∆N i ) 2 ≡ σ 2 , (∆N i ) 3 , (∆N i ) 4 , etc, where ∆N i ≡ N i − N i . In particular, the scaled variance ω, (normalized) skewness Sσ, and kurtosis κσ 2 of the particle number distribution are defined as follows: where κ n are the n-th cumulants of the N i -distribution. Such volume-independent (intensive) measures of number fluctuations (1)(2)(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 [31,32]. The equilibration time τ for fluctuations is typically larger than for mean quantities [18,19]. This is especially true for high order fluctuation measures, such as the kurtosis in the vicinity of the critical point (critical slowing down) [33,34]. 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 [35,36].
Previously, the memory effect was addressed within hydrodynamics with explicit evolution of fluctuations (see, e.g., Ref. [37] and references therein) or evolution equations for the cumulants [36]. Significant memory effects on the net proton cumulants were predicted if the sys-tem 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 freezeout 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
The UrQMD transport model [5][6][7] 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 [38,39]. 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. [40,41] for more details) 1 .
For the present work we use an EoS which was derived from the Chiral SU(3)-flavor parity-doublet Polyakovloop quark-hadron mean-field model (CMF) [44][45][46][47], both in its most recent version [48], 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 cut B = 2.6n 0 , where n 0 = 0.16 fm −3 is the nuclear matter saturation density. For n B > n cut B 1 Note, that the extension of our approach to a relativistic description is not a trivial task due to the fact that QMD is not a local mean field model. A possible solution to this problem has been presented in the RQMD approach [42,43]. However, the RQMD approach would also require the use of the scalar in addition the vector potential energy density and a method how these can be extracted from any EoS-model would have to be developed. This problem will be addressed in future works. the potential is then shifted by ∆n B = 2.6n 0 . The meanfield energy between n cut B < n B < n cut B + ∆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. [25,41,49] where this procedure was introduced.
The specific choice of the value for n cut is motivated by two factors: 1. The transition should be reachable with heavy ion collision experiments (which limits the density range to n cut B < 4 − 5n 0 ) [40]. 2. The density is higher than 2n 0 to avoid discrepancies with available constraints from heavy ion collisions and astrophysical observations [50]. 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 [51]. 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 [32]. 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.
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 lab = 1.23 ÷ 10A GeV.
In particular, we study the time dependence of the net baryon, B ≡ N B − NB, and net proton, p ≡ N p − Np, numbers inside a cubic volume of size V = 27 fm 3 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 lab below 10 AGeV one has NB N B , thus the antibaryons and antiprotons play essentially no role.
The size of the central volume is not unique and we have checked that moderate changes of the box size will only modify the extracted second order cumulants according to the fraction of the baryon number enclosed in the box α, i.e. by the factor 1 − α (see Refs. [52][53][54]) and therefore our conclusions remain independent on the box size. On the other hand a too large volume will suppress the fluctuations due to the conservation of the baryon number and a too small volume will simply give the Poissonian baseline result. As the range of the QMD interaction is determined as a few fm (given by the range parameter L in the Gaussian wave package) it is reasonable to choose a box length which is as big or larger than the interaction range but as small as possible to avoid effects from conservation, thus a box length of 3 fm was selected. 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 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 is extracted from particle and energy densities in the cell, matching to the CMF-EoS. The trajectories end at n B = n 0 . At this stage of the expansion one still has N B = n 0 V ∼ = 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 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.

III. EVENT-BY-EVENT FLUCTUATIONS
Baryon number fluctuations are sensitive probes of the interactions [8,9]. However, several additional factors affect these fluctuations, including baryon and elec-  tric charge conservation, non-equilibrium dynamics, resonance formation and decay, as well as finite size effects. The detailed microscopical UrQMD simulations naturally include these effects.
A. Coordinate space 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 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.
The enhancement of ω[B] due to the presence of the PT is most prominent at E 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.  understood as an additional acceptance effect modeled by a binomial distribution [35,55,56]. We find that ω[p] is a good qualitative proxy of ω[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.
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 fm 3 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 [57].
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 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 fm 3 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, ω[B] < 1, bell shaped curve. In the case of the PT EoS the P (B) distribution is much wider, ω[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 ≈ N B and N B ≈ 0 at E 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.

B. 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 [52][53][54]. 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 nonequilibrium 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 The fluctuation measures (1-3) for the final protons calculated in the UrQMD model at E lab = 2 AGeV for 3 million events are shown in Fig. 6 as functions of the rapidity acceptance interval ∆y. The statistical errors for the moments are estimated using the Delta method [58] using the open source package sample-moments [59].
The scaled variance ω[p] does not show large fluctuation effects. One finds ω[p] ∼ = 1 ÷ 1.04 at small ∆y ∼ = 0.2 and it decreases monotonously at larger values of ∆y. As a function of the rapidity acceptance interval ω[p] it is only weakly influenced by the presence of the PT. In addition, for the CMF EoS the values of ω[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 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 ∆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 ∆y. This behavior is considered to be the consequence of baryon conservation and excluded volume repulsion effects [60]. The qualitatively different behavior of increase of ω[p] with ∆y from unity up to large values of ω[p] ∼ = 2.3 at ∆y = 1 was reported by the HADES Collaboration in Au+Au collisions at E lab = 1.23A GeV [61].
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 [62]. 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 [63,64]. We hope to address the question of the role of nuclear clusters for proton number fluctuations at HADES energy region in future studies.