Systematic studies of correlations between different order flow harmonics in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV

The correlations between event-by-event fluctuations of anisotropic flow harmonic amplitudes have been measured in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV with the ALICE detector at the Large Hadron Collider. The results are reported in terms of multiparticle correlation observables dubbed Symmetric Cumulants. These observables are robust against biases originating from nonflow effects. The centrality dependence of correlations between the higher order harmonics (the quadrangular $v_4$ and pentagonal $v_5$ flow) and the lower order harmonics (the elliptic $v_2$ and triangular $v_3$ flow) is presented. The transverse momentum dependences of correlations between $v_3$ and $v_2$ and between $v_4$ and $v_2$ are also reported. The results are compared to calculations from viscous hydrodynamics and A Multi-Phase Transport ({AMPT}) model calculations. The comparisons to viscous hydrodynamic models demonstrate that the different order harmonic correlations respond differently to the initial conditions and the temperature dependence of the ratio of shear viscosity to entropy density ($\eta/s$). A small average value of $\eta/s$ is favored independent of the specific choice of initial conditions in the models. The calculations with the AMPT initial conditions yield results closest to the measurements. Correlations between the magnitudes of $v_2$, $v_3$ and $v_4$ show moderate $p_{\rm T}$ dependence in mid-central collisions. This might be an indication of possible viscous corrections to the equilibrium distribution at hadronic freeze-out, which might help to understand the possible contribution of bulk viscosity in the hadronic phase of the system. Together with existing measurements of individual flow harmonics, the presented results provide further constraints on the initial conditions and the transport properties of the system produced in heavy-ion collisions.


Introduction
The main emphasis of the ultra-relativistic heavy-ion collision programs at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) is to study the deconfined phase of strongly interacting QCD matter, the Quark-Gluon Plasma (QGP). The matter produced in a heavy-ion collision exhibits strong collective radial expansion [1,2]. Difference in pressure gradients and the interactions among matter constituents produced in the spatially anisotropic overlap region of the two colliding nuclei result in anisotropic transverse flow in the momentum space. The large elliptic flow discovered at RHIC energies [3][4][5][6][7] is also observed at LHC energies [8][9][10][11][12][13][14][15][16][17][18]. The measurements are well described by calculations utilizing viscous hydrodynamics [19][20][21][22][23][24]. These calculations also demonstrated that the shear viscosity to the entropy density ratio (η/s) of the QGP in heavy-ion collisions at RHIC and LHC energies is close to a universal lower bound 1/4π [25].
The temperature dependence of η/s has some generic features typical to the most known fluids. This ratio reaches its minimum value close to the phase transition region [25,26]. It was shown, using kinetic theory and quantum mechanical considerations [27], that η/s ∼ 0.1 would be the correct order of magnitude for the lowest possible shear viscosity to entropy density ratio value found in nature. Later it was demonstrated that an exact lower bound (η/s) min = 1/4π ≈ 0.08 can be conjectured using AdS/CFT correspondence [25]. Hydrodynamical simulations constrained by data support the view that η/s of the QGP is close to that limit [23]. It is argued that such a low value might imply that thermodynamic trajectories for the expanding matter would lie close to the quantum chromodynamics (QCD) critical end point, which is another subject of intensive experimental study [26,28].
Anisotropic flow [29] is quantified with n th -order flow harmonics v n and corresponding symmetry plane angles Ψ n in a Fourier decomposition of the particle azimuthal distribution in the plane transverse to the beam direction [30,31]: where E, p, p T , ϕ and η are the particle's energy, momentum, transverse momentum, azimuthal angle and pseudorapidity, respectively, and Ψ n is the azimuthal angle of the symmetry plane of the n th -order harmonic. Harmonic v n can be calculated as v n = cos[n(ϕ − Ψ n )] , where the angular brackets denote an average over all particles in all events. The anisotropic flow in heavy-ion collisions is typically understood as the hydrodynamic response of the produced matter to spatial deformations of the initial energy density profile [32]. This profile fluctuates event-by-event due to fluctuating positions of the constituents inside the colliding nuclei, which implies that v n also fluctuates [33,34]. The recognition of the importance of flow fluctuations led to the discovery of triangular and higher flow harmonics [9,35] as well as to the correlations between different v n harmonics [36,37]. The higher order harmonics are expected to be sensitive to fluctuations in the initial conditions and to the magnitude of η/s [38,39], while v n correlations have the potential to discriminate between these two respective contributions [36].
Difficulties in extracting η/s in heavy-ion collisions can be attributed mostly to the fact that it strongly depends on the specific choice of the initial conditions in the models used for comparison [19,39,40]. Viscous effects reduce the magnitude of the anisotropic flow. Furthermore, the magnitude of η/s used in hydrodynamic calculations should be considered as an average over the temperature evolution of the expanding fireball as it is known that η/s depends on temperature. In addition, part of the anisotropic flow can also originate from the hadronic phase [41][42][43]. Therefore, both the temperature dependence of η/s and the relative contributions from the partonic and hadronic phases should be understood better to quantify the η/s of the QGP.
An important input to the hydrodynamic model simulations is the initial distribution of energy density in the transverse plane (the initial density profile), which is usually estimated from the probability distribution of nucleons in the incoming nuclei. This initial energy density profile can be quantified by calculating the distribution of the spatial eccentricities ε n [35], ε n e inΦ n = −{r n e inφ }/{r n }, where the curly brackets denote the average over the transverse plane, i.e. {· · · } = dxdy e(x, y, τ 0 ) (· · · ), r is the distance to the system's center of mass, φ is azimuthal angle, e(x, y, τ 0 ) is the energy density at the initial time τ 0 , and Φ n is the participant plane angle (see Refs. [44,45]). There is experimental and theoretical evidence [9, 35,46] that the lower order harmonics, v 2 and v 3 , to a good approximation, are linearly proportional to the deformations in the initial energy density in the transverse plane (e.g. v n ∝ ε n for n = 2 or 3). Higher order (n > 3) flow harmonics can arise from initial anisotropies in the same harmonic [35,44,47,48] (linear response) or can be induced by lower order harmonics [49,50] (nonlinear response). For instance, v 4 can develop both as a linear response to ε 4 and/or as a nonlinear response to ε 2 2 [51]. Therefore, the higher harmonics (n > 3) can be understood as superpositions of linear and nonlinear responses, through which they are correlated with lower order harmonics [47,48,50,52,53]. When the order of the harmonic is large, the nonlinear response contribution in viscous hydrodynamics is dominant and increases in more peripheral collisions [50,52]. The magnitudes of the viscous corrections as a function of p T for v 4 and v 5 are sensitive to the ansatz used for the viscous distribution function, a correction for the equilibrium distribution at hadronic freeze-out [52,54]. Hence, studies of the correlations between higher order (n > 3) and lower order (v 2 or v 3 ) harmonics and their p T dependence can help to understand the viscous correction to the momentum distribution at hadronic freeze-out which is among the least understood parts of hydrodynamic calculations [45,52,55,56].
The first results for new multiparticle observables which quantify the relationship between event-byevent fluctuations of two different flow harmonics, the Symmetric Cumulants (SC), were recently reported by the ALICE Collaboration [57]. The new observables are particularly robust against few-particle nonflow correlations [8] and they provide independent, complementary information to recently analyzed symmetry plane correlators [37]. It was demonstrated that they are sensitive to the temperature dependence of η/s of the expanding medium and therefore simultaneous descriptions of correlations between different order harmonics would constrain both the initial conditions and the medium properties [57,58]. In this article, we have extended the analysis of SC observables to higher order harmonics (up to 5 th order) as well as to the measurement of the p T dependence of correlations for the lower order harmonics (v 3 -v 2 and v 4 -v 2 ). We also present a systematic comparison to hydrodynamic and AMPT model calculations. In Sec. 2 we present the analysis methods and summarize our findings from the previous work [57]. The experimental setup and measurements are described in Sec. 3. The sources of systematic uncertainties are explained in Sec. 4. The results of the measurements are presented in Sec. 5. In Sec. 6 we present comparisons to model calculations. Finally, Sec. 7 summarizes our new results.

Experimental Observables
Existing measurements for anisotropic flow observables provide an estimate of the average value of η/s of the QGP, both at RHIC and LHC energies. What remains uncertain is how the η/s of the QGP depends on temperature (T ). The temperature dependence of η/s of the QGP was discussed in [28]. The effects on hadron spectra and elliptic flow were studied in [59] for different parameterizations of η/s(T ). A more systematic study with event-by-event EKRT+viscous hydrodynamic calculations was recently initiated in [45], where the first (and only rather qualitative) possibilities were investigated (see Fig. 1 therein). The emerging picture is that the study of individual flow harmonics v n alone is unlikely to reveal the details of the temperature dependence of η/s. It was already demonstrated in [45] that different η/s(T ) parameterizations can lead to the same centrality dependence of individual flow harmonics. In Ref. [36] new flow observables were introduced which quantify the degree of correlation between amplitudes of two different harmonics v m and v n . These new observables have the potential to discriminate between the contributions to anisotropic flow development from initial conditions and from the transport properties of the QGP [36]. Therefore their measurement would provide experimental constraints on theoretical parameters used to describe the individual stages of the heavy-ion system evolution. In addition, it turned out that correlations of different flow harmonics are sensitive to the temperature dependence of η/s [57], to which individual flow harmonics are weakly sensitive [45].
For reasons discussed in [57,60], the correlations between different flow harmonics cannot be studied experimentally with the set of observables introduced in [36]. Based on [60], new flow observables obtained from multiparticle correlations, Symmetric Cumulants (SC), were introduced.
The SC observables are defined as: with the condition m = n for two positive integers m and n (for details see Sec. IV C in [60] Normalized symmetric cumulants reflect only the strength of the correlation between v m and v n , while SC(m, n) has contributions from both the correlations between the two different flow harmonics and the individual harmonics. In Eq. (4) the products in the denominator are obtained from two-particle correlations using a pseudorapidity gap of |∆η| > 1.0 which suppresses biases from few-particle nonflow correlations. For the two two-particle correlations which appear in the definition of SC(m, n) in Eq. (3) the pseudorapidity gap is not needed, since nonflow is suppressed by construction in this observable. This was verified by HIJING model simulations in [57].
The ALICE measurements [57] have revealed that fluctuations of v 2 and v 3 are anti-correlated, while fluctuations of v 2 and v 4 are correlated for all centralities [57]. It was found that the details of the centrality dependence differ in the fluctuation-dominated (most central) and the geometry-dominated (mid-central) regimes [57]. The observed centrality dependence of SC(4,2) cannot be captured by models with constant η/s, indicating that the temperature dependence of η/s plays an important role. These results were also used to discriminate between different parameterizations of initial conditions. It was demonstrated that in the fluctuation-dominated regime (central collisions), MC-Glauber initial conditions with binary collision weights are favored over wounded nucleon weights [57]. The first theoretical studies of SC observables can be found in [58,[61][62][63][64][65].

Data Analysis
The data sample of Pb-Pb collisions at the centre-of-mass energy √ s NN = 2.76 TeV analyzed in this article was recorded by ALICE during the 2010 heavy-ion run of the LHC. Detailed descriptions of the ALICE detector can be found in [66][67][68]. The Time Projection Chamber (TPC) was used to reconstruct charged particle tracks and measure their momenta with full azimuthal coverage in the pseudorapidity range |η| < 0.8. Two scintillator arrays (V0A and V0C) which cover the pseudorapidity ranges −3.7 < η < −1.7 and 2.8 < η < 5.1 were used for triggering and the determination of centrality [69]. The trigger conditions and the event selection criteria are identical to those described in [8,69]. Approximately 10 7 minimum-bias Pb-Pb events with a reconstructed primary vertex within ±10 cm from the nominal interaction point along the beam direction are selected. Only charged particles reconstructed in the TPC in |η| < 0.8 and 0.2 < p T < 5 GeV/c were included in the analysis. The charged track quality cuts described in [8] were applied to minimize contamination from secondary charged particles and fake tracks. The track reconstruction efficiency and contamination were estimated from HIJING Monte Carlo simulations [70] combined with a GEANT3 [71] detector model and were found to be independent of the collision centrality. The reconstruction efficiency increases with transverse momenta from 70% to 80% for particles with 0.2 < p T < 1 GeV/c and remains constant at (80 ± 5)% for p T > 1 GeV/c. The estimated contamination by secondary charged particles from weak decays and photon conversions is less than 6% at p T = 0.2 GeV/c and falls below 1% for p T > 1 GeV/c. The p T cut-off of 0.2 GeV/c reduces event-by-event biases due to small reconstruction efficiency at lower p T , while the high p T cut-off of 5 GeV/c reduces the effects of jets on the measured correlations. Reconstructed TPC tracks constrained to vertex are required to have at least 70 space points (out of a maximum of 159). Only tracks with a transverse distance of closest approach to the primary vertex less than 3 mm, both in the longitudinal and transverse directions, are accepted. This reduces the contamination from secondary tracks produced in the detector material, particles from weak decays, etc. Tracks with kinks (i.e. tracks that appear to change direction due to multiple scattering or K ± decays) were rejected.

Systematic Uncertainties
The systematic uncertainties are estimated by varying the event and track selection criteria. All systematic checks described here are performed independently. The SC(m, n) values resulting from each variation are compared to ones from the default event and track selection described in the previous section, and differences are taken as the systematic uncertainty due to each individual source. The contributions from different sources were added in quadrature to obtain the total systematic uncertainty.
The event centrality was determined by the V0 detectors [72] with better than 2% resolution for the whole centrality range analyzed. The systematic uncertainty from the centrality determination was evaluated by using the TPC and Silicon Pixel Detector (SPD) [73] detectors instead of the V0 detectors. The systematic uncertainty on the symmetric cumulants which arises from the centrality uncertainty is about 3% both for SC (5,2) and SC (4,3), and 8% for SC(5,3). As described in Sec. 3, the reconstructed vertex position along the beam axis (z-vertex) is required to be located within 10 cm of the nominal interaction to ensure uniform detector acceptance for tracks within |η| < 0.8. The systematic uncertainty from the z-vertex cut was estimated by reducing the z-vertex range to 8 cm and was found to be less than 3%.
The analyzed events were recorded with two settings of the magnet field polarity and the resulting data sets have almost equal numbers of events. Events with both magnet field polarities were used in the default analysis, and the systematic uncertainties were evaluated from the variation between each of the two magnetic field settings. The uncertainty due to the p T dependence of the track reconstruction efficiency was also taken into account. Magnetic field polarity variation and reconstruction efficiency effects contribute less than 2% to the systematic uncertainty.
The systematic uncertainty due to the track reconstruction procedure was estimated from comparisons between results for the so-called standalone TPC tracks with the same parameters as described in Sec. 3, and tracks from a combination of the TPC and the Inner Tracking System (ITS) detectors with tighter selection criteria. To avoid non-uniform azimuthal acceptance due to dead zones in the SPD, and to get the best transverse momentum resolution, a hybrid track selection utilizing SPD hits and/or ITS refit tracks combined with TPC information was used. Then each track reconstruction strategy was evaluated by varying the threshold on parameters used to select the tracks at the reconstruction level. A systematic difference of up to 12% was observed in SC(m, n) from the different track selections. In addition, we applied the like-sign technique to estimate nonflow contributions [8] to SC(m, n). The difference between results obtained by selecting all charged particles and results obtained after either selecting only positively or only negatively charged particles was the largest contribution to the systematic uncertainty and is about 7% for SC(4,3) and 20% for SC (5,3).  (4,2)) are taken from [57] and shown as bands. The systematic and statistical errors are combined in quadrature for these lower order harmonic correlations. The SC(4,2) and SC (3,2) are downscaled by a factor of 0.1. Systematic uncertainties are represented with boxes for higher order harmonic correlations.
Another large contribution to the systematic uncertainty originates from azimuthal non-uniformities in the reconstruction efficiency. In order to estimate its effects, we use the AMPT model (see Sec. 6) which has a uniform distribution in azimuthal angle. Detector inefficiencies were introduced to mimic the non-uniform azimuthal distribution in the data. For the observables SC(5,2), SC(5,3) and SC(4,3) the variation due to non-uniform acceptance is about 9%, 17% and 11%, respectively. Overall, the systematic uncertainties are larger for SC(5,3) and SC(5,2) than for the lower harmonics of SC(m, n). This is because v n decreases with increasing n and becomes more sensitive to azimuthal modulation due to detector imperfections.

Results
The centrality dependence of the higher order harmonic correlations (SC(4,3), SC(5,2) and SC(5,3)) are presented in Fig. 1 and compared to the lower order harmonic correlations (SC(3,2) and SC(4,2)) which were published in [57]. The correlation between v 3 and v 4 is negative, and similarly for v 3 and v 2 , while the other correlations are all positive, which reveals that v 2 and v 5 as well as v 3 and v 5 are correlated like v 2 and v 4 , while v 3 and v 4 are anti-correlated like v 3 and v 2 .
The higher order flow harmonic correlations are much smaller compared to the lower order harmonic correlations. In particular SC(5,2) is 10 times smaller than SC(4,2) and SC(4,3) is about 20 times smaller than SC (3,2).
Unlike SC(m, n), the NSC(m, n) results with the higher order flow harmonics show almost the same order of the correlation strength as the lower order flow harmonic correlations NSC (3,2) or NSC (4,2). This demonstrates the advantage of using the normalized SC observables in which the correlation strength between flow harmonics is not hindered by the differences in magnitudes of different flow harmonics. The NSC(4,3) magnitude is comparable to NSC(3,2) and one finds that a hierarchy, NSC(5,3) > NSC(4,2) > NSC(5,2), holds for centrality range 20-50% within the errors as shown in Fig. 1(b). The SC(5,2) magnitude is larger than SC(5,3), but the normalized correlation between v 5 and v 3 is stronger than the normalized correlation between v 5 and v 2 . These results indicate that the lower order harmonic correlations are larger than higher order harmonic correlations, not only because of the correlation strength itself but also because of the strength of the individual flow harmonics.  It can be seen in Fig. 1(a) that the lower order harmonic correlations as well as SC(5,2) increase nonlinearly towards peripheral collisions. In the case of SC(5,3) and SC(4,3), the centrality dependence is weaker than for the other harmonic correlations. The NSC(5,3) observable shows the strongest normalized correlation among all harmonics while NSC(5,2) shows the weakest centrality dependence. Both NSC (3,2) and NSC(4,3) are getting more anti-correlated toward peripheral collisions and have the similar magnitude.
To study the p T dependence of SC(m, n), we present the results as a function of the low p T cut-off (p T,min ), instead of using independent p T intervals; this decreases large statistical fluctuations in the results. Various minimum p T cuts from 0.2 to 1.5 GeV/c are applied. The p T dependent results for SC (3,2) and SC(4,2) as a function of minimum p T cuts are shown in Figs. 2a and 2c. The strength of SC(m, n) becomes larger as p T,min increases. The centrality dependence is stronger with higher p T,min cuts, with SC(m, n) getting much larger as centrality or p T,min increases. The NSC(3,2) and NSC(4,2) observables with different p T,min are shown in Figs. 2b and 2d. The strong p T,min dependence observed in SC(m, n) is not seen in NSC(m, n). This indicates that the p T dependence of SC(m, n) is dominated by the p T dependence of the individual flow harmonics v n . The p T,min dependence of NSC(3,2) is not clearly seen and it is consistent with no p T,min dependence within the statistical and systematic errors for the centrality range 0-30%, while showing a moderate increase of anti-correlation with increasing p T,min for the 30-50% centrality. The NSC(4,2) observable shows a moderate decreasing trend as p T,min increases. These observations are strikingly different from the p T dependence of the individual flow harmonics, where the relative flow fluctuations σ v 2 / v 2 [74] are independent of transverse momentum up to p T ∼ 8 GeV/c (see Fig. 3 in Ref. [75]).
As discussed in Sec. 2, the NSC(m,n) observables are normalized by the product v 2 m v 2 n . These products are obtained from two-particle correlations using a pseudorapidity gap of |∆η| > 1.0. In this

Model Comparisons
We have performed a systematic comparison of the centrality and transverse momentum dependence of the SC(m,n) and NSC(m,n) to the event-by-event EKRT+viscous hydrodynamics [45], VISH2+1 [76,77], and the AMPT [63,78,79] models. Comparisons for v n coefficients with the model calculations are presented in Appendix A.
In the event-by-event EKRT+viscous hydrodynamic calculations [45], the initial energy density profiles are calculated using a next-to-leading order perturbative-QCD+saturation model [80,81]. The subsequent space-time evolution is described by relativistic dissipative fluid dynamics with different parameterizations for the temperature dependence of the shear viscosity to entropy density ratio η/s(T ). This model gives a good description of the charged hadron multiplicity and the low p T region of the charged hadron spectra at RHIC and the LHC (see Figs. 11-13 in [45]). Each of the η/s(T ) parameterizations is adjusted to reproduce the measured v n from central to mid-peripheral collisions (see Fig. 15 in [45] and Appendix A).
The VISH2+1 [76,77] event-by-event calculations for relativistic heavy-ion collisions are based on (2+1)-dimensional viscous hydrodynamics which describes the QGP phase and the highly dissipative and off-equilibrium late hadronic stages with fluid dynamics. By tuning transport coefficients and decoupling temperature for a given scenario of initial conditions, it can describe the p T spectra and different flow harmonics at RHIC and the LHC [20, 76,82,83] energies. Three different types of initial conditions [58] (MC-Glauber, MC-KLN and AMPT) along with different constant η/s values have been used for our data to model comparisons. Traditionally, the Glauber model constructs the initial entropy density from the wounded nucleon and binary collision density profiles [84]. The KLN model assumes that the initial energy density is proportional to that of the initial gluons calculated from the corresponding k T factorization formula [85]. In Monte Carlo versions MC-Glauber and MC-KLN [86][87][88] of these models an additional initial state fluctuations are introduced through position fluctuations of individual nucleons inside the colliding nuclei. For the AMPT initial conditions [83,89,90], the fluctuating energy density profiles are constructed from the energy distribution of individual partons, which fluctuate in both momentum and coordinate space. Compared with the MC-Glauber and MC-KLN initial conditions, the additional Gaussian smearing in the AMPT initial conditions gives rise to non-vanishing initial local flow velocities [89].
Even though thermalization could be achieved shortly in collisions of very large nuclei and/or at extremely high energy [91], the dense matter created in heavy-ion collisions may not reach full thermal or chemical equilibrium due to its finite size and short lifetime. To address such non-equilibrium manybody dynamics, the AMPT model [78,92,93] has been developed, which includes both initial partonic and final hadronic interactions and the transition between these two phases of matter. The initial conditions in the AMPT are given by the spatial and momentum distributions of minijets and soft strings from the HIJING model [70,94]. For the data comparisons three different configurations of the AMPT model have been used: the default one and string melting with and without hadronic rescattering. The input parameters used in all configurations are: α s = 0.33, a partonic cross-section of 1.5 mb. In the default configuration, partons are recombined with their parent strings when they stop interacting. The resulting strings are later converted into hadrons using the Lund string fragmentation model [95,96]. The Lund string fragmentation parameters were set to α = 0.5 and b = 0.9 GeV −2 . In the string melting configuration, the initial strings are melted into partons whose interactions are described by the ZPC parton cascade model [97]. These partons are then combined into the final state hadrons via a quark coalescence model. In both configurations, the dynamics of the subsequent hadronic matter is described by a hadronic cascade based on A Relativistic Transport (ART) model [98] which includes resonance decays. The string melting configuration of the AMPT without hadronic rescattering was used to study the influence of the hadronic phase on the development of the anisotropic flow. Even though the string melting version of AMPT [78,99] reasonably well reproduces particle yields, p T spectra, and v 2 of low p T pions and kaons in central and mid-central Au-Au collisions at √ s NN = 200 GeV and Pb-Pb collisions at √ s NN = 2.76 TeV [79], it was observed in a recent study [100] that it fails to quantitatively reproduce the flow harmonics of identified hadrons (v 2 , v 3 , v 4 and v 5 ) at √ s NN = 2.76 TeV. It turns out that the radial flow in AMPT is 25% lower than that measured at the LHC, which is responsible for this quantitative disagreement [100]. The details of the AMPT configurations used in this article and the comparisons of p T -differential v n for pions, kaons and protons to the data can be found in [100].

Centrality Dependence of SC(m,n) and NSC(m,n)
Comparison to event-by-event EKRT+viscous hydrodynamic predictions with various parameterizations of the temperature dependence of η/s(T ) was shown in Fig. 2 of Ref. [57]. It was demonstrated that NSC (3,2) is sensitive mainly to the initial conditions, while NSC(4,2) is sensitive to both the initial conditions and the system properties, which is consistent with the predictions from [36]. The model calculations for NSC(4,2) observable show that it has better sensitivity for different η/s(T ) parameterizations but they cannot describe neither the centrality dependence nor the absolute values. The discrepancy between data and theoretical predictions indicates that the current understanding of initial conditions in models of heavy-ion collisions needs to be revisited to further constrain η/s(T ). are compared to the data in Fig. 4. As can be seen in Fig. 1 from Ref. [45], for the "param1" parameterization the phase transition from the hadronic to the QGP phase occurs at the lowest temperature, around 150 MeV. This parameterization is also characterized by a moderate slope in η/s(T ) which decreases (increases) in the hadronic (QGP) phase. The model calculations in which the temperature of the phase transition is larger than for "param1" are ruled out by the previous measurements [57]. While the correlations between v 5 and v 2 are well described at all centralities, the correlations between v 5 and v 3 are reproduced in 0-40% centrality range and deviate by about one sigma for 40-50% centrality. In the case of v 4 and v 3 , the same models underestimate the anti-correlation in the data significantly in mid-central collisions and fail similarly for the anti-correlation between v 3 and v 2 .
The comparison to the VISH2+1 calculation [58] is shown in 0.16 for AMPT initial conditions) fail to describe the centrality dependence of the SC(m,n) observables of all orders, shown in the left panels in Fig. 5. Among the calculations with small η/s (η/s = 0.08), the one with the AMPT initial conditions describes the data better than the ones with other initial conditions for all SC(m,n) observables measured, but it cannot describe the data quantitively for most of the centrality ranges.
However, NSC(4,2) is sensitive both to the initial conditions and the η/s parameterizations used in the models. Even though NSC(4,2) favors both AMPT initial conditions with η/s = 0.08 and MC-Glauber initial conditions with η/s = 0.20, SC(4,2) can only be described by models with smaller η/s. Hence the calculation with large η/s = 0.20 is ruled out. We conclude that η/s should be small and that AMPT initial conditions are favored by the data. The NSC(5,2) and NSC(5,3) observables are quite sensitive to both the initial conditions and the η/s parameterizations. The SC(4,3) results clearly favor smaller η/s values but NSC(4,3) cannot be described by these models quantitively.
The SC(m,n) and NSC(m,n) observables calculated from AMPT simulations are compared with data in Fig. 6. For SC (3,2), the calculation with the default AMPT settings is closest to the data, but none of the AMPT configurations can describe the data fully. The third version based on the string melting configuration without the hadronic rescattering phase is also shown. The hadronic rescattering stage makes both SC(3,2) and NSC (3,2) smaller in the string melting AMPT model but not enough to describe the data. Further investigations proved why the default AMPT model can describe NSC(3,2) but underestimates SC (3,2). By taking the differences in the individual flow harmonics (v 2 and v 3 ) between the model and data into account, it was possible to recover the difference in SC (3,2) between the data and the model. The discrepancy in SC (3,2) can be explained by the overestimated individual v n values as reported in [100] in all centrality ranges.
In the case of SC(4,2), the string melting configuration of the AMPT model can describe the data fairly well while the default configuration underestimates it. The NSC(4,2) observable is slightly overestimated by the string melting setting which can describe SC(4,2) but the default AMPT configuration can describe the data better. The influence of the hadronic rescattering phase on NSC(4,2) is opposite to other observables (SC(3,2), NSC (3,2) and SC (4,2)). The hadronic rescattering makes NSC(4,2) slightly smaller. It should be noted that the agreement with SC(m,n) should not be overemphasized since there are discrepancies in the individual v n between the AMPT models and the data as was demonstrated for SC (3,2). Hence the simultaneous description of SC(m,n) and NSC(m,n) should give better constraints on the parameters in AMPT models. The string melting AMPT model describes SC (5,3) and NSC (5,3) well. However, the same setting overestimates SC(5,2) and NSC (5,2). The default AMPT model can describe NSC(5,3) and NSC(5,2) fairly well as in the case of NSC (3,2) and NSC (4,2). In the case of SC(4,3), neither of the settings can describe the data but the default AMPT model comes the closest to the data. The NSC(4,3) observable is well described by the default AMPT model but cannot be reproduced by the string melting AMPT model. In summary, the default AMPT model describes well the normalized symmetric cumulants (NSC(m,n)) from lower to higher order harmonic correlations while the string melting AMPT model overestimates NSC (3,2) and NSC(5,2) and predicts a very weak correlation both for NSC (3,2) and NSC(4,3).
As discussed in Sec. 5, a hierarchy NSC(5,3) > NSC(4,2) > NSC(5,2) holds for centrality ranges > 20% within the errors. Except for the 0-10% centrality range, we found that the same hierarchy also holds in the hydrodynamic calculations and the AMPT models explored in this article. While NSC(5,2) is smaller than NSC(5,3), SC(5,2) is larger than SC(5,3). The observed inverse hierarchy, SC(5,2) > SC(5,3), can be explained by different magnitudes of the individual flow harmonics (v 2 > v 3 ). This can be attributed to the fact that flow fluctuations are stronger for v 3 than v 2 [14]. This was claimed in Ref. [58] and also seen in Ref. [101] based on the AMPT model calculations. NSC(m,n) correlators increase with larger η/s in hydrodynamic calculations in the 0-30% centrality range in the same way as the event plane correlations [102,103]. In semi-peripheral collisions (> 40%), the opposite trend is observed.
We list here the important findings from the model comparisons to the centrality dependence of SC(m,n) and NSC(m,n): (i) The NSC (3,2) observable is sensitive mainly to the initial conditions, while the other observables are sensitive to both the initial conditions and the temperature dependence of η/s.
(iii) All the VISH2+1 model calculations with large η/s fail to describe the centrality dependence of the correlations regardless of the initial conditions.
(iv) Among the VISH2+1 model calculations with small η/s (η/s = 0.08), the one with the AMPT initial conditions describes the data qualitatively, but not quantitively for most of the centrality ranges.
(v) The default AMPT model can describe the normalized symmetric cumulants (NSC(m,n)) quantitively for most centralities while the string melting AMPT model fails to describe them.
The agreement of various model calculations with the data is quantified by calculating the χ 2 /N dof : where y i ( f i ) is a measurement (model) value in a centrality bin i. The systematical and statistical errors are combined in quadrature σ i = σ 2 i,stat + σ 2 i,syst . The total number of data samples N dof in Eq. 5 is 4, which corresponds to the number of bins in the centrality range 10-50% used in χ 2 /N dof calculations. The χ 2 /N dof for model calculations which are best in describing the SC observables for each of the three different types of models is shown in Fig. 7.
The results for SC(m,n) and NSC(m,n) are presented in Fig. 7a and 7b, respectively. The χ 2 /N dof values for the individual flow harmonics v n for n=2-4 are shown in Fig. 7c. We found that in case of the calculations from VISH2+1 with AMPT initial conditions (η/s = 0.08) and the default configuration of the AMPT model, the χ 2 /N dof values for SC(m,n) are larger than those for NSC(m,n). This reflects the fact that the individual flow harmonics v n are not well described by those models compared to event-by-event EKRT+viscous hydrodynamics. This is quantified in Fig. 7c, where the χ 2 /N dof values for v n are much larger both for VISH2+1 and default AMPT calculations than event-by-event EKRT+viscous hydrodynamics. The default configuration of the AMPT model gives the best χ 2 /N dof values for NSC(m,n), especially for NSC (3,2). However the χ 2 /N dof values of this model are largest for v n among the models especially for v 2 .
The χ 2 /N dof values for v 2 and v 3 are significantly smaller than those for SC (3,2) and NSC (3,2) for all the hydrodynamic calculations. The χ 2 /N dof values for SC (4,2) and NSC(4,2) from event-by-event EKRT+viscous hydrodynamics are comparable to that for v 2 but larger than for v 4 . The χ 2 /N dof for calculations for v n with constant η/s = 0.20 ("param0") are smaller than those with temperature dependent η/s parameterization with a minimal value of η/s = 0.12 at the temperature around 150 MeV ("param1"), while an opposite trend is observed for SC(m,n), in particular for SC(4,2) and SC (5,3). This illustrates that a combination of the SC(m,n) observables with the individual flow harmonics v n may provide sensitivity to the temperature dependence of the η/s(T ) and together they allow for better constraints of the model parameters.
Even though the calculations from event-by-event EKRT+viscous hydrodynamics give the best χ 2 /N dof values for both SC(m,n) and NSC(m,n), the χ 2 /N dof values are large especially for the observables which include v 3 . Even with the best model calculations, the χ 2 /N dof value varies a lot depending on the model parameters and/or different order SC observables, which implies that the different order harmonic correlations have different sensitivity to the initial conditions and the system properties.

Transverse Momentum Dependence of Correlations between v 2 and v 3 and between v 2 and v 4
The NSC(3,2) and NSC(4,2) observables as a function of p T,min are compared to the AMPT simulations in Fig. 8 and Fig. 9, respectively. The observed p T dependence for NSC (3,2) in mid-central collisions is also seen in AMPT simulations for higher p T,min . The default configuration of the AMPT reproduces NSC (3,2), while the other AMPT configurations predict a very strong p T dependence above 1 GeV/c and cannot describe the magnitudes of both NSC(3,2) and NSC (4,2) simultaneously. In the case of NSC (3,2), the default AMPT model describes the magnitude and p T dependence well in all collision centralities except for 40-50% where the model underestimates the data and shows a stronger p T dependence than the data. As for NSC (4,2), the default AMPT configuration which describes NSC (3,2) can also reproduce the data well except for the 10-20% and 40-50% centralities. Comparison of the string melting AMPT configuration with and without hadronic rescattering suggests that a very strong p T dependence as well as the correlation strength are weakened by the hadronic rescattering. Consequently, the observed weak p T dependence may be due to hadronic rescattering. The relative contributions to the final state particle distributions from partonic and hadronic stages need further study.
The event-by-event EKRT+viscous hydrodynamic calculations are also compared to the data in Fig. 8 and Fig. 9. In the case of NSC (3,2), the hydrodynamic calculations underestimate the magnitude of the data as discussed in Sec. 6.1 and show very weak p T dependence for all centralities. The p T dependence of NSC (3,2) is well captured by the model calculations in all collision centralities except for 40-50% where the data shows stronger p T dependence than the models. The difference between the model calculations with the two different parameterizations of η/s(T ) is very small. As for NSC (4,2), the model calculations overestimate the magnitude of the data in the 5-20% centralities and underestimate it in the centrality range 30-50%. However, the p T dependence is well described by the model calculations in all centrality ranges, while the difference of the model results for the two parameterizations in most centralities is rather small.
The observed moderate p T dependence in mid-central collisions both for NSC (3,2) and NSC(4,2) might be an indication of possible viscous corrections to the equilibrium distribution at hadronic freeze-out, as predicted in [36]. The comparisons to hydrodynamic models can further help to understand the viscous corrections to the momentum distributions at hadronic freeze-out [45,52,[54][55][56].

Summary
In this article, we report the centrality dependence of correlations between the higher order harmonics (v 4 , v 5 ) and the lower order harmonics (v 2 , v 3 ) as well as the transverse momentum dependence of the correlations between v 3 and v 2 and between v 4 and v 2 . The results are presented in terms of the Symmetric Cumulants SC(m,n). It was demonstrated earlier in [57] that SC(m,n) is insensitive to nonflow effects and independent of symmetry plane correlations.
We have found that fluctuations of SC (3,2) and SC (4,3) are anti-correlated in all centralities while fluctuations of SC(4,2), SC (5,2) and SC (5,3) are correlated for all centralities. These measurements were compared to various hydrodynamic model calculations with different initial conditions as well as different parameterizations of the temperature dependence of η/s. It is found that the different order harmonic correlations have different sensitivities to the initial conditions and the system properties. Therefore they have discriminating power in separating the effects of η/s from the initial conditions on the final state particle anisotropies. The comparisons to VISH2+1 calculations show that all the models with large η/s, regardless of the initial conditions, fail to describe the centrality dependence of higher order correlations. Based on the tested model parameters, the data favor small η/s and the AMPT initial conditions.
A quite clear separation of the correlation strength for different initial conditions is observed for these higher order harmonic correlations compared to the lower order. The default configuration of the AMPT model describes well the normalized symmetric cumulants (NSC(m,n)) for most centralities and for most combinations of harmonics which were considered. Finally, we have found that v 3 and v 2 as well as v 4 and v 2 correlations have moderate p T dependence in mid-central collisions. This might be an indication of possible viscous corrections to the equilibrium distribution at hadronic freeze-out. Together with the measurements of individual harmonics the new results for SC(m,n) and NSC(m,n) can be used to further optimize model parameters and put better constraints on the initial conditions and the transport properties of nuclear matter in ultra-relativistic heavy-ion collisions.      [21] B. Schenke, S. Jeon, and C. Gale, "Elliptic and triangular flows in 3 + 1D viscous hydrodynamics with fluctuating initial conditions," J. Phys. G38 (2011) 124169.
[22] P. Bozek  A Model Comparisons of the Individual Flow Harmonics v n As discussed in Sec. 2, NSC(m, n) is expected to be insensitive to the magnitudes of v m and v n but SC(m, n) has contributions from both the correlations between the two different flow harmonics and the individual harmonics v n . Therefore, it is important to check how well the theoretical models used in Sec. 6 describe the measured v n data shown in Sec. 5. v n results presented in this section are for charged particles in the pseudorapidity range |η| < 0.8 and the transverse momentum range 0.2 < p T < 5.0 GeV/c as a function of collision centrality [11]. Results are compared to the event-by-event EKRT+viscous hydrodynamic calculations [45] for two different η/s(T ) parameterizations, labeled in the same way as in [45].
The measured v n for n = 2-4 in Pb-Pb collisions at √ s NN = 2.76 TeV are compared to the event-by-event EKRT+viscous hydrodynamic calculations [45] in Fig. A.1. In these calculations the initial conditions and η/s parameterizations are chosen to reproduce the LHC v n data. The calculations capture the centrality dependence of v n in the central and mid-central collisions within 5% for v 2 and 10% for v 3 and v 4 .
The VISH2+1 calculations with various initial conditions and η/s parameters are compared to the v n data in Fig. A.2. Neither MC-Glauber nor MC-KLN initial conditions can simultaneously describe v 2 , v 3 and v 4 . In particular, for MC-Glauber initial conditions, VISH2+1 with η/s = 0.08 can describe well v 2 from central to mid-central collisions, but overestimates v 3 and v 4 for the same centrality ranges. For MC-KLN initial conditions, VISH2+1 with η/s = 0.20 reproduces v 2 but underestimates v 3 and v 4 for the presented centrality regions. The calculations with AMPT initial conditions improves the simultaneous descriptions of v n (n = 2, 3 and 4). The overall difference to the data is quite large if all the model settings are considered, about 30% for v n (n = 2 and 3) and 50% for v 4 . The calculations with AMPT initial conditions reproduce the observed centrality dependence with an accuracy of 10-20%.   The AMPT calculations with various configurations are compared to the v n data in Fig. A.3. The string melting version of AMPT [78,99] reasonably reproduces v n as shown in Fig. A.3 within 20% for v 2 and 10% for v 3 and v 4 . The version based on the string melting configuration without the hadronic rescattering phase underestimates the data compared to the calculations with the string melting version of AMPT, which demonstrates that a large fraction of the flow is developed during the late hadronic rescattering stage in the string melting version of AMPT. The default version of AMPT underestimates v n for n = 2-4 by ≈ 20%. It should be noted that the default AMPT model can describe the normalized symmetric cumulants (NSC(m,n)) quantitively for most centralities while the string melting AMPT model fails to describe them.
Finally, few selected calculations from three theoretical models which describe the v n data best are shown in Fig. A.4. The calculations from event-by-event EKRT+viscous hydrodynamics, VISH2+1 with AMPT initial conditions (η/s = 0.08) and the string melting version of AMPT give the best description of the individual flow harmonics v n (n = 2, 3 and 4) with an accuracy of 5-20%. The centrality dependence differs in the three models as well as in the different order flow harmonics. Together with SC(m, n) and NSC(m, n), the simultaneous description of individual flow harmonics v n at all orders is necessary to further optimize model parameters and put better constraints on the initial conditions and the transport properties of nuclear matter in ultra-relativistic heavy-ion collisions.  . Results are compared with selected calculations from three different types of models which are best in describing v n coefficients.