Characterization of initial fluctuations for the hydrodynamical description of heavy ion collisions

Event-by-event fluctuations in the initial conditions for a hydrodynamical description of heavy-ion collisions are characterized. We propose a Bessel-Fourier decomposition with respect to the azimuthal angle, the radius in the transverse plane and rapidity. This allows for a complete characterization of fluctuations in all hydrodynamical fields including energy density, pressure, fluid velocity, shear stress and bulk viscous pressure. It has the advantage that fluctuations can be ordered with respect to their wave length and that they can be propagated mode-by-mode within the hydrodynamical formalism. Event ensembles can then be characterized in terms of a functional probability distribution. For the event ensemble of a Monte Carlo Glauber model, we provide evidence that the latter is close to Gaussian form, thus allowing for a particularly simple characterization of the event distribution.


I. INTRODUCTION
In recent years, data from ultrarelativistic nucleus-nucleus collisions at the CERN Large Hadron Collider [1][2][3][4] and at BNL Relativistic Heavy Ion Collider [5][6][7] have been understood as giving strong support to a dynamical picture according to which the produced soft hadronic distributions in transverse momentum, azimuthal orientation, centrality, and particle species are determined by the fluid dynamic response to fluctuating initial conditions [8][9][10][11][12]. A detailed dynamical exploration of this picture has the potential of addressing central questions in the study of hot and dense QCD matter with nucleus-nucleus collisions. In particular, one expects that fundamental transport properties of dense QCD matter, such as the ratio of shear viscosity to entropy density [13][14][15][16], can be constrained with unprecedented accuracy from the fluid dynamic propagation of fluctuations [17,18]. Moreover, to the extent to which the fluid is almost perfect and therefore almost transparent to the propagation of fluid dynamic perturbations, fluctuation analyses may provide information about the initial conditions of ultrarelativistic nucleus-nucleus collisions and their evolution towards equilibrium [19,20]. As we shall shortly recall below, and as summarized in several recent reviews [21][22][23], a large number of recent works address this program or parts of it.
To fully exploit these physics opportunities of fluctuation analyses, one may require that a fluid dynamic formulation of ultrarelativistic nucleus-nucleus collisions should be as complete and as differential as possible with respect to the characterization of fluctuating initial conditions, their fluid dynamic propagation, and their decoupling at freeze-out. In the present work, we propose to decompose fluctuating initial conditions in a complete, orthonormal basis of fluctuating modes that can be propagated individually, mode by mode, as fluid dynamic perturbations on a smooth event-averaged background. To this end, we employ in the following a Bessel-* Stefan.Floerchinger@cern.ch † Urs.Wiedemann@cern.ch Fourier expansion that-with the exception of one remarkable work [24]-has not been explored for the characterization of initial conditions so far.
On the level of single events, this can provide, for instance, a more differential understanding of how fluctuating modes that differ, e.g., with respect to wavelength are attenuated or enhanced differently during the evolution, thus providing input to the question of whether structures arising on some spatial scales in the initial conditions can leave signatures in experimental observables or whether they will remain experimentally inaccessible because they are washed out in the course of the evolution.
On the level of event ensembles, the orthonormal basis makes it possible to determine a functional probability distribution that characterizes weights and eventwise correlations of all fluctuating modes in the initial conditions. This probability distribution can actually be evolved fluid dynamically by evolving each mode. The additional control we gain by this program can help, for example, to relate subclasses of events defined by cuts on experimental data [25] to subclasses of initial conditions, thus opening further possibilities for testing the dynamical relation in between.
The present paper is devoted to a detailed discussion of the Bessel-Fourier expansion for scalar, vector, and tensor fields, the ensuing characterization of event ensembles by probability distributions formulated in this basis, and the relation of this approach to other characterizations of initial conditions for individual events and event samples. As emphasized above, one important motivation for the choice of a Bessel-Fourier expansion is that its basis modes can be propagated individually as fluid dynamic perturbations. A detailed discussion of this fluid dynamic propagation will be left to a subsequent publication, but some first results are given already in a recent letter [26], and we comment in the following on properties that make the Bessel-Fourier expansion particularly suited for such a mode-by-mode fluid dynamic propagation of fluctuations.
By far the most common characterization of fluctuating initial conditions is in terms of a cumulant expansion of the initial (entropy) density distribution [27] that underlies the characterization of spatial azimuthal anisotropies in terms of eccentricities. We discuss in Sec. II how the coefficients in a Bessel-Fourier expansion of initial conditions are related to eccentricities. Eccentricities have been determined for initial conditions from simple model distributions [17,20,[28][29][30][31][32][33][34], as well as for full dynamical models of ultrarelativistic heavy ion collisions (such as the URQMD [35], BAMPS [36], and AMTP codes [37]). Eccentricities and closely related cumulant-based formulations have also been used to characterize angular correlations between different harmonics [38][39][40], and they play currently an important role in discussing the specific initial geometry and expected fluid dynamic response of collisions between deformed nuclei (e.g., U + U) and nonidentical nuclei (e.g., Cu + Au) and of p-Pb collisions [41][42][43][44][45]. The fluid dynamic responses that result from initial conditions with characteristic eccentricities have been studied in much detail both on the level of single events or event averages [17,29,32,35], as well as on the level of event ensembles characterized by their probability distributions [18,28,30,33,36,37,46]. By demonstrating that data on soft hadronic spectra and correlations can be reproduced in viscous relativistic fluid dynamic simulations supplemented by realistic freeze-out and by constraining the transport properties of matter, these studies have established and are now further exploiting the paradigm that heavy ion collisions produce an almost perfect fluid.
Despite the obvious use and success of a dynamic framework that relates via fluid dynamic simulations a cumulant expansion of initial conditions to hadronic observables, there are questions that one may want to address within a fluid dynamic treatment of fluctuations and for which a cumulant expansion may not provide an optimal parametrization of initial conditions. In particular, any given (positive) transverse density can, in principle, be determined fully by the infinite set of its moments or cumulants. However, given a finite set of cumulants beyond the ones that determine a Gaussian, it is not possible to find a positive transverse density corresponding to them such that higher cumulants vanish. In particular, one cannot find positive transverse density configurations that correspond to a single cumulant only, as one may want to do if one is interested in studying the propagation and attenuation of single modes. Reference [27] had understood this problem and had devised a pragmatic approach to work around it by regulating the reconstructed densities to avoid negative values. However, introducing a regulator introduces further nonzero cumulants, and therefore, in principle, one still cannot formulate initial positive transverse densities to correspond to one cumulant only. Nevertheless, this approach has been very useful for understanding how specific structures in the initial conditions propagate fluid dynamically, in particular when applied to small deviations from a Gaussian transverse density distribution [27]. However, the Bessel-Fourier expansion of initial conditions that we discuss here (see Sec. IV) may be better suited for studying the fluid dynamic propagation of fluctuations individually mode by mode, because it avoids this problem. For a dynamical treatment of individual fluctuations, it is also advantageous that this is an expansion in an orthonormal basis, while the cumulant expansion is not. Moreover, as we discuss in Sec. V, the Bessel-Fourier expansion is easily extended to the characterization of initial fluctuations in vector and tensor fields and to their fluid dynamic propagation.
To the best of our knowledge, an extension of the cumulant expansion to vector and tensor fields has not been attempted so far. While essentially all currently used models of initial conditions neglect fluctuations in the initial fluid velocity and shear viscous tensor, they seem to be a natural possibility, and we regard it as an advantage to set up a formulation that treats them on an equal footing with fluctuations in the transverse density.
We have mentioned above that it can be useful to decompose initial fluctuations in an orthonormal basis. Such formulations have been explored so far, in particular, in studies that formulate fluid dynamic perturbations on top of simple, analytically given background fields [47][48][49][50]. For these special choices of the background field, the orthogonality of the basis modes is then preserved by the fluid dynamic evolution. However, such a simplification of mode-by-mode fluid dynamics can be expected only in the presence of additional symmetries. In particular, for the case of conformal symmetry, the basis functions used in Ref. [49] do not mix in the fluid dynamic evolution. For the case of translational invariance in the transverse plane as it is realized for a background field with Bjorken flow, a two-dimensional Fourier expansion of modes has this property [50]. We note that also the orthonormal modes of the Bessel-Fourier expansion discussed here will not mix during fluid dynamic evolution if embedded as fluctuations of a Bjorken background field with transverse translational invariance. This feature may be helpful, for instance, if one plans to check the numerical accuracy of the fluid dynamic simulation of fluctuations against simple analytically known limiting cases. In general, however, we want to characterize fluctuations in all fluid dynamic fields as eventwise perturbations of smooth realistic background fields that do not share such special symmetries. For this realistic case, the modes of fluid dynamic fields will mix under time evolution, and a differential understanding of how this mixing occurs may provide additional physics insights.
The present paper is organized as follows. In Sec. II we introduce the Bessel-Fourier transform for scalar fields, we explain how the coefficients of an expansion in this basis can be determined in a CPU-inexpensive way via Lemoine's method of discrete Bessel transformation, and we explain how these coefficients are related to eccentricities. In Sec. III, we illustrate first the accuracy and use of this expansion by applying it to a simple model of fluctuating initial conditions. We then turn to the question of how event ensembles can be characterized in terms of probability distributions, and we show that the latter take a particularly simple and explicit form if expressed in the expansion coefficients of the Bessel-Fourier transform. In particular, we emphasize that a Gaussian ansatz for the probability distribution, specified in terms of event-averaged two-mode correlations only, can account with high accuracy for the event distributions in a model of initial conditions that is currently used in phenomenological studies. In this sense, realistic event ensembles of initial conditions are very well approximated by simple, analytically known expressions depending on a finite number of event-averaged input data. While the Bessel-Fourier expansion of initial densities, discussed in Sec. III, does not remain positive everywhere if truncated after a finite number of modes, we show in Sec. IV that this problem does not exist if the expansion is applied to normalized density fluctuations. Because this approach underlies our dynamical treatment of fluctuations in Ref. [26], we discuss it here in some detail. Section V is finally discussing the extension of the Bessel-Fourier expansion to vector and tensor fields. Some general properties of Bessel transformations in continuous and discrete form are given in Appendixes A and B, while Appendix C discusses some properties of functional probability distributions for event samples.

II. CHARACTERIZING FLUCTUATING INITIAL CONDITIONS
The hydrodynamical description of heavy ion collisions is normally initialized on some space-time hypersurface shortly after the collision, at the end of a regime with early nonequilibrium dynamics. In Bjorken coordinates τ , r, φ, and η, related to the laboratory coordinates t, x, y, z by t = τ cosh η, x = r cos φ, y = r sin φ, z = τ sinh η, the initialization hyper surface is usually taken to correspond to fixed τ = τ 0 . In this section, we consider the initial transverse enthalpy density w(r, φ) that characterizes the matter distribution at τ = τ 0 . To keep notation simple, we do not denote explicitly the dependence of w on time or its possible dependence on the longitudinal position along the beam direction (see Sec. V for a generalization). In practice, one might want to replace w with the initial transverse energy density , entropy density s, pressure p, or some charge density associated to a single event. Our discussion focuses first on how to characterize the fluctuating density w of single arbitrary events in terms of a Bessel-Fourier transformation, and we turn to the discussion of event averages and event distributions later.

A. Radial decomposition of w: Motivation and Lemoine's method
Our starting point is the harmonic Fourier decomposition of the azimuthal dependence of w(r, φ) in terms of the harmonics We recall that the commonly used event eccentricities n,m can be defined [51] as the normalized moduli n,m = |˜ n,m |/|˜ n,0 | of the radial moments of w (m) (r), n,m = 2π dr r n+1 w (m) (r) = |˜ n,m |e i m ψ n,m .
In recent phenomenological studies, one often focusses on one radial moment per mth harmonic, selecting, for instance, the subset of eccentricities { 2,m } or { m,m } that is then denoted by the shorthand { m }. As we see in the following, this practice may be justified to some extent by the observation that the eccentricities n,m for phenomenologically relevant density profiles tend to change only gradually and smoothly with increasing n. However, while the subset of eccentricities { m } provides an incomplete characterization of w, the set of all |˜ n,m | supplemented by the angular orientations ψ n,m is complete: The shape of the transverse density w(r, φ) of a single event can, in principle, be reconstructed unambiguously from the complete set of complex-valued n,m 's. The azimuthal decomposition (1) of w(r, φ) provides a natural ordering of the azimuthal dependence that characterizes increasingly finer azimuthal structures with increasing azimuthal wave number m. In comparison, the connection between the nth moments n,m of w (m) (r) and fluctuating modes of particular radial wavelength is arguably less direct. Here we ask how to write an alternative decomposition of the radial dependence of w(r, φ) that orders fluctuating radial modes more explicitly in terms of functions of increasingly smaller radial resolution scale. Because we expect that dynamics changes the wavelength of a fluctuation only gradually and over a sufficiently long time scale, we may hope that the modes of such an alternative expansion will mix only weakly under dynamical evolution, thus facilitating studies of the relation between modes of characteristic radial wavelength in the initial distribution and measurements that are differential in transverse momenta.
A Fourier transformation of w provides arguably the simplest decomposition of a function in terms of modes of increasing resolution scale. However, in the neighborhood of r = 0, an expansion of w in Fourier modes e ikr is not possible, because these do not satisfy the boundary condition w (m) (r) ∝ r m for small r. In radial coordinates, a two-dimensional Fourier transformation is an expansion in modes ∝e ikr cos φ and the mth harmonic moment of this Fourier mode is a Bessel function, This continuous expansion becomes discrete if restricted to a finite region r ∈ [0, R] with boundary condition w (m) (r) = 0 for r > R and all m.
One can write then with the complex coefficients w (m) l given by the integral expressions 1 To minimize notation, we distinguish here the moment w (m) (r) from its Bessel transform w (m) (k) by its argument only. No confusion should arise because most of the following discussion is in terms of the coefficients w (m)

044906-3
The discrete set of wave numbers k (m) l is defined in terms of the lth zero crossings z (m) l of the Bessel function J m (z), By construction, the expansion (6) satisfies the boundary condition (5), and the terms with increasing l correspond to modes of smaller and smaller radial resolution 1/k (m) l . In this way, the characterization of the azimuthal dependence of w(r, φ) in terms of a discrete set of azimuthal harmonics (labeled by m) can be paralleled by a characterization of its radial dependence in terms of a discrete set of radial modes labeled by l. 2 Lemoine's method [52] of discrete Bessel transformation simplifies the determination of the weights w (m) l of the Bessel expansion because it allows one to replace the integral in Eq. (7) with a finite sum, Here, w (m) l are complex-valued expansion coefficients, and the approximation (9) can be improved systematically by including a larger number of terms N l . Remarkably, according to Lemoine's method [52], the determination of the coefficients w (m) l does not involve integrations but can be done by matrix multiplication of the function w (m) (r) evaluated at a discrete set of radii The coefficients w (m) l take then the form where the matrix M (m) lα is independent of the properties of w (m) (r) and reads The value of R is a parameter in the analysis that can be chosen freely as long as Eq. (5) is satisfied. In practice, it is useful to choose R as small as possible to ensure that the expansion (9) does not need to account for regions of r in which the function w vanishes. In the numerical studies for Pb-Pb collisions, discussed in later sections, we choose Here, w BG (r) denotes an appropriately normalized "background function" that depends only on the radius r and that can be defined, e.g., as the event-averaged density w(r, φ) . The motivations for this formulation and some properties are discussed in Sec. IV. R = 8 fm. Once R is fixed, one can tabulate the matrix (12) and determine the complex expansion coefficients w (m) l . If one wants to change the number N l of terms in the expansion, both the matrix M in Eq. (12) and all the coefficients w (m) l need to be reevaluated.
From the expansion coefficients w (m) l , the spatial density distribution can then be reconstructed, Here we have made the phase dependence of the complexvalued Bessel coefficients explicit, As we illustrate with examples in the next section, the reconstructed spatial transverse density becomes an increasingly better approximation of w(r, φ) if one increases the numbers N m and N l of azimuthal and radial modes included in Eq. (13). We finally note that by inserting Eq. (9) into Eq. (3), one can express the eccentricities˜ n,m in terms of a complete set of coefficients w (m) l , where For finite N l , relation (15) is approximate and can be viewed as including all contributions to˜ n,m that result from fluctuating modes of wavelength 1/k (m) N l and larger. With increasing N l , this expression becomes more and more accurate.

III. CHARACTERIZING INITIAL CONDITIONS WITH LEMOINE's METHOD: A NUMERICAL EXAMPLE
The relations (15) and (16) illustrate that the information about w contained in the Bessel coefficients w (m) l and in the eccentricities˜ n,m is complete and mathematically equivalent. It then depends on the physics problem under consideration to decide which of these two equivalent characterizations is better suited. In particular, the relation (15) makes it explicit that for any given mth moment, the nth radial moments n,m receive, in general, contributions from fluctuating modes of various different radial wavelengths 1/k (m) l , l ∈ [1, N l ]. In contrast, the expansion (9) of the fluctuations in Bessel functions is explicitly an expansion in modes of increasing radial resolution, and it is an expansion in an orthonormal basis. This can be helpful. To illustrate the use of organizing fluctuating modes of w with the help of a discrete Bessel transformation, we turn now to an explicit numerical example.

A. A simple wounded nucleon model for the initial transverse density
For illustrative purposes, we consider a simple MC Glauber model for the initial conditions in Pb-Pb collisions, similar to the one described in Ref. [34]. In the simplest version, this model determines the enthalpy density w(r, φ) as proportional to the number of wounded nucleons. Nucleons in the incoming projectiles are distributed event by event randomly in the transverse plane according to the two-dimensional projection of a standard spherically symmetric two-parameter Woods-Saxon nuclear density profile. Nucleon-nucleon correlations in the incoming projectiles are neglected. The condition for collision between nucleons i and j of different nuclei is the simple geometric one, namely that the transverse positions (x i , y i ) of both nucleons are closer than (x i − x j ) 2 + (y i − y j ) 2 σ NN /π . Here we choose for the inelastic nucleon-nucleon cross section a value corresponding to √ s NN = 2.76 TeV, namely, σ NN = 63 mb. The transverse enthalpy density in this model is then obtained by centering at the transverse position of each participating nucleon an enthalpy contribution of Gaussian shape and with width σ B , Here the factors c i give weights to the contributions from individual participating nucleons i. In the model of Ref. [34], c i = 1. Instead, we use a MC Glauber model that determines for each participating nucleon the number of collisions that this nucleon undergoes. The prefactors c i are then chosen such that the total entropy of the system scales with , where x = 0.118. This model extension is consistent with the initial conditions used in recent fluid dynamic simulations of flow [17]. It is unimportant for the arguments made in the present paper, but because the present paper serves us also to further document the input in our recent fluid dynamic study [26], and because this study is based on state-of-the-art initial conditions, we adhere to it here. The normalization N in Eq. (17) is then fixed by the total enthalpy of the system. For numerical studies, we associate to the position of each participating nucleon a Gaussian of width σ B = 0.4 fm, except where stated otherwise. We position the center of mass (the center of enthalpy) of each event at the origin of the coordinate system.

B. Reconstructing transverse density of a single event from Bessel data
We first establish the efficiency of Lemoine's method in reconstructing the transverse enthalpy w of a single event from its Bessel coefficients w (m) l . To this end, we have chosen one particular Pb-Pb collision simulated with the wounded nucleon model of Sec. III A for impact parameter b = 0. The corresponding density distribution is shown in the top left plot in Fig. 1; it is nonvanishing in a transverse extension of radius ∼6 fm, characteristic of a Pb-Pb collision, and it shows significant fluctuations. We have determined the Bessel coefficients w (m) l of this distribution as discussed in Sec. II, and we have reconstructed the density distribution w reco(N m ,N l ) (r, φ), including a varying number of modes N m , N l in the azimuthal and radial direction. As one sees clearly from the various plots in Fig. 1, with increasing number of modes, w reco(N m ,N l ) reconstructs increasingly finer details of the transverse density of this single event. A larger number of modes is needed to resolve smaller scales. In general, we observe that for a relatively small number of expansion coefficients, the main features of the transverse density can be characterized and that the approximation of the true w in terms of the reconstructed density (13) improves rapidly in accuracy with increasing number of modes N m , N l .

C. Characterizing event averages of single fluctuating modes
As seen from Fig. 1, the MC Glauber model for the initial density distribution w generates spatial distributions with significant eventwise variations. We aim at quantifying these as fluctuations around an event-averaged density distribution. To this end, we first determine the event-averaged density distribution by averaging over the densities w i of a large number of events, i ∈ [1, N events ], Here and in the following, the brackets · · · define the event average. The event-averaged density w average of a Pb-Pb collisions at impact parameter b = 0 is shown in Fig. 2. At b = 0, this average is azimuthally symmetric. As a consequence, all eccentricities of w average with m = 0 vanish and, similarly, the event-averaged Bessel coefficients satisfy The coefficients w (0) l,average quantify the shape of w average completely. As seen from Fig. 2, w (0) l,average has significant nonzero entries only for the first few radial modes l = 1, 2, 3, 4. This reflects the fact that the shape of w average is smooth and hence only long radial wavelengths (i.e., modes with small l) are needed to characterize its radial dependence.
As a first step towards quantifying fluctuations on top of the event-averaged background distribution w average , we display in Fig. 2 l around its average, illustrates the generic fact that higher modes m correspond to increasingly finer azimuthal resolution and higher Bessel modes l correspond to increasingly finer radial resolution. Figure 2 shows also the average eccentricities n,m calculated for the same sample of 1000 central Pb-Pb events. One sees that n,m tends to increase smoothly with increasing n. In principle, if one would know precisely the n-dependence of n,m for all infinitely many n's, then one could reconstruct from this information the radial dependence of w (m) (r) analogously to the reconstruction given from the Bessel coefficients in Fig. 1. However, the radial dependence of w (m) (r) is arguably much less directly characterized by the n,m than by the w (m) l . This is so, because a mode characterized by w (m) l can be associated with a characteristic radial wavelength 1/k (m) l , while a mode n,m will, in general, receive contributions from vastly different length scales.

D. Characterizing correlations between two fluctuating modes
The Bessel coefficients w (m) l that characterize the transverse density distribution w of a single event are complex valued.
In principle, a complete characterization of event samples amounts to knowing all n-mode correlators w (m 1 ) As we argue in Sec. III E, knowledge of the two-mode correlators w (m 1 ) can provide a satisfactory characterization for practical purposes. For the model studied here, we have checked numerically that correlations between different azimuthal harmonics m 1 , m 2 vanish. 3 Therefore, we focus in the following discussion on two-mode correlators that are diagonal in the azimuthal mode m, Histograms of the event distribution of w (m) l w (m) l * are shown for m = 2 and l = 1, 2 in the top row of Fig. 3. One observes a distribution that does not peak at the event average, but that is of approximately exponential shape. We discuss this shape in the following section. . The plot illustrates also that there is a significant dispersion in phase and norm around this nonvanishing event-averaged correlation.
We have inspected the event distributions of offdiagonal two-mode products w (m) for azimuthal modes 1 m 5 and for a large number of radial modes 1 l 1 l 2 9. Some results for the event-averaged mean w (m) l 1 w (m) l 2 * are shown in Fig. 4. We observe a simple and generic pattern: For fixed m, there is a significant azimuthal correlation between fluctuations of neighboring radial resolution, l 1 , l 2 = l 1 ± 1. As the difference between radial resolutions increases a bit (next-to-neighboring modes, l 2 = l 1 ± 2), the azimuthal correlation decreases, and modes with even larger differences in radial resolution l 2 = l 1 ± n, n 3 show essentially no azimuthal correlation. We have observed the same pattern for m = 1, 4, 5 and for higher radial modes l 1 where

Gaussian probability distribution
For the case of a Gaussian distribution P[w], one has (see also Appendix C) Here, N is an appropriate normalization factor. The matrix T and the averages w (m) l in (24) are determined by For collisions at vanishing impact parameter, the azimuthal symmetry of event averages implies that At finite impact parameter, also nontrivial azimuthal modes can have nonvanishing event averages, w (m) l = 0 for m = 0. However, because event-averaged distributions do not display structures at small wavelength, one expects generically that w (m) l is non-negligible only for very small l. In general, because P is real, the matrix T in Eq. (24) is Hermitian, and a nonvanishing complex phase of off-diagonal elements (T −1 ) (m 1 )(m 2 ) between the azimuthal orientations of different modes. However, the matrix T is real and symmetric if the ensemble is symmetric with respect to the transformation ϕ → −ϕ. The statistical azimuthal rotation symmetry for b = 0 implies further that the matrix T is diagonal in the azimuthal modes m, Within the MC Glauber model we observe in Fig. 4 that twomode correlations decrease quickly with increasing difference in the radial wavelengths of the two modes, that is, and we expect that this feature is shared by other models as well.

Testing the validity of the Gaussian approximation of P
The Gaussian probability distribution (24) is fully specified in terms of the event averages w (m 1 ) , and it provides a simple ansatz for the eventwise distribution of arbitrary products w (m 1 ) Here we derive within the Gaussian approximation explicit expressions for some of these eventwise distributions, and we establish for the model of Sec. III A that these eventwise distributions are correctly described by the ansatz (24). The practical interest in this statement is that to the extent to which the Gaussian approximation (24)  In the recent letter [26], we have shown that experimentally measurable flow coefficients in nucleus-nucleus collisions can be written as the fluid dynamic response to diagonal and 044906-9 off-diagonal products of two modes, l a w (m) * l a , χ ab ≡ w (m) l a w (m) * l b , and their event averages. We are therefore particularly interested in the event distributions of ξ a and χ ab around the averages w (m) Here l a and l b denote arbitrary but fixed radial wave numbers. For an arbitrary probability distribution P [{w (m) l }], the distribution in ξ a and χ ab can be calculated from These are probability distributions in the real variable ξ a and the complex variable χ ab , respectively. For the Gaussian ansatz (24), the integrals in (30) and (31) can be done analytically. For the diagonal product ξ a , one finds (32) where the inverse width σ w is given by In this section, we focus on the case of vanishing impact parameter, for which all averages w (m=2) l vanish, and Eq. (32) reduces to an exponential P ξ [ξ a ] = σ w exp(−σ w ξ a ) (ξ a ). In Fig. 3, we demonstrate that once the event average w (m) l a w (m) * l a is specified, this provides a parameter-free and accurate description of the eventwise distribution of ξ a in the model of Sec. III A. For the event distribution in the complex variable χ ab = χ r ab + i χ i ab , we find for vanishing impact parameter (i.e., for Here, T denotes the two-dimensional submatrix obtained from T l 1 l 2 , l 1 , l 2 = l a , l b . In this way, both distributions (33) and (34) are specified completely in terms of the one-mode and two-mode correlators Eqs. (25) and (26). For the model studied here, the finite number of relevant event averages is shown in Fig. 4.
To visualize the comparison of (34) to event distributions simulated in the model of Sec. III A, we show in Fig. 5 histograms of one-dimensional projections of the off-diagonal two-mode event distributions plotted in Fig. 3. These are compared to the corresponding projection of Eq. (34). We find that the Gaussian approximation (24) accounts very satisfactorily for the shape of event distributions of off-diagonal two-point correlators χ ab , too. For the case of vanishing impact parameter, these studies indicate that, in practice, a small number of two-mode correlations is sufficient to specify fully the shape of event distributions of all products of two modes around these averages.

F. Lemoine's mode decomposition at finite impact parameter: A numerical study for b = 6 fm
So far, we have focused on heavy ion collisions at vanishing impact parameter, for which event averages are azimuthally symmetric. In this section, we show that Lemoine's method applies equally well to characterizing initial conditions at finite impact parameter. To demonstrate this, we repeat in the following the study of Secs. III C-III E for Pb-Pb collisions at impact parameter b = 6 fm. Our discussion is brief, and 044906-10 we focus only on those points that arise anew at finite impact parameter.
In Fig. 6, we show the density distribution of a randomly chosen Pb-Pb collision at b = 6 fm. In comparison to the collision at vanishing impact parameter, seen in Fig. 1, the active transverse area is clearly smaller. We have checked that Lemoine's method characterizes the simulated densities with comparable accuracy irrespective of impact parameter.
In particular, the point-by-point differences between the true enthalpy density w(x, y) and the reconstruction w reco(N m ,N l ) are less than 1% of the maximal density for a reconstruction with N m = N l = 30. Figure 7 shows some elementary characterizations of event averages at finite impact parameter when the event-averaged density w average has an approximately ellipsoidal shape that breaks azimuthal symmetry. As a consequence, there are 044906-11 Also the two-mode correlators w (m 1 ) show remarkable similarities between collisions at finite impact parameter (see Fig. 8) and collisions at vanishing impact parameter (see Fig. 3). The overall normalization of all w (m 1 ) increases with the total enthalpy of the distribution, and it is thus larger for the case b = 0. Compared to the case for b = 0, however, the relative weight of w (m=2) is significantly increased in the longest wavelength modes l 1 , l 2 = 1, 2 that characterize the event-averaged azimuthal anisotropy of w average . At finite impact parameter, there can be nonvanishing eventaveraged two-mode correlators w (m 1 ) l 1 w (m 2 ) * l 2 also for modes corresponding to different azimuthal harmonics m 1 = m 2 . In particular, the event-averaged shape of w average at finite impact parameter contributes not only to the second but also to the fourth azimuthal harmonics and this leads to nonvanishing correlations w (2) l 1 w (4) * l 2 . As seen in Fig. 9, such correlations vanish for b = 0 within statistical uncertainties, but they are found at finite impact parameter in the model studied here. The strength of these correlations is weak if compared to correlations for modes at the same azimuthal harmonics. Figure 10 shows event distributions for diagonal and nondiagonal two-mode products of the types ξ a = w (m) l a w (m) * l a and χ ab = w (m) l a w (m) * l b , respectively. In the top part of Fig. 10, we compare the simulated event distribution P ξ (ξ a ) to the analytical expectation (32) for a Gaussian probability distribution (24). This comparison is analogous to the one shown in Fig. 3 for vanishing impact parameter, but it involves now a nonvanishing expectation value w (m 1 ) l a . Figure 10 Figure 10 also shows examples of event distributions for offdiagonal two-mode correlators χ ab = w (m) l a w (m) * l b . In Fig. 11, we compare one-dimensional projections of these to analytical expectations for a Gaussian probability distribution (24). For finite impact parameter, when the expectations values w (m) l do not vanish, one finds from Eq. (31) where = 0, expression (35) reduces to the simple analytical form (34). In general, however, the event distribution P χ depends on the two event averages w (m) l a , w (m) l b , and on the three independent matrix elements T ij , i, j = l a , l b . Determining the latter from (26), we confirm also at finite impact parameter that the Gaussian approximation (24) accounts very satisfactorily for the shape of event distributions of off-diagonal products χ ab . At finite, as well as at vanishing impact parameter, a very small set of event averages (25) and (26) is therefore sufficient to specify fully the shape of all two-mode event distributions.

IV. NORMALIZED DENSITY FLUCTUATION
So far, we have discussed in Sec. II general properties of the Bessel-Fourier expansion (8) for transverse scalar densities, and we have studied in Sec. III applications of this expansion to a simple model of the transverse enthalpy density at vanishing and at finite impact parameter. The full enthalpy density w is, of course, positive everywhere and for each event. However, each fluctuating mode w (m) l J m (k (m) l r) in the Bessel-Fourier expansion will take negative values in some spatial regions. Moreover, at large radial distance r, the maximal amplitudes of the oscillating modes falls off with the root of the radial distance ∝ √ 1/r only, while the enthalpy density of each  event is expected to fall off exponentially. Therefore, after truncation at a finite number of modes, the Bessel-Fourier expansion (9) of w is not guaranteed to be positive, and it may show locally negative entries, in particular, at large r. As we have seen in Sec. III, this is not a problem for characterizing the initial conditions. However, it becomes an unwanted feature if one wants to propagate single modes w (m) l J m (k (m) l r) fluid dynamically. The propagation of locally negative densities poses certainly problems in fluid dynamics, and irrespective of whether one deals with those by "ad hoc" regularizations of locally negative contributions or in another way, the effort to ensure that the physical results depend only sufficiently weakly on a chosen prescription will be an unwanted complication. One way to bypass this problem is to seek a Fourier-Bessel expansion of the enthalpy density normalized by some conveniently chosen background enthalpy w BG (r), such that for sufficiently small fluctuations the truncated expansion remains positive by construction. Here, we explore the ansatz w (m) (r) = δ m0 w BG (r) + w BG (r) for b = 0. Same as Fig. 3, but for the normalized enthalpy density w(r, φ)/w BG (r). then to r = ∞, and artificial boundary effects would disappear. A natural mapping of this kind is induced, for instance, by the background enthalpy density w BG (r) chosen such that it agrees with an appropriate event average, w BG (r) = w(r, φ) . Because w BG (r) is positive, monotonously decreasing with r and integrable in the transverse plane (the total enthalpy is finite), one can define the transformed radial coordinate such that The coordinate ρ is proportional to r for small r and it maps the interval r ∈ [0, ∞] onto ρ ∈ [0, 1]. A reformulation of the expansion (38) in this new radial coordinate is straightforward, but we do not further explore this point in the present work.

V. VECTOR AND TENSOR FLUCTUATIONS
In this section we extend the Bessel-Fourier representation of the previous section to hydrodynamical fields that transform as vectors and tensors under rotations. We also make the dependence of these fields on spatial rapidity η and time τ explicit, which has been omitted for notational simplicity so far. Again we design the expansion such that in a situation where the background field is independent of the coordinates in the transverse plane the evolution equations for different Bessel modes decouple. We discuss here fluctuations in fluid velocity; the extension to other vector fields is then straightforward.
In coordinates τ, r, φ, η it is sensible to choose the independent components of the fluid velocity as u r , u φ , and u η . The , l 1 = l 2 . Same as Fig. 3, but for the normalized enthalpy density w(r)/w BG (r). fourth component follows from the normalization condition as u τ = 1 + (u r ) 2 + r 2 (u φ ) 2 + τ 2 (u η ) 2 . For a background fluid field that satisfies rotational symmetry in the transverse plane as well as Bjorken boost invariance, the background components u φ and u η vanish, but the radial background component u r BG can be nonvanishing. We denote the fluctuating part of the velocity fields with a tilde, and we rescale all field components such that they are dimensionless in units with c = 1. Expressed with the help of transverse Cartesian coordinates v = (v x , v y ), s = (x, y), s ⊥ = (−y, x), the fluctuating part of the velocity components takes then the form We now discuss how to set up the Bessel-Fourier representation of these fields. Naively, one might think that an expansion as in Eq. (9) would do for vector valued quantities as well. However, this leads to problems, as can be seen, for example, for the m = 0 modes. At r = 0 the Bessel functions are nonzero, J 0 (0) = 1, and an expansion ofũ r as in Eq. (9) would thus contain parts that do not vanish for r → 0. This is unphysical because the divergence of the fluid velocity would have a 1/r singularity; the zeroth harmonic moment ofũ r must vanish at r = 0. However, the first harmonic moments ofũ r andũ φ can take finite values at r = 0, while the expansion (9) does not allow for that. Instead of Eq. (9), one can expand the mth moments of the velocity fieldsũ r andũ φ as linear combinations of the Bessel functions J m−1 (k (m) l r) and J m+1 (k (m) l r) that satisfy physical boundary conditions at r = 0. Here the wave numbers k (m) l are the same as in Eq. (8). The functions J m−1 (k (m) l r) and J m+1 (k (m) l r) form an appropriate orthogonal set of functions; see Appendix A . Closer inspection shows then that physical boundary conditions are realized for the linear combinations 5 wherẽ 5 In general, the mth harmonic moments ofũ r andũ φ vanish for m = 1 at r = 0, and for m = 1 they satisfy Re[ũ (1) r (r = 0)] = Im[ũ (1) r (r = 0)] and Re[ũ (1) ]. This can also be seen from Fig. 16. One checks straightforwardly that these physical boundary conditions are satisfied by the ansatz (44).
Forũ η one can use the same expansion as in the scalar case, Note that in Eqs. (45) and (46) we also expand the dependence on rapidity η into an appropriate Fourier transform. In this sense Eqs. (45) and (46) provide generalizations of Eqs. (1) and (6).
To shortly illustrate properties of the Bessel-Fourier expansion of vector fields, we have generated initial conditions with nonvanishing velocity fluctuations according to a model described in Ref. [50]. This model supplements the MC Glauber initial conditions of Sec. III with a velocity field by associating a small random transverse velocity component v to each of the participants and their individual enthalpy density distributions. For the examples considered here, we draw the random velocity components from a Gaussian distribution of width |v| = 0.1c. Figure 15 shows the fluctuations in the radial (ũ r ) and azimuthal (ũ φ ) velocity components, generated in such a model for a single event. We mention as an aside that the initial velocity fluctuations of this model have divergent and rotational (a.k.a. vorticity) components of similar size [50]. It is an open question whether such initial velocity fluctuations leave characteristic signatures in relativistic heavy ion collisions, but at least some conceivable scenarios are being explored [54,55]. However, even if velocity fluctuations should turn out to be negligible at initial time τ 0 , they will be generated at τ > τ 0 in response to fluctuations in the enthalpy density. Understanding how the Bessel-Fourier expansion extends to vector and tensor fields is therefore relevant for studying how single-density modes propagate. In particular, our first exploratory study of the dynamical evolution of single modes of the enthalpy density [26] was based on a Bessel-Fourier expansion of all vector and tensor fields at times τ > τ 0 .
Here we do not further discuss the physics of initial velocity fluctuation, but we limit our discussion to the properties of characterizing vector fields with the ansatz (44) and (45). It is a general feature of vector fields that in the limit r → 0 their first harmonic moments can take nonvanishing values, while all other harmonic moments vanish. Figure 16 shows the harmonic momentsũ r (m) ,ũ φ (m) for the velocity fluctuations of Fig. 15 and illustrates this point.
To determine the Bessel-Fourier coefficients of the expansion (45) of vector fields, one can apply again Lemoine's method of discrete Bessel transformation. The only difference to the scalar case is now that J m (k (m) l r) in (12) gets replaced by J m−1 (k (m) l r), J m+1 (k (m) l r), respectively, so that where the matrix M ±(m) lα is independent of the properties of u ± (m) (r) and reads Let us now turn to tensor-valued fields. The prime example for this is the shear stress tensor π μν . If the event-averaged background of this shear tensor has rotational symmetry in the transverse plane and Bjorken boost invariance, then it depends only on r and the only nonzero components are π τ τ BG , π τ r BG , π rτ BG , π rr BG , π φφ BG , and π ηη BG . We note that only two of these components are independent; the other are constrained by π μν = π νμ , π μ μ = 0, u μ π μν = 0. Here, the last constraint is nonlinear in the fluid dynamical fields. It is therefore not necessarily true for expectation values. It may still be reasonable, however, to assume that the relation (49)  The reason is the modified orthogonality relation for the functions J m−2 and J m+2 in Eq. (A4). Again, these relations can be inverted with Lemoine's method. We emphasize that the expressions given in this section are of practical use. In particular, the calculation of the fluid dynamical propagation of single fluctuating modes presented in Ref. [26] involves a Bessel-Fourier decomposition of all scalar, vector, and tensor fluid dynamic fields at each time step of the simulation.

VI. SUMMARY AND OUTLOOK
In summary, we have shown in the present work that a Bessel-Fourier expansion provides a convenient orthonormal basis for the characterization of fluctuating initial conditions in all fluid dynamic fields. The form of the Bessel-Fourier expansion explored in Sec. III was proposed for scalar fields already in Ref. [24], where, in particular, results closely related to Figs. 1 and 2 of the present work were presented. Here we have extended these studies to the characterization of vector and tensor fields, we have extended it to the characterization of correlations between fluctuating modes, and we have explained how the weights of these modes can be determined in practice in a CPU-inexpensive way based on Lemoine's method. Moreover, in Sec. IV, we have introduced a variant of the Bessel-Fourier expansion for normalized densities that remains by construction positive definite if truncated after a finite number of modes. As we have argued here on general grounds, and as we have demonstrated in a first fluid dynamical study of fluctuations recently [26], this property allows one to propagate single modes fluid dynamically. The Bessel-Fourier expansion, in the form given in Sec. IV is therefore a suitable starting point for the program of mode-by-mode hydrodynamics that we plan to pursue in future work.
We have also shown that the orthonormal Bessel-Fourier expansion provides for a simple and efficient characterization of the functional event-by-event probability distribution P. To illustrate this point, we have characterized P in Secs. III and IV for the MC Glauber model of fluctuating initial conditions. We have shown for this model in particular that event distributions of single modes and distributions of products of two modes are described by a Gaussian ansatz for P with high accuracy. This is important because it allows for the discussion of event distributions in terms of simple analytic expressions that depend on a finite number of event-averaged quantities only.
For a general classification of the initial conditions of ultrarelativistic nucleus-nucleus collisions, it would be interesting to understand in the future to what extent the event probability distributions P that characterize other models of fluctuating initial conditions are also well approximated by a Gaussian ansatz. We note in this context that the framework presented in Secs. III and IV need not be limited to the analysis of event-distributions of single modes and of products of two modes. In close analogy to our discussion of Eqs. (30) and (31), one can also compare for event distributions of three or more modes the true model distributions to the results of a Gaussian ansatz. One can test, of course, whether higher-mode correlators of the form (22) factorize into products of two-point correlators, as expected for a Gaussian distribution. This would establish to what extent non-Gaussianities arise in different models of fluctuating initial conditions, and it could thus contribute to a general classification of these initial conditions. We emphasize at this point that even in a situation where h l = 0 for l > N so that Eq. (B5) is exact, this is not necessarily the case for a discrete version of the inverse relation (B7). This is in contrast to other relations of similar kind such as the discrete Fourier transforms. Equation (B6) is exact and Eq. (B7) is getting better and better as N → ∞.