1 Introduction

In heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC), the fluctuations of the positions of nucleons in the overlap region are found to play an important role in controlling the shape of the initial geometry of the created matter, which subsequently controls the azimuthal anisotropy of the particles in the final state [13]. The shape of the geometry in azimuth can be characterized by a set of multi-pole components (also known as “eccentricities”) at different angular scale, calculated from the participants and binary collisions at (r,ϕ) [4]:

$$\begin{aligned} \epsilon_n = \frac{\sqrt{\langle r^2\cos n\phi\rangle^2+\langle r^2\sin n\phi\rangle^2}}{\langle r^2\rangle}, \end{aligned}$$
(1)

with a weight of δ=0.14 for binary collisions and (1−δ)/2=0.43 for participants [5], where (r,ϕ) are calculated relative to the weighted center of gravity. The orientations of the minor and major axes for each moment n are given by

$$\begin{aligned} \varPhi_n=\frac{\mathrm{atan2}(\langle r^2\sin n\phi\rangle,\langle r^2\cos n\phi\rangle)}{n}+\frac{\pi}{n} \end{aligned}$$
(2)

and

$$\begin{aligned} \varPhi_n^{*} = \varPhi_n + \pi/n \end{aligned}$$
(3)

respectively. The minor axis direction Φ n is also known as the nth-order participant plane (PP). The values of ϵ n and Φ n can be calculated easily using simple geometric models such as Monte Carlo Glauber code from [6].

When fluctuations are small and linearized hydrodynamics is applicable, each ϵ n is expected to independently drive the corresponding nth-order anisotropic flow v n along Φ n  [4]. In this case, one may rely on a simple Glauber model calculation to estimate the correlations between anisotropic flows of different order.Footnote 1 Previous studies [3, 711] show that significant correlations can exist between Φ 2 and Φ 4 due to the almond shape of the average collision geometry. However, correlations involving odd planes for n>2 are found to be generally weak except in very peripheral collisions, e.g. between Φ 2 and Φ 3 or between Φ 2 and Φ 5 [7]. Experimental studies support a strong correlation between Φ 2 and Φ 4 [12, 13], but a weak correlation between Φ 2 and Φ 3 [14, 15]. The correlations among three planes of different order have also been investigated recently [3, 10, 11], such as Φ 1+2Φ 2−3Φ 3 and Φ 1−4Φ 2+3Φ 3; they are argued to contain strong correlations between the dipole asymmetry and the triangularity. Here we propose an alternative experimental method for measuring these correlations. The expected performance of this method is evaluated based on the correlation signals from Glauber model.

2 Method

The nth-order flow has n-fold symmetry in azimuth, and the correlations between flow directions Φ n and Φ m are completely described by the differential distribution dN evts/(d(k(Φ n Φ m ))), with k being the least common multiple of n and m, i.e. k=LCM(n,m). This distribution should be an even function due to symmetry, and can be expanded into a Fourier series:

$$\begin{aligned} &{ \frac{dN_{\mathrm{evts}}}{d (k(\varPhi_n-\varPhi_m) )}\propto 1+2\sum_{j=1}^{\infty} V_{n,m}^j \cos jk(\varPhi_n- \varPhi_m) } \end{aligned}$$
(4)
$$\begin{aligned} &{V_{n,m}^j = \bigl\langle\cos jk( \varPhi_n-\varPhi_m)\bigr\rangle} \end{aligned}$$
(5)

In a real experiment, the underlying true event plane directions Φ n and Φ m are unattainable due to limited detector acceptance and finite multiplicity. They are approximated by the measured event plane angle Ψ n and Ψ m , calculated based on the azimuthal distribution of particles in the detector acceptance. The coefficients 〈cosjk(Φ n Φ m )〉 can be obtained by calculating the raw coefficients 〈cosjk(Ψ n Ψ m )〉, followed by a simple correction for finite event plane resolution:

$$\begin{aligned} &{V_{n,m}^{j} = \frac{\langle\cos jk(\varPsi_n-\varPsi_m)\rangle}{\mathrm {Res}\{jk\varPsi_n\}\mathrm{Res}\{jk\varPsi_m\}}} \end{aligned}$$
(6)
$$\begin{aligned} &{\mathrm{Res}\{jk\varPsi_n\} = \bigl\langle\cos jk(\varPsi_n- \varPhi_n)\bigr\rangle} \end{aligned}$$
(7)

The resolution factor Res{jkΨ n } can be determined using the standard two-subevent or three-subevent methods [17]. To avoid auto-correlations, the Φ n and Φ m should be measured using sub-events covering different η ranges, preferably with a gap in between.

Interestingly, some of the two plane correlators defined by Eq. (5) are related to the so called mixed harmonics, referring to v ln measured in Φ n event plane for integer l≥2, denoted as v ln {Φ n } (see Ref. [17]). For example, it is straightforward to show that \(V_{2,4}^{1}\) is simply the ratio of the integral v 4 measured in the Φ 2 plane (v 4{Φ 2}) to the integral v 4 measured in the Φ 4 plane (v 4{Φ 4}): 〈cos4(Φ 2Φ 4)〉=v 4{Φ 2}/v 4{Φ 4}. More generally, one has:

$$\begin{aligned} \bigl\langle\cos m(\varPhi_n-\varPhi_m) \bigr\rangle=\frac{v_m\{\varPhi_n\}}{v_m\{\varPhi_m\}},\quad m \bmod n =0 \end{aligned}$$
(8)

Additional examples include 〈cos6(Φ 2Φ 6)〉=v 6{Φ 2}/v 6{Φ 6} and 〈cos6(Φ 3Φ 6)〉=v 6{Φ 3}/v 6{Φ 6}.

The method described above can be generalized to correlations involving three or more event planes. As pointed out in Ref. [10], the correlations that can be measured experimentally, involve combination of l planes of different order: c 1 Φ 1+2c 2 Φ 2+⋯+lc l Φ l with c 1+2c 2+⋯+lc l =0. The correlations involving three planes of different order, e.g. Φ 1, Φ 2 and Φ 3, have the following form:

$$\begin{aligned} c_1\varPhi_{1}+2c_2 \varPhi_{2}+3c_3\varPhi_{3} =& c_22 (\varPhi_2-\varPhi_1 )+c_33 ( \varPhi_3-\varPhi_1 ), \\ =&c_2\varPhi_{2,1}+c_3\varPhi_{3,1} \end{aligned}$$
(9)

where we use the constraint c 1+2c 2+3c 3=0 and we adopt the short-hand notations: Φ n,m =k(Φ n Φ m ),Ψ n,m =k(Ψ n Ψ m ). This type of correlations can be generally determined from the underlying 2-D distribution in (Φ 2,1,Φ 3,1) via a similar Fourier expansion approach:

$$\begin{aligned} \frac{d^2N_{\mathrm{evts}}}{d\varPhi_{2,1}d\varPhi_{3,1}} \propto&1+2\sum_{j=1}^{\infty} \bigl[ V_{1,2}^{j} \cos j\varPhi_{2,1}+V_{1,3}^{j} \cos j\varPhi_{3,1} \bigr] \\ &{}+2\sum_{i,j=1}^{\infty }V_{1,2,3}^{i,\pm j} \cos (i\varPhi_{2,1}\pm j\varPhi_{3,1} ) \end{aligned}$$
(10)

This series is expected to converge quickly for non-peripheral collisions. Therefore, only the terms for i,j≤3 are required (see Fig. 5). The coefficients are:

$$\begin{aligned} V_{1,2,3}^{i,\pm j} =& \bigl\langle\cos (i\varPhi_{2,1}\pm j\varPhi_{3,1} )\bigr\rangle \\ =&\langle\cos i\varPhi_{2,1}\cos j \varPhi_{3,1}\rangle\mp \langle\sin i\varPhi_{2,1}\sin j \varPhi_{3,1}\rangle \end{aligned}$$
(11)

Under this notation, the two-plane correlator can be treated as special case: \(V_{1,2,3}^{i,0} = V_{1,2}^{i}, V_{1,2,3}^{0,j} = V_{1,3}^{j}\) and \(V_{1,2,3}^{3j,-2j} = V_{3,2}^{j}\). The average of the sine term in Eq. (11) may not be zero, if the fluctuations of Φ 2 and Φ 3 relative to Φ 1 are correlated, i.e. Φ 2 and Φ 3 prefer to appear simultaneously to one side of Φ 1. It represents a non-trivial correlation that is of great interest for understanding the nature of the fluctuations (see also [3]).

A similar resolution correction procedure can be used to connect the measured correlated with the corrected one:

$$\begin{aligned} V_{1,2,3}^{i,\pm j} = \frac{\langle\cos (i\varPsi _{2,1}\pm j\varPsi_{3,1} )\rangle}{\mathrm{Res}\{|2i\pm3j|\varPsi _{1}\}\mathrm{Res}\{2i\varPsi_{2}\}\mathrm{Res}\{3j\varPsi_{3}\}} \end{aligned}$$
(12)

where we have assumed that Ψ n is distributed randomly around Φ n , such that 〈sinjn(Ψ n Φ n )〉=0. Again, the Ψ 1, Ψ 2 and Ψ 3 should be calculated from subevents covering different η acceptances to avoid auto-correlations.

In [3], Teaney and Yan proposed to study the correlator cos(Φ 1+2Φ 2−3Φ 3) and cos(Φ 1−4Φ 2+3Φ 3). In our notations, they correspond to the cosine average of the 2-D distribution (Φ 2,1,Φ 3,1) projected along the direction (i,j)=(1,−1) and (2,−1), respectively. Our framework provides a natural way to visualize and systematize the study of these correlators.

Other triple plane correlators can be similarly analyzed, the first few are

$$\begin{aligned} &{c_1\varPhi_{1}+2c_2 \varPhi_{2}+4c_4\varPhi_{4} = c_2 \varPhi_{2,1}+c_4\varPhi_{4,1}} \end{aligned}$$
(13)
$$\begin{aligned} &{c_1\varPhi_{1}+3c_3 \varPhi_{3}+4c_4\varPhi_{4} =c_3 \varPhi_{3,1}+c_4\varPhi_{4,1}} \end{aligned}$$
(14)
$$\begin{aligned} &{2c_2\varPhi_{2}+3c_3 \varPhi_{3}+4c_4\varPhi_{4}=\frac{c_3}{2} \varPhi_{3,2}+c_4\varPhi _{4,2}} \end{aligned}$$
(15)
$$\begin{aligned} &{c_1\varPhi_{1}+2c_2 \varPhi_{2}+5c_5\varPhi_{5}=c_2 \varPhi_{2,1}+c_5\varPhi_{5,1}} \end{aligned}$$
(16)
$$\begin{aligned} &{c_1\varPhi_{1}+3c_2 \varPhi_{3}+5c_5\varPhi_{5}=c_3 \varPhi_{3,1}+c_5\varPhi_{5,1}} \end{aligned}$$
(17)

Note that c 3/2 in Eq. (15) is an integer by the requirement 2c 2+3c 3+4c 4=0. These correlators can be uniquely identified with one of the Fourier coefficients in the double differential distributions similar to Eq. (10). However, the expression of triple plane correlator in terms of the correlation between two di-plane correlators is not always possible, for example 〈cos(2Φ 2+3Φ 3−5Φ 5)〉. In this case, it can be regarded as a sum of the Fourier coefficients for triple differential distributions d 2 N evts/( 3,2 5,1 5,3) that contribute to cos(2Φ 2+3Φ 3−5Φ 5).

The measurement of correlations involving two or more event planes are feasible at the LHC due to the large detector coverage in η (−5<η<5 in ATLAS and CMS), and excellent reaction plane resolution [18, 19]. This allows the choice of many non-overlapping sub-events, each with very good η coverage for these multi-plane correlation measurements. This works as long as the true event plane angle does not rotate as a function of pseudorapidity and so far there are no experimental evidences for this rotation.

The coefficients \(V_{n,m}^{j}\) or \(V_{n,m,h}^{i,\pm j}\) are related to the previously proposed multi-particle correlators from Refs. [10, 11]. That approach effectively applies a |c n |-particle weight \({v_{n}^{\{c_{n}\}}\equiv(v_{n})_{1} (v_{n})_{2}\cdots (v_{n})_{|c_{n}|}}\) if the angle nc n Φ n appears in the correlator. This weight maximizes the correlation signal and reduce the contribution for events which have small v n . For cos(c n n +c m m ) and cos(c n n +c m m +c h h ), the weights are \({v_{n}^{\{|c_{n}|\}}v_{m}^{\{|c_{m}|\}}}\) and \({v_{n}^{\{|c_{n}|\}}v_{m}^{\{|c_{m}|\}}v_{h}^{\{|c_{h}|\}}}\), respectively; they are then divided by \(\langle v_{n}^{\{|c_{n}|\}}v_{m}^{\{|c_{m}|\}}\rangle\) and \(\langle v_{n}^{\{|c_{n}|\}}v_{m}^{\{|c_{m}|\}}v_{h}^{\{|c_{h}|\}}\rangle\) to obtain the true correlations. Note that the weighting procedure can also amplify contributions from the tail of the ϵ n distribution, especially for large values of c n , this may complicate the mapping from the measurement to correlations between ϵ n .

In contrast, all events have the same weight in our approach (including those with small ϵ n values unfortunately). In our opinion, the two methods are complimentary to each other. In fact, it is possible to construct some hybrid correlators that involves azimuthal angle of both event planes and single particles. For example, one can consider the following mixed correlator between a+b particles and event planes Ψ n ,Ψ m :

$$\begin{aligned} &{\bigl\langle\cos(nc_n\varPhi_n-mc_m \varPhi_m)\bigr\rangle_{v_n^{\{a\} }v_m^{\{b\}}}} \\ &{\quad= \frac{\langle\cos (\sum\phi _{n,m}^{a,b}+n(c_n-a)\varPsi_n+m(c_m-b)\varPsi_m )\rangle}{\langle v_n^{\{a\}} v_m^{\{b\}}\rangle\mathrm{Res}\{n(c_n-a)\varPsi_n\}\mathrm {Res}\{m(c_m-b)\varPsi_m\}}} \end{aligned}$$
(18)
$$\begin{aligned} &{\sum\phi_{n,m}^{a,b}= n( \phi_{1}+\cdots+\phi_{a})+m(\phi _{a+1}+\cdots+ \phi_{a+b})} \end{aligned}$$
(19)

where nc l mc m =0, ϕ 1,…,ϕ a+b are azimuthal angles of a+b particles, and subscript \(v_{n}^{\{a\}}v_{m}^{\{b\}}\) indicates the weighting factor introduced by those particle multiplets. Similar formula can be generalized to correlations involving more than two event planes. Three interesting examples are:

$$\begin{aligned} &{\bigl\langle\cos6(\varPhi_2- \varPhi_3)\bigr\rangle_{v_3^{\{2\}}}= \frac{\langle\cos (3\phi_1+3\phi_2-6\varPsi_2 )\rangle}{\langle (v_3)_1(v_3)_2\rangle\mathrm{Res}\{6\varPsi_2\}}} \end{aligned}$$
(20)
$$\begin{aligned} &{\bigl\langle\cos2(\varPhi_1- \varPhi_2)\bigr\rangle_{(wv_1)^{\{2\}}}= \frac{\langle\cos (\phi_1+\phi_2-2\varPsi_2 )\rangle}{\langle (wv_1)_1(wv_1)_2\rangle\mathrm{Res}\{2\varPsi_2\}}} \end{aligned}$$
(21)
$$\begin{aligned} &{\bigl\langle\cos(\varPhi_1+2 \varPhi_2-3\varPhi_3)\bigr\rangle_{wv_1}= \frac{\langle\cos (\phi_1+2\varPsi_2-3\varPsi_3 )\rangle}{\langle wv_1\rangle\mathrm{Res}\{2\varPsi_2\}\mathrm{Res}\{3\varPsi_3\}}} \\ \end{aligned}$$
(22)

where \(w( p_{\mathrm {T}} ,\eta)= p_{\mathrm {T}} -\langle p_{\mathrm {T}} ^{2}\rangle(\eta)/\langle p_{\mathrm {T}} \rangle(\eta )\) is the p T and η dependent weight (w is rapidity-even for Au+Au collisions) that designed to suppress global momentum conservation effects and maximizing the v 1 signal (or effectively increasing the resolution for Ψ 1) [10, 15, 20, 21]. Even though terms related to v 1 (ϕ 1 and/or ϕ 2) appear in Eqs. (21)–(22), the global momentum conservation effects are expected to be either negligible Eq. (21) or absent Eq. (22) for these correlators [15, 22].

The hybrid correlators are useful in real experiments where detector subsystems have limited geometrical acceptance, finite granularity (so can’t distinguish individual particles), limited p T reach or no p T information at all. In this case, it is straightforward to calibrate the event plane measurement, while the calibration procedure could be more involved for multi-particle correlations [23].

3 Expected behavior from Glauber and CGC models

3.1 Correlation between two planes

Similar to [7], we use a simple Glauber model [24] to estimate the level of the correlations between nth- and mth-order participant plane. About 2.5 Million Au+Au collisions are simulated, where each Au ion is populated randomly with nucleons with a hard-core of 0.3 fm in radii, according to the Woods–Saxon distribution with a radius of 6.38 fm and diffuseness of 0.535 fm. A nucleon-nucleon cross-section of σ=42 mb is used. The Φ n and ϵ n are defined as a combination of participants and binary collisions in the transverse plane as mentioned in the introduction. However, instead of using the minor axes of ϵ n as the proxy for the true event planes, we actually calculate the correlations between the major axes, which are related to the former by a simple phase shift δ n,m :

$$\begin{aligned} &{k\bigl(\varPhi_n^{*}-\varPhi_m^{*} \bigr) = k(\varPhi_n-\varPhi_m) +\delta_{n,m}} \end{aligned}$$
(23)
$$\begin{aligned} &{\delta_{n,m} = k(1/n-1/m)\pi} \end{aligned}$$
(24)

The corresponding Fourier coefficients are denoted as \(V_{n,m}^{j*}=\langle\cos k(\varPhi_{n}^{*}-\varPhi_{m}^{*})\rangle\), and are related to \(V_{n,m}^{j}\) as

$$\begin{aligned} V_{n,m}^{j *}=(\cos\delta)^j V_{n,m}^{j} \end{aligned}$$
(25)

The reason for doing this is a simple matter of convenience: the correlations between major axes are found to be almost always positive in the Glauber model, hence using major axes simplifies the presentation. It is interesting to note that the phase shift, when folded to [0,2π], is (δ n,m mod2π)=0 or π. The latter case leads to a sign flip: \(V_{n,m}^{j *}=-V_{n,m}^{j}\).

The top panels of Fig. 1 show the \(4(\varPhi_{2}^{*}-\varPhi_{4}^{*})\) correlations predicted by the Monte Carlo Glauber model. The correlation is weak in central collisions but is quite strong in peripheral collisions. This is understood [7, 9] due to a detailed interplay between the fluctuation and average shape for the collision geometry: the central collisions are fluctuation-dominated, hence Φ n are largely uncorrelated, while the peripheral collisions are dominated by the average geometry which has ϵ 2n components that are aligned [25]. Figure 1 also shows that the first order component captures most of the shape information in central collisions. In contrast, many components are needed to describe the tight correlation in peripheral collisions. This behavior is generally true in the Glauber model: whenever the first term \(V_{n,m}^{1 *}\) is large, more higher-order terms are needed to describe the full distribution. Note that if the participant planes are used instead, the distribution 4(Φ 2Φ 4) would show an anti-correlation: the distributions have their minima at 0, and the sign of \(V_{2,4}^{j}\) alternates between positive and negative: \(V_{2,4}^{j}=(-1)^{j}V_{2,4}^{j*}\).

Fig. 1
figure 1

(Top panels) The distribution of the angle difference between major axes for ϵ 2 and ϵ 4 in two centrality intervals, with the thin (thick) lines indicating the contribution from the first six harmonics (the sum). (Bottom panel) The Fourier spectra (Color figure online)

The calculations are extended for all types of correlations for k up to 16. Of course, additional correlations can also be calculated but the resolution is expected to deteriorate quickly for large values of k. The centrality dependence of these correlations, characterized by the first Fourier coefficient \(V_{n,m}^{1*}\), are shown in Fig. 2, where the centrality is characterized by number of participating nucleons (N part). In most cases, the correlations are found to be either consistent with zero or positive (except for \(V_{2,6}^{1*}\) in mid-central collisions). In particular, strong correlations are observed for \(4(\varPhi_{2}^{*}-\varPhi_{4}^{*})\), \(6(\varPhi _{3}^{*}-\varPhi_{6}^{*})\) and \(6(\varPhi_{2}^{*}-\varPhi_{6}^{*})\); they are presumably associated with the average geometry. The correlations are small in central collisions and over the full range for other choices of n and m, suggesting that the correlations are generally weak when they are dominated by fluctuations. We would like to draw the reader’s attention to the bottom panels, which suggest that there are significant correlations between Φ 1 and all other higher-order PPs. Recently, significant dipolar flow v 1 has been observed in Pb–Pb collisions at the LHC by the ATLAS Collaboration [18] and a theoretical group [26] based on the ALICE data [27]; large dipolar flow is also predicted in hydrodynamic [28] or transport models [21]. Therefore, it is reasonable to assume that the correlations between dipolar flow and higher-order flow is large and measurable, as long as one can find a way to determine Φ 1 without the bias of the global momentum conservation effect. One possibility might be to use the modified event plane method from Ref. [20].

Fig. 2
figure 2

Centrality dependence of first Fourier coefficient of the correlation \(V_{n,m}^{1 *}=\langle\cos k(\varPsi _{n}^{*}-\varPsi_{m}^{*})\rangle\) for various choices of n and m (Color figure online)

To get a feeling on how many \(V_{n,m}^{j*}\) terms are needed to exhaust the information encoded in distribution \(k(\varPhi_{n}^{*}-\varPhi_{m}^{*})\), in Fig. 3 we show the centrality distribution of \(V_{n,m}^{j*}\) for several values of j for j>1. Most of correlations are exhausted by including the j=1–5, except for a few cases at N part<100, such as \(4(\varPhi_{2}^{*}-\varPhi_{4}^{*})\) and \(3(\varPhi_{1}^{*}-\varPhi_{3}^{*})\).

Fig. 3
figure 3

Centrality dependence of higher-order Fourier coefficients of the correlation \(V_{n,m}^{j *}=\langle\cos jk(\varPsi_{n}^{*}-\varPsi_{m}^{*})\rangle\) for various choices of n and m for j=2–5 (from left column to right column) (Color figure online)

The results presented in Figs. 13 are obtained using the r 2 weighting (i.e. Eq. (2)) and Glauber model. Alternatively we have calculated the Φ n using the r n weighting for n>1 and r 3 weighing for n=1 [3]; we also repeated the same calculation using a CGC (Color Glass Condensate) geometry [29, 30] for both the r 2 and r n weighting. The results for these three cases are shown in Fig. 4 for j=1. The main observations are qualitatively similar to Fig. 2. However, there are some interesting changes in the magnitudes of some correlation values: the use of CGC model increases the correlation for (n,m)=(2,4) and (2,6) but decreases the correlation for (n,m)=(3,6) in mid-central collisions; the use of r n weighting also in general increases \(V_{n,m}^{j*}\) and affects the relative magnitude of \(V_{n,m}^{j*}\) between (n,m)=(2,4) and (3,6) in central collisions.

Fig. 4
figure 4

Centrality dependence of first-order Fourier coefficients of the correlation \(V_{n,m}^{1*}=\langle\cos k(\varPsi _{n}^{*}-\varPsi_{m}^{*})\rangle\) for various choices of n and m for different ways of calculating ϵ n (from left column to right column) (Color figure online)

Due to the phase shift between the two types of correlations given by Eq. (25), many positive \(V_{n,m}^{j*}\) values in Figs. 2 and 3 would imply the corresponding \(V_{n,m}^{j}\) values are negative. This happens for odd j and (n,m)=(2,3), (3,6), (2,5), (1,2), (1,4) and (1,6). If the nth-order flow direction align with Φ n as predicted by the Glauber model, one should expect the signs of the correlators between the experimental event plane in Eq. (6) to exhibit very interesting dependence on choice of (n,m). On the other hand, if the dynamic mixing between flow of different orders is important [9], then this dependence could be strongly distorted. Therefore, direct measurements of the correlations between the experimentally determined event planes of different order can help to resolve this issue.

3.2 Correlation between three planes

It is straightforward to carry out the study of correlations between three planes. The “*” notation is again used to indicate the plane calculated using the major axes. The top panels of Fig. 5 show the 2-D normalized distribution \(d^{2}N_{\mathrm{evts}}/(d\varPhi _{2,1}^{*}d\varPhi_{3,1}^{*})\) in 40–50 % centrality interval; the corresponding \(V_{1,2,3}^{i,j*}\) coefficients are shown in the bottom panels. The coefficients along i=0 or j=0 simply reflect two plane correlators, \(V_{1,2}^{j*}\) and \(V_{1,3}^{i*}\), respectively. The interesting coefficients are those for i,j≠0. A tight diagonal correlation in the top panels can be identified with a large (i,j)=(1,−1) component, which corresponds to \(\langle\cos(\varPhi _{1}^{*}+2\varPhi_{2}^{*}-3\varPhi_{3}^{*})\rangle\). This correlation is very weak in central collision and increases gradually towards peripheral collisions (see Figs. 1114), similar to the finding in [3] (there is a sign difference due to the use of major axes here). This correlation is also observed to be generally bigger for r n weighting and CGC geometry. Hence a precise determination of this correlator could allow us to distinguish between different models for initial geometry. Sizable coefficients are also observed for (i,j)=(1,−2), (2,−2) and (1,1), corresponding to \(\langle \cos(4\varPhi_{1}^{*}+2\varPhi_{2}^{*}-6\varPhi_{3}^{*})\rangle\), \(\langle\cos(2\varPhi _{1}^{*}+4\varPhi_{2}^{*}-6\varPhi_{3}^{*})\rangle\) and \(\langle\cos(5\varPhi _{1}^{*}-2\varPhi_{2}^{*}-3\varPhi_{3}^{*})\rangle\), respectively. Also note that the coefficients for (i,j)=(2,−1) and (3,−2), corresponding to \(\langle\cos(\varPhi_{1}^{*}-4\varPhi_{2}^{*}+3\varPhi_{3}^{*})\rangle\) and \(\langle \cos6(\varPhi_{3}^{*}-\varPhi_{2}^{*})\rangle\), are nearly zero consistent with the findings in Ref. [3] and Fig. 2, respectively.

Fig. 5
figure 5

(Top row) The correlation (normalized to 1) between \(\varPhi_{2,1}^{*}\) and \(\varPhi_{3,1}^{*}\) in 40–50 % centrality interval for different weighting (r 2 or r n) and initial geometry models (Glauber or CGC). (Bottom row) The corresponding Fourier spectrum \(V_{1,2,3}^{i,j*}\). The constant term (i,j)=(0,0) (or \(V_{1,2,3}^{0,0}=1\)) is omitted for clarify. The typical statistical error on the Fourier coefficients is about 0.002 (Color figure online)

More results on other types of three plane correlations are summarized in Figs. 610. It is generally observed that the Fourier components are always bigger for r n-weighting than for r 2-weighting, and they are slightly bigger for the CGC geometry than for the Glauber geometry. Some of the correlators are quite large, e.g. \(\langle\cos(2\varPhi_{1}^{*}+2\varPhi_{2}^{*}-4\varPhi_{4}^{*})\rangle\), \(\langle\cos (2\varPhi_{1}^{*}-6\varPhi_{2}^{*}+4\varPhi_{4}^{*})\rangle\), \(\langle\cos(2\varPhi _{1}^{*}+6\varPhi_{2}^{*}-8\varPhi_{4}^{*})\rangle\), \(\langle\cos(\varPhi_{1}^{*}+3\varPhi _{3}^{*}-4\varPhi_{4}^{*})\rangle\), \(\langle\cos(2\varPhi_{1}^{*}-6\varPhi_{3}^{*}+4\varPhi _{4}^{*})\rangle\),\(\langle\cos(\varPhi_{1}^{*}+4\varPhi_{2}^{*}-5\varPhi_{5}^{*})\rangle \), \(\langle\cos(3\varPhi_{1}^{*}+2\varPhi_{2}^{*}-5\varPhi_{5}^{*})\rangle\), \(\langle \cos(\varPhi_{1}^{*}-6\varPhi_{3}^{*}+5\varPhi_{5}^{*})\rangle\), \(\langle\cos(2\varPhi _{1}^{*}+3\varPhi_{3}^{*}-5\varPhi_{5}^{*})\rangle\), and \(\langle\cos(4\varPhi _{1}^{*}-9\varPhi_{3}^{*}+5\varPhi_{5}^{*})\rangle\). These correlators are shown as a function of centrality in Figs. 11, 12, 13, and 14. In general, they all increase from central to more peripheral collisions, however the rate of the change depends on the type of the correlator. The correlator \(\langle\cos(\varPhi_{1}^{*}+2\varPhi_{2}^{*}-3\varPhi_{3}^{*})\rangle\) has the largest values in most cases, except for r n weighting in central and mid-central collision, where the correlator \(\langle\cos (\varPhi_{1}^{*}+3\varPhi_{3}^{*}-4\varPhi_{4}^{*})\rangle\) has the largest values.

Fig. 6
figure 6

The correlation between \(\varPhi_{2,1}^{*}\) and \(\varPhi_{4,1}^{*}\) (top row) and the corresponding Fourier coefficients (bottom row) in 40–50 % centrality interval (Color figure online)

Fig. 7
figure 7

The correlation between \(\varPhi_{3,1}^{*}\) and \(\varPhi_{4,1}^{*}\) (top row) and the corresponding Fourier coefficients (bottom row) in 40–50 % centrality interval (Color figure online)

Fig. 8
figure 8

The correlation between \(\varPhi_{3,2}^{*}\) and \(\varPhi_{4,2}^{*}\) (top row) and the corresponding Fourier coefficients (bottom row) in 40–50 % centrality interval (Color figure online)

Fig. 9
figure 9

The correlation between \(\varPhi_{2,1}^{*}\) and \(\varPhi_{5,1}^{*}\) (top row) and the corresponding Fourier coefficients (bottom row) in 40–50 % centrality interval (Color figure online)

Fig. 10
figure 10

The correlation between \(\varPhi_{3,1}^{*}\) and \(\varPhi_{5,1}^{*}\) (top row) and the corresponding Fourier coefficients (bottom row) in 40–50 % centrality interval (Color figure online)

Fig. 11
figure 11

The centrality dependence of various significant three plane correlators for r 2 weighting with Glauber geometry (Color figure online)

Fig. 12
figure 12

The centrality dependence of various significant three plane correlators for r n weighting with Glauber geometry (Color figure online)

Fig. 13
figure 13

The centrality dependence of various significant three plane correlators for r 2 weighting with CGC geometry (Color figure online)

Fig. 14
figure 14

The centrality dependence of various significant three plane correlators for r n weighting with CGC geometry (Color figure online)

Interestingly, most of these correlators, when defined relative to the major axis, are positive (the negative values are indicated with red “x” in Figs. 510). However, some of these correlations are likely to be strongly distorted due to the mixing during the hydrodynamic evolution, especially for those involving \(\varPhi _{4}^{*}\) and \(\varPhi_{5}^{*}\). Nevertheless, it would be interesting to measure these quantities experimentally and compare with our predictions.

4 Discussion and conclusion

In summary, we discussed a method for measuring the correlations between the directions of the anisotropic flow of different orders. This method involves Fourier transforming the differential distribution of the correlations between different event planes into various Fourier components, where each component is corrected separately by an event plane resolution term. This method has the advantage of simultaneously analyzing many different correlators, especially those involving three or more different event planes, thus help identifying significant components.

The expected strength of various two- or three-plane correlators are estimated using a Monte Carlo Glauber model or CGC model. Strong positive correlations are seen between the major axes of two eccentricities in mid-central and peripheral collisions for (n,m)=(2,4), (2,6), (3,6), and those involving the first-order eccentricity. Similarly, several significant three-plane correlators have also been identified, revealing novel correlation patterns expected from the average geometry and/or initial state fluctuations. These strong correlations imply the need to measure several two- or three-plane correlators in order to describe the full distribution. A detailed comparison of the correlations calculated here with the data could shed light on the role of the initial geometry fluctuations and dynamic mixing during the hydrodynamic evolution leading to harmonic flow in the final state.

Our discussion so far has decoupled the magnitude of the flow v n from its phases Φ n . In principle, the correlation is ill-defined for events with very small v n . However these events are expected to have very broad raw correlation distribution and very poor resolution (i.e. small Res{jkΨ n }), thus their contributions to the numerator and the denominator of Eq. (6) are naturally suppressed. Nevertheless, it is possible that the strength of the event plane correlation may depend on the magnitude of the v n . This dependence can be studied by first divide the events in a given centrality bin into various classes according to e.g. their v 2 values, measure the raw correlation and resolution factors in each event class, and then use Eq. (6) to obtain the true correlation strength for each class. This may provide further insight on how the event plane correlation depends on the eccentricity of the initial geometry (e.g. ϵ 2 if events are classified according to v 2).

We acknowledge valuable discussions with Matthew Luzum and Jean-Yves Ollitrault for clarifying technical details in their multi-particle correlation framework. We thank Roy Lacey for a careful proofreading of the manuscript. This research is supported by NSF grant PHY-1019387.