Decomposition of multi-particle azimuthal correlations in Q-cumulant analysis

The method of Q-cumulants is a powerful tool to study the fine details of azimuthal anisotropies in high energy nuclear collisions. This paper presents a new method, based on mathematical induction, to evaluate the analytical form of the high-order Q-cumulants. The ability of this method is demonstrated via a toy model that uses the elliptic power distribution to simulate anisotropic emission of particles, quantified in terms of Fourier flow harmonics v_{n}. The method can help in studying the large amount of event statistics that can be collected in the future, and to allow measurements of very high central moments of the v_{2} distribution. This can in turn facilitate progress in understanding the initial geometry, input to hydrodynamic calculations of medium expansion in high energy nuclear collisions, as well as constraints on it.


Introduction
In ultra-relativistic collisions of nuclei at both the Relativistic Heavy Ion Collider (RHIC) [1][2][3][4] and the Large Hadron Collider (LHC) [5][6][7][8] a hot and dense system of strongly interacting quarks and gluons called the "quark-gluon plasma" (QGP) is created.One of the key observables used to study properties of the QGP is the anisotropic collective flow that is quantified by Fourier harmonics v n .Many methods have been developed to measure these harmonics [9][10][11].One of them, the method of cumulants based on multi-particle correlations [12,13], enables to suppress short-range correlations arising from jets and resonance decays and then reveals a genuine collectivity arising from the expansion of the QGP.The Qcumulant method [14] is an improved version of the initial cumulant method.Recently, the Q-cumulants have been employed to examine hydrodynamics predictions in Ref. [15] based on the ratio of the differences between adjacent cumulants [16].This analysis, performed using the CMS data, indicated the necessity of introducing moments higher than skewness to describe finer details of the elliptic flow v 2 distribution.In order to do so, higher-order cumulants have to be calculated.In this paper we present a way how to determine expressions needed to calculate higher-order cumulants.
The paper is organized as follows.Section 2 gives the basic quantities used in the Q-cumulant method.Section 3 gives a foundation of the method within a finite-dimensional vector space.Section 4 gives the results together with applications of the method to simple cases.Using toy models, we demonstrate the capabilities of the method in Section 5. A summary is given in Section 6. Tables that summarize coefficients needed for expressing the cumulants up to the 14-th order can be found in the Appendix.

Basics of the Q-cumulants
The general formalism of cumulants, for the purpose of flow measurements, was first used in Refs.[12][13][14].Cumulants were expressed in terms of moments of the magnitude of the corresponding flow vector.Such cumulant method is systematically biased due to trivial effects of same particle autocorrelations.To avoid the computation expensive "nested loops" in removing autocorrelations in the measurements of multi-particle azimuthal correlations, Refs.[12][13][14] have proposed an alternative method, which is based on the use of generating functions.An improved cumulant method, the Q-cumulants [14] allows, in principle, for a fast and exact calculation of all multi-particle cumulants.However, in practice the determination of the analytical expressions of multi-particle cumulants, which involve azimuthal correlations more than eight particles, is very difficult.
The Q-cumulants are built upon analytical expressions for multi-particle azimuthal correlations in terms of flow vectors Q, with: evaluated in different orders of Fourier harmonics, n.Here M is the multiplicity or number of particles in an event, and  j labels the azimuthal angle of the j-th particle measured in a fixed coordinate system in the laboratory.The method involves a decomposition of 2m-th power of the magnitude of the flow vector, |Q n | 2m : into (off diagonal) terms having 2m, 2m-1, … different indices up to the term where all 2m indices are equal (diagonal term).The first term, having 2m different indices, of decomposition is proportional to the 2m-particle azimuthal correlations m 2 : In the case of full decomposition, |Q n | 2m is expressed in readily calculable terms of powers of the flow vector given with Eq. ( 1) along with the anticipated 2m-particle azimuthal correlations introduced in Eq. (3).A thorough derivation with detailed examples is presented in Ref. [17].The analytical decomposition of |Q n | 2m for m > 4, however, becomes tedious.
In this paper we present a numerical method of decomposition of |Q n | 2m that, with the use of modern computers, enables one to easily obtain analytical expressions of multi-particle azimuthal correlations of higher orders.

Foundation of the method
The full decomposition of |Q n | 2m consists in partition of appropriate sets (or multisets in the case of repeated elements) of the azimuthal angles that lead to the partial Bell polynomials or their generalization (case of multisets).In derivation, some parts of the Bell polynomials have to be expressed through the lower-order multi-particle azimuthal correlations.This change of side in the expression makes a sign change of the corresponding term, resulting in coefficients of the partial Bell polynomials in the final form of the 2m-particle azimuthal correlations to be positive or negative integers.For example, in the case of m = 2 one obtains the four-particle azimuthal correlations might be presented as a linear combination of the corresponding Q-vector terms: with the corresponding integer coefficients (+1, -2, +1, -4, +2) in front of each of the terms.
The same holds true for all 2m-particle correlations but with different sets of integer coefficients.This fact inspired us to calculate all these integer coefficients by solving an appropriate system of algebraic equations.
The 2m-particle correlations m 2 might be considered as a member of a finite-dimensional vector space V d , where is a dimension of the vector space, and } , where x 1 , x 2 , … are unknown integer numbers that need to be obtained and ) , ( l m i f are known integer functions of multiplicity M. The left side of Eq. ( 6) can be calculated directly from Eq. (3) only in case of low multiplicity, being computationally manageable.The d unknown integer numbers,


on the right hand side of the Eq. ( 6) might be calculated numerically by solving a system of d linear algebraic equations presented in the next Section.However, in this paper we show straightforward determination of the basis by the method of mathematical induction, which enables easy calculation of the multi-particle azimuthal correlations of higher orders.First, it should be noticed, by a simple inspection of the published results [14,17], that each basis ), Re( , , Also, the basis   . Thus, the problem of determining the basis We found that the subset l = m contains different compositions (products) of the flow vectors each having a subscript (sub1, sub2,…) that corresponds to the well-known "partition of the positive integer" m.For example, in the case of the 8-particle (2m = 8) azimuthal correlations, the subset l = 4 contains these compositions of the flow vectors: where each composition contains subscripts that correspond to the partition of the integer 4.
Here we have used a convenient symbolic notation introduced in Refs.[14,17]: To obtain all basis vectors of the subset l = 4, each composition in (8) has to be combined by the ordered composition of the complex conjugate vectors as is shown in the following pattern: This pattern gives fifteen (p(4)[p(4)+1]/2 = 15) different combinations: The real part of each of these additional vectors (11), together with the basis 6 B , makes the complete basis of the 8-particle azimuthal correlations 8 B .This presentation of the composition of the flow vectors that are separated from the complex conjugate part by a vertical bar is inspired by the work of Ref. [14].
The corresponding multiplying integer functions ) , ( l m f might be also obtained via mathematical induction.These integer functions are given by: For example, in the case of 8-particle azimuthal correlations (m = 4): We list these integer functions and basis vectors for each of the 2m-particle azimuthal correlations ) in the third and fourth columns respectively of the corresponding Table 1-7.Now we have all necessary quantities to form the system of linear equations to obtain all the unknown coefficients for each of the 2m-particle azimuthal correlations: In order to set a solvable system of equations ( 14) the multiplicity has to be m M 2  (otherwise, the system of equations cannot be formed).For example, the system of equations for the 2-particle azimuthal correlations is given by ( One can notice that in the sum in the left-handed side of the Eq. ( 15) are, in fact, added a complex number and its conjugate, so the sum will be reduced to a real number.In order to set the system of equations (15) By using the same sets of azimuthal angles but now for Fourier harmonic n = 3, the system of equations (15) These two systems of linear equations are mutually different but both have the same solution: x so we fill in the second column of the Table 1 with the corresponding i x values.
The sets of the azimuthal angles might be chosen randomly, as we did in the upper example.That should secure that the calculated basis m B 2 do not have collinear vectors, otherwise the linear system of equations is not solvable.
The calculated numbers that enter the system of Eqs.(16,17) should be given with appropriate number of significant digits in order to get the obtained coefficients i x as true integers.There are no criteria to determine the appropriate number of significant digits.One can use the "trial and error" method to find out the required number.Otherwise, the obtained i x will not be integers and the rounding of i x might be problematic in the case of large basis m B 2 .The limits of the application of this method depend on the number of significant digits by which the computer operates.We have successfully applied the method up to the 14-th order of cumulants.
The two-particle azimuthal correlation formula can be formed by simply reading the corresponding values from Table 1: The similar table can be obtained for the four-particle azimuthal correlations by the solution of the corresponding linear system of equations: That might be also easily rewritten in an appropriate formula by simply reading the corresponding values of The appropriate tables for the higher order cumulants are given in the Appendix.
Many of the basis vectors are complex valued numbers like: and one should take only their real part when one writes the expression of the 2m-particle azimuthal correlations.Actually, the analytical derivation reveals that in addition to these complex valued basis vectors their complex conjugates ( ) also participate in the corresponding basis.Because of symmetry, they always enter the expression of the 2m-particle azimuthal correlations with the same coefficients and so their imaginary parts cancel.This also causes the corresponding coefficients to be even numbers.For example, ) Re( 2 , and so: is always even number.Odd coefficients might appear only in front of basis vectors that are not accompanied with their complex conjugates like Q in Eq. (19).
The coefficients have also other interesting features that should be noticed.For example, all coefficients which correspond to the same subset of basis vectors (having the same index l ) sum up to zero: (Table 3.: ), (Table 4.: , and so on).The exceptions of this rule are the two coefficients at beginning of each where d is the dimension of the vector space, e is Euler"s number, and is the incomplete gamma function of x with parameter a .
We also applied this method of derivation for the case of the weighted Q-vector evaluated in the harmonic n : where j w is the weight of the j-th particle [17].The derivation in this case is more complicated primarily because there are more basis vectors involved.However, the obtained coefficients in front of the basis vectors have very simple features: (The dimension of the vector space in this case is not equal to the previous one: Finally, having expressions for the 2m-particle azimuthal correlations, 2m one can calculate the weighted average over all events 2m (given by Eq. ( 1) in Ref. [18]).Then, one can use the recurrence relation to calculate the cumulants of any order by knowing all the cumulants of lower orders [18]: The cumulant based flow harmonics v n {2k} (k = 1, 2…) can be calculated by the following equation [16]: where coefficients a 2k are obtainable by recursion relation: which enables easy calculation of high-order v n {2k} by using any commercial program for ) obtained in the experiment [16], indicates the non-zero value of the kurtosis  40 , which is the fourth central moment of the v 2 distribution.So, to obtain the additional fourth central moments,  40 ,  22 ,  04 [16], it requires experimental values of at least three more cumulants of higher order Therefore, in order to show the validity of the above described method, the obtained expressions for the cumulants up to the 14-th order are calculated with the azimuthal angles simulated using toy models.For each event of a given v 2 a simple distribution is used to generate particle azimuthal angle.We set the input value 15 .0 2  v with no event-by-event flow fluctuations.Thus, one expects that all v 2 {2k} are equal to the input value.Fig. 1 shows that this is indeed what is seen as there is excellent agreement between the cumulant based v 2 {2k} estimations and the input v 2 value of 0.15.Due to a strong correlation between different 2m [18], the statistical uncertainties of the v 2 {2k} very slowly increases with the increase of k.This toy model therefore validates the method.
FIG. 1 Cumulant based estimations of the elliptic flow v 2 from toy Monte Carlo simulations.The input value of v 2 is 0.15.All cumulants retrieve this input value with high precision.
However, from this validation the sensitivity to flow fluctuations of the higher-order cumulants cannot be verified.Therefore, we applied the second toy model where the initial eccentricity 2  distribution is simulated using the elliptic power distribution: where  and 0  are power and ellipticity parameters respectively, which take different values, obtained by Glauber model for 5.02 GeV Pb-Pb collisions, depending on the centrality [19].The scaling factor 2  between the elliptic flow and the initial eccentricity, , is chosen to imitate the centrality dependence of the elliptic flow v 2 measured in Ref. [16].For each event of a given v 2 a simple distribution is used to v 2 {2k} k generate particle azimuthal angle.The integral in Eq. ( 27) can be carried out analytically to give [20]: However, the ROOT version [21] of the hypergeometric function in Eq. ( 28) is not defined everywhere in the interval (0, 1) of 2  .In Fig. 2 It can be seen that only the centrality range 0-15% of the 5.02 GeV Pb-Pb collisions might be well simulated by this hypergeometric function.We solved this inconvenience by applying the Pfaff transformation: For each centrality, about 1.5 x 10 7 events have been simulated using the above described toy model.Figure 3 shows v 2 {2k} values (k = 1,…,7) calculated based on the obtained expressions of the corresponding 2m correlations as a function of centrality.A gap between the v 2 {2}and higher-order cumulant values v 2 {2k} (k = 2,…,7) is present.It is due to flow fluctuations that relates higher-order cumulants based v 2 {2k} and the variance  , for k > 1 [15].As it has been seen from the experimental data [16], flow fluctuations become larger going to peripheral collisions.The elliptic flow v 2 values are reproduced using the expressions for the 2m developed in this paper.
FIG. 3 The ) as a function of centrality calculated from the data simulated by the elliptic power distribution toy model.The statistical uncertainties are negligible compared to the marker size.
As the fluctuations in the initial state are not Gaussian, the v 2 {2k} values for k > 1 will not be the same.This will produce a splitting between different v 2 {2k} values and they will be ordered as v 2 {2k} > v 2 {2(k + 1)} for any k > 1.In Fig. 3 the ordering and the fine splitting Centrality (%) between the v 2 {2k} values (k = 2,…,7) are not well visible.In order to make the splitting between the cumulants of different orders explicitly visible, we show in Fig. 4 the relative differences ) as a function of centrality.The points in Fig. 4 depicted with different symbols for different cumulant orders are calculated from the simulated data.On the other hand, the corresponding input values, obtained directly applying the elliptic power distributions, are represented by spline interpolation lines.In Fig. 4 the ordering is clearly seen, as well as the fine splitting between the cumulants of different orders.The relative difference between the cumulants decreases by about one order of magnitude for each increment of the order k.An excellent agreement between the lines and the symbols prove correctness of the obtained expressions for the 2m-particle azimuthal correlations 2m .

FIG. 4 The relative differences
) as a function of centrality.The input values, obtained directly applying the elliptic power distributions, are represented by spline interpolation lines.

Conclusions
In this paper we present a method, based on mathematical induction, that allows to calculate the high-order v 2 {2k} values from the Q-cumulants in a relatively straightforward way.The Relative Difference Centrality (%) analytical expressions for high 2m-particle correlations, 10 , 12 , and 14 , are obtained for the first time and given here in the form of table values.The validity of the proposed method is confirmed by the elliptic flow simulation with a toy model that exploits elliptic power distribution.We transformed the hypergeometric function, involved in the elliptic power distribution, by applying Pfaff transformation, and enabled its calculation in the ROOT program for all parameter values of the v 2 distribution.The ability to calculate high-order v 2 {2k} values allows the possibility to study the fine details of the v 2 distribution by extracting its skewness and higher central moments.Theoreticians can tune initial states in their hydrodynamics models such to reconstruct measured central moments.This can place stringent constraints on the initial geometry used in hydrodynamic calculations of the collective dynamics of the QGP in high energy nuclear collisions.
Coefficients, integer functions and basis vectors for calculation of the eight-particle azimuthal correlations.
Coefficients, integer functions and basis vectors for calculation of the ten-particle azimuthal correlations.
Coefficients, integer functions and basis vectors for calculation of the twelve-particle azimuthal correlations.

.
For example, the basis of 4particle azimuthal correlations contains the complete basis of the 2-particle azimuthal correlations subset l = m. i i

5 20  and 2 02
Demonstration of the method using a toy model Experimental values of the v 2 {2k} enable determination of the central moments of the v 2 distribution.As a way to obtain the lowest central moments of the v 2 distribution as are variances 2 , skewness 30 s , and co-skewness 12 s [15, 16], it requires experimental values of at least four different cumulants v 2 {2}, v 2 {4}, v 2 {6}, and v 2 {8}.However, the centrality dependence of the hydrodynamics probe h FIG. 2 Parameter values (,  0 ) of the elliptic power distribution for 5.02 GeV Pb-Pb collisions calculated in the Glauber model for different centralities, and the area of definition i

Table 1
Coefficients, integer functions and basis vectors for calculation of the two-particle azimuthal correlations.

Table 2
Coefficients, integer functions and basis vectors for calculation of the four-particle azimuthal correlations.

Table 2 :
One more interesting feature that might be used to check the correctness of the obtained expression of the 2m-particle azimuthal correlations is a sum of the absolute values of all coefficients in a Table.This sum reads: Table that correspond to the basis vectors 1 and 2 | |Q n .The coefficient in front of the basis vector "1" is m  .

Table 3
Coefficients, integer functions and basis vectors for calculation of the six-particle azimuthal correlations.