Compelling Evidence of Oscillatory Behaviour of Hadronic Multiplicities in the Shifted Gompertz Distribution

Study of charged particle multiplicity distribution in high-energy interactions of particles helps in revealing the dynamics of particle production and the underlying statistical patterns, by which these distributions follow. Several distributions derived from statistics have been employed to understand its behaviour. In one of our earlier papers, we introduced the shifted Gompertz distribution to investigate this variable and showed that the multiplicity distributions in a variety of processes at different energies can be very well described by this distribution. The fact that the shifted Gompertz distribution, which has been extensively used in diffusion theory, social networking and forecasting, has been used for the first time in high-energy physics collisions remains interesting. In this paper, we investigate the phenomenon of oscillatory behaviour of the counting statistics observed in the high-energy experimental data, resulting from different types of recurrence relations defining the probability distributions. We search for such oscillations in the multiplicity distributions well described by the shifted Gompertz distribution and look for retrieval of additional valuable information from these distributions.


Introduction
The simplest observable in high-energy interactions is a count of charged particles produced in a collision and its mean value. Its distribution measured in full or partial phase space forms both a tool for studying models and a probe for particle dynamics. A large number of statistical probability distribution functions (PDF) have been used to understand its behaviour. These include Koba, Nielsen, and Olesen (KNO) scaling [1], Poisson distribution [2,3], binomial and negative binomial [4][5][6] distributions, lognormal distribution [7], Tsallis distribution [8,9], Weibull distribution [10], modified forms of these, and several other distributions. NBD has been one of the most extensively used. It was very successful until the results from UA5 collaboration [11,12] were published. A shoulder structure was observed in the multiplicity distribution in p p collisions, showing its violations. It is also well established by various experimental results that NBD fails with increasing deviations with the growing number of charged particles produced. In order to improve the agreement with data, 2-component or 3component NBD fits [13,14] were also used.
In one of our recent papers, we introduced the shifted Gompertz distribution [15], henceforth named as SGD, to investigate the multiplicities in various leptonic and hadronic collisions over a large range of collision energies. The distribution was first introduced by Bemmaor [16] as a model of adoption of innovations. The two nonnegative fit parameters define the scale and shape of the distribution. This distribution has been widely studied in various contexts [17][18][19]. In our earlier work [15], we proposed to use the SGD for studying the charged particle multiplicities in high-energy particle collisions and showed from a detailed study for collisions in full phase space and also in limited phase space that this distribution explained the experimental data very well in high-energy particle collisions using leptons and hadrons as probes. Subsequently, we also used it to calculate the higher moments of a multiplicity distribution which also serve as a powerful tool to unfold the characteristics and correlations of particles [20]. We also used 2-component shifted Gompertz distribution, named as modified shifted Gompertz (MSGD), to successfully improve the agreement between data and fit. The details are given in our paper [15].
Wilk and Wlodarczyk, in one of their recent publications [21,22], pointed out that the 2-component or multicomponent fits improve the agreement only at large N (number of charged particles) but not at small N. They showed that the ratio data/fit deviates significantly from unity for small N. In a pursuit of retrieving additional information from measured probability of producing N particles PðNÞ, they have proposed the multiplicity distribution (MD) by a recurrence relation between the adjacent distributions PðNÞ and PðN + 1Þ. This corresponds to the assumption of a connection existing only between the production of N and N + 1 particles: The multiplicity distribution is then determined by the function form of gðNÞ, the simplest being a linear relation: where µ and ν are the parameters of the linear dependence. The more general form of recurrence relation introduced in Reference [21,22] which connects the multiplicity N + 1 with all smaller multiplicities has the form All multiplicities are then connected by means of some coefficients C j , which redefine the corresponding PðNÞ in the way such that the coefficients C j connect the probability of particle N + 1 with probabilities of all the N − j previously produced particles. These coefficients can then be directly calculated from the experimentally measured PðNÞ by exploiting the relationship. It is shown that C j shows a very distinct oscillatory behaviour which gradually diminishes with increasing N and nearly vanishes. The details are given in Section 3.
In the present work, we use shifted Gompertz distribution and its modified forms using the data at high energies from p p and pp interactions to understand the existence of such oscillatory behaviour and to check if we obtain the results consistent with the ones from [21,22].
In Section 2, we provide the essential formulae for the probability distribution function of the shifted Gompertz distribution and modified 2-component shifted Gompertz distributions, in brief. A very brief description of how the oscillations have been estimated in the multiplicity distributions by Wilk and Wlodarczyk [21,22] is included for the sake of completeness.
Section 3 presents the analysis of experimental data, the fitted shifted Gompertz distributions, the fitted modified shifted Gompertz, and the distributions giving out the oscillatory behaviour. Discussion and conclusion are presented in Section 4.

Shifted Gompertz Distribution (SGD)
Let X be any nonnegative random variable having the shifted Gompertz distribution with parameters b and ζ, where b > 0 is a scale parameter and ζ > 0 is a shape parameter. The value of the scale parameter determines the statistical dispersion of the probability distribution. The larger the value of the scale   Advances in High Energy Physics parameter, the more the distribution spread out, and the smaller the value, the more concentrated the distribution. The shifted Gompertz density function can take on different shapes depending on the values of the shape parameter ζ. It is a kind of numerical parameter which affects the shape of a distribution rather than simply shifting it or stretching or shrinking it. The multiplicity distribution is measured as the probability distribution of a number of particles being produced in a collision at a particular energy of collision and follows certain phenomenological and statistical models. The probability distribution function of X is given by The mean value ðE½XÞ of shifted Gompertz distribution is given by where γ ≈ 0:5772156 stands for the Euler constant (also referred to as the Euler-Mascheroni constant). It is well established that at high energies, the most widely adopted, negative binomial distribution [4][5][6], fails and deviates significantly for the high multiplicity tail, from the experimental data. To extend the applicability of NBD, another approach was introduced by Carruthers and Shih [4][5][6]. In this case, a weighted superposition of two independent NBDs, one corresponding to the soft events (events without minijets) and another to the semihard events (events with minijets), is obtained. These distributions combine merely two classes of events and not two different particle-production mechanisms. We used the same method to obtain the superposed distribution and call it 2-component shifted Gompertz distribution (2-component SGD) as given by equation (6). The multiplicity distribution of each component is independent SGD. The concept of superposition originates from purely phenomenological considerations. The two fragments of the distribution suggest the presence of the substructure. Each component distribution has two fit parameters, namely, scale and shape parameters. The best fit overall distribution to the experimental data, with optimised parameters, also gives an estimate of fraction, α, of the soft collisions, at a given c.m.s. energy. The dynamics of particle production is understood in terms of weighted superposition of soft and semihard contributions. Though these superimposed physical substructures are different, the weighted superposition mechanism is the same. The physical substructures are described by the same SGD multiplicity distributions and corresponding correlation functions, which are QCD inspired genuine self-similar fractal processes [4][5][6]. Same as NBD, SGD allows to describe the multiplicity distribution on purely phenomenological grounds. This may help in differentiating between different phenomenological models. Details are included in our earlier publication [15]: where α is the fraction of soft events and ðb 1 , ζ 1 Þ and ðb 2 , ζ 2 Þ are, respectively, the scale and shape parameters of the two distributions.

Modified Forms of Shifted Gompertz Distribution.
In this paper, we adopt a different approach and investigate what 3 Advances in High Energy Physics kind of changes in the structure of the multiplicity distribution described by the SGD is necessary in order to describe the same data by a single SGD, with accordingly modified parameters b and ζ. To describe data using only a single SGD, we allow the parameter b to depend on the multiplicity N, as suggested by Wilk and Włodarczyk [21,22]. To obtain an exact fit of the distribution to the experimental data, a nonmonotonic dependence of b on N is introduced. This way, the scale parameter b remains the same in nature but varies in accordance with the number of particles produced. Such a change means that we preserve the overall form of the SGD: where a 1 , d, and c are parameters. This leads to the modification of SGD (equation (4)) which describes the data very well. We call this as the modified-SGD1 (MSGD1). When another nonlinear term with a coefficient a 2 is added [21,22] to bring improved agreement with the data: we call this second modification as MSGD2. Further, we investigate the possibility of retrieving some additional information from the measured PðNÞ.

Analysis and Results
Equation (3) can be reversed, and a recurrence formula can be obtained for the coefficients C j for an experimentally measured multiplicity distribution PðNÞ, as below:  The errors on the coefficients C j are calculated from the variance: Since the coefficients C j are correlated, the last term of equation (10) introduces dependence of the error in C j on the errors of all coefficients with i < j. This leads to a cumulative effect with a large increase of errors with increasing rank j. However, despite such large errors, the values of hNiC j lie practically on the curve and the points do not oscillate in the limits of errors. Hence, the errors can be estimated reasonably well, by neglecting this cumulative effect.
In the present work, calculations are performed using the data from different experiments and following two collision types: (i) pp collisions at LHC energies ffiffi s p = 900, 2360, and 7000 GeV [23] are analysed in five rapidity windows, jηj < 0:5 up to jηj < 2:4 (ii) p p collisions at energies ffiffi s p = 200, 540, and 900 GeV [11,12] are analysed in full phase space as well as in rapidity windows, jηj < 0:5 up to jηj < 5:0 The charged hadron multiplicity experimental distributions are fitted with the SGD (equation (4)), the 2component SGD (equation (6)), MSGD1 (equation (7)), and MSGD2 (equation (8) Table 1 gives the χ 2 /ndf for all fits at different energies and rapidities. In case of 2-component SGD, the α values are taken from Reference [15].  Table 2 gives χ 2 /ndf for all fits at different energies and rapidities. It may be observed in the cases of pp collisions; MSGD2 fits the data well in comparison with other distributions, particularly at higher energies. However, for the case of p p collisions, in most of the cases, 2-component SGD improves the fits and explains the  Figure 5: Ratio R = P data ðNÞ/P SGD ðNÞ for the p p data shown in Figure 2 (red circles) and of the corrected ratio R = P data ðNÞ/P MSGD2 ðNÞ (black squares). 6 Advances in High Energy Physics data well. For a comparison between pp and p p collisions at the same ffiffi s p = 900 GeV, the trend is nearly the same and MSGD2 fit the data best.
In Figures 3 and 4, we show the ratio plots for multiplicity dependence of the ratio R = P data ðNÞ/P fit ðNÞ for the pp data shown in Figure 1 obtained from the fits SGD and MSGD2. Figure 5 shows similar ratio plots for multiplicity dependence of the ratio for the p p data shown in Figure 2. As can be seen from Figures 3 and 5, there are systematic deviations from the fits of SGD from the data at low and high multiplicities. The deviations get enhanced with increasing energy and high multiplicity values, as can also be observed in Figures 1 and  2. In addition, a structure at smaller multiplicities can also be observed. In order to understand this structure, the modified forms of SGD have been introduced as MSGD1 and MSGD2 in equations (7) and (8). The ratio R calculated with MSGD2 becomes closer to unity in all the cases, though the deviations are still present. For the possibility of retrieving some additional information from experimental multiplicity distribution, the recurrence relation given in equation (9) is used to calculate the coefficients C j . In some cases, the 2-component SGD fits exceptionally well leading to the R value around unity, as shown in Figure 4.
The coefficients C j are calculated for the pp and p p data at different centre-of-mass energies and for various pseudorapidity windows. Figure 6 shows C j for pp data in jηj < 2:4 and for the p p data for jηj < 3:0. The figures show a very distinct oscillatory behaviour in both the cases. For the case of pp interactions, oscillations occur with amplitude decreasing with the rank j at all energies. It shows that the effect of an increase in centre-of-mass collision energy has minimal effect on the amplitude and the period of the resulting oscillations. However, for the p p interactions, the trend is reversed, with the amplitude of oscillations growing with the rank j and decrease in collision energy. This intriguing property has also been observed by Rybczyński et al. [24]. The way C j oscillates between pp and p p collisions is clearly different and may be a characteristic of matter-antimatter collision. Abramovsky and Radchenko in their paper [25] have described the particle production in inelastic collisions in terms of quark and gluon strings. They have described the multiplicity distributions in terms of 2-NBD and 3-NBD and have shown how the pp and p p collisions are fundamentally different, which may lead to the observed differences. In another interesting study by Ang et al. [26], similar differences have been observed in p p (UA5) and pp (ALICE) data. Figure 7 shows the coefficients C j calculated for the pp collision data at 7000 GeV c.m. energy and for p p collisions at 200 GeV, for different pseudorapidity windows. They all show the distinct oscillatory behaviour with amplitude increasing with the pseudorapidity window for both pp and p p collisions. It is also observed that the oscillations die out with increasing rank j for all jηj bins for pp collisions, whereas for p p collisions, the oscillations grow stronger with rank j with increasing jηj bin size and somewhat random only in jηj < 5:0 bin. Similar observations are also observed by Rybczyński et al. [24] in the CMS and ALICE data [23,27]. In Figures 6 and 7, the errors on the data points are not shown for the reason that the error bars intermingle and blur the figures.
The coefficients C j are evaluated by fitting the SGD, 2component SGD, MSGD1, and MSGD2 distributions to the data. The variation of these coefficients with j is shown in Figure 8 for pp data in one pseudorapidity window for different energies. We find that C j evaluated from the SGD fit do not show this oscillatory behaviour, as compared with the data. However, with the 2-component fit, they start showing the oscillatory pattern, which further gets enhanced with MSGD1 and MSGD2 fits, following the data closely. In case of MSGD2 fits, coefficients C j follow almost exactly the oscillatory behaviour of C j obtained directly from the data at ffiffi s p = 7000 GeV. For ffiffi s p = 2360 GeV, it is MSGD1, and for ffiffi s p = 900 GeV, the 2-component SGD follow the experimental values better. Similar results are obtained by analysing the data for different pseudorapidity windows both of  Figure 6: Coefficients C j obtained from (i) the CMS data of pp collisions at different energies for one pseudorapidity window jηj < 2:4 (a) and (ii) the UA5 data of p p collisions at different energies for pseudorapidity window jηj < 3 (b). 7 Advances in High Energy Physics pp and p p collisions. However, we show the results for 7000, 2360, and 900 GeV pp collisions for only jηj < 2:4 and similarly for 900, 540, and 200 GeV p p collisions for jηj < 1:5, in Figure 9. It may be observed that none of the fits consistently follow the p p data trends. This is also seen for other η windows. To avoid too many similar figures, we present only the representative figures.
The coefficients C j evaluated from equation (9) depend on Pð0Þ. In the experimental data from complex detectors, such as CMS at the LHC, the probability Pð0Þ is very large as compared to Pð1Þ. Due to large experimental uncertainties associated with this bin, Pð0Þ is often omitted for the conventional fits to the data. However, Pð0Þ is the only bin which is very sensitive to the acceptance as explained in References [21,22]. To show the sensitivity to the value of Pð0Þ, we show in Figure 10 the coefficients calculated by using the values P ð0Þ ± δ for the pp data at ffiffi s p = 900 GeV for jηj < 2:4, where δ is the error on Pð0Þ measurement. The coefficients vary with different periods of oscillations, around the values calculated from Pð0Þ, as shown in the figure. Figure 11 shows the oscillatory behaviour when Pð0Þ is not considered; C j are calculated starting with Pð1Þ. Coefficients C j still show   Advances in High Energy Physics the oscillatory behaviour but with much reduced oscillation amplitude, with oscillations dying out quickly. In equation (9), the coefficients C j connect each probability, with every other probability. For example, PðN + 1Þ connects to PðN − jÞ, the probabilities of particles produced earlier. The most important feature of this recurrence relation is that C j can be directly calculated from the experimentally measured PðNÞ. In an interesting case study, starting with the SGD, we make changes in successive probabilities by 2%: we put Pð10Þ = P SGD ð10Þ + Δ and Pð11Þ = P SGD ð11Þ − Δ with Δ = 0:02P SGD and study the variation of C j as a function of j. The results are shown in Figure 12 for pp collisions at different energies but within the same jηj bin. Similarly, Figure 13 shows the plots for p p collisions at ffiffi s p = 900, 540, and 200 GeV for jηj < 3. The apparently insignificantly small changes in probability resulted in rather dramatic spikes occurring on the original P SGD and with rapidly falling amplitudes. This points to the sensitivity of the coefficients C j . Such a change is then provided by the MSGD, whereby spike influences then the consecutive coefficients C j and brings them to agreement with those obtained from the experimentally measured PðNÞ. With increasing value of j, smaller are the values of C j and hence weakly influencing

Conclusion
In this paper, we show and reaffirm that the MDs possess a fine structure which can be detected experimentally and analysed in terms of a suitable recurrence relation, such as the one in equation (9). The coefficients C j in the recurrence relation, which are directly connected with the combinants, give a compelling evidence that phenomenon of oscillatory behaviour of the modified combinants exists in the experimental data on multiplicities. The coefficients C j have been calculated from the shifted Gompertz distribution and its modified forms: weighted superposition of 2-component  10 Advances in High Energy Physics shifted Gompertz parametrizations and modified shifted Gompertz distributions including nonlinearity to two different orders, equations (7) and (8). The shifted Gompertz distribution, which we introduced in our publication [15], does not show any oscillatory behaviour. However, its modified forms show the oscillatory behaviour and agree with the data very well. The oscillations are large at low multiplicities for the pp data and tend to die out at large multiplicities. In case of p p collisions, the oscillations follow a reverse pattern. The behaviour of oscillations observed in present studies is very similar to what is observed in the case of negative binomial distribution (NBD), by the authors who pioneered the concept.

Data Availability
All the data used in the paper can be obtained from the references quoted or from the authors.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.