Model comparison for initial density fluctuations in high energy heavy ion collisions

Four models for the initial conditions of a fluid dynamic description of high energy heavy ion collisions are analysed and compared. We study expectation values and event-by-event fluctuations in the initial transverse energy density profiles from Pb-Pb collisions. Specifically, introducing a Fourier-Bessel mode expansion for fluctuations, we determine expectation values and two-mode correlation functions of the expansion coefficients. The analytically solveable independent point-sources model is compared to an initial state model based on Glauber theory and two models based on the Color Glass Condensate framework. We also discuss to which extent general properties of initial conditions can be understood analytically.


I. INTRODUCTION
Relativistic heavy-ion collisions arguably constitute one of the most spectacular physical experiments mankind is capable of conducting in a laboratory. At collider facilities like the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC), physicists focus beams of nuclei and make them collide at relativistic energies, producing thousands of new particles per nucleus-nucleus collision. [1][2][3][4].
The incident nuclei are Lorentz-contracted discs consisting of quarks and anti-quarks as well as gluons. The total energy density is maximal at the collision time of the nuclei. At this moment, these constituents of the nucleus strongly couple to each other and form a collective medium called the quark-qluon plasma (QGP). During LHC experiments, the energy density of this extremely dense state of matter directly after the collision is twenty times as high as that of a hadron [5][6][7][8]. A droplet of QGP quickly expands and cools down. Below a critical temperature, new hadrons are formed, which is referred to as chemical freeze-out. These particles are still interacting; only once the kinetic freeze-out has occurred, they move freely [5][6][7][8].
In recent years, we understood that the QGP can be considered as a dissipative (and almost ideal) relativistic fluid, which offers the interesting possibility to study the relationship between microscopic properties described by quantum chromodynamics (QCD) and macroscopic fluid fields like densities of energy or entropy. A recent numerical scheme for solving the relativistic fluid equations of motion is described in ref. [9]. While these simulations can compute the time-evolution of fluid dynamic fields including their fluctuations, they need informations about initial field configurations as input. These initial conditions cannot be directly measured because experimentally accessible QGP properties like hadron spectra result from an integrated history of time-evolution, dur-ing which fluctuations in the initial conditions can be intensified or attenuated [10,11]. An additional complication rises from the fact that the initial conditions and in particular their symmetry properties depend on the collision centrality, i.e. whether the two nuclei collide head-on or in a peripheral manner. The underlying geometric quantity, the impact parameter, is not directly measurable either. Most importantly, initial fields are subject to quantum fluctuations, meaning that the initial conditions vary from event to event.
In principle, given an ensemble of initial field configurations, there are two strategies to calculate final state observables. One is to time-evolve each initial field event separately in so-called event-by-event simulations. This is numerically expensive and has the drawback that only one initial state model can be studied at a time. An alternative consists of splitting up the initial fields into a background field plus fluctuations around it and to evolve the fluctuations through response functions. If these fluctuations are expanded in a basis of appropriately chosen modes, it suffices to solve -just once for a given ensemble of events -the fluid dynamic equations of motion for the background field and the response functions for the perturbations around this [12]. The response functions can actually be used to compare models of the initial state that agree in the background configuration but differ in the initial conditions for perturbations.
The present paper aims at studying the statistical properties of event-by-event initial field fluctuations in the experimentally relevant ensembles of centrality classes with randomized reaction planes. Specifically, choosing a particular mode expansion for field fluctuations, we will examine expectation values and two-mode correlation functions of the corresponding expansion coefficients. The set of basis functions constitutes Fourierand Bessel modes and turns out to be advantageous for the decomposition of profiles that exhibit on eventaverage rotation symmetry around the beam axis.
As the initial fields of a heavy-ion collision are not directly measurable, we will have to rely on initialcondition models. A large class of them are based on Glauber theory [13,14]. These models consider a nucleus-nucleus collision as a superposition of indepen-dent nucleon-nucleon interactions and promote the number of nucleons participating in the collision as well as the number of binary collisions to relevant quantities. Initial event distributions are then sampled according to the position of these variables. We will make use of the initial-condition model TrENTo, which is based on a similar ansatz and proved to reproduce a large number of LHC experiments [15].
The two-mode correlation functions of Fourier-Bessel coefficients are compared to predictions of three additional initial condition models: The independent pointsources model (IPSM) of initial conditions. This particular model assumes fluctuations to originate from independent point-shaped contributions and allows closed-form expressions for the statistical quantities of interest. As we will see, while being limited in describing finer structures in position space, the IPSM still allows to qualitatively explain many properties of two-mode correlators.
As a further model we consider the Color Glass condensate (CGC) [16][17][18][19] (we concentrate on the leading approximation for large N c ), where two-point functions of energy density have been recently derived [20]. A variation of this model called Magma [21] will also be considered.
This paper is structured as follows. In section II, we put forward our procedure of characterising fluctuations in terms of Fourier and Bessel modes. In section III we discuss general statistical properties of the event ensembles we study. In section IV present the four initial condition models that are examined in the course of this paper. For each model, its main ideas are highlighted. Special emphasis is put on how the individual models allow to compute correlation functions. A detailed comparison of two-mode fluctuations in the four models is carried out in section V and we draw conclusions in section VI.

II. MODE EXPANSION FOR INITIAL FIELD FLUCTUATIONS
It is convenient to express the fluid dynamic fields of interest at the initialization time in polar coordinates because they exhibit on average rotation symmetry around the beam axis (in a coordinate system that is conveniently centred). We will express a given event profile (r, φ) in terms of a background field plus fluctuations, and decompose the fluctuation part in terms of Fourier and Bessel modes, following the approach developed in refs. [9,22,23].
The background field will be taken to correspond to an event average. For this purpose, we introduce the following function, Throughout this paper, angle brackets · denote averages over a certain ensemble of events, for instance a centrality class. (Centrality classes will be discussed in more detail in section V A.) We will assume that (r, φ) is non-negative everywhere and in a given ensemble on average normalized to unity, In other words, (r, φ) corresponds to the transverse energy density divided by the corresponding integral over the transverse plane for a given ensemble or centralilty class. This implies that the function with defines a map from the unbound interval [0, ∞) to [0, 1). The function W (r) is typically non-vanishing for small radii and decays over a length scale specific to the given centrality class, which allows a centrality-specific parametrization of the radial dependence of a given event profile in terms of ρ. Because it is always defined on the same interval, ρ is particularly suited for the Bessel expansion to be discussed below. The functions ρ(r) and W (r) define a scalar product, which can be used conveniently to construct an orthonormal basis set. One possible choice could be a polynomial basis set constructed using the Gram-Schmidt procedure [9]. Another option is to use a set of special functions that satisfy the same orthogonality properties and convenient boundary conditions, for example the cylindrical Bessel functions J m (z). For m = 0, let z (m) l denote the l'th positive zerocrossing of J m (z), the first derivative of the cylindrical Bessel function of order m. Any choice of nodes will lead to the imposition of a specific boundary condition, in particular the latter corresponds to Neumann boundary conditions [9]. For m = 0 we also include the zero crossing of J 0 (z) at the origin in the counting, i.e. z (0) are then orthogonal to each other with respect to the scalar product in eq. (5) [24], With the case distinction in our definition of z (m) l , we include with J 0 (z (0) 1 ρ) ≡ 1 a basis function with no positive zero-crossings and non-vanishing behavior for ρ → 0. This makes our set of basis functions complete.
In summary, these considerations motivate the following mode expansion, Eq. (8) defines a mode expansion for the transverse energy density. Note that one can understand m as an azimuthal wave number and similarly l as a radial wave number with larger values corresponding to more zero crossings and therefore describing finer details in space.
The inverse relation for the coefficients (m) l is given by (9) One can easily convince oneself that for real-valued fields (r, φ) ∈ R one has Note that we have concentrated here on scalar fields, but the expansion technique can also be used in slightly modified form for vector or tensor fields such as e. g. fluid velocity and shear stress [22].

III. STATISTICAL DESCRIPTION
The discussion in this section follows partly the general principles introduced in ref. [23] for a statistical characterisation of initial conditions. Let us assume that the events have been classified into centrality classes of sufficiently small extent, for example one percent. Each of these classes contains events with random orientation of the reaction plane, so that there is a statistical azimuthal rotation symmetry.
One may then characterise the transverse density in such a class by an expectation valuē where we denote by · • an expectation value for an ensemble with statistical azimuthal rotation symmetry. As a consequence of the statistical rotation symmetry, the expectation value is actually independent of the azimuthal angle φ. The function in eq. (11) is normalized according to eq. (2). (This means that the overall normalization or integrated transverse energy must also be specified to characterize a given ensemble or centrality class.) One can use¯ (r) to define the function W (r) according to eq. (1). In terms of the Bessel-Fourier expansion in eq. (8) the expectation value or one-point function is characterised by the expectation value of weights In other words, for a rotation symmetric ensemble, only the (m = 0, l = 1) coefficient has a non-vanishing expectation value, and as a consequence of the normalisation (2) it is given by 1/π. Fluctuations around the event-averaged profile can now be characterised in terms of correlation functions such as the connected two-point correlation function Through eq. (8), this can also be written in terms of the variance of Bessel-Fourier coefficients As a consequence of the statistical azimuthal rotation symmetry, the correlator in (14) is only non-vanishing when m 1 + m 2 = 0.

A. Geometry
In the following sections we will investigate three initial state models for which the one-point function is used to specify the properties of the model. Specifically, for the independent point-sources model, the one-point function determines the probability distribution of sources and for the saturation models it determines the local saturation scale. For non-central collisions, the collision geometry will have a particular role.
One can in fact describe the collision geometry in these models by introducing a non-symmetric one-point function in a first step. It describes the expectation value of the transverse density for an ensemble with fixed reaction plane angle φ R . For such an ensemble, the expectation value of the complex Bessel-Fourier weights is actually non-trivial and of the form The coefficients¯ Note that under averaging of the reaction plane angle φ R on the interval [0, 2π) with uniform distribution, the expression in (15) reduces to the one in (12). In particular, all components with m = 0 are annihilated by this operation, while the m = 0 component is unchanged. This also implies that the m = 0 components of (15) are actually given by eq. (12), i. e.¯ (0) l = δ l,1 /π.
Let us now consider two-point correlation functions. For fixed reaction plane φ R , they are of the form (16) In particular, the right hand side features not only a connected part but also a disconnected one as a consequence of non-vanishing expectation values.
If one now performs an average over the reaction plane angle φ R , one obtains for the correlation function in a rotation symmetric ensemble The first part results from the averaging of the connected correlation function. The second term on the right hand side of (17) arises from the geometry of the collision at non-vanishing impact parameter. It has non-trivial components in particular for m = 2 (and more general even m).

B. Non-linear transformations between fields
Some initial state models are formulated for the entropy density s(x) and others for the energy density (x). In order to be able to compare them, we need to do appropriate field transformations. In a close-to-equilibrium scenario, such a transformation can be done using the thermodynamic equations of state. A difficulty arises here because the relation between the different fields is in fact non-linear and as such is difficult to implement in a stochastic theory. For our present purpose, it is convenient to expand the fields around some background configuration, e. g. for entropy density s(x) =s(x) + δs(x) and similarly for energy density ε(x) =ε(x) + δε(x). Using thermal equilibrium relations in the absence of any conserved quantum numbers besides energy and momentum, one can relate the perturbations of entropy and energy density through whereT (x) is the background temperature. Using (18) one can relate the connected correlation functions (defined as in (13)), While the fields used here are physical fields, not following the normalization condition (2), it is clear that overall normalization factors can be included easily.

IV. INITIAL CONDITION MODELS
Having put forward a mode expansion to characterise initial fluid field fluctuations, we shall now turn to some currently popular models for the initial condition of heavy-ion collisions. For each model, we will recall how one can compute initial field configurations and characterize event averages and fluctuations in terms of correlations using the mode expansion introduced in section II.
We will start with the TrENTo model, which is a Monte-Carlo implementation of a generalized Glauber model. From the numerical implementation one can obtain expectation functions and arbitrary correlation functions of transverse densities in different centrality classes.
Subsequently we will discuss the independent pointsources model (IPSM), which is a semi-analytic model based on the assumption of strongly peaked (approximately point-like) sources that are distributed in the transverse plane according to a given probability distribution.
Next we will illustrate the implementation of the colorglass condensate, starting from the two-point function of the energy density obtained by [20] (we concentrate on the limit of a large N c ) and its variation called Magma [21,25,26], in which the two-point function is simplified assuming locality in position space.

A. The TrENTo initial condition model
The reduced Thickness Event-by-event Nuclear Topology (TrENTo) initial-condition model generates eventby-event initial transverse entropy density profiles, reproducing the multiplicity distributions for a wide range of LHC experiments [15]. It constitutes a Monte Carlo model that effectively interpolates between previously existing initial condition models. TrENTo describes initial field profiles in terms of two nucleus thickness functions, T A and T B . They are modeled as superpositions of Gaussians centered around previously sampled participating nucleon positions, and similarly T B ( x). The coordinates x i denote the position of participant i. The strength w

(i)
A by which a participant contributes, is sampled from a Γ-distribution with unit mean, Here, k > 0 is a continuous shape parameter regulating the fluctuations. The distribution has a long tail for k < 1 while fluctuations are suppressed for k 1. Given the fluctuating thickness functions of the two nuclei, the TrENTo model assumes the initial entropy density profile to be proportional to the generalized mean of T A and T B , with some normalization constant N . The dimensionless parameter p ∈ R controls the mixing of the two nucleus thickness functions. Note that for p = 1, one obtains the Glauber Monte Carlo model [27]. The parameter k can be tuned to match measured multiplicity distributions once p has been chosen. The code for TrENTo is publicly available [15]. We compute initial transverse entropy density profiles on a grid of 10 × 10 fm 2 with a grid spacing of 0.2 fm and the following parameter values, • reduced thickness parameter p = 0, • fluctuation parameter k = 1.4, • nucleon width σ = 0.6 fm, • overall normalization factor N = 16, • inelastic nucleon-nucleon cross section σ NN inel = 6.4 fm 2 .
The values for p, k, σ and N have been found to best fit Pb-Pb multiplicity measurements [15]. The nucleonnucleon cross section depends on the collision energy and has been chosen such that LHC energies of √ s NN = 2.76 TeV are reproduced [1]. Impact parameters are sampled from 0-20 fm. With initial field profiles at hand, we can compute twomode correlation functions of TrENTo profiles by numerically evaluating (9) and averaging over the events of a given centrality class.
While initial fields generated by the TrENTo model correspond to entropy density profiles, we convert them to energy density profiles. This allows us to compare to models that are based on energy density fields. Since the energy density scales with the power 4/3 with the entropy density (for a thermodynamic equation of state that is approximately of the ideal gas form p ∼ T 4 ), the conversion can be done by raising each TrENTo profile to the power of 4/3. Furthermore, the event distributions of a given centrality class are scaled by a common factor such that they are on average normalized to unity, cf.
(2). In addition, each event is rotated by a random angle φ R ∈ [0, 2π] in order to realize an ensemble of events with random orientation of the reaction plane.

B. Independent point-sources model
In contrast to TrENTo, the independent point-sources model (IPSM) remarkably allows to derive analytic expressions for correlation functions [23]. Let us assume that a given event profile results from N independent and identically distributed contributions whose spatial extension is small compared to the system size. We will approximate them as point-like and write for the energy density where the positions x j are all sampled from the same probability distribution p( x), normalized to unity, d 2 x p( x) = 1. The contribution w j of each point fluctuates in strength, following a probability distributionp(w) with unit mean and standard deviation σ w . In addition, we let the contribution number N fluctuate according to a distributionp(N ) with mean µ N and standard deviation σ N .
In complete analogy to ref. [23], we can derive position space correlation functions by introducing the partition sum We obtain for the one-point function simply Similarly, the two-point function reads We have introduced here the two parameters Let the probability distribution p( x) describe an event averaged field profile for fixed reaction plane angle φ R . (We will perform the averaging over φ R later on.) We can expand then along the reaction plane angle φ R similarly as in eq. (8), The coefficients¯ (m) l are real-valued and non-vanishing only for even values of m as a consequence of the two discrete symmetries p(r, φ R − φ) = p(r, φ R + φ) and p(r, φ R + φ) = p(r, φ R + φ + π). Moreover, it follows from the event normalization (2) and the definition of W (r) in (1) that¯ Using the inverse relation (9) and the orthogonality relation (6), we can obtain two-mode correlation functions from the position space two-point function, Performing the integrals, we find The numbers b (m1,...,mn) l1,...,ln with m 1 + · · · + m n = 0 constitute integrals over Bessel functions and are defined through b (m1,...,mn) l1,...,ln Let us now also introduce an ensemble average with a randomized reaction plane angle φ R using the following notation, By doing this averaging we obtain then Randomized two-mode correlators are thus only nonvanishing if m 1 + m 2 = 0. This follows directly from the statistical azimuthal symmetry around the beam axis. Furthermore, two-mode correlators in the IPSM are real-valued. From (28) we can read off that One-mode correlators of event-plane averaged events thus vanish for modes with m = 0, which is again a consequence of azimuthal rotation symmetry. Keeping in mind that¯ (0) 1 is the only non-vanishing coefficient for m = 0 in the expansion (28), we can conclude that the two-mode correlators above are equal to their respective connected two-mode correlation function, except for m 1 = m 2 = 0 and l 1 = l 2 = 1. In the latter case, one obtains The first term in (33) corresponds to the connected two-mode correlation function for a spherically symmetric spatial distribution p(r, φ) = W (r)/π, while the second term (without the part with β) accounts for geometry. Mathematically, it arises because the operations averaging over the reaction plane angle and passing from moments to connected correlation functions do not commute with each other [23]. Specifically the contribution due to geometry reads , We have now fully characterized the one-point and two-point correlation functions within the independent point-sources model. By taking higher order functional derivatives of the partition sum (24) one can also calculate higher order correlation functions when needed.

C. Color glass condensate large-Nc model
Recently, an analytic calculation of the connected twopoint function of energy momentum tensor directly after a heavy ion collision has been reported based on the Color Glass Condensate (CGC) picture [20]. The latter is essentially a model for the field theoretic description of color fields based on the paradigm of saturation. Here we concentrate on the result of ref. [20] in the large N c limit and refer to the model as the CGC large-N c model.
Let us start with the expectation value of energy density in the McLerran-Venugopalan (MV) model [20], TheQ si denote the momentum scale that characterises the colliding nuclei and 4π∂ 2 L(0 ⊥ ) is a model-dependent constant. The latter is divergent for the MV model, however the productQ 2 s (4π∂ 2 L(0 ⊥ )) is finite, since it is proportional to the square root of the one-point function. We denote by m an infrared regulator.
Note that we have been using the symbol ε (instead of ) for the energy density, since it has not yet been normalised according to eq. (2). For this purpose, we introduce the scaled field ( x) = Λε( x), with some normalization constant Λ that also converts the units from [ε] = GeV 4 to [ ] = fm −2 . In order to determine Λ consider the normalized energy density for central collisions with the background field W (r) obtained from eq. (38).
With the saturation scale Q s,0 :=Q s (0)(4π∂ 2 L(0 ⊥ )) and W 0 := W (0) at the center of the fireball, we obtain Following [21,25,28] the saturation scaleQ si of a single nucleus can be related to a thickness function, The thickness function is defined as with ρ nucl the nuclear charge density that can be approximated with a Woods-Saxon profile [13]. Also, Q si (0) is the value of the saturation scale at the center of the nuclei.
The energy density defined in this way has to be considered at fixed impact parameter and reaction angle b and is given by The averaged energy density for each centrality class can be obtained by averaging over a suitable distribution of impact parameters. For the present paper we have used the impact parameter distribution p(b) obtained from TrENTo, since there is not a canonical choice for p(b) in the CGC large-N c model itself.
The expression for the two-point function in leading order in the large N c limit as given in ref. [20] is Here, x and y are positions in the transverse plane, r = x − y denotes their difference. The two-point function depends on the local nucleons saturation scale Q s . The latter is related to the thickness function throughQ s , the strong coupling constant g, and the model-dependent function Γ(x) as in [20], with the function Γ defined in the MV model as As pointed out before,Q s is divergent, only in combination with (4π∂ 2 L(0 ⊥ )) will it lead to a finite result. Therefore it is useful to express the saturation scale directly in terms of finite quantities like the the thickness function, The ratio Γ( x− y)/4π∂ 2 L(0 ⊥ ) is regular for small radii, since the function Γ(r) has the behaviour so that the logarithmic divergence is exactly cancelled out by 4π∂ 2 L(0 ⊥ ), leading to a finite value for small r of Q 2 s ( x, y) and consequently for the two-point function (45).
The two-point functions depend on the infrared regulator m. We verified numerically that below a certain threshold of m, the two-point functions do not depend on the infrared regulator significantly any more, so that m can be set to that threshold value. In our calculation we end up with the value of m = 5 × 10 −3 GeV.
The two-point function in this model, like the onepoint function, has to be considered at fixed impact parameter and reaction angle b. The two-point function for a given centrality can be obtained through the impact parameter distribution p(b) as a weighted average over an impact parameter window [b 1 , b 2 ]. In particular the impact parameter dependence can be introduced in a similar way as for the one-point function, namely by shifting the transverse plane dependence of the saturation scale by a b/2 term. This corresponds to the following replacement in eq. (45) The dependence of this expression on the reaction plane angle is given by the direction of b, so φ R is fixed in eq. (45). As we discussed in section III A, connected two-mode correlation functions in an ensemble with randomized reaction plane angle have two contributions, one corresponding to an azimuthal average of the correlation function, and another from the product of one-point functions.
Two-mode correlation functions can be obtained using (29) with (45) as position space two-point function and two additional integrations to account for the average over the impact parameter distribution p c (b) of a given ensemble and the reaction plane angle. It is very convenient to make use of some properties of the CGC large-N c position space two-point function (45) to reduce the number of integrals and hence the computation time. Indeed, the two-point function depends on the position variables only through as well as Here we parametrized the impact parameter as b = b(cos(φ R ), sin(φ R )) to make its two independent parameters explicit. The two-point function thus takes the form Hence, introducing additional integrations over the impact parameter and the reaction plance angle as well as performing the following change of integration variables, expression (29) simplifies to The Kronecker-delta arises from the φ 1 -integration. Just as for the IPSM, two-mode correlators thus vanish except for m 1 +m 2 = 0. Furthermore, G is upon integration over φ R an even function inφ (the relative minus sign betweeñ φ and φ R in the fifth argument of G can be handled with a change of integration variables φ R → −φ R ), which implies that two-mode correlators are real-valued. We computed two-mode correlation functions for the CGC large-N c model by numerically evaluating (55) with a non-symmetric background field from (38) and adding the contribution from geometry (37).

D. Magma model
Recently, a simplification of the CGC large-N c model has been proposed. Dubbed Magma, it considers the energy density as a superposition of localised interactions falling off like 1/r 2 [21]. The connected position space two-point correlation function of energy density in this model can be written as with the local function In (57) an infrared cutoff regulator m has been introduced. It should be chosen such that the following scale hierarchy holds where R denotes the nuclear radius. The saturation scale is proportional to the thickness function, cf. (42), and its impact parameter dependence is implemented bȳ The rescaled position-space two-point function is then given by with where we defined r ± = | x ± b/2|. Remarkably, the strong coupling constant g drops out. Following [21], we will set m = 0.14 GeV and Q s,0 = 1.24 GeV. At this point it is worth working out some general properties of two-point models of the form with a model-dependent function f . Note that for f ( x) ∝ p( x) this reproduces the connected part of the IPSM position space correlation function in eq. (26). The function f ( x) changes when one transforms fields according to eq. (19). An expression for two-mode correlation functions can be obtained from (62) using (29), The function f can in general depend on the azimuthal angle φ in a non-trivial way, therefore it is natural to consider its Fourier expansion. Introducing f m (r) = 1/(2π) dφ e −imφ f (r, φ), we obtain If the model function f depends on the impact parameter vector as well, one would have to include an additional integral db p c (b) as well as an averaging integral over the reaction plane angle. In this case, however, thanks to the delta-distribution in (62), the reaction plane angle appears only as φ − φ R . The φ R -average then assures that B m11,m2;m l1,l2 vanishes for m = 0. Hence, the twomode correlation functions vanish for m 1 +m 2 = 0. They are also real-valued because B (m1,m2;m) l1,l2 constitutes an integral over real-valued functions if m = 0.
We computed two-mode correlation functions for the Magma model by numerically evaluating (65) with a non-symmetric background field from eq.(38) and, just like for the CGC large-N c model, by adding the contribution from geometry (37).

V. COMPARISON OF INITIAL FIELD MODELS
With four initial state models at hand, we shall now compare their predicted two-mode correlations. We will first specify how we categorised TrENTo events into centrality classes. Since the other model do not have an independent way of defining a centrality class but rather relie on the distribution of the impact parameter as an external input, we adopt the same distribution as obtained in TrENTo and the corresponding centrality class definition. Next we will compare the four different models presented before in terms of their one-point functions. Subsequently, we will focus on two-mode correlators in the four models and compare them.

A. Centrality classes
The total initial transverse entropy is to a good approximation proportional to the final charged-particle multiplicity per unit rapidity [29]. This allows us to categorise TrENTo events in centrality classes according to multiplicity, the 0-1% class containing those 1% of all events with the highest multiplicity, 1-2% referring to the succeeding 1% and so on. The multiplicity of an event is highest in central collisions and decreases for increasing impact parameters. We will mainly constrain our analysis to two centrality classes: the rather central 0-1% class as well as the 20-21% class, for which we expect a nonvanishing impact parameter to play a role. Numerical data for these, as well as other centrality classes will be made available as ancillary files to this article.
Defining centrality by means of the multiplicity is useful since it is directly accessible through experiments. Another possibility is to define centrality through the impact parameter. This is the natural choice for the two CGC models, as they are parametrized by b. It is not obvious, and typically not true, that the two centrality definitions lead to identical ensembles. However, in the following, we will assume ensembles to be sufficiently similar for a comparison between the different models to be valid.
In order to define centrality classes in terms of the impact parameter, we take the impact parameter distribution p(b) obtained from TrENTo (cf. Fig. 1) and use its percentiles as impact parameter window [b min , b max ] of a given centrality class. The impact parameter distribution p c (b) for a centrality class of width 1%, which enters into (55), is then given by A more sophisticated scheme would be to determine the impact parameter distribution for each (multiplicity based) centrality class with TrENTo.

B. Expectation value of energy density
The TrENTo model and the CGC models differ in their mean energy density definition for a given centrality class, therefore the W -function (defined in eq. (1)) entering in the Bessel-Fourier expansion in eq. (8) is different for the two classes of models, as well. Recall that the TrENTo entropy density (which determines its energy density) is defined in (22) to be proportional to a generalized p-average of the reduced thickness functions. In contrast, the energy density in the CGC models is given by (44) as the product of the saturation scales Q 2 of two nuclei, the position dependence of which is taken proportional to the reduced thickness functions of the corresponding nuclei. The Magma model and the CGC share the same energy density definition. In the IPSM model the energy density can be considered as an input, therefore we choose to take TrENTo one-point functions.
In Fig. 2 we show the resulting background field profiles W (r) as function of the radius in the two classes of models, for the 0-1% and the 20-21% classes. The TrENTo profile is slightly broader than the CGC energy density on average, which is due to different definitions of energy density in terms of thickness function eq. 22 with p = 0 in TrENTo and eq. 44 for the CGC models.

C. One-point functions
The IPSM necessitates the coefficients¯ (m) l in order to predict two-mode correlations (cf. eq. (33)). They also appear in the geometry term (37), which has to be added to the CGC large-N c and Magma model to compute twomode functions with randomized reaction angles. As has been discussed in section III A, the coefficients¯ In Fig. 3 we present¯ (m) l computed from TrENTo events as a function of l and m for two centrality classes. For both centrality classes, expectation values with m odd vanish, as follows from symmetry considerations. In addition, for m = 0, only¯ (0) 1 = 1/π ≈ 0.32 is nonvanishing, as we have concluded before.
Note also that the moduli of the non-vanishing expectation values decay approximately exponentially with increasing values of l, which highlights a posteriori that our chosen set of basis functions is well suited to describe initial field profiles using a low number of expansion coefficients. We note that the decay as a function of l is much more quick in the 0-1% class than in the 20-21% class. This is only natural as the coefficients¯ (m) l quantify azimuthal variations of the background field, cf. (28). The background field, on the other hand, is symmetric for central collisions and ever more elliptically deformed for higher centrality classes.
In Fig. 4 we compare the coefficients¯   the two models. This can be seen as a consequence of the different definitions of the energy density in terms of the reduced thickness functions.
The coefficients¯ (m) l as defined in (15) and displayed in Fig. 3 can fully characterize the contribution to harmonic flow coefficients from the collision geometry. To characterize contributions from fluctuations we need also two-point functions, to which we turn next.

D. Two-mode correlation functions
We shall now turn to our key quantity of interest, namely connected two-mode correlation functions of ensembles with randomized reaction plane angles. For this discussion we introduce the following notation for the connected two-mode correlation function, and for the diagonal part It follows from the statistical azimuthal rotation symmetry that the two-mode correlation functions are realvalued and vanish except for m 1 + m 2 = 0. We can therefore constrain our analysis to real-valued correlators of the form G (m,−m) l1,l2 . We have explained in the previous section for each individual model how to retrieve two-mode correlators. However, we still need to fix the constants α and β in the IPSM for comparisons to the other models to be possible. This was carried out by fitting the IPSM expression (33) to five TrENTo correlators, namely those G (m) 1 with 0 ≤ m ≤ 4. There are two reasons for this specific choice. Firstly, we demanded correlators on the diagonal as we can expect them to be non-vanishing regardless of centrality. And secondly, we made sure to choose correlators with a small index l. The reason is that we cannot expect the IPSM with its point-shaped contributions to remain valid as we pass to finer details in position space.
a function of the radial wave numbers l 1 and l 2 for the two centrality classes 0-1% and 20-21% and comparing the four models introduced before. Note that this representation is only possible because the correlators are real-valued and depend for fixed m only on the two indices l 1 , l 2 .
In the 0-1% class, all four models show approximately diagonal two-mode correlation functions. While the offdiagonal correlators are almost rigorously zero in the IPSM, they are more pronounced in TrENTo, especially for higher values of l 1 , l 2 .
Both the IPSM and the Magma model have significantly larger diagonal values G (m) l than the TrENTo and CGC large-N c model. This is actually the reason to use two different color schemes in Fig. 5. Moreover, one observes here already that the diagonal values G (m) l decay with l in the TrENTo and CGC large-N c model, while they actually increase for the IPSM and Magma models. This is directly linked to the assumption of point-like sources. We will further comment on this below.
The off-diagonal elements in the two CGC models arerelative to the corresponding diagonal elements -smaller than they are in TrENTo. In a direct comparison of the two CGC models, the correlators of the Magma model are in this sense more diagonal than those of the CGC large-N c model.
In the IPSM, strict diagonality is a direct consequence of centrality: The coefficients¯ (m) l vanish for (m, l) = (0, 1) in central events so that only the term propor-tional to δ l1,l2 in expression (33) survives. In constrast to the IPSM, the off-diagonal elements in TrENTo are likely to result from the fact that the assumption of pointshaped sources is dropped in favour of an extended Gaussian shape. This fits with the observation that the offdiagonal elements become more pronounced for higher radial wavenumbers l, l , which probe finer structures in position space.
When passing from the 0-1% class (Fig. 5) to the 20-21% class (Fig. 6), one can equally make out deviations from the initially diagonal structure, most notably for those entries with either l 1 = 1 or l 2 = 1 of color plots with m = 2 and (less prominently) m = 4. These off-diagonal contributions, however, result from non-centrality of the collisions, expressed through the geometry term in eq. (37).
A further aspect that all color plots share, concerns the sign of the diagonal elements. Indeed we have sign G This a direct consequence of eq. (10) which implies The plots establish that two-mode correlation functions are largest on the diagonal. On the other hand, the correlator G (0) 1 is vanishing in all models. It corresponds to the variance of (0) 1 , which quantifies fluctuations of the integrated field, i.e. total transverse energy. Fluctuations of this quantity are naturally suppressed for narrow centrality classes, as we consider them here.
In Fig. 7 we compare the variation of the correlators G (m) l on the diagonal as a function of l in the four models and for different values of m. As we established before, the sign of the data series follows a (−1) m -pattern. Regardless of centrality, the diagonal elements in TrENTo and the CGC large-N c model seem to converge towards zero for large radial wave numbers l after peaking at around l ≈ 3. For m = 2, the geometry part covers the peak in the 20-21% class, making the curve decrease monotonously. In contrast, the correlators of the IPSM and Magma diverge linearly with l and with similar slopes. This is likely a consequence of the point-like approximation for correlations underlying these two models. It is well possible that the higher l modes are actually efficiently damped by a viscous fluid evolution and that the increasing behaviour for large l is therefore not visible in final state observables. This will be investigated in further work. Interestingly, the TrENTo correlators agree fairly well with those of the CGC large-N c model, in both centrality classes. For small values of l, the IPSM agrees with the previous two models as well, only from around l = 3 on can an increasing disparity be observed. Similarly, Magma data matche the other models up to l = 3 for m = 0 and l = 1 for m = 4, while diverging away for higher values of l.
The variation of the IPSM correlators shows that the assumption of point-shaped sources is fairly valid as long as one probes coarse structures in position space. However, for radial modes with large wave numbers l, the finer structure of correlations in position space is unveiled. Hence, the IPSM begins to differ from models working with finite source extensions. The close similarity between the IPSM and Magma is natural, because the latter model can be seen a generalization of the IPSM but still shares with it a Dirac-distribution-shaped contact term.
Finally in Fig. 8 we show how the diagonal parts of the two-mode correlators G (m) l (defined in eq. (68)) depend on centrality. The absolute value of the presented correlators seems to be monotonously increasing as a function of centrality. This is a natural consequence of geometry contributing more significantly with rising centrality to the correlators.
As for the agreement of the four models with each other, this can be investigated with respect to the mode numbers m and l as well as centrality. Fig. 8 suggests that the four models lead to similar results for central collisions and show increasing discrepancy as a function of centrality.
Generally, for low values of l and m, the four models lead to similar results, the discrepancy depending only slightly on centrality. Instead, the models diverge faster away from each as a function of centrality when increasing values of m and l are chosen. The m-dependence of this phenomenon seems weaker than its l-dependence. In addition, the TrENTo data matches up well with the large-N c model, and the IPSM with Magma, as we discussed before. Equally, the l-dependence supports our previous discussion of the limits of the IPSM resulting from the point-shaped nature of its sources. We note also TrENTo and IPSM show in all panels of Fig. 8 a maximum around the 80% centrality.
Interestingly the geometry contribution of the CGC large-N c model dominates the correlation function (dotted line respect to green line in Fig. 8), meanwhile the corresponding contribution in TrENTo (grey line) is much smaller than the value of the two point function (orange line). Fig. 8 thus suggests that the four models agree fairly well with each other for sufficiently low radial wave numbers l and centrality classes. In this regime, the analytically solvable independent point source model (IPSM) would then exhibit universal properties shared with other initial state models. A detailed investigation of our results shows that this universality is restricted to not too large values of the radial wave number l and more pronounced for the more central multiplicity classes. This is interesting because one expects that higher l modes are anyway damped more strongly by dissipative effects in the fluid regime.

VI. CONCLUSION AND OUTLOOK
We have studied here models for the event-by-event fluctuations in the initial energy density distributions of Pb-Pb collisions at √ s NN = 2.76 TeV in ensembles with randomized reaction plane angles. For this purpose, we used a complete Fourier-Bessel mode expansion scheme for fluctuations around a background profile specific to a given centrality class. The statistical properties of different models for the initial state can then be characterized through one-mode and two-mode correlation functions depending on azimuthal wave numbers m and radial wave numbers l.
Comparing four different initial state models that are currently discussed in the literature, we have computed two-mode correlation functions of initial energy density. While TrENTo and the two Color Glass Condensate (CGC) models demanded a numerical evaluation, the independent point-sources model (IPSM) allowed to find analytical expressions for correlators. Presenting the results exemplary for the centrality classes 0-1% and 20-21%, we were able to capture characteristics of fluctua- tions for central collisions as well as effects of the collision geometry for less central collisions.
Indeed, the IPSM qualitatively agrees with three significantly more extended initial state models, as nonvanishing two-mode correlation functions are of approximately diagonal form with off-diagonal corrections in non-central collisions. For small radial wave numbers l, the IPSM agrees even quantitatively with the other models, while differences could be explained with the pointshaped nature of sources in the IPSM for larger l.
An important outcome of this study are actually the concrete values of one-point functions and two-mode correlators for the different models. We have obtained them in numerical form and will make them available to the public as ancillary files accompanying this article. In a future publication we plan to use these initial data together with the FluiduM framework [9,12,22,23,[30][31][32][33]. This will allow to calculate from them two-particle correlation functions and harmonic flow coefficients that can be measured experimentally. We are curious to see whether any of the four models discussed here is favoured by experimental data. For a first step one could actually use the fluid parameters (overall normalization of entropy density, initialisation time, viscosities and freeze-out temperature) obtained in ref. [33] by fitting to transverse momentum spectra of identified particles. No additional parameters would be needed and one could directly see which initial state model works best.
An interesting possibility also arises from the observation that the two-mode correlation functions of all four investigated models agree reasonably well for small radial wave numbers l and for central collisions, indicating there a form of universality. It is expected that modes with higher values of l are damped more strongly by shear and bulk viscous dissipation. Also they lead to oscillating patterns on the freeze-out surface and one can therefore expect that they have a weaker contribution to final state two-particle correlation functions. The observed universality for central collisions suggests to concentrate on those in order to constrain thermodynamic and transport properties of the quark-gluon plasma, while more peripheral centrality classes might be more suitable to distinguish between different initial state models (see also [34][35][36]).
Before closing, we would like to mention possible ways to continue this line of research. One could extend the