A Fourier-cumulant analysis for multiharmonic flow fluctuation

The Fourier analysis of the final particle distribution followed by cumulant study of the Fourier coefficient event-by-event fluctuation is one of the main approaches for testing the collective evolution in the heavy-ion collision. Using a multidimensional generating function, we propose a method to extract any possible cumulant of multiharmonic flow fluctuations and classify them in terms of the order of cumulants and harmonics involved in them. In particular, we show that there are 33 distinct cumulants with order 2, 3, 4, 5 and harmonics 2, 3, 4, 5. We compute the normalized version of these cumulants from hydrodynamic simulation for Pb–Pb collisions based on T\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_\mathtt{R}$$\end{document}RENTo+VISH2+1+UrQMD. We compare the simulation with those normalized cumulants that the LHC has measured and predict the unmeasured ones. Comparing the initial and final state fluctuation normalized cumulants, we compute the linear and nonlinear hydrodynamic response couplings. We finally introduce the genuine three-particle correlation function containing information of all third-order cumulants.


Introduction
It is widely accepted in the heavy-ion community that QCD matter in the deconfined phase, quark-gluon plasma, is produced in ultrarelativistic heavy-ion collisions. The collective models equipped with QCD as an underlying physics can successfully explain the majority of experimental observations. However, measuring the accurate transport coefficients such as shear and bulk viscosity over entropy density and finding QCD critical point are still one of the main challenges in the heavy-ion community [1][2][3]. Moreover, observing similar collectivity signals in small system collisions such as p-p or p-Pb [4][5][6][7][8][9][10] raise debates about relation between the observations and the collective processes. These indicate that the ultimate goal has not been achieved, and we need to deepen a e-mail: s.f.taghavi@tum.de (corresponding author) our understanding of heavy-ion collision physics. To this end, introducing new observables to probe various aspects of the full image is crucial. The present manuscript tries to introduce a systematic procedure to generate observables related to the anisotropic flow and classify them. Some of these observables are studied before, and some others are introduced for the first time.
The anisotropic particle emission in the azimuthal direction [11][12][13] is one of the most important evidences of collectivity in ultrarelativistic heavy-ion collisions. The final particle momentum distribution in the azimuthal direction can be expanded in a Fourier series, with Fourier coefficients v n e inψ n , called flow harmonics. The imprints of different stages of a heavy-ion collision, preequilibrium, initial state, collective evolution, and freezeout are cumulatively encoded in the flow harmonics. On account of the quantum mechanical nature of partons inside the nucleons, the initial geometry has a complicated structure and more importantly fluctuates from one event to the other (event-byevent fluctuation) leading to nonvanishing and fluctuating v n and ψ n for any n [4,[14][15][16][17][18][19][20][21][22][23]. Practically, we measure the statistical properties of the v n e inψ n fluctuation in the experiment. Several attempts focus on the moments [24][25][26][27][28][29][30][31][32][33] or on cumulants of fluctuations to achieve information from flow harmonic fluctuations. Single harmonic cumulants [34,35], symmetric cumulants [36][37][38][39], generalized (or higher-order) symmetric cumulants [40,41], asymmetric cumulant [24] are some examples of these attempts. It is noteworthy that each cumulant of a probability density function (p.d.f.) contains independent information about the underlying fluctuation. Therefore, a complete study of all possible cumulants up to a given order of cumulant expansion is crucial to gain a comprehensive insight into the fluctuations. In the present work, we start with the standard defini-tion of cumulants based on generating functions for an arbitrary number of variables. We extract and classify all possible cumulants depending on their order and the involved number of harmonics. Considering the space limitation, we focus on cumulants of harmonics n = 2, 3, 4, 5 up to fifth order, which contains 33 distinct cumulants. Among them, there are already known cumulants together with new cumulants containing symmetry plane correlations. To ease the future application, we have prepared a Mathematica [42] package accessible as an ancillary file for the present manuscript or in the GitHub repository [43]. The package returns the cumulants in three different forms: in terms of v n and ψ n moment, final particle azimuthal angle correlations, and Q-vectors (see Eq. (29)) notations. For cases that the Mathematica software is not available, a list of few first multiharmonic cumulants in terms of v n and ψ n symbolic moments are tabulated in the Appendix A.
Regarding multiharmonic cumulants, other studies have been done previously [44,45] with some similarities and differences compared to the present work. In Ref. [44], the underlying p.d.f. of the fluctuations are classified into three categories, flow-amplitude (only v n fluctuations), eventplane-correlation (only ψ n fluctuations) and mixedcorrelations (both v n and ψ n fluctuations). We find that, instead of classifying the p.d.f. into three types, considering one general p.d.f. is more suitable way to extract all possible cumulants up to a given degree. The author believes that the cumulants extracted in the present paper are achievable with the method introduced earlier in Ref. [45]. In this reference, however, the explicit calculations are limited to the simple few known observables, and the dependence of the cumulant expressions to the harmonics is not presented. Our explicit calculations show that the cumulants' form is related to the involved harmonics in the cumulant (see also Refs. [44,46]). For instance, we see that the first cumulant involving simultaneous harmonics n = 2, 3, 5 appears in the third order, and the next cumulants appear in the fifth order in three different forms. For the case involving simultaneous harmonics n = 2, 4, 5, however, cumulants starts from fifth order with two different forms (see Table 1). Besides, our study is equipped with Monte Carlo hydrodynamic simulation as well (T R ENTo [47]+VISH2+1 [48,49]+UrQMD [50,51] with parameter calibrated by a Global Bayesian analysis [52]) to see which unmeasured cumulants would have larger signals in future measurements. Comparing this simulation with future measurements can also validate the heavyion collision parameters tuned with Global Bayesian analysis. We compare the simulation with normalized cumulants that are measured at the LHC, which reveals a rather good agreement.
We study the initial state fluctuation cumulants to see how the final state inherits the fluctuations from the initial state. The initial energy density anisotropy translation to the final state momentum anisotropy is formulated by a set of equations so-called hydrodynamic response [53]. By comparing the initial and final state fluctuations, we obtain the hydrodynamic response coefficients. Moreover, our systematic cumulant study guide us to extend the two-particle correlation function (2PC) notion [5,[54][55][56] to multiparticle correlation functions, qPC (see the definition in Sect. 6). We specifically study the flow-induced 3PC for harmonics 2 to 5 via our hydrodynamic simulation. This quantity contains information on all third-order cumulants.
This paper is structured as follows: The Sect. 2 is dedicated to introduce generating function method for multiharmonic cumulants. In Sect. 3, the functions available in the Mathematica package are introduced. A realistic Monte Carlo study of the first Fourier-cumulant expansion terms is presented in Sect. 4. The linear and nonlinear hydrodynamic response couplings are extracted by comparing the initial state, and final state fluctuation normalized cumulants in Sect. 5. Finally, we study genuine three-particle correlation functions in Sect. 6. In the Appendix A, a list of few first cumulants is tabulated. In Appendix B, a technical study about the statistical fluctuation of multiharmonic correlations is presented. Detail of some derivations in the text is presented in the Appendices C and D.

Multiharmonic generating function
Despite the simplicity of Eq. (1), experimental measurement of flow harmonics is a challenging task. The number of final produced particles per event in a given collision is not enough to extract statistically accurate flow harmonics values. Many events average resolves the problem of low statistics while leads to the convolution of event-by-event fluctuation of the initial state into the measurement. Strictly speaking, the measurable quantities are the moments or cumulants of the following gigantic probability density function (p.d.f.) p(v n 1 ,x , v n 1 ,y , . . . , v n k ,x , v n k ,y ), where we have used the Cartesian coordinate notation for the flow harmonics, v n,x = v n cos nψ n and v n,y = v n sin nψ n .
Noting the fact that f (ϕ) (Eq. (1)) is assigned to each initial state, the p.d.f. (2) is tightly related to the event-by-event fluctuation of the initial state, which is, accordingly, related to the underlying physics of nucleons/nuclei. Moreover, the map which connects the above p.d.f. to the initial state fluctuation is governed by the collective evolution process. Indeed, there are stochastic processes in the collective evolution that can be convoluted into the above p.d.f. as well. Besides the experiment's statistical uncertainty issues, the reaction plane angle is not experimentally measurable conveniently, which leads to loss of one more degree of freedom in p.d.f. (2). In the polar coordinate the p.d.f. (2) depends on k flow ampli-tudes v n i s and k symmetry plane angles (or event-planes) ψ n i s. After imposing the randomness of the reaction plane angle, from k variable ψ n i , only k − 1 independent variables ψ n i − ψ n 1 remain Here, ψ n 1 is chosen conventionally to compute the angle differences. One could choose any other combination of two symmetry plane angle differences. The final cumulants are independent of this choice.
In the present section, we introduce the method of generating function to extract all possible cumulants of the p.d.f. (3). To this end, we start with a standard definition of generating function of a generic multivariate p.d.f. g(x) (see for instance [57]), where x is a m-dimensional random vector. The cumulants of the distribution function g(x) (shown by K a 1 ,...,a m ) are the Taylor expansion coefficients of the cumulant generating function log G(k), The above generating function returns the cumulants of p.d.f. (2), while in heavy-ion physics, we are interested in cumulants of p.d.f. (3). To clarify the difference between cumulants of p.d.f. (2) compared to p.d.f. (3), let us start with the most simple example. Focusing on one-harmonics v n e inψ n , the p.d.f. of the flow fluctuation in the Cartesian coordinate has the form p(v n,x , v n,y ), and its cumulants can be obtained via the definition presented in Eq. (5) for k = 2. For instance, are the width of the distribution in the v n,y and the skewness in the v n,x direction, respectively. In the experiment, however, we lose one of the degrees of freedom in p(v n,x , v n,y ) due to the randomness of the reaction plane angle. After imposing event-by-event random reaction plane angle to p(v n,x , v n,y ), we obtain a rotationally symmetric p.d.f.p(v n,x , v n,y ). The information of the latter p.d.f. is encoded in a radial distribution p(v n ). Traditionally, we assume the cumulants of p(v n ) are c n {2k} [35] where the two first orders of them are given by These quantities, however, are specific combinations of twodimensional cumulants of p(v n,x , v n,y ) that survived after imposing the randomness of the reaction plane angle. For instance, it turns out [58,59], Here, we have assumed thatv n = v n,x = 0 and v n,y = 0 for simplicity. This condition is satisfied when the flow harmonics are sourced purely from fluctuations, namely ellipticity in central Pb-Pb collisions, or triangularity in round nuclei collisions. For nonvanishingv n , similar relations with more terms in the right-hand side of Eq. (8) are obtained. In any case, the cumulants c n {2k} are the maximum information we can achieve from the cumulants of the original distribution p(v n,x , v n,y ). Before proceeding, let us briefly review the generating functional method to find c n {2k} [60]. The generating function in Eq. (4) corresponds to the rotationally symmetric p(v n,x , v n,y ) reads as G(k x , k y ) = dv n,x dv n,yp (v n,x , v n,y )e iv n,x k x +iv n,y k y where in the above J 0 (k v n ) is the Bessel function of the first kind and the averaging J 0 (k v n ) p is performed with respect to p(v n ) ≡ v np (v n ). One can expand the logarithm of the generating function in terms of k, with to find the cumulants c n {2k}. We extend the above procedure into the cases with an arbitrary number of harmonics. To simplify the notation, we rewrite the p.d.f. and Considering the rotational symmetry imposed by the randomness of the reaction plane angle, one can first average the characteristic function over the azimuthal angle to find a symmetric characteristic function. Then the nonvanishing combination of cumulants of p(v n 1 ,x , v n 1 ,y , . . . , v n k ,x , v n k ,y ) can be obtained directly from the symmetric characteristic function, where and For the derivation refer to the Appendix D. Having found G(k, δφ), the cumulants of flow fluctuations can be obtained from the following expansion: where the coefficients are the cumulants (δψ n i are defined in Eq. (13)). In the righthand side of the above equation, we have represented the cumulant with its highest rank moment together with a subscript c. Here, n i = 1, 2, . . . stand for the involving harmonics, m i = 0, 1, 2, . . . are the power of flow amplitudes and α i = 0, ±1, ±2, . . . are the coefficients of the symmetry plane angle differences. The order of the cumulant is given by By definition, n i s are all distinguished while, without loss of generality, we impose a strictly ascending order convention to them, n 1 < · · · < n k . 1 The coefficient (m, α) is a numerical factor that does not depend on the moments. We fix the coefficient so that the numerical factor of the highest rank moment in cumulant c . . , m k } turns to be equal to unity. As we will see in the next section, some combinations of harmonics n i , m i , and α i lead to vanishing cumulant. A trivial example is k = 1 where the Eq. (18) reduces to c n {m} ≡ v m n c . It is known that this cumulant is nonvanishing only for even m (see Eq. (7)).
Each distinguished cumulant c . . , m k } contains a piece of independent information about the p.d.f. (3). To numerically compute these cumulants for a 1 In this manuscript, we use harmonic index a i for those repeated in the sequence, while n i is reserved for strictly ascending indices.
given p.d.f., we need to find the explicit form of them written in terms of the moments. Although the cumulants c {α 1 ,...α k−1 } n 1 ,...,n k {m 1 , . . . , m k } can be found analytically using Eq. (17), its computation is cumbersome, except for the case k = 1. By considering one flow harmonic (k = 1), the generating function reduces to that mentioned in Eq. (9). Keeping two flow harmonics, the flow fluctuation distribution contains three degrees of freedom v n 1 , v n 2 , and δψ n 1 . In this case, the function J can be written in terms of generalized Bessel function [61], and by expansion, one obtains the cumulants c However, extracting them is arithmetically more involved. For more general cases, the complexity increases, forcing us to choose a more practical and efficient way for the computation.

One package for all cumulants
Using the cumulant generating function in Eq. (17), one can find the cumulants c . . , m k } written in terms of symbolic moments of variables v n and ψ n . These results can be immediately used in theoretical studies where v n and ψ n of every single event are accurately accessible. In the experiment as well as some simulations, however, the outcomes are particles azimuthal angles. As a result, one needs to invest an extra effort to translate the azimuthal angle of final particles into the averages containing v n and ψ n . For that, we use the multiparticle techniques to compute the flow harmonic cumulants. We rewrite the moments in terms of Q-vectors (see Eq. (29)) which can be calculated from final particle azimuthal angles at every single event and employed in computing the cumulant [36,62]. This section presents a practical way to find cumulants and their statistical uncertainties in terms of symbolic moments, correlation of particle azimuthal angles, and Q-vectors.

Cumulants from generating function
As it is mentioned in the previous section, the analytical computation of a generic cumulant c α 1 ,...α k−1 n 1 ,...,n k {m 1 , . . . , m k } is cumbersome. For that reason, we do it symbolically in Mathematica. We encapsulated this computation into different functions (see below) and implemented them into a Mathematica package. One can load the package in a separate Mathematica notebook and recall the functions. This package is available as an ancillary file of this manuscript, or in the GitHub repository [43].
In the package, we first compute the Taylor expansion of the exponential function in Eq. (15) in terms k n 1 , k n 2 , . . . up to order m 1 , m 2 , . . .. Then we perform the integral in Eq. (15). After that, we replace the combinations v w 1 n 1 · · · v w k n k e β 1 δψ n 1 · · · e β k δψ n k in the expansion with a sym-bolic variable as a moment. Computing the logarithm of the result, we read the coefficients of the Taylor series for the variables k 1 , k 2 , . . . , and the Fourier series for the variables α 1 , α 2 , . . . (see Eq. (17)). We extract the cumulants in terms of our symbolic moment variables up to a numerical factor (m, α) by comparing the result with Eq. (17). Finally, we single out the highest rank moment in the result and fix the coefficient (m, α) such that the numerical factor of this moment turns to unity. We listed the available functions in a short manual at the header of the package file. One of these functions is that returns the associated cumulants written in terms of moments of variables v n and ψ n . For instance, In [3] :=c[2, 2, 0, 2, 3,v,ψ] By using cMean instead of c, the angle brackets are replaced by The above representation of cumulant is useful when flow harmonics are accurately accessible in a single event, typically in hydrodynamic simulations. If the final state is the particle azimuthal angles, one needs to employ the particle correlations, k a 1 ,...,a k ≡ e i a 1 ϕ i 1 +···+i a k ϕ i k , Here, ϕ i is the azimuthal angle of the ith particle in an event, and M is the multiplicity of the event [34][35][36]62]. A generic flow harmonic moment V a 1 ,...,a k with can be written in terms of particle correlations as V a 1 ,...,a k = k a 1 ,...,a k . The function where, From the practical point of view, computing k a 1 ,...,a k contains several nested loops depending on the value of k which are computationally expensive. A technique is introduced in Refs. [36,62] where the moments k a 1 ,...,a k are obtained in terms of Q-vectors, 2 Using this technique, only one loop over particles in an event is needed to compute the correlations.
For those who want to analyze anisotropic flow inside Mathematica, we implemented the recursive algorithm mentioned in Ref. [36] into the package to find k a 1 ,...,a k in terms of Q-vectors. Substituting v n and ψ n symbolic averages with moments written in terms of Q-vectors, we find cumulants in terms Q-vectors (so-called Q-cumulants), where M and Q are the multiplicity and Q-vector symbols, respectively. For example, which is nothing but generalized (or higher-order) symmetric cumulants for some combination of harmonics [40] (the actual form of the cumulant in terms of moments can be found in the Appendix A), 3 k,l,m {2, 2, 2}, and asymmetric cumulant [21,22,24], as special cases of what we present here. There are other cumulants, however, missed in other studies so far. For instance, is an example of a cumulant which have not been studied before.

Statistical uncertainty of cumulants
Having computed the flow harmonic cumulants, we discuss the statistical fluctuation of the cumulants. In case the flow harmonics v n e inψ n are accurately accessible in a single event, one can follow a standard procedure (see for instance Ref. [63]) to find the covariance matrix of the moments 3 The equivalence of generalized symmetric cumulants and k,l,m {2, 2, 2} is not exact. As a different approach in Ref. [40], the cumulants of flow amplitude squared fluctuations are studied, which for most of the cases it is equivalent to what is presented here. V a 1 ,...,a k (see Eq. (23)). The covariance matrix of two generic moments are simply given by where N is the number of events and the subscript i stands for a generic collective index a 1 , . . . , a k . The variance of any function of moments, f ( V 1 , V 2 , . . .) is obtained by [63], For instance, the variance of If the azimuthal angle of particles in the final state are available, the flow harmonic V a 1 ,...,a k should be replaced by particle correlations Re k a 1 ,...,a k . Although Re k a 1 ,...,a k is not an accurate estimation for V a 1 ,...,a k , in the ultimate many events average it approaches to an accurate estimation. (38), we obtain the variance of the function of correlations in terms of correlations. For instance, in this notation, we have c 2 {4} = 4 −2,−2,2,2 − 2 2 −2,2 , and accordingly, the example in Eq. (39) turns to the following form: where we have used the explicit form of covariance matrix in Eq. (38). One can immediately substitute the quantity k a 1 ,...,a k in terms of Q-vectors for practical computations by using Refs. [36,62]. The procedure explained above is implemented into the following functions in the package: The output of the above functions is N σ 2 f . The input function (func) can be any function of correlations but the correlations should be always written in the form corr[a 1 , . . . , a k ] as the output of Eq. (24). 4 For instance, to find N σ 2 c 2 {2} by using Nsigma2, we should call the function as follows: The difference between functions in Eq. (41) is in their output presentations. The output of function Nsigma2 is shown in Eq. (42), the function Nsigma2Mean output is similar to cMean function, and Nsigma2Qvec returns the variance in terms of Q-vectors similar to cQvec function (see Eq. (31)).
As an example, the statistical error of c 2 {2} = 2 −2,2 in terms of Q-vectors can be found below: Referring to Eqs. (40) and (42) (and to any other Nsigma2 outcomes), we find that the statistical uncertainty has no explicit multiplicity dependence. Let us consider two equal-size sets of events when the multiplicity of events in one set is smaller than in the other set. We expect that the statistical uncertainty of cumulants computed from the set of events with smaller multiplicity is larger than the statistical uncertainty obtained from the other set with higher multiplicity events. This apparent contradiction with intuition is due to the presence of correlations k a 1 ,...,a k a 1 ,...,a in the statistical uncertainty (for instance 2 −2,2 4 −2,−2,2,2 in Eq. (40)). These forms of correlations are the consequence of cov(V i , V j ) in Eq. (38) and have not appeared in previous studies in Refs. [36,62]. In fact, the correlations k a 1 ,...,a k a 1 ,...,a depend on multiplicity due to the presence of autocorrelations remained in them. To be more specific, we focus on the first terms in the output in Eqs. (42), The above relation is a four-particle correlation with some remaining autocorrelations. In the technical Appendix B, a more direct approach is employed to extract the statistical uncertainties. This approach is implemented in the function As seen from the above equation, on the one hand, there is no correlation with the form k a 1 ,...,a k a 1 ,...,a , and on the other hand, there are explicit M dependences as we expect. One can explicitly show that the terms inside the first bracket in Eq. (45) are exactly equal to 2 2 −2,2 . It becomes more apparent when we replace all k a 1 ,...,a in Eq. (45) with Q-vectors using the results of Refs. [36,62] where one finds that Eq. (45) in terms of Q-vectors is identical with that mentioned in Eq. (43). In general, the result of function (41b) is identical with that obtained from function Nsigma2P.
As a final remark in this section, we would like to compare our statistical uncertainty results with previous studies. By ignoring the event-by-event fluctuation, all events would be identical. Therefore, we can ignore the outer angle brackets in Eq. (45). Also, by collecting many events we can increase the accuracy such that each k a 1 ,...,a k estimates the true value of the flow harmonic V a 1 ,...,a k (see Eq. (23)). In such an ideal scenario, Eq. (45) turns into the following form: which has been obtained before in Refs. [36,64] for the case with ψ 2 = ψ 4 .

Fourier-cumulant expansion to study flow distribution
An important property of cumulant analysis is that the lowest order cumulants capture the global features of the p.d.f., and by moving toward higher-orders, more detailed features become relevant. The same is true for Fourier expansion of f (ϕ) in Eq. (1) where lower harmonics contain more coarse-grained pictures than higher harmonics. As a result, we consider a double expansion for studying p.d.f. (3): first, harmonic expansion, and second cumulant order expansion. According to this approach, c 2 {2} contains the most coarsegrained information (we ignore direct flow in the present study). The harmonic n = 2 is the first dominant flow, and c 2 {2} is the only possible second-order cumulant for n = 2. It is well-known that c 2 {2} contains information about the initial state ellipticity together with elliptic flow fluctuation. Table 1 List of all cumulants with order q = 2, 3, 4, 5 and harmonics n = 2, 3, 4, 5

Cumulant
Order Cumulant expression The first third-order cumulant is c {4} 2,4 {2, 1} which means, in addition to second harmonic n = 2, the harmonic n = 4 should be involved as well. Keeping more harmonics and higher-order cumulants, more detailed information about the p.d.f. (3) and f (ϕ) can be achieved. The number of cumulants increases rapidly by including more harmonics and keeping higher-order cumulants.
Henceforth, we study harmonics n = 2, 3, 4, 5 and cumulants up to fifth order, which includes 33 different cumulants as they are tabulated in Table 1. Apart from those cumulants that contain only one moment, there are cumulants containing more than one event-plane moment (see lines 21,24,25,26,27,29, and 32 in Table 1) which are studied here for the first time. We should point out that three different five-particle cumulants have been measured by ALICE [65,66] two of them (see [66]) contain more than one event-plane moment. The direct flow (flow harmonics with n = 1) is involved in all these cumulants. The cumulant c {4} 2,4 {4, 1} (see Ref. [44] as well) in line 22 is equivalent with v 4 {5} measured by CMS up to a normalization factor [67].

Normalized cumulants from initial and final states
To understand the observed cumulants' origin, let us recall that the central part of the flow fluctuations comes from the initial state. To study the initial state fluctuation, we need to quantify the shape of the initial state with quantities similar to those done for flow harmonics, namely eccentricities [68], where {· · · } = rdrdϕε(r, ϕ) is the averaging with respect to the energy density ε(r, ϕ). The true initial energy density shape is captured by the cumulants of the distribution ε(r, ϕ) not its moments (eccentricities) in Eq. (47) [53]. The cumulants of the energy density have been obtained by Teaney and Yan in the seminal paper [53] and follow-up papers [25,69]. Employing the convention in Ref. [69], the two-dimensional energy density cumulants in the polar coordinate are given by One notes that the eccentricities for n = 2 and 3 are cumulants as well, while for n > 3 there are contributions from lower harmonics eccentricities and radial shape of the energy density (see Eqs. (2.9) and (2.11) in Ref. [69]). The anisotropic flow in the final state is a collective response to the initial anisotropic geometry. This process is formulated via the following response relation [53], v n e i n ψ n = w n C n e i n n + nonlinear terms. (49) Up to the leading order, the physics of collectivity and freezeout is encoded in the linear coupling w n . Ignoring the nonlinear terms, the flow harmonic cumulants of the p.d.f (3) should be proportional to the cumulants of the p.d.f. p ini (C n 1 , . . . , n 2 − n 1 , . . .). The cumulants of the latter p.d.f. are obtained precisely similar to those studied so far by replacing v n e inψ n with C n e in n . 5 To be more precise, let us recall the important properties of cumulants, homogeneity. Based on this property, by scaling the random variable x i as x i → w i x i , the cumulants defined in Eq. (5) are scaled as One can examine the above relation in two explicit examples shown in Eq. (6). By ignoring the nonlinear parts in Eq. (49), we can relate the cumulants c . . , m k } obtained from eccentricities and flow harmonics as The above relation indicates that the initial state and final state fluctuation cumulants differ by a numerical factor, similar to that mentioned in Eq. (50). For comparing the initial and final state cumulants in the linear response approximation, the values of w n are needed. However, it is still possible to modify our cumulant definition such that it is independent of w n numerical value. To this end, we define normalized (standardized) cumulants, Referring to Eqs. (51) and (52), one finds that the normalized cumulants of the initial anisotropy and flow harmonics fluctuations must be precisely the same at the linear approximation. Any deviation between two normalized cumulants should be sourced from nonlinear terms. The definition of normalized cumulants in Eq. (52) is compatible with the suggestion in Refs. [25,70] made for event-plane correlations. According to the above definition, for k = 1, we have compatible with Ref. [71]. Also for nc where NSC(m, n) is the normalized symmetric cumulant. There is an alternative way of defining the normalized cumulant, footnote5 continued text. We also study C n e in n event-by-event fluctuation cumulants (the cumulants of the p.d.f. p ini (C n1 , . . . , n2 − n1 , . . .)) which are called initial state fluctuation cumulants in the manuscript.
This definition is compatible with the scalar product method mentioned in Ref. [30]. For cumulants with all m i = 1, two approaches coincide. In the present study, we use the first convention in Eq. (52) when all cumulants are normalized with the first single-harmonics cumulants. In comparison between the LHC data and simulation, we employ the alternative approach only when the first approach is not measured.

Normalized cumulants from simulation and the LHC
We study the first 29 normalized cumulants in the Fourier-cumulant expansion by using a hybrid hydrodynamic model, T R ENTo [47]+VISH2+1 [48,49]+UrQMD [50,51] (nc n {2} = 1 by definition). We simulate events for Pb-Pb collisions at √ s NN = 2.76 TeV. The global Bayesian analysis calibrates the model's free parameters (consist of the temperature dependence of shear and bulk viscosity over entropy density) to explain the measurements by ALICE experiment [52]. In Fig. 1, the normalized version of flow harmonic cumulants listed in Table 1 are plotted with red filled circles. We kept a fixed range for the vertical axes in all panels to simplify the magnitude comparison of different normalized cumulants. In the calibration, the cumulants c n {2} for n = 2, 3, 4, and c 2 {4} are used. This means that only the information of normalized cumulant nc 2 {4} (in panel (7)) has been used for calibration and the rest of the simulations are predictions.
Although the simulation is calibrated with ALICE experiment kinematics (0.2 < p T < 5 GeV and |η| < 0.8), the normalized versions of cumulants are less sensitive to these kinematics. In particular, the transverse momentum range dependence of the normalized cumulants has been examined in Refs. [22,38]. As a result, in addition to the 6 The normalization conventions in Eqs. (55) and (52) Fig. 1 and not displayed in Fig. 2 contain an independent piece of information that can be examined in the future experiment. Here, we only focus on the Pb-Pb collisions at fixed center-of-mass energy. The cumulants' system size and energy dependence would lead to interesting information about the underlying initial state fluctuation and collective evolution in the future. Referring to Fig. 1, we find that the normalized symmetric cumulants nc  (19)), and fifthorder normalized cumulants nc  (30)) are the largest cumulants up to n = 5 harmonics that have not been measured yet.
Concerning the feasibility of the cumulant measurements, there is a resent analysis by ALICE collaboration measuring higher-order normalized symmetric cumulants (shown by NSC(m, n, )) consisting flow amplitudes v 2 , v 3 , v 4 and v 2 , v 3 , v 5 [41]. To calculate these sixth-order cumulants, the moments as v 2 2 v 2 3 v 2 4 and v 2 2 v 2 3 v 2 5 are needed to be measured. Given that these cumulants are one order higher than those we presented in Table 1, we would expect that all the unmeasured normalized cumulants in Table 1 are experimentally accessible in the near future. The cumulant nc −3,5 2,3,5 {1, 1, 3} (see line25 in Table 1) contains moments with v 3 v 3 5 combination. Compared to NSC(2, 3, 5), it has one power higher for flow amplitude for v 5 and one power lower for v 3 . Considering its large signal from the simulation, we expect that an experimental observation for this cumulant as well.
So far, we have studied the normalized cumulants of Pb-Pb collisions and investigate to what extend the Monte Carlo simulation can explain the experimental data. To better 7 The measured quantity v 4 {5} in Ref. [67] is equivalent with c   Table 1 are presented  Fig. 2 The (alternative) normalized cumulants that measured by ALICE [37,38] and ATLAS [30] compared with T R ENTo+VISH2+1+UrQMD understand both initial state and collective evolution, we have also compared the initial and final state fluctuation cumulants. In the next section, we investigate it in more detail. As we have discussed earlier, up to the linear order of hydrodynamic response, the normalized cumulants of flow harmonics and energy density cumulant fluctuations must be the same. As a result, the observed difference in Fig. 1 is due to the nonlinear hydrodynamic response.

Linear and nonlinear hydrodynamic response coefficients
In Sec. 4.1, it has been discussed that by comparing the normalized cumulants of the initial and final state fluctuation, we can obtain information about the collective hydrodynamic evolution. In the present section, we quantitatively study the hydrodynamic response coefficients.
In the present study, we keep only the linear terms for n = 2, 3, and the first subleasing nonlinear terms for n = 4, 5. In particular, we consider the following explicit form of Eq. (49), v 5 e i5ψ 4 w 5 C 5 e i5 5 + w 5(23) C 2 C 3 e i2 2 +i3 3 .
We will discuss the interpretation of these two conventions later. Different approaches can be employed to extract the couplings w n , w 4(22) , and w 5(23) (or k n , k 4(22) , and k 5(23) ). We can start from a Gaussian geometrical initial energy density and deform it with one (or a few numbers) of nonvanishing C n e in n . Then we change the value of the energy density cumulant to probe the hydrodynamic response from this "single shot" simulation [25,69]. It is shown that (Marcinkiewicz theorem 8 ) reproducing a distribution (except Gaussian dis- 8 Marcinkiewicz has proved that only Gaussian distribution has a polynomial generating function with a finite number of nonzero cumulants [73]. tribution) with a finite number of cumulants leads to negative values at some parts of distribution domain. These negative values should be regulated, which consequently produces spurious cumulants [69]. Until spurious cumulants are small, we can neglect their effect in the coupling estimations. The initial state with larger deformation would have more such a problem.
As another approach, one can generate many events with complicated initial geometries and study the hydrodynamic response. For instance, to extract the response coefficient k n in Refs. [68,74], the authors multiply both sides of Eq. (49) (ignore nonlinear terms) with n e −inφ n and then average over many simulated events in a given centrality class. This approach leads to the following estimation: k est n = n v n cos(n(ψ n − φ n )) / 2 n . By writing this estimator in the Cartesian coordinate, n e inφ n = n,x + i n,y and v n e inψ n = v n,x + iv n,y , we see that this relation is the Pearson correlation between eccentricities and flow harmonics. This estimator is related to the covariance matrix of p |v ( n,x , n,y , v n,x , v n,y ). To estimate the nonlinear hydrodynamic couplings similar to what is done in Ref. [68], we need to include more harmonics in p.d.f. p |v . Using this approach, we have ignored the non-Gaussian effects encoded in p |v .
Comparing the event-by-event fluctuation of the initial and final states in a simulation is an alternative method (the one we employ here) to study the hydrodynamic response couplings. For that, we technically compare p(v n 1 , . . . , ψ n 2 − ψ n 1 , . . .) with p ini (C n 1 , . . . , n 2 − n 1 , . . .). By substituting v n e inψ n in the flow harmonic fluctuation cumulant (moments) with the corresponding expression for the hydrodynamic response Eqs. (56), we obtain a function written in terms of initial state fluctuation cumulant (moments). For instance, by employing Eq. (56c), we obtain Similarly, we work out the following flow harmonic cumulants in terms of C n e in n , We can rewrite the right-hand sides of the above equations in terms of cumulants as well.
There are six equations (Eqs. (58) and (59)) and six unknown response coupling constants that we can find numerically. We do a naïve analysis by ignoring the statistical errors in our Monte Carlo simulation to extract the unknown coefficients. At the linear level, the initial state's pressure gradient enforces that the event-plane angle ψ n to be the same as the participant plane n . As a result, we expect that all w n s are real and positive. By demanding that the linear response coupling constants w n , n = 2, . . . , 5 are positive, we obtain only one set of solutions for couplings at each centrality class. The results are depicted in Fig. 3. The couplings k 2 and k 3 (which are identical with w 2 and w 3 in our study) are explicitly reported in Ref. [74]. Considering that a different hydrodynamic model (with different tuning) has been used, and a different method (mentioned in Ref. [68]) is employed to extract the couplings, our results for w 2 and w 3 are compatible with those computed in Ref. [74]. The couplings w 4 and w 5 are computed in Ref. [25]. The values of these couplings approach zero and change the sign at mid-central collisions. In our method, there are valid positive solutions for w 4 and w 5 with no sign change. We have also found larger values for these couplings. However, our computations have rather similar behavior as those reported in Ref. [25] for the nonlinear coefficients. The observed differences between the two studies could be due to the different hydrodynamic models or different approaches of extracting the coefficients (singleshot approach has been employed in Ref. [25]). One notes that our estimations depend on the number of nonlinear response terms that we have considered in our hydrodynamic response estimation. By adding more terms, we need to employ more moments (cumulants) as input. This modification can lead to slightly different values for the couplings.
As it has been mentioned in Sect. 4.1, to understand that how much flow harmonic fluctuation is originated from the initial state, one can compute the initial state fluctuation normalized cumulants and compare them with those obtained from flow harmonic. In Fig. 1, C n e in n fluctuation normalized cumulants (blue empty circles) and n e inφ n fluctuation normalized cumulants (magenta squares) are depicted. The observed difference between n e inφ n and v n e inψ n cumulants is interpreted as the presence of nonlinear terms similar to what is mentioned in Eqs. (57). The comparison of C n e in n and v n e inψ n gives us a hint about the nonlinear terms similar to what has been shown in Eqs. (56).
We see in the figure that some of the normalized cumulants that are computed from n e inφ n have a different sign from those calculated from v n e inψ n (see for instance panels (5), (6), (22), and (30), in Fig. 1). In fact, the "wrong" sign of n e inφ n fluctuation in participant plane correlation has been observed in the previous studied [25,30,69,[75][76][77]. In particular, the quantities χ 422 and χ 523 in [77] are similar to nc 2,3,5 {1, 1, 1}, up to a normalization factor. A sign difference between initial and final state fluctuations has been observed for these quantities. The authors of Ref. [77] conclude that this sign difference is a signature of hydrodynamic response to the initial state. Up to a normalization factor, the same quantities have been studied in [25]. Again a sign change has been observed between eccentricity and flow harmonic fluctuations while the sign change is resolved by replacing n e inφ n with C n e in n . However, we cannot conclude that replacing n e inφ n with C n e in n always leads to a compatible correlation sign with final state fluctuation. For instance, one can find cases in panels (18), (23), and (28) of Fig. 1 that C n e in n correlations have opposite sign compared to flow harmonic fluctuations. Few examples can be found in [25] as well. 9 The cumulants calculated from hydrodynamic nonlinear response estimation are shown with black filled triangles in Fig. 1, where the response couplings are those displayed in Fig. 3. The figure shows that the two first cumulants (c  Fig. 1, only two ratios w 4(22) /w 4 and w 5(23) /w 5 play a role. As a result, two normalized cumulants are enough to fix these ratios. The rest of the normalized cumulants computed from the nonlinear response are, in fact, Fig. 4 The T R ENTo initial state radial shape (diamond points), and the nonlinear over linear hydrodynamic response couplings from VISH2+1+UrQMD (circle and square points) the predictions. Some of these predictions perfectly match with hydrodynamic computations. For instance, nc  (15), (22), and (30) in Fig. 1. There are cases with a poor agreement between hydrodynamic simulation and the nonlinear response estimation, namely nc 4 {4} and nc 5 {4} (panels (9) and (10) in Fig. 1). The nonlinear terms could not cure the initial and final correlation sign differences in panels (18), (23), and (28). Including more nonlinear terms would increase the accuracy of the latter cases.
Concerning the interpretation of the nonlinear hydrodynamic response coefficients, we note that the true deformation of the initial energy density is quantified by cumulants C n e in n . One expects that the final anisotropy to be proportional to the true deformation [53]. Assuming w 4(22) and w 5(23) have solely collective evolution contributions, the nonlinear couplings k 4(22) and k 5(23) in Eq. (57) have contributions from both initial shape and collective evolution. One can relate these two couplings by substituting cumulants (48) into Eqs. (56) and compare them with Eqs. (57). The linear couplings are identical in both conventions, w n = k n . However, we obtain the following relation for the nonlinear couplings: The above relations mean that the coefficients k 4(22) and k 5(23) receive contributions from the radial shape of the initial energy density as well. In Fig. 4, the coupling ratios w 4(22) /w 4 and w 5(23) /w 5 are plotted with square points (blue curves). The couplings k 4(22) /k 4 and k 5(23) /k 5 (shown by circles, black curves) are estimated by solving six equations in Eqs. (58) and (59) where energy density cumulants are replaced by eccentricities. Since the coupling ratio estimations are estimated from many events in a given centrality, we compute the right-hand side of Eqs. (60), where the initial energy density contributions are averaged over events in the given centrality class. The quantities 3 {r 2 } 2 /{r 4 } and 10 {r 2 }{r 3 }/{r 5 } (shown by diamonds, red curve) are directly computed from T R ENTo events. The estimated k 4(22) /k 4 and k 5(23) /k 5 from Eqs. (60) are shown by triangles (green curves) in Fig. 4. We see a rather perfect match between black curves and green curves, as we expected. From the figure, we see a nontrivial centrality dependence of quantities 3 {r 2 } 2 /{r 4 } and 10 {r 2 }{r 3 }/{r 5 } . This observation indicates that there are contributions from the radial shape of the initial energy density in the values of k 4(22) and k 5(23) . This contribution should be taken into account in interpreting the nonlinear couplings k 4(22) and k 5(23) as hydrodynamic response couplings.

Flow-induced genuine three-particle correlation
The observation of long-range correlations between particles at Δϕ 0 and nonzero Δη is one of the first and most important fluid(-like) signals in large (small) system collisions [5,[54][55][56]78]. In particular, function C(Δη, Δϕ), the two-particle correlation function (2PC), 10 quantifies the correlation between two particles in the final state. In the present section, we introduce a generalized version of this function and study its relation with flow harmonic cumulants. Ignoring the experimental complications, in principle, the function C(Δη, Δφ) is measured as follows: we choose all distinguished pairs of particles in an event and compute Δη and Δϕ for each pair and fill a histogram from a collection of pairs in many events. To find the connection between 2PC with cumulants, we focus only on the Δϕ part of the correlation. In other words, we choose all the distinguished particles and compute Δϕ irrespective of their position in the η direction. One can use η-gaps (see Refs. [24,79]) to decrease the contamination of nonflow effects in the correlation function estimation. For two-particle correlation, we find the well-known relation between 2PC and second-order cumulants c n {2}, Here, we have used the notation C 2 (Δϕ) instead of commonly used notation C(Δϕ) for future generalization. The advantage of measuring C 2 (Δϕ) compared to c n {2} is that it contains a cumulative information of all second order cumulants c n {2} with n > 0. We can extend the notion of two-particle correlation function into q-particle correlation function (qPC). Specifically, in the following, we focus on correlations of three particles in the final state, C 3 (Δϕ 1 , Δϕ 2 ) because it is a function of two variables and easy to visualize. For the same reason (and finding a clear connection to the flow harmonic cumulants), we ignore the η dependence similar to Eq. (61). The η-gap method can also be used in this case to decrease the nonflow effects as well.
Employing the systematic study of multiharmonic cumulants presented in this manuscript, we can find an expansion of any qPC in terms q-order cumulants, similar to what has been written in Eq. (61). The technical details of finding the relation between qPC and q order cumulants can be found in Appendix C. Here, we show the final result for 3PC where harmonics n = 2, . . . , 5 are involved,  Fig. 5. Here, we have subtracted 1/8π 3 to focus only on the nontrivial correlation.
As seen in Fig. 5, the correlation reveals repeating patterns which is a consequence of symmetries. Before explaining these symmetries in three-particle correlations, let us discuss them in a more simple case, two-particle correlation functions. To measure two-particle correlation functions, two (charged) particles in a given event are chosen and compute Δϕ = ϕ 2 − ϕ 1 . The signal distribution is obtained by mea- suring Δϕ for many events. 11 Let us call two distinguished particles in a given event as α and β with azimuthal angles ϕ α and ϕ β . Choosing all pairs of distinguished particles in an event and call them as particle 1 and 2, the particle α is labeled as particle 1 and β as particle 2 once and the particle α as particle 2 and β as particle 1 again. The former labeling leads to Δϕ = ϕ 2 − ϕ 1 = ϕ β − ϕ α and the latter to As a result, in the signal distribution both Δϕ and −Δϕ have contribution from a single configuration of particles α and β. This labeling brings us to the following symmetry, C 2 (Δϕ) = C 2 (−Δϕ). One notes that the mentioned "double counting" is also considered in c n {2} = 2 −n,n (see Eq. (22)). As a result, this symmetry is manifestly true in Eq. (61). We also have an obvious periodic symmetry C 2 (Δϕ) = C 2 (Δϕ + 2nπ) for any integer n. Using these two symmetries, we find that we have independent information in C 2 (Δϕ) in the range 0 < Δϕ < π. The correlation C 2 (Δϕ) (or C 2 (Δϕ, Δη)) is mostly reported in the range − π 2 Δϕ 3π 2 to see a clear ridge and shoulder structures (see for instance [78]). Due to the symmetries explained above, both ridge and shoulder structures are symmetric with respect to the axes pass through their peaks at Δϕ = 0 and Δϕ = π and only half of these structures contain independent information in the range − π 2 < Δϕ < 3π 2 . 11 To suppress the combinatorial backgrounds and acceptance effects, a background distribution is constructed from Δϕ = ϕ 2 − ϕ 1 when one particle is chosen from the given even and the second particle from several randomly selected other events. The ratio of signal and background distribution is proportional to C 2 (Δϕ) [80].
When we have three particles in the final state, there are 3! possible ways to label three distinguished particles, leading to six different values for (Δϕ 1 , Δϕ 2 ) with a single configuration of particles. Let us call three particles as α, β, and γ with azimuthal angles ϕ α , ϕ β , and ϕ γ . If particles α, β, and γ are labeled by 1, 2, and 3 (see Eq. (62)), we obtain Δϕ 1 = ϕ β − ϕ α , and Δϕ 2 = ϕ γ − ϕ α . The particles α, β and γ can be labeled as 2, 1, and 3. With this labeling, we have Δϕ 1 = ϕ α −ϕ β = −Δϕ 1 and Δϕ 2 = ϕ γ −ϕ β = Δϕ 2 −Δϕ 1 . Both values for (Δϕ 1 , Δϕ 2 ) and (Δϕ 1 , Δϕ 2 ) refer to the same configuration of particles and exist in the signal distribution. As a result, the final distributions have the following symmetry: Another case is labeling α, β, and γ particles with 1, 3, and 2. This labeling leads to the values Δϕ 1 = ϕ γ − ϕ α = Δϕ 2 and Δϕ 2 = ϕ β − ϕ α = Δϕ 1 . Consequently, the distribution has symmetry C 3 (Δϕ 1 , Δϕ 2 ) = C 3 (Δϕ 2 , Δϕ 1 ). There are three other permutations, but all of them lead to a combination of symmetries explained above. More than the above symmetries, for any configuration of particles, we expect the same probability for the mirrored configurations, ϕ i → −ϕ i , which leads to the symmetry C 3 (Δϕ 1 , Δϕ 2 ) = C 3 (−Δϕ 1 , −Δϕ 2 ). Finally, we have rotational symmetry ϕ i → ϕ i + 2n i π for any integer n i , which eventually leads to the periodic condition for C 3 (Δϕ 1 , Δϕ 2 ) with periodicity 2π . We summarize all the symmetries of 3PC as follows, The symmetries in Eq. (64) are responsible of repeating patterns observed in Fig. 5. Each dashed line in the figure refers to one (or combination of more than one) symmetry(ies). As seen in the figure, the range 0 < Δϕ 1 , Δϕ 2 < 2π is divided into twelve repeating triangular regions. In the experiment and the simulations, we lose the statistics twelve times by considering the full region without gaining any more information. For that reason, we confine ourselves into a "unit cell" of Δϕ 1 and Δϕ 2 containing all the nontrivial non-repeating patterns. Here, we conventionally choose the following region, as it is shown by a black triangle in Fig. 5. If any choice of particles leads to (Δϕ 1 , Δϕ 2 ) outside this triangle, the symmetries mentioned in Eq. (64) can be used to map that point inside it. The edges of the black triangle in Fig. 5 refer to the combinations of the angles, Δϕ 1 and Δϕ 2 with unequal Fig. 6 Dalitz-like plot of the "unit-cell" (black triangle in Fig. 5) of genuine three-particle correlation function lengths. It will have more symmetric visualization if we use the following variables, where 0 < δ 1 , δ 2 , δ 2 < 2π and δ 1 + δ 2 + δ 3 = 2π . Using this variables, we can plot the unit cell of C 3 (Δϕ 1 , Δϕ 2 ) in a Dalitz-like plot as it is shown in Fig. 6. Similar to C 2 (Δϕ) which has cumulative information from all second-order cumulants, experimental measurement of C 3 (Δϕ 1 , Δϕ 2 ) (or equivalently C 3 (δ 1 , δ 2 , δ 3 )) contains a cumulative information of all third-order cumulants at the same time. Measuring these correlations in large and small systems can be used as an independent method of testing the collectivity and event-by-event fluctuation in heavy-ion experiments. We leave more investigation of this measurable quantity to future studies.

Conclusion
A multidimensional generating function method was introduced to extract a large class of cumulants related to the flow harmonic fluctuations. We proposed an ordering based on the Fourier-cumulant expansion for these observables to systematically capture the most dominant features of the flow harmonic fluctuations. Using this method and reproducing the already known cumulants, we have found new cumulants consisting of symmetry plane correlations that have not been studied before. We defined the normalized cumulants to compare the cumulant's magnitude with each other and with initial state fluctuation. We employed hydrodynamic simulation for Pb-Pb collisions (T R ENTo+VISH2+1+UrQMD) calibrated by a Global Bayesian analysis to predict the unmeasured normalized cumulants' value. The observables introduced in this study can be used as inputs for calibration of the Bayesian analysis or to validate the already tuned parameters.
We extract the linear and nonlinear hydrodynamic response by comparing the initial anisotropy and flow harmonic fluctuation cumulants. This method can be extended to obtain higher-order nonlinear terms by comparing more initial and final state fluctuations cumulants.
Based on the Fourier-cumulant expansion, we also introduced a general way to find the genuine q-particle distribution function, qPC. In particular, we studied 3PC containing information of all third-order cumulants for all harmonics. The flow-induced 3PC for harmonics n = 2, · · · , 5 was presented by using T R ENTo+VISH2+1+UrQMD simulation. This measurable helps to study flow and nonflow effects in large and small systems in the future. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: List of few first order cumulants for harmonics 2 to 6
By employing the Mathematica package introduced in Sect. 3, we tabulate the cumulants containing one and two harmonics up eight order. The cumulants harmonics containing three and four harmonics are presented in up to six orders. Except for the one-harmonic, the cumulants' final expression is mostly dependent on the values of n 1 , . . . , n k . There are still repeating patterns in the expressions, like c However, we sacrifice brevity for the sake of clarity and dedicate a separate table for each combination of flow harmonics regardless of repeating patterns (Tables 3, 4, 5 and 6). Table 3 One-harmonic cumulants up to order eight Table 4 Two-harmonic cumulants up to order eight Table 5 Three-harmonic cumulants up to order six is an estimator for x n . To this end, we compute the mean value of μ n using p.d.f. in Eq. (68), By substituting Eq. (68) into the above relation, one finds μ n = x n . The equality is exact if we redo the series of measurements (x 1 , . . . , x N ) infinite times and at each time compute the summation in Eq. (69), and finally compute the average over the results. Given that measuring a quantity infinite times is not practical, we would like to estimate the actual x n by doing the summation in Eq. (69) once with a reasonable number of measurements N . With the finite number of measurements, we get the following estimation x n = μ n + σ μ n , where σ μ n is the statistical uncertainty, sourced by the finite number of measurements. To calculate the statistical error of μ n , we compute μ 2 n as follows: Consequently, we obtain the variance of the random variable as where in the last line, we substitute x n with the estimation μ n . In summary, by doing N measurements, we estimate the true value of x n with statistical uncertainty σ μ n . There are similarities and differences between the above simple example and that we are looking for in heavy-ion physics, as we will explain in the following. The features of the collective models in heavy-ion physics can be categorized into two different parts. First one is the collective evolution part, which leads to a deterministic oneparticle distribution function f (ϕ) (Eq. (1)) for a given initial state. 12 And the second one is the initial state, which is fundamentally stochastic due to the quantum nature of the nucleus wave function. For that reason, we rewrite the Eq. (1) as follows: wherev n = v n e inψ n ,v 0 = 1, {v} ≡ {v 1 ,v 2 , . . .}, and v −n =v * n . The main difference between the above equation and that in Eq. (1) is the subscript {v} in f {v} (ϕ) to label the one-particle distribution function for different events with different flow harmonics in the final state.
Our physics of interest is encoded in f {v} (ϕ), and we would like to measure its Fourier coefficientsv n experimentally. Theoretically, the Fourier coefficients can be obtained by Experimentally, we have a distribution of M particles (M is the multiplicity) in the azimuthal direction, and we need to estimate the true value ofv n from the finite number of particles. This can be done by replacing the integral in Eq (74) with the following summation, where q n is called normalized flow vector. Focusing on a single event, we can skip subscript {v} for the moment. Then the probability of finding a specific configuration of particles (ϕ 1 , . . . , ϕ M ) in the final state is given by 12 The function f (ϕ)dϕ returns the probability of finding one particle in the interval (ϕ, ϕ + dϕ) in the azimuthal direction. similar to Eq. (68). In an event, the multiplicity M is a parameter related to our model's physical parameters. For instance, in hydrodynamic models, the multiplicity is related to hydrodynamic initial time and the initial energy density deposited into a given region in the transverse direction. As a result, there are cases with low multiplicity and interesting physics which normalized flow vectors are not helpful because of high statistical uncertainty. We could overcome this problem, by using many event averages. However, as we explained earlier, two different events have two different one-particle distribution function f {v} (ϕ), and consequently two different values forv n . It means we cannot simply increase the statistics by collecting many events. However, we are still able to extract statistically stable information about flow harmonics fluctuation. Referring to Eq. (73), the function f {v} (ϕ) is fully characterized by its Fourier coefficientsv n . As a result, the fluctuating f {v} (ϕ) can be encoded in an infinite dimensional p.d.f. shown by p(v ±1 ,v ±2 , . . .). Besides, the multiplicity M is different from one event to the other generally. To consider this point, we generalize the flow harmonic p.d.f. to In the above, we have used the following short notation, Considering the flow and multiplicity fluctuations, we should modify the joint p.d.f. in Eq. (76) into the following form: which consequently leads to In the present subsection, we elaborate on the algorithm's details behind the functions Nsigma2P. The experimental estimator of the moment v a 1 · · ·v a k , as we will find out, is related to the normalized flow vectors' product q a 1 · · · q a k . By referring to Eq. (75), we start with splitting this product as the following: where in the above, we have separated summations to the terms that all k particles are auto-correlated, k−1 particles are auto-correlated, up to the case that all the auto-correlations are removed. Before proceeding, let us rewrite the above summation in a more compact form by using set partitions. Consider that the set of all partitions of the set X is shown by P X , In the analogy of P, we define the operator I and A acting on the set {i 1 , . . . , i k } and {a 1 , . . . , a k } as where B k is the Bell number corresponds to the number of all partitions of a set with k elements. Although these definitions look complicated, they help us to rewrite the Eq. (81) more compactly, Like what has been done in Eq. (70), we perform the average q a 1 · · · q a k with respect to F({ϕ}). To this end, several integrations have to be done. For example Similarly, and by employing notations in Eqs. (83) and (84), we obtain q a 1 · · · q a k = dϕ 1 dϕ 2 · · · (q a 1 · · · q a k ) F({ϕ}), In the above, the summation on I j is trivial. To find out the final expression, we need to calculate the number of terms in this summation. In fact, by noting that each i j runs from 1 to M, it is easy to see that where the (M) i = M(M − 1) · · · (M − i + 1) is the falling factorial, and I i is the number of elements in the set I j .
As a result, we eventually find where we have used the fact that I j = A j . The last expression is rather simple. To make it more clear, we elaborate on it with two examples in the following. First consider a 1 = n and a 2 = −n. In this case, which leads to q n q −n = 1 The correlations on the right-hand side are moments of p(M,v ±1 ,v ±2 , . . .). Instead of q n q −n , we could have started with q n q −n (M/M − 1) . In such a case, we would have obtained, By rearranging the above relation and using q n = Q n /M (see Eq. (30)), we recover Eq. (32). For the second example, assume a 1 = 2, a 2 = 3, and a 3 = −5. In this special case, we have which consequently leads to By replacing v 2 n from Eq. (92) and starting from  [36,62]. Given that the set I B k = {i 1 = · · · = n k } with I B k = k is the unique most populated set in the set of partitions I, only the term j = B k survives the limit M → ∞ in Eq. (89). As a result, we obtain lim M→∞ q a 1 · · · q a k = v a 1 · · ·v a k , which is manifestly correct because, in the limit M → ∞, there is no statistical fluctuation in a single event. This indicates that the true value of q a 1 · · · q a k with no statistical error is encoded in the sector I B k . Referring to Eq. (87), and keep only the term j = B k from both sides, we obtain, v a 1 · · ·v a k = 1 The meaning of Eq. (97) is that by removing the autocorrelations and averaging over many events (as it is done in Refs. [34][35][36]62]) the effect of statistical fluctuations in a single event is removed and we obtain accurate moments of The above statement is valid only when the average · · · is performed over infinitely many numbers of events. Similar to computations lead to Eq. (72), we compute the statistical fluctuation of the right-hand side of Eq. (97) with the finite number of events, N . Defining and Dϕ = dϕ 1 dϕ 2 · · · , one can compute the following integrals, The quantity k has been computed in Eq. (97) while we need to elaborate k . 13 To this end, we need to study a summation as Considering the set {i 1 , i 2 , . . . , i k+ }, we immediately find that in reorganizing the summation in Eq. (100) no partition with more than two elements is allowed. Assume a set as {i 1 , x} where the element x can be any i j with j = k + 1, . . . , . Add any third element to {i 1 , x} contradicts with one of the inequalities i 1 = · · · = i k or i k+1 = · · · = i in Eq. (100). As a result, the partitions are started from sets as Here we conventionally assumed ≥ k. The -tuple (s 1 , . . . , s ) is a permutation of (k + 1, . . . , k + ) which are ! different kinds. However, ( − k) elements of thetuple are placed in a set with a single element. Therefore, the permutation of that subset does not lead to a different partition in Eq. (101). As a result, the number of partitions in Eq. (101) is !/(k − )!. We are allowed to break each set {i r , i s r } into two, and the result will be a legal partition. After breaking all sets as {i r , i s r } we reach to Similar to i j partitions, the partitions of harmonics a i quantities start from a s 1 , . . . , a k + a s k , a k+1 , . . . , a } where (s 1 , . . . , s ) = perm(1, . . . , ) and can be continued by splitting each summation a i + a s i into two. Defining the set as set of all f k, partitions of set (103), we obtain a relation similar to Eq. (89), Now, we can define the covariance matrix as follows The statistical fluctuation of k is given by cov( k , k ). For clarity, let us elaborate on a simple example with a 1 = n and a 2 = −n. The partition set is given by Referring to Eqs. (97) and (98), one notes thatv a 1 · · ·v a k = k a 1 ,...,a k . Substituting this relation into the above equation, we obtain Eq. (45), compatible with [36,64] if we ignore the flow fluctuation and set ψ n = ψ 2n .

Appendix C: Flow fluctuation and particle distribution decomposition
This appendix introduces a generic procedure to connect qth order cumulants to the qPC. We elaborate 3PC for n = 2, . . . , 5 explicitly at the end of the appendix.

C.1. Nonflow correlations vs. flow-induced correlations
One of the main motivations of using cumulants in studying the flow harmonic fluctuations is removing the nonflow effects. Moreover, cumulants systematically classify the deviation of the flow harmonic p.d.f. (2) from Gaussianity. In this part, we would like to show how much the observed features of these two concepts are entangled. Ignoring the flow fluctuation, Eq. (76) should be replaced by [34] F(ϕ 1 , . . . , ϕ M ) = f (ϕ 1 ) · · · f (ϕ M ) +F nonflow (ϕ 1 , . . . , ϕ M ), where F nonflow (ϕ 1 , . . . , ϕ M ) is correlation developed by nonflow effects in a single event. One notes that F nonflow (ϕ 1 , . . . , ϕ M ) still contains products of two-particles or more joint distributions. For instance A standard way to suppress the nonflow effects is using cumulants. By measuring mth order cumulant we practically measure a genuine correlation among mth particles. Therefore by measuring the mth order cumulant of the distribution, we automatically remove any n-particle nonflow correlations with n < m [34,35].
In the presence of flow fluctuation, the same structure as Eq.
(109) appears even if we ignore nonflow effects. In such a case, Eq. (80) is written as the following (for the detail see Appendix C) where f (ϕ) is the many events average of the single-particle distribution Eq. (1) and F flow-fluc (ϕ 1 , . . . , ϕ M ) is the joint correlation part appears because of averaging over many fluctuating events. As a result, there are two entangled parts in genuine multiparticle correlations: nonflow effects in a single event develop the first part, and the second part fictitiously appears because of averaging over many fluctuating events.
To study lower-order cumulants, which contain nontrivial and important information about flow fluctuations, we need to use other methods to reduce the nonflow effects [24,79].
C.2. Genuine q-particle correlation functions By inserting Eq. (73) into Eq. (80), we find the following expression, 14 The above relation cannot be decomposed into the product of M distinct single-particle distribution. However, we are still able to rewrite the moment v a 1 · · ·v a M as a summation of all possible genuine correlations betweenv a i s (cumulants · · · c ), v a 1 · · ·v a M = v a 1 c · · · v a M c + v a 1v a 2 c v a 3 c · · · v a M + · · · . . .
We see from above that the apparent correlation induced by fluctuation depends on the event by event fluctuation. Considering the rotational symmetry, the constraint a 1 + · · · + a k = 0 should be included into the above summation. By this constraint, the cumulant v a 1 · · ·v a q c is equivalent to that we have already studied in Sect. 2. Defining 14 Here, we ignore the multiplicity fluctuation.
The generating function of the above distribution is written as G(P, Q) = e i L , where Similar to Eq. (122), we have used the notation P = (k n 1 ,x , . . . , k n k ,x ), Q = (k n 1 ,y , . . . , k n k ,y ).
Presenting Eq. (123) in polar coordinate and in an arbitrary reaction plane angle RP , we find where we have used polar coordinate v n e inψ n = v n,x + iv n,y together with k n,x = k n cos nφ n , k n,y = k n sin nφ n .
We eventually fix our reference frame such that φ n 1 = RP to obtain Eq. (16).