Statistics of initial density perturbations in heavy ion collisions and their fluid dynamic response

An interesting opportunity to determine thermodynamic and transport properties in more detail is to identify generic statistical properties of initial density perturbations. Here we study event-by-event fluctuations in terms of correlation functions for two models that can be solved analytically. The first assumes Gaussian fluctuations around a distribution that is fixed by the collision geometry but leads to non-Gaussian features after averaging over the reaction plane orientation at non-zero impact parameter. In this context, we derive a three-parameter extension of the commonly used Bessel-Gaussian event-by-event distribution of harmonic flow coefficients. Secondly, we study a model of N independent point sources for which connected n-point correlation functions of initial perturbations scale like 1/N^(n-1). This scaling is violated for non-central collisions in a way that can be characterized by its impact parameter dependence. We discuss to what extent these are generic properties that can be expected to hold for any model of initial conditions, and how this can improve the fluid dynamical analysis of heavy ion collisions.

In practice, testing QCD thermodynamics experimentally is complicated by the fact that data result from a convoluted time history that depends not only on hydrodynamic evolution but also on initial conditions and hadronization. In particular, essentially all flow measurements are correlation measurements and correlations can be present already in the initial conditions, or they can arise (or be attenuated) dynamically during the hydrodynamic evolution or during hadronization, respectively. The theory framework for determining the dynamical evolution is clear: once the thermodynamic information and transport properties entering the equations of dissipative fluid dynamics are specified, the propagation of a known initial condition can be controlled. A significant number of tools and techniques has been developed to this end [32,33]. However, the initial conditions are arguably less controlled so far. Understanding their dynamical origin and their statistical properties is now becoming a major focus of current research.
Earlier efforts in this direction have focussed mainly on exploring the density distributions generated in Monte Carlo models that implement variants of the optical Glauber model [13,[34][35][36][37] or supplement these with effects from parton saturation physics [19,[38][39][40], or in dynamically more complete code-based formulations [41,42]. In a different direction, the simplifying assumption that initial density perturbations follow a Gaussian distribution has served since long as a baseline for characterizing initial conditions [44][45][46]. The question arises to what extent such studies explore only a possibly limited range of the total parameter space of conceivable initial conditions, or whether it is possible to identify universal features that any phenomenologically relevant model of initial conditions is expected to satisfy on general grounds. One such general consideration that applies to heavy ion collisions is that particle production arises from a large number of essentially independent sources with identical statistical properties. It is well-known in probability theory [47] that this alone implies that n-th order cumulants (or connected n-point correlation functions) of variables that are normalized sums of these independent source contributions scale in a characteristic way ∝ 1/N n−1 with the number N of sources. Different eccentricities m {n} have been calculated for such a model to various orders in 1/N by Bhalerao and Ollitrault [48] as well as Alver et al. [10] and where shown to quantitatively reproduce results of more sophisticated Glauber models in nucleus-nucleus collisions [49]. First indications that the scaling with n is particularly relevant for initial conditions in proton-nucleus and nucleus-nucleus collisions go back to numerical findings of Bzdak, Bozek and McLerran [50]. They were sharpened subsequently due to work of Ollitrault and Yan [51] (see also Bzdak and Skokov [52]) who established a related scaling for eccentricity cumulants at vanishing impact parameter in an analytically accessible model of independent point sources (IPSM) and who showed that this reproduces with good numerical accuracy the eccentricity cumulants in other currently used models of initial conditions. The present paper aims at contributing to this important recent development. To this end, we shall show how one can solve the IPSM completely, including the set of n-point correlation functions that characterize completely the information about the radial and azimuthal dependence at zero and non-zero impact parameter. Based on this differential information, we shall provide further evidence that the IPSM shares indeed important commonalities with realistic model distributions. At finite impact parameter b, we shall find that the 1/N n−1 -scaling is broken for azimuthally averaged event samples. However, for small b, the leading b-dependence of the terms that break this scaling can be given analytically. Thus information about this b-dependence, combined with information about the 1/N n−1 -scaling for b = 0 can provide an ordering principle that applies more generally to n-point correlators at zero and non-zero impact parameter. We shall also discuss how the connected n-point correlation functions of initial fluctuations enter the calculation of measurable correlators of flow coefficients, and we shall point to possible further phenomenological applications of these insights. The main assumption underlying the scaling of connected n-point correlation functions with 1/N n−1 is that the transverse density is given by a sum of N independent and identically distributed random variables or functions of random variables. 1 As we discuss in more detail in the main text, this holds also for ensembles of non-central events, however, only if impact parameter and reaction plane orientation are kept fixed. In contrast, the phenomenologically relevant connected n-point correlation functions are defined for ensembles with random azimuthal orientation. To cope with this complication, we find it useful to work in a framework sketched in Fig. 1: we denote event averages with fixed azimuthal orientation by . . . and we construct moments and the corresponding cumulants as usual from a generating functional and its logarithm, respectively. Randomizing the azimuthal orientation φ R in the averages . . . defines the average . . . • . The scaling with 1/N n−1 is broken for the ensemble average . . . • at finite impact parameter, since the operation of averaging over φ R does not commute with the operation of passing from moments to cumulants. In other words, the cumulants with respect to the randomized ensemble do not correspond to φ R -averages of cumulants evaluated at fixed φ R .
As a significant part of this paper will study in detail the independent point-sources model, we conclude this introduction by asking to what extent the spatial dependence of correlation functions in the IPSM can be expected to have physical significance. One may argue that the long-wavelength excitations (small values of azimuthal wave numbers m and radial wave numbers l in a Bessel-Fourier expansion) do not resolve the differences between a spatially extended but short range source function and a point-like source. Since these long wavelength modes are most important for the fluid dynamic evolution (others get damped quickly by dissipative effects), one might expect that also some spacedependent features of the independent point-sources model contain realistic aspects. They are universal in the sense that a larger class of models with extended sources (and even some early non-equilibrium dynamics as long as it is local) lead to equal correlation functions for the long wavelength modes. If it could be established, such a universality for the correlations of the most important fluid dynamic modes would have profound consequences. For instance, in a mode-by-mode fluid dynamics framework one could use this knowledge of initial conditions for a detailed comparison between experimental results on correlations of harmonic flow coefficients and fluid dynamic calculations which would allow for a more detailed determination of thermodynamic and transport properties. These are some of the considerations that have prompted the following analysis.

Flow cumulants
In this section, we discuss how flow measurements are related to the n-mode correlation functions of initial density perturbations that we are going to analyze in sections 3 and 4 below. To focus on the structure of this relation, we shall defer some technical definitions to section 3. We start from a perturbative expansion of the complex-valued event-wise flow coefficients in powers of weights w (m) l that characterize these density perturbations in terms of azimuthal (m) and radial (l) wave numbers [53] Here, the indices m i are summed over the range (−∞ . . . ∞) and the indices l i are summed over the range (1, . . . , ∞). As coefficients of a Bessel-Fourier expansion, defined in eq. (3.2) below, the w The dynamical response functions S (m 1 ,...,mn)l 1 ,...,ln satisfy then S (m 1 ,...,mn)l 1 ,...,ln = (−1) m 1 +...+mn S * (−m 1 ,...,−mn)l 1 ,...,ln . For the harmonic flow coefficients one has V −m = V * m . In general, n-th order flow cumulants v m {n} n denote the connected n-point event average of flow coefficients V m . The lowest order cumulants take the explicit form [54,55] These higher order flow cumulants are measured in ion-ion and in proton-ion collisions [1,56,57]. With the help of the perturbative expansion (2.1), one can write flow cumulants as products of event averages of initial fluctuating modes w The linear response term of (2.7) can be written in terms of a connected four-point function of initial fluctuations 3 , (Summation over the indices l 1 , . . . , l 4 is implied here and in the following.) In general, the linear response contribution to v m {2n} 2n is proportional to a connected (2n)-mode correlator Flow measurements are not limited to the determination of flow cumulants. In principle, arbitrary event averages V m 1 V m 2 . . . V mn • of products of flow coefficients are experimentally accessible, see e. g. [58]. For event samples with randomized orientation of the reaction plane, the simplest generalization are 3-flow correlators To be specific, let us write here the expansion of one of them, On the right hand side we have included terms from linear dynamics as well as those quadratic corrections that contain four point functions with two opposite index pairs (m, −m). These are the leading contributions for ensembles that are close to Gaussian. Other experimentally easily accessible 3-flow correlators include V 2 V 2 V * 4 • and V 3 V 3 V * 6 • , for which similar expansions can be written down. The dynamical response functions that appear on the right hand side of (2.10) can be found also in the expansion of the flow cumulants (2.9). 4 But in the 3-flow correlators, the linear and non-linear dynamic response terms are weighted with a different set of informations about the initial conditions, namely a different set of moments w • that typically involve harmonic modes with different m.
As illustrated by the examples discussed so far, the calculation of flow correlation measurements V m 1 V m 2 V m 3 . . . • requires knowing the initial n-mode correlators w • and the dynamical response functions S (m 1 ,...,mn)l 1 ,...,ln . We note that the dynamical response functions are known in principle, in the sense that they are calculable once the thermodynamic information entering hydrodynamic evolution and the event-averaged initial enthalpy density is given. No further model dependent assumption enters their calculation. A method of how to determine them numerically was given in Ref. [53]. On the other hand, the correlators w • should be calculable in principle from a microscopic theory of thermalization dynamics. In practice, however, this program is not yet carried out, and the initial conditions are currently regarded as the most significant source of uncertainties in the calculation of flow observables. This motivates us to investigate in the following what can be said on the basis of general considerations about the structure of n-mode correlators w

Gaussian probability distributions of initial conditions
In this section, we introduce Gaussian probability distributions of fluctuations in the initial transverse enthalpy density w( x), and we discuss their implications for flow cumulants and flow probability distributions.

Gaussian model of initial fluctuations for fixed reaction plane angle φ R
We start from the general form of a Gaussian probability distribution of the enthalpy density written for fixed impact parameter and reaction plane angle φ R (see also appendix C of ref. [46]) As a Gaussian distribution it is specified completely in terms of the expectation valuew( x) and the connected two-point correlation function C( x, y), which is the inverse of M ( x, y) seen as a matrix of infinite dimension with indices x and y. For an arbitrary event, we write the enthalpy density in a Bessel-Fourier expansion Here, ρ(r) is a monotonous function that maps r ∈ (0, ∞) to ρ ∈ (0, 1). It is specified in appendix A. The real numbers z The expectation value at fixed impact parameter and reaction plane angle φ R can be written in the same Bessel-Fourier expansion, Here, the sum over m on the right hand side goes only over the even values m = ±2, ±4, . . . as it follows from the discrete symmetry thatw(r, φ) =w(r, φ + π). The function w BG (r) is defined such that the m = 0 component in the sum vanishes. The dimensionless and real coefficientsw (m) l depend on centrality and they vanish with vanishing impact parameter b, i.e., for ultra-central collisions. One can show that for small b they behave like w (m) l ∼ b |m| , see appendix C. For the two-point correlation function we write a Bessel-Fourier expansion in terms of the coefficients C . Note that Eq. (3.4) contains a factor e −i(m 1 +m 2 )φ R such that the right hand side vanishes when averaged over the reaction plane angle φ R with uniform distribution, except for m 1 + m 2 = 0. The expectation value calculated for an event sample with fixed orientation of the reaction plane reads now w 5) and the two-mode correlation function is The probability distribution in Eq. (3.1) can then be written as a function of the (complex) Bessel-Fourier coefficients w where T as a matrix with indices (m 1 , l 1 ) and (m 2 , l 2 ). Higher n-mode correlation functions can be calculated directly from p[w], but it is convenient to derive them as n-th derivatives with respect to the source terms of the partition function This equation shows nicely that the Gaussian model for a particular centrality class needs as an input besides the background density w BG (r) only the expectation valuesw (m) l that can be determined from geometrical considerations, and the two-point correlator C In particular, the three-mode correlator takes the form (3.9) and the four-point correlation function reads The connected correlation functions can be obtained from derivatives of ln Z[j]. The connected two-mode correlator equals the connected part of eq. (3.6), and the connected correlators of more than two modes vanish of course for this Gaussian distribution.

Averaging the Gaussian model of initial fluctuations over φ R
So far, we have discussed event averages for ensembles with fixed reaction plane φ R . However, essentially all measurements are for ensembles with randomized orientation of the reaction plane. One can formally introduce a distribution for an ensemble of events with random orientation by averaging over φ R , Event averages evaluated with this azimuthally symmetric probability distribution will be denoted in the following by . . . • . It is then a consequence of azimuthal symmetry that w (m) l Similarly, the 3-mode and 4-mode correlation functions for an ensemble of randomized azimuthal orientation can be obtained from equations (3.9) and (3.10) by averaging over φ R . The result is obtained from (3.9) and (3.10) by replacing on the left hand side of these equations . . . with . . . • and by replacing on the right hand side the phases is not a Gaussian distribution even if p[w] is one. As a consequence, the connected higher-mode correlators do not vanish for the azimuthally randomized average . . . • . To illustrate this point further, we write the connected 4-mode correlator that appears in the linear response term to v m {4} 4 , are defined as the connected two-mode correlators with respect to the event average for fixed φ R , see equation (3.6), while the corresponding components of the connected two-mode correlator for an azimuthally randomized event average vanish, see (3.12). We can now make the following remarks about general properties of the probability distribution p • [w]: 1. For vanishing impact parameter, the probability distribution (3.7) becomes azimuthally symmetric even without averaging over φ R . This implies Also, azimuthal symmetry of the event-averaged geometry implies that the two-point correlation function (3.4) can depend only on φ 1 − φ 2 , and hence As a consequence, the probability distribution p • [w] is Gaussian in this limit, and all connected higher-mode correlators vanish. The distribution is fully characterized by w 2. At finite impact parameter, the two-point correlation function (3.4) of fluctuations will depend in general not only on φ 1 − φ 2 , but also on φ 1 + φ 2 . Event-averages are then still symmetric under reflections on the reaction plane, Invariance of (3.4) under this reflection symmetry implies and thus the coefficients C (m 1 ,m 2 ) l 1 ,l 2 must be real. Moreover, event averages at finite impact parameter are symmetric under the rotation φ 1,2 → φ 1,2 + π. This rotation changes each term on the right hand side of (3.4) by a phase e i(m 1 +m 2 )π . Therefore, invariance under rotation by π implies with m 1 = −m 2 would vanish, then (3.4) would depend only on φ 1 − φ 2 , but not φ 1 + φ 2 . This is not the most generic case as fluctuations will depend in general on the orientation with respect to the reaction plane and therefore they will depend on 1 2 = 0 for some m 1 + m 2 even and non-zero . in the connected four-mode correlator (3.13) can be expected to take non-vanishing values at finite impact parameter. The model discussed in section 4 provides an example for which these non-vanishing terms can be calculated explicitly, see eq. (4.19).
3. The event distribution in eccentricity m can be calculated from the probability distribution p[w]. The eccentricity is (up to a small correction to normalization) linear in the (complex) Bessel-Fourier coefficients w (m) l [46], ) Since this is a linear relation, eq. (3.7) implies that at fixed reaction plane angle φ R the E m are Gaussian distributed in the complex plane. An azimuthally randomized distribution for m is obtained by (3.20) Here, the expectation value of the eccentricity at fixed φ R is 21) and the two variance parameters are This distribution is well defined for 0 < |τ m | ≤ τ m .
4. At finite impact parameter, the remaining reflection symmetries of the event-averaged enthalpy density imply thatw (m) l = 0 for m odd, (3.24) and therefore¯ For odd m = 1, 3, 5, . . ., one can then perform the integral over φ in equation (3.20) and one finds for the distribution in eccentricities is generally non-vanishing for even and non-vanishing m 1 +m 2 . It is nevertheless interesting to investigate the simplifying ad hoc assumption that C For the event distribution (3.20) in eccentricity, this corresponds to the case τ m = 0. The integral over φ R can then be done analytically and one finds This is the "Bessel-Gaussian" distribution proposed in ref. [43,44] and used by AT-LAS to compare to distributions of event-by-event flow harmonics [59].
6. Finally, for small impact parameter b one has τ m ∼ b 2m and¯ m ∼ b m , see appendix C. For b → 0 the distribution in eq. (3.20) approaches a Gaussian distribution, In the light of these remarks, the use of the Bessel-Gaussian probability distribution p BG ( m ) in (3.27) does not seem to be the best motivated choice for the comparison to model event distributions in m (and to measured event distributions in v m ). The problem with p BG [w] is two-fold. First, the derivation of p BG ( m ) from a Gaussian distribution at fixed φ R relies on the ad hoc assumption C l 1 ,l 2 δ m 1 ,−m 2 that implies that the correlation of fluctuations is independent of their orientation with respect to the reaction plane (see discussion of equation (3.18)). Moreover, for odd m = 1, 3, 5, . . ., the Gaussian model implies¯ m = 0 (see eq. (3.25)) and this calls into question the very form of (3.27).
In fact, for odd m = 1, 3, 5, . . ., the Gaussian model at fixed φ R leads without further assumption to an explicit analytical expression for p( m ) that is of Bessel-Gaussian form but that has arguments different from p BG . We emphasize in particular that the argument of I 0 in eq. (3.26) is quadratic in m while it is linear in eq. (3.27). The two distributions can also be distinguished by the cumulants, see Eq. (3.36) and the discussion thereafter. The form of (3.26) seems better motivated, as it is derived from the general form (3.1) of the Gaussian distribution without further assumptions. For the same reason, it seems preferable to use for even m = 2, 4, 5, . . . the probability distribution (3.20) that depends on three parameter. The differences between the previously used ansatz (3.27) and the expressions derived here can be traced back to our observation (3.18) that the connected two-mode correlators C (m 1 ,m 2 ) l 1 ,l 2 do not need to vanish for even and non-zero values of m 1 + m 2 .

Linear dynamic response
If we restrict the relation (2.1) between flow coefficients V m and initial amplitudes to the linear dynamic response, , then we can determine the event-by-event distribution of flow harmonics p(v m ) by paralleling exactly the calculation of eccentricities given above, We note that we have not made the assumption v m ∼ m here. In contrast, we assume that both v m and m are given as linear combinations of the Bessel-Fourier coefficients w (m) l . This is a weaker assumption since K (m) l in eq. (3.19) and S (m)l in eq. (3.30), seen as vectors with index l, do not have to be parallel.
All the remarks made above about event-by-event distributions in eccentricity carry over trivially to p(v m ) if one restricts the discussion to linear dynamic response terms. In particular, for m even, equation (3.29) depends on the three parameters κ m , κ m andv m , and p(v m ) has the same functional form as the eccentricity distribution (3.20). For odd m, reflection symmetry impliesv m = 0, and p(v m ) reduces to a two-parameter function of the form of eq. (3.26). And at vanishing impact parameter, azimuthal symmetry implies that κ m =v m = 0, and one obtains from (3.29) a Gaussian, in complete analogy to (3.28). Finally, the Bessel-Gaussian distribution proposed in ref. [44] is obtained by assuming κ m = 0 but keepingv m finite. In fig. 2 we compare these four distributions for one set of parameters κ m , κ m andv m . In this section, we have pointed out for the first time that if We also compare to the Bessel-Gaussian distribution that results from the assumption κ = 0 at finitev m . one starts from initial fluctuations that follow a Gaussian distribution at fixed orientation of the reaction plane, then one can have a non-vanishing off-diagonal variance κ m in the distribution of p(v m ). Fig. 2 serves to illustrate that such a small non-vanishing value κ m can affect the shape of event-by-event distributions in flow harmonics.
Event-by-event distributions p(v m ) of flow harmonics were measured recently in Pb+Pb collisions at the LHC for m = 2, 3, 4 and for different centrality classes [59]. These measured distributions were also characterized in terms of their variance v 2 m − v m 2 , their mean v m , and the ratio of these quantities that takes the value for Gaussian distributions. (3.33) It was found that the distributions for m = 3 and 4 are within errors consistent with (3.33), while the distribution for m = 2 is characterized by a value significantly smaller than 4 π − 1 for non-central collisions [59]. One may wonder whether the more general form of the probability distribution (3.29) derived here can lead to an improved description of these data. While a comparison to data lies outside the scope of this work, we mention in this context that the distribution obtained from (3.29) for odd m satisfies In contrast, by setting κ m = 0 in (3.29) we find a deviation from the Gaussian result (3.33) that has the opposite sign.
+O v 6 m for a Bessel-Gaussian distribution. This implies that at least one of the two basic assumptions underlying (3.29) must be wrong: the dynamical response may not be linear and/or the distribution of initial fluctuations at fixed φ R may not be Gaussian. We shall comment in the next subsection on the first possibility, before exploring in section 4 in detail the case of universal deviations from a Gaussian distribution of fluctuations.

Non-linear dynamic response
In principle, the role of the non-linear dynamic response terms in (2.1) on the event-byevent distributions of v m can still be discussed on the level of the probability distribution p(v m ) by evaluating (3.29) with a non-linear constraint in the argument. In practice, this evaluation has to be done numerically and is likely to be involved. To gain insight into the role of non-linear dynamical response terms, we therefore turn to the study of the cumulants that characterize p(v m ). Here, we make the following remarks: 1. For vanishing impact parameter, the second order cumulant flow is dominated by the linear dynamical response In contrast, the connected 4-mode correlator (3.13) vanishes for vanishing impact parameter. This implies that also the first line on the right hand side of eq. (2.7), that gives the contribution from linear response dynamics, vanishes. In other words, v m {4} 4 depends only on terms that are proportional to some power of non-linear dynamic response terms. It is easy to see from (2.7) that these terms are non-zero in general. For vanishing impact parameter, the probability distribution is Gaussian with zero mean and the correlators on the right hand side of (2.7) that involve an odd number of modes vanish. However, there are correlators involving an even number of modes, e. g. (3.39) The six-point correlation function on the right hand side is non-vanishing for Gaussian distributions. There are 15 different contractions out of which only some vanish for symmetry reasons. There is no reason that the other contributions at order w 6 should cancel on the right hand side of eq. (2.7). From this, we conclude that if one assumes a Gaussian probability distribution of fluctuations in the limit of vanishing impact parameter, the observation of a finite value for v m {4} 4 is an unambiguous sign of non-linear dynamic response.
2. For the case of an arbitrary, non-linear dynamical response at finite impact parameter, we know that the higher-order flow cumulant v m {4} 4 will depend on physics of two different origins. First, it depends on linear dynamic response to the connected 4-point correlator w (m) •,c . This connected 4-point correlation function is non-vanishing only due to deviations of p • from a Gaussian probability distribution. The second source are terms that are proportional to powers of non-linear dynamic response terms. The latter are not expected to vanish in the limit of small impact parameter b. As one approaches more and more central collisions and if p • becomes Gaussian in that limit, one expects that the non-linear response terms start to dominate at some point. The linear contribution to v m {n} n in Eq. (3.36) vanishes for b → 0 like b n·m . In order to estimate from which centrality class the non-linear terms dominate, one would have to determine their contribution quantitatively. Finally, we remark that additional terms arise on the linear level if p • is not Gaussian. An example for this is provided in the following section.

The independent point-sources model (IPSM)
Essentially all realistic microscopic models of initial conditions incorporate the plausible assumption that the initial density distribution results from the superposition of a large number of sources that are well-localized and therefore small compared to the system size. The independent point-sources model (IPSM) realizes this assumption in a setting in which all correlation functions of initial fluctuations are analytically calculable, including their radial dependence, see below. Eccentricities have been calculated in this setting to various orders in 1/N where N is the number of sources [10,48,49,51,52]. As emphasized recently by Ollitrault and Yan [51] as well as Bzdak and Skokov [52], n-mode correlators calculated in the IPSM show characteristic deviations from a Gaussian distribution even at vanishing impact parameter. Remarkably, these deviations from a Gaussian distribution display universal properties that are shared by the class of currently explored phenomenologically relevant models of initial conditions [51]. This motivates us to explore in this section the properties of the IPSM in more detail. In particular, we shall extend the discussion of the IPSM to the case of finite impact parameter when the parametric counting of n-mode correlators will be seen to be different, we shall extend the discussion from a Gaussian to an arbitrary transverse density distribution, and we shall extend it from the characterization of eccentricity moments to arbitrary correlators of the modes w (m) l evaluated for event samples at fixed and at randomly oriented reaction plane.

1/N n−1 scaling for fixed reaction plane orientation
In the IPSM, the transverse enthalpy density w( x) of a particular event is defined as a linear superposition of contributions from N point sources, 5 is the event-averaged enthalpy per unit rapidity at initial time τ 0 . The source positions x j are random two-dimensional vectors that follow the same probability distribution p( x j ) for all j. This probability distribution is normalized,   Due to the assumption that the x are independently and identically distributed, this partition function factorizes into a product of contributions from each source, Correlation functions can now be obtained as functional derivatives of Z[j], for example (4.6) Similarly, connected correlation functions can be obtained from functional derivatives of ln Z[j], for example (4.7) Observe in particular the first term in the second line of eq. (4.7). It has the form of a contact term which is due to the point-like shape of the sources. For more realistic source shape this term decays with | x − y| on a length scale that is characteristic of its size. The second term in the second line of (4.7) can be seen as a correction to the disconnected part and is closely related to the model assumption of exactly N sources. The prefactor changes when this number is allowed to fluctuate. For the further discussion, it is useful to give also the explicit form of the 3-point correlation function The corresponding connected 3-point correlation function takes the explicit form Also higher n-mode correlation functions can be given explicitly, e.g.

Bessel-Fourier coefficients in the IPSM
So far, we have derived n-point correlation functions of w( x) in position space. Similar to our discussion of the Gaussian model in section 3, it is useful to consider the Bessel-Fourier transformation (3.2). The entire information about n-point correlation functions, that are functions of n continuous variables x j , is then encoded in the n-mode correlators w that are sets of complex-valued numbers. As shown in section 2, these n-mode correlators specify the information about initial conditions that enters flow measurements.
To write the Bessel-Fourier transform, we start from equation (4.3) for the average enthalpy density in an event ensemble with fixed orientation of the reaction plane. From this, one finds the corresponding average enthalpy density for an event sample with randomized orientation of the reaction plane, The event-averaged enthalpy densityw(r, φ) can then be written as a Bessel-Fourier expansion of the form of eq. (3.3), where the coefficientsw (m) l are determined from the orthogonality relation, With the help of the orthogonality relation (4.12), we can obtain from equation (4.7) the connected 2-mode correlator  We pass now from averages . . . for event ensembles with fixed orientation of the reaction plane to averages . . . • for ensembles with randomized orientation of φ R , . . . • ≡ 1 2π dφ R . . . . We find from (4.16) Based on these calculations, we make the following comments and observations: 1. 1/N n−1 -scaling is broken at finite impact parameter for event samples . . . • with randomized orientation of the reaction plane.
In general, this follows form the fact that the operation of passing from moments to connected n-mode correlators does not commute with the operation of randomizing φ R , that means . . . •,c = . . . c,• . 2. For small impact parameter, the b-dependence of 1/N n−1 -scaling breaking terms can be given explicitly. The 1/N n−1 -scaling is restored for central collisions when the event averages . . . • and . . . become identical. As a consequence, the terms that break 1/N n−1 -scaling must vanish in the limit of vanishing impact parameter. Remarkably, as explained in appendix B, the powerlaw dependence with which these terms vanish can be determined analytically for small impact parameter. In particular, the relevant term in the connected 2-mode correlator (4.17) scales like w Here, the dimensionless scale is set by the system size that we identify roughly with the Woods-Saxon diameter D WS . The term that breaks 1/N n−1 -scaling is then of order (b/D WS ) |m|+|m | and 1/N n−1 -scaling is effectively restored if the impact parameter is sufficiently small such that (b/D WS ) |m|+|m | becomes comparable to 1/N . For instance, for m = −m = 2 and D WS ∼ 10 fm, this is the case for b = 3 fm (assuming N ∼ O(100) which is realistic for lead-lead collisions as we shall see below). This illustrates that the expansion in small b is not only of academic interest but applies to event samples for an experimentally accessible range of impact parameter. By varying centrality and thus varying b, it is possible to move from samples that satisfy 1/N n−1 -scaling to samples in which this scaling is broken with known parametric dependence.

The non-vanishing off-diagonal variance C
The properties ofw In fact, as seen in appendix B, these symbols are simply real numbers corresponding to certain integrals over Bessel functions and independent ofw (m) l and of w BG (r). They are therefore independent of the impact parameter, and they are independent of the azimuthal orientation of the collision. In n-mode correlators, they appear multiplied with a characteristic dependence in the number of sources N .

5.
Defining the Bessel expansion with ρ(r) in terms of background density coordinates has technical advantages.
In this work, we define the Bessel-Fourier decomposition (3.2) with the help of the function ρ(r) that maps r ∈ [0, ∞] monotonously to the range [0, 1] and that we define in appendix A. This is different from previous works where we used ρ(r) = r/R, R fixed, and where we denoted the weightsw This implies in particular that at vanishing impact parameter, w • is diagonal if viewed as a matrix in l, l . It follows directly that in the IPSM the linear dynamic contribution to the second order flow cumulants take the simple explict form (4.21) 6. Higher n-mode correlators of Bessel-Fourier coefficients can be given explicitly.
While expressions for higher n-mode correlators are more lengthy, the same techniques shown here for two-mode correlators allow one to find explicit expressions for correlators of more than two modes. To illustrate this point, we give here the 3-mode sions for 2-mode correlators that are more differential than the information contained in  the r.h.s. of Fig. 4) give a sense of how model-specific dependencies become quantitatively more important for higher radial and azimuthal wave numbers that resolve finer scales.

Concluding Remarks
In general, a theory of experimentally measurable flow correlation measurements V m 1 V m 2 . . . V mn • needs to provide an understanding for the statistics of initial density perturbations in heavy ion collisions and their fluid dynamical evolution. As explained in section 2, this amounts to the requirement of knowing the initial n-mode correlators w Concerning the dynamical response functions S (m 1 ,...,mn)l 1 ,...,ln , we know that they depend only on the event-averaged azimuthally randomized enthalpy density w BG of the event class, but they do not depend on finer geometric details such as the orientation of the reaction plane, and they do not depend on event-by-event fluctuations. An explicit method of how to determine them without model assumptions was give in Ref. [30,53]. Since these dynamical response functions are (at least in principle) known, the only remaining model uncertainties in the calculation of flow correlation measurements are in determining w . To the extent to which the IPSM represents universal properties shared by all realistic models of initial conditions, this remaining model dependence is removed and model-independent predictions of fluid dynamics become possible. 8 This has motivated our detailed study of the IPSM in section 4.
In section 4, we have shown how the IPSM can be solved analytically for the full set of n-mode correlators w . This allows to compare the IPSM quantitatively to other models on a level that is more differential than an analysis of cumulants of eccentricities. A short comparison to a model with MC Glauber initial conditions in section 4.3 has shown that the IPSM shares indeed universal features with other models also on this more differential level, but that the analysis of w can also serve to delineate the azimuthal and radial length scales at which the statistics of initial perturbations in different models shows deviations from a model-independent universal behavior. We further note that essentially all other models of initial perturbations are defined in terms of computer codes that implement a physics picture. This has numerous advantages but it limits the possibilities of finding beyond a purely numerical analysis ordering principles that explain the relative importance of different contributions. Comparing models to the IPSM is hence also useful since the analytical results accessible in the IPSM allow one to find interesting ordering principles.
In particular, we have further explored in section 4 the property of 1/N n−1 scaling, that is the observation that connected n-mode correlators w  [51] that v m {2n} ∝ 1/N (2n−1)/2n which explains parametrically why measurements of higher order flow cumulants v m {2n} typically do not change within experimental errors when one increases n beyond 2. We note that for central collisions, the results of section 4 allow us to extend this ordering principle to a more general class of flow measurements. For instance, 8 This is analogous to the situation in cosmology where calculations of the cosmic microwave background and large scale structure can be expected to be model independent only to the extent to which the initial conditions are constrained by general considerations (such as symmetry arguments, based on the homogeneity and isotropy of the system) and/or observations. The question in the present context is to what extent the IPSM can serve a similar role in constraining initial conditions for the phenomenology of flow measurements in heavy ion collisions. one can show that 9 Such measurements are interesting since they depend on n-mode correlators that are not tested in the measurement of flow cumulants. We note that in these expressions, the scaling in orders of 1/N applies not only to the linear response term, but also to the contribution of the non-linear dynamical response. 10 At finite impact parameter, we have pointed out that flow correlation measurements with respect to azimuthally randomized event samples cannot be ordered in powers of 1/N . Since even the most central event class contains events with finite albeit small impact parameter, this raises the question to what extent the 1/N n−1 scaling in central events can be of practical use. Here, we have shown that deviations from 1/N n−1 scaling in the IPSM show a characteristic and analytically accessible powerlaw dependence on impact parameter. This allows one to estimate the range of impact parameter for which terms that violate 1/N n−1 scaling are sufficiently small to make 1/N n−1 scaling an applicable principle. Given that the impact parameter dependence of different linear and non-linear contributions to flow correlation measurements is different in general, one may also hope that the analytical knowledge of this b-dependence can help to disentangle different dynamical contributions. However, in the present paper, we have not yet explored this possibility further.
We close by relating some of our results to the question of why p+Pb collisions at the LHC show flow cumulants v m {2}, v m {4}, (m = 2, 3) that are comparable in size and p T -dependence to corresponding measurements in Pb+Pb collisions [56,57]. This fact has been found in fluid dynamic simulations prior to the measurements [61,62], and it is currently the focus of an important topical debate, see e.g. [63][64][65]. One question in this context is whether a hydrodynamic explanation can be regarded as being generic, or whether it reproduces data only with specific model-dependent choices. Here, we observe that in the IPSM, the parameter N can be viewed as increasing monotonously with the number of participants in the nuclear overlap. The parametric estimates for flow cumulants 9 We note that for a probability distribution characterized by its moments M , the connected n-mode correlators can be written as a sum of products of moments, Here, {Pn} denotes the list of all partitions of a set of size n, B ∈ Pn is a block in a partition Pn, and |Pn | counts the number of blocks in that partition. It follows from this structure that the specific O(1/N 5 ) cumulants defined in (5.1) have a different number of non-vanishing subtraction terms on the right hand side. 10 This can be checked for the first orders of the perturbative series in equation (2.1) by direct calculation.
In general, it follows from a theorem given in reference [60].
in pPb and PbPb read then Here, we have considered only the linear dynamic response terms that we write schematically without indicating their dependence on l. Based on these parametric estimates, one can relate the strength of the dynamic response to density fluctuations in different systems. For instance, There is phenomenological support for an almost linear relation between event multiplicity and the number of participants in a pPb or PbPb collision. Relating the number of participants approximately linearly to the parameter N in the IPSM, one can then consider different limiting cases: 1. The case N pP b N P bP b that may be realized e.g. by comparing pPb and PbPb collisions of similar multiplicity. In this case, comparable flow measurements in pPb and PbPb imply comparable fluid dynamic response S PbPb 2. The case N pPb N PbPb that may be realized e.g. by comparing central pPb to central PbPb collisions. In this case, for all initial conditions for which connected n-mode correlators of initial fluctuations scale with 1/N n−1 , the dynamic flow response S (m) must be parametrically larger for larger systems to yield harmonic flow coefficients v m that are independent of system size. Comparable values for v m {4}| pPb and v m {4}| PbPb are then consistent with the intuitive expectation that the strength of flow phenomena increases with system size. 11 In both cases, we have obtained statements about the relative parametric strength of the dynamic response coefficients S (m) in different collision systems. Note that these statements can be tested in a fluid dynamic calculations involving only minimal model assumptions. The dynamical response coefficients S (m) depend on the size of the system only via their dependence on the average background enthalpy w BG (r) but they do not carry any information about finer details of the initial transverse density distribution. In the IPSM formulated in section 4, w BG (r) and the parameter N can be chosen independently, but a more complete model of the initial state and the early dynamics will relate the number of 11 The particular parametric powerlaw dependence ∝ (N PbPb /N pPb ) 3/4 given in (5.3) was obtained by requiring parametric equality of the fourth-order flow cumulants in p+Pb and Pb+Pb. If we requires parametric equality for 6-th (8-th) order flow cumulants instead, one finds a power law ∝ (N PbPb /N pPb ) α with α = 5/6 (α = 7/8).
sources N to the size and to the radial dependence of the average enthalpy density w BG (r). Since we know how to calculate without model-dependent assumptions the dependence of S (m) on w BG (r) [30,53], and since the relation between w BG (r) and N has only a relatively mild model dependence, one can therefore test whether hydrodynamic evolution is consistent with the parametric scaling of S (m) required by equation (5.3). In our view, such a test could contribute to the important question of whether fluid dynamics can account naturally for the flow coefficients measured in systems of significantly different size, or whether some elements of fine-tuning of initial fluctuations needs to be invoked. We plan to explore this point in the near future. Here, we restrict us to formulating the question with the help of the results and insights gained in section 4. This is one illustration how the knowledge about the statistics of initial density perturbations may contribute to the further understanding of flow phenomena in nucleus-nucleus and proton-nucleus collisions.

A Background density coordinates
In this appendix we discuss a special coordinate system which can be defined for a given background enthalpy density. This coordinate system is particularly well suited for the characterization of initial fluctuations and for the numerical solution of the fluid dynamic evolution equations in the background-field formalism. We start from a transverse density distribution that is azimuthal rotation and Bjorken boost invariant, Usually w BG (r) decays rather quickly with increasing r outside of some radius which is of the order of a few fm. Also, the integrated enthalpy density per unit rapidity is finite, This maps the interval r ∈ (0, ∞) to the compact interval ρ ∈ (0, 1). For small r the relation is actually linear, ρ ∼ r and for all r the function ρ(r) is monotonous. An example for a background enthalpy distribution w BG (r) and the corresponding mapping ρ(r) is shown in Fig. 5. It is also useful to note the transformation behavior but it is easy to determine them numerically from eq. (B.1) and to tabulate them when needed.

C Impact parameter dependence
In this appendix we show that the Bessel-Fourier coefficients of the expectation value of the enthalpy density at fixed reaction plane angle φ R as in Eq. We consider a collision of two (equal size) nuclei with their centers separated by the impact parameter b. We choose the coordinate origin to be in the middle of the two nucleus centers. The expectation value for enthalpy can then only depend on the distances from the two centers, r 2 A = r 2 +b 2 /4+br cos(φ−φ R ) and r 2 B = r 2 +b 2 /4−br cos(φ−φ R ), or, equivalently on u = (r 2 A + r 2 B )/2 = r 2 + b 2 /4 and v = (r 2 A − r 2 B ) = br cos(φ − φ R ). Moreover, symmetry reasoning requires that the expectation value of enthalpyw is a symmetric function of v. One can therefore writē In the last step we have expanded in the argument v as one can do at least for small impact parameter b. One can now take the Bessel-Fourier transform of this expression. One finds that the coefficientsw (m) l have contributions only from terms on the right hand side of Eq. (C.2) with n ≥ |m|. This implies Eq. (C.1).