Eventshape fluctuations and flow correlations in ultrarelativistic heavyion collisions
Abstract
I review recent measurements of a large set of flow observables associated with eventshape fluctuations and collective expansion in heavy ion collisions. First, these flow observables are classified and experiment methods are introduced. The experimental results for each type of observables are then presented and compared to theoretical calculations. A coherent picture of initial condition and collective flow based on linear and nonlinear hydrodynamic responses is derived, which qualitatively describe most experimental results. I discuss new types of fluctuation measurements that can further our understanding of the eventshape fluctuations and collective expansion dynamics.
pt
1 Introduction
Relativistic heavy ion collisions at the RHIC and the LHC create a hot and dense nuclear matter that is composed of strongly interacting quarks and gluons. This initially produced matter has an asymmetric shape in the transverse plane. Driven by the large pressure gradients arising from the strong interactions, the matter expands collectively and transfers the asymmetry in the initial geometry into azimuthal anisotropy of produced particles in momentum space [1, 2]. Hydrodynamic models are used to understand the spacetime evolution of the matter from the measured azimuthal anisotropy. The success of these models in describing the anisotropy of particle production in heavyion collisions at RHIC and the LHC [3, 4, 5, 6, 7, 8, 9] places important constraints on the transport properties, such as ratio of shear viscosity to entropy density , and initial conditions of the produced matter [10, 11, 12, 13, 14, 15].
For many years, the initiallyproduced fireball was treated as a smooth and boost invariant distribution of quarks and gluons, given by the overlap of two Woodssaxon functions that describe the distribution of nucleons in the two colliding nuclei. One important insight emerged around 2010 is the dominant role of eventbyevent (EbyE) fluctuation [2] of nuclear wavefunctions: the transverse positions of nucleons in the overlap region, as well as the parton density profiles inside those nucleons can fluctuate from collision to collision (earlier work on fluctuations can be found in Ref. [16]). As a result, each collision produces a different and lumpy fireball, each has its own shape and each features its own hydrodynamic expansion. Indeed, the azimuthal distribution of produced particles, when expanded into a Fourier series, shows significant harmonics up to at least 6 order. The relative strength of these harmonics, their centrality, transverse momentum () and particle mass dependence are well described by hydrodynamic calculations that were also performed on an EbyE basis [17, 18, 19, 20].
The rich patterns in the initial state fluctuations and the resulting hydrodynamic expansion also imply a large space of information associated with many new flow observables [21, 22, 23, 24, 25]. Initial measurements have been carried out on some of these observables, noticeably the probability distribution of individual harmonics [8] and amplitude or phase correlations between different harmonics [26, 9]. Studies of these observables are further augmented using recently proposed eventshape techniques [27, 24, 28]. These measurements already provided unprecedented insights on the nature of the initial density fluctuations and dynamics of the collective evolution. Future detailed mapping of these flow observables requires concerted efforts. The goal of this article is to provide a concise introduction of these new flow observables, review the status of current measurements, and discuss open issues and future physics opportunities in exploring these observables. More information on the current status of the harmonic flow measurements can be found in a contribution by R. Snellings in the same issue [29].
2 Eventbyevent flow observables
2.1 Eccentricities and azimuthal flow harmonics
When describing the transverse expansion dynamics, it is convenient to parameterize the probability distribution of the particle production in azimuthal angle in each event by a Fourier expansion:
(1) 
where and represent the magnitude and the phase (referred to as the event plane or EP) of the order harmonic flow, which are often represented as a twodimensional vector or in a complex form:
(2) 
The first few harmonics are referred to as dipolar, elliptic, triangular, quadrangular flow etc.
Since the number of particles in each event is finite, the values of and can not be obtained eventbyevent. Instead, they are estimated using the azimuthal distribution of particles in the event:
(3) 
where the average is over all produced particles in the event. The estimated flow vector smears around . This smearing can be removed statistically via an unfolding method [8] or corrected for on average via the eventplane method [1] or the scalarproduct method [30].
The development of harmonic flow is driven by the asymmetries in the pressure gradients of the matter, which in turn is controlled by the detailed shape configuration of the initial density profile. The shape configuration of each event is often characterized by a set of eccentricity vector , calculated from the transverse positions of the participating nucleons relative to their center of mass [2, 11]:
(4) 
where denotes an average over the transverse position of all participating nucleons, and the and angle (also known as participantplane, PP) represent the magnitude and orientation of the eccentricity vector, respectively. Hydrodynamic calculations shows that the first few flow harmonics are directly related to the eccentricities of the corresponding order, e.g. for [31, 32]. By default the radial weights are chosen as for , and for [11]. But sometimes they are also calculated using other weights, for example [2]. These alternative definitions of eccentricity capture the degree of freedom along the radial direction, and reflect important information of the system, such as the system size and density gradients [11, 32, 33]. Events with the same values in the default definition Eq. 4 may have different values when calculated with alternative weights.
2.2 Flow fluctuations and observables
In heavy ion collisions, the number of participating nucleons is finite and their positions fluctuate randomly in the transverse plane, leading to strong EbyE fluctuations of and . These fluctuations also result in nontrivial correlations between eccentricities and PP angles of different order characterized by [24, 25]. Consequently, the matter created in each collision follows a different collective expansion, with its own set of flow harmonics. Experimental observables describing harmonic flow can be generally given by the joint probability distribution (pdf) of all and :
(5) 
with each variable being a function of , etc [34]. Among these, the joint probability distribution of the EP angles
(6) 
can be reduced to the following EP correlators [21, 22, 23]:
(7) 
Due to fold symmetry of , these correlators should be invariant under a phase shift . It should also be invariant under a global rotation by any angle. The first condition requires that the EP angle of the order harmonic appears as integer multiple of , while the second condition requires the sum of the coefficients to vanish. Together they lead to the constraint in Eq. 7.
Heavy ion experiments at RHIC and LHC have recorded billions of Pb+Pb or Au+Au collisions, which are distributed according to the underlying pdf given by Eq. 5. However, it is more practical to measure the projections of the full probability distribution on a finite number of variables. These projected distributions can be generally classified into three types: 1) those involving only the flow amplitudes, 2) those involving only the flow phases or eventplane angles, and 3) those involving both amplitudes and phases. These are listed in the left column of Tab. 1.
It is not always straightforward to directly measure the flow probability distributions, instead, the information of Eq. 5 can also be obtained in terms of its moments measured by a particle azimuthal correlation of the following general form [21]:
(8)  
where . The double average in the left hand side is carried out over all particles in one event then over all events, while the single averages are carried out over the events. All sine terms drops out after the averages. The event average also removes any statistical smearing effects, leading to the second line of the equation. Note that the angle refers to the azimuthal angle of particle, while and refers to the order observed and true event plane, respectively. This equation is also conveniently expressed as:
In the absence of nonflow, the particle correlation reduces to one particular order moment (Eq. LABEL:eq:223) of the underlying flow probability distribution.
According to the traditional definition of moments or cumulants, the order moment of the underlying pdf Eq. 5 should has the form of with or . However the azimuthal correlation analysis discussed above uses a different definition , which in general does not capture the full information of the pdf. For example, since , all the odd moments of , such as or , can not be accessed directly via the multiparticle correlation method.
The particle correlation is often combined with correlations involving less number of particles to construct the corresponding particle cumulant, which removes nonflow correlations of order less than , i.e:
(10) 
For example, , and . However, since , many combinations vanish in multiparticle correlation analysis (see some concrete examples are given below).
Each of the three types of reduced pdfs in Tab. 1 has its own multiparticle correlations from Eq. LABEL:eq:223. The corresponding cumulant involves particular combination with several lower order moments. For example, ignoring nonflow, the first few moments and cumulants of accessible to the correlation analysis are :
(11) 
(12) 
with subscript “” denoting the cumulant form and with the assumption that is the same for the two particles. The single particle flow coefficients are then calculated from these cumulants, e.g:
Similarly, the lowestorder moment and cumulant for , which are accessible to the correlation analysis, involve fourparticle correlation of the following form:
Clearly, this fourparticle cumulant, first proposed in Ref. [35], reduces to zero if flow magnitudes and are uncorrelated. Similar relations can be easily derived for correlation between three or more flow magnitudes.
The correlation between eventplane angles, , can be accessed via multiparticle correlation derived from Eq. LABEL:eq:223 ( denotes the sign of ):
which involves number of particles with . The corresponding cumulant has identical expression. The righthand side of this expression, referred to as the scalar product (SP) [30, 36], is similar to that defined in Eq. 7, except for the weight. This definition was argued to be preferable over Eq. 7 as the results are independent of resolution of the event planes in the experiments [36, 37].
The last category of pdf in Tab. 1 is the mixed correlations involving both magnitudes and the phases of harmonic flow, e. g. . The cumulant form can be easily obtained from Eq. LABEL:eq:223. For example, the lowestorder moment and cumulant for accessible to the correlation analysis can be obtained from a fiveparticle correlation:
(15)  
This cumulant reduces to zero if flow amplitudes and EP angles are uncorrelated. Note that if some of the indices are the same, the cumulants takes somewhat different form. For example the lowestorder nonzero cumulants for is:
A factor of three in the second term in the righthand of the equation arises because one can swap or with .
The flow observables listed in Tab. 1 can also be accessed via the recently proposed eventshape selection method [27, 24, 28, 26]. In this method, events in a narrow centrality interval are further classified according to the observed signal ( and 3) in a forward rapidity range. This classification selects events with similar multiplicity but very different ellipticity or triangularity. The values of are then measured at midrapidity using the standard flow techniques. Since the distributions are very broad, the correlation of with and/or can be explored over a wide range. The eventshape selection techniques are also sensitive to any differential correlation between and for fixed centrality, which would otherwise be washedout when averaging over different initial configurations. One example is the strong anticorrelation between and predicted by the MC Glauber model [38, 24]. A recent transport model calculation shows that this correlation survives the collective expansion and appears as a similar anticorrelation between and [24].
pdfs  cumulants  eventshape method  
, ,2,…  
,  
…  


yes  
…  
…  Obtained recursively as above  


yes  


yes 
Last few years witnessed impressive progresses in studying these flow observables, both experimentally and theoretically. They have greatly improved our understanding of the fluctuations in the initial density profile and hydrodynamics response in the final state [39, 40, 41]. We now know that the elliptic flow and triangular flow are the dominant harmonics, and they are driven mainly by the linear response to the ellipticity and triangularity of the initially produced fireball [31, 32]:
(18) 
In contrast, the higherorder harmonics and arise from both the initial geometry and nonlinear mixing of lowerorder harmonics [32, 42, 15]. The relative contributions of linear and nonlinear effects to these higherorder harmonics can be separated cleanly using the eventshape selection techniques. The details are discussed in Secs. 3.3 and 3.4.
The presence of large EbyE fluctuations of and also has consequence on the anisotropy of jet production at high . Since the energy loss of a jet depends on the length and local energy density along its path traversing the medium, jet production rate is expected to be sensitive to fluctuations of the event shape in the initial state. Recently model calculations predicts sizable – for jet production in A+A collisions at high [43, 44, 45].
3 Results
3.1 Differential measurements of single flow harmonics
Until recently, most flow studies were aimed at eventaveraged coefficients, measured differentially as a function of , , centrality and particle species [3, 46, 47, 5, 6, 48]. Many experimental methods have been developed in these measurements, including the eventplane method {EP} [1], twoparticle correlation method {2PC} [49], scalarproduct method {SP} [30], multiparticle cumulant method {2}, {4}… [50, 51], and leeyangzero method {LYZ} [52]. These methods are all based on multiparticle correlation concept, but they are constructed to have very different sensitivity to flow fluctuations and nonflow effects. In a nutshell, the higherorder cumulant ({2} for ) and leeyangzero methods suppress both nonflow effects and flow fluctuations. The {2PC}, {SP} and {2} are closely related to the RMS value of :
(19) 
and they are sensitive to not only flow fluctuations but also nonflow effects, however the latter can be suppressed by requiring a rapidity gap between pair of particles. The most popular method, {EP}, is known to introduce nontrivial biases in the presence of flow fluctuations [53], and hence should be used with caution. For all practical purpose, it can be replaced by {SP}. More detailed discussion and comparison of these methods can be found in Ref. [54, 41].
Last few years (since 2010) also witnessed rapid development of hydrodynamic modeling of heavy ion collisions [12]. Confronted with large amount of data, theorists are able to fine tune their models in terms of both the initial conditions and the hydrodynamic response. For example, when the first measurement was obtained by the PHENIX and ALICE collaboration [3, 5], it became clear that the commonly used MCKLN initial condition was unable to simultaneously describe the and data. As more detailed differential data became available, most early initial geometry models (prior to 2010) have been ruled out. The recently developed IPGlasma model [55] takes into account gluons field fluctuations inside nucleons and the associated gluon saturation effects. As shown in Fig. 1, the IPGlasma initial condition combined with viscous hydrodynamic evolution, describe a large set of the measured spectrum. However, this description is not perfect everywhere, in particular for dipolar flow at LHC (see [12], but not shown here) and higherorder harmonics.
One challenge for the theory is the spectrum in ultracentral collisions, events in 01% centrality or less [6, 56]. The original motivation is that the initial conditions in these collisions are predominantly generated by fluctuations such that the magnitudes of first several are comparable, and that the hydrodynamic response is expected to be linear for all harmonics [31]
In a recent work, G. Denicol et.al. suggest that the nucleonnucleon correlation and effects of bulk viscosity can reduce the values relative to , thus partially but not completely resolving the hierarchy problem in Fig. 3 [60]. It is possible that the longitudinal fluctuations and eventplane decorrelation effects might be responsible for the smallness of and the strange shape, as indicated by a recent 3+1D hydro calculation that include longitudinal dynamics [58, 59].
3.2 Flow distribution
Due to large fluctuations in the initial density profile, varies strongly event to event. Calculations based on a MonteCarlo (MC) Glauber model show that, even for events in a very narrow centrality interval, can fluctuate from zero to several times its mean value, leading to a very broad probability density distribution [8]. In central and midcentral collisions where the is large and flow response is approximately linear, the fluctuations of eccentricity and flow coefficients can be approximated by a 2D Gaussian [64]:
(20) 
where the and represent the eccentricity and flow vector associated with average geometry in the reaction plane, and or reflects the width of the fluctuations. Integration of this function over the azimuthal angle gives the onedimensional (1D) probability density of in the form of the Bessel–Gaussian (BG) function [65, 64]:
(21) 
where is the modified Bessel function of the first kind. If the which is suitable for in central collision or for , the can be absorbed into the width parameter [8]:
(22) 
Thus when fluctuation is large, the value of is constrained mainly by the tail of the distribution.
Traditionally, the information of has been inferred from the multiparticle cumulants, {2}, {4} and so on. The first four of these cumulants can be expressed as [64]:
(23) 
The is calculated as the cosine average of azimuthal angle of all combination of particles, and in the absence of nonflow, it is equivalent to the corresponding moment of the distribution:
(24) 
where denote the average over all combinations in a event then over all events. For BG distribution, these cumulants have a particularly simple form [64]:
(25) 
Hence the observation that , see Figure 4 for , has been used as evidence that flow fluctuation is Gaussian. In this case, assuming that the nonflow contribution is small, the {2} and {4} can be used to extract (the component associated with average geometry) and (component associated with fluctuation).
One limitation with the cumulant framework is that it is not possible to describe fluctuation using a small set of cumulants such that higherorder cumulants systematically vanishes, as in a taylor expansion. The higherorder cumulants are obtained after cancellation between several large numbers, so it can have sizable systematic uncertainties. Furthermore, if is large, for may not be very sensitive to significant deviations from Bessel–Gaussian distribution. To illustrate this, one example BG distribution shown in Fig. 5 is divided into two equal halves, and the are calculated analytically for each half using Eqs. 23. Despite the fact that the truncated distribution in each half is nonGaussian, the higherorder cumulants for are very close to each other. Hence the similarity between , and can not be used to conclude the flow fluctuation is Gaussian
For collisions with small (p+A or peripheral A+A collisions), the fluctuations of and deviate strongly from Gaussian and maybe described by a powerlaw function [67]:
(26) 
where is proportional to . In this case, the collisions have no average geometry, nevertheless the for are generally nonzero and approximately equal to each other. The small but nonzero {4}, {6} and {4} (shown in Fig. 4) and {4} [66, 68] may be related to the nonGaussianity of the distribution.
ATLAS employed a datadriven unfolding method to obtain directly the distribution. In this method, an observed flow coefficient is calculated for each event. The distribution contains the effects of smearing due to finite number of particles and nonflow. These effects are estimated by the response function , obtained from the difference of the calculated separately from two subevents at forward and backward pseudorapidities. These two quantities are related to the observed flow vector as:
(27) 
Since the statistical smearing and most nonflow effects are not correlated between the two subevents, the distribution of observed flow vector is simply the convolution of true flow vector and the response function:
(28) 
Consequently, the truth flow distribution can be obtained by a standard unfolding procedure such a Bayesian unfolding. A detailed study based on HIJING and AMPT simulation [69] shows that appears as a random Gaussian smearing of the underlying flow distribution, justifying the unfolding procedure.
Figure 6 shows the , and from ATLAS in several centrality intervals for Pb+Pb collisions. The and distributions are well described by the BG functions. The distributions suggest a small but nonzero {4} as indicated from the deviation of the data from a pure gaussian function Eq. 22. The distributions show significant deviations from BG function in midcentral and peripheral collisions. One example is shown in the topleft panel of Fig. 7. This deviation leads to a 2% difference between {4} and {6}, but no difference between {6} and {8} (bottom panels of Fig. 7), implying that when is large, for responds slowly to deviation from B–G function. The small change of can be easily buried under the experimental systematic uncertainties.
Figure 7 also compares the obtained from particle correlations (left part of Eq. 24) with that calculated directly from the distribution (second part of Eq. 24). Excellent agreement within 1% is observed across the full centrality range. Hence cumulants calculated from distribution provides equivalent information as the traditional method, but it is more intuitive and transparent in the estimation of systematic uncertainties.
One reason for the nonGaussian behavior is that the eccentricity distribution is bounded from above by . As a result, itself deviates strongly from B–G distribution and can instead be parameterized by a “ellipticpower” function [70]. Since , the distribution is also expected to fall faster than B–G function at large value (see topleft panel of Fig. 7). But the problem is that the shape of the distributions from various initial geometry models can not be tuned to agree with distribution across the full centrality range. In peripheral collisions, for example, the always falls faster than in the tail of the distributions. There are two possible explanations: 1) Current modeling of is wrong, but response coefficient is constant across the full range [71]. In this case, the eccentricity distribution can be obtained as (see left panel of Fig. 8). 2) The modeling of is correct, but is a function of . A recent hydrodynamic model calculation [63] suggests that the increases slightly at large as shown in the middle panel of Fig. 8, hence the distribution is expect to decrease more slowly than the distribution, consistent with experimental observation [8]. Similar nonlinear behavior of is also observed in hydro calculations based on the IPGlasma initial condition as shown in the right panel of Fig. 8.
In most hydro calculation with fluctuating initial condition, a sizable spread has been observed in the correlation between EbyE and EbyE as shown in Fig. 8. This spread can be partially related to the fluctuations of eventplane angle as a function of and : , observed in hydrodynamic simulations [34, 57]. The correlation is often quantified by a linear correlation coefficient:
(29) 
or 1 indicates perfect correlation or anticorrelation. Given that can be defined with different radial weights, and that the correlation between these definitions also have significant spread, it is not surprising that the – correlation should not be perfect. However, it is not clear how much of this spread is due to the higherorder radial modes in the initial geometry or is generated dynamically in the final state (e.g. hydrodynamic noise or hadronic freezeout). This remains an open issue to be resolved in the future.
3.3 Eventplane correlations
The correlation between different harmonic planes can in general be captured by multiparticle correlation as given by Eq. LABEL:eq:223. A few such correlators have been measured by ALICE [66]. In practice, it is often desirable to cast these correlators into different forms, such as the eventplane method and the scalarproduct method used by ATLAS. Ignoring the statistical smearing and nonflow for the sake of simplicity, the correlators from these two methods reduce to:
(30)  
(31)  
Eq. 31 has been separated into two parts. The first part is the cosine average of the normalized distribution of , which is similar to Eq. 30 except for the weights. The second part, denoted as , is a factor that is always smaller than one due to CauchySchwarz inequality. If the angular correlation is independent of the magnitude of the flow, the first part of Eq. 31 should equal Eq. 30. In this case, the correlator from the SP method should always be smaller than that given by the EP method. Experimental results as shown in Fig. 10, on the other hand, suggest the opposite. This can happen only if the events with larger flow also have stronger angular correlations in order to compensate for the factor.
To verify this, we estimate the value of in the AMPT model [69], which was found to quantitatively describe the measured the EP correlators. The resulting factors, as shown in Fig. 11, are always below one as expected. Interestingly, these factors are also found to be approximately the same as the ratio of the correlators between the two methods:
(32) 
From this, we obtain the following approximate empirical relation:
(33) 
This relation suggests that the unweighted EP correlators are approximately a factor of weaker than the weighted case, but about half of this difference is cancelled out in the SP method (Eq. 31). In the future, it would be useful to calculate