Higher-order correlations between different moments of two flow amplitudes in Pb-Pb collisions at sNN=5.02 TeV

The correlations between different moments of two ﬂow amplitudes, extracted with the recently developed asymmetric cumulants, are measured in Pb-Pb collisions at √ s NN = 5 . 02 TeV recorded by the ALICE detector at the CERN Large Hadron Collider. The magnitudes of the measured observables show a dependence on the different moments as well as on the collision centrality, indicating the presence of nonlinear response in all even moments up to the eighth. Furthermore, the higher-order asymmetric cumulants show different signatures than the symmetric and lower-order asymmetric cumulants. Comparisons with state-of-the-art event generators using two different parametrizations obtained from Bayesian optimization show differences between data and simulations in many of the studied observables, indicating a need for further tuning of the models behind those event generators. These results provide new and independent constraints on the initial conditions and transport properties of the system created in heavy-ion collisions.


Introduction
An important ongoing program in high-energy nuclear physics is the exploration of the phase diagram of quantum chromodynamics (QCD) at large values of temperature and/or density.In that regime, quarkgluon plasma (QGP), an extreme state of nuclear matter in which quarks and gluons are deconfined, is formed.Ultrarelativistic heavy-ion collisions at the BNL Relativistic Heavy Ion Collider (RHIC) and the CERN Large Hadron Collider (LHC) recreate the conditions in which the QGP can be produced and, consequently, studied.Properties of the QGP are typically portrayed by specifying its transport properties (e.g., shear or bulk viscosity), or its equation of state.To date, the QGP remains the state of matter with the smallest ratio of shear viscosity to entropy density (η/s) ever discovered.Besides scrutinizing the details of individual states of nuclear matter in the QCD phase diagram, it is equally important to constrain the nature of the transition between different states, i.e., whether it is a smooth crossover, first-or second-order phase transition, etc.Another extensive effort is allocated to the search for the critical point in the QCD phase diagram, whose existence is still not confirmed experimentally.Recent reviews can be found in Refs.[1][2][3][4].
In non-central heavy-ion collisions, the measured correlations in momentum space among detected particles are dominated by correlations originating from collective anisotropic flow [5].Since anisotropic flow links the initial-state anisotropies in the collision geometry to the final-state anisotropies in momentum space, its measurements can be used to map out the whole history of the heavy-ion evolution.In particular, anisotropic flow is a sensitive probe of shear viscosity of the QGP, whose collective expansion is successfully described by its hydrodynamic response to the anisotropies in the initial-state geometry.By using a Fourier series decomposition of the azimuthal distribution f (ϕ) of reconstructed particles, anisotropic flow is quantified with flow magnitudes v n and symmetry planes Ψ n [6] Due to random fluctuations of the impact parameter vector (the vector connecting the centers of the two colliding ions) from one collision to another, this series cannot be obtained directly from the data using the simple relation v n e inΨ n = e inϕ .Instead, multiparticle azimuthal correlations, which have additional advantages and are not influenced by these fluctuations for suitably chosen harmonics n i , are utilized as experimental estimators according to the following analytic expression [7] e i(n In the above equation, angular brackets indicate an average over all distinct sets of k azimuthal angles of the tracks obtained from the same event.Further technical details about multiparticle azimuthal correlations can be found in Refs.[8][9][10][11]. With the multiparticle correlation techniques, the correlations between different flow amplitudes v n or symmetry planes Ψ n can be used to extract independent information about all stages in the heavy-ion evolution, when compared to the traditional measurements of individual flow amplitudes v n [12].Among them, the symmetric cumulants (SCs) were designed to study the genuine correlations (or cumulants) between the same lowest-order moments of two different flow amplitudes v 2 m and v 2 n .The simplest realization of SC observables in terms of flow amplitudes and the corresponding experimental estimators built from multiparticle azimuthal correlations is given by [11,13,14] SC(m, n) ≡ ⟨v 2 m v 2 n ⟩ − ⟨v 2 m ⟩⟨v 2 n ⟩ = ⟨⟨e i(mϕ 1 +nϕ 2 −mϕ 3 −nϕ 4 ) ⟩⟩ − ⟨⟨e im(ϕ 1 −ϕ 2 ) ⟩⟩⟨⟨e in(ϕ 1 −ϕ 2 ) ⟩⟩. (3) In the first line above, ⟨•⟩ indicates an average over all events.Meanwhile, ⟨⟨•⟩⟩ in the second line indicates that the averaging is performed in two steps: first over all distinct particle pairs or quadruplets in an event, and then in the second step the all-event averages are obtained from these single-event averages by weighting them with 'number of combinations' multiplicity weights M(M − 1) and M(M − 1)(M − 2)(M − 3) for two-and four-particle correlations, respectively [10].Furthermore, defined this way, SC observables are multivariate cumulants of flow amplitudes squared, and therefore they obey all nontrivial mathematical properties of multivariate cumulants for stochastic variables v 2 n by definition [15].The azimuthal correlators are then used only via Eq.( 2) to estimate each individual term in the cumulant expansion [16], i.e., the cumulant expansion is not performed directly on the azimuthal angles.
The particular importance of SCs in flow studies is that they bring new information on QGP properties, as shown in Ref. [13].From experimental measurements of various combinations of SC(m, n) and their comparison to theoretical predictions, it was concluded that SCs are sensitive to the temperature dependence of η/s.On the other hand, the analyses of single flow amplitudes can only probe the average value ⟨η/s(T )⟩ [13,14,17].Motivated by this access to new information on the system properties, studies were conducted to generalize the SCs from two to any number of harmonics using the cumulant formalism [18].The resulting three-harmonic SCs were recently measured by the ALICE Collaboration in Pb-Pb collisions recorded in 2010 [19] and 2015 [17].
In parallel to the measurements of new flow observables, techniques like Bayesian optimization, where the model parameters are systematically varied to find the best fit to all selected experimental measurements, were applied.In recent years, multiple Bayesian studies were conducted to constrain the QGP transport properties and details of the initial conditions in Pb-Pb collisions [12,[20][21][22][23][24][25][26][27].While much information on Pb-Pb collisions was obtained through Bayesian analyses, the uncertainties on the extracted parameters are still quite large.As demonstrated in recent studies [12,27], this problem can partially be solved by increasing the scope of the input experimental data, either from various beamenergy ranges or from additional independent observables sensitive to the initial conditions and system properties.In recent publications [12,27], dedicated sensitivity studies were also performed on sophisticated higher-harmonic observables for the first time, and it was concluded that the higher-order flow observables exhibit the highest sensitivity to the QGP transport properties.
The success of the SC observables in providing crucial input to Bayesian analyses has prompted a natural quest towards further generalizing the SCs to a new and independent set of multiharmonic observables involving different moments of the flow amplitudes without any influence from the symmetry planes.Among examples of such estimators are the asymmetric cumulants (ACs) [16].As they probe the genuine correlations between different moments of the flow amplitudes, they can, by definition, extract new information (about event-by-event flow fluctuations, non-linear response, etc.) which is inaccessible to the SCs themselves.
The structure of this article is as follows.Section 2 reviews the theoretical definitions of ACs and the corresponding experimental estimators.The technical details of the data analysis are introduced in Sec. 3. In Sec. 4, the first experimental measurements of the ACs and their comparisons with state-of-the-art model predictions are presented.The article is summarized in Sec. 5, while additional support material can be found in the Appendix.

Definitions and scope
This article presents the experimental measurements of the AC observables in their simplest realization which involves only two flow amplitudes, evaluated for three different pairs of harmonics chosen to be the same as the ones studied in Ref. [12].For a given pair of different harmonics (m, n), measuring the genuine correlations between v 2a m and v 2b n with AC observables (where a and b are non-equal positive integers) leads to independent information, due to the uniqueness of each multivariate cumulant.Therefore, this analysis focuses on the following new AC a,b (m, n) flow observables: AC 2,1 (m, n), AC 3,1 (m, n), and AC 4,1 (m, n), to study correlations of the second, third, and fourth moments of v 2 m with the first moment of v 2 n .Two remarks have to be made: the lowest-order AC, AC 1,1 (m, n), corresponds to the well-known SC(m, n) (Eq.( 3)), showing that ACs are a natural generalization of SCs.Furthermore, AC observables are invariant under a simultaneous permutation of the indices and the harmonics, e.g., AC a,b (m, n) = AC b,a (n, m).These properties will be exploited in Sec. 4 to discuss the comparisons between AC a,b (m, n) and AC a,b (n, m).
The expressions for the ACs in terms of flow amplitudes are derived in Ref. [16] and are recalled as Due to their cumulant nature, much information can already be extracted from the SCs and ACs themselves.One example is the sign of the measured observables.In Ref. [13], it was shown how the event-by-event flow fluctuation patterns and the signatures of SC observables are related.This idea was also applied in the interpretation of the three-harmonic SC presented in Ref. [19].For ACs, however, the defining equations are more involved and the direct physical interpretation of their signature, positive or negative, is not straightforward, even in the simplest case of AC 2,1 .Another example is the observation of zero magnitude, indicating either an absence of genuine correlations between all involved variables or a trivial consequence of underlying symmetries.In the case of cumulants, however, this has to be interpreted carefully because in the absence of underlying symmetries either all the cumulants at given order are identically zero, or none of them are identically zero [15].Therefore, the nonvanishing results for SC (2,4), obtained previously in Refs.[13,14], imply that all of AC a,b (2, 4) must be nonvanishing as well, and the fact that some of them are consistent with zero would indicate only the limited sensitivity in the available datasets.In addition, a multivariate cumulant is nonvanishing if and only if all variables in its definition are correlated to each other [15].
As was done for the SCs, the ACs can also be normalized to remove the sensitivity of the observable to the strength of the individual harmonics and extract only the contribution of the genuine correlations between the various moments of v 2 m and v 2 n .This is done according to the standard procedure from Refs.[16,28], using Equations ( 4) to ( 6) cannot be used experimentally, as only the azimuthal angles of the final-state particles can be measured.Reliable experimental estimators for the various ACs were presented in Appendix C.2. of Ref. [16].As an example, the expression corresponding to Eq. ( 4) is Each of the azimuthal correlators in the above expression is isotropic (i.e., invariant under the transformation ϕ → ϕ + α where α is an arbitrary angle), which ensures that ACs are not affected by random event-by-event fluctuations of the impact parameter vector.In addition, in each azimuthal correlator, each harmonic appears an equal number of times with positive and negative signs, which ensures that any dependence on symmetry planes is canceled out exactly (see Eq. ( 2)).

Data analysis
The analysis is performed with data from Pb-Pb collisions at the center-of-mass energy per nucleon pair √ s NN = 5.02 TeV recorded by the ALICE detector in 2015 during the LHC Run 2 period.
The detector closest to the beam pipe is the Inner Tracking System (ITS) [34,35], which, besides the TPC, is also used for track reconstruction.It covers the full range in azimuthal angle and consists of six silicon layers, with three different detector technologies.The two innermost layers, the Silicon Pixel Detectors (SPD), have a high spatial granularity making them ideal for reconstructing both primary and secondary vertices.
Minimum bias (MB) events are triggered by a coincident signal in both the V0A and V0C.The selected MB Pb-Pb events all have a reconstructed primary vertex within ±8.0 cm from the nominal interaction point along the beam direction.Background events like beam-gas collisions and pileup were removed using the correlation between the information from the V0 detector and the SPD, as well as the correlation between the number of reconstructed tracks and the information from the ITS and TPC.A further reduction of the pileup events was obtained using the correlation between the information provided by the SPD and by the Time-Of-Flight detector [36].Furthermore, the data taken with each of the two possible directions of the magnetic field in the central barrel are both considered for the analysis.Finally, an additional event criterion is applied after the track selection to reject all remaining events with less than ten reconstructed tracks.This requirement is motivated by the calculation of the multiparticle azimuthal correlators themselves, and their corresponding event-by-event multiplicity weights.For a k-particle correlator, at least k tracks are needed to have a valid expression for the correlator and corresponding multiplicity weight.This criterion is set to ten tracks, as this is the minimum number of tracks needed to calculate the largest correlator in this analysis, AC 4,1 (m, n) (see Appendix C in Ref. [16] for the full estimator).After all event selections are applied, approximately 50 million events in the centrality range of 5% to 60% are used in the analysis.Furthermore, this range is divided into the following centrality classes: the two central classes 5-10% and 10-20%, and the four semicentral classes 20-30%, 30-40%, 40-50%, and 50-60%.
Tracks are reconstructed using the combined information given by the ITS and the TPC.These tracks are selected to be within |η| < 0.8 and 0.2 < p T < 5 GeV/c, where p T is the transverse momentum.The low p T cutoff is implemented for reducing biases due to a smaller reconstruction efficiency at very low transverse momentum, while the high p T cutoff is chosen for suppressing the contribution to the anisotropies in the measured azimuthal correlations from jets.In order to avoid contributions from secondary particles, each track is required to have a distance of closest approach (DCA) to the primary vertex of less than 2 cm in the longitudinal direction and less than 0.0105 + 0.0350/p 1.1 T cm in the transverse direction, meaning a value ranging from 0.016 cm at p T = 5.0 GeV/c to 0.22 cm at p T = 0.2 GeV/c.Additionally, each track is required to have a minimum of 70 TPC space points out of the maximum of 159, and a χ 2 value per space point from the track fit to be within the range of 0.1 to 2.5.The tracks are required to have a minimum of two hits in the ITS.Requirements were also added to reject tracks with a kink topology (i.e., tracks that appear to change direction due to multiple scatterings and/or weak decays).
Corrections for both non-uniform reconstruction efficiency (NUE) as a function of transverse momentum and non-uniform acceptance (NUA) as a function of azimuthal angle are applied in the same way as in Refs.[11,19,37].For the NUE corrections, a HIJING (Heavy-Ion Jet INteraction Generator) [38] simulation with a GEANT3 [39] detector model is used to correct for the effects of track reconstruction efficiency and contamination from secondary particles by constructing a p T -dependent track weighting correction.The NUA corrections are obtained in a data-driven manner, based on the method presented in Ref. [11].Finally, the nonflow contribution from the two-particle correlations appearing in the denominator of the NAC is suppressed with the use of a pseudorapidity gap |∆η| > 1.As the ACs themselves were shown to be robust against nonflow [16], no gap is used in their computation.
Obtaining the statistical uncertainty through standard error propagation is impractical due to both nontrivial functional dependence and many covariance terms.Therefore, in this analysis, the statistical uncertainties were determined using the standard subsampling procedure (see Appendix F in Ref. [16]).
The systematic uncertainties are estimated by independent variations of each selection criterion.The significance of each variation is determined according to the Barlow test [40].For each variation, the results for AC a,b (m, n) are compared with the ones obtained with the default selection, taking into account the correlations between their statistical uncertainties.If the deviation is larger than one σ , where σ is the uncertainty on the difference, the variation is deemed significant according to the Barlow test and is included into the quadratic sum to get the total systematic uncertainties.Each selection criterion variation is described below and the trend of their relative deviations with respect to the default values observed in semicentral collisions (collisions in the centrality range 20-60%) is reported in parentheses.Due to the small size of the signal, the variations measured in central events (centrality up to 20%) and for higher-order observables can fluctuate more significantly than the ones indicated here.
The effects of the centrality determination are estimated by changing the default V0 estimator to the SPD (20%).The quality of the event selection is further tested by varying the longitudinal position of the primary vertex from ± 8 cm to ± 6 cm (10%), and by applying a tighter pileup rejection factor (4%).Furthermore, the impact of the two configurations of the magnetic field polarity in the solenoid magnet of ALICE is investigated by performing the same analysis on the data sets taken for each orientation separately (5%).The minimum number of space points in the TPC required for the track reconstruction is increased from 70 to 80 (2%), and the quality of the reconstruction fit per TPC cluster is reduced from 0.1 < χ 2 < 2.5 to 0.1 < χ 2 < 2.3 (4%).The longitudinal DCA is also tightened from 2 cm to 1 cm.It has to be noticed that this systematic variation mainly affects the higher-order combinations (3-60%).Finally, a change of the default tracking requirements to the so-called hybrid tracks, which are TPC-ITS tracks with looser DCA requirements in the transverse plane and along the beam axis, leads to an effect on all studied observables (0-30%).

Results and discussion
The results obtained in this analysis are first presented together, in order to gain insights on the strength of the different multiharmonic correlations and their dependence on centrality and cumulant order.Then, specific comparisons between the data and the model predictions are presented for a few selected ACs and NACs.
The centrality dependence of the thirty-four different AC a,b (m, n) and NAC a,b (m, n) is shown in Fig. 1.However, due to their large uncertainties, NAC 4,1 (m, n) and NAC 1,4 (m, n) are not presented here.As the results for AC a,b (2, 3) and AC a,b (2, 4) (Fig. 1 (a) and 1 (c), respectively) present similar behaviors, albeit with different signs, they are here discussed together.In the following lines, the notation (m, n) thus stands for both pairs of harmonics (2,3) and (2,4).A dependence of the magnitude of the genuine correlations on the cumulant order is observed.The results for SC(m, n) ≡ AC 1,1 (m, n) present the largest signals, followed by AC 2,1 (m, n), while AC 3,1 (m, n) and AC 4,1 (m, n) are several orders of magnitude smaller.This decreasing effect is more pronounced for the mirror combinations AC b,a (m, n) when compared to their respective counterparts.The splitting between the results for AC a,b (m, n) and its mirror AC b,a (m, n) is also present for NAC a,b (2, 3) and NAC a,b (2, 4) (Fig. 1 (b) and 1 (d), respectively).While the results for NAC 2,1 (2, 3) and NAC 3,1 (2, 3) demonstrate a trend similar as NSC(2, 3), albeit with smaller magnitudes, the ones for NAC 2,1 (2, 4) and NAC 3,1 (2,4) are in good agreement with NSC(2, 4) and deviate only in more peripheral collisions.In both panels, the mirror combinations (open markers) show an increase of the magnitude with the cumulant order.Finally, Figs. 1 (e) and 1 (f) present the centrality dependence of the different orders of AC a,b (3,4) and NAC a,b (3,4), respectively.The results for AC a,b (3,4) present the same ordering of the magnitudes with the order of the AC itself, as observed in the two other pairs of harmonics.However, no splitting between AC a,1 (3,4) and AC 1,a (3,4) can be seen at any order a = 1, . . ., 4 within uncertainties.
The measurements are compared with results from simulated events generated with the state-of-the-art chain of models T R ENTo+VISH(2+1)+UrQMD [41][42][43][44], which will be referred to here as T R ENTo+iEBE-VISHNU.This model chain consists of the T R ENTo model [45] for the initial conditions, which are connected by free streaming to a 2+1-dimensional hydrodynamical model VISH(2+1) [41].The description of the heavy-ion evolution after the collision is continued after particlization with a hadronic cascade (UrQMD) model [43,44].The initial conditions, the specific shear viscosity η/s(T ), the specific bulk viscosity ζ /s(T ), and other free parameters of this hybrid model, for instance, the Gaussian nucleon width and the minimum inter-nucleon distance, are obtained from a global Bayesian analysis.The datato-model comparisons are performed with two sets of best-fit parametrizations given by maximum a posteriori (MAP) from two independent Bayesian analysis groups.As the values for these two sets of best-fit parameters agree within their respective uncertainties from the Bayesian analysis, this work allows the observation of the deviations between the two parametrizations in the predictions of observables that were not used in the Bayesian optimization.The extraction of the Duke 2019 parametrization [22] used a series of ALICE measurements (two-particle correlations, charged-particle multiplicities, etc.) as inputs into the Bayesian analysis.However, most of the used data are from Pb-Pb collisions at √ s NN = 2.76 TeV and only a subset of available flow observables were used.The Jyväskylä 2022 parametrization [12] utilized for the first time higher-order flow observables and data from Pb-Pb collisions at both √ s NN = 2.76 TeV and √ s NN = 5.02 TeV.New observables, such as SC [13,14,17,19] and flow harmonic mode couplings [37], recently measured by the ALICE Collaboration, are more sensitive to the hydrodynamic transport coefficients η/s and ζ /s.This resulted in a reduction of the uncertainties of the extracted transport coefficients [12].In these hydrodynamical model calculations, all the kinematic selections, such as p T and η ranges, are matched with the ones reported in this article.In addition, the observables are extracted from the simulations using the same experimental analysis techniques.
Figure 2 shows the comparisons of the experimental data for NAC a,b (2, 3) with the predictions from hydrodynamic simulation using the Duke 2019 and Jyväskylä 2022 parametrizations.As discussed above for the data, the same splitting between NAC a,1 (2, 3) and NAC 1,a (2, 3) (a = 2, 3) can be observed in both model calculations.The negative sign is visible both for NSC(2, 3) (Fig. 2 (a)) and for some AC observables in the same pair of flow harmonics.The notable exception here is NAC 3,1 (2, 3), which is close to zero for the data and predicted to be positive in semicentral collisions by both model parametrizations  The data points for NSC(2, 3) are taken from Ref. [17].
(see Fig. 2 (d) for a close-up of the results).Furthermore, no major difference can be seen between the predictions from Duke 2019 and Jyväskylä 2022 for any of the observables in Fig. 2. As the linear response v n = k n ε n (n = 2, 3), where k n are the linear response coefficients and ε n the eccentricities in the initial state, dominates in central collisions, the genuine correlations between v 2 and v 3 are expected to be sensitive only to the initial conditions in this regime.Therefore, a possible explanation of the similar predictions is the use of T R ENTo for both parametrizations.Future studies would thus benefit from the comparison between different initial state models.It also has to be noted that both calculations reproduce the experimental data only qualitatively, and thus the use of such new measurements as input in future Bayesian studies would allow further constraints on the descriptions of the early stages of heavy-ion collisions.
The comparison of the results for AC a,b (2, 4) with the calculations from both model parametrizations is shown in Fig. 3. Overall, the behavior of the experimental data is better reproduced by the Jyväskylä 2022 parametrization, quantitatively for SC and qualitatively for the higher orders.The model predictions are also closer to the experimental values for AC 1,a (2, 4) than for AC a,1 (2,4), for a = 2, 3, 4. The observable AC a,1 (2,4) presents the same positive sign as SC (2,4).This is the case as well for AC 1,2 (2, 4) and AC 1,3 (2,4) in semicentral collisions (see Fig. 3 (e) for a close-up of the latter observable).However, both experimental and theoretical results for AC 1,4 (2, 4) are in agreement with zero within uncertainties up to centralities of 30% for the data and 45% for the model (Fig. 3 (f)).Because of previous nonvanishing results for SC (2,4), this indicates that the used sample of Run 2 data is not sufficient to extract nonvanishing values for these particular AC observables.Past studies have shown the existence of a non-linear response between the fourth-order flow amplitude, v 4 , and the initial-state ellipticity term ε 2 2 , which is dominant in non-central collisions [46,47].Furthermore, higher-order harmonics are more sensitive to small variations of η/s and ζ /s than the lower-order coefficients [47][48][49].This latter point was further discussed in recent Bayesian estimations [12], where higher-order cumulants of flow amplitudes exhibit better sensitivity to QGP properties when compared to the traditional lower-order flow observables, generalizing the findings from Refs.[13,14,49].The remaining tensions between data and predictions observed in the current analysis for AC 1,a (2,4) are thus a sign that the measurements of the genuine correlations between v 2 2 and v 2a 4 can bring additional input on the higher-order terms of the nonlinear response of v 4 .Comparisons between the initial and final state predictions for these observables would be one of the crucial steps in that direction, as any discrepancies between those would shed more light on the hydrodynamic description of v 4 in the models.
The experimental results for NAC a,b (3,4) are compared with the theoretical calculations from the two model parametrizations, as shown in Fig. 4. Similarly to the results for the data, the predictions for NAC 2,1 (3,4) and NAC 1,2 (3,4) are in good agreement with each other within uncertainties (Fig. 4 (b)), for both the Duke 2019 and Jyväskylä 2022 parametrizations.Additionally, previous studies [46] have shown that the measurements of v 4 can be explained by a contribution from ε 2 2 and none from ε 3 .The observed agreement may then originate from the interplay between the different stages of the evolution of the heavy-ion collisions.For instance, an agreement between the initial-state predictions for NAC 2,1 (3,4) and NAC 1,2 (3,4) could indicate an absence of non-linear response between the higher orders of ε 3 and v 4 .Similarly, different initial-state predictions could mean that the hydrodynamic evolution differently impacts the correlations between the higher moments of ε 3 with v 4 and the ones between ε 3 and the higher moments of v 4 .However, these non-linear hydrodynamic effects would be such that the two observables converge onto each other in the final state.Future studies of both the initial and final state, for instance, comparisons of the AC results from this study with initial-state predictions, would help in understanding the origin of the different behaviors -splitting or no splitting -observed in the various combinations of NAC.Furthermore, it is noted that there is good agreement between the predictions obtained from the Duke 2019 and Jyväskylä 2022 parametrizations, and that they both capture quantitatively the experimental data for centralities higher than 15%.However, none of the parametrizations capture the sign change between central and semicentral events observed in the experimental data.In Ref. [17], this effect was already observed for NSC (3,4), where the sign change could be reproduced by models using AMPT initial conditions but not by the ones using T R ENTo.It is therefore expected that future precision measurements of NAC a,b (3,4) and their use as input data in Bayesian studies would help in tuning the different details of system evolution as a function of the collision centrality in the models.The main findings of this study on ACs and NACs are summarized here as follows: (i) The magnitude of the NAC observables varies with different moments as well as with the collision centrality, indicating that non-linear responses are also reflected in higher moments.
(ii) While the NAC 2,1 (2, 3) observable shows a negative decreasing magnitude toward peripheral col-  lisions, the higher-moment correlations, NAC 3,1 (2,3), are close to zero except in the most central collisions.The model calculations show a similar trend, but underestimate the data.
The AC 1,3 (2,4) observable is close to zero, while the AC 1,4 (2, 4) is compatible with zero within uncertainties.As non-zero measurements for SC(2, 4) [13] indicate that v 2 2 and v 2 4 are positively correlated, this does not mean that the higher-order fluctuations of v 2 2 and v 2 4 are uncorrelated, but that the used data sample is not enough to extract their nonvanishing values.
The agreement between model and data for each cumulant order is now quantified with the χ 2 test (χ 2 /N dof ), performed as in Ref. [14] where y i ( f i ) is a measurement (model) value in a centrality bin i.The systematic and statistical uncertainties from the data are combined in quadrature σ 2 i = σ 2 i,stat +σ 2 i,syst +σ 2 f i ,stat , together with the statistical errors of the model calculations.The quantity N dof corresponds to the number of bins in the centrality range used in χ 2 calculations.The results are shown in Fig. 6 for both model parametrizations indicated on the upper x-axis.A grey box indicates that the experimental data have large uncertainties, which may influence the interpretation of the obtained χ 2 value.For instance, an agreement between these large uncertainties of the data and the predictions could bias the χ 2 test, leading to an artificially small value.[22] and Jyväskylä 2022 [12] parametrizations.The statistical uncertainties in the calculations are indicated by the thicknesses of the colored bands.The data points for NSC (3,4) are taken from Ref. [17].
Such bias would then be resolved with a larger data sample.As they were not measured in the current analysis, NAC 1,3 (2,4) and NAC 1,3 (3,4) are shown with a white filling.A better agreement between the data and a model parametrization is represented by lower χ 2 values.In Fig. 6, it is observed that the model description is better for higher-order harmonics such as all AC a,b (3, 4), and higher-order moments such as AC 4,1 (m, n), as well as most of NAC 2,1 (m, n) and NAC 3,1 (m, n).Generally, the Jyväskylä 2022 parametrization describes the data better than the Duke 2019 parametrization.This result is expected as a larger set of data and higher-order observables are included in Jyväskylä 2022.However, the Duke 2019 still provides better predictions than the newer parametrization in the case of SC(2, 3), AC 3,1 (2,3) and NSC(m, n) for all three pairs of harmonics.It is also known that Jyväskylä 2022 has a deviation in its prediction of v 2 at √ s NN = 5.02 TeV [12], which may explain the differences observed in the descriptions of ACs involving v 2 but not of the NACs.Ultimately, the χ 2 test indicates that the independent AC observables are needed as input for future Bayesian analyses in order to reduce the uncertainty on the initial conditions and/or on the hydrodynamic transport coefficients such as η/s(T ) and ζ /s(T ).

Summary
In conclusion, the first measurements of the correlations between different moments of two flow amplitudes in Pb-Pb collisions at √ s NN = 5.02 TeV, using asymmetric cumulants, are shown.This work is a generalization of the previous measurements obtained with symmetric cumulants, proven to have an increased sensitivity to the initial conditions and medium properties.Due to the cumulant uniqueness, each AC has access to new and independent information unreachable by measuring the SC only.[22] and Jyväskylä 2022 [12] parametrizations.The last point of NAC 1,2 (2,4) and the results for NAC 1,3 (2,3) are not shown for visibility purposes.The statistical uncertainties in the calculations are indicated by the thicknesses of the colored bands.
The data points for a = 1 are taken from Ref. [17].
Therefore, this article provides new measurements and resulting additional constraints for 34 new observables.A dependence of the magnitude of the correlations on the moments and the centrality of the collision is observed, showing an influence of the non-linear response between the higher-moments as well.Furthermore, it has been seen that an increasing (decreasing) dependence with the moments for a given pair (m, n) in NAC a,1 (m, n) generally leads to the decreasing (increasing) dependence for the mirror combination NAC 1,a (m, n).The comparison of the experimental results with the predictions from two different parametrizations of the T R ENTo+iEBE-VISHNU model chain obtained with Bayesian analyses shows a need for further improvements of the initial conditions and non-linear response.Due to the high sensitivity of these measurements to the model parameters [12,27], the results can improve the current understanding of the detailed nucleon structure in nuclei [50,51], the nuclear skin effect [52], and the pre-equilibrium phase [53][54][55][56][57][58].Incorporating these details into the models used by the Duke and Jyväskylä groups will lead to a better quantitative understanding of the properties of QCD matter.NACs and all pairs of harmonics presented in this paper.The observables with large uncertainties are indicated with a grey filling, while an absence of observable (see Fig. 1) is shown with a white filling.

A Additional observables
The comparisons of the model predictions to the experimental data for NAC  [12], could be at play in SC(2, 3) and cancel out in NSC(2, 3).
In      [22] and Jyväskylä 2022 [12] parametrizations.The statistical uncertainties in the calculations are indicated by the thicknesses of the colored bands.

Figure 1 :
Figure 1: Centrality dependence of the different orders of AC a,b (m, n) (left) and NAC a,b (m, n) (right) for four pairs of indices (a, b) and three pairs of harmonics (m, n) in Pb-Pb collisions at √ s NN = 5.02 TeV.The mirror combinations of AC a,b (m, n) (closed markers), i.e., AC b,a (m, n), are indicated with open markers and similar color scheme.The statistical (systematic) uncertainties are shown with lines (boxes).The data points are shifted horizontally for visibility.

Figure 5 :
Figure5: Centrality dependence of the different orders a = 1, . . ., 3 of NAC a,1 (m, n) and NAC 1,a (m, n) compared with the theoretical predictions from T R ENTo+iEBE-VISHNU for the Duke 2019[22] and Jyväskylä 2022[12] parametrizations.The last point of NAC 1,2 (2, 4) and the results for NAC 1,3(2,3) are not shown for visibility purposes.The statistical uncertainties in the calculations are indicated by the thicknesses of the colored bands.The data points for a = 1 are taken from Ref.[17].

Figure 6 :
Figure 6: Values for χ 2 test calculated between the data and two different model calculations for all ACs andNACs and all pairs of harmonics presented in this paper.The observables with large uncertainties are indicated with a grey filling, while an absence of observable (see Fig.1) is shown with a white filling.
the case of NAC a,b (2, 4) (Fig. A.2), a slightly better agreement with the data from Jyväskylä 2022 over Duke 2019 parametrizations can be observed.As stated above, the non-linear response between v 4 and ε 2 2 makes AC a,b (2, 4), and thus NAC a,b (2, 4), especially sensitive to the hydrodynamical description of the collision.
(2,3),3), AC a,b(2,4), and NAC a,b(3,4)are presented in the main part of this article.In this Appendix, the comparisons for the other observables of this analysis are discussed.The complete list of combinations can be found in Figs.1 and 6.The predictions from Duke 2019 and Jyväskylä 2022 parametrizations and the experimental data for AC a,b (2, 3) and NAC a,b(2,4)are presented in Figs.A.1 and A.2, respectively.The same splitting with the cumulant order as observed in Fig.3for NAC a,b (2, 3) and AC a,b (2, 4) can be seen here.In both cases, NAC a,1(2,4)show increasingly larger magnitudes than NAC a,1 (4, 2), while the opposite stands true for the unnormalized AC.The observables AC a,b (2, 3) and AC a,b (2, 4) (Figs.A.1 and 3) show similar agreements between predictions and data.The Jyväskylä 2022 parametrization has overall a better agreement with the data than the Duke 2019 one, and the higher-order ACs are overall better reproduced than SC(m, n) and AC 2,1 (m, n).It can be noted, however, that while the Jyväskylä 2022 parametrization reproduces quantitatively SC(2,4), whereas the Duke 2019 one does not; the opposite behavior is observed for SC(2,3).On the other hand, no big difference in the predictions of NAC a,1 (2, 4) and NAC a,1 (4, 2) between the Duke 2019 and Jyväskylä 2022 parametrizations can be observed.Possible effects, such as the ones at the origin of the 10% discrepancy in the description of v 2 at 5.02 TeV by the Jyväskylä 2022 parametrization