A look at multiparticle production via modified combinants

As shown recently, one can obtain additional information from the measured charged particle multiplicity distributions, P(N), by investigating the so-called modified combinants, Cj\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_j$$\end{document}, extracted from them. This information is encoded in the observed specific oscillatory behaviour of Cj\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_j$$\end{document}, which phenomenologically can be described only by some combinations of compound distributions based on the Binomial Distribution. So far this idea has been checked in pp and e+e-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^+e^-$$\end{document} processes (where observed oscillations are spectacularly strong). In this paper we continue observation of multiparticle production from the modified combinants perspective by investigating dependencies of the observed oscillatory patterns on type of colliding particles, their energies and the phase space where they are observed. We also offer some tentative explanation based on different types of compound distributions and stochastic branching processes.


Introduction
Multiplicity distributions (MDs) of high energy collisions have been extensively studied in the field of multiparticle production. It is one of the first observables to be determined in new high-energy experiments. This is partly due to the ease with which such information can be obtained, and also because MDs contain useful information on the underlying production processes. Due to the inability of perturbaa e-mail: ang.h.w@u.nus.edu b e-mail: phycahp@nus.edu.sg c e-mail: mghaffar@mun.ca d e-mail: maciej.rybczynski@ujk.edu.pl (corresponding author) e e-mail: grzegorz.wilk@ncbj.gov.pl f e-mail: zbigniew.wlodarczyk@ujk.edu.pl tive Quantum Chromodynamics (pQCD) to provide a complete theoretical account for the observed MDs incorporating both the hard and soft processes, various phenomenological approaches had to be adopted. These can range from dynamical approaches in the form of coloured string interactions [1] and dual-parton model [2], to geometrical approaches [3,4] resulting in the fireball model [5], stochastic approaches [6][7][8] modelling high energy collision as branching [6][7][8] or clans [9].
The myriad of stochastic models since proposed have described the experimental data well with reasonable χ 2 /do f values with the Negative Binomial Distribution (NBD) and its variants being the most ubiquitous [10]. However, as has been proposed recently [11][12][13][14], a good fit to the MD from a statistical distribution is only one aspect of a full description of the multi-faceted set of information derivable from the MDs. A more stringent requirement before any phenomenological model is considered viable is to also reproduce the oscillatory behaviour seen in the so called modified combinants, C j , which can be derived from experimental data. In fact, this phenomenon is observed not only in pp collisions discussed in [11][12][13][14] but also, as demonstrated recently in [15], in e + e − annihilation processes. Such oscillations may be therefore indicative of additional information on the multiparticle production process, so far undisclosed. Specifically, the periodicity of the oscillations of modified combinants derived from experimental data is suggestive.
It is in this spirit that this study sets forth to understand the effects of the collision systems and various experimental observables on the period and extent of oscillations in C j . In Sect. 2, the concept of modified combinant will be reviewed in light of its connection to the earlier concept of combinant [16][17][18]. From this link, an attempt is made on the potential interpretation of modified combinant applied in the context of multi-particle production. Section 3 discusses the problem of dependence on collision system whereas Sect. 4 discusses the effect various experimental variables have on the modified combinant oscillations and summarises the key points observed.
Our concluding remarks are contained in Sect. 5 together with a tentative proposal of employing the characteristics of oscillations in experimental modified combinants to distinguish between different collision types. Some explanatory material is presented in appendices: Appendix A presents the relationship between C j and the K q and F q moments that are more familiar to the particle physics community whereas Appendix B shows the possible origin of the observed oscillations of C j based on the stochastic approach to the particle production processes.

Modified combinant and combinant
Statistical distributions describing charged particle multiplicity are normally expressed in terms of their generating function, G(z) = ∞ N =0 P(N )z N , or in terms of their probability function P(N ). One other way to characterise a statistical distribution is a recurrent form involving only adjacent values of P(N ) for the production of N and (N + 1) particles, Cast in this form, every P(N ) value is assumed to be determined only by the next lower P(N −1) value. In other words, the link to other P(N − j)'s for j > 1 is indirect. In addition, the eventual algebraic form of the P(N ) is determined by the function g(N ). In its simplest form, one can assume g(N ) to be linear in N , such that Some prominent distributions have been defined in this form. For example, when β = 0 one gets the Poisson Distribution (PD). The Binomial Distribution (BD) arises for β < 0 and β > 0 results in the Negative Binomial Distribution (NBD).
While conceptualising a phenomenological model, the form of g(N ) can be modified accordingly to describe the experimental data, cf., for example, [19,20]. However, the direct dependence of P(N + 1) on only P(N ), as seen in Eq. 1, seems unnecessarily restrictive. This constraint can be further relaxed, by writing the probability function connecting all smaller values of P(N − j) as follows [16], The coefficients C j are known as the modified combinants and forms the core of this study. They are related to the combinants C * j first defined for the study of boson production models [17,18] by the following relation [11]: Combinants were first introduced to quantify the extent any distributions deviate from a Poisson distribution. For the Poisson distribution C 0 = 1 and C j>1 = 0. In this way, any non-zero C j at higher orders indicate a deviation from the Poisson distribution. From Eq. (3), two obvious interpretations for C j follow. First, there is a one-to-one map between C * j to C j via Eq. (4). Modified combinants can be interpreted as a proxy to the extent of deviation from a Poisson distribution at different higher orders. Secondly, C j 's are the normalized weights in the series for the value of (N + 1)P(N + 1). This can be interpreted as the "memory" which P(N + 1) has of the P(N − j) term. In other words, the modified combinants are the weights in which all earlier P(N − j) values has on the current probability. In this interpretation the links between P(N + 1) to all P(N − j) values are clearly established.
One further notes that since C j 's are expressed in terms of the probability function in Eq. (3), it may be reasonable to attempt casting the modified combinant in terms of the generating function G(z) = N P(N )z N . Such an expression is immensely useful should a theoretical distribution avail itself to describe experimental data. In this case, C j can be expressed as follows: Modified combinants for some prominent distributions are shown in Table 1. Note that the generating functions of NBD and BD are in fact some quasi-power functions of z and as such can be written in the form of the corresponding Tsallis distribution [21], Table 1 Distributions P(n) used in this work: Poisson (PD), Negative Binomial (NBD) and Binomial (BD), their generating functions G(z) and modified combinants C j emerging from them whereas for q → 1 in both cases we obtain G(z) for PD. Equations (7) and (8) allow to write C j for all three distributions differentiated by the above choice of the paramater q in one formula, Note that while for the PD and NBD coefficients C j are monotonic and positive functions of rank j, they strongly oscillate for the BD. This feature will be very important in all our further analysis here.
To understand the effects of various experimental variables on oscillations of modified combinants, a mathematical expression is required for calculating the value of C j given P(N ). From Eq. (3), it follows that Note that Eq. (10) will require P(0) > 0 which is often the case as most experimental data on non-single diffraction collision exhibits enhanced void probability [12,22]. In the event that the void probability is not made available, it will be inferred from the normalization of probability.

Dependence of C j oscillations on collision systems
We shall start with a reminder of two distinct observed patterns of modified combinants, one observed in e + e − annihilation [13,15] (cf. Fig. 1) and another observed in pp scattering [12,13] (cf. Fig. 2). In the first case we use the additivity property of modified combinants, i.e. that for a random variable composed of independent random variables, with its generating function given by the product of their generating functions, G(x) = j G j (x), the corresponding modified combinants are given by the sum of the independent components. For the e + e − data we shall use then the generating function G(z) of the multiplicity distribution P(N ) in which N consists of both the particles from the BD (N B D ) and from the NBD (N N B D ): Bottom panel: the modified combinants C j deduced from these data on P(N ). They can be fitted by C j obtained from the same generating function with the same parameters as used for fitting P(N ) In this case generating function is and multiplicity distribution can be written as and the respective modified combinants are Figure 1 shows the results of fits to both the experimentally measured [23] multiplicity distributions and the corresponding modified combinants C j calculated from these data (cf. [15] for details).
In the case of pp collision the satisfactory agreement in fitting observed oscillatory pattern is obtained by using the sum of two Compound Binomial Distributions of BD with NBD type,

Fig. 2
Top panel: Multiplicity distributions P(N ) measured in pp collisions by ALICE [24]. Bottom panel: The corresponding modified combinants C j . Data are fitted using sum of two compound distributions (BD+NBD) given by Eqs. (16) and (15) with parameters: with the generating function of each component equal to As seen in Fig. 2, one gains satisfactory control over both the periods and amplitudes of the oscillations, as well as their behavior as a function of the rank j. More importantly one can reproduce the enhancement of void probability of P(0) > P(1) in addition to fitting both the P(N ) and C j . The results presented in Figs. 1 and 2 suggest the possibility that the enhanced oscillatory behavior is, perhaps, a trait of the annihilation type of the process considered.
To check this we turned to pp processes measured by UA5 [25]. Figure 3 demonstrates that the outcome is rather intriguing and brings in new questions. At 900 GeV one observes oscillatory pattern which follows that observed in annihilation process e + e − , and which can be fitted by the same kind of P(N ). However, the observed oscillatory pattern changes dramatically at 200 GeV and resembles that seen before in the pp collisions. It can still be fitted using generating function G(z) given by Eq. (12) but with Bino-mial Distribution replaced by compound distribution CBD of the Binomial Distribution with Poisson distribution, i.e., by where generating function for Compound Binomial Distribution (CBD) is given by Such replacement allows to preserve oscillating power of BD but, at the same time, to gain better control over the period of oscillations which is detemined by the mean multiplicity λ in the PD [12]. Note that the BD used at 900 GeV can be considered as such compound distribution but with the PD replaced by δ N ,1 . It means therefore that, in order to fit the annihilation data at lower energies, one has to somehow smear out this delta-like behavior. In fact, one could as well use instead of the PD a NBD with large k and p such that λ = kp/(1 − p).
We close this Section by noting that the use of G(z) in the form of Eq. (12) corresponds to a QCD-based approach based on stochastic branching processes used in [15], the so-called Generalized Multiplicity distribution (GMD), with initial number of gluons given by a BD. The links between adopting a stochastic branching approach in the study of QCD phenomena has its roots in [26]. In fact, similar approach was also formulated on general grounds in [27] where it was shown that the stochastic birth process with immigration and with initial conditions given by BD results in the so-called Modified Negative Binomial Distribution (MNBD) (both approaches are presented in more detail in Appendix B). With more general choice of initial conditions, i.e., by replacing BD by some compound distribution CD based on BD, one can, as presented here, describe also pp processes. However, in the case of pp collisions this CD is more complicated (we have now K = 3 in our BD, which could, perhaps, correspond to 3 valence quarks; additionally, to describe P(N ) we need in this case at least two such components).

Dependence of C j oscillations on phase space being tested
In addition to dependence on the collision system discussed above there are data [24,25,[28][29][30][31] (see also [14]) which allows to investigate the possible oscillatory behavior of C j in different pseudorapidity windows |η|, for different transverse momentum cuts p T and for different collision energies √ s. We shall study them in this section.  [25]. Right panels: The corresponding modified combinants C j . Data at 900 GeV are fitted by the distribution obtained from the generating function given by the product

Dependence on pseudorapidity window
The dependence of the extent of oscillations on the pseudorapidity window from which the experimental data was obtained is the most obvious. Intuitively, one would expect experimental data collected from a larger pseudorapidity phase space to be more representative of the collective behaviour of the underlying collision (e.g. e + e − , pp or pp) and the associated secondary particles. Figure 4 shows example of the observed differences between the oscillations in C j derived from different rapidity windows by ALICE Collaboration [24]. The first observation is that oscillations, which are almost non-existent at small pseudorapidity window (|η| < 1.5) are becoming very strong at the maximal pseudorapidity window (|η| < 3.4. There is also a change in the period of oscillations (where present) with a change in pseudorapidity window. In general, the period decreases from around 18 for |η| < 2.4, to approximately 11 for |η| < 1.5. The amplitude of oscillations for any smaller pseudorapidity window is too weak to discern the period. Nevertheless, the oscillations for the data from the ALICE Collaboration are relatively smooth within the pseudorapidity phase space.
Data from the ALICE Collaboration had been taken over a larger pseudorapidity window, up to η < 3.4. This allows the investigation of behaviour of C j oscillations beyond the limited window |η| ≤ 2.4 available in by CMS data (this is due to challenges surrounding the drastic drop in reconstruction efficiencies at |η| > 2.4 [28]). The bottom panel of Fig. 4 has been plotted using ALICE data from η < 2.4 to η < 3.4 for better clarity. The trend of increasing period with larger pseudorapidity window continues beyond |η| < 2.4. However, the rate of amplitude decay slows significantly between |η| < 2.4 and |η| < 3.0 and reverses at |η| < 3.4. From this observation, it is inferred that the amplitude stops its decay and reversed somewhere between 3.0 < |η| < 3.4.

Dependence on p T
In earlier study presented in [14] it was noted that the C j obtained from data obtained for p T > 100 MeV/c cut by ATLAS [29] exhibit minimal oscillation for |η| < 2.5, which are completely absent for data with p T > 500 MeV/c cut. This observation suggests that the p T phase space plays a role in the extent of C j oscillations as well. For this subsection, we will consider pp collision data obtained from the ATLAS and CMS collaborations across different p T cuts at √ s = 7 TeV. The CMS collaboration performs an extrapolation down to p T = 0 MeV/c for their MD data. This allows further exploration of the behaviour of the C j oscillations over the complete p T phase space. The resulting uncertainty due to the extrapolation is less than 1%, inclusive of systematic uncertainty [28].
Top panel of Fig. 5 presents results for data at √ s = 7 TeV, from CMS at |η| < 2.4 and from ATLAS at |η| < 2.5. The small difference in the pseudorapidity window over which they are obtained is considered insignificant, as can be seen in the close tracking of the data points from CMS and ATLAS for p T > 500 MeV/c. Note that the C j oscillations are the strongest at p T > 0 MeV/c (from CMS) while having

Fig. 4 Top panel:
The plots of C j oscillations using pp experimental data at √ s = 7 TeV derived from ALICE Collaboration over a pseudorapidity range up to |η| < 2.4 [24]. The magnitude and period is comparable to C j derived from the CMS Collaboration at the same energy and pseudorapidity. Bottom panel: C j plots from ALICE Collaboration [31] obtained for pseudorapidity up to |η| < 3.4 plotted separately for clarity. Note the increase in oscillatory magnitude at |η| < 3.4 minimal oscillations at p T > 500 MeV/c (both CMS and ATLAS). Due to the lack of availability of data points with consecutive integral N from ATLAS at p T > 500 MeV/c, the plot has to be truncated at j = 30. Unfortunately, no p T data is available from earlier experiments. The dearth of such data prohibits further investigation into the effects on oscillations between various p T cuts in pp collisions. Nevertheless, even these limited results can be very helpful in understanding the message of C j . They are very similar to what is presented in the bottom panel of Fig. 5 which shows schematic view of modified combinants C j for separate components from the two component compound distribution given by Eqs. (15) and (16) with parameters fitting experimental P(N ) shown in Fig. 2. This comparison seems to suggest that particles The plot of C j vs j with CMS data at 7 TeV with |η| < 2.4 [28] and ATLAS data [30] with η < 2.5. CMS has extrapolated its data all the way to p T > 0 MeV/c in the cited reference. This allows us to compare it with the data obtained experimentally with p T > 500 MeV/c, also from CMS. The C j derived from ATLAS data tracks that of CMS closely. Bottom panel: Schematic view of modified combinants C j for separate components from the two component compound distribution given by Eqs. (15) and (16) with parameters fitting experimental P(N ) shown in Fig. 2 with large transverse momenta mainly come from the first component (with smaller mean multiplicity) in our two component compound distribution. In other words, top panel of Fig. 5 seems to show that reducing the p T phase space eliminates (at least to some extent) one of the components.

Dependence on √ s
The reason why data from √ s = 7 TeV has been extensively exploited in the earlier parts of this work is due to the fact that oscillatory behaviour are more apparent at higher collision energies. Hints of this potential dependence on collision energy can first be observed in Fig. 3 between pp collisions at The modified combinants derived form CMS [28] and ALICE [31] across centre-of-mass energies are plotted in Fig. 6 on the top and bottom panel respectively. The difference between the data sources is that ALICE provides data up to √ s = 8 TeV while that from CMS is up to √ s = 7 TeV. To facilitate comparison, only data at |η| < 2.4 is used, on considerations that it shows the most distinct oscillatory behaviour without the amplitude blowing up. Note that CMS does not provide data obtained from wider pseudorapidity windows, which makes comparison difficult.
For the C j from CMS, there is no clear effect on the amplitude with increasing collision energies. The C j made with data from lower collision energy of √ s = 0.9 TeV appeared to have a slightly higher initial amplitude but decayed at similar rates to that from √ s = 7 TeV. There is also an increase in the period of oscillation at higher energies. On the other hand, the graph derived from ALICE data seemed to show a more distinct difference in the amplitude between data from √ s = 0.9 TeV and that from √ s = 7 TeV with a slower rate of decay. The shorter period at lower energy is also observed here, and is consistent up to √ s = 8 TeV. However, there are some observed differences in the details between both plots in Fig. 6. Some examples include the higher amplitudes of the oscillating C j for ALICE than than CMS at the same √ s, the location of the amplitudes with respect to q etc. These discrepancies can be traced back to different P(N ) values obtained between the two experiments due to slightly different methods in which measured data is being treated between the two experiments. A more detailed comparison can be found in [31]. However, this difference should not mask the trend in C j oscillations with increasing energies, which is the main point behind the plot.

Summary and discussion of results
The pseudorapidity window within which the data has been obtained appears to have the most significant effect on the oscillatory period for the corresponding derived value of C j . This feature can be clearly observed in the plot across various pseudorapidity windows from ALICE data in Fig. 4. There is direct correspondence between the size of the pseudorapidity window to the oscillation period. While C j up to |η| < 1.0 barely exhibits any oscillations, the Top panel of Fig. 4 shows an increase in period from 11 at |η| < 1.5 to 18 for |η| < 2.4. With reference to the bottom panel in Fig. 4 for large pseudorapidity windows, we see that increasing the size of the window results in a corresponding increase in oscillatory period, from 18 at |η| < 2.4 up to 23 at |η| < 3.4. Data from UA5 paints a different story. The C j oscillates with period 2 at √ s = 900 GeV at |η| < 3.0 and above. Coupled with the modified combinants derived from e + e − [15], this seems to suggest that C j from matter-antimatter ( pp and e + e − ) collisions oscillates more violently at comparatively lower energies than their pp counterparts. This may be a feature useful in distinguishing between the two types of collision data.
The second effect of larger pseudorapidity window is on the amplitude of the oscillations of C j . Referring to Fig. 4, the amplitude of oscillation increases from just below 1.5 for data from η < 1.5, to around 1.8 for C j derived form η < 2.4.
In Fig. 7 it is observed that both in UA5 and ALICE data the amplitudes of oscillations increase as a power-law from η < 3.0 onwards. The increase is more prominent for higher energies and for pp data from UA5. Note that when we use G(z) as given by Eq. (12) then amplitude of oscillations is given by [ p/(1− p)] j . If the modified combinants were to be interpreted as weights of the various P(N )'s, as discussed in Sect. 2, the oscillations in the weights are more pronounced and periodic in a larger pseudorapidity phase space.
Another aspect which the oscillatory behaviour can be discussed is in terms of the p T phase space. In the top panel of Fig. 5, results from both CMS and ATLAS data shows an unambiguous relation between p T phase space and oscillatory extent of C j . The extrapolation of the CMS data from p T > 0 MeV/c for √ s ≤ 7 TeV allows us access to full p T phase space for LHC Run 1 energies. By comparing the derived C j from both CMS, ALICE and ATLAS, it is clear The same for ALICE data at √ s = 900 GeV and for |η| < 3.0. In this case the data points are fitted against a line y(x) = 2.09 × 1.01 x . In both cases, the oscillation amplitude increases in a power-law fashion as function of x = j that like pseudorapidity, the larger the p T phase space, the larger the extent of oscillations. The comparison of these results with view of C j from separate components of distribution used to fit experimental P(N ) shown in Fig. 2 which seems to suggest that particles with large transverse momenta mainly come from the first component is very instructive and suggest further investigations which, however, go beyond the goals of this work.
On the other hand, varying the collision energies does not produce such drastic changes in the extent of oscillations as compared to pseudorapidity and p T cuts. In Fig. 6, we see that the effects of an increase in collision energy has minimal effects to the amplitude decay and the period of oscillatory behaviour. Both the amplitude and period of oscillations do not change significantly from √ s = 0.9 to 7 TeV for CMS, and up to √ s = 8 TeV for ALICE. Note that usually the oscillatory behavior of C j (as well as the lack of oscillations) is observed in the ideal cases, i.e., for P(N ) described by analytical formulas. Experimentally C j 's are obtained from the measured multiplicity distributions, which are recorded with some acceptance in limited phase space and which contain both the systematic and statistical uncertainties. However, as was shown in [15], experimental acceptance do not generate oscillations of C j . For example, acceptance procedure applied to NBD gives again the NBD with the same k but with the modified p, which is now equal to p = pα/(1 − p + pα), where α denotes the probability of the detection of a particle in the selected phase space. For a distribution described by NBD the acceptance does not alter the smooth decrease of C j . In addition, the influence of statistical and systematic uncertainties on C j was discussed in details in earlier works [12,14]. It was found that at sufficiently-high statistics, modified combinants C j become relatively insensitive to statistical uncertainties, although the effects of systematic uncertainties of the measurements still remain.
However, it turns out that, notwithstanding this sensitivity, the oscillatory signal observed in the modified combinants derived from ATLAS, ALICE, CMS and UA5 data remains statistically significant. Therefore, the regularity and periodicity of the observed oscillations cannot be results of random fluctuations but instead, justify detailed and careful analysis of oscillations in modified combinants in the study of multiplicity distributions.

Conclusion
The utility of a phenomenological approach to analysis of multiplicity distributions stems from the lack of a comprehensive theoretical explanation transcending the hard and soft regimes of QCD. If enlarging the pseudorapidity phase space results in more distinct oscillatory behaviour, then the C j oscillations could find their origins in soft hadronic collisions.
This paper discusses dependence of C j oscillations on collision systems and the impact of varying pseudorapidity, p T cuts and collision energies on the oscillatory behaviour of C j . It is clear that pseudorapidity has the greatest impact on the oscillatory behaviour among the experimental variables considered. The general trend inherent in the data shows increased oscillatory behaviour with an increase in the extent of phase space under considerations. Sampling within a larger extent of experimental phase space allows the collection of information from a larger domain. This in turn implies more representative data to be collected when the extent of phase space is large.
The way the C j oscillates between pp and pp collisions is clearly different, in terms of the order of magnitude as well as the period. For pp data from ALICE, the C j oscillates with a period of 20. This is close to the earlier discussion in Sect. 4.1 with C j oscillating at a period of 18 at √ s = 7 TeV, |η| < 2.4. In the case of pp, C j oscillates with a period of 2. Such a short period is reminiscent of our earlier work [15] exploring C j oscillations derived from e + e − collisions at √ s = 91 GeV. Based on these two observations, it seemed that at sufficiently wide pseudorapidity window, C j from particle-antiparticle collisions at different energies oscillates with period 2, while that from particle-particle collisions do not exhibit such regularity. Such power-law increase in amplitude may potentially be a characteristic of matter-antimatter collision, including that from e + e − .
Another distinguishing feature between pp and pp collisions is the order of magnitude over which the oscillations take place. At √ s = 900 GeV, C j from pp goes up to a magnitude of 10 20 while that for pp stays below 10. Should more data between the two types of collisions become available in the future, such figures can be tabulated to explore the dependence of the scaling coefficients on energy and pseudorapidity.
The relationship between C j and F q and K q moments as discussed in the Appendix A may offer some clues as to why C j derived from experimental MD data oscillates. The H q = K q /F q moments, with its roots in gluodynamics [32,33], were conceived of and observed to undergo oscillations in earlier studies. On the other hand, F q has shown to be a valuable a tool in the study of intermittent behaviour [34] in multiparticle production. Any attempts at a physical interpretation of C j can be considered in analogy to the relationship between H q and F q . However, before that, the exact physical interpretation of C j still remains open and is subject for further investigation.
Finally, we will refer to the imposing question: what lesson can be learned from the behavior of modified combinants C j deduced from the measured multiplicity distributions P (N ) in what concerns the the dynamics of the multiparticle pro-duction process. First, it seems that the oscillations of the C j are closely related to the need to use some specific form of multiplicity distribution (MD) in the description of these processes. It must be a compound distribution based on BD (CBD), which gives oscillations, with some other MD, which controls their period and amplitudes. In fact, as shown in [12] the successful use of simple sum of 3 NBDs presented in [14] is possible only because such sum acts effectively as a kind of BD. In our investigations we were usually using MD which were either compound distributions of BD with NBD (either the sum of two such compound distributions to get perfect agreement with data) or MD for the sum of multiplicities from BD (or CBD) and NBD. In all cases BD is crucial to describe the oscillatory behaviour of modified combinants. This result, if taken seriously, imposes certain restrictions on the selection of the appropriate multi-particle production model. In Appendix B we present a summary of two potential candidates for such model, both based on some specific stochastic approach, one of which was used in this work. A broader discussion on this topic, in particular what other classes of models can meet the criteria required here, would require a separate work. moments K q . Both can be defined by the multiplicity distributions, P(N ), and modified cumulants, C j . The basic quantity to start with is the generating function from which P(N ) emerges as and combinants C * j are defined as note that according to Eq. (4) N C j = ( j + 1)C * j+1 . Similarly, the respective derivatives taken at z = 1 define factorial moments, F q , and cumulant factorial moments, K q , Continuing this presentation, note that similarly as P(N ) defines F q , Note that K q , share the additive property of C j . As an example, for a random variable made up of a sum of other random variables each described by a generating function G j (z), the generating function of the sum is given by G(z) = j G j (z).
In this case, the value of K q of the sum is the sum of the K q values of the individual components, similar to how the modified combinants behave. While culmulants are suited to study the densely populated region of phase space, modified combinants are better suited for the sparsely populated regions. This can be seen from Eq. (10), which only requires a finite sum of P(N − j) terms in the calculation of C j .
Recurrence relation given by Eq. (10) follows naturally from definition of C j . Using Leibniz's formula for the j th derivative of the quotient of two functions, x = G (z) G(z) , which is just Eq. (10) used before. On a separate note, a variant of the unnormalized factorial moment F q has proved useful in the study of intermittent behaviours in high energy collisions [34]. It has been shown that if intermittent behaviours do indeed persist in the detected multiplicity spectra, the multiparticle production mechanism takes the form of a cascading process [37] via relations in the scaled factorial moments.

Appendix B: The possible origin of observed oscillations of C j
In [27] as a model for the particle production was considered the so called birth process with immigration. The production process proceeds via emission of particles from an incident colliding particle (by a kind of bremsstrahlung process) which can further produce another particles (via the birth process). This specific branching process is defined by the following evolution equation: where P(n; t) is the distribution of the number of particles at t (the parameter describing the evolution of a particle system from the initial state, t = 0 to the final state corresponding to the maximum value t = T , with T being some energy dependent parameter chosen to reproduce the energy dependence of the observed mean multiplicity), λ 0 is the immigration rate in an infinitesimal interval (t, t +dt) and λ 2 is the production rate of the birth process in the interval (t, t + dt).
In [15] we have used specific, QCD based, realization of such approach based on the stochastic branching equa-tion (describing the total multiplicity distribution of partons inside a jet, [6]), In both approaches, defined by Eqs. (B.1) and (B.2), one has to define initial condition. For a set number of initial particles, P(n; t = 0) = δ n,k , one gets G(z; t = 0) = z k (this is the case of the GMD discussed in [15]). For initial condition for P(n; t = 0) chosen in a form of binomial distribution, with two new parameters, α representing the production rate of additional particles present at t = 0, and K denoting their maximal number, one gets boundary condition G(z; t = 0) = ∞ n=0 P(n; t = 0)z n = [1 + α(z − 1)] K . (B.5) (used in [27]), which leads to the following generating function: where κ = exp (λ 2 T ) − 1 and ξ = λ 0 /λ 2 . Note that this is simply just a product of generating functions for the BD and NBD, G(z) = G B D z; p = 1 − κ(1 − α); K · G N B D z; p = κ 1 + κ ; k = K + ξ , (B.7) and the respective modified combinants are given by Eq. (14).
In the case when initially some complex objects (firebals, cluster, jets and so on) are produced and subsequently each of them produces secondary particles, the corresponding multiplicity distribution is described by an appropriate compound distribution. If initial objects (sources) are produced according to a BD and subsequeny production process is defined by Eqs. (B.1) or (B.2) with initial condition P(n; t = 0) = δ n,0 (i.e., multiplicity distribution from sources is given by the NBD), we have final compound distribution defined by the generating function H (z) = G B D G N B D (z; p = κ, k = ξ); p = α, K (B.8) as given by Eq. (16). In this case combinants C j oscillate with period equal to ∼ 2m, where m = κξ/(1 − κ) denote mean multiplicity from single source.