Event-by-Event flow

Event-shape fluctuations and flow correlations in ultra-relativistic heavy-ion collisions


I review recent measurements of a large set of flow observables associated with event-shape 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 non-linear hydrodynamic responses is derived, which qualitatively describe most experimental results. I discuss new types of fluctuation measurements that can further our understanding of the event-shape fluctuations and collective expansion dynamics.



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 space-time evolution of the matter from the measured azimuthal anisotropy. The success of these models in describing the anisotropy of particle production in heavy-ion 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 initially-produced fireball was treated as a smooth and boost invariant distribution of quarks and gluons, given by the overlap of two Woods-saxon functions that describe the distribution of nucleons in the two colliding nuclei. One important insight emerged around 2010 is the dominant role of event-by-event (EbyE) fluctuation [2] of nuclear wave-functions: 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 event-shape 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 Event-by-event 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:


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 two-dimensional vector or in a complex form:


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 event-by-event. Instead, they are estimated using the azimuthal distribution of particles in the event:


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 event-plane method [1] or the scalar-product 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]:


where denotes an average over the transverse position of all participating nucleons, and the and angle (also known as participant-plane, 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 non-trivial 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 :


with each variable being a function of , etc [34]. Among these, the joint probability distribution of the EP angles


can be reduced to the following EP correlators [21, 22, 23]:


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 event-plane 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]:


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 non-flow, 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 multi-particle correlation method.

The -particle correlation is often combined with correlations involving less number of particles to construct the corresponding -particle cumulant, which removes non-flow correlations of order less than , i.e:


For example, , and . However, since , many combinations vanish in multi-particle correlation analysis (see some concrete examples are given below).

Each of the three types of reduced pdfs in Tab. 1 has its own multi-particle correlations from Eq. LABEL:eq:223. The corresponding cumulant involves particular combination with several lower order moments. For example, ignoring non-flow, the first few moments and cumulants of accessible to the correlation analysis are :


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 lowest-order moment and cumulant for , which are accessible to the correlation analysis, involve four-particle correlation of the following form:

Clearly, this four-particle cumulant, first proposed in Ref. [35], reduces to zero if flow magnitudes and are un-correlated. Similar relations can be easily derived for correlation between three or more flow magnitudes.

The correlation between event-plane angles, , can be accessed via multi-particle correlation derived from Eq. LABEL:eq:223 ( denotes the sign of ):

which involves number of particles with . The corresponding cumulant has identical expression. The right-hand 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 lowest-order moment and cumulant for accessible to the correlation analysis can be obtained from a five-particle correlation:


This cumulant reduces to zero if flow amplitudes and EP angles are un-correlated. Note that if some of the indices are the same, the cumulants takes somewhat different form. For example the lowest-order non-zero cumulants for is:

A factor of three in the second term in the right-hand 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 event-shape 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 mid-rapidity 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 event-shape selection techniques are also sensitive to any differential correlation between and for fixed centrality, which would otherwise be washed-out when averaging over different initial configurations. One example is the strong anti-correlation 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 anti-correlation between and  [24].

pdfs cumulants event-shape method
, ,2,…
Obtained recursively as above
Table 1: Event-by-event flow observables in terms of probability density distributions (left column) and lowest-order cumulants accessible to the correlation analyses (middle column). Most of them can also be accessed via event-shape selection methods (right column). Note that if some of the indices are the same, the cumulants take somewhat different form, e.g. Eq. LABEL:eq:pmix2.

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]:


In contrast, the higher-order harmonics and arise from both the initial geometry and non-linear mixing of lower-order harmonics [32, 42, 15]. The relative contributions of linear and non-linear effects to these higher-order harmonics can be separated cleanly using the event-shape 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 event-averaged 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 event-plane method {EP} [1], two-particle correlation method {2PC} [49], scalar-product method {SP} [30], multi-particle cumulant method {2}, {4}… [50, 51], and lee-yang-zero method {LYZ} [52]. These methods are all based on multi-particle correlation concept, but they are constructed to have very different sensitivity to flow fluctuations and non-flow effects. In a nut-shell, the higher-order cumulant ({2} for ) and lee-yang-zero methods suppress both non-flow effects and flow fluctuations. The {2PC}, {SP} and {2} are closely related to the RMS value of :


and they are sensitive to not only flow fluctuations but also non-flow 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 non-trivial 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 MC-KLN 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 IP-Glasma model [55] takes into account gluons field fluctuations inside nucleons and the associated gluon saturation effects. As shown in Fig. 1, the IP-Glasma 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 higher-order harmonics.

Figure 1: Calulation of for –5 from viscous hydro-calculation based on IP-Glamas initial condition [12], compared with data rom RHIC (left column) with constant (top-left) and temperature dependent (bottom-left) and from LHC (right column) with initial flow (top-right) and temperature dependent (bottom-right). The is the mean value, while are the RMS values.

One challenge for the theory is the spectrum in ultra-central collisions, events in 0-1% 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] 2. Hence these collisions are expected to provide better constraints on the mechanism of the density fluctuations. The precision data from CMS and ATLAS show several features that are not described by models: 1) is comparable or larger than as shown in Fig. 2. This is challenging since naively one expect , but the should suffer larger viscous correction. 2) The has very different shape comparing to other harmonics (Fig. 3). It peaks at much lower , around 1.5 GeV compare to 3-4 GeV for other harmonics. 3) The is also observed to break the factorization relation at 20-30% level, while such breaking is much less in other centrality and for higher-order harmonics in all centrality [6, 56]. The current hydro calculations describe the factorization data for but over-predict those for the  [57].

Figure 2: The -integrated vs. in ultra-central Pb+Pb collisions from CMS for 0-0.2% centrality (left) and ATLAS for 0-1% centrality (right), compared with hydrodynamic model calculations with different initial conditions. Figures are taken from Ref. [56] and Ref. [10], respectively.
Figure 3: The for –6 in ultra-central Pb+Pb collisions from CMS for 0-0.2% centrality [56] and from 3+1D ideal hydrodynamic model calculations with initial conditions from AMPT model [58, 59].

In a recent work, G. Denicol et.al. suggest that the nucleon-nucleon 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 event-plane 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].

The data is precise enough to also improve the present modeling of other aspects of space-time evolution of the collisions, such as the early-time dynamics and initial flow, the temperature dependence of (see also Fig. 3), and non-linear hydrodynamic response for or  [12, 61, 62, 63].

3.2 Flow distribution

Due to large fluctuations in the initial density profile, varies strongly event to event. Calculations based on a Monte-Carlo (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 mid-central collisions where the is large and flow response is approximately linear, the fluctuations of eccentricity and flow coefficients can be approximated by a 2-D Gaussian [64]:


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 one-dimensional (1D) probability density of in the form of the Bessel–Gaussian (B-G) function [65, 64]:


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]:


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 multi-particle cumulants, {2}, {4} and so on. The first four of these cumulants can be expressed as [64]:


The is calculated as the cosine average of azimuthal angle of all combination of particles, and in the absence of non-flow, it is equivalent to the corresponding moment of the distribution:


where denote the average over all combinations in a event then over all events. For B-G distribution, these cumulants have a particularly simple form [64]:


Hence the observation that , see Figure 4 for , has been used as evidence that flow fluctuation is Gaussian. In this case, assuming that the non-flow contribution is small, the {2} and {4} can be used to extract (the component associated with average geometry) and (component associated with fluctuation).

Figure 4: ALICE measurements of (left) centrality dependence of , and estimated with multi-particle cumulants, and (right) transverse momentum dependence of and estimated with four-particle cumulants [66].

One limitation with the cumulant framework is that it is not possible to describe fluctuation using a small set of cumulants such that higher-order cumulants systematically vanishes, as in a taylor expansion. The higher-order 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 B-G 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 non-Gaussian, the higher-order cumulants for are very close to each other. Hence the similarity between , and can not be used to conclude the flow fluctuation is Gaussian 3

Figure 5: The distribution of Bessel-Gaussian function requiring . This distribution is divided into two halves with equal integral, labeled by A and B, respectively. The table at the right panel summarizes the values of cumulants for the full distribution and each half in units of , calculated via Eqs. 23.

For collisions with small (p+A or peripheral A+A collisions), the fluctuations of and deviate strongly from Gaussian and maybe described by a power-law function [67]:


where is proportional to . In this case, the collisions have no average geometry, nevertheless the for are generally non-zero and approximately equal to each other. The small but non-zero {4}, {6} and {4} (shown in Fig. 4) and {4} [66, 68] may be related to the non-Gaussianity of the distribution.

Figure 6: The event-by-event distributions of , and  [8]. Curves are fits to Bessel-Gaussian function Eq. 21 but with (or equivalently Eq. 22).

ATLAS employed a data-driven 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 non-flow. These effects are estimated by the response function , obtained from the difference of the calculated separately from two sub-events at forward and backward pseudorapidities. These two quantities are related to the observed flow vector as:


Since the statistical smearing and most non-flow 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:


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 B-G functions. The distributions suggest a small but non-zero {4} as indicated from the deviation of the data from a pure gaussian function Eq. 22. The distributions show significant deviations from B-G function in mid-central and peripheral collisions. One example is shown in the top-left 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: (Top-left) The distribution and associated fit by Bessel-Gaussian function. (Top-right) Centrality dependence of {8} obtained from multi-particle correlation method compared with that directly calculated from the . (Bottom) The ratios of the {6}/{4} and {8}/{4} from the two methods. Plots taken from Ref. [8, 68].

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 non-Gaussian 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 “elliptic-power” function [70]. Since , the distribution is also expected to fall faster than B–G function at large value (see top-left 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 non-linear behavior of is also observed in hydro calculations based on the IP-Glasma initial condition as shown in the right panel of Fig. 8.

Figure 8: (left) The data fit by Bessel-Gaussian and elliptic power function [71]. (Middle) vs. from viscous hydro calculation with EbyE fluctuating initial condition [63]. (Right) The scaled distribution compared to viscous hydro calculation with EbyE IP-Glasma initial conditions [12, 18].

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 event-plane angle as a function of and : , observed in hydrodynamic simulations [34, 57]. The correlation is often quantified by a linear correlation coefficient:


or -1 indicates perfect correlation or anti-correlation. 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 higher-order 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.

Figure 9: The pion vs. for , 3 and 4 in 0-5% central Pb+Pb collisions from EbyE viscous hydro calculations [13].

3.3 Event-plane correlations

The correlation between different harmonic planes can in general be captured by multi-particle 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 event-plane method and the scalar-product method used by ATLAS. Ignoring the statistical smearing and non-flow for the sake of simplicity, the correlators from these two methods reduce to:


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 Cauchy-Schwarz 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:


From this, we obtain the following approximate empirical relation:


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