Intriguing feature of multiplicity distributions

Multiplicity distributions, P(N), provide valuable information on the mechanism of the production process. We argue that the observed P(N) contain more information (located in the small N region) than expected and used so far. We demonstrate that it can be retrieved by analysing specific combinations of the experimentally measured values of P(N) which we call {it modified combinants, Cj, and which show distinct oscillatory behavior, not observed in the usual phenomenological forms of the P(N) used to fit data. We discuss the possible sources of these oscillations and their impact on our understanding of the multiparticle production mechanism.


Introduction
The multiplicity distribution, P(N), is an important characteristic of the multiparticle production process, one of the first observables measured in any multiparticle production experiment [1]. The way in which the consecutive P(N) are connected reflects the dynamics of the multiparticle production process. In the simplest case one assumes that the multiplicity N is directly influenced only by its neighbouring multiplicities (N ± 1) in the way dictated by the simple recurrence relation: (N + 1)P(N + 1) = g(N)P(N) where g(N) = α + βN, The most popular forms of P(N) emerging from this recurrence relation are (p denotes probability of particle emission): Usually the first choice of P(N) in fitting data is a single NBD. However, with growing energy and number of produced secondaries it increasingly deviates from the data for large N (see Fig. 1) and is therefore replaced either by combinations of two [2], three [3], or multi-component NBDs [4], or  Figure 1. (a) Charged hadron multiplicity distributions for the pseudorapidity range |η| < 2 at √ s = 7 TeV, as given by the CMS experiment [10] (points), compared with the NBD for parameters N = 25.5 and k = 1.45 (dashed line) and with the 2-component NBD (solid line) with parameters from [11]. (b) Multiplicity dependence of the ratio R = data/ f it for the NBD (circles) and for the 2-component NBD for the same data as in panel (a) (squares). (c) The multiplicity dependence of the modified probability of particle emission p in the MNBD as given by Eq. (5) resulting in flat R = R(N) = 1. (d) Coefficients C j emerging from the CMS data used in panel (a) compared with the NB and 2-NBD fits shown there and with the C j obtained from the MNB with modifications proposed in [9] and shown in panel (c). by some other forms of P(N) [1,[5][6][7][8]. However, as seen in Fig. 1 (a) − (b), such a procedure only improves the agreement at large N, the ratio R = data/ f it deviates dramatically from unity at small N for all fits. This observation, when taken seriously, suggests that there must be some additional information hidden in the small N region, not investigated yet [9]. In [9] we retrieved it using a single NBD form of P(N) in which we allowed for the multiplicity dependence of the particle emission ratio p = m/(m + k) in Eq. (4). It turns out that for with parameters c = 20.252, a 1 = 0.044, a 2 = 1.04·10 −9 and b = 11, one gets the desired flat behavior of R as a function of multiplicity N, now R = 1 for all N [9]. Such a choice corresponds to a rather complicated, nonlinear and non monotonic spout-like form of g(N) in the recurrence relation Eq. (1) and to a non monotonic, depending on multiplicity, probability of particle emission, p = p(N), with a sharp minimum around N = 10, after which p(N) grows steadily, see Fig. 1 (c).

Modified combinants C j
The above example shows that there is room for change in P(N) resulting in agreement with data over the whole region of N. However, the recurrence relation (1) is too restricted to be helpful in this case and in [9] we proposed a more general form of the recurrence relation, that used in counting statistics when dealing with cascade stochastic processes [12]. Contrary to Eq. (1), it now connects all multiplicities by means of some coefficients C j which define the corresponding P(N) in the following way: The coefficients C j contain the memory of particle N + 1 about all the N − j previously produced particles. They can be directly calculated from the experimentally measured P(N) by reversing Eq. (6) and putting it in the form of the following recurrence formula [9]: The coefficients C j can therefore replace the ratio R = data/ f it in quality assessment of P(N) used to fit data. The result was striking, as can be seen in Fig. 1 (d), where the coefficients C j obtained from the data used in Fig. 1 (a) show very distinct oscillatory behavior (with a period roughly equal to 16), gradually disappearing with N. It turns out that it can be reproduced only by the MNB model [9] mentioned before, which makes R(N) = 1 for all N [9]. As shown in [9,13,18] such oscillations of C j are seen for different pseudorapidity windows, in data from all LHC experiments and energies. The only condition is that the statistics of the experiment must be high enough, for with small statistics the oscillations become too fuzzy to be recognized [13,18] (the simplest way to obtain oscillatory behaviour of C j for an otherwise smooth distribution P(N) is to distort P(N) slightly at some point, this distortion then propagates further [9,18]). The reason why the single NBD is not able to reproduce data is that in this case all C j > 0 [9]: The oscillations can occur only for combinations of NBD [9,13]. However, the parameters of the 2-NBD fit used in Fig. 1 (a) result in a very small trace of oscillations of C j and we were not able to find parameters of a 2-NBD allowing for a reasonable description of P(N) and C j at the same time.
The best result obtained so far is the 3-NBD fit proposed in [3] and based on the claim that there is a place in data for a third component aiming to describe the low N events (see [3] for details, it agrees with our observations mentioned before). Fig. 2 shows the results obtained for the parameters from [3]. Note that the agreement of P(N) with data and the behavior of both the ratio R = data/ f it and the coefficients C j improved substantially. However, as one can see in Fig. 2 (b), the low N region of P(N) still shows some deviations, which, albeit rather small, result in R departing from unity (downwards) at small N and in C j missing the data for large j. Contrary to the case of the NBD, the modified combinants for the BD, cf. Eq. (2), oscillate rapidly, with a period equal to 2. In Fig. 3 (a) one can see that the amplitude of these oscillations depends on the emission probability p, in this case the C j increase with rank j for p > 0.5 and decrease for p < 0.5 (this is, however, not generally true as we shall see later, cf., for example Fig. 3 (b)). However, their general shape lacks the distinctively fading down feature of the C j observed experimentally. This means that the BD used alone cannot explain the data (see also [14]) 1 . . Results for P(N), (a); its enlarged part for N < 50, (b); ratio R = data/ f it, (c), and for C j , (d), obtained using a 3-component P(N) composed of 3 NBD as proposed in [3] (with parameters the same as in [3]).  It turns out that the coefficients C j defined by the recurrence relation (7) are closely related to the so called combinants C ⋆ j which are defined in terms of the generating function, G(z) = ∞ N=0 P(N)z N , as while preparing the final data. However, such a statement has so far not been substantiated by any known experimental analysis of the procedure used, furthermore, the peculiarities seen in the ratio R (tightly connected with the oscillations of C j ) were not addressed as well. Therefore we assumed that this is a real new effect, connected with some dynamical features of the production mechanism (in fact in [16,17] the cascade stochastic processes leading to Eq. (7) were successfully applied to multiparticle phenomenology). and which were introduced in [19] (see also [1,[21][22][23][24][25][26]). Namely, Therefore, henceforth we shall call the C j modified combinants. Note that the recurrence relation, Eq. (6) can be written in terms of C ⋆ j as: and that the C j can be expressed by the generating function G(z) of P(N) as This relation will be used in what follows when calculating C j from multiplicity distributions P(N) defined by some generating function G(z).
Because a single distribution of the NBD or BD type cannot describe data we shall check the idea of compound distributions (CD) applicable when the production process consists of a number M of some objects (clusters/fireballs/etc.) produced according to some distribution f (M) (defined by a generating function F(z)), which subsequently decay independently into a number of secondaries, n i=1,...,M , following some other (always the same for all M) distribution, g(n) (defined by a generating function G(z)) [27] 2 . The resultant multiplicity distribution, h(N) = f ⊗ g, where N = M i=0 n i , is a compound distribution of f and g with generating function 3 .
Let us take, as an example, f as a Binomial Distribution with generating function F(z) = (pz + 1 − p) K (for which the C j oscillate with a period of 2, cf. Fig. 3(a)), and g as a Poisson distribution with generating function G(z) = exp[λ(z − 1)] (for which C 0 = 2 and C j>0 = 0, cf. Eq. (13)). The generating function of the resulting Compound Binomial Distribution (CBD) is then The analytical forms of the corresponding C j and P(N) are presented in [18] (as Eqs. (136)-(138)). Fig. 3 (b) shows C j for the CBD with K = 3 and λ = 10 calculated for three different values of p in the BD: p = 0.54, 0.62, 0.66. Note that, in general, the period of the oscillations is now equal to 2λ (i.e., in Fig. 3 (b) where λ = 10 it is equal to 20). This example shows that the choice of a BD as the basis of the CD used is crucial to obtain the oscillatory behavior of the C j (for example, a compound distribution formed from a NBD and some other NBD provides smooth C j ). 2 In fact the NBD is a compound Poisson distribution with the number of clusters given by a Poissonian distribution and the particles inside the clusters distributed according to a logarithmic distribution [29]. In [28] we proposed a specific compound distribution to explain the Bose-Einstein correlation phenomenon. It consisted of a combination of k elementary emitting cells (EEC) producing particles according to a geometrical distribution. For k = const the resultant P(N) was of the NBD type, for k distributed according to a BD it was a modified NBD. However, applying it to the present situation we could not find a set of parameters providing both the observed P(N) and oscillating C j . 3 Note that for the class of distributions of M that satisfy the recurrence relation Eq. (1) the compound distribution h = f ⊗g is also given by the so-called Panjer's recurrence relation [30], , with initial value h(0) = f (0). It could be considered as a generalization of Eq. (6) used to define modified combinants with coefficients C (P) j depending additionally on N. However, Eq. (6) is not limited to the class of distributions satisfying Eq.(1) but is valid for any distribution P(N), therefore this recursion relation is not suitable for us. Unfortunately, such a single component CBD (with P(N) = h(N; p, K, λ) depending on three parameters: p, K and λ), does not describe the experimental P(N) and C j . We return therefore to the idea of using a multicomponent version of the CBD presenting two examples. The first is 3-component CBD defined as: The results of using Eq. and As one can see in Figs. 4 (c) and (d), using Eq. (19) (with parameters: K 1 = K 2 = 3, p 1 = 0.7, p 2 = 0.67, k 1 = 4, k 2 = 2.3, m 1 = 6, m 2 = 19.0 and w 1 = w 2 = 0.5) improves substantially the behavior of C j . This means that to describe data one has to use some multicomponent compound distributions based on the BD (as responsible for the oscillations in C j ) and some other distribution providing damping of the oscillations for large N (here the NBD). p+p, √ s = 7 TeV CMS, |η| < 2 ALICE, |η| < 2 ALICE, |η| < 3 ALICE, −3.4 < η < 5 Figure 5. Comparison of N C j and C j obtained from the P(N) measured by CMS [10] and ALICE [31] for three different rapidity windows. See text for details.
With such knowledge we proceeded to a description of the recent ALICE data [31] for multiplicity distributions (NSD events at 7 TeV), at three different rapidity windows: |η| < 2, |η| < 3 and −3.4 < η < 5. In Fig. 5 we show on the left panel the results for N C j obtained from the measured P(N); for comparison the previously used CMS data [10] for |η| < 2 are also shown (and they agree with the data from ALICE). The most intriguing features observed is the rather dramatic increase of both the period of the oscillations and their amplitude with the width of the rapidity window used to collect the data and, most noticeably, the previously observed fading down of their amplitude is now replaced by an (almost) constant behavior (for |η| < 3) and by a rather dramatic increase (for −3.4 < η < 5). Because, roughly, N ∼ ∆η, one would expect that at least part of the increase of amplitude could come from the increase of N with ∆η. The right panel of Fig. 5 shows that this can only partially be true, the previously observed effects remain, albeit they are perhaps not so dramatic as before.
We close our presentation by showing in Fig. 6 that the results presented in Fig. 5 can be nicely fitted, including (not shown there) P(N) from [31], using the two component CBD (BD+NBD) discussed above (Eqs. (18) and (19) with the parameters listed in Table 1. When performing these fits we were only concerned to find such values of the parameters as would reproduce the results in the best possible way (although no χ 2 estimations were used). This means that, at the moment, the values presented in Table 1 must be taken cum grano salis. In general one observes a slow increase of p and p ′ (but with some nonmonotonicity seen for p 2 ), already observed in Fig. 3. A more detailed analysis would need much more involved investigations than presented here at the moment and is currently un- Table 1. Parameters w i , p i , K i , k i and m i of the 2-component P(N), Eqs. (18) and (19), used to fit the data in Fig.  6. For completeness p ′ i = m i /(m 1 + k i ) from Eq. (18) are also included.
√ s = 7 TeV ALICE, −3.4 < η < 5 Figure 6. Multiplicity distributions P(N) (left panels) measured by ALICE [31] and modified combinants C j emerging from them (right panels) for three different rapidity windows fitted using the two compound distribution (BD+NBD) given by Eqs. (18) and (19) with parameters listed in Table 1 (see text for details). der investigation. At the moment we must admit that the problem of the physical meaning of these fits (in other words: what information on the mechanism of production of particles they convey) remains still an open one.

Some explanatory remarks
Note that in all the above discussions we always remained on the same phenomenological level, either modifying the emission probability p in the NBD, or combining it with some other distribution. So far we have not looked for physical justifications of the methods used but concentrated on reproducing P(N) and the ratio R = data/ f it or the modified combinants C j . The only more theoretically oriented attempt is Ref. [32] describing the pion multiplicity using the combinants C ⋆ j and discussing two scenarios. The first one had N sources emitting bosons without any restrictions on their number, and resulted in a NBD and smooth and diminishing combinants. The second one had M sources, each emitting only a limited number of bosons (one or two in [32]), and this resulted in a BD and oscillating combinants. We expect therefore that our C j would follow the same behavior. Note also that the NBD belongs to the class of the so-called infinite divisible distributions whereas the BD does not. In [26] it is claimed that combinants of all ranks are all non-negative if and only if the probability distribution is infinitely divisible, therefore we would expect that our C j share this property. However, for a uniform distribution in the interval (0, K) which is not infinitely divisible, the resulting C j are strictly positive, in fact C j = 2/(K + 1), which invalidates the above statement (or, at least, weakens it considerably). The true origin of the oscillations still remains not fully specified.
Let us finally note that multiplicity distributions P(N) are usually studied by analyzing factorial moments cumulant factorial moments, or their ratios [26,33], which are very sensitive to the details of the multiplicity distribution. The advantage of their use is that they seem to be well described by perturbative QCD considerations, especially their oscillations in sign as a function of the rank q [33]. On the other hand, K q can be expressed as an infinite series of the modified combinants, C j and, conversely, C j can be expressed in terms of K q [26], Note that the combinants can be, by analogy to factorial cumulants, understood as exclusive correlation integrals [1,20]. However, they differ in the region of phase space they are most suitable to study: whereas cumulants are particularly well suited for the study of densely populated phase-space bins, combinants are better suited for the study of sparsely populated regions and their calculation requires only a finite number of P(N), with N < j, which compensates the drawback caused by the requirement that one must have P(0) > 0. Additionally, the advantage of combinants is that, being finite combinations of the probability ratios P(N)/P(0), they do not suffer from a bias (empty-bin effect) present at high resolution in factorial moments and cumulants [1]. Combinants share with cumulants the property of additivity, i.e., 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 combinants are given by the sum of the independent components [26]. Because the C j are directly connected with the C ⋆ j (cf. Eq. (11)) they also share all their properties mentioned above. However, whether they share their oscillating pattern is still to be checked.
To summarize, we argue that only compound distributions based on the BD (like) and the NBD (like) components can fit adequately observed oscillations of modified combinants C j . The question of which particular theoretical mechanism is at work remains, however, still open and one may expect a number of particular models to emerge here.