Primordial non-Gaussianity in heavy-ion collisions

TL;DR

This paper demonstrates that cumulants of anisotropic flow in heavy-ion collisions reveal the non-Gaussian nature of the initial energy density fluctuations, with experimental data indicating a positively skewed distribution.

Contribution

It introduces a perturbative approach to connect flow cumulants with initial density fluctuations, highlighting the non-Gaussian features of the early-stage energy distribution.

Findings

Flow cumulants probe non-Gaussian initial conditions.

Negative $sc(3,2)$ indicates positive skewness in energy distribution.

Perturbative expansion links flow correlations to density fluctuations.

Abstract

We show that experimental data on cumulants of anisotropic flow in heavy-ion collisions probe the non-Gaussian statistics of the energy density field created right after the collision. We carry out a perturbative expansion of the initial anisotropies of the system in terms of its density fluctuations. We argue that the correlation between the magnitudes of elliptic flow and triangular flow, dubbed $sc(3,2)$, is generically of the same sign and order of magnitude as the kurtosis of triangular flow in a hydrodynamic picture. The experimental observation that these quantities are negative implies that the distribution of energy around a given point has positive skew.

Figures

Equations

$${\int }_{\mathbf{s}}\mathbf{s}⟨\rho (\mathbf{s})⟩=0,$$

$${\mathbf{s}}_0\equiv \frac{{\int }_{\mathbf{s}}\mathbf{s}\rho (\mathbf{s})}{{\int }_{\mathbf{s}}\rho (\mathbf{s})}=\frac{{\int }_{\mathbf{s}}\mathbf{s}\delta \rho (\mathbf{s})}{{\int }_{\mathbf{s}}\rho (\mathbf{s})}.$$

$${\epsilon }_n=\frac{{\int }_{\mathbf{s}}(\mathbf{s}-{\mathbf{s}}_0{)}^n\rho (\mathbf{s})}{{\int }_{\mathbf{s}}|\mathbf{s}-{\mathbf{s}}_0{|}^n\rho (\mathbf{s})},$$

$$\delta f$$

$$⟨f⟩$$

$$⟨E⟩={\int }_{\mathbf{s}}⟨\rho (\mathbf{s})⟩.$$

$${\epsilon }_2$$

$${\epsilon }_3$$

$$c_n\{2\}\equiv ⟨{\epsilon }_n{\epsilon }_n^{\ast }⟩,$$

$$c_n\{4\}$$

$$=$$

$$SC(3,2)$$

$$=$$

$$C_2({\mathbf{s}}_1,{\mathbf{s}}_2)\equiv ⟨\delta \rho ({\mathbf{s}}_1)\delta \rho ({\mathbf{s}}_2)⟩,$$

$$C_3({\mathbf{s}}_1,{\mathbf{s}}_2,{\mathbf{s}}_3)\equiv ⟨\delta \rho ({\mathbf{s}}_1)\delta \rho ({\mathbf{s}}_2)\delta \rho ({\mathbf{s}}_3)⟩.$$

$$C_4({\mathbf{s}}_1,{\mathbf{s}}_2,{\mathbf{s}}_3,{\mathbf{s}}_4)$$

$$-$$

$$⟨\delta f\delta g⟩=\frac{1}{⟨E{⟩}^2}{\int }_{{\mathbf{s}}_1,{\mathbf{s}}_2}f({\mathbf{s}}_1)g({\mathbf{s}}_2)C_2({\mathbf{s}}_1,{\mathbf{s}}_2).$$

$$⟨\delta f\delta g\delta h⟩=\frac{1}{⟨E{⟩}^3}{\int }_{{\mathbf{s}}_1,{\mathbf{s}}_2,{\mathbf{s}}_3}f({\mathbf{s}}_1)g({\mathbf{s}}_2)h({\mathbf{s}}_3)C_3({\mathbf{s}}_1,{\mathbf{s}}_2,{\mathbf{s}}_3).$$

$$⟨\delta f\delta g\delta h\delta k{⟩}_c$$

$$⟨\delta f\delta g\delta h\delta k{⟩}_c=\frac{1}{⟨E{⟩}^4}{\int }_{{\mathbf{s}}_1,{\mathbf{s}}_2,{\mathbf{s}}_3,{\mathbf{s}}_4}f({\mathbf{s}}_1)g({\mathbf{s}}_2)h({\mathbf{s}}_3)k({\mathbf{s}}_4)C_4({\mathbf{s}}_1,{\mathbf{s}}_2,{\mathbf{s}}_3,{\mathbf{s}}_4).$$

$$⟨\delta f\delta g\delta h\delta k\delta l⟩$$

$$⟨\delta f\delta g\delta h\delta k\delta l\delta m⟩$$

$$c_n\{2\}=\frac{⟨\delta {\mathbf{s}}^n\delta {\mathbf{s}}^{\ast n}⟩}{⟨|\mathbf{s}{|}^n{⟩}^2}.$$

$$c_n\{4{\}}_K$$

$$SC(3,2{)}_K$$

$$c_n\{4{\}}_S$$

$$SC(3,2{)}_S$$

$$c_n\{4{\}}_V$$

$$SC(3,2{)}_V$$

$$c_n\{4\}$$

$$SC(3,2)$$

$$\rho (\mathbf{s})=\sum _{j=1}^N\delta (\mathbf{s}-{\mathbf{s}}_j),$$

$$⟨\delta f\delta g⟩$$

Full text

I introduction

Particles emitted in heavy-ion collisions display specific correlation patterns, which depend little on rapidity, but strongly on azimuthal angle Chatrchyan:2012wg . These correlations are analyzed through the observation of elliptic flow, $v_{2}$ Aamodt:2010pa , and triangular flow, $v_{3}$ Alver:2010gr , which are the second and third Fourier coefficients, respectively, of the azimuthal distribution of particles Luzum:2011mm . In a hydrodynamic framework, the values of $v_{2}$ and $v_{3}$ originate from an almost linear response Niemi:2012aj to the elliptic and triangular anisotropies Teaney:2010vd , dubbed $\varepsilon_{2}$ and $\varepsilon_{3}$, respectively, that characterize the energy density profile of the system created in the collision process. Since the energy density profile fluctuates event to event, $\varepsilon_{2}$ and $\varepsilon_{3}$ also fluctuate. Their probability distribution, and consequently that of the final-state $v_{n}$ coefficients, carries information about the initial density fluctuations.

In this article, we use experimental data on $v_{2}$ and $v_{3}$ to gather evidence that primordial energy-density fluctuation in nucleus-nucleus collisions do not follow Gaussian statistics. Non-Gaussianities are generic in microscopic systems, where they appear as corrections to the central limit theorem Alver:2008zza ; Bhalerao:2011bp ; Yan:2013laa . Unlike the situation in the early Universe, where primordial non-Gaussianities are compatible with zero Ade:2013ydc , it is natural to expect that they are sizable in the early stages of heavy-ion collisions.

Since the relation between the initial anisotropies and the initial density is not linear, non-Gaussian fluctuations of $\varepsilon_{2}$ and $\varepsilon_{3}$ are nontrivially related to non-Gaussian fluctuations of the initial energy density. Our aim in this paper is to derive expressions relating them. We show that the correlation between the magnitudes of $v_{2}$ and $v_{3}$, as measured by the symmetric cumulant $SC(3,2)$ Bilandzic:2013kga ; ALICE:2016kpq ; STAR:2018fpo , and the kurtosis of the distribution of $v_{3}$ Abbasi:2017ajp as measured by $v_{3}\{4\}$ ALICE:2011ab , have a common origin. Their negative values can be ascribed to the positive skewness of the distribution of energy in a given transverse area. We also show that, paradoxically, the skewness of elliptic flow fluctuations Giacalone:2016eyu ; Sirunyan:2017fts ; Acharya:2018lmh is likely to be less sensitive to the skewness of density fluctuations.

In order to relate the fluctuations of the initial anistotropies to the fluctuations of the initial density profile, we carry out a systematic perturbative expansion in powers of the density fluctuations Blaizot:2014nia ; Gronqvist:2016hym , which is introduced in Sec. II. This is a general valid approach for large systems such as nucleus-nucleus collisions. Small systems (proton-nucleus and proton-proton collisions Albacete:2017ajt ; Sirunyan:2017uyl ; Aaboud:2018syf ; Acharya:2019vdf ) are not considered here. For simplicity, the calculations in the main body of this article are carried out for central collisions (i.e., at zero impact parameter). As we shall argue, the geometry of non-central collisions only matters when studying the fluctuations of $\varepsilon_{2}$, whose discussion is relegated to Appendix A. In Sec. III, we recall the definitions of cumulants of $\varepsilon_{n}$, and the definitions of cumulants of the initial energy density field. In Sec. IV, we evaluate the cumulants of order 4 of the distribution of $\varepsilon_{n}$, including the mixed cumulant $SC(3,2)$, to leading order in the perturbative expansion. In Sec. V, we test these perturbative results using Monte Carlo calculations. We define several measures of primordial non-Gaussianity which allow to compare models of initial conditions with experimental data, and we discuss the implications of current measurements.

II Perturbative expansion of $\varepsilon_{n}$

We first introduce a few notations. We denote the initial energy density profile of an event at mid-rapidity by $\rho({\bf s})$, where ${\bf s}$ labels a point in the transverse plane.111We shall not deal with the longitudinal dynamics of the system, which is not relevant for our discussion. We are interested in the event-by-event fluctuations of $\rho({\bf s})$ in a situation where macroscopic quantities (typically, the centers of colliding nuclei) are fixed. We denote the statistical average of $\rho({\bf s})$ by $\langle\rho({\bf s})\rangle$, and the the event-to-event fluctuation by $\delta\rho({\bf s})\equiv\rho({\bf s})-\langle\rho({\bf s})\rangle$. We choose a coordinate frame where the origin of the transverse plane lies at the center of the average energy density, so that

[TABLE]

where we use the short hand $\int_{{\bf s}}=\int{\rm d}x{\rm d}y$ for the integration over the transverse plane. Due to fluctuations, the center of the distribution $\rho({\bf s})$, which we denote by ${\bf s}_{0}$, fluctuates event to event. It is defined by:

[TABLE]

The initial eccentricity $\varepsilon_{2}$ and the initial triangularity $\varepsilon_{3}$ are then defined by

[TABLE]

where we have used the complex notation ${\bf s}=x+iy$, and the recentering correction ${\bf s}_{0}$ ensures that anisotropies are evaluated in a centered frame Teaney:2010vd ; Alver:2006wh . Under the transformation ${\bf s}\to{\bf s}e^{i\alpha}$, $\varepsilon_{2}$ and $\varepsilon_{3}$ are multiplied by $e^{2i\alpha}$ and $e^{3i\alpha}$, respectively, thus $\varepsilon_{n}$ is a measure of the deformation of $\rho({\bf s})$ in the $n^{\rm th}$ Fourier harmonic.

We now carry out a perturbative expansion of $\varepsilon_{n}$. We decompose $\rho=\langle\rho\rangle+\delta\rho$, and expand the right hand side of Eq. (3) in powers of $\delta\rho$,222Note that even though $\delta\rho({\bf s})$ may be of the same order, or even larger than the average density $\langle\rho({\bf s})\rangle$ locally, its contribution to $\varepsilon_{n}$ remains small for a large system Blaizot:2014nia . keeping in mind that ${\bf s}_{0}$ is itself of order $\delta\rho$. For simplicity, we consider a central collision where the average density $\langle\rho({\bf s})\rangle$ is azimuthally symmetric, so that $\varepsilon_{n}$ is solely due to fluctuations. The noncentral case is discussed in Appendix A. In central collisions, there is no mean anisotropy, and the dominant contribution to $\varepsilon_{n}$ is of order $\delta\rho$. In order to evaluate the non-Gaussian fluctuations to leading order, we shall need to include the next-to-leading correction, of order $(\delta\rho)^{2}$. We show in Appendix B that the next-to-next-to-leading correction, of order $(\delta\rho)^{3}$, does not contribute to the quantities calculated in this article.

In order to achieve compact expressions, we introduce the following shorthand notations, for any function $f({\bf s})$ Gronqvist:2016hym :

[TABLE]

where $\langle E\rangle$ is the average energy per longitudinal length:

[TABLE]

Thus $\langle f\rangle$ is the average value of $f({\bf s})$ with the weight $\langle\rho({\bf s})\rangle$. With these notations, the center of the distribution is ${\bf s}_{0}\simeq\delta{\bf s}$, where we have made the approximation $\rho({\bf s})\simeq\langle\rho({\bf s})\rangle$ in the denominator of Eq. (2), which amounts to expanding ${\bf s}_{0}$ to leading order in $\delta\rho$. Expanding Eq. (3) to second order in $\delta\rho$, we then obtain:

[TABLE]

In each line, there is one term of order $\delta\rho$, which is the dominant contribution to the anisotropy, and was the only term kept in Blaizot:2014nia . In addition, there are two next-to-leading corrections of order $(\delta\rho)^{2}$, which originate from the non-linear relation between $\varepsilon_{n}$ and $\delta\rho$ Gronqvist:2016hym . The first correction is the effect of the fluctuation in the fireball size (denominator of Eq. (3)), while the last term originates from the recentering correction. Note that the recentering correction was neglected in Gronqvist:2016hym . We retain it in this paper because, as will be shown below, it plays a crucial role for the symmetric cumulant $SC(3,2)$.

III Definitions of cumulants

III.1 Cumulants of initial anisotropies

Anisotropic flow is not measured on an event-by-event basis. Experimentally-measured quantities are multi-particle correlations averaged over events, which yield average values of products of $v_{n}$ Bhalerao:2014xra . If $v_{n}$ is proportional to $\varepsilon_{n}$ in every event Gardim:2011xv , then, the relevant quantities for phenomenology are average values of products of $\varepsilon_{n}$, i.e., moments. Since $\varepsilon_{n}$ is a complex number whose phase is uniformly distributed by azimuthal symmetry, moments which are not invariant under rotations, such as $\langle\varepsilon_{n}\rangle$ or $\langle(\varepsilon_{n})^{2}\rangle$, vanish by symmetry. The lowest-order non-trivial moment is the rms value of $\varepsilon_{n}$, denoted by $c_{n}\{2\}$:

[TABLE]

where $\varepsilon_{n}^{*}$ denotes the complex conjugate, and angular brackets an average over events.333Note that we use the same angular brackets to denote an average over events, or an average value taken with the mean density profile, as in Eq. (4). There should be no confusion depending on the context.

Moments of order $\geq 4$ (note that all odd moments vanish by azimuthal symmetry) are usually combined with lower-order moments to form cumulants. Cumulants were originally introduced to suppress nonflow correlations in experimental analyses Borghini:2000sa , but they turn out to be useful observables also for studying flow fluctuations Voloshin:2007pc ; Yan:2013laa . In this article, we study cumulants of order 4, also called kurtosises, which are obtained from moments of order 4 by subtracting all pairwise correlations. There are cumulants involving only one harmonic $n$, denoted by $c_{n}\{4\}$ Borghini:2001vi :

[TABLE]

In the second equality, we have used azimuthal symmetry to simplify the expression. There is also a cumulant involving both $\varepsilon_{2}$ and $\varepsilon_{3}$, called symmetric cumulant, and denoted by $SC(3,2)$.444This quantity was denoted by $SC(3,2)_{\varepsilon}$ in ALICE:2016kpq to avoid confusion with the measured quantity, which involves $v_{n}$ rather than $\varepsilon_{n}$. It is defined in a similar way Bilandzic:2013kga :

[TABLE]

where, in the last equality, we have again eliminated the terms that vanish by symmetry. In this article, we refer to $SC(3,2)$ as to the “mixed kurtosis”, because it involves two different Fourier harmonics.

III.2 Cumulants of initial density

If one replaces $\varepsilon_{n}$ by its perturbative expansion, Eq. (7), moments and cumulants involve average values of products of $\delta\rho$, that is, correlation functions of the field $\rho({\bf s})$, which we now define.

The 2-point function, $C_{2}({\bf s}_{1},{\bf s}_{2})$, characterizes the variance of the fluctuations. It is defined by:

[TABLE]

where we used $\langle\delta\rho({\bf s})\rangle=0$. Similarly, the connected 3-point function, which characterizes the skewness of density fluctuations, is defined by:

[TABLE]

In this article, we shall also need the connected 4-point function, or kurtosis of density fluctuations, which is defined by subtracting out all pairwise correlations:

[TABLE]

For Gaussian density fluctuations Floerchinger:2014fta , $C_{3}$ and $C_{4}$ vanish. However, the positivity condition $\rho(x)\geq 0$ naturally generates positive cumulants to all orders Gronqvist:2016hym . One typically expects that the functions $C_{n}$ are short ranged (they vanish unless all the arguments are close to one another) and positive. However, long-range correlations can be induced by energy conservation, i.e., if one requires that all events have the exact same energy Blaizot:2014nia ; Gronqvist:2016hym .

III.3 $n$-point averages and orders of magnitude

In this article, we shall need to evaluate weighted products of $\delta\rho$, averaged over events, up to order $(\delta\rho)^{6}$. We now explain how such products are evaluated. For any function $f$, the one-point average $\langle\delta f\rangle$ vanishes by definition of $\delta\rho$. The first non-trivial average is the two-point average, which is obtained by multiplying Eq. (15) with functions of transverse coordinates and integrating over the transverse plane:

[TABLE]

The three-point average is obtained in a similar way using Eq. (16):

[TABLE]

Four-point averages are decomposed according to Wick’s theorem into products of two-point averages and a connected part, which we denote by the subscript $c$:

[TABLE]

With this notation, Eqs. (17) and (21) give immediately:

[TABLE]

For short-range correlations, one can make the approximation that the weighting functions $f,g,h,k$ in Eqs. (21), (22) and (25) are evaluated at the same point Gronqvist:2016hym , and integrate $C_{n}$ over relative positions Giacalone:2019kgg . Thus $\langle\delta f\delta g\rangle$ is just a weighted integral of $f({\bf s})g({\bf s})$, with a positive weight, and so on.

Let us now discuss orders of magnitude. If the weighting function $f$ is of order unity, then, $\delta f$ is typically much smaller than unity. Hence, the $n$-point averages are strongly ordered as a function of $n$. One must however distinguish odd and even moments. If $n$ is even, $(\delta f)^{n}$ is typically of the same order of magnitude as its average value over events $\langle(\delta f)^{n}\rangle$. For odd moments, there are cancellations upon averaging over events, so that $\langle(\delta f)^{2n-1}\rangle$ is of the same order as $\langle(\delta f)^{2n}\rangle$ for $n\geq 2$. For short-range correlations, the connected part $\langle(\delta f)^{n}\rangle_{c}$ is much smaller than $\langle(\delta f)^{n}\rangle$ for $n\geq 4$.

In this article, we shall need to evaluate 5-point and 6-point averages, but only to leading order. These leading-order contributions are the disconnected terms given by the Wick decomposition. The 5-point average is decomposed as a sum of products of 2-point averages and 3-point averages:

[TABLE]

while the dominant contribution to a 6-point average is from products of 2-point averages:

[TABLE]

IV Perturbative evaluation of cumulants

IV.1 General expressions

The perturbative expression of $c_{n}\{2\}$ in terms of the two-point function was derived in Blaizot:2014nia to leading order $(\delta\rho)^{2}$. It is obtained by inserting the first term in the right-hand side of Eqs. (7) into Eq. (9):

[TABLE]

In this article, we shall evaluate cumulants of order 4 to the first non-trivial order.

We start with the leading contribution to $\varepsilon_{n}$, corresponding to the first term in the right-hand side of Eq. (7), which is proportional to $\delta\rho$. One naively expects that its contribution to cumulants of order 4 is of order $(\delta\rho)^{4}$. However, this dominant contribution vanishes after averaging over events and subtracting the pairwise correlations in Eqs. (10) and (13). Only the subleading connected part, defined by Eq. (24), remains:

[TABLE]

The numerator in these equations is proportional to the kurtosis of the density field $C_{4}$, according to Eq. (25), which is the reason why we label it with the subscript $K$. Note that $c_{n}\{4\}_{K}$ and $SC(3,2)_{K}$ are positive for short-range correlations. Indeed, for the reasons explained in Sec. III.3, the numerators in the right-hand side of Eq. (31) are integrals of $|{\bf s}|^{4n}$ and $|{\bf s}|^{10}$, respectively, with a positive weight.

If the anisotropies depended linearly on the initial density $\rho$, the kurtosises of the initial anisotropy would be proportional to $C_{4}$, the kurtosis of density fluctuations. However, the nonlinearity, which results in the second-order terms in Eq. (7), generates other contributions to the kurtosises which are of the same order as $c_{n}\{4\}_{K}$ and $SC(3,2)_{K}$. These contributions originate from terms of order $(\delta\rho)^{5}$ and $(\delta\rho)^{6}$ in the expansion of the product of the four $\varepsilon_{n}$ factors in Eqs. (10) and (13).

We now evaluate the contribution of order $(\delta\rho)^{5}$, which is obtained by keeping the correction of order $(\delta\rho)^{2}$ for just one of the four factors, and only the terms of order $\delta\rho$ for the other three factors. The corresponding contribution to the product of $\varepsilon_{n}$ is of the type $\delta_{2}\delta_{2}\delta_{1}\delta_{1}\delta_{1}$, where we use the shorthand notation $\delta_{1}$ for the leading contribution to $\varepsilon_{n}$ or $\varepsilon_{n}^{*}$, and $\delta_{2}\delta_{2}$ for one of the two subleading contributions. When averaging over events, this 5-point average is split according to Wick’s theorem, Eq. (27). The contractions of the type $\langle\delta_{2}\delta_{2}\rangle\langle\delta_{1}\delta_{1}\delta_{1}\rangle$ or $\langle\delta_{2}\delta_{2}\delta_{1}\rangle\langle\delta_{1}\delta_{1}\rangle$ cancel out when subtracting the disconnected parts in Eqs. (10) and (13). Only the contractions of the type $\langle\delta_{2}\delta_{1}\rangle\langle\delta_{2}\delta_{1}\delta_{1}\rangle$ remain. Using azimuthal symmetry to simplify the expressions, one obtains:

[TABLE]

where we use the subscript $S$ because this contribution involves three-point averages, which are proportional to the skewness of the density field $C_{3}$, according to Eq. (22). The only contribution to $c_{n}\{4\}_{S}$ is from the second term in the right-hand side of Eq. (7), i.e., the fluctuation in the size. The recentering correction does not contribute to $c_{n}\{4\}_{S}$, but it is the origin of the last term in $SC(3,2)_{S}$. Note that $c_{n}\{4\}_{S}$ and $SC(3,2)_{S}$ are always negative for short-range correlations, contrary to $c_{n}\{4\}_{K}$ and $SC(3,2)_{K}$ which are positive.

Finally, we consider contributions of order $(\delta\rho)^{6}$, whose average over events are evaluated using Eq. (29). Such contributions can arise in two different ways. The first possibility is to have a correction of order $(\delta\rho)^{3}$ to one of the four $\varepsilon_{n}$ factors, so that the product is of the type $(\delta_{3}\delta_{3}\delta_{3})\delta_{1}\delta_{1}\delta_{1}$. We show in Appendix B that these terms do not contribute to the kurtosises. The other contributions of order $(\delta\rho)^{6}$ arise from having a correction of order $(\delta\rho)^{2}$ to two of the four $\varepsilon_{n}$ factors, so that the product is of the type $(\delta_{2}\delta_{2})(\delta^{\prime}_{2}\delta^{\prime}_{2})\delta_{1}\delta_{1}$. The contractions $\langle\delta_{2}\delta_{2}\rangle\langle\delta^{\prime}_{2}\delta^{\prime}_{2}\rangle\langle\delta_{1}\delta_{1}\rangle$ or $\langle\delta_{2}\delta^{\prime}_{2}\rangle\langle\delta_{2}\delta^{\prime}_{2}\rangle\langle\delta_{1}\delta_{1}\rangle$ cancel out in the subtraction of Eqs. (10) and (13), while the contractions $\langle\delta_{1}\delta^{\prime}_{2}\rangle\langle\delta_{1}\delta^{\prime}_{2}\rangle\langle\delta_{2}\delta_{2}\rangle$ or $\langle\delta_{1}\delta_{2}\rangle\langle\delta_{1}\delta_{2}\rangle\langle\delta^{\prime}_{2}\delta^{\prime}_{2}\rangle$ vanish by azimuthal symmetry. Eventually, only the contractions $\langle\delta_{1}\delta_{2}\rangle\langle\delta_{1}\delta^{\prime}_{2}\rangle\langle\delta_{2}\delta^{\prime}_{2}\rangle$ contribute to the cumulant. One obtains

[TABLE]

where the subscript $V$ means that this contribution only involves two-point averages (21), i.e., the variance $C_{2}$ of the density field. The only contribution to $c_{n}\{4\}_{V}$ comes from the size fluctuation (second term in the right-hand side of Eqs. (7)). The first contribution to $SC(3,2)_{V}$ is also from the size fluctuation, while the second is from the recentering correction to $\varepsilon_{3}$ (last term in the second line of Eq. (7)). Both $c_{n}\{4\}_{V}$ and $SC(3,2)_{V}$ are positive.

The full expressions of the cumulants to leading order in the perturbative expansion are obtained by summing the contributions of Eqs. (31), (33), (37):

[TABLE]

Note that the recentering correction does not contribute to $c_{n}\{4\}$ (which is the reason why it was neglected in Gronqvist:2016hym ), but it contributes to $SC(3,2)$.

As we shall illustrate in the specific case of identical pointlike sources Bhalerao:2006tp , $V$ and $K$ contributions are positive, while the $S$ contribution is negative. Therefore, the experimental observation that both $c_{3}\{4\}$ and $SC(3,2)$ are negative implies that the $S$-type contribution, which is induced by the 3-point function of the density field, dominates. We come back to this point in Sec. V.2.

IV.2 Identical pointlike sources

We now simplify the results of Sec. IV.1 in the specific case where the energy density is a sum of Dirac peaks Gronqvist:2016hym :

[TABLE]

where $N$ is a fixed number, and ${\bf s}_{j}$ are independent variables, whose probability distribution is defined by the average density $\langle\rho({\bf s})\rangle$. In this case, Eqs. (21), (22) and (25) reduce to Alver:2008zza :

[TABLE]

where we have used the shorthand $\hat{f}\equiv f-\langle f\rangle$. For pointlike sources, the limit of large $N$ is also the limit of small fluctuations. The perturbative expansion of Sec. II is equivalent to an expansion in $1/N$ Alver:2008zza .

For pointlike sources, Eq. (30) simplifies to:

[TABLE]

where $r\equiv|{\bf s}|$. Adding the contributions of Eqs. (31), (33) and (37), one obtains:

[TABLE]

The expression of $c_{n}\{4\}$ coincides with that derived in Alver:2008zza for $n=2$, and then extended in Bhalerao:2011bp to $n=3$. The result for $SC(3,2)$ is new. These equations illustrate that these quantities depend in a non-trivial way on the initial density profile. In particular, there are positive and negative terms, which are typically of the same order of magnitude. More importantly, the above results illustrate the symmetry between $c_{n}\{4\}$ and $SC(3,2)$. In the next section, we introduce scaled observables which allow to compare their magnitudes.

V Measures of non-Gaussianity

The natural dimensionless observable is, rather than the kurtosis itself, the kurtosis divided by the square of the variance, called standardized kurtosis. By scaling $c_{n}\{4\}$ and $SC(3,2)$, one obtains:

[TABLE]

One advantage of these quantities is that if $v_{n}$ is proportional to $\varepsilon_{n}$ in every event, the proportionality constant cancels out in the ratio. Therefore, the standardized kurtosis of the distribution of $v_{n}$ is that of the distribution of $\varepsilon_{n}$, which allows one to compare directly models of initial conditions to experimental data. Another advantage is that since they are “standardized”, we can perform meaningful comparisons between their magnitudes.

In Sec. V.1, we check the validity of the perturbative results of Sec. IV.2 using Monte Carlo simulations. In Sec. V.2, we compare several measures of primordial non-Gaussianity, and discuss the implications of existing data.

V.1 Monte Carlo simulations

We carry out Monte Carlo simulations using the identical source model of Sec. IV.2, in which the average density profile is a symmetric Gaussian in the transverse plane. In this case, using Eqs. (47) and (48), one obtains

[TABLE]

Note that the standardized kurtosises are proportional to $1/N$. For large $N$, fluctuations are approximately Gaussian according to the central limit theorem, which is the reason why the standardized kurtosis is small. While the $1/N$ behaviour is generic, the dimensionless coefficient in front is sensitive to the density profile. This sensitivity is enhanced by the fact that there are positive and negative terms of the same order of magnitude in Eq. (48): The sum is typically smaller than any individual term.

Figure 1 displays our Monte Carlo results for the standardized kurtosis, together with the analytic results (52). Since Eq. (52) is the leading-order result in $1/N$, one expects that it is valid for large $N$. This is confirmed by the results of Fig. 1, where numerical results converge to the curves as $N$ increases. The open symbols in Fig. 1 are obtained by switching off the recentering correction, i.e., by setting ${\bf s}_{0}=0$ in Eq. (3). As expected from the perturbative calculation of Sec. IV, this has little effect on $c_{2}\{4\}$ and $c_{3}\{4\}$. On the other hand, it has a dramatic effect on $sc(3,2)$, which changes sign when the recentering correction is switched off. In the perturbative calculation, the recentering correction corresponds to the fourth and sixth terms of the second line of Eq. (48), with factors $-6$ and $9$ in front. If one omits this term, one obtains $sc(3,2)=9/(4N)$ for a Gaussian profile. This value is displayed as a dashed line in Fig. 1. It agrees with the corresponding Monte Carlo results for large $N$, as expected.

Our perturbative calculations are carried out for central collisions, where the mean density is azimuthally symmetric. In order to test how results are modified for noncentral collisions, we carry out Monte Carlo calculations for an asymmetric Gaussian profile, whose widths along $x$ and $y$ differ. We define the mean eccentricity by:

[TABLE]

Results are presented in Fig. 2 as a function of $\bar{\varepsilon}_{2}$. The kurtosis of $\varepsilon_{3}$ (panel (a)) and the mixed kurtosis (panel (b)) depend only mildly on $\bar{\varepsilon}_{2}$. On the other hand, the kurtosis of $\varepsilon_{2}$ (panel (c)) increases by orders of magnitude. The reason is that it is driven by the mean eccentricity, not by fluctuations. The value $-1$ corresponds to the limit where fluctuations are negligible, and $\varepsilon_{2}\simeq\bar{\varepsilon}_{2}$. Therefore, the standardized kurtosis (50) is not an appropriate measure of non-Gaussian fluctuations for elliptic flow. Fluctuations of $\varepsilon_{2}$ in non-central collisions are discussed thoroughly in Appendix A, where we show that a measure of non-Gaussian $\varepsilon_{2}$ fluctuations, of the same order as the standardized kurtosis, is the scaled skewness defined by:

[TABLE]

While this quantity is typically negative, its negative sign cannot be attributed to the positive skewness of the underlying density field, unlike the negative signs of $c_{3}\{4\}$ and $SC(3,2)$.

V.2 Experimental data

Figure 3 displays three measurements of non-Gaussianity of $\varepsilon_{n}$ fluctuations: the scaled skewness of elliptic flow fluctuations, the kurtosis of triangular flow Abbasi:2017ajp , and the mixed kurtosis. Based on the perturbative expansion, one typically expects these quantities to be much smaller than unity, and roughly of the same order or magnitude. Their magnitude should increase with the centrality percentile. Based on the results obtained for a Gaussian density profile, Eqs. (78) and (52), one also expects them to be negative. These trends are confirmed by existing data to some extent. The scaled skewness has the expected order of magnitude and centrality dependence. The centrality dependence of the two standardized kurtosises is somewhat weaker than one would expect based on the $1/N$ law, taking for $N$ the number of wounded nucleons as estimated in a Glauber model Miller:2007ri . This is confirmed by Monte Carlo calculations of the kurtosis of $\varepsilon_{3}$ Bhalerao:2011ry ; Giacalone:2017uqx which give values in rough agreement with data, but with a stronger centrality dependence. The striking observation in Fig. 3 is that the mixed kurtosis is larger in absolute value than the kurtosis of $v_{3}$, while it is a factor 10 smaller for a Gaussian density profile, Eq. (52). Thus, the measured value of $sc(3,2)$ seems unusually large. (Similar values were recently reported by the ATLAS collaboration Aaboud:2019sma .) This trend is confirmed if one compares with standard models of initial conditions ALICE:2016kpq or full hydrodynamic calculations Gardim:2016nrrNiemi:2015qia; Acharya:2017gsw , which predict smaller values. However, values of $SC(3,2)$ similar to those observed in ALICE have been reported Alba:2017hhe ; Zhao:2017yhj ; Moreland:2018gsh in recent hydrodynamic calculations using the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo model of initial conditions Moreland:2014oya . We stress that the negative sign of the kurtosises implies that the three-point function $C_{3}$ of the density field is positive, as it gives the only negative contribution to the skewness (see Sec. IV.1).

VI Conclusions

We have presented three different measures of non-Gaussian anisotropy fluctuations: The scaled skewness of $v_{2}$ fluctuations, the kurtosis of $v_{3}$ fluctuations, and the mixed kurtosis between $v_{2}$ and $v_{3}$. In hydrodynamics, these observables can be directly evaluated from the model of initial conditions, using the proportionality between $v_{n}$ and $\varepsilon_{n}$.555It has been found in extensive hydrodynamic calculations Abbasi:2017ajp that this proportionality breaks down, and that the kurtosis of $v_{3}$ is less negative than the kurtosis of $\varepsilon_{3}$. Such departures from linear hydrodynamic response deserve further investigations. We have argued that these three quantities, which measure deviations to the central limit theorem, are typically much smaller than unity, of the same order of magnitude and with a similar centrality dependence — their magnitude becomes larger as the centrality percentile increases. These trends are to some extent confirmed by existing data.

We have carried out a perturbative expansion of anisotropies in terms of the fluctuations of the initial energy density field. We have applied this perturbative scheme only to a few observables, but it is systematic and could be applied to any initial-state quantity, in particular participant-plane correlations Jia:2012ju .

We have evaluated the kurtosises to leading order in perturbation theory. Even to leading order, we find that they involve the two-, three- and four-point functions of the density field, which makes their interpretation difficult. Nonetheless, they receive a negative contributions from the three-point function only. Hence, the observation that the kurtosises are both negative implies that the three-point function of the density field is positive, and large enough that it overcomes the other contributions. The scaled skewness of elliptic flow fluctuations involves the two- and three-point functions of the density field, but the sign of each contribution depends on the density profile. Therefore, it represents a less stringent probe of primordial non-Gaussianity.

Acknowledgments

RSB would like to acknowledge the hospitality of the IPhT, Saclay, France where a part of this work was done and the support of the CNRS LIA (Laboratoire International Associé) THEP (Theoretical High Energy Physics) and the INFRE-HEPNET (IndoFrench Network on High Energy Physics) of CEFIPRA/IFCPAR (Indo-French Center for the Promotion of Advanced Research). RSB also acknowledges the support of the Department of Atomic Energy, India for the award of the Raja Ramanna Fellowship.

Appendix A $\varepsilon_{2}$ fluctuations in non-central collisions

The overlap area of two identical nuclei in a noncentral collision is almond shaped, hence, the mean energy density profile $\langle\rho({\bf s})\rangle$ has a large elliptic deformation. In this situation, $c_{2}\{4\}$ no longer represents the kurtosis of $\varepsilon_{2}$ fluctuations, and signatures of primordial non-Gaussianity are less easy to identify.

In this Appendix, we derive the general expression of $c_{2}\{4\}$ for non-central collisions. We thus generalize the result obtained by Alver et al. (Eq. (B37) of Alver:2008zza ) to a continuous density profile. We also provide a transparent physical interpretation of the various terms in their result. We then identify a signature of non-Gaussian $\varepsilon_{2}$ fluctuations, which we relate to experimentally-measured quantities.

We denote by $\bar{\varepsilon}_{2}$ the value of $\varepsilon_{2}$ obtained by replacing $\rho({\bf s})$ with $\langle\rho({\bf s})\rangle$ in Eq. (3):

[TABLE]

We choose a coordinate frame where the impact parameter is along the $x$ axis, so that the system has $y\to-y$ symmetry and $\bar{\varepsilon}_{2}$ is real (it is negative with the convention chosen in Eq. (57)). We decompose $\varepsilon_{2}$ into a (real) mean value and a (complex) fluctuation:

[TABLE]

The magnitude of $\varepsilon_{2}$ fluctuations is characterized by its rms value $\sigma$:

[TABLE]

whose perturbative expression was derived in Bhalerao:2019uzw . In this Appendix, we evaluate the asymmetry of the fluctuations, which we define by

[TABLE]

and the genuine non-Gaussian fluctuations, characterized by the skewness

[TABLE]

and the kurtosis

[TABLE]

Inserting Eq. (58) into Eq. (10) and expanding in powers of $\delta\varepsilon_{2}$, one obtains Bhalerao:2018anl :

[TABLE]

For central collisions, $\langle\varepsilon_{2}\rangle$ vanishes and only the kurtosis remains, so that $K$ in this equation corresponds to $c_{2}\{4\}$ calculated in Sec. IV. We consider that $K$ is not sensitive to the almond shape of the profile, and that it has the same value as for central collisions. Note that $A$ and $S$ vanish by symmetry for central collisions. For noncentral collisions, simplifications occur if the fluctuations are both symmetric ($A=0$) and Gaussian ($S=K=0$). One then obtains Voloshin:2007pc $c_{2}\{4\}=-\langle\varepsilon_{2}\rangle^{4}$.

In this Appendix, we refine this result by deriving the expressions of $A$ and $S$ to leading order in perturbation theory. In the identical source model, $A$, $S$ and $K$ are proportional to $1/N$, $1/N^{2}$ and $1/N^{3}$, respectively. Thus, the various terms in Eq. (63) correspond to the successive orders in the $1/N$ expansion of Alver:2008zza .666More precisely, the $1/N$ term is the sum of two contributions, where the first contribution is the asymmetry $A$, and the second contribution comes from the difference between $\langle\varepsilon_{2}\rangle$ and $\bar{\varepsilon}_{2}$ Bhalerao:2019uzw .

We shall neglect the recentering correction to $\varepsilon_{2}$. Its contribution to the quantities we evaluate is not strictly zero for non-central collisions Bhalerao:2019uzw , but it vanishes for identical pointlike sources, and we anticipate that it is always small. The presence of a mean eccentricity modifies Eq. (7) to Bhalerao:2019uzw :

[TABLE]

There are two terms of order $\delta\rho$ and two terms of order $(\delta\rho)^{2}$. Note that the average over events $\langle\varepsilon_{2}\rangle$ gets a contribution from second-order terms, and therefore differs from $\bar{\varepsilon}_{2}$ Bhalerao:2006tp ; Bhalerao:2019uzw . We evaluate the asymmetry defined by Eq. (60), which can be rewritten as:

[TABLE]

We evaluate $A$ to the first non-trivial order $(\delta\rho)^{2}$. Terms which are already of this order in Eq. (64) cancel out in the difference of Eq. (65). Therefore, only terms of order $\delta\rho$ in Eq. (64) contribute to $A$. One immediately obtains

[TABLE]

The numerator is a sum of two-point averages, which is proportional to the two-point function of the density field according to Eq. (21). For identical pointlike sources, Eq. (66) can be simplified along the lines of Sec. IV.2:

[TABLE]

It is typically of order $(\bar{\varepsilon}_{2})^{2}/N$ Alver:2008zza .

We now evaluate the skewness $S$, which is the connected part of the moment of order 3, to leading order in the fluctuations. The terms of order $(\delta\rho)^{3}$ and $(\delta\rho)^{4}$ contribute to the same order after averaging over fluctuations. We denote by $S_{S}$ the contribution of order $(\delta\rho)^{3}$ to the skewness. It is obtained by keeping only terms of order $\delta\rho$ in Eq. (64). We neglect terms such as $\langle(\delta{\bf s}^{2})^{3}\rangle$ which are of order $(\bar{\varepsilon}_{2})^{3}$ and only retain the terms of order $\bar{\varepsilon}_{2}$:

[TABLE]

The numerator is a sum of 3-point averages, which are evaluated using Eq. (22). It is directly proportional to the 3-point function of the density field $C_{3}$, as the contribution (33) to the kurtosises.

For identical pointlike sources, Eq. (68) gives:

[TABLE]

The two terms have opposite signs. For an asymmetric Gaussian density profile, they mutually cancel. The sign of this contribution depends on the details of the density profile, and there is no definite effect of the skewness of the density field on the skewness of $\varepsilon_{2}$ fluctuations. This is at variance with the corresponding contribution to the kurtosises, Eq. (33), which is always negative.

We now evaluate the contribution of order $(\delta\rho)^{4}$ to the skewness, which we denote by $S_{V}$. It is obtained by taking into account only the terms of order $(\delta\rho)^{2}$ in one of the three $\varepsilon_{2}$ factors in Eq. (61), and only the terms of order $\delta\rho$ in the two remaining factors. One is led to evaluate contributions of the type $\langle\delta_{2}\delta_{2}\delta_{1}\delta_{1}\rangle$, where we write schematically the order $(\delta\rho)^{2}$ term as $\delta_{2}\delta_{2}$ and the order $\delta\rho$ terms as $\delta_{1}$. This four-point average is decomposed according to Wick’s theorem. The contractions of the type $\langle\delta_{2}\delta_{2}\rangle\langle\delta_{1}\delta_{1}\rangle$ cancel out when subtracting disconnected parts. Only the contractions of the type $\langle\delta_{2}\delta_{1}\rangle\langle\delta_{2}\delta_{1}\rangle$ remain. One obtains:

[TABLE]

The numerator only contains two-point averages. Therefore, this contribution only involves the 1- and 2-point functions of the density field, as the contribution (37) to the kurtosises. For identical pointlike sources, this expression gives:

[TABLE]

The total skewness is obtained by adding the contributions from Eqs. (68) and (70):

[TABLE]

For pointlike sources, using Eqs. (69) and (71), one recovers the $1/N^{2}$ contribution to $c_{2}\{4\}^{4}$, which was derived in Alver:2008zza (Eq. (B37)). We have seen that $S_{S}$ can have either sign depending on the profile. $S_{V}$ again contains two terms of opposite signs, but the first typically dominates, resulting in a negative $S_{V}$. For an asymmetric Gaussian profile, one obtains

[TABLE]

Eventually, we relate the skewness, which is the leading contribution to non-Gaussian $\varepsilon_{2}$ fluctuations in non-central collisions, to experimental observables. We construct a dimensionless observable, $\Sigma$, which is proportional to $S$ and of the same order of magnitude as the standardized kurtosises (50), i.e., $1/N$ (Eq. (52)). $S$ can be extracted from the small splitting between $\varepsilon_{2}\{6\}$ and $\varepsilon_{2}\{4\}$ Giacalone:2016eyu :

[TABLE]

where we have used approximate azimuthal symmetry which implies $\langle\delta\varepsilon_{x}^{3}\rangle=3\langle\delta\varepsilon_{x}\delta\varepsilon_{y}^{2}\rangle$. The quantity (74) is of order $1/(N^{2}\bar{\varepsilon}_{2})$. We next use the approximate equality $\varepsilon_{2}\{4\}\simeq\varepsilon_{2}\{6\}\simeq\bar{\varepsilon}_{2}$ to write:

[TABLE]

This quantity is now of order $1/N^{2}$. Eventually, we make it dimensionless by dividing by the variance of $\varepsilon_{2}$ fluctuations $\sigma^{2}$:

[TABLE]

This quantity is of order $1/N$. If $v_{2}$ is proportional to $\varepsilon_{2}$, it can be extracted from data using

[TABLE]

We dub this quantity the scaled skewness, to distinguish it from the standardized skewness Giacalone:2016eyu ; Sirunyan:2017fts ; Acharya:2018lmh , which is of different order, $\bar{\varepsilon}_{2}/N^{1/2}$ Bhalerao:2018anl . For a Gaussian density profile, $\sigma^{2}=2/N$. Using Eqs. (76) and (73), one obtains:

[TABLE]

Figure 4 displays a comparison between Monte Carlo calculations and this perturbative estimate, for two values of the mean eccentricity $\bar{\varepsilon}_{2}$. One sees that the scaled skewness is independent of $\bar{\varepsilon}_{2}$, and has limited sensitivity to the recentering correction, which we have neglected throughout this Appendix. Agreement with the perturbative result (78) is good for large $N$ as expected. It is better for the smaller value of $\bar{\varepsilon}_{2}$.

Appendix B Expansion of $\varepsilon_{n}$ to order $(\delta\rho)^{3}$

Equations (7) give the expansion of Eq. (3) up to order $(\delta\rho)^{2}$. If one includes the next order, $(\delta\rho)^{3}$, one obtains:

[TABLE]

Note that there are more terms of order $(\delta\rho)^{3}$ for $\varepsilon_{3}$ than for $\varepsilon_{2}$. In the last line of Eq. (85), $\delta 1$ denotes the relative fluctuation of the energy, $\delta 1=\int_{\bf s}\delta\rho/\langle E\rangle$, according to the notation of Eq. (4). This quantity appears when expanding Eq. (2) up to order $(\delta\rho)^{2}$.

We now explain why the terms of order $(\delta\rho)^{3}$ do not contribute to the kurtosises defined by Eq. (10) and (13). When evaluating a product of four $\varepsilon_{n}$ factors up to order $(\delta\rho)^{6}$, one must include terms where one of the $\varepsilon_{n}$ factors is expanded up to order $(\delta\rho)^{3}$, while all the other factors are evaluated to leading order $\delta\rho$. Such terms are of the form $(\delta_{3}\delta_{3}\delta_{3})\delta_{1}\delta_{1}\delta_{1}$, and their average over events is evaluated using Eq. (22). Contractions of the type $\langle\delta_{3}\delta_{3}\rangle\langle\delta_{3}\delta_{1}\rangle\langle\delta_{1}\delta_{1}\rangle$ cancel out when subtracting disconnected terms. Only contractions of the type $\langle\delta_{3}\delta_{1}\rangle\langle\delta_{3}\delta_{1}\rangle\langle\delta_{3}\delta_{1}\rangle$ contribute to the kurtosises. However, one easily sees that there is always one factor $\langle\delta_{3}\delta_{1}\rangle$ which vanishes by azimuthal symmetry, so that these contractions also vanish.