Benchmark values for the net proton number fluctuations
Boris Tomasik Ivan Melo Lukas Laffers Marcus Bleicher




§ INTRODUCTION


Statistical physics arguments lead to the conclusions about strongly interacting matter, 
that in the vicinity of the critical point of the phase diagram the susceptibilities 
of higher orders should strongly depend on temperature and the baryon chemical potential μ_B 
<cit.>.
The susceptibilities with respect to μ_B represent themselves in the cumulants of the net-baryon number 
distribution <cit.>. 
In principle, they can be even calculated from the first principles by simulating the QCD partition function on the 
lattice <cit.>. 
These arguments and calculations are intrinsically connected with the grand-canonical ensemble of systems. 

Unfortunately, the baryon number fluctuations, as they are elegantly treated in the grand-canonical formalism, cannot be directly measured 
in an experiment. Three major reasons make this impossible: 


  * In real collisions, baryon number is always conserved. Hence, if the whole phase-space would be covered by detectors that would 
register all baryons and antibaryons with perfect efficiency, there would be no fluctuations. On the other hand,
if the detectors cover only a part of the phase-space, then fluctuations may arise just from the binomial distribution of the 
appearance of a particle in the acceptance window <cit.>. 

  * It is experimentally impossible to investigate a set of heavy-ion 
collisions with always exactly the same number of nucleons participating 
in the collisions. This brings in another source of baryon number fluctuations. It has to be taken into account also in theoretical 
simulations which are compared to experimental data. Another possibility is to filter out such fluctuations from the experimental data 
my means of unfolding <cit.>. 

  * Not all baryons (and antibaryons) can be detected. A typical detector does not register neutrons. Therefore, the net-proton  number  
is usually measured as a good proxy for the net baryon number. The argument in favour of this proxy is usually based on the high 
rate of the isospin-changing reactions that can change protons into neutrons and vice versa <cit.>.   Note, however, 
that such arguments may be of limited validity at the lowest collision energies within the RHIC Beam Energy Scan programme. 


The less-than-perfect acceptance of the detectors also acts as a source of fluctuations. This problem has been addressed in the 
literature and we will not analyse it here  <cit.>. 
Instead, we assume that the detector efficiency is 100%.

In this study we investigate the influence of the effects mentioned above on the net-baryon and the net-proton number 
fluctuations. Two new points are introduced in our treatment. 


  * We simulate the participating baryons and the produced BB̅ pairs with rapidity distributions that depend on the collision energy. 
This, combined with the rapidity window fixed for all collision energies (since it is really fixed by the actual detector), leads to the
collision energy dependence of the fluctuations. The rapidity distributions are constructed with an eye on the experimental data <cit.>. 
The distribution of the  wounded nucleons is different from the distribution of the produced BB̅ pairs. 

  * Since for lower collision energies the isospin randomisation may become ineffective, we also explore the possibility 
that the  wounded nucleons remember their isospin. 


We study the fluctuations as functions of the size of the acceptance window in rapidity but also as a function 
of the position of the acceptance window. The latter is relevant as a baseline study for the recent proposal 
to look for the position of the critical point by measuring fluctuations at slightly forward rapidities, where the 
net baryon density may be higher than at midrapidity <cit.>. This is doable for us owing to the implemented rapidity 
distributions of the wounded nucleons.





§ THE MODEL


Our results will be based on Monte Carlo simulations that will assign the baryons and antibaryons their rapidity. 
First of all, we respect the baryon number conservation. Only protons, neutrons, and their antiparticles 
are included into our simulations. From the colliding gold nuclei only a part of the incoming nucleons 
participate in the collision. 

The number of wounded nucleons is determined with the help of GLISSANDO <cit.>. We first create a sample 
of minimum bias collisions for each of the investigated collision energies. 
Then we select centrality classes based
on the relative deposited strength (RDS) which should be proportional to the produced multiplicity

    M ∝1-α/2 N_w + α N_bin   ,

where the parameter α increases with the collision energy. It has been set by fitting the centrality dependence of the 
multiplicity in Au+Au collisions at √(s_NN) = 19.6 GeV and 200 GeV and interpolating for other energies 
with the logarithmic dependence

    α(√(s_NN)) = α_0 + α_1ln√(s_NN)   .

This method is  closest to what is done in real experiment. 

One of the features of our model, which can be turned on and off, is the isospin memory of wounded nucleons. 
This may be (partially) present in  collisions at the lowest energies. The isospin of a wounded nucleon is set in such a way 
that one never can obtain more protons (or neutrons) than there are in the incoming nucleus. This means that if 
N_p wounded protons and N_n wounded neutrons have been generated, then the probability that the next 
wounded nucleon will be a proton is (Z-N_p)/(A-N_p-N_n), where A and Z refer to gold nucleus. The number of 
wounded protons is then distributed according to hypergeometric distribution. 

In addition to the wounded nucleons there are also produced pairs of nucleons and antinucleons. The number of pairs fluctuates 
according to Poisson distribution, with its mean proportional to the number of wounded nucleons

    μ_BB̅= dN_p̅/dy  y_m N_w/⟨ N_w ⟩   ,

where dN_p̅/dy is the rapidity density of antiprotons measured at given energy and centrality <cit.>, 
y_m is the width of 
the whole rapidity distribution of the produced BB̅ pairs, N_w is the number of wounded nucleons in the particular 
event, and ⟨ N_w ⟩ is the mean number of wounded nucleons at given collision energy and centrality, as determined 
by GLISSANDO. 
The isospin is assigned with equal probabilities for (anti)protons and (anti)neutrons. 
Antinucleons are, of course, only present among the produced BB̅ pairs. 

Not the whole rapidity region is covered by the detectors, thus only a part of the (anti)nucleons with rapidities in the experimentally accepted
window is assumed to be registered. 

It is a crucial feature of this model, that the relative abundance of wounded and produced nucleons depends on the collision energy, centrality,
and the acceptance window in rapidity.

Finally, let us specify the rapidity distributions. Wounded nucleons partially remember their original rapidity, so we parametrize
their rapidity distribution with a double Gaussian function

    dN_w/dy = N_w/2√(2πσ_y^2){exp ( - (y + y_m)^2/2σ_y^2 )
    +
    exp ( - (y - y_m)^2/2σ_y^2 )
    }   .

We set the width σ_y = 0.8, and N_w is the total number of wounded nucleons in the given event. 
The shifts y_m are determined from data. We use the fact that wounded protons make up the difference between 
observed protons and antiprotons. Thus the number of the observed net protons in a given rapidity window (-y_b,y_b)
can be determined

    N_p-p̅ = Z/A∫_-y_b^y_bdN_w/dy  dy   .

We use data for y_b = 0.25 measured at different collision energies by the STAR collaboration <cit.>. 

For the produced nucleon-antinucleon pairs we assume that they are produced mainly around midrapidity. To soften sharp edges 
of their distribution we use the Woods-Saxon distribution

    dN_BB̅/dy = N_BB̅C/1 + exp ( |y| - y_m/a )   ,

where C is the proper normalisation constant C = (2a ln(e^y_m/a + 1))^-1, and we choose a = σ_y /10.
The number of pairs N_BB̅ is determined from data: note that this is the only source of antiprotons in our model. 
Thus we set the parameters so that we can reproduce the measured number of antiprotons N_p̅ in the interval 
|y| < y_b = 0.25 <cit.>

    N_p̅ = 1/2∫_-y_b^y_bdN_BB̅/dy  dy   ,

where the factor 1/2 in front of the integral stands for taking only antiprotons and no antineutrons. 

The pairs of nucleons and antinucleons are generated so that their mutual distances in rapidity are distributed exponentially with 
the correlation length 1. The correlation length may have influence on the fluctuations <cit.> and we plan to investigate this effect 
in the future. 

With this model we determine the central moments μ_i up to the fourth order and the volume-independent ratios of cumulants
χ_i

    Sσ    =    χ_3/χ_2 = μ_3/μ_2
    κσ    =    χ_4/χ_2 = μ_4/μ_2 - 3 μ_2 
    κσ^4/n̅    =    χ_4/χ_1 = μ_4 - 3 μ_2^2/μ_1   ,

where we have used σ for the width, S for the skewness, and κ for the kurtosis of the multiplicity distribution. 






§ A WARM-UP EXERCISE: BARYON NUMBER CONSERVATION AND MODIFICATIONS


Before we embark on more realistic simulation, let us test the model and present the mere effect of baryon number 
conservation and limited acceptance. 
We inspect a sample of 5.10^7 events with fixed number of wounded nucleons N_w = 338 and other 
parameters set to reproduce Au+Au collisions at √(s_NN) = 19.6 GeV: y_m = 1.019 and 
N_BB̅ = 16.946.

In Figure <ref> 

	

we plot the central moments and volume independent cumulant ratios of the net baryon number 
distribution, as functions of the width of the acceptance window Δ y around midrapidity. 
We show the central moments since they are the primary quantities which could be determined experimentally
and the other ratios are based on them. Thus one can relate the behaviour of the Δ y dependences
of the scaled skewness and the scaled kurtosis to the behaviour of the central moments. 
As Δ y increases, the mean value μ_1 saturates at N_w while the higher moments go to 
0 as there is no room for fluctuations due to baryon number conservation. The dependences are actually 
described by the  binomial distribution where the probability p equals μ_1(Δ y)/N_w. This is 
demonstrated in the plots for the scaled skewness Sσ and the scaled kurtosis κσ^2. The simulated points
are compared to the curves obtained for the binomial distribution

    Sσ    =     1 - 2p = 1 - 2μ_1(Δ y)/N_w
    κσ^2     =     1 - 6p + 6 p^2 = 1 - 6 μ_1(Δ y)/N_w + 6 ( μ_1(Δ y)/N_w )^2   .

The observed dependences result just from the rapidity distribution of baryons and antibaryons 
combined with the binomial distribution of the number of accepted nucleons. 

The third central moment and consequently the skewness is  positive for small Δ y and then turns negative. This means, that 
in small rapidity windows numbers larger than the mean are rather populated. As the rapidity window gets larger, 
the total baryon number limits the possible observed numbers and they are more often populated below the mean.

Let us next look at how the results are modified, if now instead of baryons one measures the fluctuations of net protons. 

	

They are shown in Fig. <ref>.
The rough features of the baryon number distribution are inherited here. However, 
just by looking at the net protons—whose number is not conserved—there are qualitative differences  from 
the net baryons and the binomial distribution is no longer applicable. This is due to two reasons: firstly, the produced nucleon-antinucleon 
pairs populate protons and neutrons  with equal probability. Secondly, also the wounded nucleons may change isospin. The result 
is that now the total number of net protons is not fixed and fluctuates from event to event. This is seen in the behaviour 
of the higher central moments which saturate at a non-vanishing values for large rapidity windows. The asymptotic values of 
Sσ and κσ^2 are modified, as well. 



	

The next complication comes from the fluctuation of the number of participating nucleons. Centrality classes in an experiment 
are usually determined from the produced multiplicity. That usually scales with the relative deposited strength as defined in 
eq. (<ref>) which can be evaluated in GLISSANDO. 
Within the 0–5% percentile with highest RDS among all events, which corresponds to the experimental class of 
0–5% centrality, the number of wounded nucleons fluctuates. Figure <ref> shows how this influences the net-proton 
number fluctuations. In general, all central moments characterising fluctuations grow considerably with respect to the 
simulation with fixed N_w. 
Nevertheless, some qualitative features of their dependence on Δ y remain, e.g., the third moment flips the sign 
at the same Δ y as it did for fixed N_w.
These features are then reflected also in Sσ and κσ^2. 








§ RESULTS



Figures <ref>-<ref> demonstrate the influence of the three basic effects on the observed 
fluctuations: baryon number conservation (Fig. <ref>), measurement of only protons (Fig. <ref>), 
and fluctuations of the interacting volume (Fig. <ref>). From Fig. <ref> one would conclude that the latter 
influence is tremendous. However, no real detector in the current experiments covers 10 units of rapidity. Therefore, 
in Figure <ref> we zoom into a narrower Δ y-interval of 2 units, which may realistically be covered
by the central detector, e.g. in the STAR experiment at RHIC. 

	

Here, the only large difference  among the central moments shows up for the third moment: results from Glauber MC simulations
lead to larger μ_3.
This is then translated in about a constant increase of Sσ by 0.25 in the whole interval of Δ y. Up to 
Δ y = 1, the values of κσ^2 for fluctuating and for fixed N_w are practically identical. 
However, they start to depart from each other strongly beyond this limit. 

An interesting suggestion <cit.> has been presented also at this conference <cit.> that the critical point might be sought 
by inspecting the net-proton number fluctuations at different rapidities as an alternative to different collision energies. 
Our model includes the feature, that rapidity distributions of protons and antiprotons are different. Hence, we can make 
predictions for the rapidity dependence of the net-proton number fluctuations in absence of any critical behaviour.  
This is shown in Figure <ref>.

	

Simulations were done for a set of parameters that correspond to √(s_NN) = 19.6 GeV. Results have been extracted 
for a rapidity window of Δ y = 0.5, which is put in different positions along the rapidity distribution. Figure <ref>
shows the dependences on the central value of the rapidity window. We compare the results obtained for fixed N_w = 338 
(2 ×10^7 events in the sample) with events from 0–5% central collisions according to GLISSANDO Glauber MC (1.2× 10^6 events).
The central moments calculated for fixed N_w generally decrease as we go towards higher y. In contrast to that, the curves
for μ_1, μ_2, and μ_4 simulated with Glauber MC show peaks around y = y_m, i.e., where the antiprotons 
die out. Note that the antiproton rapidity distribution is cut more sharply than that of wounded protons. This results in a dip in 
both Sσ and κσ^2as functions of y around y =1.

Since y_m increases with the collision energy, we would also expect that the position of the dip 
would move to higher  values of y_m. 

	

This is confirmed in Figure <ref>, where we show Sσ, κσ^2, and χ_4/χ_1 as functions 
of y for three collision energies. 
The parameters used in the simulations for all energies are summarised in Table <ref>.



In the Figure we only show results from Glauber MC version of our model. Unfortunately, the 
statistics of 1.2× 10^6 events is not sufficient to avoid large statistical errors in the curves for the fourth order. Nevertheless,
one can see that this simple baseline model predicts rather strong dependence on y of the fourth order cumulant ratios which 
for √(s_NN) = 200 GeV start
around 1 and fall down all the way to 0 at y = 1.2. The decrease is less severe when the collision 
energy decreases. 

In order to present the centrality dependence of the fluctuations we turn back to √(s_NN) = 19.6 GeV. 

	

Selected results 
are shown in Fig. <ref>: Sσ and κσ^2 as functions of the number of wounded nucleons N_w. 
The number N_w is either fixed or we use for plotting the mean number from GLISSANDO for the selected centrality class. 
With these results we also show the influence of wounded nucleon isospin randomisation or isospin memory.  
We investigated the influence of this effect on the simulations with fixed N_w. The centrality dependence appears flat 
with randomised isospins. The effect of isospin memory is most pronounced in most central collisions since it is here that the 
limit on incoming proton number kicks in most strongly. As a result, the third and the fourth cumulant are slightly lowered in central 
collisions. The isospin memory has also been assumed in the simulations with the Glauber MC model. Although the statistics of 
5× 10^5 events for each centrality still leads to large statistical uncertainties, we see that these fluctuations increase 
Sσ by about 0.1 and κσ^2 by about 0.2. 


Finally, we investigated the dependence of the moments on the collision energy in Figure <ref>.

	

(Note that model parameters for all energies have been summarised in Table <ref>.)
We show results for 0–5% central  collisions. Since we are showing results for the central rapidity window
of Δ y = 0.5, the baryon-antibaryon asymmetry decreases as the collision energy goes up. The relative 
contribution of produced BB̅ pairs is increasing. This is demonstrated by √(s_NN) dependence of 
μ_1. In general, all higher central moments decrease, as well, if we look at the results obtained for fluctuating 
N_w. For μ_2 and μ_4 we see that the curves for fixed N_w (i.e. which do not include all baseline effects)
increase together with √(s_NN). This behaviour is reflected in Sσ which starts  below 0.8 at low collision 
energies and approaches 1 as √(s_NN)→ 200 GeV for Glauber MC simulations, but is strongly decreased in 
fixed N_w results.
The scaled kurtosis κσ^2, on the other hand, does not change when N_w fluctuations are turned on 
and off, since μ_2 and μ_4 behave similarly, here. The scaled kurtosis starts at very low values for 
√(s_NN)=7.7 GeV and grows together with √(s_NN). Unfortunately, this is exactly opposite to the 
observed experimental result <cit.> which shows a huge increase of κσ^2 as 
√(s_NN)→ 7.7 GeV. Hence, our result just highlights the uncommonness of the measured value.





§ CONCLUSIONS


We have set up a benchmark model which includes many non-critical effects that influence the fluctuations of net-proton number. 
Generally, it is always necessary to look at such trivial predictions when interpreting the measured data. 

In particular, our model allowed to look at the dependence of the moments of the net-proton number distribution on the rapidity 
of the centre of the acceptance window. Our predictions, presented in Figures <ref> and <ref>, are thus complementary 
to the proposal to use the rapidity dependence for looking for the critical point of the phase diagram <cit.>.

We have also looked at the effect of the wounded nucleons remembering their original isospin. Such a mechanism may be present
in nuclear collisions at √(s_NN) of the order of a few GeV, because there the isospin randomisation may become ineffective.

Nevertheless, the presented baseline effects are not sufficient to explain all observed features of the data; in particular the huge 
enhancement of κσ^2 towards √(s_NN)=7.7 GeV is not reproduced by this model. 
One has to, however, always keep them in mind when interpreting any measured result. 






  §.§.§ Acknowledgements

This work was supported by the grant 17-04505S of the Czech Science Foundation (GAČR)
and the grant VEGA  1/0348/18 (Slovakia).
International collaboration was supported by collaboration grant in framework of the German-Slovak PPP programme and 
the COST Action CA15213 THOR.



99

Hatta:2003wn
  Y. Hatta and M. A. Stephanov,
  
  Phys. Rev. Lett.  91 (2003) 102003
   Erratum: [Phys. Rev. Lett.  91 (2003) 129901]
  
  
  
  

Athanasiou:2010kw
  C. Athanasiou, K. Rajagopal and M. Stephanov,
  
  Phys. Rev. D 82 (2010) 074008
  
  
  
  

Gavai:2010zn
  R. V. Gavai and S. Gupta,
  
  Phys. Lett. B 696 (2011) 459.

 
  
  

Stephanov:2011pb
  M. A. Stephanov,
  
  Phys. Rev. Lett.  107 (2011) 052301.

 
  
  

Gupta:2011wh
  S. Gupta, X. Luo, B. Mohanty, H. G. Ritter and N. Xu,
  
  Science 332 (2011) 1525

 
  
  

Stephanov:2008qz
  M. A. Stephanov,
  
  Phys. Rev. Lett.  102 (2009) 032301.

 
  
  

Bazavov:2017tot
  A. Bazavov et al. [HotQCD Collaboration],
  
  Phys. Rev. D 96 (2017) no.7,  074510.

 
  
  

Borsanyi:2018grb
  S. Borsányi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K. Szabo, A. Pasztor, I. Portillo and C. Ratti,
  
  JHEP 1810 (2018) 205
  
  
  
  

Braun-Munzinger:2016yjz
  P. Braun-Munzinger, A. Rustamov and J. Stachel,
  
  Nucl. Phys. A 960 (2017) 114
  
  
  
  
  
  Bzdak:2017ltv
  A. Bzdak and V. Koch,
  
  Phys. Rev. C 96 (2017) no.5,  054905.

 
  
  
  
 Luo:2013bmi
  X. Luo, J. Xu, B. Mohanty and N. Xu,
  
  J. Phys. G 40 (2013) 105104.
 
 
  
  


Sugiura:2019toh
  T. Sugiura, T. Nonaka and S. Esumi,
  
  arXiv:1903.02314 [nucl-th].
  

KA
  M. Kitazawa and M. Asakawa,
  
  Phys. Rev. C 85 (2012) 021901.
  
  
  
  

KA2
M. Kitazawa and M. Asakawa,
  
  Phys. Rev. C 86 (2012) 024904,
   Erratum: [Phys. Rev. C 86 (2012) 069902].
  
  
  
  

Bzdak:2013pha
  A. Bzdak and V. Koch,
  
  Phys. Rev. C 91 (2015) no.2,  027901
  doi:10.1103/PhysRevC.91.027901
  [arXiv:1312.4574 [nucl-th]].
  
  
  
  Luo:2014rea
  X. Luo,
  
  Phys. Rev. C 91 (2015) no.3,  034907
   Erratum: [Phys. Rev. C 94 (2016) no.5,  059901]
  doi:10.1103/PhysRevC.91.034907, 10.1103/PhysRevC.94.059901
  [arXiv:1410.3914 [physics.data-an]].
  
  
  
  Kitazawa:2016awu
  M. Kitazawa,
  
  Phys. Rev. C 93 (2016) no.4,  044911
  doi:10.1103/PhysRevC.93.044911
  [arXiv:1602.01234 [nucl-th]].
  
  
  
  Bzdak:2016qdc
  A. Bzdak, R. Holzmann and V. Koch,
  
  Phys. Rev. C 94 (2016) no.6,  064907
  doi:10.1103/PhysRevC.94.064907
  [arXiv:1603.09057 [nucl-th]].
  
  

Brewer
J. Brewer, S. Mukherjee, K. Rajagopal and Y. Yin,
  
  Phys. Rev. C 98 (2018) no.6,  061901.
  
  
  
  
  
glis
  P. Bożek, W. Broniowski, M. Rybczynski and G. Stefanek,
  
  arXiv:1901.04484 [nucl-th].
  
  
glis2
  M. Rybczynski, G. Stefanek, W. Broniowski and P. Bożek,
  
  Comput. Phys. Commun.  185 (2014) 1759.


  
  
  

glis1
  W. Broniowski, M. Rybczynski and P. Bożek,
  
  Comput. Phys. Commun.  180 (2009) 69.


  
  


star2009
B.I. Abelev et al. [STAR collaboration], Phys. Rev. C 79 (2009) 034909

starbulk 
L. Adamczyk et al. [STAR collaboration], Phys. Rev. C 96 (2017) 044904

pruneau
C. A. Pruneau,
  
  arXiv:1903.04591 [nucl-th].
  

jasmine
J. Brewer, talk at this conference.

STARfluc X. Luo, Nucl. Phys. A 956 (2016) 75.



