A Monte Carlo study of multiplicity fluctuations in proton-proton collisions at $\sqrt{s}=$~7~TeV

With large volumes of data available at LHC, it has possible to study the multiplicity distributions. It is interesting as well to check how well event generators can describes the properties and the behavior of multi-particle production processes. In this paper, we analyse the oscillatory behavior of modified combinants in proton-proton collisions at centre of mass energy of 7 TeV.


I. INTRODUCTION
Multiplicity distributions (MDs) of charged particles produced in high-energy nuclear collisions have been extensively studied in the field of multi-particle production.The determination of multiplicity distribution is among the initial observations in new high-energy experiments, primarily because it is relatively easy to obtain such information.Furthermore, MDs provide valuable insights into the underlying production processes.Since perturbative QCD fails to fully explain the observed MDs, a range of phenomenological approaches have been employed.These approaches include dynamical methods like colored string interactions [1] and the dual-parton model [2], as well as geometrical approaches [3,4] leading to the fireball model [5], and stochastic approaches [6][7][8] that model high-energy collisions as branchings [6][7][8] or clans [9].
Charged particle multiplicity distribution, P (N ), is usually fitted with a single negative binomial distribution (NBD) [10]: NBD has two free parameters: p describing probability of particle emission and parameter k ≥ 1 influencing shape of the distribution.Nevertheless, as the energy and the number of charged secondaries, denoted as N , increase, the negative binomial distribution tends to deviate from the observed data for large values of N , as discussed in [11].In these cases, alternative approaches are adopted, including combinations of two [12,13], three [14], or multi-component NBDs [15], or even different forms of P (N ) distributions [1,10,[16][17][18].However, it should be noted that such adjustments primarily improve the agreement for large N , while the ratio R = data/fit exhibits significant deviations from unity at small N across all fitting scenarios [11,19].
Such a observation suggests that there is additional information in the measured multiplicity distribution that * maciej.rybczynski@ujk.edu.pl is not covered by the following recurrence relation: (2) Three commonly encountered forms of P (N ) resulting from the recurrence relation ( 2) are as follows: the binomial distribution, where α = Kp/(1 − p) and β = −α/K; the Poisson distribution, where α = λ and β = 0; and the negative binomial distribution, where α = kp and β = α/k.Here, the parameter p again represents the probability of particle emission.In our previous work [11], we introduced a more generalized form of the recurrence relation that is applicable in counting statistics when considering multiplication effects in point processes [20].Unlike Eq. ( 2), this new relation connects all multiplicities using coefficients C j , which determine the corresponding P (N ) in the following manner: The modified combinants, C j , can be obtained from experimental data and exhibit a pronounced oscillatory pattern.This behavior is not only observed in proton-proton collisions, as discussed in previous works such as [11,[21][22][23], but has also been recently demonstrated in e + e − annihilation processes, as shown in [24].These oscillations suggest the presence of additional information regarding the multi-particle production process that remains undisclosed.The periodic nature of the oscillations observed in the modified combinants derived from experimental data is particularly indicative in this regard.Nevertheless the probability that such oscillations are statistically insignificant is very small (∼10 −16 , see Ref. [21] for more details) the sensitivity to experimental procedures are still under debate.
The aim of this paper is to show that the observed oscillations have a physical origin and are not the result of experimental procedures.We focus on the analysis of the Monte Carlo simulated events and comparison with existing experimental data.The paper is organized as follows.In Sec.II we discuss the methodology of event generation and analysis of model data.In Sec.III we provide a concise description of the results we obtained for proton-proton interactions.Finally, in Sec.IV we made several comments referring to the oscillatory behavior of the higher-order moments of multiplicity distributions observed both in experimental data and the models.

II. EVENT GENERATION AND ANALYSIS METHODOLOGY
PYTHIA [25] is a widely used Monte Carlo event generator program designed to generate events in highenergy physics.It serves as a tool for simulating collisions at high-energies involving elementary particles like e + , e − , protons, anti-protons, as well as heavy-ions, and various combinations thereof.The program encompasses a wide range of physics aspects, such as total and partial cross sections, interactions at both hard and soft scales, parton distributions, initial-and final-state parton showers, matching and merging of matrix elements with showers, multi-parton interactions, as well as processes related to hadronization/fragmentation and particle decays.
The Energy Parton Off-shell Splitting (EPOS) transport model [26] is a Monte Carlo event generator program designed for simulation high-energy particle collisions.It provides a framework to study various aspects of particle interactions and the resulting hadron production in both nucleus-nucleus and proton-proton collisions.EPOS considers each nucleus-nucleus or proton-proton collision as a collection of many elementary collisions happening simultaneously.These collisions involve the exchange of "parton ladders", which represent the evolution of partons from the projectile and target sides towards the central region (small x).The evolution of partons in EPOS is governed by an evolution equation, typically based on the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) formalism.The intermediate gluons in the parton ladder are treated as kink singularities within the framework of relativistic strings.These strings represent flux tubes that connect the interacting partons.These flux tubes eventually decays by producing quarkantiquark pairs, giving rise to fragments that are identified as hadrons [27].
The Ultra-relativistic Quantum Molecular Dynamics (UrQMD) model, as outlined in [28,29], is a microscopic transport model that employs the covariant propagation of hadrons along classical trajectories, coupled with stochastic binary scatterings, color string formation, and resonance decay.It provides a comprehensive framework to study the dynamics and interactions of particles across a wide energy range.It operates as a Monte Carlo solution to a complex system of coupled partial integro-differential equations, which describe the time evolution of phase space densities for various particle species.In the UrQMD model, baryon-baryon colli-sions at lower energies consider the exchange of electric and baryonic charge, strangeness, and four-momentum in the t-channel.On the other hand, meson-baryon and meson-meson interactions are treated through the formation and subsequent decay of resonances, following the s-channel reaction mechanism.At higher energies, a vast array of particle species can be generated, and the model accounts for subsequent re-scatterings.The UrQMD model allows generate all types of particles in hadron-hadron collisions and enables their further interaction with one another [27].For the analysis described in this paper, we used UrQMD set to LHC mode.It means that no hydrodynamics functions were activated.This compilation mode of the UrQMD model, in short words, is prepared to do calculations over high values of multiplicity.No essential changes to the model parameters were introduced.
In this study we have used PYTHIA 8.308 [25], EPOS LHC [26] and UrQMD 3.4 set to LHC mode [28] to generate proton-proton interactions at √ s = 7 TeV in accordance to data on charged particles multiplicity distributions obtained by the ALICE experiment at CERN LHC [30].In PYTHIA simulation we have implemented the inelastic component of the total crosssection for soft-QCD processes with the parameter Soft-QCD:inelastic=on.The remaining set of PYTHIA parameters we left with its default values.In the case of EPOS LHC and UrQMD we used default values of the parameters.To match with the experimental conditions, charged particle multiplicities have been chosen in the trigger conditions and acceptance of the ALICE detector, defined in [30].Namely, the generated events of collisions (EOCs) were divided into two classes: inelastic (INEL) class and non-single diffractive (NSD) class.The generated EOC belongs to INEL class if there is at least one charged particle in either the −3.7 < η < −1.7, |η| < 2.98, or 2.8 < η < 5.1 pseudorapidity interval corresponding to the acceptances of the V0-C, SPD, and V0-A ALICE sub-detectors, respectively.The NSD class requires charged particles to the detected in both −3.7 < η < −1.7 and 2.8 < η < 5.1 pseudorapidity intervals [30].

III. RESULTS
Multiplicity distributions P (N ) of charged particles in simulated events and the modified combinants C j that results from them are shown in Fig. 1 for INEL class and in Fig. 2 for NSD class.Modified combinants in all models exhibit oscillating behavior, however their amplitudes and periods of oscillations are different.Growth of amplitudes with rank j can be described as N C j ∼ a j with a = 1.042, 0. The corresponding modified combinants Cj emerging from them.PYTHIA 8 [25] with Soft-QCD:inelastic processes (solid lines), EPOS LHC [26] (dotted lines) and UrQMD 3.4 set to LHC mode [28] (dashed lines).For all models the applied kinematic cuts as described in the ALICE paper [30] (INEL class).
oscillations are smaller and equal to 16, 4 and 3, respectively.Remarkable oscillations of modified combinants C j and multiplicity distribution P (N ) given by PYTHIA model are compared with experimental data in Fig. 3.The model and ALICE data show substantial discrepancies at small multiplicities, the experimental results for P (N ) cannot be described exactly.In particular, for the void probability P (0), we observe large difference between ALICE data and PYTHIA prediction.Since this find reflection in behavior of modified combinants.
Comparing with ALICE data, in PYTHIA model the pe- The most commonly used form of P (N ), the NBD form given by Eq. ( 1) does not describe experimental data nor the multiplicity distributions given by models.To describe P (N ) using NBD, the negative binomial distribution parameter p must depend on N [11].The probability of particle emission, p, which is constant in the standard NBD, is dependent on the multiplicity N in the way presented in Fig. 4.Both in experimental data and model the non-monotonic form of p is clearly visible (upper panel) Multiplicity distribution P (N ) of charged particles produced in proton-proton non-single diffractive interactions at √ s = 7 TeV as measured by ALICE experiment [30] (NSD class).(lower panel) The corresponding modified combinants Cj emerging from them.PYTHIA 8 [25] (dotted lines) with SoftQCD:inelastic processes and all kinematic cuts as described in the ALICE paper.
and can be described by with parameters listed in Table I.Nevertheless the behaviors are roughly similar, we observe significant difference in position N = D exp (−E 2 /2) of minimum of p.The probability of particle emission p = min for multiplicity N 23 in PYTHIA model and N 14 for ALICE experimental data.
Multiplicity distributions measured by ALICE can be successfully described by a two-component compound binomial distribution where h(N ) is the compound binomial distribution (BD + NBD) given by the generating function Comparison with PYTHIA simulation is shown in Fig. 5 for parameters given in Table II.

IV. DISCUSSION
It is worth pointing out that modified combinants evaluated from models exhibit oscillatory behavior, though the oscillation period differs from experimental data.Modified combinants C j can be expressed by the generating function: of the count probability P (N ) as: The generating function can be shown to be a sum over the "averaged" connected correlation function g (n) of all orders, (n): where m is the number of particles in a cell of the phasespace volume [31].The modified combinants C j can be expressed as an infinite series of the g (n) :   [25] with SoftQCD:inelastic processes and kinematic cuts as described in the ALICE paper [30] (NSD class).(lower panel) The corresponding modified combinants calculated from AL-ICE data [30] (NSD class).Solid lines in both panels show fits obtained with parameters listed in Table II.
Note that the correlation functions g (n) are associated with widely used cumulant factorial moments (see Appendix for more details).
Higher-order correlations, characterized by an n−body correlation function g (n) , are of general interest and have been investigated in many fields of physics, including astronomy, particle physics and quantum optics.In 1963,  Glauber predicted that the maximal value of the samepoint normalized n−body correlation function g (n) calculated for thermal light is directly related to the order of the function by a simple relationship n! [32].This n! dependence is a consequence of Wick's theorem [33], which enables higher-order correlations to be expressed using products of one-body correlation functions.The applicability of Wick's theorem is not limited to correlation functions for light, it has also been applied in many other fields; for example, it is commonly used in radio-astronomy, nuclear physics, and generally in quantum field theory.The validity of Wick's theorem has been firstly demonstrated with thermal photons, and recently proved to higher-order correlations for massive particles [34].
For correlation function with real positive parameter k we have: as for NBD, where m is the average multiplicity and the two-body correlation function determines the value of NBD shape parameter, 1/k = g (2) .
To assure oscillating behavior of modified combinants we choose correlation function in the form: which leads to to the following formula for modified combinants: with ı denoting imaginary unit.
It is remarkable that k parameter only affects oscillation amplitude, see Fig. 6.The period of oscillations (∼ 2πm) is determined by the value of m parameter, see Fig. 7.
The average number of particles in a cell given by the value of m parameter is not derived from any first principle, although some suggestions have been made to equate it with the average number of partons in the QCD cascade.Considering hadrons as a dense system of partons we expect that m ∼ Q 2 S , where Q S denotes the saturation scale [35,36].
In the PYTHIA model, the saturation scale Q S can be connected with p T 0 parameter [25].The period of oscillations depends on p T 0 , but ∼ p 1/3 T 0 dependence is very weak.By setting MultiPartonInteraction:pT0Ref=1.56 instead of its default value 2.28 we found the period of oscillations decreasing by 15% (see Fig. 8), but a more thorough re-tune is necessary in order to simultaneously obtain the correct multiplicity distribution.
Our results can, hopefully, lead to wider theoretical investigations and provide a better understanding of multi-particle production processes in hadronic collisions.

Appendix: Relationship of correlation functions with cumulants
Usually information contained in P (N ) is obtained by examining their corresponding cumulant factorial mo-ments [37] K q = F q − q−1 i=1 q − 1 i − 1 K q−i F i , (A.1)Note that the moments K q require knowledge of all P (N ) while, according to Eq. ( 4), calculation of C j (and corresponding correlation function g (n) ) requires only a finite number of probabilities P (N < j) which may be advantageous in applications.
FIG. 3.(upper panel) Multiplicity distribution P (N ) of charged particles produced in proton-proton non-single diffractive interactions at √ s = 7 TeV as measured by ALICE experiment[30] (NSD class).(lower panel) The corresponding modified combinants Cj emerging from them.PYTHIA 8[25] (dotted lines) with SoftQCD:inelastic processes and all kinematic cuts as described in the ALICE paper.

FIG. 4 .
FIG. 4. Values of NBD p parameter chosen to obtain NBD fits.See text for details.

FIG. 5 .
FIG. 5. (upper panel) Modified combinants Cj calculated from multiplicity distribution of charged particles generated in proton-proton interactions at√ s = 7 TeV using PYTHIA 8[25] with SoftQCD:inelastic processes and kinematic cuts as described in the ALICE paper[30] (NSD class).(lower panel) The corresponding modified combinants calculated from AL-ICE data[30] (NSD class).Solid lines in both panels show fits obtained with parameters listed in TableII.

FIG. 6 .
FIG. 6. Modified combinants Cj calculated for m = 4 and different values of k parameter.

FIG. 7 .
FIG. 7. Modified combinants Cj calculated for different values of m parameter.The values of k parameter were adjusted in order to obtain the same amplitude of oscillations.

3 )
FIG.8.Modified combinants Cj calculated from multiplicity distribution of charged particles generated in proton-proton interactions at √ s = 7 TeV using PYTHIA 8[25] with Soft-QCD:inelastic processes and kinematic cuts as described in the ALICE paper[30] (NSD class).Simulation performed for two values of pT0Ref parameter.Using red solid line we show the results for pT0Ref =2.28, which is default PYTHIA value.Black dashed line show the result obtained for pT0Ref =1.56.

TABLE II .
Parameters wi, pi, Ki, ki and mi of the 2-component P (N ) given by Eq. (7), used to fit the ALICE data and corresponding PYTHIA simulation, √ s = 7 TeV and |η| < 3.