1 Introduction

One of the most important discoveries in the physics of heavy-ion collisions at ultrarelativistic energies is the observation of a deconfined state of nuclear matter dubbed quark–gluon plasma (QGP). This extreme state is produced during the heavy-ion collision evolution, and its properties resemble the properties of a perfect liquid. Unprecedentedly large data sets collected at the LHC enable the most quantitative description of the QGP to date. Given the complexity of the system produced in heavy-ion collisions, an important program in the field is the development of observables that provide new and independent information inaccessible with previous measurements [1,2,3,4,5,6,7,8].

The intersecting volume of two heavy ions is anisotropic in coordinate space, either due to collision geometry (particularly in non-central collisions with large values of impact parameter) or due to fluctuations of positions of participating nucleons (most significant in central head-on collisions). Anisotropic pressure gradients, which develop in this volume containing the strongly interacting nuclear matter, transfer the initial-state spatial anisotropies into final-state anisotropies in momentum space. This phenomenon is known as anisotropic flow and it is a sensitive probe of all stages in the heavy-ion collision evolution [9]. Anisotropic flow measurements are used to constrain the transport properties of the QGP, for instance ratios of shear and bulk viscosities to entropy density  [4, 7, 10,11,12,13]. The anisotropic emission of particles in the plane transverse to the beam direction is quantified with amplitudes \(v_n\) and symmetry planes \(\Psi _n\) by using the Fourier series decomposition of the azimuthal angle (\(\varphi \)) distribution of produced particles [14]

$$\begin{aligned} f(\varphi )=\frac{1}{2\pi }\left[ 1+2\sum _{n=1}^{\infty }v_n\cos [n(\varphi -\Psi _n)]\right] . \end{aligned}$$
(1)

A detailed discussion of properties of \(v_n\) and \(\Psi _n\) can be found in Ref. [15]. The symmetry plane \(\Psi _n\) has a simple geometrical interpretation when the anisotropic distribution can be parameterized only with one harmonic n, since then it can be shown that \(f(\Psi _n+\varphi ) = f(\Psi _n-\varphi )\), i.e. symmetry plane \(\Psi _n\) is the plane for which it is equally probable for a particle to be emitted above or below it.

Historically, the emphasis was on studying the amplitudes \(v_n\), but the symmetry planes also carry a very important information about different stages in heavy-ion collision evolution. Unlike the flow amplitudes \(v_n\), a single symmetry plane \(\Psi _n\) cannot be estimated directly in an experiment using correlation techniques – the simplest available observables are symmetry plane correlations (SPC), for instance \(\langle \cos 4(\Psi _4-\Psi _2)\rangle \) [15, 16]. Such correlations are the subject of this study.

In the early anisotropic flow analyses, the goal was to measure \(v_n\) with respect to the reaction plane (a plane spanned by the beam axis and impact parameter vector), and it was assumed that all symmetry planes are approximately the same and equal to the orientation of the reaction plane. Therefore, the first flow measurements were exclusively of flow amplitudes \(v_n\). The first experimental results for SPC can be traced back to the E877 experiment [17]. These initial measurements were performed by the standard event plane method with the subevent technique [18]. The first measurements of SPC involving two symmetry planes in the RHIC era were obtained by PHENIX in Refs. [19, 20]. An alternative approach was pursued by NA49 and STAR using 3-particle mixed-harmonic correlations, which by definition have contributions from SPC [21, 22]. In the first flow studies at LHC energies, the ALICE Collaboration demonstrated in Ref. [23] that the symmetry planes \(\Psi _2\) and \(\Psi _3\) fluctuate independently in all considered centralities. Finally, the most detailed experimental results to date were published by the ATLAS Collaboration in Ref. [24], where also for the first time the strength of correlations among three symmetry planes was presented. ATLAS systematically studied the centrality dependence of SPC both in the initial and final state using the analysis technique from Refs. [16, 25]. It was concluded that SPC originate both from correlated fluctuations in the initial geometry and from the non-linear mixing between different flow harmonics in the final state. Subsequent experimental publications which used SPC to constrain the details of the non-linear hydrodynamic response can be found in Refs. [26,27,28,29,30].

In theoretical studies, SPC can be obtained directly both in coordinate and in momentum space [16, 24, 25, 31,32,33,34,35,36,37,38]. State-of-the-art modeling of heavy-ion collisions covers all stages of its evolution starting from the initial conditions to the final free streaming of produced particles. The SPC in the initial state can be obtained event-by-event directly from the underlying model of the collision geometry using for instance energy density distribution or nucleon positions, while in the final state SPC are the event-by-event output of the model used to describe all subsequent stages in the evolution. Therefore, in theoretical studies it is not, in general, necessary to build an estimator for SPC from the azimuthal angles of final-state particles, like it is done in an experiment. In order to ease the comparison between theoretical and experimental results, azimuthal correlators were used to indirectly estimate SPC also in Refs. [10, 39,40,41,42,43,44]. Other types of theoretical studies involving symmetry planes can be found in Refs. [45,46,47,48,49,50].

Several experimental difficulties associated with the SPC render their measurements particularly challenging. Even in the simplest realisation, it is necessary to construct non-trivial estimators for SPC to resolve these issues. Unlike the flow amplitudes \(v_n\), each symmetry plane \(\Psi _n\) taken individually is not invariant under rotations of the coordinate system in the laboratory frame in which azimuthal angles are measured (see Eq. 1). Therefore, the simplest rotationally-invariant physical observable involving symmetry planes is the difference of two symmetry planes. In an actual experiment such rotations are unavoidable as a direct consequence of random event-by-event fluctuations of the direction of the impact parameter vector. Only symmetry planes that are different, apart for trivial periodicity, carry independent information, and therefore any dependence on periodicity must be removed from all SPC observables by definition. The widely used technique to suppress systematic biases from short-range nonflow correlations by introducing pseudorapidity gaps in the measured azimuthal correlators which are used to estimate SPC is not applicable due to decorrelations of symmetry planes as a function of pseudorapidity [50,51,52,53,54,55,56]. Moreover, it has been shown recently that the effect of flow magnitude correlations, which have been either completely [24] or partially [27] neglected in the existing measurements, may overshadow the correlations of symmetry planes in the analysis with the Scalar Product (SP) method [15]. The new and improved estimator for SPC from Ref. [15], which overcomes these limitations, is introduced next.

The starting point is the following relation between \(v_n\) and \(\Psi _n\), and multiparticle azimuthal correlations [15, 39, 57]:

$$\begin{aligned} v_{n_1}^{a_1}\cdots v_{n_k}^{a_k}\,e^{i(a_1 n_1\Psi _{n_1}+\cdots +a_k n_k\Psi _{n_k})} = \left<e^{i(n_1\varphi _1+\cdots +n_l\varphi _l)}\right> . \end{aligned}$$
(2)

In this equation, angular brackets indicate an average over the azimuthal angles of all distinct sets of l particles measured in the same event.

The coefficients \(a_i\) are positive integers which ensure that all harmonics \(n_i\) and symmetry planes \(\Psi _{n_i}\) are unique on the left-hand side in the above expression. These coefficients can be understood in the following way: \(a_i\) counts how many times a harmonic \(n_i\) appears in the azimuthal correlator on the right-hand side of Eq. (2) (harmonics with positive and negative signs are counted separately). The total number of particles, i.e. the order of the multiparticle azimuthal correlator, is given by \(\sum _i a_i\). The index k on the left-hand side labels only unique harmonics in the original set \(n_1, n_2, \ldots , n_l\), therefore \(k\le l\). As an example, for the correlator \(\left<e^{i(2\varphi _1-\varphi _2-\varphi _3)}\right>\) it follows that \(n_1=2, a_1=1, n_2=n_3=-1, a_2=2\). The advantage of this generalized notation is that now \(n_i\) and \(a_i\) decouple naturally either into a subscript or into an exponent when associated with flow amplitudes \(v_{n_i}\) in Eq. (2), which enables their distinct physical interpretation. Finally, solely from the definition of the Fourier series in Eq. (1) one can prove that \(v_{-n} = v_n\) and \(\Psi _{-n} = \Psi _n\), which is used in the rest of the paper. Due to this property, the final a coefficient for harmonic n in Eq. (2) is a sum \(a_n + a_{-n}\).

Taking into account all these technical considerations, the simplest definition of SPC observables is provided by the following expression [24, 25, 39]:

$$\begin{aligned} \left<e^{i(a_1n_1\Psi _{n_1} + \cdots + a_kn_k\Psi _{n_k})}\right>,\quad \sum _i^k a_i n_i = 0, \end{aligned}$$
(3)

where all \(a_i\) are positive and all \(n_i\) are unique integers. Angular brackets \(\langle \ \rangle \) indicate here an average over all events. Defined this way, SPC observables are rotationally invariant and therefore invariant with respect to random event-by-event fluctuations of the impact parameter vector, while the periodicity of each individual symmetry plane is accounted for by definition. Experimentally, Eq. (2) is used as a starting point for an estimator for SPC. However, to isolate the true SPC part, the prefactor \(v_{n_1}^{a_1}\cdots v_{n_k}^{a_k}\) has to be divided out. The importance of this technical detail was neglected in all previously used SPC estimators.

The new and improved SPC estimator, named the Gaussian Estimator (GE), was developed recently in Ref. [15]. Its key improvement amounts to using the following expression to estimate SPC:

$$\begin{aligned} \begin{aligned}&\langle \cos \left( a_{1} n_1 \Psi _{n_1}+\cdots + a_{k} n_k \Psi _{n_k}\right) \rangle _{\text {GE}} \\&\quad = \sqrt{\frac{\pi }{4}}\;\;\;\frac{\langle v_{n_1}^{a_{1}} \;\cdots \; v_{n_k}^{a_{k}}\cos \left( a_{1} n_1 \Psi _{n_1}+\cdots +a_{k} n_k \Psi _{n_k}\right) \rangle }{\sqrt{\langle v_{n_1}^{2a_{1}} \;\cdots \;v_{n_k}^{2a_{k}} \rangle }}, \end{aligned} \end{aligned}$$
(4)

which was derived by approximating multi-harmonic flow fluctuations with a two-dimensional Gaussian distribution. Both the numerator and denominator on the right-hand side in the above expression can be estimated by using Eq. (2) with suitably chosen harmonics \(n_i\). Further explanations of the technical details of the GE based on the example \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle \) are provided in Appendix A. The main improvement of this new estimator can be found in the denominator where the GE has the joined multivariate moment of different flow amplitudes, \(\langle v_{n_1}^{2a_{1}} \;\cdots \;v_{n_k}^{2a_{k}}\rangle \). This is in contrast to the previously used SP estimator, defined as [43]

$$\begin{aligned}{} & {} \langle \cos \left( a_{1} n_1 \Psi _{n_1}+\cdots +a_{k} n_k \Psi _{n_k}\right) \rangle _{\text {SP}}\nonumber \\{} & {} \quad =\frac{\langle v_{n_1}^{a_{1}} \;\cdots \; v_{n_k}^{a_{k}} \cos \left( a_{1} n_1 \Psi _{n_1}+\cdots +a_{k} n_k \Psi _{n_k}\right) \rangle }{\sqrt{\langle v_{n_1}^{2a_{1}} \rangle \cdots \langle v_{n_k}^{2a_{k}} \rangle }}, \end{aligned}$$
(5)

which uses instead \(\langle v_{n_1}^{2a_{1}}\rangle \;\cdots \;\langle v_{n_k}^{2a_{k}} \rangle \) in the denominator and therefore assumes that event-by-event fluctuations of flow amplitudes are mutually independent. This assumption is in contradiction with recent experimental results which reported strong and non-trivial correlated fluctuations of different flow amplitudes, both at RHIC and LHC energies, and across different collisions systems [26, 58,59,60,61]. These shortcomings of the previous SPC results are the main motivation for the current work. As it was pointed out in Ref. [15], the GE does not account for cross-correlations between flow amplitudes and symmetry planes. However, the study in Ref. [15] showed that the contribution by these cross-correlations is minor when compared to the correlations between flow amplitudes.

The rest of the article is organized as follows. In Sect. 2 the ALICE detector is introduced, together with the analyzed data set and analysis details, such as the event and track selection criteria. In Sect. 3 the SPC results using the GE are presented, comparisons with previous experimental results are discussed, and confrontation with state-of-the-art theoretical models is displayed. The article concludes in Sect. 4 with the summary. A more detailed discussion about the technical details of the GE can be found Appendix A.

2 Data analysis

The data set consists of Pb–Pb collisions at a center-of-mass energy per nucleon pair \(\sqrt{s_{\text {NN}}}=2.76~\text {TeV}\) recorded by ALICE in 2010. A detailed description of the apparatus and its performance is given in Refs. [62, 63]. The silicon pixel detector (SPD), which comprises the two innermost layers of the inner tracking system (ITS) [64, 65], and both V0 detectors  [66] were used for triggering. The latter consists of two arrays of scintillator counters, the V0A and V0C, covering a pseudorapidity range of \(2.8<\eta <5.1\) and \(-3.7<\eta <-1.7\), respectively. The SPD covers pseudorapidities of \(|\eta | < 2.0\) for its inner and \(|\eta | < 1.4\) for its outer layer. Minimum bias collisions were selected by requiring a signal in at least two out of the three following: two chips in the outer layer of the SPD, the V0A, and the V0C. For this analysis, only events with a primary vertex within \(\pm 10~\text {cm}\) of the nominal interaction point along the beam axis were used. The centrality of the collisions [67] was estimated with the SPD. Backgrounds events due to beam–gas interactions and parasitic beam–beam interactions were removed by using V0 and Zero Degree Calorimeter [68] timing information. Overall, after the event selection the used data set consists of \(7.36\times 10^6\) reconstructed collisions for the centrality range 0–50%.

The reconstruction of charged particle trajectories was performed using only information from the time projection chamber (TPC) [69] due to its uniform acceptance in azimuth. This analysis used tracks with transverse momenta \(0.2<p_{\text {T}}<5.0~\text {GeV}/c\) and in a pseudorapidity range of \(|\eta |<0.8\), while covering the full azimuth. The lower boundary of the transverse momentum selection ensured a large and stable tracking efficiency in the TPC, while the upper cutoff decreases the contribution from jets which in general have larger momenta. The charged tracks were accepted for the analysis if they had a minimum of 70 out of a maximum of 159 space points in the TPC. The \(\chi ^2\) per space point from the track fit was set to be within \(0.1<\chi ^2/ \textrm{NDF}<4.0\). The distance of closest approach (DCA) of the extrapolated tracks to the primary vertex was required to be at maximum 2.4 cm in the transverse plane and 3.2 cm in the beam direction. Daughter tracks with a reconstructed secondary weak-decay kink topology (i.e. tracks with an abrupt change of direction) were discarded. The contamination from secondaries as well as the reconstruction efficiency with this track selection can be found in Ref. [70].

The \(p_{\text {T}}\)-dependent reconstruction efficiency was corrected using particle weights according to Ref. [57]. These weights were obtained with the HIJING (Heavy-Ion Jet INteraction Generator) Monte Carlo generator [71] by comparison of generated and reconstructed tracks. For the latter, a GEANT3 [72] detector simulation and event reconstruction was used in addition to HIJING. At the same time, weights to correct for non-uniform acceptance in azimuthal angle did not have to be applied due to the uniform acceptance of the TPC over the whole azimuth in the analyzed data set. Nonflow contributions, i.e. correlations between a few particles unrelated to collective anisotropic flow, were investigated with HIJING for the numerator and denominator of the GE in Eq. (4) separately. For all SPC combinations, both the numerator and denominator were found to be consistent with zero in all considered centrality ranges, demonstrating that the analyzed SPC observables are not influenced by most important sources of nonflow correlations such as jets or resonance decays.

The statistical uncertainties of the measured SPC were obtained via propagation of uncertainties of the numerator and denominator in Eq. (4). Systematic uncertainties were evaluated by varying the default event and track selections. All variations were performed one at a time and only those with a difference larger than \(2\sigma \), where \(\sigma \) is the uncertainty of the difference, with respect to the default selection were taken into account for the final systematic uncertainty. All individual systematic variations were considered independent and combined in quadrature to obtain the total systematic uncertainty. Regarding the event selection criteria, the position of the primary vertex along the beam line was varied to \(\pm 6~\text {cm}\) and \(\pm 8~\text {cm}\), where a relative effect on the measured observables of up to 5% was found. A systematic uncertainty of up to 6% from the centrality estimation was determined by using the V0 instead of the SPD. To evaluate the uncertainty due to the track selection, the number of TPC clusters used in the track reconstruction was varied to a required minimum of 80, 90 and 100 compared to the default 70. This resulted in a systematic uncertainty of up to 6%. The sensitivity of the results to the track quality was checked by varying the \(\chi ^2/ \textrm{NDF}\) to \(0.3<\chi ^2/ \textrm{NDF}<4.0\) and \(0.1<\chi ^2/ \textrm{NDF}<3.5\), which led to an additional uncertainty of up to 6%. Two variations were performed regarding the DCA by changing the upper limit in the transverse direction to 1 cm and in the longitudinal direction to 2 cm. The variation of the DCA changes the contribution from secondaries in the analysis as these particles usually have a larger DCA than primary particles. The DCA variation in the transverse plane led to a systematic uncertainty of about 3–10%, while the check along the beam axis had a relative variation of about 4%. Additionally, an independent analysis was performed by using a different track reconstruction procedure, which employs combined information from both the TPC and the ITS. This led to an uncertainty in the range of 5–10%.

3 Results

Fig. 1
figure 1

Comparison of the extracted correlations between different combinations of two symmetry planes (a) and between three and four planes (b) using the GE in Eq. (4). Statistical (systematic) uncertainties are shown as lines (boxes)

The centrality dependence of the correlations between different combinations of two and three symmetry planes, as well as the first measurement of a correlation between four planes, are presented in Fig. 1. In the case of two symmetry planes shown in Fig. 1a, the strongest correlation is observed for \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle _\text {GE}\), while the correlation strength gets weaker for \(\langle \cos \left[ 6\left( \Psi _6-\Psi _2\right) \right] \rangle _\text {GE}\) and \(\langle \cos \left[ 6\left( \Psi _6-\Psi _3\right) \right] \rangle _\text {GE}\). The results for \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\) are compatible with zero within uncertainties. A hierarchy, \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle _\text {GE}\) > \(\langle \cos \left[ 6\left( \Psi _6-\Psi _3\right) \right] \rangle _\text {GE}\) > \(\langle \cos \left[ 6\left( \Psi _6-\Psi _2\right) \right] \rangle _\text {GE}\), holds for the centrality range 5–50%, with an exception of \(\langle \cos \left[ 6\left( \Psi _6-\Psi _3\right) \right] \rangle _\text {GE}\) and \(\langle \cos \left[ 6\left( \Psi _6-\Psi _2\right) \right] \rangle _\text {GE}\) being comparable at centralities above 20%. The details of the centrality dependence vary for the different combinations of symmetry planes. While \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle _\text {GE}\) and \(\langle \cos \left[ 6\left( \Psi _6-\Psi _2\right) \right] \rangle _\text {GE}\) are increasing non-linearly from central to semicentral collisions, \(\langle \cos \left[ 6\left( \Psi _6-\Psi _3\right) \right] \rangle _\text {GE}\) shows a weak centrality dependence. The observed zero signal for \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\) indicates that no correlation is present within the current uncertainties for the final-state planes \(\Psi _2\) and \(\Psi _3\), while \(v_2\) and \(v_3\) are anti-correlated [26, 44, 58, 59, 73,74,75]. This result justifies the necessity of measuring separately correlations of symmetry planes and flow magnitudes, because these measurements can be used to independently constrain properties of the matter produced in heavy-ion collisions.

The different magnitudes of correlations are also observed for three symmetry planes as shown in Fig. 1b. The magnitude and details of the centrality dependence vary for different combinations of flow harmonics. The \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\)Five observable exhibits the strongest correlations and \(\langle \cos \left[ 2\Psi _2-6\Psi _3+4\Psi _4\right] \rangle _\text {GE}\) shows the weaker signal. The SPC \(\langle \cos \left[ 2\Psi _2-6\Psi _3+4\Psi _4\right] \rangle _\text {GE}\) is the only correlator with a negative sign, which will be discussed later on in more detail. The \(\langle \cos \left[ 8\Psi _2-3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\) observable is consistent with zero within uncertainties, similar to \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\). The correlation between four planes, \(\langle \cos \left[ 2\Psi _2-3\Psi _3-4\Psi _4+5\Psi _5\right] \rangle _\text {GE}\) shows the strongest centrality dependence among all harmonic combinations and increases towards peripheral collisions.

The magnitudes of SPC are ordered approximately based on the corresponding order of the particle correlations. The two largest SPCs, \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle _\text {GE}\) and \(\langle \cos \left[ 2\Psi _2+3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\), are both measured with three-particle correlators. In contrast, the smallest ones are \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\) and \(\langle \cos \left[ 8\Psi _2-3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\), which are five- and six-particles correlations, respectively. One possible explanation is the following: the flow vector fluctuations encoded in the observed correlations are mainly attributed to the fluctuation of the initial state. Also, the initial state fluctuation is attributed to the fluctuation of a finite amount of “sources” produced at the degrees of freedom collision points, namely protons and neutrons, in the collision region. The central limit theorem (CLT) states that for independent random variables (here, the position of sources), the sample average tends toward a Gaussian distribution when the number of sampling increases. A clear example of such behavior was studied for initial ellipticity in Ref. [76], where it was shown how the ellipticity fluctuation distribution changes from elliptic-power distribution with large skewness to a Gaussian distribution at a large number of sources. The order of particle correlations corresponds to the order of the cumulants of the underlying flow vector fluctuation. To see the clear connection, correlations should be written in a Cartesian notation rather than polar notation (see Refs. [77,78,79] for the relation between skewness and Kurtosis of flow vector distribution to the particle correlations). Only the second-order cumulant, namely the width of the distribution, is nonvanishing for a Gaussian distribution. As a result, higher-order cumulants (skewness, kurtosis, etc.) are small for distributions close to Gaussian. These studies are done for flow amplitudes with only one harmonic, but the logic is true for more than one harmonic as well. The observed ordering of magnitudes in Fig. 1 indicates that the contribution of higher-order cumulants is smaller compared to lower ones in general, meaning the lowest-order cumulants have the dominant role in deviation from Gaussianity. A crossing between \(\langle \cos \left[ 6\left( \Psi _6-\Psi _3\right) \right] \rangle _\text {GE}\) (a three-particle correlation) and \(\langle \cos \left[ 6\left( \Psi _6-\Psi _2\right) \right] \rangle _\text {GE}\) (a four-particle correlation) is observed with centralities above 25% where the number of final state particles is lower. The same is true for \(\langle \cos \left[ 2\Psi _2+4\Psi _4-6\Psi _6\right] \rangle _\text {GE}\) (a three-particle correlation) and \(\langle \cos \left[ 2\Psi _2-3\Psi _3-4\Psi _4+5\Psi _5\right] \rangle _\text {GE}\) (a four-particle correlation). The effect of non-Gaussianity is expected to be more dominant in this centrality region since the system size is smaller and less number of sources are expected. At a finite number of sources, the actual ordering of the correlation magnitudes depends on the details of the underlying source fluctuation that needs a separate study.

Fig. 2
figure 2

Experimental data of correlations between two symmetry planes obtained with the GE compared with measurements from ATLAS [24] and ALICE [27] using the SP method. Statistical and systematic uncertainties are represented by lines and boxes, respectively

Fig. 3
figure 3

Correlations between three symmetry planes obtained with the GE compared with measurements from ATLAS [24] and ALICE [27] using the SP method. Statistical and systematic uncertainties are shown as lines and boxes, respectively

In Figs. 2 and 3 the experimental data for SPC estimated with the GE are compared with the results obtained from ATLAS [24] and ALICE [27] using the SP method. While the analysis of the SP method by ALICE used the same kinematic range as the work presented in this article, the analysis by ATLAS was performed in a wider range of \(0.5~\text {GeV}/c < p_\text {T}\) and \(|\eta |<2.5\). Despite this difference in kinematic regions, the SPC extracted by the SP method from ALICE and ATLAS agree within uncertainties. In general, the obtained data from the GE are significantly smaller than the estimates performed with the SP method for centralities larger than 10%. This difference is mainly attributed to the fact that correlations between flow amplitudes were not removed in the SP method as it was demonstrated in Ref. [15]. For the SPC \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle \) and \(\langle \cos \left[ 6\left( \Psi _6-\Psi _2\right) \right] \rangle \) shown in Fig. 2, the GE and SP method are compatible only in 0–5% centrality, while for \(\langle \cos \left[ 6\left( \Psi _6-\Psi _3\right) \right] \rangle \) the GE differs in all centrality intervals when compared to the SP method by ATLAS. For centralities larger than 5%, a clear splitting between all of the previously mentioned SPC is visible with significantly smaller values obtained by the GE. For \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle \) the experimental data of the GE are compatible with zero within the uncertainties in all considered centrality intervals. In contrast to that, the results of the SP method show a small, but non-zero value. However, the results obtained with the GE show larger uncertainties when compared to the SP for this particular SPC. Future studies with larger data sets will show whether the SPC \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle \) remains compatible with zero within uncertainties when using the GE or if a small non-zero correlation exists which cannot be resolved within the present uncertainties. In the latter case, the results of the GE will nonetheless lead to significantly smaller values than reported by the SP method.

Similarly, the experimental results of the GE and the SP method are compared to each other for SPC between three planes. The results are presented in Fig. 3. For the combinations \(\langle \cos \left[ 2\Psi _2+3\Psi _3-5\Psi _5\right] \rangle \), \(\langle \cos \left[ 2\Psi _2+4\Psi _4-6\Psi _6\right] \rangle \) and \(\langle \cos \left[ 2\Psi _2-6\Psi _3+4\Psi _4\right] \rangle \) the GE again leads to significantly smaller values than the SP method for centralities larger than 10%. For \(\langle \cos \left[ 2\Psi _2+3\Psi _3-5\Psi _5\right] \rangle \) it has to be noted that the observables previously employed by ALICE [27] and ATLAS [24] differ in the denominator. ATLAS uses a fully factorized denominator \(\langle v_2^2\rangle \langle v_3^2 \rangle \langle v_5^2\rangle \) as in the definition of the SP method (5), while the denominator in the ALICE measurement is only partially factorized \(\langle v_2^2v_3^2 \rangle \langle v_5^2\rangle \) and thus is not defined exactly as in Eq. (5). To ease the notation in Fig. 3 we still label \(\langle \cos \left[ 2\Psi _2+3\Psi _3-5\Psi _5\right] \rangle \) measured by ALICE [27] as SP method. The SPC \(\langle \cos \left[ 8\Psi _2-3\Psi _3-5\Psi _5\right] \rangle \) is the only combination where the estimates by the GE and the SP method are compatible with each other within uncertainties in all considered centralities, as the results from the SP method are already close to zero. The difference in physical interpretation between the two SPC involving \(\Psi _2\), \(\Psi _3\) and \(\Psi _5\) is discussed later.

Fig. 4
figure 4

Experimental data for correlations between two symmetry planes compared to theoretical predictions in the initial and final state obtained with T\(_\text {R}\)ENToand T\(_\text {R}\)ENTo+VISH(2+1)+UrQMD [80,81,82,83,84], respectively. For \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\) (c), the initial state predictions calculated via eccentricities and energy density cumulants, \(\langle \cos [6(\phi _2-\phi _3)]\rangle _{\text {GE}}\) and \(\langle \cos [6(\Phi _2-\Phi _3)]\rangle _{\text {GE}}\), fully overlap. Statistical (systematic) uncertainties of the ALICE data are shown as lines (boxes). The statistical uncertainties of the models are represented by the colored bands

The new measurements of SPC with the GE are compared with Monte Carlo simulations with the T\(_\text {R}\)ENTo+VISH(2+1) +UrQMD event generator [80,81,82,83,84]. In this article, the maximum a posteriori (MAP) estimation obtained in the Bayesian analysis in Ref. [7] is used for the parameters of the model. In inferring the MAP parameterization, a series of ALICE measurements (two- and four-particle correlations, charged particle multiplicities, etc.) were used as inputs into the Bayesian analysis, while the SPC are not included in these studies. Including new observables (e.g. SPC) in the Bayesian analysis can lead to an improvement in the uncertainty of the inferred parameter and resolving the discrepancies [12, 13]. If the discrepancy between model and data persists even after including new observables as input, the model itself needs to be revised.

In addition to the model predictions of final-state SPC, initial-state participant plane correlations are studied with T\(_\text {R}\)ENTo. The participant plane of order n takes the same role in the initial state as the symmetry plane in the final state. The correlations between participant planes are extracted from the initial state where flow vectors \(v_ne^{in\Psi _n}\) are replaced first by eccentricities [85], and second by cumulants of the initial energy density [38, 86, 87]. The eccentricities are defined as

$$\begin{aligned} \epsilon _n e^{in\phi _n} = -\frac{\{r^n e^{in\varphi }\}}{\{r^n\}},\quad n>1, \end{aligned}$$
(6)

where \(\{\cdots \}=\int r\text {d}r\text {d}\varphi \,\varepsilon (r,\varphi )\) stands for the average with respect to the initial energy density \(\varepsilon (r,\varphi )\) in the transverse direction and \((r,\varphi )\) are the polar coordinates in the transverse plane. Eccentricities are the moments of the initial energy density distribution. The cumulants of the initial energy density distribution, \(c_n e^{in\Phi _n}\), are obtained as a combination of eccentricities and the radial moments of the energy density, \(\{r^n\}\). In fact, cumulants are a better measure to study the deformation of a distribution close to a Gaussian. Borrowing a motivating example from Ref. [86], a Gaussian distribution \(e^{-x^2/2\sigma _x^2-y^2/2\sigma _y^2}\) has infinitely many non-vanishing moments, while only its second order cumulants are non-zero. Following the convention of Ref. [86], the first two cumulants and eccentricities are equivalent, \(c_n e^{in\Phi _n}= \epsilon _n e^{in\phi _n}\) for \(n=2,3\). Higher order cumulants have non-trivial relations to eccentricities. Here, only the fourth harmonic is shown as an example:

$$\begin{aligned} c_4 e^{i4\Phi _4}=\epsilon _4 e^{i4\phi _4}+3\left( \frac{\{r^2\}^2}{\{r^4\}}\right) \epsilon _2^2 e^{i4\phi _2}. \end{aligned}$$
(7)

More details can be found in Refs. [38, 86, 87].

The comparison with initial and final state SPC demonstrates how much of the observed correlation is inherited from the initial state. This is due to the linear and non-linear response of the medium [88]. For the second- and third-order anisotropies, the linear response is expected to dominate i.e. \(v_2e^{i2\Psi _2} = \omega _2 c_2e^{i2\Phi _2}\) and \(v_3e^{i3\Psi _3} = \omega _3 c_3e^{i3\Phi _3}\), especially in central and semicentral collisions [73, 85, 89]. The \(\omega _i\) describe the linear hydrodynamic coupling constants. For higher orders, non-linear contributions will play a significant role, e.g. in case of the fourth order as

$$\begin{aligned} v_4e^{i4\Psi _4} = \omega _4 c_4e^{i4\Phi _4} + \omega _{422}c_2^2e^{i4\Phi _2}, \end{aligned}$$
(8)

where \(\omega _{422}\) is the non-linear coupling between the second- and fourth-order anisotropies [38, 86, 87]. As an example of how this impacts the SPC, one can build the quantity \(v_2^2v_4e^{i4(\Psi _4-\Psi _2)}\). The real part of its phase corresponds to the SPC \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle \). Using the linear and non-linear response, one can translate this into the initial state as:

$$\begin{aligned} v_2^2v_4e^{i4(\Psi _4-\Psi _2)} = \omega _2 \omega _4 c_2^2 c_4 e^{i4\left( \Phi _4-\Phi _2\right) } + \omega _{422} \omega _2^2 c_2^2. \end{aligned}$$
(9)

The latter equation shows that the initial and final state SPC \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle \) will be equal to each other in the limit of pure linear response, while they will deviate in case of a non-zero non-linear coupling \(\omega _{422}\).

Figure 4 shows the comparison between the model calculations and the experimental data for the correlation between two symmetry planes. For the SPC \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle _\text {GE}\) the initial-state participant plane correlations given by the energy density cumulants overlap with the final-state prediction of the SPC up to 10% in centrality, indicating a vanishing non-linear coupling in the regime dominated by flow fluctuations. For higher centralities, the two curves increasingly deviate, showing the presence of a non-zero non-linear coupling between the second- and fourth-order flow vectors. In particular, it can be observed that the final-state prediction increasingly deviates from the data with increasing centrality. This is expected, since in this regime strong correlations between the second and fourth harmonics can originate from the initial ellipsoidal geometry. The non-linear coupling constant between initial-state ellipticity and \(v_4 e^{i 4\Psi _4}\) from T\(_\text {R}\)ENTo to iEBE-VISHNU was studied in Ref. [90]. It was demonstrated that this coupling is very small up to 10% and it is positive up to 70% centrality. For \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\) the experimental data show a flat centrality behavior and are compatible with zero within the uncertainties. The model predicts small values for \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\) with a flat centrality behavior although the predictions for the final state slightly overestimate the data. However, the initial-state correlations decrease to more negative values for centralities above 30%. For \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\), a linear response is expected to be an accurate approximation for harmonics \(n=2, 3\) in central and semicentral collisions [73, 85, 89]. As such, one would expect that the initial state participant plane correlations should be the same as the final state symmetry plane correlations between the second- and third-order harmonics. One possible explanation is that higher-order terms beyond linear response are responsible for decreasing the correlation during the hydrodynamic evolution, and the final value accidentally lands on very small numbers. In this respect, more rigorous study is needed in the future. While the model captures the qualitative behavior of the experimental signal in the two previous cases, a quantitative agreement is not observed in every case, particularly not for SPC involving \(\Psi _6\) or correlations between four symmetry planes. For the SPC \(\langle \cos \left[ 6\left( \Psi _6-\Psi _2\right) \right] \rangle _\text {GE}\) and \(\langle \cos \left[ 6\left( \Psi _6-\Psi _3\right) \right] \rangle _\text {GE}\), a large deviation between the data and the final-state model prediction can be observed. This deviation could be related to the complex dynamics of the sixth-order symmetry plane, which involves multiple non-linear responses to the lower order symmetry planes \(\Psi _2\), \(\Psi _3\) and \(\Psi _4\). For the SPC \(\langle \cos \left[ 6\left( \Psi _6-\Psi _3\right) \right] \rangle _\text {GE}\) it is in particular interesting that the initial-state correlation becomes stronger with increasing centrality, while the final-state correlation becomes weaker with increasing centrality in the model. In contrast to that, the experimental data shows only a very weak centrality dependence.

The sign change between correlations obtained from eccentricities (green short-dashed curves) and initial energy density cumulants (blue long-dashed curves) in Fig. 4 was pointed out in Ref. [87]. The actual shape of the initial state is captured by \(c_ne^{in\Phi _n}\). As a result, the linear hydrodynamic response approximation is more accurate when employing cumulants. As an example, SPC \(\langle \cos \left[ 4\left( \Psi _4-\Psi _2\right) \right] \rangle _\text {GE}\) in Fig. 4 panel (a) is considered. Referring to Eq. (7), the difference between \(\epsilon _4 e^{in\phi _4}\) and \(c_4 e^{in\Phi _4}\) is proportional to \(\{r^2\}^2/\{r^4\}\), which merely depends on the radial shape of the initial state. The numerator of the correlation \(\langle \cos (4\phi _2-4\phi _4)\rangle _{\text {GE}}\) is proportional to the real part of \(\langle \epsilon _2^2\epsilon _4 e^{4i\phi _2-4i\phi _4}\rangle \). Substituting from Eq. (7) (\(c_2e^{2i\Phi _2}\) and \(\epsilon _2e^{2i\phi _2}\) are identical), one finds \(\langle \epsilon _2^2\epsilon _4 e^{4i\phi _2-4i\phi _4}\rangle = \langle c_2^2c_4 e^{4i\Phi _2-4i\Phi _4}\rangle -3\langle (\{r^2\}^2/\{r^4\}) c_2^4\rangle \). In case the second term on the right is small, SPC of cumulants are expected to be close to those calculated from eccentricities. However, noting that \(3\langle (\{r^2\}^2/\{r^4\}) c_2^4\rangle \) is real and positive, the contribution of this term is such that it leads to a negative sign for \(\langle \epsilon _2^2\epsilon _4 e^{4i\phi _2-4i\phi _4}\rangle \). Therefore, one concludes that, regarding the hydrodynamic response, using energy-density cumulants for the initial state is more appropriate. A more detailed study can be found in Ref. [90].

Fig. 5
figure 5

Experimental data for correlations between three (ad) and four (e) symmetry planes compared with theoretical predictions in the initial and final state obtained with T\(_\text {R}\)ENToand T\(_\text {R}\)ENTo+VISH(2+1)+UrQMD [80,81,82,83,84], respectively. Statistical (systematic) uncertainties of the ALICE data are shown as lines (boxes). The statistical uncertainties of the models are represented by the colored bands

In Fig. 5, the same comparison of model predictions with respect to the experimental data was performed for SPC involving three or four symmetry planes. For the SPC \(\langle \cos \left[ 2\Psi _2+3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\), a non-zero correlation signal was extracted, which existed to a larger extent already in the initial state. The final-state predictions of T\(_\text {R}\)ENTo+iEBE-VISHNU are in good agreement with the experimental data. Further considering the same three symmetry planes \(\Psi _2\), \(\Psi _3\) and \(\Psi _5\), the combination \(\langle \cos \left[ 8\Psi _2-3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\) results in an experimental signal that is compatible with zero within uncertainties for all reported centrality intervals. The deviation between \(\langle \cos \left[ 2\Psi _2+3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\) and \(\langle \cos \left[ 8\Psi _2-3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\) can be attributed to different contributions from the initial-state in the non-linear response of the two observables. While the non-linear response term for \(\langle \cos \left[ 2\Psi _2+3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\) does not contain contributions from any participant plane correlations, the non-linear part of \(\langle \cos \left[ 8\Psi _2-3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\) picks up such an additional correlation from the initial state. In particular, the SPC \(\langle \cos \left[ 8\Psi _2-3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\) has a non-linear coupling to the initial-state correlation between the second- and third-order participant planes. Similar to previously presented examples of SPC that involve the sixth-order symmetry plane \(\Psi _6\), the final-state model prediction for \(\langle \cos \left[ 2\Psi _2+4\Psi _4-6\Psi _6\right] \rangle _\text {GE}\) shows a large deviation compared to the measurements. This deviation is again attributed to the complexity of \(\Psi _6\), which makes it particularly sensitive to the model parameters. One finds that \(\langle \cos \left[ 2\Psi _2+3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\) has the strongest signal in data while \(\langle \cos \left[ 2\Psi _2+4\Psi _4-6\Psi _6\right] \rangle _\text {GE}\) is larger than \(\langle \cos \left[ 2\Psi _2+3\Psi _3-5\Psi _5\right] \rangle _\text {GE}\) in the model. The results of \(\langle \cos \left[ 2\Psi _2-6\Psi _3+4\Psi _4\right] \rangle _\text {GE}\) are the only measurement with a negative signal for the SPC. While the model predicts this behavior as well, it can be observed that the initial-state correlations are strictly non-negative. Thus, the sign-change between the initial and final state can be linked to the hydrodynamic evolution of the system. Lastly, Fig. 5 shows the comparison with the first experimental measurement of correlations between four symmetry planes. While the final-state model prediction captures the qualitative behavior, a quantitative agreement between experimental data and model is not observed. A recent Bayesian analysis of the QGP hydrodynamic properties [13] has shown that higher-order correlations as well as measurements involving higher order flow harmonics are more sensitive to changes in the QGP parameters. Thus, in particular the measurements of SPC including \(\Psi _6\) and the correlations between four symmetry planes are expected to give more stringent constraints in future analyses of this kind.

4 Conclusion

Utilizing the recently introduced Gaussian Estimator, the first measurements of symmetry plane correlations, which are not influenced by correlations between different flow amplitudes, are presented in Pb–Pb collisions at \(\sqrt{s_\text {NN}}=2.76\) TeV. Correlations between two, three and, for the first time, four symmetry planes were shown. The data show a clear order for the different SPC, which can be related to the cumulants of the underlying flow vector fluctuations. The measurements using the new GE show significantly smaller symmetry plane correlations than previously reported by the SP method. This observation is in qualitative agreement with the study in Ref. [15] which reported the bias of the SP method to larger values due to the influence of correlations between the flow amplitudes. In contrast to the SP method, the results of \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\) are consistent with zero with the current uncertainties. Within the uncertainties, this shows that \(\Psi _2\) and \(\Psi _3\) are fully uncorrelated, which was qualitatively reported by a previous ALICE study [23]. Future studies using Run 2 Pb–Pb data as well as the upcoming Run 3 Pb–Pb campaign will show whether a small non-zero correlation exists between \(\Psi _2\) and \(\Psi _3\), which cannot be resolved within the present uncertainties. Employing state-of-the-art hydrodynamic model calculations, one could see that the predictions and the measurements are not in quantitative agreement for all SPC. The most significant deviations are observed in the correlations \(\langle \cos \left[ 6\left( \Psi _6-\Psi _3\right) \right] \rangle _\text {GE}\), \(\langle \cos \left[ 6\left( \Psi _6-\Psi _2\right) \right] \rangle _\text {GE}\), \(\langle \cos \left[ 2\Psi _2+4\Psi _4-6\Psi _6\right] \rangle _\text {GE}\) and \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\)FourFive. Future studies have to address how the initial state correlations between the second- and third-order participant planes are suppressed in the final state SPC \(\langle \cos \left[ 6\left( \Psi _2-\Psi _3\right) \right] \rangle _\text {GE}\) as one would expect linear response to be a good approximation for the second- and third-order flow harmonics. Since the measured SPC contain independent information about flow vector fluctuations, they will provide useful inputs for future Bayesian analyses aiming at extracting the properties of the QGP.