Azimuthal Anisotropy Distributions in High-Energy Collisions

Elliptic flow in ultrarelativistic heavy-ion collisions results from the hydrodynamic response to the spatial anisotropy of the initial density profile. A long-standing problem in the interpretation of flow data is that uncertainties in the initial anisotropy are mingled with uncertainties in the response. We argue that the non-Gaussianity of flow fluctuations in small systems can be used to disentangle the initial state from the response. We apply this method to recent measurements of anisotropic flow in Pb+Pb and p+Pb collisions at the LHC. The response coefficient is found to decrease mildly as the system becomes smaller. This mild decrease is consistent with a low value of the ratio of viscosity over entropy.


I. INTRODUCTION
The large magnitude of elliptic flow, v 2 , in relativistic heavy-ion collisions at RHIC [1] and LHC [2] has long been recognized as a signature of hydrodynamic behavior of the strongly-interacting quark-gluon plasma [3]. v 2 is understood as the hydrodynamic response to the initial anisotropy, ε 2 , of the initial density profile [4]. However, the magnitude of this anisotropy is poorly constrained theoretically [5,6]. This uncertainty hinders the extraction of the properties of the quark-gluon plasma from experimental data [7,8].
The statistical properties of anisotropic flow are now precisely known [9]. The ATLAS collaboration has analyzed the full probability distribution of v 2 , v 3 and v 4 in Pb+Pb collisions for several centrality windows [10]. In p+Pb collisions, information is less detailed, but the first moments of the distribution of v 2 have been measured [11]. Our goal is to make use of these measurements to separate the initial state from the response without assuming any particular model of the initial conditionsby only using a simple functional form which goes to zero at the geometric limits of ε n = 0 and 1.
In theory, one can describe the particles emitted from a collision with an underlying probability distribution [12]. Anisotropic flow, v n , is defined as the n th Fourier coefficient of the azimuthal probability distribution P(ϕ): where we have used a complex notation [13,14]. Note that the underlying probability distribution P(ϕ) and V n fluctuate event to event, but they are both theoretical quantities which cannot be measured on an event-byevent basis. The particles that are detected in an event represent a finite sample of P(ϕ), and the measurement of the probability distribution of v n involves a nontrivial unfolding of statistical fluctuations [10].

II. DISTRIBUTION OF εn
We assume that the fluctuations of v n for n = 2, 3 are due to fluctuations of the initial anisotropy ε n in the corresponding harmonic, defined by [15] E n = ε n e inΦn ≡ − r n e inφ ρ(r, φ)rdrdφ r n ρ(r, φ)rdrdφ .
where ρ(r, φ) is the energy density near midrapidity shortly after the collision, and (r, φ) are polar coordinates in the transverse plane, in a centered coordinate system. We assume for the moment that v n in a given event is determined by linear response to the initial anisotropy, v n = κ n ε n , where κ n is a response coefficient which does not fluctuate event to event. Event-by-event hydrodynamic calculations [16] show that this is a very good approximation for n = 2, 3. Within this approximation, it has already been shown that one can rule out particular models of the initial density using either a combined analysis [17,18] of elliptic flow and triangular flow [19] data, or the relative magnitude of elliptic flow fluctuations [20][21][22]. Our goal is to show that one can extract both κ n and the distribution of ε n from data. We hope to show that this is true even if we relax the linear assumption. We make use of the recent observation that the distribution of ε n is to a large extent universal [23,24] and can be characterized by two parameters.
Both the magnitude and direction of E n fluctuate event to event. The simplest parametrization of these fluctuations is a two-dimensional Gaussian probability distribution which, upon integration over azimuthal angle, yields the Bessel-Gaussian distribution [25]: where ε 0 is the mean anisotropy in the reaction plane, which vanishes by symmetry for odd n, and σ is the typical magnitude of eccentricity fluctuations around this mean anisotropy. In a previous publication [24], we have introduced an alternative parametrization, the Elliptic Power distribution: where α describes the fluctuations and is approximately proportional to the number of sources in an independentsource model [4]. When ε 0 1 and α 1, Eq. (4) reduces to Eq. (3) with σ ≈ 1/ √ 2α. Its support is the unit disk: it naturally takes into account the condition |ε n | ≤ 1 which follows from the definition, Eq. (2). For this reason, it is a better parametrization than the Bessel-Gaussian, in particular for large anisotropies. Eq. (4) has been shown to fit various initial-state models [24]. Note that ε 0 is not strictly equal to the mean reaction plane eccentricity for the Elliptic Power distribution, but the difference is small for Pb+Pb collisions [24].
When the anisotropy is solely due to fluctuations, ε 0 = 0, the Bessel-Gaussian reduces to a Gaussian distribution: and the Elliptic Power distribution reduces to the Power distribution [23]:

III. DISTRIBUTION OF vn
The probability distribution of anisotropic flow, P (v n ), is obtained from the distribution of the initial anisotropy Assuming v n = κ n ε n , this becomes: In this case the distribution is rescaled by the response coefficient κ n . Figure 1 displays the probability distribution of v 2 and v 3 in various centrality windows [10] together with fits using rescaled Bessel-Gaussian and Elliptic Power distributions for v 2 , and rescaled Gaussian and Power distributions for v 3 . Note that σ and α need not be the same for v 2 and v 3 [24]. Both parametrizations give very good fits to v 2 and v 3 data for the most central bins shown on the figure. 1 As the centrality percentile increases, however, the quality of the Bessel-Gaussian fit becomes increasingly worse, which is reflected by the large χ 2 of the fit, and also clearly seen in the tail of the distribution: it systematically overestimates the distribution for large anisotropies. On the other hand, the Elliptic Power fit is excellent for all centralities. In particular, it falls off more steeply for large v n , in close agreement with the data. Note that the Bessel-Gaussian distribution Eq. (3) is scale invariant: rescaling it by κ n amounts to multiplying both ε 0 and σ by κ n , so that the fit is degenerate: only the products κ n ε 0 and κ n σ can be determined. Therefore the Bessel-Gaussian fit to ATLAS data is in practice a 2parameter fit for v 2 , and a 1-parameter fit for v 3 . On the other hand, the Elliptic Power fit is not degenerate because of the non-Gaussian cut-off at ε n = 1, and returns both the response coefficient κ n and the parameters pertaining to the shape, namely α and ε 0 (for v 2 ). However, the fit parameters are still correlated in the sense that the combinations κ n / √ α and κ n ε 0 (for v 2 ) have much smaller errors than each individual parameter.

IV. EXPERIMENTAL ERRORS
Our ability to separate the response from the initial eccentricity thus lies in the difference between the Bessel-Gaussian and the Elliptic Power fits, that is, in the non-Gaussianity of flow fluctuations. Since the difference is small, errors must be carefully evaluated.
The ATLAS collaboration reports the statistical error, the systematic error on the mean v n , and the systematic error on the relative standard deviation σ vn / v n . The first systematic error is an error on the scale of the distribution, while the second is an error on its shape. The error on the scale directly translates into an error of the response coefficient κ n , of the same relative magnitude. Since our analysis uses the deviations from a Gaussian shape, the dominant source of error is -by far-the error on the shape. In order to estimate the corresponding error on our fit parameters, we distort the distribution of v n in such a way that the mean v n is unchanged, and σ vn / v n is increased or decreased by the experimental uncertainty. This is done in practice by shifting the values of the v n bins according to v n → v n + δ(v n ), where δ(v n ) is a small non-linear shift. We choose the ansatz where v max is the tail of the v n distribution, λ is chosen in such a way that v n is unchanged, and is chosen in such a way that σ vn / v n is increased or decreased by the systematic error. This non-linear transformation leaves the minimum and maximum values of v n invariant. Figure 2 displays the value of the response coefficients κ 2 and κ 3 as a function of the centrality percentile. They are smaller than unity, with κ 3 < κ 2 , in line with expectations from hydrodynamic calculations [16], and decrease mildly as a function of the centrality percentile, which is the general behavior expected from viscous corrections to local equilibrium [26,27]. We estimate that the lowp T cut of ATLAS at 0.5 GeV increases κ 2 by a factor 1.4 to 1.5. The systematic error for κ 3 is very large and therefore not shown: for most bins, the upper error bar goes all the way to infinity. Now, if one takes the limit κ 3 → ∞ while keeping the rms v 3 constant, α in Eq. The other parameters of the fit to v n distributions, namely, ε 0 and α, characterize the shape of the distribution. They are displayed in Fig. 3 as a function of the collision centrality. ε 0 increases smoothly with the centrality percentile: extrapolation to the most central collisions (where the fit does not converge) gives ε 0 = 0, as required by azimuthal symmetry. ε 2 for that particular model. Specifically, the dashed line at the edge of the band is the value returned by a 2-parameter Elliptic Power fit to the distribution of ε 2 . The full line at the other edge of the band is the value that the fit to the v 2 distribution would return if v 2 ∝ ε 2 , with ε 2 given by that model. If one assumes linear response, ATLAS data deviate from both models.

VI. CUMULANTS
For p+Pb collisions, the full distribution of v 2 has not been measured, but only its first cumulants [33,34] 35,36]. Assuming linear response to the initial eccentricity, each measured cumulant is proportional to the corresponding cumulant of the initial eccentricity [37], v 2 {k} = κ 2 ε 2 {k}, for k = 2, 4, 6... The eccentricity in p+Pb collisions is solely due to fluctuations [38,39], therefore Eqs. (5) and (6) apply. While cumulants of order 4 and higher vanish for the Gaussian distribution Eq. (5), the Power distribution Eq. (6) always gives ε n {4} > 0 [23]. We again use this non-Gaussianity to disentangle the initial state from the response: We extract α from the measured ratio v n {4}/v n {2} ε n {4}/ε n {2} = (1+ α 2 ) −1/4 [23]. The rms anisotropy is then obtained as ε n {2} = 1/ √ 1 + α [23]. One finally obtains for the Power distribution: The values of κ 2 extracted from CMS p+Pb data [11] using this equation are also displayed in Fig. 2. We multiply them by a factor 1.19 to correct for the different low-p T cut (0.3 GeV/c) assuming a linear dependence of v 2 on p T . We plot p+Pb data at the equivalent centralities, determined according to the number of charged tracks [11]. General arguments have been put forward which suggest that the hydrodynamic response should be identical for p+Pb and Pb+Pb at the same equivalent centrality [40]. Once rescaled, the p+Pb slope is in line with Pb+Pb results, albeit somewhat steeper. Note that the fit parameters can also be obtained from cumulants for v 3 in Pb+Pb collisions using Eq. (9). For v 2 , there is a third parameter ε 0 , therefore one needs a third cumulant v 2 {6}. α and ε 0 , which control the shape of the distribution and its non-Gaussian features, can be extracted from the ratios v 2 {6}/v 2 {4} and v 2 {4}/v 2 {2} using the Elliptic Power distribution (Eq. (A5) of Ref. [24]). Note that while the Bessel-Gaussian Eq. (3) gives ε n {4} = ε n {6} = ε 0 [41], the Elliptic Power distribution always gives ε n {6} < ε n {4}. We have checked that α and ε 0 thus extracted from cumulant ratios are essentially identical to those obtained by fitting the distribution of v 2 . This approach has the advantage that cumulants can be analyzed without any unfolding procedure [34] but v 2 {2} may suffer from nonflow effects.

VII. DEVIATIONS FROM LINEAR SCALING
We now discuss the effect of deviations from linear eccentricity scaling of anisotropic flow. Because of such deviations, the shape of the v n distribution is not exactly the same as that of the ε n distribution, as already noted in event-by-event hydrodynamic calculations [42]. There are two distinct types of non-linearities: v n can be a function of ε n which is not exactly linear, or v n can depend on properties of the initial state other than ε n . We study these effects in turn.
Adding a quadratic term would be equivalent to rotating the distribution 90 degrees or changing the sign of v n . Thus the first significant non-linear term is the cubic: Several hydrodynamic calculations show evidence that κ > 0 [17,43,44], but no quantitative analysis has been done yet. One typically expects κ to depend mildly on centrality. When fitting the experimental v 2 distributions with the added parameter of the cubic term, κ had large errors but was in the range from 0 to 0.15. Thus we fixed κ at 0.10 and plotted the κ 2 values also in Fig. 2.
The effect of the non-linear response is to reduce the linear response coefficient κ 2 , essentially by a constant factor. In the case where the distribution of ε 2 is the Power distribution (6) and in the limit α 1, the relative change of κ 2 is to leading order. 2 With the Elliptic-Power distribution, the relative effect is also ∼ −2κ as can be seen in Fig. 2. Note that this non-linear correction to the response is much larger than one would naively expect from Eq. (10): the relative magnitude of the cubic term is κ ε 2 2 , yet it produces an effect of order κ . The reason is that the non-linear response contributes to the non-Gaussianity of flow fluctuations.
We now discuss deviations from linearity due to the fact that v n is not entirely determined by ε n [15]. One can generally decompose the flow as V n = κ n E n + X n , where X n is uncorrelated with the initial eccentricity E n . There can be various contributions to X n from non-linear coupling between different harmonics [45] or radial modulations of the initial density [15]. In order to estimate their effect on the hydrodynamic response, we further assume that X n is a Gaussian noise. Then, the distribution of V n is a rescaled Elliptic-Power distribution, convoluted with a Gaussian: the deviation from linearity here results in a Gaussian smearing of the distribution.
A quantitative measure of the magnitude of X n is the Pearson correlation coefficient r n between the anisotropic flow and the initial anisotropy, defined as where angular brackets denote an average value over events in a centrality class. Our analysis assumes the maximum correlation, |r n | = 1. Event-by-event hydrodynamic calculations show that there are small deviations around eccentricity scaling [46]. Ideal hydrodynamics [47] gives |r 2 | ∼ 0.95 for elliptic flow. However, the correlation between v n and ε n has been shown to be significantly larger in viscous hydrodynamics [16], and a value |r 2 | = 0.99 seems reasonable, but there is to date no quantitative estimate of |r 2 | as defined in Eq. (12). The effect on the fit parameters can be obtained using the fact that the rms flow v n {2} = |V n | 2 is increased by the noise, while higher-order cumulants v n {4} and v n {6} (see below Sec. VI) are unchanged. We find that a decrease of |r 2 | by 1% results in an decrease of the extracted κ 2 by 6% to 9%, depending on the centrality, the effect being maximum in the 20-30% centrality range. The value of |r 2 | found in ideal hydrodynamic calculations [47] depends mildly on centrality and is closest to 1 also in the 20-30% centrality range. Therefore one can conjecture -this should eventually be confirmed by detailed calculations-that the effect of the noise X n is to reduce the extracted response essentially by a constant factor, independent of centrality.
Note that the cubic response in Eq. (10) does not contribute to r 2 to first order in κ , so that the two effects are in practice well separated.
The conclusion is that deviations from linear eccentricity scaling all make κ 2 smaller, by a factor which can be significant, but depends little on centrality. This is of crucial importance for the extraction of the viscosity over entropy ratio (see below). The decrease in κ 2 makes ε 0 larger and α smaller (see Fig. 3), thereby improving compatibility with existing initial-state models.

VIII. VISCOUS HYDRO
We now compare our result for κ 2 with hydrodynamic calculations. In ideal hydrodynamics, scale invariance implies that the response coefficient κ 2 is independent of the system size, i.e., independent of centrality. Deviations from thermal equilibrium generally result in a reduction of the flow which is stronger for peripheral collisions [26,27]. In a hydrodynamic calculation [8], such deviations are due to the shear viscosity [7] and, to a lesser extent, to the freeze-out procedure at the end of the hydrodynamic expansion. Therefore the dependence of κ 2 on centrality in Fig. 2 can be used to estimate the shear viscosity over entropy ratio η/s of the quarkgluon plasma. We use the same hydrodynamic code as in Ref. [45] to estimate κ 2 . The resulting values are significantly smaller than the data in Fig. 2. Since we have shown that deviations from linear eccentricity scaling reduce κ 2 without altering its centrality dependence, we compensate for this effect, and the low p T cut of the AT-LAS data, by multiplying our hydrodynamic result by a constant, while tuning the viscosity so as to match the centrality dependence of κ 2 . The κ 2 lines in Fig. 2 are obtained with η/s = 0.19. This extracted value of η/s is consistent with that reported in the literature [8], using specific models of the initial state.

IX. SUMMARY
We have shown that a rescaled Elliptic Power distribution fits the measured distributions of elliptic and triangular flows in Pb+Pb collisions at the LHC. These distributions become increasingly non-Gaussian as the anisotropy increases. We have used this non-Gaussianity to disentangle for the first time the initial anisotropy from the response without assuming any particular model of initial conditions -just using a simple functional form which meets the geometrical constrants of eccentricity. This is another aspect of the analogy between heavyion physics and cosmology [48,49], where initial quantum fluctuations give rise to correlations, and the non-Gaussian statistics of these correlations can be used to unravel the properties of the initial state [50][51][52]. The non-Gaussianity is stronger for smaller systems, which is an incentive to analyze flow in smaller collision systems.
We have found that the hydrodynamic response to ellipticity has the expected overall magnitude and centrality dependence: it decreases mildly with centrality percentage. A somewhat similar slope is found for p+Pb collisions. This decrease can be attributed to the viscous suppression of v 2 . Comparison with hydrodynamic calculations supports a low value of the viscosity over entropy ratio, η/s ∼ 0.19.
The present study can be improved by constraining the cubic response coefficient κ in Eq. (10) as well as the Pearson coefficient due to other non-linear terms in Eq. (12). This could be done in future hydrodynamic calculations. Taking into account these nonlinear terms will decrease the magnitude of the response and therefore improve the agreement with hydrodynamic calculations. However, we have argued that this decrease is essentially a constant factor, independent of centrality, so that our estimate of η/s is likely to be robust. Our study is a first step toward the extraction of the viscosity over entropy ratio of the quark-gluon plasma from experimental data, without any prior knowledge of the initial state.