Event-by-event cumulants of azimuthal angles

We further develop the recently proposed event-by-event cumulants of azimuthal angles. The role of reflection symmetry, permutation symmetry, frame independence, and relabeling of particle indices in the cumulant expansion is discussed in detail. We argue that mathematical and statistical properties of cumulants are preserved if cumulants of azimuthal angles are defined event-by-event in terms of single-event averages of azimuthal angles, while they are violated in the traditional approach in which cumulants are defined in terms of all-event averages. We derive for the first time the example analytic solutions for the contribution of combinatorial background in the measured two- and three-particle correlations. We demonstrate that these solutions for the combinatorial background are universal, as they can be written generically in terms of multiplicity-dependent combinatorial weights and marginal probability density functions of starting multivariate distribution. The new general results between multiparticle azimuthal correlators and flow amplitudes and symmetry planes are presented.


Introduction
In ultrarelativistic heavy-ion collisions, the properties of an extreme state of nuclear matter, dubbed quark-gluon plasma, can be studied in a controlled environment [1][2][3]. Among various experimental observables used to probe its properties, in the present work we focus exclusively on multivariate cumulants [4]. They have been regularly applied at both the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) experiments for the measurement of anisotropic flow phenomena [5,6] and in femtoscopic studies [7,8]. In the context of flow analyses, the cumulant expansion has been performed traditionally on azimuthal angles [9][10][11][12], and only recently also directly on flow amplitudes [13][14][15]. Both a e-mail: ante.bilandzic@tum.de (corresponding author) the statistical and systematic uncertainties of these measurements are suppressed with inverse powers of multiplicity, which means that multivariate cumulants are a precision tool only in heavy-ion datasets characterized by large values of multiplicity. On the other hand, their measurements in the collisions of smaller objects, such as proton-proton collisions, are unreliable, and there is still no consensus in the field on how to interpret them in this environment.
The mathematical properties of cumulants are rigorous and well established (for a recent concise review and summary, see Ref. [15]). The non-trivial physics information which can be extracted from cumulants stems from the observation that a cumulant is identically zero if one of the variables in it is statistically independent of the others [4]. However, this statement is not a mathematical equivalence, and therefore Theorem 1 from Ref. [4] cannot be used to conclude that variables in the cumulant definition are statistically independent of each other if the cumulant is zero. On the contrary, a cumulant is not zero if and only if all variables in it are statistically connected [4], i.e. there exists a genuine multivariate correlation among all variables in the cumulant which cannot be written as a superposition of correlations among a smaller number of variables. By definition, each higher-order cumulant extracts a piece of new and independent information that is not accessible at lower orders.
The traditional cumulant expansion over azimuthal angles yields the final expressions that contain only the isotropic all-event averages of azimuthal correlators [9,10]. All nonisotropic averages are identically zero due to random fluctuations of the impact parameter vector for detectors which have uniform azimuthal acceptance [9,10]. However, it was recently realized that these results for cumulants of azimuthal angles in which the non-isotropic azimuthal correlators are missing no longer satisfy the fundamental properties of cumulants [13,15]. The reason for this failure is that mathematical and statistical properties of cumulants are in general preserved only if all terms in the cumulant expansion are kept. As soon as there are underlying symmetries due to which some terms in the cumulant expansion are lost, their mathematical and statistical properties are invalidated, while their physical interpretation can be completely different.
An alternative definition of cumulants of azimuthal angles was recently proposed in Ref. [15]. Instead of defining cumulants in terms of all-event averages of azimuthal angles, cumulants are defined in terms of single-event averages. These new event-by-event cumulants of azimuthal angles satisfy all mathematical and statistical properties of cumulants. However, their direct application and measurement are plagued by contributions from the combinatorial background, which are very difficult to remove completely. The standard and approximate mixed-event technique used in high-energy physics to remove a combinatorial background is not applicable for azimuthal angles.
The main motivation for the current draft is an attempt to make further progress with these recently introduced eventby-event cumulants of azimuthal angles. As our main new result, we show how the contributions from the combinatorial background in the measured two-and three-particle correlations can be solved analytically by using universal combinatorial weights which depend only on multiplicity and using marginal probability density functions (p.d.f.). These example solutions can be straightforwardly generalized to correlations involving more than three particles and for cases when multiple particles were misidentified in the calculated correlators. The generalization will be presented in our future work.
The rest of the paper is organized as follows: In Sect. 2, we discuss the role of various underlying symmetries which can appear in cumulant expansion. In Sect. 3, we discuss how the presence of a combinatorial background affects measurements obtained with correlation techniques. The new set of general results between multiparticle azimuthal correlators and flow amplitudes and symmetry planes are presented in Appendix A.

Role of symmetries
In this section we discuss the role of a few symmetries which can appear in the cumulant expansion. As our main conclusion, we demonstrate that genuine multivariate correlations can be reliably extracted with cumulants only if there are no underlying symmetries due to which some terms in the cumulant expansion are identically zero. This topic is of great relevance particularly for anisotropic flow analyses in high-energy physics, where most of the cumulants used now suffer from this problem. In what follows next, we will frequently use the phrase "p.d.f. does not factorize," with the following meaning: The starting nvariate (or joint) p.d.f. f (x 1 , x 2 , . . . , x n ) cannot be decom-posed into the product of marginal p.d.f.'s corresponding to some set partition of x 1 , x 2 , . . . , x n . For instance, a factorizable three-variate p.d.f. would permit the decompositions where f x i x j and f x i are two-and single-variate marginal p.d.f.'s of the starting three-variate p.d.f. f (x 1 , x 2 , x 3 ). We use indices for p.d.f.'s to indicate that their functional forms can be different, instead of naming them differently. For ease of notation, we drop indices only when we refer to the starting n-variate p.d.f. of all x 1 , x 2 , . . . , x n variables in question. Theorem 1 from Ref. [4] states that the n-variate cumulant is identically zero if the underlying n-variate p.d.f. is factorizable. Physically, this means that there are no genuine correlations among all n particles in the system, i.e. n-particle correlation is just a trivial superposition of a correlation involving less than n particles. Mathematically, it is possible to find a set partition of x 1 , x 2 , . . . , x n in which subsets are statistically independent of each other. We now demonstrate that, in general, the opposite is not true, i.e. the cumulant can be identically zero because of underlying symmetries and not because variables on which the cumulant expansion has been performed are statistically independent.

Reflection symmetry
We start the discussion with the following well-known argu- , is a two-variate p.d.f. which does not factorize, the corresponding twovariate cumulant is not zero. Based on this observation, we conclude that the two variables x 1 and x 2 are not statistically independent, i.e. there exists a genuine twobody correlation between them. This conclusion fails, however, if in addition the starting p.d.f. f (x 1 , x 2 ) has the following reflection symmetry: Because of this symmetry, the corresponding two-variate cumulant, x 1 x 2 − x 1 x 2 , is now identically zero, simply because both x 2 ) dx 1 dx 2 are identically zero. But this result does not imply that x 1 and x 2 are statistically independentthe starting p.d.f. is still not factorizable. This simple example clearly illustrates the failure of cumulants in the presence of underlying symmetries [16].
We now generalize this argument to higher orders. If f (x 1 , x 2 , . . . , x n ), x 1 , x 2 , . . . , x n ∈ (−∞, ∞), is an nvariate p.d.f. which does not factorize, and if we again assume that there is a reflection symmetry in one variable , it still follows that identically x 1 = 0, since The above derivation holds true if we add some of the remaining variables x j , x k , . . . into the correlator, i.e. we have Each term in the cumulant expansion is a possible partition of the set x 1 , x 2 , . . . , x n , in which in addition the averages have been taken over each subset. Therefore, in each term in the cumulant expansion there is always exactly one average which involves x 1 , which is identically zero if the underlying n-variate p.d.f. has a reflection symmetry . .) in that variable. From this we can conclude that due to reflection symmetry in any of the variables, cumulants at any order can be trivial and identical to zero. This does not imply that variables x 1 , x 2 , . . . , x n are statistically independent; instead, this is an extreme example in which the cumulant expansion is invalidated at all orders due to underlying symmetries.
To resolve between the two possibilities, the following additional cross-check needs to be performed in practice: If the underlying n-variate p.d.f. factorizes, cumulants are identically zero, irrespective of the choice of sample space for n variables x 1 , x 2 , . . . , x n . In this case, cumulants will be identically zero for any choice of boundaries in x 1 ∈ (x 1,min , x 1,max ), x 2 ∈ (x 2,min , x 2,max ), etc., and it can be safely concluded that there are no genuine multivariate correlations among all n variables. On the other hand, if cumulants are identically zero due to reflection symmetry, that will occur only for some specific choice of boundaries in

Permutation symmetry
Next, we consider the role of permutation symmetry in the starting n-variate p.d.f. f (x 1 , x 2 , . . . , x n ). For the sake of clarity, we present the argument for the simplest two-variate case, since the generalization to higher orders is trivial.
If for the starting two-variate p.d.f. we have that f (x, y) = f (y, x), and if the sampling spaces of x and y are identical, it follows immediately that marginal p.d.f.'s are the same, since This conclusion is used frequently later in the paper.

Frame independence
A mandatory requirement which any physics observable must fulfill is its independence from the rotations of a coordinate system in the laboratory frame. We can inspect the role of rotations already at the level of two-particle cumulants of azimuthal angles. In the traditional approach, two-particle cumulants of azimuthal angles are defined as [9,10]: The double angular brackets indicate that in the first step, averaging is performed over all distinct combinations of two azimuthal angles within the same event, and then these results are averaged over all events. If we now rotate randomly by angle α from one event to another the coordinate system in which azimuthal angles are measured, then both e inϕ 1 and e −inϕ 2 are trivially averaged to zero [9,10]. The reason is the following mathematical identity: In practice, the above identity is a direct consequence of random event-by-event fluctuations of the impact parameter vector. Only the first term in Eq. (4) is invariant with respect to the rotations. In this traditional approach, the rotational invariance of cumulants was achieved by simply averaging out all non-isotropic terms. However, by following such a procedure, especially at higher orders, most of the terms in the cumulant expansion are lost, and the final expressions are no longer valid cumulants of azimuthal angles. If we change the definition and instead define cumulants of azimuthal angles event-by-event [15] the rotational invariance is satisfied by definition, and all terms in the cumulant expansion are kept. We now have that This argument generalizes trivially to higher orders. As our new general result, we argue that all terms in the cumulant expansion can be kept, and simultaneously the rotational invariance of cumulants can be satisfied, only if cumulants of azimuthal angles are defined event-by-event in terms of single-event averages. This conceptual difference in the definition of cumulants of azimuthal angles (all-event averages vs. single-event averages) yields different final results and therefore has non-trivial consequences.

Relabeling
Another common practice in the field is to group together in the final expressions for cumulants all azimuthal correlators which are identical, apart from relabeling of indices of azimuthal angles. For instance, in the derivation of the traditional expression for a four-particle cumulant, there is the following intermediate result [9,10] All two-particle correlators in the above equation are estimated from the same sample of azimuthal angles. Therefore, they are mathematically identical, and to reflect that, they are grouped together to obtain the final expression We now argue that such a practice conflicts with the strict mathematical properties of cumulants. In general, the nvariate cumulant κ ν 1 ,...,ν n is by definition a sum in which each individual summand can be thought of as one possible partition of the set consisting of n starting variables {x 1 , ..., x n } [4]. After relabeling variable indices in the final expressions for cumulants, this fundamental property of cumulants is lost. This problem is closely related to the general problem of combinatorial background, which we tackle in the next section. Relabeling discussed here amounts essentially to ignoring all the non-trivial effects of the combinatorial background. The formalism of cumulants can be easily applied only in cases when combinatorial background plays no role (e.g. in the study of event-by-event fluctuations of different flow magnitudes). However, in cases when the combinatorial background is not under control, the usage of cumulants is not that straightforward.

Combinatorial background
In this section, we investigate the role of combinatorial background in analyses that rely on the use of correlation techniques and cumulants. We keep the discussion general, and therefore our main conclusions are applicable to diverse fields of interest (anisotropic flow, femtoscopy, etc.).

Trivial multiplicity scaling
One of the well-known results in high-energy physics is that in the absence of collective phenomena, the two-particle correlators exhibit trivial multiplicity dependence, the so-called ∼ 1/M scaling. In fact, the breakdown of this scaling is typically proposed as a very robust signature of the collective behavior of nuclear matter that is produced.
We now decipher this trivial multiplicity scaling and demonstrate that it does not originate solely from few-particle correlations, but rather from the interplay between combinatorial background and few-particle correlations. This is a subtle yet important difference. The importance of the combinatorial background has been regularly ignored in this context, but we now demonstrate that the contributions from fewparticle correlations do not exhibit this trivial multiplicity scaling in the measured two-particle correlators per se-the universal nature of ∼ 1/M scaling for vastly different physical sources of few-particle correlations instead stems naturally from the presence of the combinatorial background.
In general, in the measurements dominated by few-particle correlations, the nonflow contribution, δ k , in the k-particle correlator, k , exhibits the following universal scaling as a function of multiplicity M [9,10]: The derivation of this scaling is purely probabilistic and can be derived with the following simple argument. For the sake of simplicity, we outline the argument for four-particle correlators, but it can be trivially generalized to any higher-order correlator. In an event with multiplicity M, we consider four particles which are correlated through direct few-particle correlations which are not of collective origin (for instance, these four particles originate from the same resonance decay). The probability of selecting precisely these four particles into the four-particle correlator, out of M particles available in an event, is given by After we have taken the first particle into the four-particle correlator to be the one from the specific four-particle resonance decay, we have three particles left from the same resonance decay out of M − 1 in total; after we have selected the next one, there are two particles left from the same resonance decay out of M − 2 in total, etc. If we would have instead considered the five-particle resonance decays but still used four-particle correlators in the measurements, only the numerical factors in the numerator of Eq. (11) would change, and the multiplicity dependence in the denominator would remain the same. As a consequence, the universal multiplicity scaling in Eq. (11) is always present if the measured fourparticle correlators are dominated by contributions from fewparticle correlations, and does not depend on the details of few-particle correlations.
Only if there are collective effects which induce correlations among all produced particles is the universal multiplicity scaling in Eq. (10) broken. By following the same probabilistic argument, we have instead that the contribution of collective correlations among all M particles in an event, in the measured k-particle correlator, is We note that there is no room for ambiguity here. For instance, in the context of anisotropic flow analyses, theoretically we can have the fine-tuned possibility that a flow harmonic itself has the following multiplicity dependence: If we now sample all azimuthal angles individually from the single-variate Fourier-like p.d.f. parameterized solely with that harmonic (we use the standard techniques and notation from Ref. [12], see also Sect. 3.1.1), the two-particle azimuthal correlator yields which is the same result as a trivial multiplicity scaling in two-particle correlators in Eq. (10). The relation between any azimuthal correlator and flow harmonics can be obtained from the general procedure outlined in Appendix A. To distinguish these two vastly different possibilities, we must look at higher orders. For a four-particle azimuthal correlator, 4 , a contribution from few-particle correlations gives the following generic scaling 4 = δ 4 ∼ 1/M 3 (see Eq. (10)), while the non-trivial dependence of flow harmonics on multiplicity in Eq. (13) gives 4 ∼ 1/M 2 . Therefore, we can establish the following argument: If we observe 2 ∼ 1/M, we cannot conclude based solely on that measurement whether the two-particle correlator is dominated by few-particle correlations or by correlations which are of collective origin. However, at higher orders (k > 2, k is even) we have the splitting between two possibilities: Fig. 1 In this study, the elliptic flow v 2 itself exhibits ∼ 1/ √ M multiplicity scaling, which yields ∼ 1/M scaling in the two-particle correlators (blue markers), which is typically attributed to few-particle correlations. We cannot conclude based solely on the measurement of the two-particle correlator 2 whether this is a contribution from collective or non-collective effects. However, the four-particle correlator 4 (red markers) can discriminate between these two possibilities, because it follows what would be the collective flow scaling 1/M 2 in this example (red line) and not a universal few-particle nonflow scaling 1/M 3 (black line). To increase visibility, both horizontal and vertical axes are plotted on log scale 1. k ∼ 1/M k−1 ⇒ the k-particle correlator is dominated by an interplay between few-particle correlations and combinatorial background; 2. k ∼ 1/M k/2 ⇒ the k-particle correlator is dominated by collective effects.
We have restricted the above argument to even k, because for odd k there is an additional contribution from symmetry planes (see Appendix A), which is not relevant for the discussion presented here. We support this conclusion with the following toy Monte Carlo study.

The 1/M scaling in 2-p correlators: flow or nonflow?
As the starting point in this toy Monte Carlo study, we sample the multiplicity M of an event. From the sampled multiplicity, we determine the elliptic flow in that event via v 2 (M) = 1/ √ M, so that by definition in this study, 2 = v 2 2 ∼ 1/M, and 4 = v 4 2 ∼ 1/M 2 . In the next step, we sample particle azimuthal angles from Fourier-like single-particle p.d.f. f (ϕ) = 1 2π (1+2v 2 (M) cos(2(ϕ−Ψ ))), where Ψ is an event-by-event randomized symmetry plane (this is the standard procedure in the field and it resembles the random event-by-event fluctuations of impact parameter vector). Finally, from the sampled azimuthal angles we calculate the correlators 2 and 4 with the generic framework [12] as a function of multiplicity. The results of this study are presented in Fig. 1.

Example solution for combinatorial background for two-particle correlations
While the discussion in the previous section on multiplicity scaling in the measured correlators was only qualitative and purely phenomenological, we now provide the robust mathematical treatment. In particular, we establish the mathematical procedure to obtain the analytic solution for the combinatorial background for the example process in which two particles are emitted pair-wise from the most general joint two-variate p.d. f. f xy (x, y). In the next section we present the solution for three-particle correlations in order to demonstrate that the same mathematical procedure is applicable at higher orders. For the sake of clarity, we consider only the most important cases of particle misidentification in the measured correlators. The general solutions which cover all possible cases of particle misidentification will be presented in our future work. The starting problem can be formulated mathematically as follows: as a "signal" and the p.d.f. w(x, y) after randomization as a "signal + background." The general two-particle correlator which corresponds only to the signal is x y S = x y f xy (x, y) dxdy, while the one that corresponds to both signal and background is x y S+B = x y w(x, y) dxdy.
We now demonstrate that the relation between f xy (x, y) and w(x, y) is universal, in the sense that whatever the starting functional form of f xy (x, y) is, the same generic equation determines the functional form of w(x, y). That equation involves only combinatorial weights that depend only on multiplicity M and marginal p.d. f.'s of f xy (x, y). For instance, the marginal p.d.f. f x (x) is the probability of obtaining x whatever value of y, and it can be obtained from the following definition f x (x) ≡ f xy (x, y) dy. Further discussion about fundamental properties of marginal p.d.f.'s can be found in Ref. [16].
When particles are sampled pair-wise, and when only the most important case of particle misidentification is considered, without loss of generality we can divide the final dataset of M particles after randomization into the following four disjoint subsets: 1. Diagonal: Both particles originate from the same sampling. Only these two particles can be physically correlated, i.e. this is the signal and is therefore described with the starting p.d. f. f xy (x, y). Its combinatorial weight is

Misidentified:
To account for the fact that identical particles are indistinguishable, we must allow the possibility that the particle which was originally sampled as y was reconstructed as x, and vice versa. For instance, if two pions are emitted from two different processes, after randomization we do not even know in principle which pion originated from which process. This case is governed by two possibilities: f x (x) f x (y) and f y (y) f y (x). In both cases, the combinatorial weight is the same and it reads In the above expressions, f x and f y are marginal p.d.f.'s of x and y, respectively. In total, there are four disjoint subsets, and as a cross-check we can immediately see that for their respective combinatorial weights it follows that We can now write the example solution for the p.d.f. of a randomized sample The combinatorial weights p A , p B and p C are universalthey depend only on multiplicity M and not on any details of the starting p.d. f. f xy (x, y). The universal multiplicity dependence of these two-particle combinatorial weights is presented in Fig. 2. We stress that the example analytic result in Eq. (16) can be further generalized by including the contributions when misidentified particles corresponding to the Fig. 2 Multiplicity dependence of three distinct combinatorial weights, p A , p B and p C (see the main text for their definitions and explanation), which always appear in two-particle correlations. The contribution from genuine two-particle correlations (signal) is weighted with p A (green line). For large M this contribution is suppressed when compared with two distinct contributions from a combinatorial background (red lines). The results from this figure cover any case in which particles are produced in pairs from the general two-variate p.d.f. f xy (x, y) same initial particle type are coupled directly to each other in the two-particle correlators.
The practical use case of the above result amounts to the fact that we can now decompose and study separately the individual contributions to the measured two-particle correlators in the randomized sample and deduce their relative importance. For instance, the standard two-particle azimuthal correlation used in flow analyses can be calculated in the randomized sample as follows: where the p.d.f. w(x, y) is given by Eq. (16), while the azimuthal angles ϕ 1 and ϕ 2 can in general be taken from different sample spaces. The above equation describes analytically what is measured experimentally when we use twoparticle azimuthal correlations, for the case when particles are emitted pair-wise, when the final dataset is randomized, and when only the most important cases of particle misidentification are taken into account. Given the nature of the example expression for the p.d. f. w(x, y) in Eq. (16), we are naturally led to decompose distinct individual contributions to the two-particle correlator as follows: where and f ϕ 1 (ϕ 1 ) ≡ f ϕ 1 ϕ 2 (ϕ 1 , ϕ 2 ) dϕ 2 and analogously f ϕ 2 (ϕ 2 ) ≡ f ϕ 1 ϕ 2 (ϕ 1 , ϕ 2 ) dϕ 1 . We illustrate and support this general decomposition with the toy Monte Carlo study in the next section.

Toy Monte Carlo study for the two-particle case
We test the validity of Eq. (16) with the following clear-cut toy Monte Carlo study, which can be solved analytically. The starting normalized two-variate p.d.f. is where M is a parameter and corresponds to multiplicity. The stochastic variables are azimuthal angles ϕ 1 and ϕ 2 , whose sample space is [0, 2π). One can easily check that for any value of multiplicity M. With straightforward calculus, and after recalling that the combinatorial weight of the signal contribution for the twoparticle case is p A = 2! M 2 M(M−1) , we have analytically obtained for flow harmonic n = 2 for the signal contribution: and similarly for the three distinct background contributions 2 w 2 , 2 w 3 and 2 w 4 . All results are shown in Fig. 3. The blue markers show the experimental result, i.e. the average 2 = cos[2(ϕ 1 − ϕ 2 )] measured in the randomized data sample, after particle azimuthal angles were sampled pairwise from the p.d.f. in Eq. (20). The analytic expression for the signal contribution from Eq. (22) is shown with the green curve. Red curves represent analytic results for three distinct cases of combinatorial background contribution, 2 w 2 , 2 w 3 and 2 w 4 , respectively. Finally, the blue curve is a superposition of the signal (green curve) and distinct background contributions (three red curves). In this section, we address the combinatorial background for three-particle correlations. The problem we try to solve can be mathematically formulated as follows: Physically, f xyz (x, y, z) describes the "signal," while w(x, y, z) describes the "signal + background." For instance, if we are interested in the average threeparticle correlation x yz , then the signal contribution is given by x yz S = x yz f xyz (x, y, z) dxdydz, while the contribution from both signal and background is given by x yz S+B = x yz w(x, y, z) dxdydz. Given that we have only randomized the initial sample, there should be no systematic biases involved, i.e. we seek the universal relation between w(x, y, z) and f xyz (x, y, z). The derivation of that universal relation is the main result of this section-its generic structure deciphers the role of combinatorial background after we have randomized the final data sample for any starting p.d. f. f xyz (x, y, z).
We keep the discussion both general and as close to the experiment as possible. For instance, we introduce particle misidentification by blinding the particle labels x, y, z in the final data sample. This corresponds to the real-case scenario when, for instance, three pions are emitted from three different physical processes. After randomization, we cannot deduce even in principle which pion originated from which process, since identical particles are indistinguishable. Just as in the two-particle example, we consider here only this most important case of particle misidentification. The more general solutions which cover all possible cases of particle misidentification will be presented in our future work.
If particles are sampled in triples from f xyz (x, y, z), and if we consider only the most important case of particle misiden- Fig. 4 Multiplicity dependence of six combinatorial weights for threeparticle correlations: p A , p B1 , p B2 , p C1 , p C2 and p C3 (see the main text for their definitions and further explanation). The contribution from genuine three-particle correlations (signal) is weighted with p A (green line), and already for M ≈ 50 this contribution becomes negligible. The contribution from genuine two-particle correlations is also suppressed for large multiplicities (solid red line). The remaining weights, p C1 , p C2 and p C3 , all correspond to independent particle emission, and their relative contribution increases as multiplicity increases (c) f x f y f z . Finally, this case corresponds to the situation when we have three different particles from three different samplings in the three-particle correlator. The com- In total, in this example there are 20 disjoint cases to deal with: one from the signal and 19 from the combinatorial background. We can check immediately that the corresponding probabilities sum up correctly: From the above decomposition we can now read off and write down our final solution for this particular example. If the starting p.d.f. is f xyz (x, y, z), the p.d.f. w(x, y, z) after randomization is given by the following universal result: This result is universal because the six combinatorial weights p A , p B 1 , p B 2 , p C 1 , p C 2 and p C 3 depend only on multiplicity, while all marginal p.d.f.'s are always by definition calculated in the same way from the starting p.d.f. y, z) dydz, etc.). Therefore, whatever the starting p.d. f. f xyz (x, y, z) is, one needs to calculate all marginal p.d.f.'s and insert them in the above equation, and the result for the p.d. f. w(x, y, z) of the randomized sample follows. We stress again that the example analytic result in Eq. (24) can be further generalized by including the contributions when misidentified particles corresponding to the same initial particle type are coupled directly to each other in the three-particle correlators. Figure 4 shows the multiplicity dependence of six combinatorial weights for three-particle correlations: p A , p B 1 , p B 2 , p C 1 , p C 2 and p C 3 . As expected, all weights that correspond either to genuine three-particle correlation ( p A ) or genuine two-particle correlations ( p B 1 and p B 2 ) diminish rather quickly as multiplicity increases.
In the next section, we use the clear-cut toy Monte Carlo study to confirm the validity of our main result in Eq. (24).

Toy Monte Carlo study for the three-particle case
In this section, we set up a toy Monte Carlo study to confirm the validity of our main result in Eq. (24). We organize the study in such a way that it does not lead to simplifications at any step. For instance, we choose random observables x, y and z to have different functional forms of single-variate marginal p.d.f.'s, defined over different sample spaces.
As defined above, the three stochastic variables are x, y and z, while M is a parameter that corresponds to multiplicity. The normalization constraint is satisfied for any M, i.e.
measured two-and three-particle correlations were derived. The main result is the demonstration that these solutions for the combinatorial background are universal-they can be written generically in terms of multiplicity-dependent combinatorial weights and marginal probability density functions of starting multivariate distribution. Finally, the new general results between expectation values of multiparticle azimuthal correlators and flow amplitudes and symmetry planes were presented.
In the derivation, we have used result 1.320 from Ref. [18] in combination with the orthogonality relations of trigonometric functions. As an example, we can now easily read off the following specific cases: For instance E[cos n 1 ϕ i cos n 2 ϕ j ] = v n 1 v n 2 cos Ψ n 1 cos Ψ n 2 , ∀n 1 , n 2 , i, j (i = j) , (A.12) E[sin n 1 ϕ i sin n 2 ϕ j ] = v n 1 v n 2 sin Ψ n 1 sin Ψ n 2 , ∀n 1 , n 2 , i, j (i = j) , (A. 13) We now demonstrate that the expectation values of all highlevel observables can be reduced to these fundamental ones.

Compound observables
The expectation value of any compound observable can now be trivially deduced from the fundamental expectation values presented in the previous section. For instance, for the single-event average two-particle correlation, 2 , we use the following statistics built from azimuthal angles ϕ 1 , . . . , ϕ M (M is the multiplicity of an event (A.14) Its expectation value is