Three-component multiplicity distribution, oscillation of combinants and properties of clans in $pp$ collisions at the LHC

Multiplicity distributions of charged particles measured by the ATLAS Collaboration in proton-proton collisions at $\sqrt{s}=8$ and 13 TeV are analyzed in the framework of a weighted superposition of three negative binomial distribution functions. The examination of the experimental data confirms the existence of a narrow peak at low multiplicities found earlier at $\sqrt{s}=0.9$ and 7 TeV. The peak is described by the third separate component of the total distribution. The energy dependence of the multiplicity characteristics of the three-component scenario is studied. It is demonstrated that the third multiplicity component might be responsible for oscillations of combinants observed lately in the analysis of the multiplicity distributions measured in proton-proton collisions at the LHC. Some consequences of the three-component description of the data for the clan parameters are analyzed.

The multiplicity of produced particles is an important characteristic of the final state in high energy protonproton collisions. The shapes of the charged-particle multiplicity distributions (MDs) measured with high accuracy can provide valuable information on various processes in these interactions. The distributions reflect correlations in the system encoded in an integrated form. At lower collision energies, the two-particle correlations are dominant and a single negative binomial distribution (NBD) [1] gives a satisfactory description of the MD up to the ISR energies.
At higher collision energies, a shoulder structure of the MD at high multiplicities [2] and oscillation of the ratio H q = K q /F q of the cumulant to factorial moments [3] have been found. No pattern of this kind appears in a phenomenological fit by a single NBD [4]. It has been proposed [5] to describe the observed shoulder structure as a weighted superposition of soft events (without minijets) and semi-hard events (with mini-jets), each of the NBD type. The suggestion was emphasised by the interesting and remarkable observation that the oscillations of H q vs. rank q obtained from data is reproduced by the K q over F q ratio calculated from the superposition of two NBDs. The appearance of the second multiplicity component under the tail of the MD is connected with cumulants each of which involves an infinite "cumulative" sum over all multiplicity probabilities P (n).
The extrapolation of the weighted superposition of the two defined classes of events from the TeV to the multi-TeV energy domain and a study of the properties of particle clans raised intriguing questions concerning the onset of a third class of events in the tail of the distribution, both in the full phase space [6] and in the limited windows in pseudorapidity [7]. Later study [8] of the MDs * Electronic address: zborovsky@ujf.cas.cz measured at √ s = 7 TeV indicated existence of the third multiplicity component at low n. In this region, near the maximum of P (n), there is a suitable tool, the method of combinants [9], to be used in searching for new phenomena in the behavior of the MDs. The combinants C(i) have similar additive properties as cumulants which are expressible as an infinite sum over all probabilities. The latter are therefore convenient for an analysis of the total shape of the MDs. On the other hand, the combinants can be expressed as a finite combination of the ratios P (n)/P (0). The combinants are thus extremely suitable for study of the final set of the multiplicity frequencies which make up the experimental data at low multiplicities.
In this paper we study in terms of a three-component description the data on MD of charged particles produced in pp collisions at the LHC at new high energies, √ s = 8 and 13 TeV. Each component is represented by a single NBD and parametrized by two parameters,n and k. The total distribution is weighted sum of three NBDs. The paper extends an earlier analysis [8] performed with data on MD measured at √ s = 0.9 and 7 TeV. The study of the ATLAS measurements [10][11][12][13] show that the two-component description of the MD is unsatisfactory. Within the multi-component NBD parametrization, the structure of the data indicates the necessity of a third component in the region of low multiplicities. We demonstrate that, unlike the two-component fits to the MDs measured at the LHC, the phenomenological fits with the third NBD component at low n can give oscillations of the combinants multiplied by their rank discovered lately [14] in multiplicity data measured by the CMS [15] and ALICE [16] Collaborations. Although the sensitivity of the oscillations to the systematic uncertainties of the measurements is large, the natural applicability and power of the method of combinants in the analysis of the fine structure of the MDs observed in pp collisions at the LHC is of interest in the search for and study of new phenomena emerging in the global characteristics of particle production.
The second part of the paper is dedicated to a study of the obtained results within the clan structure analysis of data on MD. The qualifying assumption is that each of the components of the weighted superposition used to parametrize the experimental data has the NBD form. The emergence of the peak in the MD at low n, described by the third NBD, has a serious impact on the clan parameters of the first and the second component of the total distribution. The parameters changed dramatically in comparison with the two-NBD analyses considered usually in different scenarios of soft and semi-hard events. The most appealing observation is the increase of the average number of clans in the high-n (semi-hard) component with the center-of-mass energy √ s. This is in sharp contrast with two-NBD parametrizations of the multiplicity data. The contribution of the third NBD to the low multiplicity part of the charged particle distribution has a serious consequence for the average number of particles per clan in the first and the second component. It turns out that the dominant component with largest probability is characterized with few clans containing many particles. We illustrate that the second, high-n component consists of large number of clans containing less particles. The properties of clans of the three-NBD parametrization of the MDs measured at the LHC and some corresponding mechanisms of particle production are discussed.

B. Weighted superposition of NBDs
The data on MD provide information on various processes underlying the production mechanism contain information as regards correlations in the system and can serve as a test for probing the dynamics of the interaction. The particle production is based on quantum chromodynamics (QCD). The multi-particle processes responsible for the final multiplicity involve also soft scales, including hadronization, which remain beyond the reach of perturbative QCD. The complex phenomena that influence MD are therefore hard to describe in detail and one has to rely on phenomenology. One of the most successful distribution functions used in describing the probability distributions of produced particles in pp/pp collisions is the two-parameter NBD, P (n,n, k) = Γ(n + k) Γ(k)Γ(n + 1) wheren is the average multiplicity and k characterizes the width of the distribution. The history and genesis of NBD can be found e.g. in [17]. An important application of this distribution in the particle production is connected with cascading mechanisms. A frequently used and popular interpretation of NBD comes from the clan model [21] where a particle emits additional particles in a self-similar branching pattern. The clans are produced independently and contain particles of the same ancestry. The Poisson distribution is the distribution of clans consisting of single particles. It is obtained for k = ∞. A superposition of NBDs was exploited by different approaches to multiple production in hadron collisions. It is based on the decomposition where P (n,n i , k i ) is given by (1). In this paper we study high statistic multiplicity data [11][12][13] obtained by the ATLAS Collaboration at the LHC in the framework of a weighted superposition of three NBDs. The aim is to achieve a detailed description of the data including the shape of the distributions in the region of maximal values of P n . Specifically, we consider the function (2) for N = 3 and extract values of the eight independent parameters (n i , k i , α 2 and α 3 ) corresponding to the ATLAS data with different transverse momentum cuts at the energies √ s = 8 and 13 TeV. The results are compared with our previous analysis [8] at √ s = 0.9 and 7 TeV and with a two-NBD parametrization of the same data.

C. Analysis of new ATLAS data at midrapidity
The multiplicity distributions of charged particles produced in pp collisions have been measured at the LHC by different experiments, in different kinematic regions, at different energies and for different classes of events. The experimental data accumulated by the ATLAS Collaboration are exceptional, for they have small systematic errors and are based on the analysis of an extremely large number of events. The ATLAS data at √ s = 8 TeV [11] and √ s = 13 TeV [12,13] exploit a similar methodology to that used at lower centre-of-mass energies [10]. The analyzed charged particle MDs were obtained in the same phase-space regions with application of the same multiplicity cuts in the same pseudorapidity window |η| < 2.5. The recorded number of events for the conditions p T > 100 MeV/c and n ch > 1 exceeds 9 million for both energies, √ s = 8 and 13 TeV. A similar large number of events was considered in the data sample with the higher transverse momentum cut p T > 500 MeV/c and n ch > 0. The measurements in both regions give a severe restriction on models of multiparticle production in pp collisions at high energies. This concerns in particular the weighted superposition of two NBDs, which is usually attributed to the classification of events into soft and semi-hard events with respect to the momentum transfer in parton-parton scatterings.

MDs at low transverse momentum
The most inclusive phase-space region covered by the ATLAS measurements corresponds to the conditions p T > 100 MeV/c and n ch > 1. We have fitted the data [11,13] on MD measured in the interval |η| < 2.5 FIG. 1: MD of charged particles measured by the ATLAS Collaboration [11,13] in the pseudorapidity interval |η| < 2.5 for pT > 100 MeV/c, n > 1 at a √ s = 8 TeV and b √ s = 13 TeV. The error bars include both the statistical and the systematic uncertainties summed in quadrature. The solid lines represent a superposition of three NBDs with the parameters from Tables I and II, respectively. The dash-dot, dash and dash-dot-dot lines show single components of the total MD in this region with a weighted superposition of three NBDs. The decomposition of the total MD at the energies √ s = 8 TeV and 13 TeV is shown in Fig. 1a, b, respectively. The experimental data are indicated by symbols and the fitted three-component function (2) is depicted by the solid line. The dash-dot, dash and dashdot-dot lines represent single NBD components of the total distribution. The corresponding parameters and values of χ 2 are quoted in Tables I and II. Terms of the fitting procedure are explained in Appendix A.
One can see from Fig. 1 that the decomposition into three components is similar at both energies and reveals  identical features to those found [8] in the same class of events at lower √ s. The dominant component with the largest probability gives the main contribution α 1n1 to the total average multiplicity. The other two components contribute to the high-and low-multiplicity region, respectively. The average multiplicities of the first and the second component,n 1 andn 2 , increase with energy. This results in broadening of the total distribution. The average multiplicityn 3 ≃ 11 of the third component is nearly energy independent. Within the errors quoted in Table I, the probabilities α i show a weak energy dependency as well. The values of the parameters k i increase with decreasing probabilities α i . The first NBD with the largest probability α 1 is characterized by the smallest parameter k 1 . The NBD component under the peak of the distribution at low multiplicities is narrow with large value of k 3 . We have fitted the ATLAS data with a weighted superposition of two NBDs. The obtained values of the parameters for the two-component hypothesis are given in Tables I and II. Figure 2 shows the relative residues of data with respect to the two-NBD fits. The points correspond to the measurements of the ATLAS Collaboration in the pseudorapidity interval |η| < 2.5 for p T > 100 MeV/c,  Tables I and II, respectively. The symbols correspond to data on MDs measured by the ATLAS Collaboration [11,13] in the interval |η| < 2.5 with the cut pT > 100 MeV/c, n > 1 at a √ s = 8 TeV and b √ s = 13 TeV. The error bars include both the statistical and the systematic uncertainties summed in quadrature. The solid lines represent a three-NBD superposition with the parameters from Tables I and II, respectively. The inserts show the detailed behaviour of the residues at low n n > 1. The inserts show the detailed structure of the residues at low multiplicities. The solid lines are given by the three-component description of the ATLAS data with the parameters quoted in Tables I and II. One can see that the two-NBD approximation of the measured MDs is unsatisfactory at both energies. The corresponding values of χ 2 /dof are too large, especially at the energy √ s = 13 TeV. The high statistic ATLAS data measured with relatively small systematic uncertainties manifest a distinct peak around n ∼ 11. The description of the peak clearly seen in the residues in Fig. 2 was obtained by the third negative binomial component withn 3 ≃ 11, similar to the demonstration in [8] at lower √ s. The third component is depicted by the dash-dot-dot lines in Fig. 1.

MDs with the cut pT >500 MeV/c
The ATLAS Collaboration presented data on MD produced in pp collisions in the separate phase-space region defined by the conditions p T > 500 MeV/c and n ch > 0 at the energies √ s = 8 TeV [11] and 13 TeV [12]. The distributions were measured in the window |η| < 2.5. They differ from earlier analyses of primary charged particles in that charged particles with a life time 30 ps < τ < 300 ps were considered as secondary particles and thus excluded.
The ATLAS data at √ s = 8 and 13 TeV together with three-NBD fits are shown in Fig. 3a, b, respectively. The symbols and the lines have the same meaning as in Fig. 1. The corresponding parameters and values of χ 2 /dof are given in Tables III and IV. One can see some similarities when comparing the description of the most inclusive and the p T -cut data shown in Figs. 1 and 3, respectively. The component with the largest probability α 1 has the smallest parameter k 1 . The parameters k i increase as the probabilities α i decrease. The average multiplicityn 3 ≃ 3 of the third component is nearly energy independent. The probability α 3 of the third component is non-zero also in the p T -cut data sample. On the other hand, there are differences between the description of data shown in Figs. 1 and 3, respectively. With the imposed p T cut, the value ofn 3 becomes significantly smaller relative to the average multiplicities of the first and the second NBD component. The multiplicityn 1 does not increase with energy but shows signs of saturation. The contribution α 2n2 of the second component to the total average multiplicity increases rapidly with √ s and becomes dominant at √ s = 13 TeV. A decreasing tendency of the probabilities α i with the index i (seen in minimal p T -biased data) is well visible for p T > 500 MeV/c at √ s = 13 TeV only. A similar description applies to the hierarchy of the parameters k −1 i characterizing the widths of the single NBD components.
The normalized residues relative to the two-NBD parametrization (see Tables III and IV) of the ATLAS data measured in the pseudorapidity window |η| < 2.5 with p T > 500 MeV/c, n > 0 at the energies √ s = 8 TeV [11] and 13 TeV [12] are depicted in Fig. 3c, d, respectively. Both data show sizeable discrepancies with respect to the weighted superposition of two NBDs. The values of χ 2 /dof = 43.8/34 at √ s = 8 TeV and χ 2 /dof = 548/76 at √ s = 13 TeV for the two-NBD fits are too large, especially at the higher energy.
The residues manifest a clear peak visible at low multiplicities around n ≃ 3. The description of the peak is obtained by the third negative binomial component with n 3 ≃ 3 shown by the dash-dot-dot lines in Fig. 3a [11,12] in the pseudorapidity interval |η| < 2.5 for pT > 500 MeV/c, n > 0 at a √ s = 8 TeV and b √ s = 13 TeV. Normalized residues of the data on MDs relative to weighted superposition of two NBDs (P 2NBD n ) with the parameters listed in Table III and Table IV    13 TeV confirms the existence of the peak [8] emerging near n ≃ 3 at √ s = 7 TeV. Such an effect is negligible in the ATLAS measurements with p T > 500 MeV/c at √ s = 0.9 TeV (see Fig. 6(b) in [8]).

D. Energy dependence of three-NBD description
The ATLAS data on MDs in two p T -cut regions show a distinct peak at multiplicities where the soft production processes dominate. The peaky structure can be well described by a separate component within the threecomponent parametrization of the experimental data. A description of the MDs by superposition of three NBDs reveals some properties which are compared below for both analyzed kinematic regions. It concerns the energy dependence of the corresponding parameters obtained at lower energies ( √ s = 0.9 and 7 TeV) [8] and in this analysis ( √ s = 8 and 13 TeV). Figure 4 shows the probabilities α i of single NBD components in dependence on √ s. The parameters α i were obtained from fits to the ATLAS data on multiplicities in the window |η| < 2.5 for (a) an almost inclusive (p T > 100 MeV/c) and (b) a transverse momentum cut (p T > 500 MeV/c) data sample. The symbols with error bars and shaded rectangles correspond to the parameter values obtained by minimization of Eqs. (A1) and (A4), respectively. Despite considerable errors, there are some trends visible in the behavior of the probabilities. The values of α i in Fig. 4a show a weak or nearly no energy dependence in the multi-TeV energy region. Errors of their determination allow to trace a hierarchy in relation to the index i which is rather stable against √ s. The same ordering of α i is obvious in Fig. 4b at √ s = 13 TeV only. The ordering of α 2 and α 3 for the data with p T > 500 MeV/c is due to large errors unclear at √ s =7 and 8 TeV. The probability of the third NBD component is non-zero in all analyzed cases. Figure 5 shows the √ s dependence of the average multiplicitiesn i of single NBDs in the window |η| < 2.5 for (a) p T > 100 MeV/c and (b) p T > 500 MeV/c cuts. The values ofn i extracted from the ATLAS data [10][11][12][13] are depicted in the log-log plot to check their power behavior. As seen from Fig. 5a, the average multiplicitiesn 1 andn 2 demonstrate a power increase with the energy √ s. Both dependences are well parametrized in , which is the maximal rapidity of pions in the pp c.m. system. A similar description seems to be valid for the √ s dependence of n 2 shown in Fig. 5b, thought for the multiplicityn 1 with p T > 500 MeV/c it is problematic to draw such a conclusion. The third NBD component at low n reveals remarkable properties. Its average multiplicityn 3 is nearly energy independent for both p T -cut data samples.
The inverse values of the parameters k i of the weighted superposition of three NBDs are depicted as a function of √ s in Fig. 6. They were found from an analysis of the same data as the parameters shown in Figs. 4 and 5. As  Fig. 6a is observed for the data sample with p T > 500 MeV/c at √ s = 13 TeV only. The ordering of the parameters at lower energies is hard to tell and their energy dependence is unclear from Fig. 6b.

E. Combinants of the MDs
A general form of the MD is useful to make a characterization in terms of its deviations from the Poisson distribution. The natural logarithm of the Poisson generating functional is given by a linear dependence on its argument. The higher order terms of the power series expansion of the logarithm of a general generating func-tional denote the deviation from the Poisson distribution. The expansion coefficients C(i), named combinants [9], characterize the MD in terms of the generating function The combinants possess the additivity property and are expressible as a finite combination of the ratios P (n)/P (0). Some of them are permitted to have negative values [9]. Individual interpretations of the coefficients depend on different models of multiparticle production. The quantity of interest, named here "cumulative combinant", is iC(i). In the model of independent boson production with geometric distributions of particles in each production mode (e.g. in the mode with an average number of bosons n i with momenta in the interval (p i , p i + dp i )), the cumulative combinant (i + 1)C(i + 1) is the mean number of modes which have an occupation number greater than i [23]. The cumulative combinants can be rewritten in the form (i+1)C(i+1) = N C i where N = G ′ (z = 1) is the average multiplicity and the coefficients C i are commonly used in the relation Such a kind of recurrence formula occurs in cascadestochastic processes [24,25]. The relation was applied to the parametrization of the MD of charged particles in hadron-hadron collisions at high energies [26,27].
The quantities N C i calculated by Eq. (4) from MDs measured by the CMS [15] and ALICE [16] Collaborations in pp collisions possess remarkable oscillating properties [14]. Unlike this, the two-NBD superposition used to fit the data is not able to account for the oscillating structure of N C i .
Next we show that oscillations of the cumulative combinants N C i = (i+1)C(i+1) can arise from the three-NBD superposition fitted to the data with the third component, which accounts for the peak at low n. We calculate the coefficients N C i from Eq. (4) using the values of P (n) obtained from the corresponding fits. Figure 7 shows the quantities N C i in dependence on the rank i calculated from a weighted superposition of (a) three and (b) two NBDs used to fit the data on MD measured by the ATLAS Collaboration in the interval |η| < 2.5 with the transverse momentum cut p T > 100 MeV/c. The parameters of the fitted functions are taken from Table 1 of [8] (for √ s =0.9 and 7 TeV) and from Tables I and II of the present paper (for √ s =8 and 13 TeV). One can see that two-NBD fits to the ATLAS data give a smooth behavior of the cumulative combinants in dependence on their rank. In contrast to this, the three-NBD fits to the ATLAS data are characterized by the oscillations of N C i at low i. The shape of the first oscillation is nearly energy independent. The position of the first minimum corresponds FIG. 7: Energy dependence of the cumulative combinants N Ci calculated from a weighted superposition of a three and b two NBDs fitted to the charged-particle MDs [10,11,13] measured by the ATLAS Collaboration in the interval |η| < 2.5 for pT > 100 MeV/c to the energy independent value of the average multiplicityn 3 ∼ 11 of the third NBD component of the total distribution. As discussed in [14], the oscillations can be connected with memory (or correlations) in the produced multiparticle system. The correlation of particles in the minima is weaker. The first minimum corresponds to the production of particles within the third independent NBD component which makes the average correlation of the total system weaker. The consecutive oscillations represent non-monotonic loss of memory (or correlations) away from a considered multiplicity n.
The cumulative combinants obtained from a weighted sum of three and two NBDs used to fit the ATLAS data in the interval |η| < 2.5 with the cut p T > 500 MeV/c are shown in Fig. 8a, b, respectively. One can see that the oscillations have disappeared for the data sample with FIG. 8: Energy dependence of the cumulative combinants N Ci calculated from a weighted superposition of a three and b two NBDs fitted to the charged-particle MDs [10][11][12] measured by the ATLAS Collaboration in the interval |η| < 2.5 for pT > 500 MeV/c FIG. 9: The cumulative combinants N Ci calculated from a weighted superposition of a three and b two NBDs fitted to the charged-particle MDs [15] measured by the CMS Collaboration at √ s =7 TeV for pT > 0 in different pseudorapidity intervals p T > 500 MeV/c even for the three-NBD description. This is despite the fact that a three-NBD fit to the data is still needed to account for the peak at low multiplicities (see Fig. 3c, d). The disappearance of the oscillations is caused by the small value ofn 3 ∼ 3 (cf. the diminishing of the oscillations for small windows in Fig. 9a wherē n 3 = 2.7 for |η| <0.5).
The data on MD of charged particles measured by the CMS Collaboration in the central interaction region [15] allow one to study the cumulative combinants in dependence on the pseudorapidity window. As shown in [14], the amplitude and periodicity of the oscillations increase with the window size. For small windows they practically vanish. Figure 9a demonstrates that similar properties of the cumulative combinants can be obtained from three-NBD fits [8] to the CMS data. As seen from Fig. 9b, the two-NBD fits [8] to the CMS data do not give any oscillations. The cumulative combinants calculated from weighted superposition of two NBDs used to fit the data are decreasing functions of the rank i. In contrast, the combinants calculated from three-NBD fits to the data reveal oscillations at low i. The magnitude of the oscillations decreases with the decreasing window size. The larger periodicity is in the larger windows. The first min-imum corresponds to production of particles within the third-NBD component, which makes the average correlation of the total system weaker. Its position (for larger windows) corresponds to the pseudorapidity dependence of the average multiplicityn 3 (see Table 4 of [8]) of the third-NBD component of the total distribution. Let us note that, for the CMS data, the multiplicity n = 0 was excluded from the fits. Therefore the cumulative combinants shown in Fig. 9 are calculated with the parameters quoted in Tables 4 and 5 of [8] without any correction to the value of P (0) obtained from the fits.
We have studied the sensitivity of the oscillations of the cumulative combinants to changes of the probabilities α i of single NBD components of the total distribution. For this purpose, we used auxiliary parameters ǫ i by which the corresponding probabilities α i are gradually turned down. In order to keep the total probability equal to unity, the attenuation of α 1 with ǫ 1 < 1 is compensated by the proportional amplification of α 2 and α 3 , which is realised as follows: Similar equations were used for weakening of the second and the third component by ǫ 2 and ǫ 3 , respectively. Figure 10a demonstrates the increase of the amplitude of oscillations when the probability α 1 of the main NBD component is reduced by 5 or 10%. The locations of the minima and maxima of the oscillations remain intact as the other parameters were unchanged. Figure 10b depicts the sensitivity of the cumulative combinants to the probability α 2 of the NBD component under the tail of the total distribution. One can see that the first wave up to n ∼ 24 is not affected, even for ǫ 2 = 0. The behavior of the cumulative combinants in dependence on ǫ 3 is shown in Fig. 10c. One can see that the amplitudes of the oscillations diminish with decreasing probability α ′ 3 of the third NBD component at low multiplicities. The oscillations disappear completely in the case of α ′ 3 = 0. The uncertainties of the quantities N C i and, more specifically, correlation of their size with the pattern of their oscillations are sensitive to statistical fluctuations and systematic uncertainties of the raw data to such an extent that the wavy structure in N C i may be rendered insignificant given sufficiently large experimental uncertainties, as demonstrated above by varying the fitted three-NBD parameters. This will be further confirmed by direct extraction of modified versions of the cumulative combinants from experimental data on MDs in Sect. G.
The natural applicability of the method of combinants to the study of fine structure of the MDs appears to be useful. This together with the relation of the combinants to the finite set of multiplicities which make up the data seems to be an appropriate tool in justification of the existence of the third multiplicity component emerging in pp collisions at the LHC at low multiplicities. The correlation structures of multi-particle states created in high energy collisions are often studied in the framework of the clan concept introduced to interpret NBD occurring in different reactions over a wide range of energies. In hadron-hadron collisions, the properties of clans have been investigated and discussed within a twocomponent model used to describe the MD of the produced particles by a superposition of two NBDs [5,18]. The first component, connected with soft events, reveals invariant properties as a function of the c.m. energy [19]. The second component, under the shoulder of the MD at high n, is associated with semi-hard processes including the production of mini-jets. It is believed that mini-jets with low p T are an important part of particle production at high energies [20].
The clan model is based on the decomposition of NBD to Poisson distributed clusters or independent bunches of produced particles [21,22]. According to the common definition, a single cluster (clan) contains at least one particle. The number of particles inside the clan follows a logarithmic distribution. Among possible interpretations of the logarithmic distribution let us mention the time average of time-dependent cascading [21] and a model of cascade processes with collapsing of clans [29]. The average number of clans,N , and the average number of particles per clan,n c , are expressed in terms of the NBD parametersn and k as follows: For a weighted superposition of NBDs, the clan characteristics of a single NBD depend on the number of components used for a description of the experimental data. The analysis of the data in terms of two and three NBDs is instructive, because it allows one to study changes of clan parameters with the emergence of a distinct peak observed in MD at the LHC at low n. Figure 11a shows the √ s dependence of the average number of clans calculated from a weighted superposition of two NBDs. The parameters were calculated from fits to the ATLAS data measured in the interval |η| < 2.5 for p T > 100 MeV/c. The symbols with error bars and shaded rectangles are given by parameter values using Eqs. (A1) and (A4), respectively. As one can see, the average number of clans of the first component,N 1 , is approximately constant with energy (except the last shaded rectangle corresponding to large value of χ 2 /dof = 776/81). A similar observation [28] follows from the measurements of the CMS Collaboration [15] in the interval |η| < 2.4. The most peculiar feature of the two-NBD description is the decrease of the average number of clansN 2 of the semi-hard component with energy (upper panel in Fig. 11a). This property, confirmed by the analysis [28] of CMS data, was discussed [18] within different scenarios for soft and semi-hard events.
The √ s dependence of the average number of clans is Energy dependence of the average number of clansNi calculated from a weighted superposition of a two and b three NBDs fitted to the charged-particle MDs [10,11,13] measured by the ATLAS Collaboration in the interval |η| < 2.5 for pT > 100 MeV/c. The symbols with error bars and shaded rectangles correspond to the parameter values obtained by minimization of Eqs. (A1) and (A4), respectively quite different in the three-NBD model. As seen from Fig. 11b, the values ofN 1 for the dominant component with largest probability α 1 are nearly constant at high energies. The average number of clans in the second component,N 2 , increases with energy and becomes much larger thanN 1 at high √ s. This is a new feature, substantially different from the parametrization of the data by two NBDs (cf. the upper panel in Fig. 11a). Another observation is that the values ofN 3 corresponding to the third NBD are relatively high and nearly constant. In a certain sense, the energy independence ofN 1 and N 3 substitutes the approximate constancy of the average number of clans of the soft component in the two-NBD scenario (cf. the lower panel in Fig. 11a).
FIG. 12: Energy dependence of the average number of particles per clannci calculated from weighted superposition of a two and b three NBDs fitted to the charged-particle MDs [10,11,13] measured by the ATLAS Collaboration in the interval |η| < 2.5 for pT > 100 MeV/c. The symbols with error bars and shaded rectangles correspond to the parameter values obtained by minimization of Eqs. (A1) and (A4), respectively Figure 12 shows the average number of particles per clan,n ci , calculated from a weighted superposition of (a) two and (b) three NBDs fitted to the ATLAS data in the interval |η| < 2.5 for p T > 100 MeV/c. In the two-NBD scenario (Fig. 12a),n c1 for the soft component does not change much with energy. The value ofn c2 grows with energy and becomes larger thann c1 . Such a trend is compensated with the decrease of the average number of the semi-hard clans with √ s, as is illustrated in the upper panel in Fig. 11a.
The situation becomes completely different for the three-NBD description. As demonstrated in Fig. 12b, the average number of particles inside clans of the first (dominant) NBD component,n c1 , increases rapidly with √ s. The values ofn c2 of the second component under the tail of MD reveal growth with energy as well, but not as rapid as in Fig. 12a. The substantial difference is that, contrary to the two-NBD description,n c1 >n c2 in the multi-TeV energy domain. The average number of particles inside clans of the third NBD component,n c3 , is only slightly larger than unity and depends weakly on √ s. Based on the performed analysis of the ATLAS data with weighted superposition of three NBDs, we draw the following conclusions concerning the clan properties of the first and the second NBD component of the total distribution. The first component, corresponding to the dominant class of events, contains large clans consisting of many particles. The number of particles inside the clans grows strongly with energy. The average number of such clans is small (N 1 ≃ 4) and independent of √ s. The small and constant number of clans suggests that there are few (or no) "mini-jet clans" in these events.
The multiplicity component under the tail of the distribution corresponds to a class of events in which much more clans are produced than in the first, dominant one. The average number of clans of the second component increases with energy and reaches large values (N 2 ∼ 15) at high √ s. The corresponding number of particles per clann c2 is somewhat smaller thann c1 in the high energy region. It grows with √ s as well. The large, increasing number of clans and the growing number of particles per clan suggest that there are many clans "mini-jets" in this class of events. The properties of the second multiplicity component under the tail of the MD are characteristic for semi-hard processes with abundant production of minijets.
We have studied the pseudorapidity dependence of the clan parameters using the data [15] measured by the CMS Collaboration in the intervals |η| < η c = 0.5, 1.0, 1.5, 2.0, and 2.4 at √ s = 7 TeV. Figure 13 shows (a) the average number of clansN i and (b) the average number of particles per clann ci calculated from two-NBD description [8] of the CMS data in dependence on the width of the pseudorapidity window. The behavior of the parameters calculated from three-NBD fits [8] to the same data is depicted in Fig. 14.
As one can see from Figs. 13a and 14a, the η c dependence of the average number of clans of the first component,N 1 , is nearly the same for two-and three-NBD parametrization. Despite large errors, a difference is visible in the behavior ofN 2 for the semi-hard events with mini-jets, especially in the large pseudorapidity windows. The values ofN 2 increased slightly for the three-NBD description and do not show a sign of saturation for larger values of η c . The average number of clans grows with η c for all three components. Another difference between the parametrization of the CMS data with two and three NBDs is seen in the η c dependence of the average number of particles per clan. There is some indication that clans of the first component contain more particles than the clans of the second one (n c1 >n c2 ) when fitting the data with weighted superposition of three NBDs in larger windows. This is in accord with the clan analysis of the ATLAS data shown in Fig. 12.
The above observations are consistent with the interpretation of 1/k i as an aggregation parameter [18] 1 where P i (N, m) is the probability that m particles belong to N clans. As seen from Figs. 11b and 12b, clans of the second component are more numerous and contain fewer particles in comparison with the clans of the first com-ponent. This means that there is much less aggregation of particles in the semi-hard events with mini-jets than in the first, dominant class of events. Such a property is reflected by the relation 1/k 2 < 1/k 1 observed in central pseudorapidity intervals at √ s = 7 TeV (Fig. 18(a) in [8]) and at different energies in the window |η| < 2.5 (Fig. 6a).
We have studied the influence of transverse momentum on clan structure analysis using the ATLAS data on MDs [10][11][12][13] with different p T cuts. Figure 15 shows the ratios of the clan parameters extracted from fits to the data with p T > 500 MeV/c and p T > 100 MeV/c at √ s = 13 TeV. The ratios for the two-and three-NBD superposition are depicted in the upper and lower panels, respectively. As one can see from Fig. 15a, the average number of clans in the third component shows the greatest suppression with the p T cut. It means that these clans, consisting of very few particles, are produced mostly with low transverse momenta. For the three-NBD description, the average number of clans of the first componentN 1 is reduced the least, slightly less thanN 2 . This trend is quite opposite in the two-NBD model, as depicted in the upper panel in Fig. 15a. The suppression of the average number of particles per clan with p T is illustrated in Fig. 15b. The value ofn c3 is reduced minimally. It is natural because, by definition, a clan contains a minimum of one particle and because of the very small number of particles in these clans (n c3 ∼ 1). At the same timen c1 is diminished considerably with the harder p T cut for the three-NBD description. The opposite follows from the two-NBD scenario. There is a very small suppression of particles inside clans of the soft multiplicity component, as depicted for the first NBD in the upper panel in Fig. 15b. Let us interpret the results obtained by the description of the ATLAS data on MD with different p T cuts when using weighted superposition of three NBDs. The relatively small suppression of the average number of clans N 1 and considerable reduction of the average number of particles per clann c1 suggest that clans of the first, dominant component contain particles with a wide spectrum of transverse momenta. There are many particles in these clans (see Fig. 12b), some of them with larger, some of them with smaller p T . The imposed cut p T > 500 MeV/c reduces the number of particlesn c1 inside the clans by almost half, as those with small p T are lost. At the same time,N 1 is least suppressed with the p T cut. This observation suggests that the large clans of the first multiplicity component consisting mostly of soft particles must contain also some particles with considerably high p T . It is therefore reasonable to assume that there are semihard processes with relatively large momentum transfer in the dominant class of events.
The properties of the clan structure of the second multiplicity component are different. As depicted in the lower panel in Fig. 15b, the number of particles per clan n c2 is less suppressed with the cut p T > 500 MeV/c than n c1 . The relatively small reduction ofn c2 means that fewer particles are affected by the p T cut. It suggests that there are fewer particles with low p T in the clans of the second component in comparison with the clans of the first, dominant one. At the same time (see lower panel in Fig. 15a), the average number of clans,N 2 , is reduced similarly and even slightly more thanN 1 . The considerable suppression ofN 2 and small reduction ofn c2 indicate that the transverse momenta of particles inside clans of the second component have a quite narrow distribution concentrated around somep T . Thisp T is higher than the average <p T > of particles of the first component, most of which carry low transverse momenta. Moreover, considering thatn c2 is obviously smaller thann c1 for p T > 100 MeV/c at √ s = 13 TeV (see Fig. 12b), one can make the following conclusions. The clans of the second NBD component contain fewer particles than the large clans of the first component. The particles of the smaller clans carry on average higher p T than most of the particles in clans of the dominant component. Their transverse momenta are concentrated around somep T . The average number of these clans,N 2 , increases with energy and becomes much larger thanN 1 at high √ s. The mentioned properties of the clan structure of the second semi-hard component resemble some basic features of mini-jets: clustering of particles of close transverse momenta with increasing frequency at high energies. The application of the clan structure analysis to MDs measured by the ATLAS Collaboration with different limitations on transverse momenta reveals properties which differentiate between two-and three-NBD description of the data. It concerns both the average number of clans and the average number of particles per clan, which are reduced with varying intensity when applying the higher p T cut. As seen from Fig. 15, there are opposite trends of the suppression for both models in the first and the second multiplicity component. In order to clarify what is the cause of such a difference, one can look at the p T dependence of the corresponding NBD parameters. Figure 16a shows the ratios of the average number of particles,n i , obtained from the fits to the ATLAS data measured in the interval |η| < 2.5 for p T > 500 MeV/c and p T > 100 MeV/c. The ratios of the NBD parameters 1/k i extracted from the data are depicted in Fig. 16b. As concerns the first and the second multiplicity component, the ratios ofn i follow a similar trend for two-and three-NBD parametrization of the data. The average multiplicitiesn 1 andn 2 are suppressed almost proportionally for both models. Contrary to this, the corresponding ratios of the aggregation parameters, 1/k i , behave completely opposite for parametrization of the ATLAS data with two and three NBDs. The aggregation of particles into clans increases strongly with the cut p T > 500 MeV/c for the soft component in the two-NBD model. At the same time, the aggregation parameter 1/k 2 changes very little for the semi-hard component. In the three-NBD model, the aggregation of particles into clans increases slightly with p T for the first, dominant component. The aggregation is slightly stronger for the semi-hard component with mini-jets. As concerns the third NBD, it is hard to tell because of the large error. The observations indicate that it is just another aggregation of particles into clans, which results in different properties of clans in two-and three-NBD models used for the description of the ATLAS data.

G. Discussion
A detailed examination of the sensitivity of cumulative combinants to the systematic uncertainties of measurement and to the unfolding procedures applied to raw data requires information on the response matrix for the considered experiments as well as knowledge of the experimental methods used by obtaining the final MD.
Though a study of the raw data is not the subject of the present paper, we discuss here an estimation of errors of the cumulative combinants based on the published data. The combinants are expressible in terms of the first i probability ratios, P n /P 0 , n = 1, .., i, which depend on P 0 . For the ATLAS data, the analysed events are required to satisfy special criteria in order to minimize the systematic uncertainties. One of them is the minimal number of charged tracks that depends on the particular phase space region. The most inclusive phase space region covered by the measurements corresponds to the conditions n ch ≥ 2 and p T > 100 MeV/c. The experimental information on P 0 and P 1 is missing in these data. For the CMS data at √ s = 7 TeV, the measured value of P 0 is very large. Due to the experimental difficulties connected with the large error of this probability and because of the rise of the MDs in the zeroth bin, P 0 is usually omitted in the NBD fits to the data.
In order to avoid inaccuracies as to the assumptions with regard to the multiplicity probabilities and their errors in the first two bins, let us consider a modification of combinants based on the relation 4. Instead of the full MD {P 0 , P 1 , P 2 , ...} one can apply this recurrence relation to the truncated sets {P 2 , P 3 , P 4 , ...} or {P 1 , P 2 , P 3 , ...}. The corresponding modified quantities (i+1)C 2 (i+1) and (i+1)C 1 (i+1) are defined by the formulas and respectively. The lower index 2 (1) at C means that the truncated probability set begins with P 2 (P 1 ). As illustrated below, the quantities (i+1)C 2 (i+1) ((i+1)C 1 (i+1)) defined by Eqs. (8) ( (9)) reveal similar oscillating properties to the cumulative combinants calculated by Eq. (4) from a superposition of NBDs with parameters obtained from the corresponding fits. Figure 17 shows the coefficients (i + 1)C 2 (i + 1) in dependence on the rank i. Their values were calculated from the data on MD (n ch ≥ 2) measured by the ATLAS Collaboration in the interval |η| < 2.5 with the transverse momentum cut p T > 100 MeV/c at √ s =13, 8, 7 and 0.9 TeV. The error bars shown in the figure were obtained from experimental errors of the data on P n . The latter include both the statistical and the systematic uncertainties summed in quadrature.
In determination of the errors of (i + 1)C 2 (i + 1) we have proceeded in the following way (see Appendix A). Taking into account the positive correlation between adjacent multiplicity bins and the anti-correlation between the opposite sides of the MD distribution maximum, we have considered the distributions shifted to the right and left with respect to the error variation in single multiplicity bins. The right shifted distribution is constrained by the upper values of the errors of P n for n > m and by their lower values for n < m where m is the multiplicity with the maximal value of P m . The left shifted distribution was considered in the analogous way, exchanging mutually upper and lower errors at both sides of the maximum. The right shifted distribution of P n is responsible for the upper error bars of (i + 1)C 2 (i + 1) at maxima and the lower error bars near the minimum visible in Fig. 17. The left shifted P n governs the error bars in the opposite direction. For higher values of the rank i, one can see a scattering of points around the curves calculated from three-and two-NBD fits to the experimental MDs. This corresponds to a statistical spread of the mean values of P n with respect to the fits. As the statistical uncertainties of the data are much smaller than the spread of the mean values in the ATLAS data, the error bars shown in Fig. 17 are practically hidden within the symbols for increasing i.
Because of relatively small systematic uncertainties of the measurements, one can see a sizeable wavy structure of (i + 1)C 2 (i + 1) in Fig. 17, which clearly discriminates between the two-and three-NBD parametrization of the ATLAS data with p T > 100 MeV/c. The difference is seen in the region of low i and corresponds to a distinct peak near the maximum of the measured MD. Figure 18 shows the coefficients (i + 1)C 1 (i + 1) calculated from data on MD (n ch ≥ 1) measured by the ATLAS Collaboration in the window |η| < 2.5 with the cut p T > 500 MeV/c at √ s =13, 8, 7 and 0.9 TeV. The errors of the coefficients were obtained from experimental uncertainties of P n in the same way as in Fig. 17. Except for a few points at low i in Fig. 18a, the errors are smaller than the size of the displayed symbols. One can see that there are no oscillations of (i + 1)C 1 (i + 1) for FIG. 18: a The coefficients (i + 1)C1(i + 1) calculated from data on MD measured by the ATLAS Collaboration [10][11][12] in the interval |η| the ATLAS data with p T > 500 MeV/c at all considered energies. The quantities calculated from the data on MD reveal approximately the same behavior as the cumulative combinants obtained from weighted sum of three and two NBDs used to fit the data (see Fig. 8). The disappearance of oscillations is connected with the small average multiplicityn 3 ∼ 3 of the third NBD component for the data with higher p T cut (cf. the diminishing of the oscillations in small windows in Fig. 9a wheren 3 = 2.7 for |η| < 0.5). The situation is somewhat different for the CMS measurements with p T > 0 [15] in sufficiently large pseudorapidity intervals. Figure 19a shows the quantities (i + 1)C 1 (i + 1) in dependence on the rank i calculated from the n ch distributions in the window |η| < 2.4 at √ s =7 TeV. The black symbols correspond to the mean values of P n . The error bars were obtained in the same way as in the case of the ATLAS data. The right shifted distribution of P n results in a large oscillation of errors that increases strongly in the region of higher rank i. The right shift determines the upper errors at maxima and the lower errors at minima of the first oscillations and creates a repeating pattern. The left shift of the distribution is mostly responsible for errors in the opposite direction. One can see from Fig. 19a that the structure of the error bars agrees with the initial wavy character of the black symbols quite well. The errors grow very rapidly at extrema with the increasing number of oscillations. Such a behavior is caused by the systematic uncertainties of the CMS data, which are much larger than in the ATLAS measurements. The systematic uncertainties govern the errors of residues of the MD with respect to the two-NBD fits [8]. They form a characteristic envelope around the mean values of P n that follows the fine structure of the residual mean values quite accurately.
We have examined how a reduction of the systematic uncertainties of MD can affect the information on the oscillatory behavior of the coefficients (i+1)C 1 (i+1). Assuming that the systematic errors of P n could be reduced by 50%, the oscillating structure of cumulative combinants including errors is expected to be similar to Fig. 19b. Such an increase in accuracy would help in discriminating between two-and three-component description of the data considerably. The oscillating pattern depicted in Fig. 19a can help to discriminate between the two-and three-NBD hypotheses also on the basis of the following considerations. Let us assume that the distribution of mean values of P n is described exactly by two NBDs, e.g. with parameters quoted in [8]. In this model situation we take the statistical and systematic uncertainties provided by the CMS Collaboration and construct the quantities (i+1)C 1 (i+1). The result for the interval |η| < 2.4 at √ s = 7 TeV and p T > 0 is shown in Fig. 20a. No characteristic pattern of oscillation can be seen there. If we assume, however, that the mean values of P n follow exactly a three-NBD curve, with parameters taken e.g. from [8], and take the experimental errors from the corresponding measurement, the uncertainties of (i+1)C 1 (i+1) constructed as above look like Fig. 20b. One can see that the oscillating pattern based on the three-NBD hypothesis shows a very similar structure, as depicted in Fig. 19a. The later corresponds to the modified cumulative combinants and their errors calculated from data published by the CMS Collaboration.
We examined the sensitivity of the oscillation pattern of errors to the change of the truncation condition from n ≥ 1 to n ≥ 2. For that purpose we used data on MD measured by the CMS Collaboration [15] for p T > 0 in the interval |η| < 2.4 at √ s = 7 TeV and calculated the coefficients (i+1)C 2 (i+1). The result depicted in Fig. 21 demonstrates that the error pattern seen in Fig. 19 is not destroyed with the increase of the truncation point. A comparison of both figures shows that omitting of P 1 from the construction of the modified combinants results in a reduction of the amplitude of their errors and is accompanied by a slight decrease of their periodicity. Let us note that the structure illustrated in Fig. 21b is similar to that depicted by full symbols in Fig. 17b.
A detailed analysis of the cumulative combinants requires study of raw data and knowledge of proper response of detectors. The discussed estimation of their uncertainties indicates that oscillations of the cumula-tive combinants seem to be a true property of n ch distributions. The ATLAS data are most conclusive in that respect. The CMS data are less decisive as to the oscillations, although some supplementary support for such a hypothesis is seen in larger pseudorapidity intervals.
In the following part we discuss some ideas concerning the development of the clan structure in the context of the present study of the clan parameters obtained from the three-NBD description of the analyzed data. The emerging picture is consistent with a stochastic-physical FIG. 21: a The coefficients (i + 1)C2(i + 1) calculated from data on MD measured by the CMS Collaboration [15] in the interval |η| < 2.4 at √ s = 7 TeV for pT > 0. The error bars correspond to the statistical and the systematic uncertainties of Pn summed in quadrature. b The coefficients (i + 1)C2(i + 1) calculated from the same data with the systematic errors reduced by 50%. The solid (dashed) lines were calculated from three-(two-) NBD fits to the measured MD using Eq. (8) scenario [29] of clans' evolution involving the ingredients of QCD branching processes. In this scenario, clans, during their (QCD) evolution, can perform collapses onto single particle states. The states evolve further by subsequent branching. The collapses of clans are most probable in the early stage of their development. We consider that such collective behavior reflects early neutralization of color configurations and simulates rapid freeze-out of QCD degrees of freedom inside clans. The processes of a clan's collapsing define directions along which the final particles originated from their ancestors are spreaded. Those clans that experienced (one or more) collapses can be identified (under additional conditions) with mini-jets. The "mini-jet clans" are expected to have different properties in comparison with clans that evolved without collapses. It concerns the transverse momentum distribution of particles which, due to collapses, should be concentrated around somep T . Moreover, it is natural to expect that the number of clans mini-jets should increase with energy and should become larger than the number of clans of which the evolution happened without collapses. Due to collapses, clans identified with mini-jets should contain less particles than clans that did not suffer any collapse.
As can be seen from Figs. 11b and 12b, this picture is in accord with the values of the clan parameters extracted from three-NBD description of the ATLAS data. The first multiplicity component consists of large clans containing many particles. The average number of these clans is low and constant with energy. These properties are in line with the conjecture that the occurrence of the mini-jets in the dominant component is small. The second NBD component describes a high-multiplicity tail of the total distribution. As usually considered, the component is enriched with mini-jets. Their symptoms concern clustering of particles into smaller clans with increasing intensity at high energies. This assumption complies with the clan properties extracted from the experimental data. First, the clans of the second component contain fewer particles in comparison with the clans of the first one (n c2 <n c1 ). The value ofn c2 increases with energy. Second, the average number of the smaller clans,N 2 , increases with energy and becomes much larger thanN 1 at high √ s.
The properties of clans in the third component are different. There is only one or very few particles in these clans, as illustrated by empty squares in Fig. 12b. It points to little or nearly no branching in such a case. The parameter 1/k 3 is small (see Fig. 6a), reflecting a minimal aggregation of particles in comparison with the first and the second component. A characteristic feature of the three-NBD fits to the ATLAS data is the approximate constancy of the average multiplicityn 3 with energy. Moreover, the average number of clans,N 3 ≃ 9, depends very weakly on √ s as well (see Fig. 11b). These observations suggest that there are different mechanisms responsible for particle production in the first two components of the MD and the third one emerging at low n.
The performed study of the clan parameters indicates that the peak around maximum of the MD is influenced by a mechanism which depends weakly on energy and is characterized by minimum branching and small transverse momenta. This is typical for soft particle production. It is unlikely that the peak is a remnant of diffractive events. In the central rapidity intervals studied at the LHC, the single-diffractive events influence 0-and 1-bin in multiplicity. The double-diffractive events influence 1-and 2-bin. It is therefore hard to assume that the peak atn 3 ≃ 11 in the central region (|η| < 2.5) is caused by diffraction. The relevant mechanism could be the production of particles from fragmentation of (color) strings stretched between the leading partons of the colliding protons. The fragmentation of the longitudinally stretched strings is a soft process and results in narrow (nearly Poissonian) MD of particles with a limited range of transverse momenta.
The properties of the clan parameters in the first and the second NBD component point to a quite different kind of mechanism responsible for particle production. The large amount of branching in both of them can be understood by a two-stage mechanism; via collisions of constituents of the colliding protons with production of massive secondary objects (fireballs, clans' ancestors, etc.) in the first stage, and the successive branching of these objects into observable particles in the second one.
The events corresponding to the first, dominant multiplicity component are characterized by the production of a small number of clans containing many particles. The independence of their number,N 1 , from energy corresponds to the occurrence of a small number of mini-jets. Most particles in the clans carry relatively low transverse momenta. However, the weak suppression ofN 1 with p T indicates that the clans contain also particles with relatively high p T . The above-mentioned properties suggest that there are semi-hard processes with a small number of mini-jets followed by intensive branching and successive production of many low-p T particles in the dominant class of events. The second multiplicity component, under the tail of the MD, describes production of particles in semi-hard processes with extensive production of mini-jets. The corresponding NBD consists of many clans containing less particles on average than in the first component. We consider these clans to represent to a certain extent mini-jets, i.e. clusters of directed particles of a common ancestor. Such assumption is supported by the simultaneous increase of the average number of clans N 2 and the average number of particles per clann c2 with energy.

H. Summary
We have analysed charged-particle MDs using new high-statistics data [11][12][13] measured by the ATLAS Collaboration in pp collisions at energies √ s =8 and 13 TeV. The data include measurements in the low-p T regime (p T > 100 MeV/c) and with the cut p T > 500 MeV/c. The analysis extends our previous study [8] of the MDs obtained at the LHC at lower energies. The results of our investigations show that the new ATLAS data confirm the existence of a distinct peak in MD at low multiplicities up to the energy √ s =13 TeV. The description of the data within a two-component superposition of NBDs is unsatisfactory. The ATLAS data can be well parametrized by a weighted sum of three NBDs. The probabilities of the single NBD components show a weak or nearly no energy dependence in the most inclusive phase-space region (p T > 100 MeV/c). The third NBD component accounts for the peak-like structure of the experimental MD at low n. Its average multiplicityn 3 is approximately energy independent and decreases with the imposed transverse momentum cut. The parameters k −1 i , characterizing widths of the single NBDs, demonstrate a clear hierarchy with the index i for the p T > 100 MeV/c data sample. The third component of the total distribution is narrow, well described by the NBD with a large value of k 3 . The probability α 3 of the low-multiplicity component does not diminish when a cut on the transverse momentum is imposed.
We have studied the properties of the cumulative combinants N C i calculated from a weighted superposition of three-and two-NBD parametrization of the ATLAS and CMS data on MDs. It was demonstrated that the third NBD component at low multiplicities with sufficiently largen 3 can be responsible for the oscillating structure of the combinants. The oscillations are clearly visible when they are calculated from the three-NBD fits to the ATLAS data measured in the pseudorapidity interval |η| <2.5 with p T > 100 MeV/c. It was demonstrated that the oscillatory behaviour of N C i obtained from the three-NBD fits [8] to the CMS data on MDs at √ s =7 TeV in different pseudorapidity windows reveal similar properties as those found in [14]. They are the decrease of the magnitude and periodicity of the oscillations with the decreasing window size. The position of the first minimum of the oscillations corresponds to the average multiplicityn 3 of the third NBD component. The result suggests that the average correlation of the system is weaker if just this multiplicity contributes to the build-up of the probability P (n) for a higher multiplicity n. To the contrary, the cumulative combinants N C i calculated from a weighted superposition of two NBDs used to fit the data give no oscillations and reveal a monotonic decrease only.
We conclude that a detailed parametrization of the LHC data on MD which would give oscillations of the cumulative combinants N C i = (i+1)C(i+1) requires third component within weighted sum of NBDs. This is analogous to the oscillations of the ratio H q = K q /F q of the cumulant to the factorial moments which were described by two-NBD fits used to account for the shoulder of MDs observed already at lower energies. The oscillating structures corresponding to the low-and high-multiplicity regions cannot be obtained from the two-and one-NBD fits to the experimental data, respectively.
A modification of cumulative combinants based on a truncation of the MD in the first two multiplicity bins was examined. The modified combinants, (i+1)C 2 (i+1) and (i+1)C 1 (i+1), were obtained from experimental data on MD using Eqs. (8) and (9), respectively. It was illustrated that they reveal similar properties as the cumulative combinants calculated by Eq. (4) from a superposi-tion of three NBDs with the parameters obtained from the corresponding fits. The results of the analysis of the ATLAS data measured in the interval |η| <2.5 with the transverse momentum cut p T >100 MeV/c at √ s = 13, 8, 7 and 0.9 TeV show that the quantity (i+1)C 2 (i+1) reveals an oscillating structure in the region of low i. The wavy structure is clearly visible because of relatively small systematic uncertainties of the ATLAS measurements. The oscillation is a consequence of the distinct peak at the maximum of the MD. The oscillating wave at low i allows to discriminate between the two-and three-NBD parametrization of the ATLAS data.
We have computed the quantities (i + 1)C 1 (i + 1) and (i+1)C 2 (i+1) using data on MD measured by the CMS Collaboration in the pseudorapidity interval |η| < 2.4 at √ s=7 TeV. The results of the analysis show that the CMS data is inconclusive as regards the discussed oscillations (see Figs. 19a and 21a).
We have analyzed the clan structure of the MD of charged particles produced in pp collisions. The clan parameters were studied in the framework of the weighted superposition of three NBDs used to describe the ATLAS data on MD in the interval |η| <2.5 and p T > 100 MeV/c. The examination of the experimental data shows considerable differences in the behavior of the average number of clansN i and the average number of particles per clan n ci in comparison with two-component parametrization of the same data. The differences concerning the average values of the parameters are summarized as follows.
The three-NBD model gives:  3. more particles per clan in the first component than in the semi-hard one with mini-jets (n c1 >n c2 ).
The two-NBD model gives: It was shown that the above properties of clans are consistent with the results of the analysis performed with data on MD measured by the CMS Collaboration in different pseudorapidity intervals at √ s = 7 TeV. As concerns the third multiplicity component, the three-NBD model gives approximately constant value of the average number of clansN 3 with energy. These clans contain one or very few particles on average.
We have studied the influence of the transverse momentum cut on the clan parameters using the ATLAS data with p T > 500 MeV/c. The average number of clans and the average number of particles per clan are reduced in all three multiplicity components with the p T cut. It was shown that the reduction reveals a different behavior from the two-NBD model. The size of the changes suggests that there are semi-hard processes with simultaneous production of many low-p T particles in the first, dominant class of events. The values of the clans' parameters point to intensive branching and a small number of clans mini-jets in these events. On the other hand, the semi-hard component under the tail of the MD reveals the existence of many particle clusters with properties expected for mini-jets. The transverse momenta of particles in these clusters are on average higher than for the most particles belonging to the clans of the first component. The suppression of the clan parameters with the cut p T > 500 MeV/c indicates that the transverse momenta are concentrated around somep T . A novel feature is growth of the average number of the clans,N 2 , with √ s. As concerns the third multiplicity component, the present study shows that it consists mostly of particles with low transverse momenta. The biggest suppression ofn 3 with p T suggests that there is a soft mechanism responsible for particle production under the peak of the MD at low n.
Based on the performed analysis we argue that there is strong evidence of a new component in the MD of charged particles measured by the ATLAS Collaboration in pp collisions up to the energy √ s =13 TeV. The data are well described by a superposition of three NBD functions in a wide range of √ s. The third component at low n gives oscillations of the cumulative combinants N C i = (i + 1)C(i + 1) observed for the first time in [14]. The analysis of the clan properties applied to three-NBD fits to the ATLAS data on multiplicities reveals new features in the LHC energy domain. Further study of the clan characteristics and the fine structure of the MD with future high quality data is a challenging problem, which can give important information on production mechanism and correlation structures of the multi-particle states created in high energy protonproton collisions.

Appendix A
Here we summarize the details of the least-squares criteria and the conditions used in the minimization procedure. The value of χ 2 was computed as where M is the number of bins considered in the experiment. The bins are defined using the information accessible in the following form: the probability P ex m in the point n m and the bin interval (a m , b m ). In general, the values of n m are not given as integers. In the case when the mth bin consists of a single multiplicity n, n m = n and P ph m = P (n) with P (n) given by (2). For the bins containing more multiplicities, such as in the tails of the distributions, n m = n and we define  Tables I -IV. The obtained values of χ 2 /dof are quoted in the separate rows. The number of degrees of freedom, dof , is stated explicitly as the number of experimental points minus the number of free parameters.
In order to estimate how covariations in the data on MD can affect the results obtained by minimization of the expression (A1), we have performed calculations using the generalized formula The weight matrix V is the inverse of the covariance error matrix Σ = V −1 . The covariance matrix of the chargedparticle MD takes into account correlations between multiplicities as well correlations introduced by corrections applied on the data. The published MDs are obtained from data which suffer from low track reconstruction efficiencies ε trk , resolution effects etc., so that experimental errors in the resulting MDs are strongly correlated across adjacent multiplicity bins. The track efficiency uncertainty is dominant one for the n ch distributions. It results in positive correlations between the adjacent multiplicity bins and in anti-correlation between the low-and high-multiplicity tails. The variations of ε trk induce shifts of the maximum of the MD which, together with a proper normalization, increases the right tail of the distribution with a simultaneous decrease of the left tail and vice versa. As the leading uncertainty is due to the track efficiency systematics, the shifts correspond to the up/down systematic uncertainties of P n . Because systematic errors largely exceed the statistical ones, the effect of covariations between the latter is likely to be negligible compared to the systematic covariations. Taking all this into account, given the asymmetry of the systematic errors, we write the covariance matrix in the following way. Its diagonal part is constructed in the form The off diagonal part (i = j) is considered as follows The symbols σ st i and σ + i , σ − i > 0 mean the statistical and the systematic (up and down) errors in the ith bin, respectively. The value of m denotes multiplicity, which corresponds to the maximum of the measured distribution. The correlation factors ρ i,j = ρ c g i g j are governed by the overall correlation coefficient ρ c . As the track efficiency uncertainty is dominant one for the n ch distributions and because it results in positive correlations between the adjacent multiplicities, we consider ρ c > 0. Due to the alternating signs of Σ in different multiplicity quadrants defined in (A6), we introduce the functions which model the necessary interpolation of the shifted distributions in the vicinity of the point m. We used ∆ = 1 so that g i has a perceptible impact on the nearest neighbouring bins next to this point only. The definition of the matrix Σ satisfies the two requirements: positive covariations between adjacent multiplicity bins, and longrange negative covariations between bins on the opposite sides of the distribution maximum. The considered structure of the covariance matrix leads to a certain limitation on the parameter ρ c . The restriction follows from the requirement of the positive definiteness of the matrix Σ and from data on bin-to-bin errors in single measurements. The extra strong correlations (ρ c → 1) do not meet the criteria for valid correlation matrices [30]. Our calculations were performed with 80% correlations (ρ c = 0.8) for the ATLAS data with p T > 100 MeV/c and for the CMS data at √ s =7 TeV. We have used a smaller value of ρ c = 0.6 for the ATLAS data with the cut p T > 500 MeV/c. This is motivated by an increase of the track reconstruction efficiency parameter V: The parameters of the three-NBD superposition extracted from fits to data on MDs [10,11,13] measured by the ATLAS Collaboration in the window |η| < 2.5 with the cut pT > 100 MeV/c, n > 1. The parameter values were obtained by minimization of Eq. (A4) √ s i αini ki ε trk ≃ 0.75 [12] for the p T -cut data, which results in effective decrease of correlations between multiplicity bins. The change of ρ c affects the uncertainties of the obtained parameters but their values remain within the specified errors. We have checked the positive definiteness of the corresponding error matrices Σ in all considered cases. The error matrices were inverted to obtain the weight matrices V . The parameters of the weighted superposition of the three NBDs obtained by minimization of Eq. (A4) are shown in Tables V and VI. The shaded rectangles depicted in the respective figures represent the parameter values within the stated errors. The CMS data measured in the intervals |η| < η c = 0.5, 1.0 and 1.5 were fitted with the fixed value of k 3 = ∞. Except for some results of the fit to the ATLAS data with p T > 500 MeV/c at √ s = 7 TeV, one can see reasonable agreement between the parameters obtained by Eq. (A1) (full symbols with error bars) and by Eq. (A4) (shaded rectangles) for all fittings when using the superposition of three NBDs. The fits with the three-component function (2) seem to be rather stable with respect to covariations of the systematic uncertainties reflected by Eqs. (A5) and (A6) for the error matrix Σ.