Measurements of mixed harmonic cumulants in Pb-Pb collisions at $\mathbf{\sqrt{{\textit s}_{\rm NN}}}=5.02$ TeV

Correlations between moments of different flow coefficients are measured in Pb-Pb collisions at $\sqrt{s_{\rm NN}} = 5.02$ TeV recorded with the ALICE detector. These new measurements are based on multiparticle mixed harmonic cumulants calculated using charged particles in the pseudorapidity region $|\eta|<0.8$ with the transverse momentum range $0.2<p_{\rm T}<5.0$ GeV/$c$. The centrality dependence of correlations between two flow coefficients as well as the correlations between three flow coefficients, both in terms of their second moments, are shown. In addition, a collection of mixed harmonic cumulants involving higher moments of $v_2$ and $v_3$ is measured for the first time, where the characteristic signature of negative, positive and negative signs of four-, six- and eight-particle cumulants are observed, respectively. The measurements are compared to the hydrodynamic calculations using iEBE-VISHNU with AMPT and TRENTo initial conditions. It is shown that the measurements carried out using the LHC Run 2 data in 2015 have the precision to explore the details of initial-state fluctuations and probe the nonlinear hydrodynamic response of $v_2$ and $v_3$ to their corresponding initial anisotropy coefficients $\varepsilon_2$ and $\varepsilon_3$. These new studies on correlations between three flow coefficients as well as correlations between higher moments of two different flow coefficients will pave the way to tighten constraints on initial-state models and help to extract precise information on the dynamic evolution of the hot and dense matter created in heavy-ion collisions at the LHC.


Introduction
One of the fundamental questions in the phenomenology of quantum chromodynamics is what are the properties of matter at extreme densities and temperatures where quarks and gluons are in a state of matter called the quark-gluon plasma (QGP) [1,2]. High-energy heavy-ion collisions at the Relativistic Heavy Ion Collider (RHIC) at BNL and the Large Hadron Collider (LHC) at CERN create such a state of strongly interacting matter allowing us to study its properties in the laboratory. Anisotropic flow is a key phenomenon that provides important information about the transport properties of the created QGP matter. Due to large pressure gradients, the anisotropy of the overlapping region between two colliding nuclei causes an anisotropic distribution of the emitted particles in the final state. This anisotropic particle distribution can be quantified by anisotropic flow [3,4] which is characterized by the single-particle azimuthal distribution, Here ϕ is the azimuthal angle of the emitted particle, v n and Ψ n are the n-th order flow coefficient and flow symmetry plane, respectively. Both v n and Ψ n define the n-th order flow-vector as − → V n = v n e inΨ n . The size and direction of − → V n related to the initial anisotropy − → ε n vector is defined by the moments of the shape of the transverse positions (r, φ ) of the participating nucleons, − → ε n = ε n e −inΦ n = − r n e −inφ r n , (n > 1) (2) where ε n and Φ n are the magnitude and orientation of − → ε n , respectively, and stands for the average over all participating nucleons in the initial state. For lower orders, n = 2 and 3, originally a linear response of v n to ε n was expected, with v n = κ n ε n [5,6] where κ n is a parameter that encodes the transport properties of the produced QGP. Later on, it was noticed in models that, already in semi-peripheral collisions, the correlation between the initial ε 2 (ε 3 ) and the final-state v 2 (v 3 ) is not completely linear, with a nonnegligible spread in the correlation between v n and ε n [7]. Such a nonlinear response of lower-order v n should be related to the dynamic evolution of the system, but it was briefly investigated in previous studies [7][8][9]. For the higher orders, n ≥ 4, − → V n receives a significant nonlinear hydrodynamic response from − → ε 2,3 in non-central collisions, which was studied in great detail [10][11][12][13][14][15][16][17].
One can describe the distribution of final-state anisotropies using a joint probability density function (p.d. f .) in terms of v n and Ψ n as P(v m , v n , ..., Ψ m , Ψ n , ...). This is sensitive to the spatial anisotropy ε n , its event-by-event fluctuations, the correlations between different orders of anisotropy coefficients and initial participant planes Φ n carried by P(ε m , ε n , ..., Φ m , Φ n , ...) and it also reflects the early state dynamics and the transport properties of the QGP. Although ideally one would like to measure P(v m , v n , ..., Ψ m , Ψ n , ...), this is not straightforward to achieve in experiments, but what can be measured are the projections of the full p.d. f . on a finite number of variables [18]. Most of these projected distributions could be classified into the following types: (1) v n fluctuations P(v n ) for both integrated and differential v n measurements, The v n coefficients were measured up to the ninth order with an unprecedented degree of precision [16]. The full p.d. f . of single v n coefficients P(v n ) was either measured with a Bayesian unfolding procedure [19,20] (for n = 2, 3 and 4) or constructed via the measured moments (for n = 2) [21]. It was found that the P(v n ) distribution, which originates from the p.d. f . of initial-state ε n distribution P(ε n ), is described better by an elliptic-power function than a Bessel-Gaussian function [21]. It was also realized that during the expansion the produced particles might not share a common flow symmetry plane at different transverse momenta, p T , and pseudorapidity, η [22,23]. These transverse momentum and pseudorapidity dependent flow vectors fluctuate event-by-event, which also breaks the factorization of two-particle correlations V (p t T , p a T ) into the product of flow coefficients v n (p t T ) · v n (p a T ) [24][25][26]. Such phenomena were predicted by hydrodynamic calculations and are found to be sensitive to the initialstate density fluctuations and/or to the specific shear viscosity of the expanding medium [22,23,27]. In addition, analyses of correlations between different order flow vectors [13,28,29] show promise to shed additional light on the initial-state conditions. The correlations between different order symmetry planes were initially investigated in the observable v 2n/Ψ n [30][31][32][33]. This was followed by measurements of nonlinear flow modes of higher harmonics by ALICE [14][15][16]33] as well as event-plane correlations by ATLAS [28].
The correlation observables involving only anisotropic flow coefficients v m and v n were at first measured with event-shape engineering studies [13] proceeded by investigations using symmetric cumulants [34], To study such correlations without the dependence on individual flow coefficients, the normalized symmetric cumulant NSC(m, n) was further proposed [29]. It was found that NSC (3,2), which studies the correlations between v 2 2 and v 2 3 , is very sensitive to the initial conditions and can be used as a good tool to probe initial state ε 2 2 and ε 2 3 correlations. On the other hand, NSC(4, 2) and also NSC involving higher order flow coefficients, are sensitive to both initial conditions and the QGP properties. Thus, these NSC measurements have the potential to distinguish between various models of QGP evolution in hydrodynamic and transport models [9,[34][35][36][37].
It is evident that the study of correlations between various moments of different flow coefficients will deepen our knowledge of the joint p.d. f . for flow magnitudes and angles. However, only correlations involving the second moments of two flow coefficients, v 2 n and v 2 m , have been measured utilizing SC(m, n) while the rest have not yet been explored in experiments. In this Letter, an additional step has been made in this direction by using mixed harmonic cumulants (MHC) [38] to investigate correlations involving more than two different flow coefficients and to study the relationship between higher moments of different flow coefficients in heavy-ion collisions at the LHC. These new measurements establish a milestone for the study of the underlying p.d. f . from P(v n ) to P(v n , v m ...), and significantly improve the overall understanding of the initial conditions and the transport properties of the created QGP at the LHC.

Observables and methods
The multiparticle cumulant of mixed harmonics that involves only flow coefficients, named MHC, was introduced in Ref. [38]. It is defined as an m-observable cumulant [39] in terms of azimuthal angles. By construction, lower order correlations have been subtracted to form genuine multiparticle correlations. Thus, MHC is expected to be insensitive to non-flow effects. This was confirmed in the study of MHC using the HIJING model [40], which does not generate collective flow phenomena [38]. For MHC involving only two flow coefficients of second-order, it is identical with the previously defined fourparticle symmetric cumulants, i.e. MHC(v 2 m , v 2 n ) = SC(m, n). The six-particle cumulant MHC involving v 4 2 and v 2 3 is Here the double angular brackets indicate the averaging procedure performed first over all possible combinations of m-particle tuples that form the m-particle correlation and subsequently the weighted average of all events is calculated with the number of combinations used as an event weight [41]. In the above expressions, lower order (i.e. two-and four-particle) correlations were removed from the six-particle correlation, which results in a genuine six-particle correlation between v 4 2 and v 2 3 . One can rewrite the expression in terms of flow coefficients v 2 and v 3 , To study genuine multiparticle correlations that are independent of the magnitude of the flow coefficients, the normalized mixed harmonic cumulants nMHC involving two flow coefficients v m and v n are constructed according to: Here m k = n l to ensure that nMHC(v k m , v l n ) does not contain flow symmetry plane correlations. The expression of Eq. 10 is also independent of the magnitudes of v m and v n and can therefore be used to quantitatively compare genuine correlations between v k m and v l n determined from experimental data to those determined from the model calculations. Analogously, for MHC involving three flow coefficients without flow symmetry plane correlations, we define the corresponding nMHC, Here m k = n l = p q , and the sum of any two of m k , n l and p q is not equal the third term to avoid flow symmetry plane correlations. Since systematic studies of v n coefficients were carried out for n = 1 − 9 in a previous work [16], this Letter focuses only on the normalized measurements to avoid repeating the earlier discussions on the v n coefficients themselves. Additionally, multi-particle correlations with sub-event method have been used in the normalizations [14,43], to suppress potential non-flow contamination.
In general, one should be able to construct arbitrary mixed harmonic cumulants to any order. However, due to the limited amount of data available, mixed harmonic cumulants higher than the eighth order will not be examined here. All of the previously mentioned two-and multiparticle azimuthal correlations can be measured by using the latest development of the generic algorithm for multiparticle azimuthal correlations [38].

Data sets and systematic uncertainty
This analysis uses data sample from Pb-Pb collisions at √ s NN = 5.02 TeV recorded with the ALICE detector [44,45] during the LHC Run 2 (year 2015) data-taking period. Minimum bias events were triggered by a coincidence signal in the two scintillator arrays of the V0 detector, V0A and V0C, which cover the pseudorapidity ranges of 2.8 < η < 5.1 and −3.7 < η < −1.7, respectively [46]. Only events with a reconstructed primary vertex within ±10 cm from the nominal interaction point along the beam direction were used in this analysis. Removal of background events from, e.g., beam interactions with the residual gas molecules in the beam pipe and pileup events was performed based on the information from the Silicon Pixel Detector (SPD) and the V0 detector. A sample of 55 × 10 6 Pb-Pb collisions, which passed these event selection criteria, were used for the analysis.
Charged tracks were reconstructed using the Inner Tracking System (ITS) [47] and the Time Projection Chamber (TPC) [48]. The selected tracks are required to have at least 70 TPC space points (out of a maximum of 159), and the average χ 2 per degree of freedom of the track fit to the TPC space points is required to be lower than two. Additionally, a minimum of two hits are required in the ITS to improve the momentum resolution. A selection requiring the pseudorapidity to be within |η| < 0.8 is applied. Tracks with a transverse momentum p T < 0.2 GeV/c or p T > 5.0 GeV/c were rejected due to the magnetic field and to reduce the contribution from jets, respectively [49]. In addition, a criterion on the maximum distance of closest approach of the track to the collision point of less than 2 cm in the longitudinal direction and less than a p T -dependent selection in the transverse direction, ranging from 0.2 cm at p T = 0.2 GeV/c to 0.016 cm at p T = 5.0 GeV/c, was applied. This results in a residual contamination from secondary particles from weak decays and from interactions in the detector material of 1-3%, which is negligible in the final systematic uncertainty. These selection criteria result in a transverse momentum dependent efficiency of track reconstruction of about 80%.
Numerous potential sources of systematic uncertainty were investigated in the analysis, including variations of the event and track selection and the uncertainty associated with possible remaining non-flow effects. These are the azimuthal angle correlations not associated with the common flow symmetry planes, including contributions from jets, resonance decays, and are estimated using the HIJING model and found to be negligible for all of the presented observables. The variation of the results with the choice of collision centrality was calculated by alternatively using the SPD to estimate the event multiplicity and was found to contribute less than 5% for all observables. Results with opposite polarities of the magnetic field within the ALICE detector and with narrowing the nominal ±10 cm range of the reconstructed vertex along the beam direction from the center of the ALICE detector to 9, 8, and 7 cm showed a difference of 0-5.4% compared to results with the default selection criteria. The contribution from pileup events was investigated by varying the selections on the correlations between multiplicities from the V0 and SPD, and was found to be negligible. The sensitivity to the track selection criteria was explored by varying the number of TPC space points and by comparing the results to those obtained with tracks with different requirements on hits in the ITS. The effect of varying the number of TPC space points from 70 to 80 and 90 resulted in a negligible systematic uncertainty. Using different track requirements led to a difference with respect to the default selection criteria of less than 4.3% except for nMHC where it was about 16%. The systematic uncertainty evaluated for each above-mentioned source found to be statistically significant according to the recommendation in [50] were added in quadrature to obtain the measurements' total systematic uncertainty.

Results
The centrality dependence of mixed harmonic cumulants with two and three flow coefficients are measured in Pb-Pb collisions at Fig. 1   ) are observed for all centralities, which means that v 2 2 and v 2 4 are correlated while v 2 2 and v 2 3 are anticorrelated. This indicates that finding v 2 larger than v 2 in an event enhances the probability of finding v 4 larger than v 4 and v 3 smaller than v 3 in that event. For nMHC(v 2 3 , v 2 4 ), a similar centrality dependence as for nMHC(v 2 2 , v 2 3 ) is seen for centralities above 20-30% where the nonlinear hydrodynamic response of − → V 4 plays a significant role [14,16]. These new measurements are compared in Fig. 1 to the previously published results at √ s NN = 2.76 TeV, which were named SC(m, n) in Ref. [29,51], shown with open markers. The results of nMHC(v 2 2 , v 2 3 ), nMHC(v 2 2 , v 2 4 ) are compatible within uncertainties at the two different energies, which indicates a weak dependence on the collision energy of these two observables. However, there are differences for nMHC(v 2 3 , v 2 4 ) between the two studied energies, which increase towards central collisions. In particular, the measurement at 5.02 TeV changes sign from negative to positive in central collisions, while it remains negative at the lower energy. A similar study of multiparticle cumulants in the most central collisions was investigated in great detail in Ref. [52], where a significant effect from centrality fluctuations was found in Pb-Pb collisions at 5.02 TeV. Moreover, the amplitude of the centrality fluctuations depends on how the centrality was determined.
Besides the measurements of correlations between two flow coefficients, the new measurement of correlations between three flow coefficients, , is shown with green diamonds in Fig. 1.
, which has been recently measured at a lower energy [53]. By construction, the lower order few-particle correlations have been subtracted from the higher order correlations in the nMHC. Thus, it is expected that nMHC(v 2 2 , v 2 3 , v 2 4 ) should be consistent with zero, if the correlations between three flow coefficients are purely driven by the correlations between two flow coefficients. It is seen in Fig. 1  In order to gain more information on the initial conditions and transport properties of the created QGP at the LHC, the results are compared with those from hydrodynamic model calculations. Results from the hybrid iEBE-VISHNU model with TRENTo initial conditions with specific shear viscosity η/s(T ) and bulk viscosity ζ /s(T ) extracted from the best fit of a Bayesian analysis [55] as well as calculations with AMPT-initial conditions with η/s = 0.08 and no bulk viscosity [54] are compared to the data. Both calculations can quantitatively describe the flow coefficients from inclusive and identified hadrons [54,56] and also provide a reasonable description of more complicated flow observables, e.g., nonlinear modes of higher-order flow [15,16]. Besides these two calculations, the iEBE-VISHNU model with AMPT-initial conditions with η/s = 0.20 and no bulk viscosity is also used. This model does not describe the particle spectra nor the flow coefficients and thus should not be compared with the experimental data. In the remaining text, the hydrodynamic calculations using AMPT initial conditions and η/s = 0.08 will be refereed to as "AMPT calculations". However, the comparison of hydrodynamic calculation from the same initial state model but with different η/s values can be very useful to study the sensitivity of various nMHC to the η/s of the QGP.
Comparisons of the measured nMHC(v 2 2 , v 2 3 ) and nMHC(v 2 2 , v 2 4 ) to hydrodynamic calculations are shown in Fig. 2. In general, the hydrodynamic calculations with both AMPT and TRENTo initial conditions, shown as blue and red shadowed bands, respectively, predict qualitatively the centrality dependence of nMHC. In addition, as v 2 and v 3 are linearly correlated with the initial ε 2 and ε 3 in central and semicentral collisions, compatible results of the final-state nMHC(v 2 2 , v 2 3 ) calculations and the nMHC(ε 2 2 , ε 2 3 ) calculations from the initial-state models are expected [38]. This is indeed shown by the shaded areas and  the dashed lines in Fig. 2 (left). In the same figure, there is also no difference between the calculations using AMPT-initial conditions with different η/s values. This suggests that for the presented centrality ranges, v 2 2 (v 2 3 ) is linearly correlated with the initial ε 2 2 (ε 2 3 ). Thus, the nMHC(v 2 2 , v 2 3 ) measurements shown in Fig. 2 (left) can be used to directly constrain the correlations between the initial anisotropy coefficients ε 2 2 and ε 2 3 without much consideration of the exact value of the transport coefficients in the hydrodynamic models. For nMHC(v 2 2 , v 2 4 ) results shown in Fig. 2 (right), both calculations underestimate the data; the TRENTo calculation fits the data better in central collisions, while the AMPT calculation works slightly better for centralities above 20%. The initial-state calculations of nMHC(ε 2 2 , ε 2 4 ) are significantly lower than the final-state nMHC(v 2 2 , v 2 4 ) calculations, which suggests that the correlation between v 2 2 and v 2 4 is not driven solely by the initial correlation between ε 2 2 and ε 2 4 , but it is mainly developed at later stages of the system's dynamic evolution, especially the nonlinear response contribution to v 4 . , v 2 4 ) measurement. In general, both models generate the same trend of centrality dependence as is seen in data. Notably, the AMPT calculations also predict the sign change in central collisions, while the TRENTo calculations remain negative for the entire centrality range. It has also been seen in Ref. [54] that the AMPT calculations always predict a positive correlation in the most central collisions at 2.76 and 5.02 TeV, while the TRENTo calculations are always negative at both collision energies. Although the hydrodynamic calculations of nMHC(v 2 3 , v 2 4 ) from AMPT and TRENTo initial conditions are almost compatible for non-central collisions, the initial correlations between ε 2 3 and ε 2 4 , quantified by nMHC(ε 2 3 , ε 2 4 ), are utterly different from the two initial-state models, and are far away from the final-state nMHC(v 2 3 , v 2 4 ) calculations. It could be attributed to a significant nonlinear hydrodynamic response in v 4 from ε 2 2 . This nonlinear contribution is strongly anti-correlated with ε 2 3 and plays a dominant role in the final nMHC(v 2 3 , v 2 4 ) results for non-central collisions. On the other hand, in the same centrality region, the linear response of v 4 to ε 4 is rather weak [14], and the contributions from correlations between the initial ε 2 3 and ε 2 4 in the final nMHC(v 2 3 , v 2 4 ) appear to be minor. To extend the discussion from correlations of two flow coefficients to three flow coefficients, the mea-surement of nMHC(v 2 2 , v 2 3 , v 2 4 ) and its comparison to hydrodynamic calculations with both AMPT and TRENTo initial conditions are presented in Fig. 3 (right). In general, the agreement between the initial nMHC(ε 2 2 , ε 2 3 , ε 2 4 ) correlations and the final-state nMHC(v 2 2 , v 2 3 , v 2 4 ) calculations worsens as the collision centrality becomes more peripheral, which can be expected due to the increasing contribution from the nonlinear hydrodynamic response in v 4 . Figure 3 (right) also shows clearly that the calculation with AMPT initial conditions and η/s = 0.08 describes the data reasonably well. The calculation with η/s = 0.20 is two times larger than the one with η/s = 0.08. Such a difference is more significant compared to what has been seen in the correlations between two harmonics, where no obvious difference is observed for nMHC(v 2 2 , v 2 3 ) and nMHC(v 2 3 , v 2 4 ) and only a relatively small difference is seen for nMHC(v 2 2 , v 2 4 ). This demonstrates the novelty of the new correlations between three flow coefficients constraining the transport properties of the QGP. However, despite the fact that the hydrodynamic calculations using TRENTo initial conditions are consistent with the measured nMHC(v 2 2 , v 2 3 ) and nMHC(v 2 3 , v 2 4 ), and also provide a reasonable description of the nMHC(v 2 2 , v 2 4 ) measurement, they significantly underestimate the data by roughly a factor of two. Considering an apparent discrepancy between the data and TRENTo calculations, there is little doubt that the hydrodynamic framework and its corresponding parameters can be better tuned in a future Bayesian analysis if this new nMHC(v 2 2 , v 2 3 , v 2 4 ) measurement is used as an input. The first measurement of correlations between three harmonics provides additional independent constraints on the theoretical models beyond those provided by the correlations of two harmonics that have been studied before. Centrality percentile With the recently proposed observable nMHC, one can study not only the correlations between two or three different flow coefficients, in terms of their second moments, but also the correlations between the k th order moment of v m and the l th order moment of v n where k ≥ 2, l ≥ 2. It is particularly interesting to study the correlations between various moments of v 2 and v 3 because in central and semi-central collisions, both v 2 and v 3 are linearly correlated to their corresponding initial eccentricities ε 2 and ε 3 [8,9]. Thus the measurement of nMHC(v k 2 , v l 3 ) in central and semi-central collisions might provide a direct constraint on the initial correlation between ε k 2 and ε l 3 . This information is extremely important for the understanding of the initial conditions of heavy-ion collisions but it has never been measured before. Conversely, the potential nonlinearity of v 2 and v 3 , more pronounced in peripheral collisions, strongly depends on the dynamical evolution of the created QGP. The study of nMHC(v k 2 , v l 3 ) will enable a new way to study this effect on v 2 and v 3 . The first measurements of nMHC(v k 2 , v l 3 ) are presented in Fig. 4. In addition to the negative value of nMHC(v 2 2 , v 2 3 ), which has been studied before, positive correlations are observed for both six-particle mixed harmonic cumulants nMHC(v 4 2 , v 2 3 ) and nMHC(v 2 2 , v 4 3 ). The eight-particle mixed harmonic cumulants nMHC(v 6 2 , v 2 3 ), nMHC(v 4 2 , v 4 3 ) and nMHC(v 2 2 , v 6 3 ) are all negative. Such characteristic negative, positive and negative signs of four-, six-, and eight-particle mixed harmonic cumulants, respectively, are very similar to the previously measured pattern for two-, four-, six,-and eight-particle single harmonic cumulants in Pb-Pb collisions [21], which show positive, negative, positive, and negative signs, respectively. These findings agree qualitatively with the initial-state predictions based on the MC-Glauber [38], AMPT, and TRENTo models [54]. It should be pointed out that the measured negative nMHC(v 2 2 , v 2 3 ) shown above could only confirm the negative correlations of (v 2 2 , v 2 3 ), while the results presented in Fig. 4 illustrate further the positive correlations of (v 4 2 , v 2 3 ) and (v 2 2 , v 4 3 ) as well as the negative correlations of (v 6 2 . Moreover, one can see the following hierarchy, |nMHC(v 6 2 This agrees qualitatively with the predictions based on initial-state models [38,54]. Furthermore, the calculations based on the HIJING model [40], which does not generate anisotropic flow in the created system, are consistent with zero [38], and thus do not reproduce the characteristic signs of the multiparticle mixed harmonic cumulants observed in experiments. Future studies with other non-flow models, i.e. PYTHIA [57,58], could confirm if the aforementioned characteristic signs of the multiparticle mixed harmonic cumulants can be regarded as a flow signature and thus could be used for searching for collective flow in small collision systems like pp or pA collisions [59][60][61].
As mentioned above, for non-peripheral collisions, both v 2 and v 3 are expected to be linearly correlated with the initial eccentricity ε 2 and ε 3 . Thus, the final-state result of nMHC(v k 2 , v l 3 ) could reflect the initial correlation between ε k 2 and ε l 3 . This behavior is observed in the case of nMHC(v 2 2 , v 2 3 ), where good agreement with nMHC(ε 2 2 , ε 2 3 ) was found. Moving to higher moments of v 2 and/or v 3 , one can further probe the nonlinearity of v 2 (v 3 ) to ε 2 (ε 3 ) by seeing if the agreement between the initial and final-state correlations persists, because of the better sensitivity of higher moments to the nonlinear hydrodynamic response. Figure 5 presents the comparison of data with iEBE-VISHNU calculations with , independent of whether the AMPT or TRENTo initial-state model are used. For AMPT calculations, there is no difference between the results using η/s = 0.08 or 0.20, which confirms that the precision measurement of nMHC(v 4 2 , v 2 3 ) can offer an additional approach to constrain initial-state models. However, the situation is different in the case of nMHC(v 2 2 , v 4 3 ), shown in Fig. 5 (right). The hydrodynamic calculations with both AMPT and TRENTo initial conditions are compatible with the measurement within the considerable uncertainty, but there is an apparent discrepancy between nMHC(v 2 2 , v 4 3 ) and nMHC(ε 2 2 , ε 4 3 ). This does not agree with the naive expectation of both v 2 and v 3 being linearly correlated with their respective initial eccentricities ε 2 and ε 3 , which might be because, generally, the linearity of v 3 to ε 3 is worse than that of v 2 to ε 2 as shown by hydrodynamic calculations [7]. When one examines higher-order moments, the linear response of v 4 2 remains and thus nMHC(ε 4 2 , ε 2 3 ) = nMHC(v 4 2 , v 2 3 ) is observed. However, the nonlinearity of v 4 3 becomes non-negligible in non-peripheral collisions, which creates the discrepancy between the initial nMHC(ε 2 2 , ε 4 3 ) and final state nMHC(v 2 2 , v 4 3 ) correlations. This hypothesis is further confirmed in Figs. 6 and 7 where eight-particle cumulants are reported, which involve even higher moments of v 2 and/or v 3 . Firstly, hydrodynamic calculations with AMPT initial conditions quantitatively predict the new measurements of nMHC(v 6 2 , while the TRENTo calculations show compatible results except for nMHC(v 4 2 , v 4 3 ), where the calculations overestimate the data by a factor of two. Secondly, although in hydrodynamic calculations with two different initial conditions there is a similar centrality dependence of nMHC(ε 6 2 , ε 2 3 ) and nMHC(v 6 2 , v 2 3 ), a clear difference between the two is observed already in semi-central collisions. This difference becomes much larger when fourth and sixth order moments of v 3 are involved, with no obvious agreement between the initial and final-state calculations; this is especially shown by the TRENTo calculations in Figs. 6 (right) and 7. It is expected that the effect of the nonlinearity of v 2 and v 3 will be enhanced when studying nMHC(v k 2 , v l 3 ) with higher moments (e.g. k, l ≥ 4). The resulting nMHC(v k 2 , v l 3 ) in the final state, instead of being determined solely by the initial correlations of nMHC(ε k 2 , ε l 3 ), receive non- negligible contributions from the nonlinearities of higher moments of v 2 and v 3 , developed during the dynamic evolution of the system. Thus, the new measurements of nMHC(v k 2 , v l 3 ) presented in this Letter provide direct access to the initial correlations between ε k 2 and ε l 3 when lower moments of v 2 and v 3 are involved, while enabling a new possibility to study the nonlinearities of v 2 and v 3 when higher moments are involved.

Summary
The normalized mixed harmonic cumulants nMHC between two and three flow coefficients as well as between higher moments of two flow coefficients, were measured in √ s NN = 5.02 TeV Pb-Pb collisions with ALICE. It is found that nMHC(v 2 2 , v 2 4 ) is positive, while nMHC(v 2 2 , v 2 3 ) and nMHC(v 2 3 , v 2 4 ) are negative. In addition, the first measurement of three harmonic correlations nMHC(v 2 2 , v 2 3 , v 2 4 ) is closer to nMHC(v 2 2 , v 2 4 ) and nMHC(v 2 3 , v 2 4 ) in central collisions, and then becomes closer to nMHC(v 2 2 , v 2 3 ) and nMHC(v 2 3 , v 2 4 ) for more peripheral collisions. These measurements compared with iEBE-VISHNU hydrodynamic calculations using AMPT and TRENTo initial conditions exhibit different sensitivities to the initial conditions and the specific shear viscosity of the QGP. Thus the measurements presented in this Letter can be used to more tightly constrain theoretical models. Furthermore, the correlations between higher moments of v 2 and v 3 were investigated for the first time. The four-, six-and eight-particle mixed cumulants of nMHC show the characteristic signature of negative, positive, and negative signs, respectively, similar to the multiparticle cumulants of single harmonics. The comparison with hydrodynamic calculations reveals that the correlations involving higher-order moments could significantly enhance the contributions that arise from nonlinearities of v 2 and v 3 to the initial eccentricity ε 2 , triangularity ε 3 , respectively. Such contributions mainly develop during the expansion of the system and reflect the time evolution of the shear and bulk viscosities of the QGP. These new measurements of correlations between different moments of two and three flow coefficients, together with comparisons to state-of-the-art hydrodynamic calculations, provide further information on the initial conditions and considerably tighten the constraints on the evolution of the QGP created in heavy-ion collisions at the LHC.

Acknowledgments
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The