Correlated stopping, proton clusters and higher order proton cumulants

We investigate possible effects of correlations between stopped nucleons on higher order proton cumulants at low energy heavy-ion collisions. We find that fluctuations of the number of wounded nucleons $N_{\mathrm{part}}$ lead to rather nontrivial dependence of the correlations on the centrality; however, this effect is too small to explain the large and positive four-proton correlations found in the preliminary data collected by the STAR collaboration at $\sqrt{s}=7.7$ GeV. We further demonstrate that, by taking into account additional proton clustering, we are able to qualitatively reproduce the preliminary experimental data. We speculate that this clustering may originate either from collective/multi-collision stopping which is expected to be effective at lower energies or from a possible first-order phase transition, or from (attractive) final state interactions. To test these ideas we propose to measure a mixed multi-particle correlation between stopped protons and a produced particle (e.g. pion, antiproton).

The LQCD results show that, at the physical pion mass, hot nuclear matter exhibits residual properties of both dynamical chiral symmetry breaking and confinement at finite temperature. In the LQCD calculations, it was demonstrated that the transition between hadrons and quark/gluon degrees of freedom is of crossover type [25]. At finite baryon densities, progress of LQCD calculations is impeded by the notorious sign problem. Several attempts to circumvent the sign problem were carried out [26][27][28][29][30][31]. Some of these studies indicate the existence of the expected critical point (CP) at a finite value of the baryon chemical potential [26,29].
In experiment, the fluctuations of the conserved charges are believed to be a promising probe of the QCD critical point. Based on universality arguments, it was predicted that the higher order cumulants of baryon/charge number fluctuations are very sensitive to the correlation length, ξ, and thus convey the information about the underlying behaviour of the critical mode [32][33][34][35]. Quantitatively, it was shown that the singular parts of the cumulants of the net baryon/charge distribution scale with the correlation length according to χ sing n ∝ ξ nβδ ν −3 where β, δ and ν are critical exponents of the three-dimensional Ising universality class. However, this sensitivity of the higher order cumulants to the critical dynamics does not come for free: they probe the tails of the probability distribution which are also susceptible to various non-critical effects including baryon number conservation [36], volume or number of wounded nucleon fluctuations [37][38][39], detector efficiency and acceptance [40][41][42][43][44][45], hadronic rescattering [46], deuteron formation [47], non-equilibrium effects [48,49], non-critical correlations between centrality trigger and the observable, etc.
Recent experimental results collected in the Beam Energy Scan, a dedicated experimental program at RHIC, demonstrated a non-monotonic dependence of the fourth order cumulant or kurtosis of net proton fluctuations on the energy of collisions [23]. Experiments also showed, that, at low energies or, equivalently, higher baryon densities, the kurtosis increases with decreasing collision energy. In the central rapidity region, where most of the measurements are performed, the baryon stopping is what makes the high baryon densities possible. Indeed at low energies, where proton-antiproton pair-production is negligible, all the observed protons originate from the target and projectile. Therefore, the dynamics of baryon stopping is yet another source of fluctuations which can distort the measurement and can be potentially responsible for the non-monotonic behavior of the kurtosis [50].
An analysis of the genuine multi-particle correlations [51] found that these correlations have a range of the order or larger than one unit of rapidity, δy ≥ 1. Typically, long-range correlations in rapidity suggest the correlations to be formed early in the evolution of the fireball, although this argument is less stringent at the lower energies. Therefore, it would be worthwhile to study the effect of initial conditions, such as so-called volume or rather participant fluctuations as well as effects of stopping.
It is the purpose of this paper to address some of these issues. We will explore the influence of wounded nucleon [52] fluctuations as well as fluctuations due to multi-collision nucleon stopping, effectively resulting in proton clusters, on the correlation functions at low energies. We show that for the integrated n-particle correlation functions, C n , the effect of wounded nucleon fluctuations is small for n > 2 and negligible for C 4 . However, for C 2 , we find that this contribution is quite significant. Additionally, we show that taking into account multi-proton clusters reproduces the right order of magnitude of the experimental values for C n .
This manuscript is organized as follows: In Sec. II, we briefly review the definition of the correlation functions previously discussed in Ref. [51]. In Sec. III, we explore the effect of the wounded nucleon fluctuations. Additional sources of correlations due to proton clustering is discussed in Sec. IV. We conclude with Sec. V.

II. NOTATION
Let us start by defining our notation. The proton multiplicity distribution measured in a given rapidity bin, ∆y, will be denoted by P (N ). The corresponding generating function H(z) is given by so that the factorial moments F k of the proton multiplicity distribution can be obtained by taking the appropriate number of derivatives of H(z): We note that the first factorial moment corresponds to the average number of particles, F 1 = N , the second factorial moment F 2 = N (N −1) gives the number of pairs, F 3 the number of triplets etc. The integrated k-particle correlation function C k , also known as the factorial cumulant κ k [41], is given by with C 1 = N . As an example, consider the two-particle rapidity correlation function this can be indeed recognized as the proper definition of the integrated two-particle correlation function. Here C 2 (y 1 , y 2 ) is the differential correlation function in rapidity, and ρ 2 (y 1, y 2 ) and ρ 1 (y) are the two-particle and the single-particle rapidity distributions, respectively. It is convenient to define the reduced correlation functions c k which, following Ref. [51], we shall refer to as couplings Finally, the proton cumulants, K n , as recently measured by the STAR Collaboration 1 , are related to the integrated correlation functions C n through where δN = N − N . 2 As was recently emphasized in Refs. [41,51], the cumulants mix the correlation functions of different orders. In Ref. [51], the correlation functions C n and the couplings c n were extracted from the preliminary STAR data [53,54]. Here we highlight the most important conclusions from this analysis. For peripheral collisions at √ s = 7.7 GeV, the couplings c k scale like 1/N k−1 ; this is consistent with the production from independent sources. At N part ∼ 200, the correlations C 3 and C 4 change signs and reach large values for the most central collisions both with large error bars. C 3 and C 4 decrease with the increasing energy and at √ s = 19.6 GeV they are consistent with zero for N part ≈ 350 whereas C 2 does not vary significantly and approximately equals 7C 2 ∼ −15 at all energies. After these preliminaries let us now discuss the effect of volume or rather participant fluctuations on the correlation functions.

III. CORRELATIONS FROM FLUCTUATIONS OF THE NUMBER OF WOUNDED NUCLEONS
Fluctuations of the number of wounded or participating nucleons N part , which in the context of fluctuation studies are often referred to as volume fluctuation [37,55,56], constitute one of the more obvious sources for correlations and fluctuations of the final protons. Also, at low energies, where baryon pair-production is negligible, practically all observed protons at mid-rapidity originate from stopped target and projectile nucleons. Since pions are abundant even at 7.7 GeV, isospin exchange reactions should be very fast and, thus, the observed protons may originate equally likely from stopped protons and neutrons [46,57]. This allows us to formulate the following minimal model which takes into account fluctuations of the number of wounded nucleons, baryon number conservation [36] and fast isospin exchange: Based on the Glauber model (see, e.g., Ref. [58]), a certain centrality selection determines the distribution of participating nucleons P (N part ). Given the number of participating nucleons N part in an event, the number of protons N observed in ∆y then follow a binomial distribution B(N ; N part ; p), where p = N / N part is the probability for any wounded nucleon to end up in the rapidity interval ∆y as a proton. Here, N is the observed mean number of protons in the rapidity interval ∆y, and N part is the average number of wounded nucleons for a given centrality selection. Obviously this model assumes that each nucleon stops independently from the other; this finds some support at higher energies [59]. Thus, the probability P (N ) to observe N protons in ∆y is given by and Given the above generating function, the correlation functions C k and couplings c k as defined in Sec. II can be 1 Note, that the STAR collaboration denotes cumulants by Cn which we reserve for the correlations.
2 For completeness, let us add that that Kn can be obtained as well from the above generating function . computed (see also the Appendix): where In absence of fluctuations in the number of participants, that is P (N part ) = δ Npart, Npart in Eq. (8), the couplings, c k , reduce to In this case, the couplings, c n , as functions of the order, n, alternate in sign and are suppressed by powers of the mean number of participants, ∼ 1/ N part n−1 . This behavior is qualitatively consistent with the analysis of the preliminary STAR data for peripheral collisions, N part < 200 [51]. As expected, the couplings, Eqs. (10-12), do not depend on the binomial probability p = N / N part , because (binomial) efficiency corrections do not alter the reduced correlations functions (see e.g. Ref. [60]).

A. Monte Carlo calculation
To calculate the contribution of N part fluctuations to the multi-particle correlations C k and the couplings c k we need to define the centrality of the collision. Following the STAR procedure, see e.g. Ref. [53], we use the tightest centrality cuts, that is, we calculate c k and C k at a given number of produced charged particles (except protons) N ch in |η| < 1 [53].
In our analysis we first calculate N part using a standard Glauber model, see, e.g. Ref. [58]. We used σ in = 31 mb for √ s = 7.7 GeV. 3 Next, for each N part we sampled charged particles from the Poisson distribution 4 with the average N ch in |η| < 1 given by Here we take a = 0.75 for √ s = 7.7 GeV. We verified that small variations in the value of a do not change our conclusions (see Sec. V for further discussion). Let us comment on (N part /2) 0.1 . It is well known that in the midrapidity region, the number of charged particles grows a bit faster than N part , see e.g. Refs. [62,63]. For example, in Au+Au collisions at RHIC energies, N ch /N part grows roughly by a factor of 1.7 between proton-proton (N part = 2) and central Au+Au collisions (N part 350). This feature can be easily understood in the wounded constituent quark (quark-diquark) model [64,65], where nucleons undergoing several collisions generate more particles than protons in p+p collisions.
After generating a sufficient number of events, for each value of N ch we calculate N part , N 2 part etc. and evaluate the couplings c k following Eqs. (10)(11)(12). Our results for c 2 , c 3 , and c 4 are shown in Fig. 1 by the solid (red) curves. The (blue) dashed curves represent the results without volume (N part ) fluctuations (no VF), see Eq. (14).
We observe that N part fluctuations lead to rather nontrivial effects in very central collisions. The coupling −c 2 changes the trend from decreasing to increasing with growing N ch (mind the minus sign in front of c 2 plotted in in Fig. 1), and both c 3 and c 4 change signs. Interestingly, similar trends were observed for c k extracted from the preliminary STAR data, see Ref. [51], although the effects for c 3 and c 4 observed in Fig. 1 are at least an order of magnitude too small for c 3 , and roughly three orders of magnitude too small for c 4 , see Eq. (7). Moreover, in our results the signs change only for very central collisions whereas in the analysis of the preliminary data this change is present at about N part ∼ 200. Finally, as we shall discuss below, the sharp wiggles observed in c 3 and c 4 disappear once one averages the couplings over a centrality region of 5%, as it is done in the STAR analysis [53].
In order to illustrate the contribution from N part fluctuations let us factor out the leading term c k ∼ 1/ N part k−1 from c 2 , c 3 and c 4 in Eqs. (10)(11)(12), that is define R n according to so that R 2 , R 3 and R 4 represent the ratios of the contributions from N part fluctuations over those arising without N part fluctuations in Eq. (14). In Fig. 2 we plot these ratios as functions of N part | N ch . Even though we apply the tightest centrality cuts, (we fix the number of charged particles with the finest possible bin width) we find corrections of 50% or more for off-central collisions and much larger modification in the most central collisions.
Finally, let us calculate the integrated correlation functions C k = N k c k ; they are directly related to the cumulants measured by STAR, see Eq. (6). To proceed we need to determine the dependence of the average number of protons, N , on N part . From the preliminary STAR data [53] we get where b = 40 for √ s = 7.7 GeV. Our conclusions are not sensitive to small variations of b and changing the exponent from 1.25 to 1. The results are presented in Fig. 3 by the solid curves. The dashed curves correspond to calculations without volume (N part ) fluctuations (no VF). The symbols represent the correlations after averaging over bins in centrality of 5%, i.e. 0 − 5%, 5 − 10% etc. Only the five most central points are shown. For less central collisions, the centrality averaging does not alter our results and points fall right on the solid lines. Clearly, the contribution originating from N part fluctuations is important for the two particle correlation, C 2 ; there is also some but less significant effect of N part fluctuations on the three particle correlation C 3 in central collisions. On the other hand, when compared to the STAR data, fluctuations of wounded nucleons are all but irrelevant for the four particle correlation, C 4 . In our model calculation, C 4 is negative for off-central collisions and it gets positive for large N part . After averaging over centrality bins, the model predicts around −0.3 for C 4 while the analysis of the preliminary STAR data gives ∼ 170. Also, as already mentioned, the strong oscillations exhibited in C 3 and C 4 at large N part disappear after averaging over centrality bins. Obviously our model of independent stopping together with baryon number conservation clearly fails to explain the preliminary STAR data, reported in Ref. [51] (see Fig. 1   Before we close this section, let us make a few more remarks. First, the results without the number of wounded nucleon fluctuations presented in this section can be verified analytically. At a fixed N part , Eq. (9) reduces to and using Eq. (3) we obtain Since p < 1 this explains the relative magnitude of the correlation functions. Next, in our analysis we assumed that each nucleon is stopped in ∆y with the same probability p. This is rather unphysical since nucleons that collide once are expected to have significantly smaller p than nucleon from the centers which collide several times. However, as long as we have independent stopping of the nucleons, individual stopping probabilities do not really change our conclusions. Suppose that each nucleon is characterized by its own stopping probability, p (i) , i = 1, ..., N part . Neglecting N part fluctuations we obtain at a given N part which obviously reduces to Eq. (20) if p i = p. Calculating C k we observe that it is enough to replace N part p n → i p n (i) in Eq. (21) and thus the signs of C k do not change. We conclude that this effect cannot lead to a large and positive C 4 as seen in the STAR data. The corollary of this section is the following. The two-particle correlations obtained in our model of independent nucleon stopping together with baryon-number conservation and fast isospin equilibration are of the same magnitude as in the preliminary STAR data. Also, the model produces a non-negligible three-particle correlation. On the other hand, the model four-particle correlation comes out to be almost three orders of magnitude smaller than in the preliminary STAR data at √ s = 7.7 GeV. 6 The large discrepancy of the four-particle correlation between our model and the preliminary data suggests that additional strong effects must be at play. In the following section, we will explore what happens if we relax the assumption of independent stopping of nucleons and allow for proton clustering.

IV. PROTON CLUSTERS
In this section go we go beyond the independent stopping assumption and consider a possible clustering of protons. This clustering can be attributed to, e.g., the collective stopping or to a first order phase transition. We note, that little is known about the stopping of nucleons at lower energies, and, therefore, this discussion is somewhat speculative and should be considered merely as a motivation to explore the stopping mechanism in more detail.
In our view, it is not at all obvious that the standard Glauber model still provides the correct stopping framework at energies as low as √ s = 7.7 GeV. Indeed, it seems plausible that a nucleon passing through a center of a nuclei may actually loose enough energy and start undergoing elastic scatterings leading to an inter-nucleus cascade. This may lead to some multi-particle correlations between stopped protons. For example, in the extreme case when two protons scatter elastically they go back-to-back and they either both end up in our rapidity bin or they both go outside. This obviously results in two-particle rapidity correlations since protons tend to come in pairs. In central collisions more protons can be involved in this multi-collision (quasi-) elastic stopping, and thus, higher-order proton multiplets are not a priori excluded. We note that this mechanism is expected to turn on at relatively low collision energy and in central collisions, where protons collide several times, and thus may lose sufficient energy to engage in subsequent elastic scatterings. In other words, the stopped protons from the centers effectively form proton clusters. Of course cluster formation may also arise from a first order phase transition (see, e.g., [66][67][68]) or final state interactions and it will require more detailed study to disentangle the various possibilities.
We explored various scenarios where protons come in pairs, triplets and higher-order multiplets. We found that it is not an easy task to reproduce C 2 < 0, C 3 < 0 and C 4 > 0 with 6C 3 and C 4 of the order of 100 for √ s = 7.7 GeV, as seen in Ref. [51]. To discuss this in more detail, let us focus on central √ s = 7.7 GeV collisions, were the signal is strongest. In this case N 40, 7C 2 −15, 6C 3 −60, and C 4 170. Before we discuss two scenarios which give results in the right ballpark, let us put the difficulty of the task in perspective. Consider a system of clusters distributed according to a Poisson distribution which decay into a fixed number of protons, m. In this case, the generating function is given by where N cl is the average number of clusters. The resulting integrated correlation functions C k are given by C k = N cl m!/(m − k)! so that for four-particle clusters (m = 4) we have C 4 = 24 N cl . Thus, in order to get C 4 170 we need on average 6−8 clusters, about 25 of the 40 observed protons would have to originate from clusters, a rather large fraction indeed. Of course the same cluster model would also lead to positive two-and three-particle correlations, contrary to what is seen in the data. Let us now turn to discuss two scenarios which lead to a correlation structure in qualitative agreement with the preliminary STAR data. Consider central Au+Au collisions and suppose that nucleons at the surface (and in general those who do not scatter sufficiently often) are stopped according to a simple string picture, and, thus, follow the previously discussed model of independent binomial stopping. Next, suppose that nucleons from the centers of target and projectile, which scatter several times, loose most of the energy and engage in elastic scattering and, thus, fall into ∆y in pairs. In this case the generating function is given by where M is the number of pairs. Let us clarify the above formula. N part − 2M nucleons either stop in ∆y with the probability p 1 or they do not with the probability 1 − p 1 , and the generating function for each nucleon is simply given by (1 − p 1 )z 0 + p 1 z 1 . In the second term we have M pairs which either fall into ∆y with probability p 2 or they do not with the probability 1 − p 2 . In this case we either get 0 or 2 protons (from each pair) and the pair generating function is given by (1 − p 2 )z 0 + p 2 z 2 . We also expect that p 2 is much larger than p 1 since pairs loose energy due to scatterings and are likely to fall into the midrapidity bin of interest ∆y (for STAR ∆y ∈ [−0.5, 0.5]). The multi-particle correlations C k are given by the appropriate derivatives (at z = 1) of resulting in Taking N part = 350, N = 0.12N part = 42 and M = 8 (8 pairs of protons) we obtain the relation between p 1 and p 2 . In Fig. 4 (left) we plot 7C 2 , 6C 3 and C 4 as a function of p 2 . We observe that for p 2 > 0.5 both C 3 and C 4 have the right signs and can reach substantial values. The right panel of Fig. 4 shows the results of an analogues calculation where protons come in quartets instead of pairs. In this case and N = p 1 (N part − 4M ) + 4p 4 M, where M in this case is the number of proton quartets. In this calculation we use M = 4 so that the number of correlated protons is the same as in the previous case. We observe that the signal for p 4 > 0.7 is much larger, and all C n agree qualitatively with the STAR data. We have also verified that the signal increases even further if protons are clustered in even larger multiplets.

V. DISCUSSION AND CONCLUSIONS
Let first summarize the main findings of this paper.
• We have studied the proton correlations at low energies where proton-antiproton pair production can be neglected. To this end we developed a minimal model which is based on independent stopping of nucleons, baryon number conservation and fast isospin-exchange. We find that this model qualitatively reproduces the twoproton correlations seen in the preliminary STAR data, while it underpredicts the magnitude of the four-proton correlations by almost three orders of magnitude. • Fluctuations of the number of participating nucleons, though significant even for the tightest centrality cuts, are nowhere near large enough to explain the observed four-proton correlations.
• The observed large four-particle correlations as well as the signs and rough magnitudes of the two-and threeparticle correlations can be reproduced if one assumes that about 40% of the observed protons originate proton quartets. Given that at lower energies the incoming nucleons in the center are likely to loose most of their energy, such a scenario may be not as far fetched as one would think initially.
Before we conclude, we want to mention that: • Although our minimal model fails to reproduce the measured data, we consider it as a better baseline for low energies than the Poisson distribution, which is typically used. It respects baryon number conservation, and contains the effect of participant fluctuations.
• Although we achieved a qualitative agreement with the observed two-particle correlations, the details do matter. This is demonstrated in Fig. 5. There we vary the parameter a of Eq. (15) from its standard value of a = 0.75 to a = 0.6 and a = 0.9 and plot the resulting couplings c 2 , c 3 , and c 4 . While the three-and four-particle couplings are hardly changed, the two-particle couplings exhibit strong sensitivity on the dependence of the number of charged particles on N part . This dependence can only and should be further constrained by the charged particle distribution, which, however, is not yet publicly available for the STAR measurements.
• Collective stopping or a possible first-order phase transition may lead to clustering of protons at midrapidity. As we pointed out above, by taking this clustering into account it is possible to qualitatively understand the STAR data for the multi-particle correlation functions.
• It would be interesting to measure mixed correlations of fourth order. For example a correlation C 3,1 between three protons and one produced particle, such as anti-proton or pion may be useful (see the appendix of Ref. [51] for the relevant formulas). If these mixed correlations, or rather couplings c 3,1 , to avoid trivial sensitivity to the number of particles, are as large as the four-proton coupling c 4 it would rule out collective stopping and provide evidence for a possible first-order phase transition, as speculated in Sec. IV.
• In this paper we have concentrated on the lowest energy, √ s = 7.7 GeV. Our basic model equally applies to higher energies as long as one can ignore pair production, which is probably still the case at √ s = 19.6 GeV. There, the observed correlations are more or less in agreement with our basic model, subject to the aforementioned uncertainty related to the distribution of the number of charged particles.
In conclusion, the observed large four-proton correlation exhibited in the preliminary STAR data at √ s = 7.7 GeV cannot be explained by a combination of independent baryon stopping, participant fluctuations, and baryon-number conservation. However, more speculative, though not necessarily unrealistic, scenarios involving the "collective" stopping of multiple baryons, may be able to explain the large part of the observed correlation structure. Therefore, it is important to test these ideas by, e.g., measuring mixed correlations. However, given the present analysis and the vary large positive value for the four-proton correlation, C 4 , observed in the preliminary STAR data, it is very difficult to imagine a scenario which explains the data without employing some kind of cluster formation. If these clusters originate from collective stopping, a possible first order phase transition [66][67][68][69], or some final state effect, remains to be seen, however.
Appendix A: General independent source model Here we generalize our Eqs. (10,11,12) for any independent sources of protons (see also [70]). Suppose we have v independent sources of protons distributed according to G(v) and each source produces protons according to p(n i ). The distribution of measured protons in a given rapidity bin ∆y is given by n1,...,nv p(n 1 )p(n 2 ) · · · p(n v )δ n1+...+nv−N . (A1) We note that this formula is quite general and the only assumption we make is that each source decides on its own how many particles it produces. The generating function reads where H 1 (z) is the generating function for a single source.
To calculate correlation functions we use where C 1 (z) characterizes correlations from a single source. We have and the generating function for a single source is given by Substituting the above values to Eqs. (A5,A6,A7) and putting v = N part , we obtain Eqs. (10,11,12).