Higher order Symmetric Cumulants

TL;DR

This paper introduces higher order symmetric cumulants as new observables to study complex correlations among flow amplitudes in high-energy nuclear collisions, providing deeper insights into the Quark-Gluon Plasma.

Contribution

It generalizes symmetric cumulants with higher order observables and demonstrates their feasibility and potential for revealing new information about initial conditions and QGP properties.

Findings

Feasibility demonstrated with realistic models

Predictions for centrality dependence in Pb--Pb collisions

New observables contain unique information beyond pairwise correlations

Abstract

We present the generalization of recently introduced observables for the studies of correlated fluctuations of different anisotropic flow amplitudes, dubbed Symmetric Cumulants. We introduce a new set of higher order observables and outline a unique way how the genuine multi-harmonic correlations can be extracted from multi-particle azimuthal correlators. We argue that correlations among flow amplitudes can be studied reliably with the general mathematical formalism of cumulants only if that formalism is applied directly on the flow amplitudes. We have tested all the desired properties of new observables with the carefully designed Toy Monte Carlo studies. By using the realistic iEBE-VISHNU model, we have demonstrated that their measurements are feasible and we have provided the predictions for the centrality dependence in Pb--Pb collisions at LHC energies. A separate study was performed for their values in the coordinate space. The new observables contain information which is inaccessible to individual flow amplitudes and correlated fluctuations of only two flow amplitudes, and therefore they provide further and independent constraints for the initial conditions and the properties of Quark-Gluon Plasma in high-energy nuclear collisions.

Figures

Equations

$$f(\phi )=\frac{1}{2\pi }\left[1+2\sum _{n=1}^{\infty }{v}_n\cos [n(\phi -{\Psi }_n)]\right].$$

$$⟨e^{i(n_1{\phi }_1+\cdots +n_k{\phi }_k)}⟩=v_{n_1}\cdots v_{n_k}{e}^{i(n_1{\Psi }_{n_1}+\cdots +n_k{\Psi }_{n_k})},$$

$${⟨⟨\cos [n({\phi }_1+{\phi }_2-{\phi }_3-{\phi }_4)]⟩⟩}_c$$

$${⟨⟨\cos (m{\phi }_1+n{\phi }_2-m{\phi }_3-n{\phi }_4)⟩⟩}_c$$

$${⟨⟨\cos [n({\phi }_1+{\phi }_2-{\phi }_3-{\phi }_4)]⟩⟩}_c$$

$${⟨⟨\cos (m{\phi }_1+n{\phi }_2-m{\phi }_3-n{\phi }_4)⟩⟩}_c$$

$$\mathrm{SC}(k,l,m)\equiv ⟨v_k^2{v}_l^2{v}_m^2⟩-⟨v_k^2{v}_l^2⟩⟨v_m^2⟩-⟨v_k^2{v}_m^2⟩⟨v_l^2⟩-⟨v_l^2{v}_m^2⟩⟨v_k^2⟩+2⟨v_k^2⟩⟨v_l^2⟩⟨v_m^2⟩.$$

$$\mathrm{SC}(k,l,m)$$

$$⟨\cos [n(3{\phi }_1+2{\phi }_2+{\phi }_3-3{\phi }_4-2{\phi }_5-{\phi }_6)]⟩=v_n^2{v}_{2n}^2{v}_{3n}^2.$$

$${⟨\cos [n(3{\phi }_1+2{\phi }_2+{\phi }_3-3{\phi }_4-2{\phi }_5-{\phi }_6)]⟩}_c$$

$${⟨v_n^2{v}_{2n}^2{v}_{3n}^2⟩}_c$$

$$f(\phi )=\frac{1}{2\pi }\left[1+2v_2\cos \left(2\phi \right)+2v_3\cos \left(3\phi \right)+2v_4\cos \left(4\phi \right)\right].$$

$$M_{\text{final}}=2M_{\text{initial}}.$$

$${\delta }_3^{\mathrm{SC}}=\frac{\alpha }{M^5}+\frac{\beta }{M^4}+\frac{\gamma }{M^3},$$

$$\begin{array}{rl}\alpha & =27.01\pm 9.79,\\ \beta & =-0.0947\pm 0.1743,\\ \gamma & =(-1.383\pm 5.876)\cdot 10^{-4},\end{array}$$

$${ϵ}_n{e}^{in{\Phi }_n}\equiv -\frac{\int r^n{e}^{in\phi }\rho (r,\phi )rdrd\phi }{\int r^n\rho (r,\phi )rdrd\phi },n=2,3,\dots ,$$

$$\text{NSC}(n,m,p)=\frac{\text{SC}(n,m,p)}{⟨v_n^2⟩⟨v_m^2⟩⟨v_p^2⟩},{\text{NSC}}_{ϵ}(n,m,p)=\frac{{\text{SC}}_{ϵ}(n,m,p)}{⟨{ϵ}_n^2⟩⟨{ϵ}_m^2⟩⟨{ϵ}_p^2⟩}.$$

$$v_4{e}^{4i{\Psi }_4}$$

$$v_5{e}^{5i{\Psi }_5}$$

$$⟨e^{in({\phi }_1-{\phi }_2)}⟩=⟨\cos [n({\phi }_1-{\phi }_2)]⟩=v_n^2,$$

$$f({\phi }_1,{\phi }_2)=f_{{\phi }_1}({\phi }_1)f_{{\phi }_2}({\phi }_2).$$

$$⟨e^{in({\phi }_1-{\phi }_2)}⟩=v_n^2+{\delta }_2,$$

$$⟨e^{in({\phi }_1-{\phi }_2)}⟩={\delta }_2,$$

$$⟨e^{i(n_1{\phi }_1+\dots +n_{k/2}{\phi }_{k/2}-n_1{\phi }_{k/2+1}-\dots -n_{k/2}{\phi }_k)}⟩$$

$$\mathrm{SC}(m,n{)}_{\mathrm{correct}}$$

$$\mathrm{SC}(m,n{)}_{\mathrm{wrong}}$$

$$f({\phi }_1,\dots ,{\phi }_k)\ne f({\phi }_1)\cdots f({\phi }_k).$$

$${\delta }_k\sim \frac{1}{M^{k-1}}.$$

$${\delta }_k^{\mathrm{SC}}\sim {\delta }_{2k}+\cdots +\prod _{i=1}^k{\delta }_2\sim \sum _{i=k}^{2k-1}\frac{1}{M^i}.$$

$${\delta }_3^{\mathrm{SC}}=\frac{\alpha }{M^5}+\frac{\beta }{M^4}+\frac{\gamma }{M^3},$$

$$⟨⟨\cdots ⟩⟩=\frac{{\sum }_i{w}_i{⟨\cdots ⟩}_i}{{\sum }_i{w}_i},$$

$$\mathrm{SC}(k,l,m)=⟨v_k^2⟩⟨v_l^2⟩⟨v_m^2⟩-⟨v_k^2⟩⟨v_l^2⟩⟨v_m^2⟩-⟨v_k^2⟩⟨v_m^2⟩⟨v_l^2⟩-⟨v_l^2⟩⟨v_m^2⟩⟨v_k^2⟩+2⟨v_k^2⟩⟨v_l^2⟩⟨v_m^2⟩=0.$$

Full text

I Introduction

In the explorations of the phase diagram of strong nuclear force, some of the most intriguing questions are associated with the phase dubbed Quark-Gluon Plasma (QGP), in which under large temperatures and/or baryon densities quarks and gluons are deconfined. Properties of this extreme state of matter have been investigated in high-energy physics in the past 20 years with the plethora of different observables and across different collision systems, at SPS, RHIC and LHC experiments. One of the most important programs in these studies were the analyses of anisotropic flow phenomenon Ollitrault:1992bk ; Voloshin:2008dg , which primarily have been carried out with the two- and multi-particle correlation techniques. The anisotropic flow measurements proved to be particularly successful in the studies of transport properties of QGP. For instance, they were used to constrain the ratio of shear viscosity to entropy density ($\eta/s$) of QGP to the very small values and therefore helped to establish a perfect liquid paradigm about the QGP properties Braun-Munzinger:2015hba .

Traditionally, anisotropic flow is quantified with the amplitudes $v_{n}$ and symmetry planes $\Psi_{n}$ in the Fourier series parametrization of anisotropic particle emission as a function of azimuthal angle $\varphi$ in the plane transverse to the beam direction Voloshin:1994mz :

[TABLE]

Correlation techniques have been utilized in the past to provide estimates for the average values of various moments of individual flow amplitudes, $\left<v_{n}^{k}\right>$ ($k>1$), where each moment by definition carries a piece of independent information about the event-by-event fluctuations, i.e. the stochastic nature, of anisotropic flow. Flow fluctuations are unavoidable in ultrarelativistic collisions as they originate from the non-negligible fluctuations of positions of the nucleons inside the colliding nuclei, as well as from quantum fluctuations of the quark and gluon fields inside those nucleons. As a consequence, even the nuclear collisions with the same impact parameter (the vector connecting centers of two colliding nuclei) can generate different anisotropic flow event-by-event. Therefore, a great deal of additional insights about QGP properties can be extracted from flow fluctuations. It was demonstrated recently that more detailed information of QGP properties, like temperature dependence of $\eta/s$, cannot be constrained with the measurements of individual flow amplitudes due to their insensitivity. Instead, a set of more sensitive flow observables, focusing on correlations of fluctuations of two different flow amplitudes $v_{n}$ and $v_{m}$, have been proposed Niemi:2012aj ; Bilandzic:2013kga . In this paper, we generalize this idea and introduce a new set of independent flow observables, which quantify the correlations of flow fluctuations involving more than two flow amplitudes.

By using solely the orthogonality relations of trigonometric functions, one can from Eq. (1) derive that $v_{n}=\left<\cos[n(\varphi\!-\!\Psi_{n})]\right>$, where the average $\left<\cdots\right>$ goes over all azimuthal angles $\varphi$ of particles reconstructed in an event. However, this result has little importance in the measurements of flow amplitudes $v_{n}$, since symmetry planes $\Psi_{n}$ cannot be measured reliably in each event. Since azimuthal angles $\varphi$, on the other hand, can be measured with the great precision, we can estimate instead the flow amplitudes $v_{n}$ and the symmetry planes $\Psi_{n}$ by using the correlation techniques Wang:1991qh ; Jiang:1992bw . The cornerstone of this approach is the following analytic result derived recently Bhalerao:2011yg

[TABLE]

where the true value of the average $\left<\cdots\right>$ can be estimated experimentally by averaging out all distinct tuples of $k$ different azimuthal angles $\varphi$ reconstructed in the same event. With the advent of a generic framework for flow analyses with multi-particle azimuthal correlations Bilandzic:2013kga , the LHS of Eq. (2) can be evaluated fast and exactly for any number of particles $k$, and for any choice of non-zero integers for harmonics $n_{1},\ldots,n_{k}$. This new technology was successfully used in the first measurements of Symmetric Cumulants (SC) in ALICE:2016kpq and can be used also for their generalized version, which we present in this paper.

Before elaborating further on all the technical points which are relevant for the generalization of SC observables, we briefly review the existing studies of correlated fluctuations of different flow amplitudes, with the special focus on the ones performed with SC, both by theorists and experimentalists. We traced back the first theoretical study in which correlations of amplitude fluctuations of two different flow harmonics were used to extract the independent information about the QGP properties to Niemi:2012aj . The observables utilized in that work are available only when it is feasible to estimate flow amplitudes event-by-event, which is the case only in the theoretical studies. Alternatively, experimentalists use the correlation techniques to estimate the correlated fluctuations of amplitudes of two different flow harmonics, via the SC observables, in the way it was first proposed in the Sec. IVC of Bilandzic:2013kga . Complementary, such correlated fluctuations can be also probed with the Event Shape Engineering (ESE) technique Schukraft:2012ah . Even though the SC are relatively novel observables, a lot of theoretical studies utilizing them have been already performed: the standard centrality dependence with state-of-the-art models was obtained in Bhalerao:2014xra ; Zhou:2015eya ; Gardim:2016nrr ; Zhao:2017yhj ; Eskola:2017imo ; Alba:2017hhe ; Gardim:2017ruc ; Giacalone:2018wpp ; Bernhard:2018hnz ; Moreland:2018gsh ; Schenke:2019ruo ; Moreland:2019szz ; Wen:2020amn ; the relation between SC and symmetry plane correlations was studied in Giacalone:2016afq ; the more differential studies (including transverse momentum and pseudorapidity dependence, or using subevents) were performed in Zhu:2016puf ; Ke:2016jrd ; Nasim:2016rfv ; Jia:2017hbm ; the use of SC to constrain the details of energy density fluctuations in heavy-ion collisions was recently discussed in Bhalerao:2019fzp ; extensive coordinate space study was reported in Broniowski:2016pvx ; the study in the collisions of smaller systems was carried out in Dusling:2017dqg ; Dusling:2017aot ; Albacete:2017ajt ; Nie:2018xog ; Sievert:2019zjr ; Rybczynski:2019adt ; Zhao:2020pty . The complementary theoretical studies of $v_{n}-v_{m}$ correlations, without using SC, has been performed in Qian:2016pau ; Qian:2017ier . The nontrivial technical details about the implementation of multi-particle cumulants, with the special mentioning of SC, was recently rediscussed from scratch in DiFrancesco:2016srj . How SC can emerge in the broader context of flow fluctuation studies was briefly mentioned in Mehrabpour:2018kjs .

Measurements of correlations between the amplitude fluctuations of two different flow harmonics in heavy-ion collisions with the ESE technique was first reported by ATLAS in Jia:2014jca ; Aad:2015lwa ; Radhakrishnan:2016tsc . On the other hand, analogous measurements by using SC observables were first reported by ALICE in Zhou:2015slf ; ALICE:2016kpq ; Kim:2017rfn . After these initial studies, the measurements of SC have been successfully extended to different collision systems and energies in Guilbaud:2017alf ; Sirunyan:2017uyl ; STAR:2018fpo ; Aaboud:2018syf ; Gajdosova:2018pqo ; Acharya:2019vdf ; Aaboud:2019sma ; Sirunyan:2019sef . The detailed differential analyses and the extension to correlations of subleading flow harmonics were published in Acharya:2017gsw . Feasibility of measurements of SC in the future experiments was recently addressed in Citron:2018lsq .

The rest of the paper is organized as follows. In Sec. II the generalization of SC observables is motivated and presented. In Sec. III the detailed Monte Carlo studies for the higher order SC are presented, and first predictions provided for heavy-ion collisions at LHC energies. We have summarized our findings in Sec. IV. In a series of appendices, we have placed all technical steps which were omitted in the main part of the paper.

II Generalization to higher orders

The possibility of higher order SC was first mentioned in Table 1 of Jia:2014jca , while the first systematic studies of their properties were performed independently in Sec. 3.2 of Deniz:2017 . In this section, we elaborate on the generalization of SC observables for the concrete case of three flow amplitudes, but our results and conclusions can be straightforwardly generalized to any number of amplitudes by following the detailed list of requirements presented in Appendix A. The few final solutions for generalized SC involving more than three amplitudes we have outlined in Appendix B. We have cross-checked all our statements about higher order SC with the set of carefully designed Toy Monte Carlo studies in Appendix C. In Appendix F we have outlined the argument that there are no built-in mathematical correlations in SC observables which might originate as a consequence of Parseval theorem. In all our derivations, we use as the starting point the general mathematical formalism of cumulants published in Kubo .

II.1 Physics motivation

We now briefly motivate this generalization with the following simple observation: In mid-central heavy-ion collisions the correlated fluctuations of even amplitudes $v_{2},v_{4},v_{6},\ldots$ can solely originate from the fluctuations of magnitude in the ellipsoidal shape Kolb:2003zi . This sort of fluctuations will contribute only to SC(2,4), SC(2,4,6), etc., but not to SC(2,3,4). The genuine correlation among three amplitudes, $v_{2}$, $v_{3}$ and $v_{4}$ can develop only due to fluctuations in the elliptical shape itself. The non-vanishing result for SC(2,3,4) implies that there are additional sources of fluctuations in the system which couple all three amplitudes $v_{2},v_{3},v_{4}$. In this sense, SC(2,3,4) can separate the following two sources of fluctuations: a) fluctuations in the shape of ellipsoidal; b) magnitude fluctuations of the persistent ellipsoidal shape. On the other hand, SC(2,4) cannot disentangle these two different sources of fluctuations. This argument is supported by concrete calculus in the simple mathematical model presented in Appendix E.

II.2 Technical details

Before going to the generalization of SC, we recall here a few important points about the use of multi-particle azimuthal correlations in flow analyses. In general, these correlations are susceptible to the non-trivial collective phenomena (e.g. anisotropic flow), but also to the phenomena involving only a few particles (e.g. momentum conservation, jet fragmentation, resonance decays, etc.). The latter, called nonflow, is considered as systematic bias in flow analyses. One way to separate the two contributions is to compute the genuine multi-particle correlations or cumulants.

In the general mathematical formalism of cumulants Kubo , for any set of $k$ ($k>1$) stochastic observables, $X_{1},\ldots,X_{k}$, there exists a unique term, the $k$-particle cumulant, sensitive only to the genuine $k$-particle correlation among all the observables. This mathematical formalism has been introduced for anisotropic flow analyses in Borghini:2000sa ; Borghini:2001vi , and further improved and generalized in Bilandzic:2010jr ; Bilandzic:2013kga . As one has the freedom to choose the experimental observables of interest (e.g. azimuthal angles, flow amplitudes $v_{n}$, multiplicities of particles in a given momentum bin DiFrancesco:2016srj , etc.), different variants of cumulants can appear in practice. In flow analyses, the traditional choice is to identify $X_{i}\equiv e^{in_{i}\varphi_{i}}$, with $\varphi_{i}$ the azimuthal angles of reconstructed particles and $n_{i}$ the non-zero flow harmonics, as the observable of interest. This specific choice allows the identification of the averages in the cumulant expansion (e.g. Eq. (2.8) in Kubo ) with the single-event averages of azimuthal correlators (Eq. (2)). If the detector has a uniform azimuthal acceptance, only the isotropic correlators with $\sum_{i}n_{i}$ = 0 are not averaged to zero in the generalization of the single-event averages $\left<\cdots\right>$ to all-event averages $\left<\left<\ldots\right>\right>$. Finally, in the resulting expression one groups together all terms which differ only by the re-labelling of azimuthal angles $\varphi_{i}$ Borghini:2000sa ; Borghini:2001vi .

For instance, the flow amplitude $v_{n}$ estimated with four-particle cumulant, $v_{n}\{4\}$, can be obtained by identifying $X_{1}\equiv e^{in\varphi_{1}}$, $X_{2}\equiv e^{in\varphi_{2}}$, $X_{3}\equiv e^{-in\varphi_{3}}$ and $X_{4}\equiv e^{-in\varphi_{4}}$. Following the steps described above, one obtains

[TABLE]

A generalization of this idea for the case of non-identical harmonics leads to the definition of new observables dubbed SC Bilandzic:2013kga . With the more general choice $X_{1}\equiv e^{in\varphi_{1}}$, $X_{2}\equiv e^{im\varphi_{2}}$, $X_{3}\equiv e^{-in\varphi_{3}}$, and $X_{4}\equiv e^{-im\varphi_{4}}$, where $n$ and $m$ are two different positive integers, the general formalism of four-particle cumulants from Kubo translates now into:

[TABLE]

It follows that, in both Eqs. (3) and (4) above, the final expressions for the cumulants depend only on the flow amplitudes, since by using Eq. (2) we obtain immediately:

[TABLE]

The success of observables $v_{n}\{4\}$ and SC$(m,n)$ lies in the fact that they suppress much better the unwanted nonflow correlations, than for instance the starting azimuthal correlators $\left<\left<\cos[n(\varphi_{1}\!+\!\varphi_{2}\!-\!\varphi_{3}-\!\varphi_{4})]\right>\right>$ and $\left<\left<\cos(m\varphi_{1}\!+\!n\varphi_{2}\!-\!m\varphi_{3}\!-\!n\varphi_{4})\right>\right>$, and therefore they are much more reliable estimators of anisotropic flow properties. We will elaborate more on systematic biases due to nonflow later in the paper.

After these general considerations, to which we refer from this point onward as traditional or standard cumulant expansion in flow analyses, we proceed with the generalization of SC. As a concrete example, we provide the generalization of SC observables for the case of three flow amplitudes, and discuss in detail all consequences. As opposed to the traditional cumulant expansion in Eqs. (3) and (4), we now directly define SC$(k,l,m)$ by using the general cumulant expansion for three different observables (Eq. (2.8) in Kubo ), in which as fundamental observables, instead of azimuthal angles, we choose directly the amplitudes of flow harmonics. This difference in the choice of fundamental observables on which the cumulant expansion is performed is the landmark of our new approach. It follows immediately:

[TABLE]

In theoretical studies, in which flow amplitudes can be computed on an event-by-event basis, this expression suffices. In Appendix A we demonstrate that SC$(k,l,m)$ satisfies the fundamental cumulant properties, while the remaining requirements needed for generalization are cross-checked in Appendix C.

In an experimental analysis, it is impossible to estimate reliably flow amplitudes event-by-event, and all-event averages of azimuthal correlators need to be used to estimate them. This connection is provided by the analytic result in Eq. (2). Solely by using that result, it follows immediately that the SC$(k,l,m)$ observable defined in Eq. (6) can be estimated experimentally with:

[TABLE]

Since each harmonic in each azimuthal correlator above appears with an equal number of positive and negative signs, any dependence on the symmetry planes $\Psi_{n}$ is canceled out by definition in each of the ten correlators (Requirement 4 in Appendix A, and proof therein). All correlators in Eq. (7) are invariant under the shift $\varphi_{i}\rightarrow\varphi_{i}+\alpha$ of all azimuthal angles, where $\alpha$ is arbitrary, which proves their isotropy (Requirement 5 in Appendix A).

We have arrived at the final result in Eq. (7) by following our new approach in which the cumulant expansion has been performed directly in Eq. (6) on flow amplitudes. While we have started with the cumulant expansion of order 3 in Eq. (6), the final expression in Eq. (7) depends on 6 azimuthal angles. We now make a comparison with the traditional approach in which the cumulant expansion is performed directly on 6 azimuthal angles $\varphi$, and discuss the differences.

II.3 Difference between old and new cumulant expansion in flow analyses

We now scrutinize the cumulant expansion for one concrete example of six-particle azimuthal correlators, by following the traditional procedure, in which the starting observables are azimuthal angles. We start first with the analytic expression (derived solely using Eq. (2)):

[TABLE]

The corresponding cumulant in the traditional approach is:

[TABLE]

In terms of flow amplitudes, this expression evaluates into:

[TABLE]

The result in Eq. (10) is not a valid cumulant of flow amplitudes $v_{n}$, $v_{2n}$ and $v_{3n}$, since it has the extra terms in the 2nd line which depend also on the symmetry planes. Because of this contribution, this expression does not reduce to 0 for three independent flow amplitudes $v_{n}$, $v_{2n}$ and $v_{3n}$, which violates the elementary Theorem 1 on multivariate cumulants from Kubo .

Based on this concrete example, we conclude that one cannot start in general to calculate cumulants for one set of observables (e.g. azimuthal angles $\varphi_{1},\varphi_{2},\ldots$), and then interpret the final results to be a cumulant of some other observables (e.g. $v_{n}$’s and/or $\Psi_{n}$’s). At best one can state that the cumulants of multi-variate distribution of azimuthal angles $P(\varphi_{1},\varphi_{2},\ldots)$ can be parameterized with $v_{n}$’s and $\Psi_{n}$’s (since they are related via Fourier series expansion in Eq. (1) and analytic expression in Eq. (2), but one cannot state that the final results are direct cumulants of $v_{n}$’s and/or $\Psi_{n}$’s. In particular, the cumulants of flow degrees of freedom $v_{n}$’s and $\Psi_{n}$’s can be obtained only from the underlying p.d.f. $P(v_{1},v_{2},\ldots,\Psi_{1},\Psi_{2},\ldots)$ which governs the stochastic nature of $v_{n}$’s and $\Psi_{n}$’s observables, and not from $P(\varphi_{1},\varphi_{2},\ldots)$, which governs the stochastic nature of particle azimuthal angles.

In our new approach, we use azimuthal angles merely to estimate flow amplitudes via the mathematical identity in Eq. (2), but only after the cumulant expansion has been already performed directly on flow amplitudes $v_{m},v_{n},\ldots$. We do not perform the cumulant expansion on azimuthal angles directly, as it was done in flow analyses using cumulants so far (e.g. in Sec. II of Borghini:2001vi ). This subtle difference is one of the main points in our paper. This is a necessary change if one wants to estimate reliably genuine correlations of flow amplitudes, with the measured multi-particle azimuthal correlators.

This change also implies different expressions for cumulants of just one flow amplitude. For instance, the variance (second cumulant) of $v_{n}^{2}$ is $\left<v_{n}^{4}\right>-\left<v_{n}^{2}\right>^{2}$, which is not the same as the usual 4-particle cumulant which gives $\left<v_{n}^{4}\right>-2\left<v_{n}^{2}\right>^{2}$.

We now support our conclusions further by presenting the Monte Carlo studies.

II.4 Comparison of old and new cumulant expansion with Monte Carlo studies

To compare results for cumulants calculated by treating azimuthal angles as fundamental observables (tagged ‘old’ in this section) and calculated with our new approach by treating flow amplitudes as fundamental observables (tagged ‘new’), we use as a concrete example the respective results for the higher order SC built out of the flow amplitudes $v_{1},v_{2}$ and $v_{3}$. In the traditional approach, we use result in Eq. (10) by setting $n=1$, and in the new approach we use Eq. (6) by setting $k=1,l=2,m=3$. Despite this specific choice, the outcome of this study is more general given the generic nature of Eq. (10), which covers all cases when the order of one harmonic is the sum of orders of two remaining harmonics.

The comparison of the centrality dependence of these two expressions for the realistic VISHNU model (the details of VISHNU setup are presented later in Sec. III) is presented in Fig. 1(a), while in Fig. 1(b) we have presented only the centrality dependence of the 2nd line in Eq. (10), i.e. the difference between ‘new’ and ‘old’ approach. Clearly, in this concrete and realistic example, this difference is not negligible in mid-central and in peripheral collisions.

The remaining question is to demonstrate unambiguously which of the two expressions above is a valid cumulant of flow amplitudes. To answer this question we have set up the following clear-cut Toy Monte Carlo study, which demonstrates that only the Eq. (6) is the valid cumulant of three flow amplitudes $v_{1}$, $v_{2}$ and $v_{3}$: a) The amplitudes $v_{1}$, $v_{2}$ and $v_{3}$ are sampled in each event uniformly and independently in the intervals (0.03,0.1), (0.04,0.1), and (0.05,0.1), respectively; b) The symmetry planes $\Psi_{1}$ and $\Psi_{2}$ are sampled in each event uniformly and independently in the interval (0,2$\pi$); c) The symmetry plane $\Psi_{3}$ is determined in each event as: $\Psi_{3}\equiv\frac{1}{3}(\frac{\pi}{4}+2\Psi_{2}+\Psi_{1})$. With such a setup, in each event $\cos(3\Psi_{3}-2\Psi_{2}-\Psi_{1})=\cos\frac{\pi}{4}=\frac{\sqrt{2}}{2}$ (and similarly for the sine term), and therefore the 2nd line in Eq. (10) is non-vanishing, which means that in this study Eqs. (10) and (6) will yield systematically different results. However, since the three flow amplitudes $v_{1}$, $v_{2}$ and $v_{3}$ are sampled independently, the corresponding cumulant must be 0 (Theorem 1 in Kubo ). After we have carried out this study for 10000 events, we have obtained the results presented in Fig. 2. Since only $\left<v_{1}^{2}v_{2}^{2}v_{3}^{2}\right>_{c,{\rm new}}$ is consistent with 0 (see Fig. 2(b)), we conclude that this observable is a valid cumulant of three flow amplitudes $v_{1}$, $v_{2}$ and $v_{3}$. This simple Toy Monte Carlo demonstrates also that another observable, $\left<v_{1}^{2}v_{2}^{2}v_{3}^{2}\right>_{c,{\rm old}}$, is not a valid cumulant of the three flow amplitudes $v_{1}$, $v_{2}$ and $v_{3}$ due to the presence of 2nd line in Eq. (10) involving symmetry planes, which spoils the cumulant property.

Before presenting our first predictions for higher order SC in heavy-ion collisions at LHC energies, we discuss the systematic bias from nonflow contribution.

II.5 Nonflow estimation with Toy Monte Carlo studies

In this section, we discuss the generic scaling of the nonflow contribution in higher order SC. We begin with the implementation of the Fourier distribution $f(\varphi)$ parameterized as in Eq. (1), which we use to sample the azimuthal angle of each simulated particle. For simplicity, we define $f(\varphi)$ with three independent parameters, the flow amplitudes $v_{2}$, $v_{3}$ and $v_{4}$:

[TABLE]

The setup of our Toy Monte Carlo simulations goes as follows: For each one of the input number of events $N$, we set the values of the multiplicity $M$ and the flow harmonics $v_{2}$, $v_{3}$ and $v_{4}$ according to the requirements of the current analysis. Examples of such conditions can be that these input parameters are kept constant for all events or that they are uniformly sampled event-by-event in given ranges. We indicate this second case with the notation ($\cdot$,$\cdot$). After insertion of the harmonics, Eq. (11) is used to sample the azimuthal angles of the $M$ particles. We then compute all the needed azimuthal correlators with the generic framework introduced in Bilandzic:2013kga with the possibility left open to choose which event weight to use during the transition from single- to all-event averages. Finally our SC observable is obtained using Eq. (7).

Now that our Toy Monte Carlo setup is in place, we can use it to check if our SC observable has the needed properties. We simulate one example of simple nonflow correlations and look at the scaling of our expression as a function of the multiplicity. Our nonflow is described as follows: for each one of the $N=10^{8}$ events we generate, we sample a fixed initial number of particles $M_{\textrm{initial}}$ amongst the possibilities: 25, 38, 50, 75, 100, 150, 200, 250 and 500. The flow harmonics are all set to zero. We then introduce strong two-particle correlations by taking each particle two times in the computation of the two- and multi-particle correlations. This means our final multiplicity is given by

[TABLE]

As detailed in the list of requirements in Appendix A, the nonflow scaling of our three-harmonic SC can be described with the following expression:

[TABLE]

where $M$ corresponds to the final multiplicity introduced in Eq. (12). The fit is done over all simulated multiplicities and can be found in Fig. 3. The obtained fit parameters are as follows:

[TABLE]

with a goodness of the fit of $\chi^{2}/ndf=0.7$. We see that the results are well described with the chosen fit. Parameters $\beta$ and $\gamma$ are both consistent with zero, meaning the dominant contribution to the nonflow comes from the six-particle correlator. This is the same leading order behavior as in the corresponding traditional cumulant expansion for this particular SC. Therefore SC($k$,$l$,$m$) is not sensitive to lower order correlations.

This simulation has been done in a totally controlled environment which is not a true representation of what happens in a heavy-ion collision. We check now SC properties in realistic Monte Carlo simulations in the next section, and provide first predictions.

III Predictions from realistic Monte Carlo studies

In this section we discuss the physics of our new observables. We provide studies and obtain separate predictions both for the coordinate and momentum space, by using two realistic Monte Carlo models: iEBE-VISHNU and HIJING.

III.1 Predictions from iEBE-VISHNU

To inspect to what extent the generalized SC can capture the collective behavior of the heavy-ion collision evolution, we use iEBE-VISHNU Shen:2014vra , a heavy-ion collision event generator based on the hydrodynamic calculations. In this event generator, after preparing the initial state using Monte Carlo approach, the evolution of the energy density is calculated via 2+1 causal hydrodynamics together with equation of state calculated from a combination of lattice QCD and hadronic resonance gas model (s95p-v1 Huovinen:2009yb ). After the hydrodynamic evolution, the Cooper-Frye formula Cooper:1974mv is used to convert each fluid patch to a distribution of hadrons. Using iEBE-VISHNU, one can study the hadronic gas evolution after hadronization by using Ultra-relativistic Quantum Molecular Dynamics (UrQMD) transport model Schenke:2012wb . The evolution stops if no collision or decay happen in the system. In the present study, we have not used UrQMD to decrease the simulation time.

In this paper, we study Pb–Pb collision with center-of-mass energy per nucleon pair $\sqrt{s_{\rm NN}}=2.76\,$TeV. For the initial state, MC-Glauber with $0.118$ as wounded nucleon/binary collision mixing ratio is used. For the hydrodynamic evolution, the DNMR formalism Denicol:2010xn ; Denicol:2012cn with fixed shear viscosity over entropy density ($\eta/s=0.08$) and zero bulk viscosity is exploited. The hydrodynamic initial time is set to $0.6\;$fm/$c$. The events are divided into 16 centrality classes between 0–80% with equally sized bins. For each centrality, we have generated 14000 events. Let us point out that in the present study we have taken into account the particles $\pi^{\pm}$, $K^{\pm}$, $p$ and $\bar{p}$ in the final particle distribution, since they are the most abundant particles in the final distribution. In iEBE-VISHNU, the reaction plane angle, i.e. the angle between the impact parameter vector and a reference frame, is fixed to be zero for all events. One notes that here we are dealing with a 2+1 dimensional hydrodynamic calculations in which the boost invariance is considered in the third (longitudinal) direction. For that reason, there is no pseudorapidity dependence in the present simulation.

It is worth mentioning that the aim of the simulation in this paper is not to present a precise justification of the experimental observation. In fact, we would try to demonstrate that in the presence of flow the generalized SC have non-trivial values and that their measurements are feasible for Pb–Pb collisions at LHC energies in terms of required statistics. Nevertheless, to show that our simulations can be considered as at least a qualitative prediction for future experimental observation, we present our Monte Carlo simulation results of few well-studied two-harmonic SC and compare them to the published data from ALICE in Appendix D.

Regarding the $p_{T}$ range, the range $0.28<p_{T}<4$ GeV has been chosen in the present study. We should say that in ALICE $p_{T}$ is in the range $0.2<p_{T}<5$ GeV. It turns out that SC is very sensitive to the lower limit of $p_{T}$ range where we expect to have a considerable amount of particles. For the range $0.28<p_{T}<4$ GeV, which is more close to ALICE $p_{T}$ range, we have a qualitative agreement between simulation and data for SC of two flow amplitudes (see Appendix D). The reason we have used $0.28$ (not $0.3$ or $0.2$) for the lower $p_{T}$ range is that the $p_{T}$ dependent output of VISHNU is reported in a fixed discrete range between $0<p_{T}<4$ GeV. We should point out that our computations with lower $p_{T}=0.28$ GeV show a reasonable agreement with Ref. Zhu:2016puf in which VISHNU simulation with MC-Glauber model for the initial state and constant $\eta/s=0.08$ has been studied.

According to the collective picture of the produced matter in heavy-ion collision, the anisotropic flow corresponds to the anisotropy in the initial state, quantified by eccentricities $\epsilon_{n}e^{in\Phi_{n}}$. In fact, eccentricities are defined as the moments of initial energy density Teaney:2010vd ,

[TABLE]

where $(r,\varphi)$ are the polar coordinates in the two dimensional transverse space, and $\rho(r,\varphi)$ is the energy density in this space. Similar to the anisotropic flow, the SC of the eccentricity distribution shown by $\text{SC}_{\epsilon}(n,m,p)$ can be studied by replacing $v_{n}$ with $\epsilon_{n}$ in Eq. (5).

In addition, the flow fluctuations are the manifestation of initial state fluctuations after the collective evolution. One notes that the difference between initial state models (Glauber, CGC, TrENTo, etc.) is not only encoded in different cumulants of eccentricities $\epsilon_{2}$ or $\epsilon_{3}$ but also in the other details of the initial state fluctuations, for instance in the correlations between eccentricities. In other words, to discriminate between different models, we need to study different aspects of the eccentricity fluctuations and their correlations. The higher order SC observables introduced in this paper should be considered as an example of such a study. Modeling collective evolution as a linear/non-linear response to the initial state helps our understanding of both the initial state and the collective evolution models. Using it, one can explain the experimentally observed flow fluctuations, in the sense that we can qualitatively inspect how much of the observed value is a manifestation of the initial state and how much of it is the consequence of the collective evolution. Even more than that, it is also possible to examine the strength (couplings) of the linear/non-linear terms with this method. A more detailed investigation in this context is postponed to future studies.

Using our Monte Carlo simulations, we have shown SC(2,3,4), SC(2,3,5), SC(3,4,5) and SC(2,4,6) in Fig. 4. Also the same quantities which are obtained from the initial state are presented in Fig. 5. One can see from the Figs. 4 and 5 that except for SC(3,4,5) in the final state (Figs. 4(c)) all the other presented cumulants show significant non-vanishing values. Moreover, one can observe that the cumulants obtained from the initial state (Fig. 5) are monotonically increasing, while the same cumulants for the distribution after hydrodynamic evolution change the slope after centrality classes $45\%$ to $55\%$ and approach to zero. This can be due to the fact that, in the more peripheral collisions, the collective evolution duration is shorter, and small system size cannot transport the initially produced correlation into the final momentum distribution. Moreover, one can see a sign flip in $\text{SC}(2,3,4)$ and a suppression in $\text{SC}(3,4,5)$ in Fig. 4, compared to the initial state generalized SC in Fig. 5. We will return to these points after defining the normalized generalized SC in the following.

As it can be seen from Fig. 4 and Fig. 5, due to the different scale of the initial and final distributions, we are not able to compare $\text{SC}(n,m,p)$ and $\text{SC}_{\epsilon}(n,m,p)$. In order to clarify to what extent the observed values of SC are related to the correlations in the initial state, we will exploit the normalized generalized SC,

[TABLE]

With such a study, for instance, one can answer what is the influence of the collective stage in the heavy-ion evolution. If NSC calculated in terms of $\epsilon_{n}$ and $v_{n}$ are the same, that means that the anisotropies in the initial state are the dominant source of anisotropies in the final state, and therefore NSC can be used to tune the details of initial conditions only. Using normalized (generalized) SC, we also avoid the sensitivity to the $p_{T}$ range in the final particle distribution. In fact, NSC clearly eliminates any dependence of multiplying amplitude of flow harmonics, which was obtained independently from correlators depending only on the same harmonics. Therefore, the independent information contained only in the correlations can be extracted best only from the normalized SC. While this is rather straightforward to achieve in models having only flow correlations, it is much more of a challenge in experimental analyses, due to difficulties in suppressing nonflow contribution in the denominator in Eq. (16). In ALICE:2016kpq , where normalized SC were measured for the first time, the nonflow in denominator was suppressed by introducing pseudorapidity gaps in two-particle correlations.

The generalized NSC are depicted in Fig. 6. The similarity and discrepancy between $\text{NSC}(n,m,p)$ and $\text{NSC}_{\epsilon}(n,m,p)$ can be explained qualitatively by considering the linear and non-linear hydrodynamic response. As a matter of fact, the linear response is approximately true for $\epsilon_{2}$ and $\epsilon_{3}$ in the central and mid-central collisions Gardim:2011xv . It means the event with larger ellipticity $\epsilon_{2}$ has larger elliptic flow $v_{2}$, and the dependence is approximately linear. The same relation holds between triangularity $\epsilon_{3}$ and triangular flow $v_{3}$. However, this is not the case for higher harmonics. For instance, it has been shown that for $v_{4}e^{4i\Psi_{4}}$ and $v_{5}e^{5i\Psi_{5}}$ we have the following relations Gardim:2011xv ; Teaney:2012ke ; Teaney:2013dta ; Giacalone:2018syb ,

[TABLE]

where $k_{n}$ and $k^{\prime}_{n}$ are coefficients related to the hydrodynamic response.

Interestingly, compared to $\text{NSC}_{\epsilon}(2,3,4)$ a sign change can be seen in NSC(2,3,4) (Fig. 6(a)) which is similar to what has been observed for NSC(3,4) and $\text{NSC}_{\epsilon}(3,4)$ in Fig. 13(c) in Appendix D. In the present case, we are dealing with genuine three-harmonic observable. However, the main difference between generalized SC of the initial state and the final state comes from the contribution of the non-linear term $\epsilon_{2}^{2}$ in the $v_{4}$. In fact, the term $\epsilon_{2}^{2}$ and its anti-correlation with $\epsilon_{3}$ should be responsible for this sign change similar to the NSC(3,4) and $\text{NSC}_{\epsilon}(3,4)$ case (see Appendix D). The same logic can explain the suppression of the cumulant $\text{NSC}(3,4,5)$ in Fig. 6(c) too. We know that there is a non-linear contribution with the term $\epsilon_{2}\epsilon_{3}$ in $v_{5}$ (see Eq. (17)). As a result, the terms $\epsilon_{2}^{2}$ and $\epsilon_{2}\epsilon_{3}$ can explain the small value of SC(3,4,5) (or suppression of $\text{NSC}(3,4,5)$ in comparison with $\text{NSC}_{\epsilon}(3,4,5)$). However, in $\text{NSC}(2,3,5)$ only the term $\epsilon_{2}\epsilon_{3}$ plays the role. As can be seen in Fig. 6(b), the effect of the term $\epsilon_{2}\epsilon_{3}$ is small in $\text{NSC}(2,3,5)$. It is worth to note that in even simpler cumulants $\text{NSC}(2,3)$ and $\text{NSC}_{\epsilon}(2,3)$ (see Fig. 13(b) in Appendix D) we do not have such an agreement in centrality classes above $40\%$. We think there must be a more complex reason such as the presence of extra non-linear terms behind the approximate agreement between $\text{NSC}(2,3,5)$ and $\text{NSC}_{\epsilon}(2,3,5)$ in a wide range of centrality classes. Finally, compared to $\text{NSC}_{\epsilon}(2,4,6)$, we observe a considerable enhancement in $\text{NSC}(2,4,6)$ in Fig. 6(d). This enhancement is even larger than what has been observed in Fig. 13(b) for SC(2,4). The reason would be due to the fact that in $\text{NSC}(2,4,6)$ the terms $\epsilon_{2}^{2}$ in $v_{4}$ and the $\epsilon_{2}^{3}$ in $v_{6}$ are responsible for this enhancement. It seems the term $\epsilon_{3}^{2}$ in $v_{6}$ and its anti-correlation with $\epsilon_{2}$ does not have enough power to compete with $\epsilon_{2}$, $\epsilon_{2}^{2}$ and $\epsilon_{2}^{3}$ trivial correlations. It is important to point out that we have explained the generalized SC from the initial state and hydrodynamic response only qualitatively. In this context, a further rigorous study is needed to be done in the future.

III.2 Estimating nonflow contribution with HIJING

As detailed in Appendix A, one of requirements for our observables to be called SC is that it should be robust against nonflow. We have already studied in a Toy Monte Carlo simulation the case of strong two-particle correlation for different multiplicities in Sec. II.5. We now investigate nonflow contribution further and introduce HIJING, which stands for Heavy-Ion Jet INteraction Generator Gyulassy:1994ew . It is a Monte Carlo model used to study particle and jets production in nuclear collisions. It contains models to describe mechanisms like jet production, jet fragmentation, nuclear shadowing, etc. The correlations these mechanisms introduce involve generally only few particles and should not be included in the analysis of collective effects like anisotropic flow. Since HIJING has all the phenomena produced in a heavy-ion collision except flow itself, we can use it to test the robustness of our SC observable against few-particle nonflow correlations.

In general, when one uses an azimuthal correlator in an expression like

[TABLE]

one assumes in the derivation that the underlying two-variate p.d.f. $f(\varphi_{1},\varphi_{2})$ fully factorizes into the product of two marginal p.d.f.’s, i.e.:

[TABLE]

If all azimuthal correlations are induced solely by anisotropic flow, such factorization is exactly satisfied. In the case of few-particle nonflow correlations the above factorization is broken, and flow amplitudes $v_{n}$ are now systematically biased when estimated with azimuthal correlators, i.e. at leading order we expect to have:

[TABLE]

where $\delta_{2}$ denotes the systematic bias in 2-particle azimuthal correlations due to nonflow. We have written the above argument for the two-particle azimuthal correlations, but the argument can be trivially generalized for any larger number of particles (see the discussion on nonflow in Appendix A). It is impossible in practice to calculate and to quantify the systematic bias $\delta_{2}$ stemming from vastly different sources of nonflow correlations. Since HIJING has all relevant sources of correlations, except anisotropic flow, it is a great model to estimate solely the value of this systematic bias $\delta_{2}$, i.e. in HIJING we expect the following relation to hold:

[TABLE]

and analogously for higher order azimuthal correlators. If the given multi-particle azimuthal correlator, or some compound observable like SC which is estimated from azimuthal correlators, when measured in the HIJING dataset is consistent with 0, that means that nonflow contribution in that observable is negligible, and therefore it can be used as a reliable estimator of anisotropic flow properties.

In this section we show the predictions from HIJING for the centrality dependence of two different combinations of harmonics: SC(2,3,4) and SC(2,3,5). The data used here correspond to Pb–Pb collisions taken at the center of mass energy of $\sqrt{s_{\text{NN}}}$ = 2.76 TeV. Two kinetic cuts have been applied: 0.2 $<p_{T}<$ 5 GeV and -0.8 $<\eta<$ 0.8. The results obtained with HIJING are shown on Fig. 7 (a) for SC(2,3,4) and on Fig. 7 (b) for SC(2,3,5), alongside with the VISHNU predictions for the same combinations of harmonics. We can see that in both cases, our new SC observable in HIJING is compatible with 0 for head-on and mid-central collisions, meaning our observable is robust against nonflow. The comparison also shows that flow has generally a more important contribution than nonflow. This means that any observation of nonzero values of SC in real data could be attributed to collective effects.

IV Summary

In summary, we have presented the generalization of recently introduced flow observables, dubbed Symmetric Cumulants. We have presented the unique way how the genuine multi-harmonic correlations can be estimated with multi-particle azimuthal correlators. All desired properties of higher order Symmetric Cumulants were tested with the carefully designed Toy Monte Carlo studies. By using the realistic iEBE-VISHNU model, we have demonstrated that their measurements are feasible and we have provided the first predictions for their centrality dependence in Pb–Pb collisions at LHC energies. A separate study has been presented for their values in the coordinate space. Study based on HIJING demonstrated that these new observables are robust against systematic biases due to nonflow correlations. These generalized, higher order observables, can be therefore used to extract reliably new and independent information about the initial conditions and the properties of QGP in high-energy nuclear collisions.

Acknowledgements.

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 759257). This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe.”

Appendix A List of requirements for the generalization of Symmetric Cumulants

In this self-contained appendix we establish how SC need to be generalized for the case of correlated fluctuations involving more than two flow amplitudes. Out of plenty of ways to perform such a generalization, we narrow down the one which uniquely satisfies the following list of requirements, which must hold for any number of flow amplitudes:

No built-in trivial contribution. For constant flow amplitudes $v_{k}$, $v_{l}$, $v_{m}$, …, the observable SC($k,l,m,\ldots$) is identically zero. 2. 2.

Genuine multi-harmonic correlations (cumulants). If fluctuations only of a subset of amplitudes $v_{k}$, $v_{l}$, $v_{m}$, …, are correlated, SC($k,l,m,\ldots$) is identically zero. This requirement ensures that the SC($k,l,m,\ldots$) observable is a well defined cumulant, i.e. it isolates only the genuine correlation among all flow amplitudes in question. This cumulant property secures that the SC extracts unique and independent information at each order: SC constructed from $n$ distinct flow amplitudes provides information which cannot be accessed from SC constructed from smaller number of $n-1$, $n-2$, …, flow amplitudes. This requirement also ensures in particular that if fluctuations of amplitudes $v_{k}$, $v_{l}$, $v_{m}$, …, are all mutually uncorrelated, SC($k,l,m,\ldots$) is identically zero. This requirement also ensures that all amplitudes in the SC definition must be different (e.g. if one starts with the definition of two-harmonic SC$(k,l)\equiv\left<v_{k}^{2}v_{l}^{2}\right>-\left<v_{k}^{2}\right>\left<v_{l}^{2}\right>$ and identifies in the resulting expression $l=k$, then SC($k,l$) is trivially non-zero as it reduces to the variance $\left<v_{k}^{4}\right>-\left<v_{k}^{2}\right>^{2}$ of one flow amplitude). For other formal mathematical properties which any cumulants must satisfy we refer to Kubo . 3. 3.

Symmetry. Observable SC$(k,l,m,\ldots)$ must be the same by definition for any permutation of amplitudes $v_{k},v_{l},v_{m},\ldots$. Example: $\mathrm{SC}(k,l)=\mathrm{SC}(l,k)$, $\mathrm{SC}(k,l,m)=\mathrm{SC}(k,m,l)=\mathrm{SC}(l,k,m)=\mathrm{SC}(l,m,k)=\mathrm{SC}(m,l,k)=\mathrm{SC}(m,k,l)$, etc. This requirement distinguishes SC from other possibilities of building cumulants from flow amplitudes. As a consequence of this requirement, all amplitudes in the definition of SC must be raised to the same power. 4. 4.

Cleanliness. Any dependence on the symmetry planes must be canceled out by definition in SC. This is an important requirement, since in general in Fourier series parameterizations of anisotropic distributions (see Eq. (1)), both flow amplitudes $v_{n}$ and symmetry planes $\Psi_{n}$ appear at equal footing, since they are two independent degrees of freedom to quantify anisotropic flow phenomenon. This requirement restricts tremendously the number of azimuthal correlators which are suitable to build the SC. There are two cases of interest in this context which need to be clarified:

1. (a)

In all azimuthal correlators from which SC are built, one always have equal number of positive and negative harmonics (see e.g. Eqs. (7), (33) or (36-40)). Therefore, if there are $k$ azimuthal angles in the multi-particle correlator, we can always label the corresponding harmonics as $n_{1},n_{2},\ldots,n_{k/2}$ and $-n_{1},-n_{2},\ldots,-n_{k/2}$. It follows directly from the general result in Eq. (2):

[TABLE]

In the 1st line above we have used the fact that for any Fourier harmonic $n$ we have $v_{n}=v_{-n}$ and $\Psi_{n}=\Psi_{-n}$, which can be derived straight from the formal properties of Fourier series. Therefore, for the case when symmetry planes are all different and do not fluctuate within the same event, by definition SC observables do not have any contribution from symmetry planes; 2. (b)

Experimentally, if the symmetry planes $\Psi_{n}$ change with $\eta$ or $p_{T}$, which is known as decorrelation or factorization breaking Jia:2014vja ; Pang:2014pxa ; Khachatryan:2015oea ; Bozek:2017qir , the differential analysis of SC observables might still exhibit dependence on symmetry planes. In order to avoid this problem, we advocate the differential SC analyses in the narrow phase-space windows of interest, e.g. narrow $\Delta\eta$ bins, when all particles from the multiplet must be taken from the same $\Delta\eta$ bin. 5. 5.

Isotropy. All azimuthal correlators used in estimating SC experimentally must be isotropic since these are the only correlators which for a detector with uniform acceptance give non-trivial contributions, after the single-event averages of azimuthal correlators have been extended to all-event averages. Isotropic azimuthal correlators are the correlators for which the harmonics satisfy the relation: $\sum_{i}n_{i}$ = 0 Bhalerao:2011yg . Such correlators are invariant under the arbitrary rotation of the coordinate system in which azimuthal angles are being measured. 6. 6.

Uniqueness. No other combination of multi-particle azimuthal correlators shall yield to the same final expression for flow amplitudes. For instance, out of two possibilities to estimate SC$(m,n)$:

[TABLE]

only the definition in Eq. (23) leads to the desired properties of SC. While it is clear, by using Eq. (2), that definition in Eq. (23) evaluates to the desired expression $\left<v_{m}^{2}v_{n}^{2}\right>-\left<v_{m}^{2}\right>\left<v_{n}^{2}\right>$, it is not so obvious that the alternative definition in Eq. (24) evaluates to something else. We clarify this technical point further in Appendix C. This property implies that SC must be estimated experimentally in terms of all-event averages of each individual azimuthal correlator. 7. 7.

Robustness against nonflow. SC must be consistent with zero for a system containing only a few-particle nonflow correlations. Mathematically, nonflow correlations are all correlations for which the factorization of joint multi-variate p.d.f. of $k$-azimuthal angles cannot be written as a product of single-particle p.d.f.’s., i.e. for nonflow correlations:

[TABLE]

If the particle emission is dominated by anisotropic flow, in the above expression we have equality. For $k$-particle azimuthal correlator, the nonflow contribution $\delta_{k}$ exhibits the universal generic scaling as a function of multiplicity $M$:

[TABLE]

By using the above result, it follows that in general one obtains the generic scaling also for the nonflow contribution $\delta^{\rm SC}_{k}$ in $k$-harmonic SC observable:

[TABLE]

The leading order nonflow contribution of $k$-harmonic SC observable is the same as for $2k$-azimuthal correlator. For instance, the scaling of nonflow contribution to SC($k,l,m$) can be described well for any choice of harmonics $k$, $l$ and $m$ with the following relation:

[TABLE]

where constants $\alpha$, $\beta$ and $\gamma$ quantify various sources of nonflow contribution and can be obtained from the fit of multiplicity dependence of corresponding SC observable. Any violation from the few-particle nonflow scaling in Eq. (27) indicates that SC is dominated by contribution from collective correlations, and robust against nonflow. In practice, the nonflow contribution can be estimated by evaluating the SC observable over the data simulated with HIJING event generator (see Sec. III.2), which contains only nonflow correlations. 8. 8.

Event weights. When estimating SC, it is essential that all azimuthal correlators are expressed as all-event averages $\left<\left<\cdots\right>\right>$. This raises a non-trivial question of what is the correct weight to use to weight the event-by-event averages $\left<\cdots\right>$ of azimuthal correlators. Namely, in the following generic expression:

[TABLE]

where the two sums go over all events, there are various possibilities to choose event weights $w_{i}$ from. Different choices will in general yield both to the different statistical and systematical properties of the final SC results. In the present work, we have only tested that the standard choice, namely ‘number of combinations’ event weight Bilandzic:2010jr , yields to the smallest statistical spread of SC (see Appendix C). There is no performance penalty due to double counting (e.g. by treating $\cos[n(\varphi_{1}\!-\!\varphi_{2})]$ and $\cos[n(\varphi_{2}\!-\!\varphi_{1})]$ separately, despite the fact that cosine is an even function) when using this weight, since all expressions for multi-particle azimuthal correlations (which by definition include double counting) are analytically expressed in terms of $Q$-vectors Bilandzic:2013kga . For subevent cumulants, the optimal weight is ‘number of unique combinations’. For the specific case when the event weights depends only on multiplicity, one can completely circumvent this problem in the following two independent ways:

1. (a)

One bins all available data in bins of multiplicity, where the bin width is 1. The final SC results in each multiplicity bin are independent and therefore can be trivially combined with the standard $1/\sigma^{2}$ weighting, where $\sigma$ is the spread of SC results in a given unit multiplicity bin; 2. (b)

For a given subsample of the whole dataset (e.g. events corresponding to the specific centrality), one randomly selects for the analysis in each event the same number of particles (typically corresponding to the smallest multiplicity in that subsample).

In both cases a) and b), the weight $w_{i}$ in Eq. (29) is a constant function of multiplicity by construction, and any constant weight is equivalent to the unit weight. We leave the further study of event weights for our future work. This non-trivial problem was addressed recently also in Gardim:2016nrr .

As an example, we now demonstrate how the Requirements 1–3 are satisfied for the SC$(k,l,m)$ defined in Eq. (6). It is trivial to verify that Eq. (6) is identically zero if all three amplitudes are constant (Requirement 1). When all three amplitudes fluctuate independently, all averages in Eq. (6) can be factorized completely and it follows:

[TABLE]

When all three amplitudes fluctuate, but only two fluctuations are correlated (e.g. among harmonics $k$ and $l$), one can only partially factorize some of the averages in Eq. (6), but the1 resulting observable is nevertheless still identically zero:

[TABLE]

The same conclusion follows if only correlations among the other two amplitudes, namely $v_{k}$ and $v_{m}$, or $v_{m}$ and $v_{l}$, are correlated. The results in Eqs. (30) and (31) demonstrate that $\mathrm{SC}(k,l,m)$ observable is well defined cumulant of flow amplitudes (Requirement 2). Trivially, $\mathrm{SC}(k,l,m)$ is symmetric in all three harmonics (Requirement 3). The remaining requirements are discussed in detail in Appendix C with Monte Carlo studies.

Appendix B Definitions for higher order Symmetric Cumulants

In this appendix, we write down explicitly the equations of few higher order SC, which satisfy all requirements set upon their generalization in Appendix A. We provide first the theoretical definitions in terms of flow amplitudes, and then the corresponding relations in terms of azimuthal correlators, which can be used to estimate SC experimentally.

The four-harmonic SC$(k,l,m,n)$ we define directly as:

[TABLE]

Experimentally, SC$(k,l,m,n)$ can be estimated with the following combination of azimuthal correlators:

[TABLE]

The five-harmonic SC$(k,l,m,n,o)$ we define directly as:

[TABLE]

Experimentally, SC$(k,l,m,n,o)$ can be estimated with the following expression:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Analogously, one can proceed with all higher order SC.

Appendix C Checking the other requirements for the Symmetric Cumulants

In this appendix, we develop the argumentation about the remaining requirements outlined in Appendix A. In what follows, all simulations are done using the Toy Monte Carlo setup described in Sec. II.5.

We first look at the effects of the correlations between less than three flow amplitudes. We consider $N=10^{8}$ events in total and set the values of $v_{2}$, $v_{3}$ and $v_{4}$ according to the current test for five different multiplicities: 100, 250, 500, 750 and 1000. As we want to study the multiplicity fluctuations and the event weights separately, we keep constant the multiplicity. This implies that any event weight is equal to unit weight. The configurations of flow amplitudes used in our basic tests are the following: $v_{2}$, $v_{3}$ and $v_{4}$ are set to zero (Fig. 8 (a)), constant nonzero amplitudes fixed to $v_{2}=0.15$, $v_{3}=0.13$, $v_{4}=0.1$, respectively (Fig. 8 (b)), uncorrelated amplitudes uniformly sampled event-by-event in the interval (0.05, 0.09) (Fig. 8 (c)) and finally $v_{2}$ and $v_{3}$ uniformly sampled in (0.05, 0.09) and $v_{4}$ correlated to $v_{2}$ with the relation $v_{4}=v_{2}-(0.005,0.025)$ (Fig. 8 (d)). All four studies lead to the same conclusion that SC(2,3,4) is compatible with zero for all multiplicities. This means our definition in Eq. (7) does not contain any built-in bias in absence of fluctuations in the amplitudes, but also that it is insensitive to any subset of amplitudes where the fluctuations of at least one of the $v_{n}$ are not correlated with the others.

We now look at an example of genuine correlations among the three flow amplitudes. We study the case where $v_{2}$ is uniformly sampled in (0.03, 0.1), $v_{3}=v_{2}-(0,0.02)$ and $v_{4}=v_{2}-(0.005,0.025)$. These values are input in our Toy Monte Carlo setup for $M$ = 100, 250, 500, 750 and 1000, and for $N=10^{8}$. The theoretical value of SC(2,3,4) is $5.13992\cdot 10^{-9}$. We can see on Fig. 9(a) that for all considered values of multiplicities, SC(2,3,4) obtained through simulations is compatible with the theoretical result within the statistical errors. This result, combined with the ones obtained in the previous tests (Fig. 8), implies that our definition of SC given by Eq. (7) is sensitive only to the genuine correlation among all three flow amplitudes in question.

Now that we have demonstrated that our SC is sensitive only to three genuinely correlated flow amplitudes, we provide an argumentation about the uniqueness of its expression (Requirement 6 in Appendix A). We first show which expression must be chosen to solve the ambiguity between the two definitions in the case of SC($m$,$n$) using a mathematical derivation. As the exact computations are tedious for three different amplitudes, we then look at results for SC($k$,$l$,$m$) obtained with our Toy Monte Carlo setup and with HIJING.

We start by recalling here the definitions of the two-, three- and four-particle azimuthal correlators for a number of particle per event $M$ and for unit particle weights Bilandzic:2010jr :

[TABLE]

where we have the condition that $m+n+o=0$ for $\langle 3\rangle_{m,n,o}$. The computation of the single-event average of the first term of Eq. (24), $\langle 2\rangle_{m,-m}\langle 2\rangle_{n,-n}$, gives us:

[TABLE]

For $M\gg 3$, we can write that $M-k\simeq M$ for $k=1,2,3$. This is the case in heavy-ion collisions where a large amount of particles are created in each event. After averaging over all the events, we obtain finally that

[TABLE]

We can therefore conclude that in the definitions in Eqs. (23) and (24) the first terms are not equal. The reason is that in the alternative expression in Eq. (24) the first term contains self-correlations between the particles. Such self-correlations must be removed exactly in correlations techniques as they trivially and non-negligibly bias the final results.

In order to look at the ambiguity in the correct definition of SC($k$,$l$,$m$), we introduce two other expressions leading to the same final theoretical result given by Eq. (7). These are:

[TABLE]

In what follows, we will refer to Eq. (7) as the usual expression of SC($k,l,m$), to Eq. (46) as its alternative expression and to Eq. (47) as its mixed expression. Like for the case of SC($m$,$n$) in Eq. (24), the alternative definition contains only two-particle correlators for every harmonic. The mixed formula is called this way as it is a combination of two- and four-particle correlators. Since the exact derivation made for two amplitudes becomes too tedious in the case for three differents amplitudes, we will try to establish which one is the unique expression to use for SC($k$,$l$,$m$) both by using the Toy Monte Carlo setup described in Sec. II.5 and by using HIJING.

We first look at the Toy Monte Carlo simulation. We sample uniformly event-by-event, for $N=10^{7}$ events, $v_{2}$ and $v_{4}$ in the interval (0.05, 0.09) and set that $v_{3}=v_{2}+0.02$, before computing the three expressions described above. Like in our basic tests, the multiplicity is kept constant. Since not all three amplitudes are correlated, we expect that SC(2,3,4) is zero. The results of the simulations are shown on Fig. 9(b). The usual expression is zero for all considered multiplicities. However, this is not the case for the alternative and the mixed definitions which show a clear multiplicity dependence. Similarly as in the expression presented above, one can expect the terms of the form $\langle\langle\cos(a\varphi_{1}+b\varphi_{2}-a\varphi_{3}-b\varphi_{4})\rangle\langle\cos(c(\varphi_{1}-\varphi_{2}))\rangle\rangle$ for any $a,b,c=k,l,m$ and $\langle\langle\cos(k(\varphi_{1}-\varphi_{2}))\rangle\langle\cos(l(\varphi_{1}-\varphi_{2}))\rangle\langle\cos(m(\varphi_{1}-\varphi_{2}))\rangle\rangle$ to introduce self-correlations which are not visible in the definition of SC($k$,$l$,$m$) with flow amplitudes.

We now look at the results obtained with HIJING (the description of the HIJING setup and resulting dataset can be found in Sec. (III.2)). Figure 10 shows the centrality dependence of SC(2,3,4) and SC(2,3,5) in the case of the usual expression in Eq. (7), the alternative expression in Eq. (46) and the mixed expression in Eq. (47). Both alternative and mixed expressions show large statistical errors, which prevent us to rule out the compatibility of these results with zero. One possible reason of the size of these errors can be that the multiplicity dependence of the terms with the self-correlations is also present for their statistical errors. However, the usual expression has smaller statistical errors and, therefore, seems to be more consistent to zero for all centrality bins and less sensitive to nonflow.

Finally, we conclude this appendix with the discussion on event weights (Requirement 8 in Appendix A). We therefore simulate $N=10^{8}$ events with $v_{2}=0.09$, $v_{3}=0.07$ and $v_{4}=0.05$ with the setup described in Sec. II.5. Since all event weights are equal to unity if the multiplicity is kept constant, we sample $M$ uniformly in the interval (50, 500). We consider three different event weights. The first one is the weights used in ALICE:2016kpq and which was named the number of combinations. In our case, this weight is equal to $M(M-1)$ for the two-particle correlators, $M(M-1)(M-2)(M-3)$ for the four-particle correlators and finally $M(M-1)(M-2)(M-3)(M-4)(M-5)$ for the six-particle correlator. The other two event weights we introduce are the multiplicity of the event itself and the unit weight. On Fig. 11 we see that the number of combinations weight has the smallest statistical spread. This result is consistent with what has been done in the analyses of individual flow harmonics with correlation techniques. This implies that the number of combinations is still the event weight to use in the case of the three-harmonic SC observable.

Appendix D Vishnu predictions for SC($m,n$)

In this appendix, we reproduce few well-studied two-harmonic SC with our Monte Carlo simulations. In Fig. 12, we have shown centrality dependence of $\text{SC}(2,3)$ and $\text{SC}(2,4)$ obtained from iEBE-VISHNU in three different $p_{T}$ ranges, and we have compared it with the experimental results from ALICE at center-of-mass energy per nucleon pair $\sqrt{s_{\rm NN}}=2.76$ TeV ALICE:2016kpq . The ALICE $p_{T}$ is in the range $0.2<p_{T}<5$ GeV, and the range $0.28<p_{T}<4$ GeV in our simulation shows reasonable agreement with the ALICE data. Also we have shown few examples of normalized SC in Fig. 13(b-d). One can see in Fig. 13(b) that $\text{NSC}(2,3)$ and $\text{NSC}_{\epsilon}(2,3)$ are in good agreement in the centrality classes below $40\%$. It is due to the fact that the linear response is approximately true for $\epsilon_{2}$ and $\epsilon_{3}$ in this range Gardim:2011xv . As matter of fact, the large discrepancy between $\text{NSC}(2,4)$ and $\text{NSC}_{\epsilon}(2,4)$ (Fig. 13(c)) can be explained by non-linear term in $v_{4}$ equation. This effect has more non-trivial impact in the case of $\text{NSC}(3,4)$ and $\text{NSC}_{\epsilon}(3,4)$ (Fig. 13(d)) where the normalized SC of the initial state shows a sign flip after the non-linear hydrodynamic response Zhu:2016puf . In this case, the $\epsilon_{2}$ contribution inside $v_{4}$ is anti-correlated with $\epsilon_{3}$. Therefore the non-linear term leads to the suppression of $\text{NSC}(3,4)$ compared to $\text{NSC}_{\epsilon}(3,4)$ (see Ref. Zhu:2016puf ). This is fascinating because the same feature has been observed for the generalized SC in our study.

Appendix E Ellipse-like distributions

In this appendix, we demonstrate mathematically how the higher order SC observables can provide further and independent constraints. For clarity sake, we consider a simple Toy Monte Carlo model of ellipse-like distributions. Despite its simplicity, and since in mid-central heavy-ion collisions due to collision geometry the overlapping volume containing the strongly interacting nuclear matter is to leading order ellipsoidal, this toy model has some relevance also for the real-case scenarios. When the polar coordinates are measured from the center of ellipse, the ellipse is determined by the following relation:

[TABLE]

where $\varepsilon$ is the eccentricity defined as

[TABLE]

and $a(b)$ is semimajor(semiminor) axis. Since we take $b\leq a$, it follows from (49) that eccentricity satisfies

[TABLE]

For the distribution $f(\varphi)$ in Eq. (48) one can calculate the Fourier harmonics $v_{n}$ by using the following standard expression (due to symmetry of the problem all symmetry planes $\Psi_{n}$ are the same, and we can always span the coordinate system in such a way that $\Psi_{n}$ coincides with the orientation of semiminor axis $b$, which means that $\Psi_{n}=0,\ \forall n$, in this simple model):

[TABLE]

The integral in (51) cannot be expressed in terms of elementary functions. For instance, for the overall normalization $v_{0}$ and for the harmonic $v_{2}$ we have obtained the following results:

[TABLE]

where ${}_{2}F_{1}\left(a,b;c;z\right)$ is the hypergeometric function whose series expansion reads

[TABLE]

We now inspect the leading order behavior of these expressions. Since the eccentricity satisfies relation (50), we can expand the results for harmonics obtained from (51) as a series in $\varepsilon$. After straightforward calculation, we have obtained the following results for the first non-vanishing harmonics,

[TABLE]

We stress out that all odd harmonics vanish for ellipse-like distributions when the polar coordinates are measured from the center due to the symmetry with respect to $y$-axis.

Using all above results, with the straightforward calculus we have obtained for the normalized ellipse-like distribution (normalization amounts simply to dividing with the overall normalization $v_{0}$ the starting expression in Eq. (48) and all final results for harmonics $v_{n}$ in Eq. (55)) the following relations among the even flow amplitudes $v_{2},v_{4},v_{6},\ldots$

[TABLE]

We conclude that when the distribution of particles in heavy-ion collisions is ellipsoidal, there is always a geometric correlation among all even harmonics, determined by relations in Eqs. (56). If the fluctuations occur only in the magnitude of elliptical shape (e.g. only $\varepsilon$ and $b$ in Eq. (48) fluctuate event-by-event), the shape remains in each event elliptical, and therefore odd harmonics cannot develop Kolb:2003zi . This sort of fluctuations will contribute only to SC(2,4), SC(2,4,6), etc., but not to SC(2,3,4). The non-vanishing result for SC(2,3,4) implies that there are additional sources of fluctuations in the system which couple all three amplitudes $v_{2},v_{3},v_{4}$. In particular, SC(2,3,4) can separate the fluctuations in the shape of ellipsoidal from the magnitude fluctuations of the persistent ellipsoidal shape. On the other hand, SC(2,4) cannot disentangle these two different sources of fluctuations.

Appendix F Comment on Parseval theorem

One may wonder to what extent correlated fluctuations of different flow harmonics might originate from some non-physical built-in mathematical property, solely from definitions. To clarify this point we start with the standard parametrization of Fourier series in Eq. (1). We can attach to the function $f(\varphi)$ a physical interpretation in the context of flow analyses only if we also promote $f(\varphi)$ to the level of p.d.f. This is achieved by imposing the following two constraints:

[TABLE]

From the purely mathematical point of view, flow amplitudes $v_{n}$ in Eq. (1) always satisfy the following relation, so-called Parseval theorem:

[TABLE]

In the combination with the probability constraint in Eq. (57), this relation could then lead always to some trivial built-in correlation among flow amplitudes. For instance, if anisotropic flow is quantified only with two amplitudes $v_{m}$ and $v_{n}$ and if the RHS in Eq. (58) would be constrained to some constant value, we would then have:

[TABLE]

i.e. amplitudes $v_{m}$ and $v_{n}$ would be always trivially anticorrelated, due to Parseval identity (on the other hand, correlated amplitudes can be simply generated with the different relative signature in an analogous expression, i.e. $v_{m}^{2}-v_{n}^{2}=\mathrm{const}$). However, that is not the case, since the knowledge of integral of a function has a priori no relation with the knowledge of integral of that function squared. That means that even though in each high-energy nuclear collision the condition $\int_{0}^{2\pi}f(\varphi)\,d\varphi=1$ is strictly satisfied, $\int_{0}^{2\pi}|f(\varphi)|^{2}\,d\varphi$ can nevertheless take any value. From another perspective, even if Parseval identity in Eq. (58) would enforce some correlation among flow amplitudes $v_{m}^{2}$ and $v_{n}^{2}$, it would not be able to explain why then some amplitudes would be correlated, and some anticorrelated, as was observed experimentally for $v_{2}^{2}$ and $v_{4}^{2}$, and $v_{2}^{2}$ and $v_{3}^{2}$, respectively ALICE:2016kpq . We conclude that there is no trivial built-in correlation among different flow amplitudes which is originating solely from the fact that flow amplitudes are the degrees of freedom in the Fourier series expansion, to which also the probabilistic interpretation was attached via Eq. (57).