Higher harmonic non-linear flow modes of charged hadrons in Pb-Pb collisions at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV

Anisotropic flow coefficients, vn, non-linear flow mode coefficients, χn,mk, and correlations among different symmetry planes, ρn,mk are measured in Pb-Pb collisions at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV. Results obtained with multi-particle correlations are reported for the transverse momentum interval 0.2 < pT< 5.0 GeV/c within the pseudorapidity interval 0.4 < |η| < 0.8 as a function of collision centrality. The vn coefficients and χn,mk and ρn,mk are presented up to the ninth and seventh harmonic order, respectively. Calculations suggest that the correlations measured in different symmetry planes and the non-linear flow mode coefficients are dependent on the shear and bulk viscosity to entropy ratios of the medium created in heavy-ion collisions. The comparison between these measurements and those at lower energies and calculations from hydrodynamic models places strong constraints on the initial conditions and transport properties of the system.


Introduction
One of the primary goals in the ultra-relativistic heavy-ion collision programs at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) is to study the nuclear matter at extreme conditions. The pressure gradients in the strongly interacting matter, known as the Quark-Gluon Plasma (QGP), are believed to drive the hydrodynamic expansion observed through anisotropy in multi-particle correlations in high energy collisions at RHIC and the LHC [1,2]. The anisotropic expansion of the medium, commonly referred to as anisotropic flow [1], can be characterized by a Fourier decomposition of the azimuthal particle distribution with respect to the reaction plane [3,4] dN dϕ ∝ 1 + 2 ∞ n=1 v n cos(n(ϕ − ψ RP )), (1.1) where the flow coefficient v n is the magnitude of the n-th order flow, and the reaction plane ψ RP defined by the beam direction and impact parameter which is defined as the distance between the centers of two colliding nuclei. Due to fluctuations in the initial state energy density profile, it is useful to define symmetry planes of different orders, where the n-th order plane ψ n defines the orientation of the n-th order complex flow vector V n ≡ v n e inψn . The expansion of the azimuthal distribution about ψ n also yields finite values of odd coefficients [5,6]. Anisotropic flow measurements through two-and multi-particle azimuthal correlations [6][7][8][9][10][11][12][13] have provided important information on the medium response and in particular its transport coefficients such as the shear viscosity to entropy density ratio (η/s), bulk viscosity to entropy density ratio (ζ/s) and the equation of state [14]. Studies have shown the relativistic hydrodynamic nature of the medium [1,2,[15][16][17][18][19][20][21][22], with η/s close to the AdS/CFT minimum 1/(4π) [23].

JHEP05(2020)085
The initial state eccentricity, determined from the energy density profile, is obtained from the definition [5] ε n e inΦn = −{r n e inϕ }/{r n }, n ≥ 2, (1.2) 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 the azimuthal angle, e(x, y, τ 0 ) is the energy density at the initial time τ 0 , and Φ n is the participant plane angle, defining the spatial symmetry of the nuclear constituents in the participant region (see refs. [24,25]). Hydrodynamic models demonstrate that the second and the third harmonic flow coefficients exhibit an almost linear dependence on the initial eccentricity coefficients ε n [26]. Considering that the anisotropic expansion is a result of a hydrodynamic evolution governed by η/s, a measurement of the second and third harmonics combined with hydrodynamic calculations can constrain the properties of the medium. Several estimates for the limits of η/s were determined through measurements of elliptic flow coefficient v 2 [27 -32] and their comparison with hydrodynamic calculations. Consequently, the early constraints placed the value of η/s between 0.08 to 0.16 [33][34][35]. However, the limited sensitivity of the elliptic flow to η/s and the large uncertainty in the initial state anisotropy inhibit a precise determination of the value of η/s [34,[36][37][38], and its temperature dependence, which was recently shown to be explorable during the second run of LHC [39,40]. In addition, part of the anisotropic flow can also originate from the hadronic phase [41][42][43]. It has been shown in [43,44] that the inclusion of the temperature dependent bulk viscosity ζ/s(T) in hydrodynamic simulation lead to a better description of the average transverse momentum of charged hadrons and on the elliptic flow coefficient.
The effects of bulk viscosity should be considered when extracting any transport coefficient from the data [45][46][47].
Flow harmonics of order n ≥ 3 reveal finer details of initial conditions [6,8,10,11,13], enabling to constrain η/s better [39,40,48,49]. Higher flow harmonics n > 3 do not exhibit a linear response to the initial anisotropy [26] as a finite contribution is induced by the initial state anisotropy of the lower orders [50,51]. For example, the fourth order flow vector V 4 gains contributions not only from the fourth order flow (linear flow mode), but also from the second order flow (non-linear flow mode). Starting from the V n estimators studied in [50], the flow can be expressed as a vector sum of the linear and non-linear modes where χ n,mk is called non-linear flow mode coefficient, characterizing the non-linear flow mode induced by the lower order harmonics. The high order linear component is denoted by V nL , while the many higher order linear couplings are depicted by E(. . . ) for V 8 . The -2 -

JHEP05(2020)085
V nL is linearly related to a cumulant-defined anisotropy [52] ε 4 e i4Φ 4 = ε 4 e i4Φ 4 + 3 r 2 2 r 4 ε 2 e i4Φ 2 (1.4) as opposed to the relation v n ∝ ε n , where v n is the magnitude of the total contribution and ε n is given by eq. (1.2). In earlier measurements performed by ALICE [53], the non-linear flow mode coefficients were measured up to the sixth harmonic order in Pb-Pb collisions at √ s NN = 2.76 TeV. It was indicated that the coefficients χ 5,23 and χ 6,33 are sensitive not only to η/s, but also to the distinctive energy density profiles generated by different initial conditions. It was observed that the hydrodynamic models with their respective initial conditions Monte-Carlo (MC)-Glauber [54,55], MC-KLN [33,56], and IP-Glasma [57]), are unable to reproduce these measurements, which indicates that the model tuning and η/s parameterization require further work. In this paper, measurements of high order flow coefficients in Pb-Pb collisions at √ s NN = 5.02 TeV are presented. The flow coefficients v n are measured up to the ninth harmonic, v 9 , extending the previous measurements of v 2 -v 6 [58]. The data recorded during the 2015 heavy-ion run of the LHC allow the measurements of non-linear flow mode and correlations between symmetry planes to be extended. A total of six non-linear flow mode coefficients are measured, including the non-linear flow mode coefficient χ 7,223 , for which the sensitivity to η/s is expected to be significantly stronger than for the lower oddharmonic coefficient χ 5,23 [37,59]. The results are compared with those in Pb-Pb collisions at √ s NN = 2.76 TeV [53] and various state of the art hydrodynamical calculations.

Formalism and observables
In order to separate the linear and non-linear contributions from eq. (1.3), one assumes the respective contributions to be uncorrelated [60]. For example for the fourth order V 4 , by mean-squaring the equations in eq. (1.3) and setting (V Here denotes an average over all events and * the complex conjugate. The magnitudes of the linear and non-linear flow coefficients are denoted with v 4L and v 4,NL , respectively. The observables of the non-linear response mode are constructed from the projections of flow vectors on to the symmetry planes of lower harmonics [61,62]. For n = 4, the magnitude of the non-linear response mode is given by where v 2 4,22 ≡ v 2 4,NL ≡ χ 2 4,22 |V 2 | 4 . The right-hand side approximation holds if the low (n = 2, 3) and high order flow is weakly correlated. Only the fourth harmonic is shown here and the complete list of other harmonics are provided in appendix A.

JHEP05(2020)085
The contributions from short-range correlations unrelated to the common symmetry plane, commonly referred to as "non-flow", are suppressed by using the subevent method where the event is divided into two subevents separated by a pseudorapidity gap [4]. The underlying multi-particle correlation coefficient for subevent 1/2 as determined using eq. (2.2), 1 and a similar treatment is applied for the subevent B, for which v B 4,22 is obtained by swapping B for A in the aforegiven expression. The final result is then the average of the results from subevents A and B.
The symmetry-plane correlations are defined as the ratio between the magnitude of the non-linear flow modes and flow coefficients [63]. For n = 4, one obtains A value close to zero indicates weakly correlated symmetry planes, while a value reaching one implies a strong correlation. The correlations between symmetry planes reflect those of the corresponding initial state participant planes [53,64], therefore providing valuable information on the evolution of the QGP. Correlations between symmetry planes have been previously studied using the event-plane method [64,65], event plane describing an experimentally approximated symmetry plane. However, these results depend on the event-plane resolution [66], which complicates the comparison between data and theoretical calculations. Recently, the ALICE Collaboration has measured symmetry-plane correlations [53]. It was found that the correlations of symmetry planes of higher harmonics with second and third order symmetry planes increased for less central collisions. Furthermore, the comparison with hydrodynamic calculations revealed the importance of final-state collective dynamics in addition to the initial-state density fluctuations [33] as it is known that the observation of correlated final state symmetry planes implies the existence of fluctuations in the initial state eccentricity vectors. The fourth non-linear flow mode coefficient, with the aforementioned assumptions, is given by [59] (2.4)

Experimental setup and data analysis
The data sample consists of about 42 million minimum bias Pb-Pb collisions at √ s NN = 5.02 TeV, recorded by ALICE [67,68] during the 2015 heavy-ion run at the LHC. Detailed descriptions of the detector can be found in [67,69,70]. The trigger plus crossing of beam is provided by signals from the two scintillator arrays, V0A and V0C [67,71], covering the pseudorapidity intervals 2.8 < η < 5.1 and −3.7 < η < −1.7, respectively. A primary vertex position less than 10 cm in beam direction from the nominal interaction point is required. Pile-up events are removed by correlating the V0 multiplicity with the multiplicity from the first Silicon Pixel Detector (SPD) [67,72] layer. To further remove pile-up events,

JHEP05(2020)085
information from the Time-of-Flight (TOF) [73] detector is used: the multiplicity estimates from the SPD are correlated with those imposed with a TOF readout requirement. The centrality of the collision is determined using information from the V0 arrays. Further details on the centrality determination in ALICE are given in [74]. Only events in the centrality range 0% to 60% are used in the analysis. The track reconstruction is based on combined information from the Time Projection Chamber (TPC) [67,75] and the Inner Tracking System (ITS) [67,72]. To avoid contributions from secondary particles, the tracks are required to have a distance of closest approach to the primary vertex of less than 3.2 cm and 2.4 cm in the longitudinal and transverse directions, respectively. Such a loose Distance of Closest Approach (DCA) track cut is chosen to improve the uniformity of the ϕ-distribution for the Q n -vector calculation. Furthermore, each track is required to have at least 70 TPC space points out of the maximum 159, and the average χ 2 per degree of freedom of the track fit to the TPC space points to be less than 2. Minimum 2 hits are required in the ITS. In order to counteract the effects of track reconstruction efficiency and contamination from secondary particles [76], a HIJING simulation [77,78] with GEANT3 [79] detector model is employed to construct a p T -dependent track weighting correction. The track reconstruction efficiency is approximately 65% at p T = 0.2 GeV/c and 80% at p T > 1.0 GeV/c, while the contamination from secondaries is less than 10% and 5%, respectively. Only particle tracks within the transverse momentum interval 0.2 < p T < 5.0 GeV/c and pseudorapidity range 0.4 < |η| < 0.8 are considered. A pseudorapidity gap |∆η| > 0.8 is used to suppress the non-flow. The observables in this analysis are measured with multi-particle correlations obtained using the generic framework for anisotropic flow analysis [80].

Systematic uncertainties
The systematic uncertainties are estimated by varying criteria for selecting the events and tracks. The systematic evaluation is done by independently varying the selection criteria, and the results given by this variation are then compared to the default criteria given in section 3. The total uncertainty is obtained by assuming that the individual sources are uncorrelated, which are then quadratically summed.
Summaries of the relative systematic uncertainties are given in tables 1-4. Uncertainties stemming from the event selection criteria are estimated by changing the rejection based on the vertex position from 10 cm to 8 cm and by adjusting the pile-up rejection criteria. It is found that the contribution to the uncertainty is generally negligible, below 1%. An alternative centrality determination is employed using the event multiplicity estimates from the SPD layers. The uncertainty related to the centrality determination is less than 2% for all observables, except for v 7 to v 9 for which the uncertainty increases to 10%.
The ALICE detector can be operated with either positive or negative solenoid magnetic field polarity. The polarity of the field affects the direction of the charged particle curvature, while also subjecting the structural materials of the detector itself to either a positive or negative magnetic field. The default data set is composed of events recorded with both polarities. The results produced with exclusively either negative or positive magnetic field configurations deviate from the default by up to 15% in case of flow coefficients, and 28% for Event Selection z-vertex cut < 0. of the higher interaction rates is estimated by comparing results from different regions of the TPC, one for η > 0 and the other η < 0. The maximum uncertainty is evaluated less than 15%. The track reconstruction related uncertainty, referred to as tracking mode, is evaluated by comparing the results obtained with tracks for which the requirement for the number of hits in the ITS layers is changed. In this case, the uncertainty is generally less than 15%, and a maximum 20% is evaluated for ρ 7,223 . Furthermore, the track selection criteria is tightened by increasing the minimum number of the TPC space points from 70 to 90, resulting in uncertainties around 1% to 3%.

Results
In this section, the measurements of the flow coefficients, the non-linear modes, symmetryplane corre-lations and the non-linear flow mode coefficients are presented. They are compared with hydrodynamic calculations with various settings [25,57,81,82]. The first calculation is based on an event-by-event viscous hydrodynamic model with EKRT initial conditions [25,81] using a value of η/s = 0.20 (param0) and a temperature dependent η/s(T ) (param1). For both configurations, ζ/s is set to zero. The visualization of the model parameters can be found in figure 1. The second calculation employs the iEBE-VISHNU hybrid model [83] with AMPT [84][85][86] and T R ENTo [87] initial conditions. The η/s = 0.08 and ζ/s = 0 are taken for param2, while the η/s(T ) and ζ/s(T ) (param3), extracted using Bayesian analysis [45] (except for the normalization factors) from a fit to the final multiplicities of the charged hadron spectra in Pb-Pb collisions at √ s NN = 5.02, are used for the T R ENTo initial conditions. The third calculation uses the MUSIC model [88] with IP-Glasma [89] initial conditions with a value of η/s = 0.095 and ζ/s(T ) (param4). Each of the η/s(T ) parameterizations is adjusted to reproduce the measured charged hadron multiplicity, the low-p T region of the charged-hadron spectra, and v n from central to midperipheral collisions up to the fourth harmonic at RHIC and the LHC [25,44,57,84,[90][91][92].
The model configurations are summarized in table 5.
In figure 2, the measurements of the flow coefficients from v 2 to v 9 are presented. The first two coefficients up to v 6 have been extensively measured at RHIC and LHC [6][7][8][9][10][11][12][13], and more recently also v 7 [49]. The present analysis reports the first results on higher harmonic coefficients from v 7 to v 9 in ALICE, where v 8 and v 9 are measured for the first time at the LHC energies. The coefficients exhibit a modest centrality dependence, and their magnitude is similar to that of v 7 within statistical and systematic uncertainties. The measurements up to v 6 are compatible with those published previously [58]. Figure 2 also shows the comparison between the measured v n and model calculations. The hydrodynamic calculations qualitatively reproduce the v n measurements, and the overall model depiction is very good for v 2 and v 3 . For n ≥ 4 however, the calculations show noticeable overestimations, especially in mid-peripheral collisions. For v 5 and v 6 , the data are well described by EKRT+param0, showing a better agreement than the temperature dependent EKRT+param1. The data are also described by AMPT+param2, for which the agreement for v 5 and v 6 is good in mid-central and mid-peripheral collisions. IP-Glasma+ param4 and T R ENTo+param3 overestimate the measurements by a factor of 1.5∼2. Values of v 7 are well estimated by AMPT+param2, and v 8 by both AMPT+param2 and T R ENTo+param3 within uncertainties. In other cases, the data are overestimated by the other models.
To study the dependence on the harmonic order of the anisotropy coefficients [97], figure 3 shows values of different coefficients as a function of n for all centralities. This presentation is particularly well suited in visualizing the harmonic dependence, and the acoustic scaling [97] observed across the harmonic orders. The decrease in v n with increasing harmonic order up to n = 7 indicates viscous damping [97]. This means that the higher frequency waveform propagating through the medium should get more damped is the same for param3 and param4 and taken from eq. 5 in [45] motivated by refs. [43,[93][94][95]. For the parameters with T R ENTo initial condition, the ones based on identified yields are taken from table 4 in [45]. The ζ/s normalization factor used with IP-Glasma (T R ENTo) initial conditions is 0.9 (1.25). The models with ζ/s = 0 are not shown on the right.  until freeze-out takes place. In [98,99] the viscosity effect is explained as the main contributor to the observed damping. It is speculated, that another driving factor is the phase of the oscillation itself, which also contributes to the magnitude at the time of freeze-out. The measurements show that there is a hint of v 9 > v 8 , as also predicted in the acoustic model [97].  Figure 3. v n as a function of the harmonic order n for various centrality intervals. figure 4 show the flow modes of v 6 and v 7 . Only the non-linear flow modes of v 6 and v 7 are presented. The v 6,222 increases from zero to approximately half of the total v 6 in mid-central collisions, while the other mode, v 6,33 , has a much weaker centrality dependence. The relatively large magnitude of these flow modes imply strong contributions from the second and third order harmonics. Finally, v 6,24 follows the trend of the total magnitude. The magnitude of v 6,24 comes close to the total, which in turn suggests not only strong contributions from the second harmonic order, but also the fourth one. The v 6,24 induced by the fourth order is seen to be the dominant contribution to the sixth order from 10% centrality classes and higher. For the seventh order v 7 , there are three non-linear contributions, of which v 7,223 is measured. The centrality dependence is similar as with the v 6 coefficient, and there is a similar general trend as for the lower order harmonics among the non-linear flow modes.

Panels (c) and (d) of
The coefficients ρ n,mk , quantifying the correlations amongst different symmetry planes, are presented as a function of centrality in figure 5. Except for ρ 6,33 , all coefficients indicate an increase in correlation between symmetry planes with increasing centrality class of the collision. The measurements generally agree with the ones obtained at the lower energy. The ρ 6,222 is the only coefficient for which an energy dependence can be observed. The hydrodynamic calculations reproduce the measurements within the large theoretical uncertainties. For ρ 4,22 , ρ 5,23 , and ρ 6,222 , T R ENTo+param3 however underestimates the data in mid-central collisions. in [53]. For χ 4,22 and χ 5,23 , the centrality dependence and overall magnitude agree well with the results from the lower beam energy. The centrality dependence of the new data is similar to the previous results: a larger value in more central collisions, decreasing close to unity towards 50% centrality.
All of the non-linear flow mode coefficients for the sixth harmonic agree with the previous measurements. The centrality dependence of χ 6,222 is similar to the ones of the lower order coefficients, and the overall magnitude similar to χ 4,22 . As for χ 6,33 , no clear centrality dependence is observed within the current experimental uncertainties. Whereas the previous measurements are unable to distinguish between the magnitudes of χ 6,222 and χ 6,33 , the current results show that χ 6,222 > χ 6,33 across the whole centrality interval. For χ 7,223 , the overall magnitude is larger than for the other non-linear flow mode coefficients.
The hydrodynamic calculations for the non-linear flow mode coefficients show slightly more variation compared to the symmetry-plane correlations. As seen from the panels of figure 6, one observes the reproduction of the data points by EKRT+param0 up to  the modes of the sixth harmonic, and T R ENTo+param3 in all harmonics. The EKRT+ param1 calculations slightly overestimate the centrality dependence of the non-linear flow mode coefficients. It can be seen that the parameterizations of the EKRT presented here imply χ n,mk across all harmonic orders to have sensitivity to η/s, whereas in the previous calculations in [53], weak η/s dependence was found for χ 4,22 and χ 6,222 . The fifth order coefficient χ 5,23 is expected to be quite sensitive to η/s in central collisions as can be seen from the difference of the predicted values from EKRT+param0 and EKRT+param1. The AMPT+param2 calculations underestimate the magnitude of some of the measured nonlinear flow mode coefficients in various centrality classes, especially χ 5,23 , χ 7,223 as well as χ 6,24 . The IP-Glasma+param4 calculations overestimate the measurements in some centrality intervals. The agreement between data and the model calculations is quantified by calculating the reduced χ 2 /N dof defined as where y i and f i are the values for data and calculations, respectively, and σ 2 i = σ 2 i,stat + σ 2 i,syst + σ 2 f i ,stat is the quadratic uncertainty in terms of statistical measurement σ i,stat , model uncertainties σ f i ,stat , and systematic uncertainties σ i,syst in centrality bin i. Here N dof represents the number of data points across the centrality interval.
The χ 2 /N dof for the flow coefficients are presented in figure 7, panel (a). It is observed that IP-Glasma+param4 gives the best description of v 2 and v 3 compared to the other models, indicated by the overall low value of χ 2 /N dof . However, the overall performance -14 -  Figure 7. Overview of various model comparisons with data, quantified by χ 2 /N dof . Lower χ 2 /N dof represents a better overall description for a given observable.
of IP-Glasma+param4 is considerably worse at n ≥ 4, for which the data are overestimated, as seen in figure 2. For v 4 to v 6 , EKRT+param0 gives the lowest value of χ 2 /N dof . In the case of EKRT+param1, the χ 2 /N dof is slightly worse than EKRT+param0. The χ 2 /N dof of T R ENTo+param3 is very close to that of IP-Glasma+param4, indicating a comparable description of data between the two model configurations. At low harmonic orders, T R ENTo+param3 performs slightly worse than IP-Glasma+param4. For n ≥ 4, description of the data between these two models are comparable except for n = 8, where T R ENTo+param3 clearly has a better magnitude and centrality depiction. Notably this can be seen for v 8 where the χ 2 /N dof value is the lowest across all models. Finally, the performance of AMPT+param2 can be considered good within the reported χ 2 /N dof values. It is noted that the magnitude of v 7 is best depicted by AMPT+param2 amongst the three models used.
The performance of the models with respect to the symmetry-plane correlations is quantified in panel (b) of figure 7. IP-Glasma+param4 has by far the best estimates of ρ n,mk for ρ 4,22 and ρ 5,23 . For other models, the model depiction is comparable. In low harmonic orders, EKRT+param0 shows good agreement with the data, as well as AMPT+ param2, which has the best agreement in higher harmonic orders. For T R ENTo+param3, the agreement is slightly worse for ρ 5,23 and ρ 6,222 .
The panel (c) of figure 7 shows the χ 2 /N dof for non-linear flow mode coefficients. As seen also in figure 6, T R ENTo+param3 consistently provides the most successful overall description of the data. For other models the data are more frequently over-or underestimated. T R ENTo+param3 estimates χ n,mk better than it does the v n coefficients, for which significant overestimation was present at almost every harmonic order (see figure 2). For EKRT+param0 the agreement is good, but the calculation over-or underestimates in some cases especially in most central or mid-peripheral collisions. Most of the observables are better described by the calculations using EKRT+param0 with a const η/s = 0.2 as compared to results from EKRT+param1 which uses a temperature dependent η/s value.

JHEP05(2020)085
AMPT+param2 performs worse for low-order harmonics as it overpredicts the data in central and mid-central collisions. Of the five configurations, IP-Glasma+param4 describes the data worse in all harmonic orders.
The deviation of the calculated results from the measured value of each observable is of the same order of magnitude for the different models. Even where the model results show gross agreement with overall features in data, the values of χ 2 /N dof vary considerably from one harmonic order to another. Considering the χ 2 /N dof to be a goodness-of-fit estimate to validate any model, these variations suggest that the sensitivity of the different observables on the initial conditions, η/s, and ζ/s are reflected differently in the model calculations. Since the current uncertainties in the model calculations are large for higher order harmonics, the absolute χ 2 test should not be over-interpreted. Both, improved statistical uncertainties in the model calculations and different values of input parameters would be beneficial. However, large sets of calculations in many parameter spaces require substantial computing power. In order to constrain the model parameters Bayesian analysis can provide a plausible approach as demonstrated in [45,47]. At present it is limited to low harmonic-order observables, and the extracted parameters have large uncertainties. Extending the Bayesian analysis to include the results in this paper will help reduce the uncertainties of the model parameters.

Summary
The measurements of anisotropic flow coefficients (v n ), non-linear flow mode coefficients (χ n,mk ), and correlations among different symmetry planes (ρ n,mk ) in Pb-Pb collisions at √ s NN = 5.02 TeV are presented. The anisotropic flow coefficients are measured up to v 9 , where v 8 and v 9 are measured for the first time at LHC energies. It is observed that v n decreases as n increases, observing n-ordered damping up to n = 7. The v n is found to be enhanced for n > 7. The non-linear contribution becomes dominant towards peripheral collisions in all harmonic orders. The strength of the non-linear flow mode and the symmetry-plane correlations depends also on harmonic orders. The non-linear flow mode coefficients show a clear centrality and harmonic order dependencies and the strongest non-linear mode coefficients is observed for the fifth and seventh harmonic orders. These results are compared with various hydrodynamic model calculations with different initial conditions, as well as different parameterizations of η/s and ζ/s. None of the models presented in this paper simultaneously describe the v n coefficients, χ n,mk , or ρ n,mk . Based on the model and data comparisons, among all the models, the event-byevent viscous hydrodynamic model with EKRT initial conditions and a constant η/s = 0.2 is observed to describe the data best, as far as the harmonics up to the sixth order are concerned. As a result further tuning is required to find the accurate parameterization of η/s and ζ/s. It is found that the different order harmonic observables presented in this paper have different sensitivities to the initial conditions and the system properties. These results allow further model parameters to be optimized and the initial conditions and the transport properties of nuclear matter in ultra-relativistic heavy-ion collisions to be better constrained.

A List of observables
In this section the complete list of the measured observables is presented. By root-meansquaring the equations in eq. (1.3), one obtains a starting point for the definitions presented in this section. Provided that the linear and non-linear parts are uncorrelated, the following harmonic projections are obtained The higher order superpositions in eq. (1.3) include the coupling constants for the higher order linear responses. In a more complete treatment [100], the extraction of the higher order non-linear flow mode coefficients are performed by correlating the corresponding superpositions with those of the relevant harmonics, effectively resulting in a more general projection. The results agree with the expressions in eq. (2.4), and generate additional high order linear coupling coefficients .

(A.4)
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.