Baryon clustering at the critical line and near the hypothetical critical point in heavy-ion collisions

We study clustering of baryons at the freeze-out point of relativistic heavy-ion collisions. Using a Walecka-Serot model for the nucleon-nucleon (NN) interaction we analyze how the modified/critical $\sigma$ mode---responsible for the NN attraction---allows for clustering of nucleons when the system is close to a possible critical point of QCD. We investigate clusters of few nucleons, and also the internal cluster configuration when the system is long lived. For realistic heavy-ion collisions we study to how extend such clusters can be formed in a finite time, and perform the statistical analysis of cumulants and higher-order moments (skewness and kurtosis) for collisions at the Beam Energy Scan of RHIC.

We start by contrasting known facts about high-and low-energy heavy-ion collisions, after which we will define the phenomena to be discussed in this work.
By high-energy collisions we mean those at Relativistic Heavy-Ion Collider (RHIC) full energy √ s N N ≈ 200 GeV and at the Large Hadron Collider (LHC) √ s N N ≈ 2 − 8 TeV. In these cases the particle yields are very accurately described by the so-called "resonance gas model", assuming that all interactions between hadrons can be effectively treated as an ideal gas of all known resonances, as suggested by the Beth-Uhlenbeck formula. As shown e.g. in Ref. [1], the yields per degree of freedom are on the same thermal exponent, from the lightest speciespions, kaons, etc-up to the baryons, hyperons and their antiparticles, all the way up to light nuclei such as 4 He. This trend is observed to hold over about 9 decades. The lesson is that, at chemical freeze-out temperatures T ch 150 MeV and at near-zero baryonic chemical potential µ B 0 the fireball is very well thermally equilibrated. Such high degree of equilibration undoubtedly is related with kinetic properties of the strongly coupled quark-gluon plasma (QGP), preceding the freeze-out stage of these collisions.
At low non-relativistic collision energy, √ s N N < 1 GeV, creating nuclear matter with the temperature T 10 MeV, one observes the so-called "multifragmentation" phenomena, production of large variety of nuclear fragments, with a wide power-like distributions. It is attributed to a nearby presence of a critical point, separating liquid nuclear matter from a gas-like phase. For a review see e.g. [2]. The production of various nuclear clusters is not in equilibrium, and is very sensitive to the relation between the temperature and time available for cluster formation. This regime can be compared to that in atomic physics studying various out-of-equilibrium situations, for example "snow production" machines, operating in between water and ice phases. Heavy-ion collisions at intermediate energies have been studied in 1980s, both at CERN Super Proton Synchrotron and the Brookhaven National Laboratory Alternating Gradient Synchrotron, but not in sufficient detail. Many models predict that baryon-rich matter will also have the first order transition line, ending in certain critical point. Its search using enhanced fluctuations was proposed in Refs. [3,4]. The Beam Energy Scan (BES) towards the lowest energies possible at RHIC is currently under way. Significant modification of the baryon number distributions, such as its large kurtosis is indeed observed at the low energy end [5], perhaps indicating outof equilibrium fluctuations related with criticality. Using STAR detector at RHIC in a fixed-target mode is in the plans.
The topic of this paper is the baryon clustering phenomenon happening at the so-called freeze-out stages of heavy-ion collisions, in this intermediate baryon-rich domain. We will show how relatively small modification of nuclear forces at distances r = 1 − 2 fm can dramatically change the binding of baryonic clusters, as well as the kinetics of their production.
The paper is structured as follows. We start in Sec. II with a motivation for the modified inter-nucleon potential, which are quantified in Sec. III. Then we present preliminary studies of how such modifications change the binding of clusters in Sec. IV, using two opposite limits: the uncorrelated Gaussian-shape clusters in Sec. IV A, and fully correlated clusters forming certain classical shapes in Sec. IV B. Afterwards we form mean-field clusters consisting of bound nucleons only: the corresponding approximation follows the theory of globular clusters in galaxies which is briefly described in App. D. Connection to globular clusters is a new element of this work, which very instructively shows how a dynamical system can be in out-of-equilibrium partially-clustered stage, for billions of years. We use a similar approximation to evaluate the properties of mean-field self-consistent clusters of bound nucleons in Sec. IV C. Before we describe our main body of simulations, we classify the observables to be used in Sec. V. The bulk of our studies introduced in Sec. VI is done using classical dynamical approaches, such as molecular dynamics (MD) complemented by Langevin (MD+L) forces representing the effects of the mesonic heat bath. We proceed from small number of nucleons and finite clusters to rather large ones, with N ∼ 100 particles, see Sec. VII. Finally, in Sec. VIII we calculate the resulting skewness and kurtosis in a setting modeling experimental conditions of the RHIC BES program, and indeed find that its growth can be caused by the modified inter-baryon forces.
The remainder of this introduction contains a summary of the ideas motivating this work.
One important notion is the very high sensitivity of the dynamical clustering to the details of the inter-nucleon effective potential. Since the time of Yukawa's suggestion, nuclear forces are traditionally described in terms of certain meson exchanges. Furthermore, as all nuclear physicists know, any model of nuclear forces needs special tuning, needed to reproduce two delicate phenomena: (i) strong cancellation between repulsion and attraction in the mean potential energy; and (ii) partial cancellation of the remainder in the mean potential energy by quantum kinetic energy. The final result should be that neutron systems, and in fact many species of light nuclei, are not bound. Even infinite nuclear matter, with equal number of protons and neutrons (and QED effects switched off) is only slightly bound.
Because of these cancellations, a small modification of the inter-nucleon potential can induce quite significant changes in binding, even up to an order of magnitude. This is of crucial importance, because the temperatures of the hadronic phase we discuss is ranging from the critical temperature T c ≈ 120−155 MeV down to the kinetic freeze-out temperature of baryons T kin ≈ 80 − 100 MeV. Such temperatures may appear large compared to the usual nuclear potential depth ∼ 50 MeV and bindingper-nucleon ∼ 10 MeV. And yet, even with such conditions we do find significant clusters of trapped baryons. We therefore suggest to look not only at higher-order moments of the net-baryon distribution, but also out-ofequilibrium production of light nuclei.
Why do we think that inter-nucleon effective potentials might be modified in the conditions discussed, from wellknown forces in cold nuclear matter?
One generic reason-suggested many times before-is that in the baryon-rich end of the phase diagram certain modification of meson masses and couplings should be much larger than in the (well studied) small-µ B mesondominated regime. In the spirit of the resonance gas model, one may argue that there are much more baryonic resonances than mesonic ones. Studies of dilepton spectral density [6] and related ρ-meson modifications [7] have indeed shown such baryonic dominance. It is furthermore quite reasonable to think that what happens with ρ should happen with other wide resonances, the σ in particular.
Another generic reason, emphasized in Ref. [3] and also widely known, is the possible existence of the (hypothetical) QCD critical point, as the endpoint of the first-order phase transition line. On general theoretical grounds we know that second-order phase transitions have massless modes, which lead to the phenomenon of critical opalescence at scales much larger than the microscopic scales of matter. If exchanges of such long-range critical modes do appear in the inter-nucleon potential-even with relatively small coupling-we will find a significant enhancement of both the binding of certain nuclear clusters, and the kinetic clustering rates.
Finally, as multiple studies on the kinetics near the phase transitions indicate, the so-called "critical slowing down" phenomenon prevents complete equilibration, and opens the door to multiple out-of-equilibrium scenarios, some with significant cluster production.

II. BARYONIC DENSITIES AND FORCES AT FREEZE-OUT
We already mentioned the "resonance gas model", which is so successful for predicting hadronic yields for high-energy heavy-ion collisions. It is based on the standard statistical expression for the equilibrium particle densities at number of baryons of the type i where γ i , V tot , T ch are the statistical weight, total effective volume of the chemical freeze-out surface, and the corresponding chemical freeze-out temperature. We put a bar on the chemical potential indicating that we include the mean value of the inter-baryon potential in it, There are certain important distinctions between highenergy collisions and the conditions we are going to study. First of all, in the former case µ B ≈ 0 and baryons/antibaryons are both very much suppressed by the Boltzmann factor, since m i /T ch ∼ 10. Second, at T ch ≈ T c ≈ 155 MeV, excitation of baryonic resonances N * , ∆ * and their strange counterparts is very significant. For example, the population of the S = 3/2, I = 3/2 ∆ resonance, relative to that of the nucleon, is about 4 exp [(m ∆ − m N )/T ch ] ≈ 0.7. On the other hand, in the time between the chemical and kinetic freeze-out, with T kin ≈ 80 − 100 MeV, most of them decay into a baryon and one (or more) mesons, providing large "feed-down corrections" to nucleon yields.
For the conditions of the BES, on the other hand, the chemical potential is in the rangeμ = 500−700 MeV, and the Boltzmann factors exp[−(m i −μ)/T ch ] are not so punishingly small. Furthermore, the number of antibaryons is negligible, and we will not discuss them in the following. The chemical freeze-out temperature is lower, and thus a fraction of excited baryonic resonances is much smaller. In the following we will (maybe crudely) ignore their existence and feed-down. Another way to explain this approach is to assume that effectively the baryonic resonances have the same effective potentials as the nucleons. We thus normalize our calculations to the total final nucleon number observed experimentally.
Let us finally comment on the distinctions between our molecular dynamics computations and those for lowenergy heavy-ion collisions. If the temperature T ∼ 10 MeV, the thermal kinetic energy is comparable with the Fermi energy of matter at nuclear densities, and therefore quantum effects play a significant role and needs to be taken care of, by some kind of approximation. The freezeout temperatures we deal with are significantly higher, many states are excited and the role of Fermi repulsion is significantly reduced. In essence, the baryon component of the system can be approximated by a classical gas. Nevertheless, we will study some effective quantum corrections in App. E.

III. MODIFIED BARYONIC POTENTIALS
Let us now proceed to the discussion of in-matter forces between the baryons, starting with the so called "mass shifts" issue, which is somewhat controversial. On one hand, a significant part of the nucleon mass is believed to be due to "constituent quark masses" induced by chiral symmetry breaking. If so, in view that the freeze-out is not far from the restoration line of the chiral symmetry, it was predicted by many phenomenological and theoretical models that there should be a significant downward shift of such contributions to the effective quark mass. On the other hand, as we mentioned already, the successful thermodynamical description of the particle yields at chemical freeze-out uses the"resonance gas model" without any modifications of the particle masses.
Furthermore, the range of the inter-nucleon potentials is defined by masses of the corresponding mesons. For one of them, the vector meson ρ, we have direct access to its spectral density via the dilepton production, and its significant widening has indeed been observed [6]. For the ω meson no changes are observed, which is expected, since due to its longer lifetime most of them decay outside of the fireball. The σ meson, wide even in vacuum [8], is often represented as a correlated ππ pair, and is perhaps getting even wider in matter. The effective potential, convoluting Yukawa potential with its spectral density, is expected to become longer-range or even infinite-range at the critical point.
Unfortunately, lattice QCD at the moment can only extrapolate to µ B /T < 2 or so, which is far from the regime we are interested in. Some hints can perhaps be gained from the lattice study by the Graz group [9], which performed restoration of chiral symmetry "by surgery" i.e. simply removing the lowest Dirac eigenstates from the hadronic mass evaluation. What is observed is that the chiral partners (such as the nucleon P = +1 and N * P = −1, ρ and a 1 , etc.) modify their masses in the opposite directions, meeting somewhere in between. Perhaps such effects cancel each other in the calculations of the total baryon and meson yields. If so, note that the chiral partner of the σ is the pion. Moving towards it means reducing its mass, maybe to a half of it, or even all the way to zero (close to the second-order phase transition).
Completing the motivation, we now explain the reader the simplified form of nuclear forces we will be using. It follows from the popular relativistic model by Serot and Walecka [10]. One important simplifying characteristic is that it only includes the isoscalar mesons, scalar σ and vector ω, so there is no difference between the interaction of protons and neutrons. We will also ignore electromagnetism, as the clusters studied are not so large as to make it important. The Lagrangian density of their model is where the Abelian field strength of the vector field F µν ≡ ∂ µ V ν − ∂ ν V µ is the same as in electrodynamics. There are thus three fields, Dirac nucleons ψ, vector ω-mesons V µ and scalar σ-mesons φ, interacting with each other in relativistically invariant way. Their masses are considered to be an input. For definiteness we use m σ = 500 MeV, m ω = 782 MeV and m N = 938 MeV. The resulting static potential between nucleons is where the coupling values selected by Serot and Walecka [10] are  (4). Dotted line: Same as before but with the repulsive strength increased by a factor 1.4. Dashed line: Bonn N N potential [11] in channel 1 S0, taken from [12].
The ω coupling is stronger, thus dominating at small distances. Note further that these two terms nearly cancel each other, leaving us with a relatively shallow potential, |V A | < 100 MeV ∼ m N /10, see Figs. 1 and 24. It is also important to notice that the couplings are selected not to fit the binary scattering phase-shifts and deuteron binding, as done for all other phenomenological potentials, but from the fit to nuclear matter in the mean-field approximation. The details of that are further delegated to Sec. IV C.
For our studies of the baryonic clustering in this work we will use the Serot-Walecka model in four different versions of the mesonic masses: (A) The unmodified Walecka potential (3) with the parameters computed at mean field quoted in (4).
(A ) Walecka potential with increased repulsion g 2 ω → 1.4g 2 ω to make it closer to the phenomenological Bonn potential.
(B1) One in which the σ mass squared decreases "half way" (that is m 2 σ → m 2 σ /2), presumed to hold at the critical line for µ B < µ c . The "minimal modification" version changes the coupling as well g 2 s → g 2 s /2, keeping the mean potential energy constant.
(B2) This version is the same as (B1) except that the scalar coupling is not modified. The mean potential from σ thus is a factor 2 larger than in B1.
(C) An admixture of the (B2) potential with the one with very light critical mode σ, m 2 σ → m 2 σ /6 (denoted as V crit ),  (4)) is shown by the black solid line, the versions B1 and B2 correspond to the dashed brown and dotted blue dashed lines, respectively. The version C potential, with x = 0.1 is represented by the red dash-dotted line.
In Fig. 2 we show the corresponding potentials, multiplied for convenience by r 3 (note that 4πr 3 /3 times the density of other baryons tell us effectively how many "partners" a given baryon has). As one can see, these four models show progressively increasing depth and range of the attractive potential.

IV. PRELIMINARY STUDIES OF CLUSTER BINDING
Before we discuss our dynamical out-of-equilibrium studies of multi-baryon systems, it is instructive to report some simplified approaches. We considered either N = 4 − 13 nucleons, or clusters of certain fixed size, and use all versions of the modified potentials described above. In Sec. IV A we consider a limit in which there are no correlations between locations of the nucleons, so that the N -body distribution is simply factorizable into a product with the same Gaussian-like spatial distribution. In Sec. IV B we turn to the opposite limit, in which the nucleons are set to specific locations, defined by symmetry considerations, which in turn depend on the particle number, and study the dependence of the total energy on the scale parameter. Finally, in subsection IV C we calculate properties of self-consistent mean field clusters, formed of only bound nucleons.

A. Clusters made of uncorrelated nucleons
Before we study clustering rates (correlation growth), it is instructive to illustrate the effect of different potentials defined above in a simple model. Let us consider a Gaussian-shaped cluster, with the nuclear matter density n 0 = 0.16 fm −3 at its core, and the r.m.s. size R = 2 fm. The integral N = d 3 r n(r) ≈ 20, so this is a crude model of a medium-size nucleus.
Using the Thomas-Fermi expression for local Fermi momentum with γ degrees of freedom p F (r) = [6π 2 n(r)/γ] 1/3 one can calculate the kinetic energy per nucleon. For γ = 4 it is independent or R. For pure neutron matter, with γ = 2, it is K/N ≈ 27.1 MeV. Now, ignoring pair correlations n( r 1 , r 2 ) → n( r 1 )n( r 2 ) (mean-field approximation), one can calculate the potential energy corresponding to forces defined in the preceding section.
The results are P = −15.1, 3.7, 6.7, −52.5 and −70.9 MeV per nucleon, for models A, A , B1, B2 and C, respectively. In Fig. 3 we summarize our results for the total energy per nucleon, for the different potentials used in this work. We also consider different r.m.s. of the nuclear density, with a total number of nucleons of 2.5 (circles), 8.5 (squares) and 20 (triangles). The total energy per nucleon (K + P )/N for Walecka model A (whose parameters are chosen at mean field) is, within the accuracy of about 2 MeV, equal to zero for nuclear matter (p, n equal mixture), and +10 MeV for pure neutrons.
The model A is chosen to be much more shallow than the previous Walecka potential, and it is not able to bind this kind of clusters in the mean-field approximation. Similarly, the potential B 1 has comparable effect, with a total energy per nucleon of several dozens of MeV. However models B2 and C lead to an large binding. As we will detail later, the addition of binary correlation function can increase this binding even more. The main lesson from this initial calculation is that significant binding (non-negligible compared to the temperature T ≈ 100 MeV) can be produced at mean field only for significantly modified potentials (models B2 and C).

B. Clusters made of strongly correlated nucleons
For vanishingly small temperatures and small values of the particle number N , the geometry of the classical lowest energy states is suggested by symmetry. In this section we present some expectations as functions of N . Later, we will study not only the near-freeze-out T ∼ 100 MeV cases, but also cool the systems down to T ≈ 1 MeV and even T ≈ 10 −3 MeV and test that the symmetric configurations considered in this section are indeed obtained from the MD+L simulations.
For definiteness, the potential used in this section is V (r) = V A (r), which is enough to bind the nucleons as we will not account for thermal effects (T = 0).
The smallest number of particles we consider is four, N = 4, which form a tetrahedron. As it is known from studies of few-body nuclei, such correlation between four nucleons is indeed rather strong inside the 4 He, and persists in "alpha-particle nuclei" such as 12 C, 16 O. All 6 pair distances between the 4 nucleons are in this case the same, denoted by a. In general, a is defined as the minimum distance between 2 nucleons in equilibriumwhich is not necessarily the minimum of the inter-nucleon potential. The energy per nucleon V in this simplest case is just Octahedron has N = 6 particles and 15 pairs: 12 of them of distance a and 3 of distance √ 2a . The energy per nucleon is in this case The next cluster to consider, of N = 8 particles is the hexahedron (or cube). It has 12 distances a, 12 distances √ 2a and 4 of distances √ 3a, 28 in total, The largest particular cluster we discuss is the icosahedron with 12 vertices, to which we added one particle at the center, making N = 13. It has 78 distances: 12 distances at a, 30 distances at 2 − 2/ √ 5a, 30 at 2 + 2/ √ 5a, and 6 at 2a, The energy per particle V N (r) for all four clusters, as a function of the distance r, is shown in Fig. 4. One can see that augmenting N this potential energy increases, eventually exceeding the range of temperatures in the problem T = (100 − 150) MeV by a significant factor (even assuming that the potential is not modified by the temperature). Previous experience of working with strongly coupled Coulomb plasmas, see Ref. [13] and references therein, tells us that for such range of V N /T the factorized mean field theory is completely inadequate, and the correlations are significant. At the same time, this range of the ratio is also too small to cause solidification of the system, keeping the system in the strongly correlated but still liquid phase. The value of the minimal distance between 2 nucleons in equilibrium was denoted by a, and it can be obtained by minimizing the potential energy per particle in Fig. 4 for each N . For future reference, we summarize these distances and the associated potential energy in Table I. Notice that a coincides with the minimum of the potential V A (r) only for the case N = 4. Finally note that suggested by a totally different minimization problem (Thomson problem in electrostatics [14]), we have tried a different configuration for N = 8, the square antiprism, whose energy per particle is denoted as V 8s . We indeed find a lower potential energy than the cubic configuration, providing an example where the expectation based on symmetry (provided in this case by the Platonic solids) does not provide the optimal configuration. We will come back to these geometries when applying our MD+L simulations to cold systems.   principle be self-consistent, in analogy to globular clusters in galaxies. Let us assume homogeneous matter at rest, with certain mean density (1) and the mean potentialV , and on top of it a cluster, as a deviation from the mean. It is cause by a deviation of the mean potential δV (r) = V (r) −V . In thermal equilibrium it will add an extra density of baryons, Furthermore, following the setting of the globular clusters in the galaxies described in App. D, we will consider times at which all unbound particles has already left the cluster, and in the phase space integral we include only bound particles. This means in the momentum integral we only integrate over the region where To make cluster self-consistent, this extra potential δV (x) should be created by the extra density itself. We write this condition in the integral form equivalent to the Poisson Eq. (D8) for the Newtonian potential in App. D. The two equations (14,16) together make a system of equations which needs to be solved. One simplification is to ignore +1 in (14), that is proceed from Fermi to Boltzmann statistics. Note further, that when δV /T is small, one can expand the bracket to the first order in it, and then take the momentum integral using the binding condition (15). The resulting contribution is δn ∼ δV 5/2 . The exact integral without expansion can also be done analytically, leading to the following function of z ≡ δV /T , given with its (rather well convergent) series (see also Eq. (D7)). In practice we adopt the following procedure: start with a certain ansatz for δV (r), e.g. Gaussian with two parameters, the amplitude and the radius. Then, via N (z) function, calculate numerically the r.h.s. of Eq. (16), and tune the parameters to minimize the difference between the l.h.s., the obtained δV , and the input one. Of course, inside a given variational ansatz one cannot get a very good match of the shape, but the overall difference was kept at a reasonable level, of the order of 15-20 %.
We found it instructive to keep the radius of the cluster fixed, say r.m.s. radius R = 2.2 fm, and modify only the potential depth. For different potentials defined above, we find the best depth of the potential: the resulting number of nucleons in the cluster and the mean potential per nucleon in it, see Table II. One can see that while the original Walecka potential require quite deep potential and large number of baryons, the modified potentials B2, C expected near the critical point can, due to its longer range, bind a smaller number of nucleons.  N is the integrated number of baryons in the cluster, V /N and V (r = 0) are the mean potential per baryon, and the potential depth at the center.

V. OBSERVABLES
In this section we include some generic discussion of the observables involved.
The thermodynamical susceptibilities in equilibrium -derivatives of log Z over various chemical potentials of three light quarks-are usually recombined into We would call those global observables, because they correspond to mean correlation functions of fully integrated quark densities. Many of these quantities, up to N = N B + N Q + N s = 6, are currently calculated on the lattice, see Ref. [15]. For their comparison to the heavyion data on event-by-event fluctuations see e.g. Ref. [16]. At the opposite end are what we will call local observables, related to unintegrated local densities. For example, bi-local distribution function n( r 1 , r 2 ), which is usually defined in a "uncorrelated plus a correlation" form. In homogeneous matter it is defined as Similar definitions can be given for the N -point correlators. Obviously, the local observables include the full information about the correlations in the system. However, for N > 2 they are multi-dimensional functions, which is difficult to work with. Say, for N = 4 and homogeneous matter, there are 6 relative distances, and it is not practical to calculate 6-dimensional histograms. Furthermore, as we will see, in bound clusters there are strong velocity-position correlations, so that in classical approaches we adopt below one has to work with the phase space distributions, e.g. 6-dimensional one body distribution f ( r, v). Their local-in-phase-space correlators obviously are even of higher dimensions.
As a result, one needs to invent/use certain observables in between global and local ones. Experimentalists naturally use what we would call semi-global observables, in which integral is done over the detector acceptance. For example, it can be a certain range of longitudinal rapidities y ∈ [−Y, Y ] and transverse momenta p ⊥ ∈ [p ⊥,min , p ⊥,max ]. Typically, the included kinematic range is comparable to the excluded one, colloquially known as 50-50 percent setting, maximizing the fluctuations. One can measure distributions in the number of net protons P (N p ), or electric net charge P (Q), or net strangeness P (S), deduce the corresponding moments, cumulants, etc., or correlations between these charges.
As will be discussed later, for the net-proton case, the kinematical cuts imposed in experimental analyses reduce the measured multiplicities by a factor around 5 − 15% of the total multiplicity (not really following 50-50 setting). Such a reduced multiplicity allows to reach Poissonian fluctuations of protons and antiprotons, thus observing, for high energies, the Skellam expectations.
Another natural set of observables, which we would call semi-local ones, in which the densities (in coordinate space or the phase space) are integrated, but over the same small volume V (or analogous small region in the phase space). In studies of clustering we will do, the effect is of course maximal when V is of the order of the volume of the clusters produced.
The last set of observables can in fact be directly observed in experiment, via physical clusters in the final state. One well known indicator of the baryon clustering is the deuterium d production. The so-called coalescence models assume that d yield is proportional to where W d is the so-called Wigner function related to the deuteron wave function. In this case the microscopic volume V is that of the deuterons or other light nuclei, such as t, 3 He, 4 He (and hypernuclei, and their antiparticles), currently observed.
With out-of-equilibrium production of 4 He in mind, we will use below a 4-particle observable, a normalized sum of 6 inter-particle distances.

VI. MOLECULAR DYNAMICS+LANGEVIN SIMULATIONS
We use molecular dynamics (MD) simulations to study the agglomeration of nucleons and the time scales required for cluster formation. Previous applications of MD in nuclear matter to study clustering can be seen in Refs. [17,18]. We start by testing our code checking total energy and momentum conservation for finite systems, and then proceed to relatively large number of particles. These are contained in a cubic box with periodic boundary conditions and "reflections" on all sides, simulating infinite homogeneous matter. We reproduced a number of correlation functions for gas and liquid argon in a comparable regime, see App. C. We also apply the same approach to a (modified) Walecka potential to access the average properties of cold nuclear matter, introducing an effective quantum localization potential, see App. E. We relegate that study to an appendix because it turns out not to be important for systems at temperatures around the freeze-out one.
Furthermore, we modify the MD code for a nuclear system at fixed temperature, using Langevin dynamics. The corresponding stochastic forces can be thought of as interactions of ambient heat bath made of multiple mesons.
Presenting the results, we begin with systems with a small number of nucleons, starting with N = 4 nucleons. Using different temperatures, we check that they group into an average tetrahedral shape minimizing their energy per nucleon by sitting at mutual distances close to the minimum of the potential. Then, we will consider a larger number of nucleons and analyze their clustering rate. We will study the nuclear density profile of these clusters and their higher-order correlation functions.

A. Setting
A system with a small number of nucleons is useful to check and validate our MD+L code. Equilibrium configurations can be easily found for such systems. In finite systems we find no extra complications due to periodic boundary conditions, such as breaking of the periodicity of the pair-wise potential. Nevertheless, to avoid the particles to escape from the region of interest we sometimes implement a confining potential, which is written in terms of the Woods-Saxon potential where V 0 is the strength of the potential, R is the radius of the volume, and a is the skin depth. Such potential does not appreciable modify the dynamics in the region | x| 0. The temperature of the system is fixed by the light degrees of freedom (pions and kaons), which we encode in the nucleon Langevin dynamics. Therefore, we introduce a stochastic force to the nucleons as well as a drag force proportional and opposed to the nucleon momentum, with λ is the drag coefficient and ξ is the random noise following a white Gaussian distribution, with a, b = 1, 2, 3 and we made use of the fluctuationdissipation theorem to relate the drag coefficient with the variance of the fluctuation noise. A reasonable value for λ is taken from the baryon diffusion coefficient where the latter is extracted from URASiMA simulations for similar conditions of density and temperature as those used here for the hadronic evolution until freeze-out [19], which is found to be around D B 0.5 fm. Incidentally this number is not to far to the often quoted estimate using strongly-coupled QGP from holography [20] D B (2πT ) −1 for temperatures around T ch = 120 MeV.
The final value used in our simulations will be λ = 0.256 fm −1 . The precise number is not important as long as it allows for a rapid thermalization of the system.

B. Few-nucleon configurations
It is instructive to remind the different distribution of distances for the first Platonic polyhedra discussed in Sec. IV B. We summarize them in Table III.  III: Summary of the distances of edges of some polyhedra. a denotes the length of the minimal edge. Included for completeness the cube configuration does not appear as an equilibrium configuration, rather the square antiprism. We denote φ± ≡ 2 ± 2/ √ 5. We first apply this MD scheme to a system of N = 4 particles and V (x ij ) = V A (x ij ) i.e. the unmodified Walecka potential. We first run a simulation at very low temperatures T = 10 −3 MeV to match the analysis in Sec. IV B. The initial sampling of velocities is done at higher temperatures so that the particles are given some time to acquire their equilibrium configuration. The evolution of the potential, kinetic and total energies is seen in the top panel of Fig. 5. While the kinetic energy is negligible, the potential energy per nucleon takes the equilibrium value V N = −62.47 MeV (as predicted from Table I).

N = 4: Tetrahedron
It is easy to see that the geometrical configuration is the expected tetrahedron shape (Fig. 5, bottom), whose center of mass is evolving with time but the relative distances are preserved. We perform a (time) distribution of the distances between pairs of nucleons in the top panel of Fig. 6. The probability distribution function (PDF) shows a single peak at 0.873 fm, which is the expected value quoted in Table I, and corresponds to the minimum of the Walecka potential V A . The cumulative distribution function (CDF) jumps from 0 to 1 precisely at this distance (bottom panel of the same figure). These distribution functions-not particularly informative in this case-will become useful for the cases with larger N . An increase of the temperature produces a broadening of the PDF (although the tetrahedral shape is still preserved for small T ). We present the same distributions for T = 10 MeV in Fig. 7. The kinetic thermal energy is the responsible of making the average distances increase with temperature (in this case the average distance is computed as 1.03 fm), eventually preventing any kind of clustering among nucleons when the temperature dominates over the attractive N N potential.
We make notice that at temperatures of T = 120 MeV we obtain no bound system for N = 4 with the original Walecka potential V A . The clustering of four nucleons (and the eventual formation of 4 He) requires a deeper potential for the freeze-out temperatures of baryon-rich HICs.

N = 6: Octahedron
The case with N = 6 nucleons is still relatively simple to predict that the octahedron configuration will be the equilibrium shape. For T = 10 −3 MeV a fast equilibration is reached (see top panel of Fig. 8), sitting until t = 50 fm in a metastable minimum of the potential (until the last particle is finally captured by the cluster). The final potential energy per nucleon is equal to V N = −95.78 MeV, in agreement with Table I. A snapshot of the spatial configuration is presented in the lower panel of Fig. 8.
The distribution function of mutual distances is presented in the top panel of Fig. 9. It is possible to verify that the geometry is consistent with the expectations of an octahedron. This polyhedron has 2 different sets of relative distances, one at some distance a and another at √ 2a with relative strength 12-to-3. This is precisely what we observe in the histogram, where the ratio between the area under the peaks is exactly 4. This can alternatively be checked in the cumulative distribution function of the bottom panel of the same figure. The steps in this function are located at 0.848 fm, and 1.199 fm, which correspond to the 2 distances between nucleons in the octahedron configuration. The minimum distance coincides with the expectations in Table I, and the second one is a factor √ 2 larger.
In Fig. 10 we can observe how already at T = 1 MeV the two peaks are smeared out due to the thermal motion of the nucleons. Nevertheless, it is still possible to identify the octahedron configuration.

N = 8
For N = 8 we notice that the naive expectation of a cubic geometry was already ruled out in Sec. IV B in favor of a square antiprism configuration. The later configuration has a lower potential energy for N = 8 nucleons. The distribution of mutual distances is rather different from the cubic configuration case as seen in Table III. A calculation at finite temperature T = 1 MeV seems to be roughly consistent with this expectation. The PDF shown in Fig 11 is clearly inconsistent with a cubic configuration after comparing to the numbers in Table III. To test the square antiprism configuration we run a calculation at T = 10 −3 MeV. The resulting PDF shown in Fig. 12 shows that this distribution is much richer and not consistent with this geometry. The potential energy per particle at T = 10 −3 MeV is V N = −119.45 MeV, also not consistent with neither cube of square antiprism (see Table I). We were not able to identify the precise geometrical shape (shown in the bottom panel of Fig. 12), but we have classified 7 different distances with relative weights 2:4:1:2:1:2:2.  Table III.
We conclude the study of small cluster by considering N = 12 + 1 nucleons at T = 10 −3 MeV, where the expected configuration is an icosahedron plus one nucleon at the center. It is easy to see by naked eye that the geometrical configuration resembles this expectation. In the top panel of Fig. 13 we present a snapshot of the spatial configuration at some time after equilibration. In the middle panel we also present the distribution of (78) mutual distances. We observe 4 different sets of distances, and with a relative weight (see cumulative distribution function in the middle of the same figure) in excellent agreement with the expectations of Table III. Finally, the minimum distance (position of the first peak of the distribution) is 0.782 fm, and the potential energy per nucleon obtained is V N = −177.32 MeV, both in agreement with the values in Table I.

C. Clustering at freeze-out temperatures
In this section we describe simulations following the scheme presented in the previous section (MD+Langevin with modified Walecka potentials). The number of nucleons is large N = 128, and the temperature is fixed at the typical freeze-out temperatures T kin = 120 MeV [21].
In this section we use the potential V A to see how a deep potential can bind nucleons and eventually produce a large cluster. In this example we look for a clear example of such a cluster, and to apply various systematic procedures to analyze its internal structure. The initial state and a configuration after its equilibration are shown in Fig. 14. The initial geometry is spherical with a density of n = 0.16 fm −3 .
When equilibrium is reached we obtain a big cluster which includes all N = 128 particles. In Fig. 15 we show the time evolution of the kinetic, potential and total energies (top panel) as well as the temperature evolution (bottom panel). We observe that the total energy of the system is dominated by large negative potential energy, so to see one cluster structure is not surprising. The temperature at equilibrium (plateaux formed after t ∼ 50 fm in the bottom panel) fluctuates around the value T kin = 120 MeV. The sudden kicks in temperature and the steps in energies occur when one more particle is captured by the cluster, and falls to the deep potential well.
Following the mean-field approach, and the King's solution a decreasing distribution of particles is expected as a function of the radial distance. We want to analyze the internal arrangement of nucleons in this cluster by finding the radial distribution of nucleons starting at r = 0 (defined as the centroid of the cluster). As the centroid evolves in time, we monitor its position at each time step, and perform the radial distribution of the nucleons. To have independent events in the distribution, and to avoid spurious correlations we choose time steps  well separated to perform the average. We measure the number of nucleons per unit volume/radial distance, and plot it versus the distance from the centroid.
The distributions dN/dV and dN/dr are represented in Fig. 16 showing a non monotonous structure suggesting a shell like organization with accumulations of nucleons every 0.3 fm in the radial direction.
For clarity, let us note that this study is only done for investigative purposes. The time scales considered in the plot above, up to t ∼ 300 fm/c, are much longer than those available in heavy-ion collisions, t ∼ 10 fm/c. Furthermore, this analysis was done in static, rather than exploding, heat bath. So, by no means we suggest that such clusters are actually produced in experiment: at best we hope to find evidences of the very beginning of the clustering process.

VII. BARYONIC CLUSTERS NEAR THE CHIRAL TRANSITION
In this study we continue the study of big cluster formation, and their time scales as the clustering process becomes more and more important. This will happen when the original parameters of the Walecka potential are modified as a consequence of the changes in the properties of the σ mode.
We will compare the potentials V B1 , V B2 , V C , each one thought to be acting closer and closer to the chiral transition.

A. Formation of clusters
All simulations begin with randomly placed nucleons. Naturally the cluster formation starts with small clusters, which then assemble into larger and larger ones. We decided to follow the process by defining variables in which one can separate clustered and non-clustered baryons in the most direct way, and then histogram the distributions at different times.
We performed a number of such studies, demonstrating here one example, for 4-particle variable. The variable S (from sum) is defined as the normalized sum of all mutual distances between particles in the system, where i, j = 1, ..., N run over all nucleons, x ij is the vector joining pairs, and N d = N (N − 1)/2 is the number of mutual distances between different nucleons. As one can see from an example shown in Fig. 17 (note the scales), for potentials A, B1, B2 we observe that the entropy wins over the energy. With time the distribution slowly become wider due to the diffusion of baryons. In contrast to that, the potential C with longer-range attraction shows the opposite trend, the potential wins over the entropy, leading to a rather robust clustering. An example of the time evolution for the C(x = 1) potential is shown in Fig. 18. The clear separation of the distribution into two peak structure, in this one particular event, corresponds to a formation of two clusters (in this event, those have sizes of 9 and 22, with only one particle evaporating out). The first peak corresponds to intra-cluster distances in both clusters, whereas the second peaks reflect inter-cluster distances.

B. Time scales
We consider a system of N = 128 nucleons at temperature T = 80 MeV, with an initial nuclear density n = 0.13 fm −3 and finite size. In Fig. 19 we show the time dependence of the energies per particle (top panel) and the temperature (bottom panel). After a fast thermal equilibration the temperature is approximately constant, while the total energy is not conserved in the evolution as dissipation occurs due to the drag force.
The potential V B1 is able to produce full clustering after long times. From the example in Fig. 19, the full equilibration time is of the order of ∼ 800 fm/c, and clustering has taken place. Individual particles can escape the cluster thanks to their kinetic energy, however, we avoid the lose of particles with the external trapped potential in Eq. (21).
We can define an equilibration time by noticing that the total energy has an approximate exponential decay exp − t teq . We obtain t eq = 187 fm in this particular example.
Although the V B1 is enough to form a big cluster in several hundreds fm/c, these scales are totally irrelevant for HICs. A slightly critical potential V B2 produces the clustering in a much faster way ∼ 40 fm/c. We present the time dependence of the energies and temperature in Fig. 20.
In this case, the exponential decay is much less evident. We find an initial regime of ∼ 10 fm/c where the energy is approximately constant. Between 10 fm and 17 fm we find a good exponential decay with an equilibration time of t eq = 3 fm/c. After this transient exponential decay the relaxation is much softer, reaching equilibration in around 40 fm/c. The time scales for the clustering with this potential are much closer to those in heavy-ion collision, so it seems reasonable to consider this mechanism as potentially important close to the critical point (where the equilibration time is even reduced using a deeper potential like V C ).
Nevertheless, it seems clear that the time for full clustering is still large to take place completely in heavy-ion collisions. We only hope to have a potential effect close to the critical point where the signatures of initial clustering might certainly occur (perhaps clusters of few nucleons as 4 He). Starting from a system away from the critical point, we will calibrate our model with noncritical potential V A to experimental data at energies where Poissonian fluctuations are observed. Then, we will modify our potential to increase criticality and compute observables like higher-order moments of the (net-)proton distribution.

VIII. BARYONIC CLUSTERS IN BES CONDITIONS
In this section we apply our model to heavy-ion collisions in the condition of BES. Rather than providing a quantitative result, for what one would need a more sophisticated evolution code, we contempt ourselves to show that the effect is consistent with what is observed experimentally using the closest experimental conditions we are able to implement.
At high collision energies above √ s N N = 19.7 GeV STAR data shows approximate Skellam distribution for net protons [5,22] consistent with thermal equilibrium fluctuations. This will be our baseline energy. For collisions at this energy and below we can safely neglect antiprotons. In Table IV we show the ratio of proton/antiproton in the kinematic cut |y| < 0.5, 0.4 < p ⊥ < 0.8 (GeV/c) for the most central collisions at different energies considered in this work.
In our simulations, Poisson statistics is achieved when measuring the distribution of protons in a sub-volume or the order % of the initial volume. This is consistent with the fact that experimental net-proton distribution in the narrower p T cut is 5% of the total net-proton multiplicity [23], and matches very well Poisson expectations.
Our first task is to achieve similar multiplicities in a noncritical scenario, where Poissonian fluctuations dominate, which we identify with experimental data at IV: Proton-to-antiproton ratio for |y| < 0.5 and 0.4 < p ⊥ < 0.8 (GeV/c) at centrality bin 0-5% for collision energies √ sNN ≤ 19.6 GeV. From Ref. [22].
√ sNN (GeV) 7.7 11.5 19.6 proton/antiproton 114.4 ± 0.6 30.64 ± 0.07 9.89 ± 0.01 √ s N N = 19.6 GeV. For that energy we will be running the potential V A , i.e. Walecka potential with an additional repulsion. The kinetic freeze-out temperature for low energies is roughly T kin = 120 MeV. In our code, the MD will simulate a few Fermi/c between hadronization and freeze-out, so we set the temperature to T = 150 MeV (although the calculation is not very sensitive to this parameter). The baryon density is close to n kin ∼ 0.12 fm −3 at freeze-out, but at earlier stages it can take a few times this value [24]. We use a value of n = 0.3 fm −3 . For numerical convenience we use a reduced number of protons N = 32 and then scale up the different cumulants as suggested in Ref. [25], to be able to compare with the experimental cumulants. In particular, we note that scaled skewness Sσ and kurtosis κσ 2 do not depend on volume, so we can compare them without the need of scaling. The number of events for each potential is N ev = 10 5 .
We summarize the parameters used in our MD+L simulations for this section in Table V. These parameters will be common for all potentials, as we would like to isolate the only effect of the potentials. In addition, we checked that the parameters from thermal fits do not change too much within these energies (the most sensible parameter would be the baryochemical potential).  We run our MD+L a total time of ∆t = 5 fm/c, corresponding to an approximate time between hadronization and kinetic freeze-out. While this time can be extended, perhaps up to a factor 2, we prefer to be conservative not to overestimate the effect of clustering (as we have seen larger times help to create more bound clusters).
The calculation is performed in a non expanding medium, i.e. without radial flow implemented. This is convenient for the use of nonrelativistic dynamics at all times. A final boost in rapidity and transverse momentum will take care of the mapping of the particles into the appropriate kinematic domain, consistent with experiment.
A first simulation using V A generates the expected Gaussian distributions of rapidity and transverse momentum for nonrelativistic dynamics. The fitted p ⊥ distribution is perfectly consistent with the temperature used (providing an intermediate check of the code). Then, we perform a mapping of these distributions to mimic the experimental findings and take into account radial flow. Otherwise, the distribution of particles cannot match experiment. First, the experimental p ⊥ distribution of netprotons is well fitted to a double exponential [21]. However, we have found that within the same kinematic cut a Gaussian form is already very approximate. Therefore we opt by simply rescaling the transverse momenta of our particles by a single factor of 1.7 to match experimental distribution.
The rapidity distribution is a relative narrow Gaussian one centered at zero. Therefore, most of the particles live at midrapidity, which is not consistent with experiment, where one expects ∼ 10% of particles at midrapidity. This is normal, as in our simulation we do not account for any longitudinal boost. Then, we also transform our Gaussian distribution into a uniform distribution in rapidity between the kinematic limits for that energy. The uniform rapidity distribution is quite a crude approximation, but it does not involve any additional parameters.
Once fixed our final distributions, we proceed to test the proton distribution by performing the kinematical cuts considered in the literature. In what follows we will denote Cut 1 as the one with rapidity |y| < 0.5 and transverse momentum 0.5GeV/c < p ⊥ < 0.8GeV/c; whereas Cut 2 extends the p ⊥ coverage up to 2 GeV/c, thanks to the time of flight detector for the particle identification [5].
We summarized our results for the proton cumulants in Fig. 21 compared to experimental data for the two mentioned cuts. In our simulation we simply perform the analysis within the kinematic cuts including the N ev = 10 5 events. The statistical uncertainty is coming from the Delta theorem, as explained in Ref. [26]. We need to make a final scaling from our N = 32 to match the absolute number of protons observed in experiment. For this we choose the average number of proton C 1 in Cut 1 and scale up our value to the experimental results. As it is well known, all the cumulants scale with volume, so all the other cumulants (up to order fourth) are scaled up by the same amount.
We observe that despite our crude model, we can match in a good degree of accuracy the cumulants for proton number, in both kinematical cuts for √ s N N ≤ 19.6 GeV. Therefore, the noncritical scenario is reasonably under control. As a last check, after multiplication of dN/dy by the same scaling factor we obtain dN/dy = 30.5 at midrapidity. The experimental value given for the most central collisions is dN/dy = 34.2 ± 4.5.
It is also possible to compute ratios of these cumulant to obtain the skewness Sσ = C 3 /C2 and kurtosis κσ 2 = C 4 /C 2 and compare with experiment. As being directly related to cumulants these moments for proton . Experimental data for Cut 1 is taken from [22] and for Cut 2 from [5]. distribution will be in accordance within the same levels at those in Fig. 21.
A step forward would correspond to compare our skewness and kurtosis with the same quantities for the netproton distribution. However, at this energy there is an additional 10 % systematic error because antiprotons (not accounted in our simulation) are still important (see Table IV).
Once the obtained multiplicities and cumulants for the proton distribution are calibrated to those as in experiment, we simply repeat the calculation with different potentials at our disposal. Each of them are supposed to encode the modification of the NN potential closer and closer to the QCD critical point. Notice that a precise matching between the experimental collision energy and the particular NN potential is far from clear, so we cannot compare directly to our results with real data. We can nevertheless observe the qualitative effect on the skewness and kurtosis after increasing the criticality of our model. We use Models A, B1, B2 and Model C with x = 0.1, 0.5, 1.
In Fig. 22 we our present our results for Sσ and κσ 2 in our simulations. From top to bottom we show the theoretical skewness as a function of the NN potential, the experimental skewness as a function of the collision energy, and the same dependences for the kurtosis, respectively. In all cases we consider both Cut 1 and Cut 2. As mentioned, a direct comparison is not possible due to the difficulty of matching a given potential to a precise collision energy. However, we base our study in the idea that lowering the collision energy from high energies, should necessarily approach the expanding system to the critical region, until some particular value of √ s N N . In our setup this is achieved by increasing the attraction of the NN potential towards a more critical one. One important result is that the increase of the kurtosis is consistent with the presence of a critical point. Therefore, clustering of nucleons close to T c , and the increase of NN correlations, translated into an increase of higher moments.
In this document we have proven that the nuclear clustering is a solid phenomenon with relatively conservative assumptions if the system is left for large amount of time. In this respect, we do not find realistic to experimentally find big clusters of nucleons due to this phenomenon, as the required time for this is much longer than the hadronic phase. However, early signatures of clustering can be reflected in higher-order cumulants of (net-)protons, and we showed that these signatures are compatible with what has been preliminarily observed in experimental data.

A. Light-nuclei clusters at freeze-out
The results of the previous sections show that not only a strong correlations among nucleons, but also a clustering of few of them might be possible during a time interval of several fm/c. A mechanism producing such effects between nucleons is required, according to Ref. [27], to explain STAR experimental data. In that reference the third and fourth order cumulants cannot be explained by a model with nucleon stopping and baryon global conservation alone. The conclusion of [27] is that some sort of clustering is needed to describe the data. In this work we provide such a natural mechanism for clustering, if the N N interaction is attractive enough close to the critical point.
While for the calculation of the higher-order moments of the proton distribution we have included the contribution of all protons-within the corresponding kinematic cuts-we will now extract the nuclear clusters which may give rise to light nuclei at the end of the evolution. We will denote these clusters as "pre-nuclei" as they are products of the nucleon coalescence at freeze-out. If such "pre-nuclei" are able to survive as bound objects until the final stage of the fireball at very low temperatures, then it will form states like 3 H, 3 He, 4 He...resulting in an excess of light nuclei over the expected thermal production.
Let us focus on 4 He, which has a binding energy of 28.3 MeV [28]. Freeze-out temperatures are larger than this energy, but the modification of the nuclear potential can provide extra binding to it. We will look for candidates of 4 He nuclei at the final time of our simulation. We just need to find 4 nucleons close in the phase space at the moment of the freeze-out. If a pair of nucleons are separated by a large distance in the phase space, then they are assumed not to belong to the same cluster. We run our code using several versions of the N N potentials and identify configurations of 4 nucleons, or "pre-4 He". We apply the following criteria: 1. We only search clusters of 4 nucleons. If any nucleon also belongs to a different cluster, the whole set is ruled out (as the nucleons belong to a bigger nuclei).
2. The relative position between pairs of nucleons should be small. The rms of 4 He is 1.67 fm, and the rms for proton is 0.87 fm. Assuming a tetrahedron configuration (see Sec. VI), one obtains that the distance between the center of 2 nucleons should be 1.69 fm. Giving some freedom to this value due to the thermal motion (deformation of the tetrahedron), we assign a maximal distance of ∆r = 2 fm.
3. The momenta should also be similar. Taking a typical thermal momentum for 4 He of √ mT = 0.77 GeV, this gives a momentum of 0.11 GeV to each of the Cartesian components of each individual nucleon. We impose the condition that any component of the relative momentum cannot be larger than ∆p = 0.22 GeV.
These two numbers satisfy ∆p∆r = O(1), so this choice seems reasonable. In Fig. 23 we present the number of clusters of 4 nucleons ("pre-4 He") per event (we use N ev = 10 5 events for each N N potential). Similarly to the calculation of the moments of the proton distribution, we associate the noncritical potential V A to the experimental collision energy of √ s N N = 19.7 GeV. Then, using this potential in our MD+Langevin code we count the number of such clusters defined by the previous criteria.
The results illustrate the effect of larger clustering formation with the N N potential used. Going from right to left in Fig. 23 we find that the number of "pre-4 He" increases with the attraction of the nuclear potential. Surprisingly, the most attractive potential V C presents a decrease of the number of clusters. The explanation is that because the huge attraction the nucleons belong to bigger clusters, i.e. it is more difficult to find 4 nucleons isolated from the rest.
In spite of the qualitative motivation we can check that the numbers are not unrealistic for the V A potential (identified with the collision energy of √ s N N = 19.7 GeV). The multiplicity of 4 He is not measured at this energy by STAR. However, other light nuclei have been measured by the NA49 experiment [29]. For a close energy of √ s N N = 17.3 GeV (E beam = 158AGeV) we know that after increasing the mass number in one unit (from A = 1 to A = 2, and from A = 2 to to A = 3), the dN/dy at midrapidity decreases a factor of 100 [29]. Assuming that this scaling holds also up to A = 4 we then expect one nucleus of 4 He for each 10 6 protons. Using that in our simulation we have N ev = 10 5 events with N = 32 nucleons we expect a total of 3 nuclei of 4 He for the V A case. We have numerically obtained 11 "pre-4 He". Given the number of simplifications made in our study, this number seems reasonable due to the following argument. Notice that our 4-nucleon clusters might not necessarily become 4 He at the end of the evolution. Assuming a rather sharp freeze-out process, one would need to project the Wigner function of these "pre-nuclei" configurations to the actual wave function of 4 He. This study-which would give a more precise prediction for the produced light nuclei-is left for a future work, and here we restrict ourselves to show the qualititative increase of clusters with the reduction of the σ mass close to T c .

IX. SUMMARY AND OUTLOOK
In this work we have studied baryonic clustering at the freeze-out conditions corresponding to baryon-rich heavy-ion collisions. More specifically, we have observed that both the clustering rate and the properties of the resulting clusters are very sensitive to the magnitude of the effective inter-nucleon potential, and suggest that detailed studies of the baryon distributions will be able to fix such potentials, and ultimately tell us whether the QCD critical point exists or not.
In Sec. III we have defined a set of inter-nucleon effective potentials, which are modifications of the Walecka-Serot model, some with the addition of a long-range component related to massless critical mode at the (hypothetical) critical point. Then in Sec. IV we performed some initial studies of baryonic clusters which such potentials can support. The main tool we used is classical molecular dynamics, complemented by Langevin stochastic terms accounting for the mesonic heat bath, and also by additional repulsive potential modeling quantum Fermi effects for the case of cold infinite nuclear matter in App. E.
If the matter is not exploding and the system evolves long enough, we do observe that the initial stage, with random baryon positions, is always clustering, in one or few large clusters. If the time is not so long, corresponding to ∆t ∼ 5 fm/c available for the hadronic phase of heavy-ion collisions, the degree of clustering is very strongly dependent on the version of the potential used. Our main result is thus the high sensitivity of this phenomenon to the inter-nucleon potential.
We also tried to imitate an experimental fireball, mapping it to an expanding system. We also impose similar cuts to the experimental acceptance of STAR papers, and calculated the baryon number distribution. We do observe an increase of kurtosis, by about a factor of 3, from the original Walecka potential to our most attractive version.
Our main qualitative conclusion is that while the evolution time available is insufficient to produced fully developed "nucleosynthesis" with large clusters, one definitely should find the baryon distribution in the final state far from thermal equilibrium. Indeed, the confidence in this statement is also provided by similar studies in atomic systems and globular clustering in galaxies (briefly outlined in the corresponding appendices). We therefore suggest to look at possible deviations from thermal equilibrium in the yields of light nuclei, such as d,t, 3 He, 4 He.
While in this paper we cannot directly compare our results to the STAR BES data, we do focus on one important finding: a growth of the kurtosis of the proton distribution near mid-rapidity, at the lowest collision energies [5].
Although the specific critical enhancement of the multi-particle fluctuations remain the major goal of this program, one needs to also study other phenomena which can lead to those. In this paper we focused on the clustering of baryons due to their attractive interaction. As we detailed above, significant clustering should in fact occur due to the usual nuclear forces.
Appendix A: Mean-field approach to the Serot-Walecka model In this appendix we remind the reader a simplified form of nuclear forces, following a model by Serot and Walecka [10]. One important simplifying characteristic is that it only includes the isoscalar mesons: scalar σ and vector ω, so there is no difference between coupling to protons and neutrons.
Its Lagrangian density is shown in Eq. (2), and the inter-nucleon potential is written in Eq. (3), which we reproduce here again for convenience, with parameters in (4) chosen by mean-field calculations.
In Fig. 24 we illustrate the partial cancellation of the attractive and repulsive terms of this potential.
Considering the case of infinite homogeneous matter of density n and ignoring correlations between the nucleons, one gets the mean potential energy If matter is cold, T = 0, the baryons are in a form of degenerate Fermi gas of quasiparticles with dispersion relation Note that if one expands the square root, the leading term of the mean potential g ω V 0 − g σ φ 0 will be the same as the one from the usual non-relativistic theory, but the kinetic energy term would be k 2 /2M * rather than k 2 /2m N . The total energy of the gas is where the statistical weight γ = 4 for symmetric nuclear matter, and 2 for neutron matter in neutron stars. Two densities, the vector and scalar, can now be written as integrals over the Fermi sphere Note that the latter has scalar mass M * in the numerator and the energy in the denominator, which is needed because Lorentz invariant integration measure is d 3 k/E k .
At this stage all is fixed except the scalar mean field (or alternatively M * ): this is a parameter of the homogeneous-field trial function, which as any variational parameter, it should be found from minimization of the ground state energy. This leads to the following equation for M * to be solved numerically. As shown in the original work [10], such mean-field result can be fitted to reproduce the nuclear matter density and nuclear binding.
For finite spherical nuclei the procedure includes the solution of the mesonic equations of motion supplemented by Thomas-Fermi-like treatment of baryons. For heavy nuclei the results are rather good. While this model is only a stripped-down version of nuclear forces and the mean field is only the first of various approximations used for nuclear matter description, we will use it below due to its simplicity. In particular, this model only includes isoscalar exchanges, which means that pp and pn forces are the same. As a result, the only place where isospin matters is in the quantum kinetic energy, since it depends on the number of species. We are however fully aware of the fact that Walecka model parameters are only good for mean field treatment, and the resulting forces do not describe elastic N N scatterings or the deuteron binding. To improve on this it is possible to increasing the repulsion, via higher ω coupling g 2 ω → 1.4g 2 ω , so that the resulting potential gets very similar to the Bonn potential [11] (see right panel of Fig. 1). This is what we will call "modified Walecka potential", which will be used in this work and denoted by V A .
Appendix B: σ-meson dependence of the N N potential from the functional renormalization group.
In Sec. III we have analyzed the modification of the attractive part of the N N potential due to the σ mass modification. In this work we have not dealt with the precise dependence of this mass with the temperature/density. This would imply an additional uncertainty dependent on the model used e.g. linear sigma model, quark-meson model, Nambu-Jona-Lasinio model... On the other hand a more rigorous treatment would involve the modification of the whole spectral function of this state.
In this appendix we will illustrate how these two issues can be addressed using results of the σ-meson properties in the N f = 2 quark-meson model, approached by the application of the functional renormalization group (FRG) [30,31].
The version of the quark-meson model presented in [30,31] contains quarks, antiquarks, pion and σ degrees of freedom. After the evolution of the FRG equations one is able to obtain the medium-modified properties of these states in the infrared limit. In particular, the spectral function of the σ meson can be obtained at different temperatures and chemical potentials.
In Ref. [30] the critical point of the quark-meson model is located arount T 10 MeV, which seems to be quite low from the phenomenological point of view (this critical temperature is supposed to increase when extending the calculation to N f = 3 flavors and after introducing the effects of the Polyakov loop potential). Therefore, we will consider two cases: at T = µ = 0 where the σ screening mass is very close to our vacuum mass m σ 500 MeV for the potential V A ; and T = 150 MeV where the σ screening mass drops to values around m σ 280 MeV.
The attractive part of the static N N potential is computed as a Fourier transform of the σ-meson exchange diagram, where p · x ≡ p 0 t − p · r, and the σ retarded propagator is used in the Lehmann representation to account for the complete spectral function ρ σ (ω, p), The data from Ref. [31] is given between ω ∈ (−1, 1) GeV and p ∈ (−1, 1) GeV. For the case at T = µ = 0 MeV the σ-mass pole lies away from the real energy axis (i.e. its real part is above the π − π unitary threshold). Therefore, the mass appears in the spectral function as a broad pole, which can be numerically integrated in energy and momentum. However, at T = 150 MeV, µ = 0 MeV the σ mass goes down below the unitary threshold and the pole is located on the real axis. In such a case the mass appears in the spectral function as a Dirac delta, which we need to add by hand to the spectral function as the discretized data cannot capture it. Following the conventions in [31] we add to the ρ σ (ω, p) a term like: where m σ is the pole mass of the σ and Z −1 is the pole weighting factor (we refer to [31] to see how to compute these from the spectral function).
In Fig. 25 we compare the results coming from the Fourier transform of the spectral function of the σ meson, and the simple potential as given in our Eq. (3) (the repulsive part due to the ω meson is kept the same). Even at the quantitative level the two sets of potentials look similar. The main difference of the σ potential occurs at small distances, but in this limit the full potential is dominated by the ω repulsion. In the results using the σ spectral function we observe some spurious oscillations around zero, which are nothing but the Gibbs effect  coming from the inverse Fourier transform of the spectral function performed within a compact support of energy and momentum.

Appendix C: Kinetics and clustering in atomic systems
The simplest atomic systems are those of the noble gases, with spherical atoms and forces depending solely on distances. For a large enough atomic weight, one can neglect quantum effects. For all these reasons, the object of choice is argon, with its A = 40 (for the most abundant argon isotope) being ten times heavier than 4 He. By tradition, theoretical studies of it use the simple potential Its minimum is at V (2 1/6 σ) = − . Following one of the classic MD simulations from 1960's [32], one can use parameter values σ = 3.4Å, = 120 K. The shapes of this potential and that of the nuclear forces (Walecka model) are compared in Fig. 26. It shows that Lennard-Jones potential is much more narrow. The ratio of the potential to the temperature are similar to the problem we study, provided the temperature of argon is T ∼ 100 K.
The work [32] focused on one temperature T = 94 K and one density ρ = 1.37 g/cm 3 , which is well in the liquid phase. We minimally modify our MD (without Langevin dynamics) to run an isocanonical simulation (by rescaling of instantaneous temperature). We Increasing the density one crosses the phase transition to a solid like phase. A new simulation with n * = 1.1 gives the radial distribution function in the lower panel of Fig. 27. The amount of very pronounced peaks is a signature of the solid (crystalline) structure of the system. In this case the distribution of peaks can be identified with a face-centered cubic distribution (which is the configuration used to initialize our simulation).
The standard MD simulation, unlike Monte Carlo ones, have not just static (fixed time) but also the timedependent information, such as velocity-velocity and other correlation functions, Using standard Green-Kubo formulas one can calculate diffusion constant, viscosity, etc. We however would not go into vast literature on the kinetic properties, except to note that liquid argon, like other liquids, has a second order critical point, and studies of the singularities of kinetic coefficients there remain to be better understood.
Finally, we would like to mention instead a particular large-scale MD simulation [33], using as many as a billion atoms, and focusing on transition from homogeneous particle distributions to liquid phase, at supersaturated conditions. As it is well known, the process can be divided into two stages: (i) creation of critical clusters, with i * particles in them; and (ii) their subsequent linear growth as a function of time with certain rate. Large scale of the simulation had allowed to cover a range of temperatures and densities, in which the clustering rates change over many decades, and cluster sizes grow to well over 100 particles. However, what is most important, is that in all cases the critical clusters are relatively small, ranging from i * ∼ 12 to about 100 atoms. Therefore, the classical theory of nucleation-treating these clusters as macroscopic drops with a surface and volume free energies-needs to be corrected. After the actual energies of these clusters are used, the corrected theory was shown to work well. Equilibrium configuration of small and medium size cluster in Lennard-Jones interaction has been studied e.g. in [34].
which is a shif ted Maxwell-Boltzmann distribution with temperature T = σ 2 .
Step three is a calculation of the corresponding density of stars, which includes the integration over the velocity. Note that it is limited by the escape velocity defined via the potential, so the density obtained is the function of the potential ψ = Φ − Φ 0 , This complicated function is plotted in Fig. 28 , and one can see that it is a monotonously rising one. The density is the source of the potential itself, so now we come across the main dynamical equation to be solved, the Poisson equationn for the potential. In case of spherical symmetry it is 1 r 2 d dr r 2 dψ dr + 4πG N ρ K (ψ(r)) = 0 , which can be solved numerically starting from the center. The value ψ(0) is the single input parameter, the derivative needs to be vanishing at the center ψ (0) = 0. Solution can be followed until the point where ψ = 0: and as it is clear from the expression above for the density, at that point the density vanished as well since the integration region till the escape velocity shrinks to zero. Substituting the resulting ψ(r) into the universal ρ(ψ) one finally obtains the spatial distribution of the stars in the cluster.
Appendix E: Cold nuclear matter and quantum Fermi repulsion To account for quantum repulsion in the simulations, the simplest thing one can do is adding the Fermi energy, evaluated in spirit of the Thomas-Fermi approach from the density profile, to the classical kinetic and potential energy. Full account for both quantum and thermal fluctuations can be done in approaches called "quantum open systems", see e.g. Ref. [35] for its application to motion and heavy quarks and quarkonia, as well as general references.
Strictly speaking, a complete account for quantum effects would requires going from classical molecular dynamics to full path integrals. As it is well known, while for distinguishable particles and bosons it can be considered to be just a technical complication, for fermions the amplitude needs to be anti-symmetrized, with brings in the notorious sign problem. The effective Fermi repulsion, acting as a kind of repulsive potential, generates correlations between particles which depend strongly on their mutual distance.
In 1980's O. Zhirov and one of us studied paths of fermions moving in one dimension. This case is special because one can always enumerate fermions along the line, and thus pretend that the "exchange" never happens. For a small time step t a the one-particle amplitude is and for two particles it can be written as with the "Pauli potential" defined as Note that when two particles get close, the exponent becomes close to 1, the argument of the log near zero and the potential gets very high. So, a node of the amplitude can be viewed as a repulsive potential. This is going in the right direction: indeed, fermions must have a larger energy than distinguishable particles in the same setting. If particles never jump over the note-generated barrier, their order along the line remains preserved, and if the Pauli potential is included, the simulation can be done by traditional Monte Carlo. We checked it for several (n = 3 − 5) particles put in a harmonic potential: for distinguishable particles the ground state energy is ω(1/2 + 1/2 + ...) but for fermions it should be ω(1/2 + 3/2 + 5/2 + ...) since each must be put into the next available level. So, our algorithm with this "Pauli potential" worked correctly. The work was concluded in Ref. [36]. Description of the method and its usage is also described in Ref. [37], in which many tests have been successfully performed.
The next step forward, allowing to use this idea in any dimension, was made by Ceperley [38]. It has been applied to fermionic problems, including liquid 3 He. The main idea can be explained if one considers various paths of one fermion, keeping all other fermion paths frozen. The 1-dimensional node of the amplitude gets promoted to a "nodal surface", which surrounds each fermion, keeping it inside a "nodal cell". Paths which are not allowed to leave the nodal cell are called "restricted": sum over the restricted paths obviously has no sign change.
The nodal surface model corresponds to a certain constant potential well with the location of of the wall depending on those of other particles. The radius of this surface can be tuned to reproduce the Fermi energy of an ideal gas.
Instead of a sharp wall we decided to include a more smooth localization potential, of the form where the exponent is chosen rather arbitrary as long as Pauli repulsion as short distances is achieved.
To normalize this effective potential we attempted to simulate properties of cold homogeneous nuclear matter by our molecular dynamics scheme, where V represents the pair-wise potential, sum of the localization potential plus one among the different possibilities described in the main text. We use the Walecka potential V A (r) with increased repulsion, whic is closer to the NN phenomenological potentials for nuclear matter. To simulate an infinite system be work on a cubic box with periodic boundary conditions. In such a box the particle density is fixed to the nuclear density at saturation n 0 = 0.16 fm −3 with N = 64 the total number of nucleons. To account for the interactions of the particles in the box and those outside, we use the method of images, where in the sum of Eq. (E5) we consider the contributions from all j-particles within a number of copies of the box in each spatial direction (positive and negative). The number of images (or copies of the elementary box) per each direction is set to 2. After a transient regime, the MD simulation reaches an equilibrium state with constant potential and kinetic energies (with statistical fluctuations of O(1/ √ N )). For infinite nuclear matter at saturation an average Fermi momentum of p F ∼ 260 MeV translates into a kinetic energy per nucleon of K/N ≈ 25 MeV. Lacking of quantum dynamics in the classical MD we achieve this value of K/N by forcing a isokinetic simulation by rescaling the velocity of each particle by K/K inst , where K inst is the instantaneous value of the kinetic energy at a given time step. The expected energy per nucleon at saturation E/N = −16 MeV provides the additional constraint that helps us to fix the remaining parameter of the simulation, the strength of the localization potential, to a = 0.75 fm 3 . The resulting energies versus time are given in Fig 30. After the equilibration time (∼ 100 fm/c) we can measure the average total energy (binding energy) per nucleon. We obtain −16.6 MeV, a fair value for our illustrative purposes. For dedicated computations a more precise value of a can be extracted, using more nucleons in the simulation in order to reduce the statistical fluctuations of E/N (going as 1/ √ N ). We find a rather homogeneous system at equilibrium with evidences of a slight grouping of nucleons. In the upper panel of Fig. 31 we show the initial configuration of nucleons at random positions in a volume of (7.37 fm) 3 . In the lower panel we show the spatial configuration of the nucleon for an arbitrary time well after the equilibration time. Quantum effects via localization potential-important for a T 0 calculation-will be absent around the freezeout temperatures, where kinetic energy is expected to be dominated by thermal fluctuations.