A method for studying initial geometry fluctuations via event plane correlations in heavy ion collisions

A method is proposed to measure the relative azimuthal angle distributions involving two or more event planes of different order in heavy ion collisions using a Fourier analysis technique. The analysis procedure is demonstrated for correlations involving two and three event planes (Phi_n, Phi_m and Phi_h). The Fourier coefficients of these distributions are found to coincide with previously proposed correlators, such as cos(6Phi_2-6Phi_3) and cos(Phi_1+2Phi_2-3Phi_3) etc, hence the method provides a natural framework for studying these correlators at the same time. Using a Monte Carlo Glauber model to simulate Au+Au collisions with fluctuating initial geometry, we are able to identify several new two- or three-plane correlators that have sizable magnitudes and should be measured experimentally.


I. 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 [1][2][3]. 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]: ǫ n = r 2 cos nφ 2 + r 2 sin nφ 2 r 2 , 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 Φ n = atan2( r 2 sin nφ , r 2 cos nφ ) n + π n (2) and Φ * n = Φ n + π/n respectively. The minor axis direction Φ n is also known as the n th -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 n th -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 1 . Previous studies [3,[7][8][9][10][11] 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.

II. METHOD
The n th -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: 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 cos jk(Φ n − Φ m ) can be obtained by calculating the raw coefficients cos jk(Ψ n − Ψ m ) , followed by a simple correction for finite event plane resolution: 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 1 2,4 is simply the ratio Additional examples include cos 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: where we use the constraint c 1 + 2c 2 + 3c 3 = 0 and we adopt the short-hand notations: . 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: This series is expected to converge quickly for nonperipheral collisions. Therefore, only the terms for i, j ≤ 3 are required (see Fig. 5). The coefficients are: V i,±j 1,2,3 = cos (iΦ 2,1 ± jΦ 3,1 ) = cos iΦ 2,1 cos jΦ 3,1 ∓ sin iΦ 2,1 sin jΦ 3,1 (11) Under this notation, the two-plane correlator can be treated as special case: 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: (12) where we have assumed that Ψ n is distributed randomly around Φ n , such that sin jn (Ψ n − Φ n ) = 0. Again, the Ψ 1 , Ψ 2 and Ψ 3 should be calculated from subevents covering different η acceptances to avoid auto-correlations.
Other triple plane correlators can be similarly analyzed, the first few are 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 /(dΦ 3,2 dΦ 5,1 dΦ 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 j n,m or V i,±j n,m,h are related to the previously proposed multi-particle correlators from Ref. [10,11].
That approach effectively applies a |c n |-particle weight v 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Φ n + c m mΦ m ) and cos(c n nΦ n + c m mΦ m + c h hΦ h ), the weights are v 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 : where nc l − mc m = 0, φ 1 , ..., φ a+b are azimuthal angles of a + b particles, and subscript v {a} n v {b} m 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: 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 Eq. 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 ac-ceptance, 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].

III. EXPECTED BEHAVIOR FROM GLAUBER AND CGC MODELS
A. Correlation between two planes Similar to [7], we use a simple Glauber model [24] to estimate the level of the correlations between n th -and m th -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 nucleonnucleon 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 : The corresponding Fourier coefficients are denoted as V j * n,m = cos k(Φ * n − Φ * m ) , and are related to V j n,m as 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 mod 2π) = 0 or π. The latter case leads to a sign flip: V j * n,m = −V j n,m .  The top panels of Fig. 1 show the 4(Φ * 2 − Φ * 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 1 * n,m 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 anticorrelation: the distributions have their minima at 0, and the sign of V j 2,4 alternates between positive and negative: V j 2,4 = (−1) j V j * 2,4 . 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 1 * n,m , 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 1 * 2,6 in mid-central collisions). In particular, strong correlations are observed for 4(Φ * 2 − Φ * 4 ), 6(Φ * 3 − Φ * 6 ) and 6(Φ * 2 − Φ * 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].
To get a feeling on how many V j * n,m terms are needed to exhaust the information encoded in distribution k(Φ * n − Φ * m ), in Fig. 3 we show the centrality distribution of V j * n,m 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(Φ * 2 − Φ * 4 ) and 3(Φ * 1 − Φ * 3 ). The results presented in Figs. 1-3 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] 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 j * n,m and affects the relative magnitude of V j * n,m between (n, m) =(2,4) and (3,6) in central collisions.
Due to the phase shift between the two types of correlations given by Eq. 25, many positive V j * n,m values in Fig. 2 and 3 would imply the corresponding V j n,m 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 n th -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.
More results on other types of three plane correlations are summarized in Figs. 6-10. 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. cos(2Φ , and cos(4Φ * 1 − 9Φ * 3 + 5Φ * 5 ) . These correlators are shown as a function of centrality in Figs. 11-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 cos(Φ * 1 + 2Φ * 2 − 3Φ * 3 ) has the largest values in most cases, except for r n weighting in central and mid-central collision, where the correlator has the largest values. Interestingly, most of these correlators, when defined relative to the major axis, are positive (the negative values are indicated with red "x" in Figs. 5-10). However, some of these correlations are likely to be strongly distorted due to the mixing during the hydrodynamic evolution, especially for those involving Φ * 4 and Φ * 5 . Nevertheless, it would be interesting to measure these quantities experimentally and compare with our predictions.

IV. 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 ).