Morphology of High-Multiplicity Events in Heavy Ion Collisions

We discuss opportunities that may arise from subjecting high-multiplicity events in relativistic heavy ion collisions to an analysis similar to the one used in cosmology for the study of fluctuations of the Cosmic Microwave Background (CMB). To this end, we discuss examples of how pertinent features of heavy ion collisions including global characteristics, signatures of collective flow and event-wise fluctuations are visually represented in a Mollweide projection commonly used in CMB analysis, and how they are statistically analyzed in an expansion over spherical harmonic functions. If applied to the characterization of purely azimuthal dependent phenomena such as collective flow, the expansion coefficients of spherical harmonics are seen to contain redundancies compared to the set of harmonic flow coefficients commonly used in heavy ion collisions. Our exploratory study indicates, however, that these redundancies may offer novel opportunities for a detailed characterization of those event-wise fluctuations that remain after subtraction of the dominant collective flow signatures. By construction, the proposed approach allows also for the characterization of more complex collective phenomena like higher-order flow and other sources of fluctuations, and it may be extended to the characterization of phenomena of non-collective origin such as jets.


I. INTRODUCTION
The bulk of low-momentum particles produced in heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) and at the Large Hadronic Collider (LHC) appears to emerge from a common flow field that is largely driven by the density gradients produced in the earliest stage of the collision [1][2][3][4]. This fluid dynamic paradigm of heavy ion collisions is supported by a large set of data from the LHC [5][6][7] and RHIC [8][9][10][11][12][13] which includes the dependence of particle production on the azimuthal angle relative to the reaction plane and its dependence on transverse momentum,pseudorapidity η, particle species and centrality (impact parameter) of the collision. Complementary support comes from recent indications that also event-by-event fluctuations of the energy density distribution produced in the earliest stage of a heavy ion collision, are propagated fluid dynamically. In particular, the measurement of non-vanishing odd harmonic coefficients, v n , in the azimuthal dependence of particle distributions is an unambiguous signal for the presence of sizable fluctuations in the initial stage of a heavy ion collision [14][15][16][17][18][19][20].
These findings point to a remarkable set of commonalities between the physics of the "Big Bang" determining the evolution of the Universe, and the physics of the "Little Bangs" formed in heavy ion collisions, that creates hot and dense matter under conditions akin to those prevailing in the early universe around the first microsecond, where the transition from the Quark Gluon Plasma (QGP) to a confined hadron gas (HG) occurred.
In both cases, the physical system is viewed at an ini-tial time as exhibiting a phase space distribution with a high degree of symmetry, overlaid with distributions of localized fluctuations. In both cases, these initial fluctuations are propagated fluid dynamically and are experimentally accessible, in the first case as fluctuations of the (photon) Cosmic Microwave Background (CMB) [21][22][23][24][25][26][27] and, in the second case, as fluctuations and characteristic variations in the density of particles produced in very energetic heavy ion collisions [28,29]. Also, in both cases, the decoupling of the different particle species from the common fluid dynamical system provides important possibilities for experimental verification of the collective dynamics. In the case of Big Bang physics, it sets the time scale for the decoupling of photons and determines the primordial abundances of light nuclei. In the case of heavy ion physics, it sets, i.a., the hierarchy of temperatures for chemical and kinematic freeze-out, and determines the relative abundances of hadronic resonances. Furthermore, the study of fluctuations is motivated mainly by the idea that details of the fluid dynamical evolution of a system make it possible to constrain its material properties. In the case of Big Bang physics, this has allowed to constrain the composition of the Universe in terms of visible matter, dark matter and dark energy [23]. In the case of heavy ion collisions, this has allowed to establish that the produced matter exhibits the properties of an almost perfect fluid with a shear viscosity to entropy density ratio that is close to minimal [1-4, 30, 31]. The question arises whether these qualitative analogies between the physics of the Early Universe and the physics of heavy ion collisions can provide conceptual or technical advances in either field. Can heavy ion physics learn from the techniques developed to analyze Big Bang physics, or vice versa, can it contribute to advances in modern cosmology? Synergies between both fields may arise on different levels. For instance, cosmology has studied in significant detail the question of how fluid dynamic perturbations are propagated in an expanding system, and there is at least some understanding of how non-linearities and turbulent phenomena relevant for structure formation build up, and which experimental measures may give access to them. In the context of heavy ion collisions, first analyses of the propagation of fluid dynamic perturbations have been undertaken with analytical and numerical techniques [32][33][34][35][36][37] , but these are arguably at the very beginning, and further explorations in close analogy to developments in cosmology may be beneficial.
Similarities between the Big Bang and the Little Bangs may also lead to synergies on the level of the techniques for the analysis of data on fluctuations. The main purpose of the present paper is to initiate such an effort. To this end, we shall use here a standard tool [38] for analyzing the CMB anisotropy, and apply it to the study of sets of simulated particle densities from ultra-relativistic heavy ion collisions.
We note that studies of the physical system produced in the Big Bang and the systems produced in heavy ion collisions share not only remarkable commonalities but have also significant differences. For an analysis of fluctuations and their separation from other sources of signals, some of the obvious differences (such as the vastly different scales involved in both problems, the different fundamental forces that govern the dynamics, and the difference in formulating the dynamics in general relativity or special relativity) may be less important. Potentially relevant, however, is the difference between a study of one 'big' event (even if consisting of a large but finite sample of statistically independent patches) and the study of an -in principle arbitrarily abundant -sample of mesoscopic events.
The paper is organized as follows. In section 2 we introduce the spherical harmonic approach that underlies CMB data analysis. In sections 3 and 4, we present and discuss an application of this approach to simulated heavy ion data, and in section 5 we apply the method to the particular case of elliptic flow, a standard measurement in heavy ion collisions. In section 6 we summarize our findings and we present a brief outlook. Finally, in Appendix A we have presented mathematical details of the method and generalization of the model with constant flow amplitude to more complex one, when the amplitude of the flow depends on pseudorapidity .

II. FORMULATION OF THE PROBLEM
In cosmology, the CMB anisotropy is analyzed by means of two-dimensional maps that essentially cover the full sky (like those from the COBE, WMAP and PLANCK experiments [21][22][23]), or only patches of the sky (Boomerang, MAXIMA, CBI, ACT etc. [24][25][26][27]). Analysis of such maps reveal the small temperature fluctuations and anisotropies that the early Universe has left on top of an otherwise smooth and constant background. By expressing the measured signal in terms of an expansion in spherical harmonics, useful quantities such as the power spectrum, alignment between different multipoles, statistical anisotropy and Gaussianity of the CMB can be constructed and compared with theory even from a single map (realization) of the CMB sky [21]. In the analysis of CMB data, the methods focus on the separation of the different components of the observed signal from that of the primordial CMB. The by-product of this separation yields additional maps for different kinds of foreground effects (e.g. synchrotron, free-free and dust emission) and the point-sources catalog [39].
The analysis of heavy ion data from nuclear collisions at the LHC and the RHIC faces similar challenges. On one hand, several independent statistical methods are in use [40,41] to analyze two-and multi-particle correlations in terms of flow coefficient v n , taken to be sensitive to the fluid dynamical properties of the high density partonic system in its early stages of expansion. On the other hand, the same particle correlations are known to be sensitive to other important dynamical features of heavy ion collisions, such as jets or resonance decays. Typical tasks in heavy ion collisions are then to either clean the flow signal as far as possible from confounding dynamical features, or to analyze a different class of interesting physical phenomena (such as the 'point-sources' produced by jets) without biasing its analysis by underlying collective effects. The question arises, what one could learn for this class of problems if the investigation of morphological features commonly done in the CMB data analysis was applied directly to heavy ion collisions?
Each heavy ion collision results in a distribution f (η, φ) of particles as a function of azimuthal angle φ and pseudorapidity η. Pseudorapidity is related to the polar angle θ of the produced particle with respect to the beam direction via η = − ln(tan( θ 2 )), and φ characterizes the angle in the transverse plane orthogonal to the beam. We simulate such distributions with HIJING [42] (Heavy-Ion Jet INteraction Generator), which is an event generator that reproduces a number of main features of heavy ion physics in the energy range of RHIC and LHC. It is build on PYTHIA [43] routines for hard interactions and JETSET [44] routines for string fragmentation. The initial geometry and matter distribution is determined via a Glauber model [45,46]. Various phenomenologically relevant nuclear effects are also implemented, including e.g. the nuclear modification of parton distribution functions, or parametrization of the expected parton energy loss. However, HIJING does not directly model a collective dynamics resulting in flow. As this is relevant for our study, we include flow effects a posteriori by modulating the azimuthal distribution of particles from HIJING.
Particle production in heavy ion collisions is known to display a regular η-dependence that peaks around η = 0 and exhibits a forward-backwards symmetry in the case of collisions between identical nuclei. Fig. 1 shows results for a semi-peripheral Pb-Pb collision at the LHC with a total of 12316 particles produced over more than 16 units of pseudo-rapidity. The pseudo-rapidity distribution of this randomly chosen simulated single event shows visible but relatively small fluctuations around a smooth event-averaged distribution. In preparation of the following analysis, we have replotted this distribution as a function of the polar angle θ. We see that the θ-distribution is approximately flat over an extended range π 4 < θ < 3π 4 which corresponds to a sizeable window around mid-rapidity (−0.88 < η < 0.88). A pre- sentation in the variable θ thus highlights the region of approximately one unit around mid-rapidity for which most differential data are available. As remarked already, HIJING does not generate realistic flow signals. For our study, we superimpose such flow signals a posteriori on HIJING events by redistributing the generated particles such that event-averaged ensembles of collisions show a second order harmonic oscillation of the azimuthal angle φ with some prescribed magnitude v 2 . This procedure results in events that carry azimuthal oscillations supplemented by event-by-event fluctuations. Fig. 2 shows a typical event realization for an average flow amplitude v 2 = 0.10. We caution the reader that experimental data on flow reveal characteristic dependencies of v 2 on transverse momentum, rapidity and particle identity that are not included in this ad hoc simulation. The simplified procedure adopted here will not account for all phenomenologically relevant properties of particle flow, but it will be sufficient for the purpose of our study.

III. SPHERICAL HARMONIC DECOMPOSITION AND ANALYSIS OF THE SYMMETRIES OF THE MODEL
The harmonic analysis of particle production in heavy ion collisions focus typically on the Fourier decomposition of the φ-dependence (see Fig. 2) over a narrow region in pseudorapidity around η = 0. Here, we contrast this procedure with the one employed in the CMB data analysis that starts from an expansion of the entire twodimensional map f (θ, φ) in terms of spherical harmonics where a l,m = |a l,m | exp(−imΦ l,m ) are the coefficients of decomposition with amplitudes |a l,m | and phases Φ l,m for each component l, m, P m l (cos θ) are the associated Legendre polynomials, and is the total power spectrum.
As it is seen from Eq. (1), the constant term with l = 0 which corresponds to the average density of the particles on the sphere is not included in our analysis. Thus, f (θ, φ) represents the fluctuations of the number density of particles with respect to the mean value and it can take positive or negative values. The parameter l max in Eq. (1) characterizes the angular resolution of the map. In principle, infinitely precise pointing resolution requires (l max → ∞) which is computationally prohibitive. In the CMB analysis, one choses the parameter l max according to the angular resolution of the experimental data which is dictated by the segmentation of the detector. For instance, l max = 100 corresponds to an angular resolution ∆ = π/l max ≃ 1.8 o .In heavy ion physics, other physics considerations such as the typical angular separation of nearest charged tracks in the ∆η and ∆φ may motivate the choice of l max . In the following, we use a resolution l max = 100 that is sufficiently high to reveal fluctuations on small scale. We work for illustrative purposes with the standard GLESP [38] pixelization where l max = 100 corresponds to a number of pixels in the θdirection N θ = 2l max + 1 and in the azimuthal direction N φ = 2N θ .

A. Mollweide projection of heavy ion events
To familiarize ourselves with the main features of a CMB-like harmonic analysis of heavy ion collisions, we start from a 'baseline' distribution f (θ, φ) of a heavy ion collision in which features of global collective dynamics are absent. To this end, we consider as baseline a HI-JING event without superimposed flow signal (v 2 = 0). Fig. 3 shows the corresponding fluctuations of the distribution f (θ, φ) in the so-called Mollweide projection that is heavily used in CMB analysis. In this representation the polar angle θ runs vertically from the 'north pole' to the 'south pole' and the azimuthal angle φ runs along the 'equator' from the center to left . Let us pause to discuss in more detail the information shown in this map.
We have seen already in the context of Fig. 1 that particle distributions are approximately flat in a wide window of the polar angle π 4 < θ < 3π 4 . Accordingly, the particle density and fluctuations in θ remain approximately unchanged in a broad band around the equator of the Mollweide projection (see upper part of Fig. 3). Larger particle densities are found at the poles reflecting the large number of particles emitted at small angles with respect to the beam direction (see bottom panel of Fig. 1). This enhanced particle yield at the poles is a kinematically trivial but potentially confounding factor for studies that aim at establishing event-by-event global or local signals on top of an event-average background that varies significantly in η or θ. Thus, irrespective of the choice of coordinates, the question arises to what extend potentially interesting information about fluctuations and azimuthal asymmetries can be disentangled from the strongly varying background in η or θ. The approximately uniform distribution of fluctuation in the lower panel of Fig. 3 demonstrates that removing from the expansion of f (θ, φ) the m = 0 mode of the harmonic decomposition goes a long way towards achieving this aim. We note that the m = 0 spherical harmonic components do not depend on φ but only on θ. Therefore, removing this component will remove the dominant dependence on the average pseudorapidity particle distribution without modifying any flow signal. In general, a dependence of particle distributions on pseudo-rapidity will remain after subtraction of the dominant m = 0 component, of course. In appendix A, we explain how an arbitrary η-dependence can be dealt with in the present formalism. For the sake of a particularly transparent presentation, we focus in the main text on the simpler case in which the η-dependence can be treated satisfactorily by subtracting the m = 0 mode.
To further quantify this separation of background from fluctuations, we discuss now in more detail the spherical harmonic analysis of f (θ, φ). Using Eq. (1) and, for instance, the GLESP transform method, we can decompose this distribution into its multipole components (see Appendix A for details). The result of this operation is the set of amplitude coefficients a l,m . With a suitably high choice of l max all morphological features of the initial angular distribution of counts are preserved. For the decomposition we have used l max = 100 in accordance with the original binning of the fluctuations shown in Fig. 3. The corresponding power spectrum C(l) is shown in the top panel of Fig. 4.
In the case of a perfect reflection symmetry of particle production around mid-rapidity, θ − π 2 ↔ − θ − π 2 , the odd components C S (l) of the total power spectrum vanish. Event-by-event fluctuations, however, will break this reflection symmetry. As a consequence of such fluctuations, the top panel of Fig. 4 shows very small but nonvanishing values C S (l) ∈ 10 −5 ; 10 −4 for essentially all odd integers l. In contrast, the even modes of C S (l) are for small even integers l up to a factor 1000 larger since they are dominated by the strong dependence of particle production on pseudorapidity. Removing the dominant m = 0 mode from this power spectrum, we see that all even and odd modes in C S (l) have now approximately equal signal strength, and that the odd modes do not change in strength. This quantifies the extent to which removing the m = 0 mode removes the trivial background dependence without affecting information on fluctuations. We conclude that this procedure allows one to represent a 'baseline' event, in which fluctuations are superimposed on top of a varying background, in a way in which the varying background is efficiently removed and fluctuations dominate the visual representation (see Fig. 3) and the statistical information (see Fig. 4). We shall discuss in the following section how this changes if the event contains, for example, a collective flow component.
We finally mention yet another way of checking that (3) The first part of Eq. (3) , C(l), corresponds to the most symmetric m = 0 azimuthal modes, while the second part, D(l), shows the power spectrum of fluctuation above the |m| = 0 threshold. It would be worth to note, that both these components have significantly different statistical properties, which can be characterized in terms of the probability distribution functions. The bottom panel of Fig. 4 shows the difference in the statistical properties of the signal when all m modes are included (black) and the case |m| ≥ 1 (red). Normalized to the total number of pixels ( N tot ≃ 8l 2 max ), the probability distribution function, H(s), is defined as the number of pixels with corresponding amplitude of fluctuation within an interval [s−δs, s+δs]. Here, δs = (s max −s min )/N bin , s max and s min are the absolute maxima and minima of the signal in the map with m ≥ 1 (red), and N bin = 200 is a number of intervals ( bins). As it is seen from Fig. 4 the probability distribution function for signal with removed the m = 0 mode, is quite close to a Gaussian distribution, while the one including all m modes (black) is far from a symmetric Gaussian distribution with respect to s = 0.
To the extent to which the m = 0-mode dominates the distribution, any distortion of the morphology of the map will provide relatively small corrections to the total power spectrum C S , but it may result in strong modifications of the power spectrum D(l) (see Eq. (3)), and of the corresponding multipole coefficients a l,m . Note that any features of the a l,m and of the corresponding power spectrum can be detected only if they exceed the level of statistical event-by-event fluctuations. This is why the above mentioned 'Gaussianisation' trend of the probability distribution function for the sub-dominant component power spectrum, D(l), reflects the level of detectability of particular features with small amplitudes above the statistical level.

IV. SINGLE EVENT vn-FLOW
Due to Lorentz contraction the colliding ultrarelativistic nuclei appear pancake-like just before a heavy ion collision in the laboratory frame of reference. For non-central collisions, the overlap between the two nuclei has an approximately almond-shaped cross section. This initial and anisotropic collision geometry results in large pressure gradient differences along the principal axes of the distribution and gives rise to a fluid dynamical flow characterized by fundamental properties of the fluid like the reinteraction time and residual interaction between the constituents. This flow, in turn, leads to an anisotropic distribution of particle momenta and numbers in the transverse plane which may be quantified as a function of the azimuthal angle, φ [40,41].
Recent experimental results from both RHIC and the LHC suggest that the particle flow is established at the quark-gluon level over a characteristic time scale of about 1 − 2f m/c and that the flow is quite sensitive to detailed features of the system, such as the viscosity [30,31]. Recently, the presence of higher order flow moments was understood as resulting from fluctuations around the almond shape in the initial collision geometry [14][15][16][17][18][19][20].
In the present study, the angular distribution of particles produced in each heavy ion collision will be viewed as a 'stochastic' map This expression maps the random distribution of particles without flow, f (θ, φ) ≡ d 2 N /dφdη| vn=0 onto a distribution S(θ, φ) = d 2 N/dφdη that shows azimuthal dependencies with harmonic flow amplitudes v n in azimuthal orientations Ψ n . The model distributions considered in the following were obtained in line with this expression by modulating, a posteriori, an azimuthal uniform distribution f (θ, φ) simulated in HIJING or a simple random generator with prescribed flow amplitudes v n with angular orientations Ψ n 1 . The case n = 2 corresponds to so called elliptic flow, to which we restrict ourselves in the present study. For simplicity, the model used in this section assigns the same flow value to all regions of η and transverse momentum p t , and to all particle species (e.g. mesons and baryons). This could be extended easily to account for dependencies in these parameters, or to treat the flow coefficients as random functions with statistical properties that differ from those of f (θ, φ). The arguments in the present study will not require such a more detailed modeling.
A typical experimental task is to extract from a given particle distribution d 2 N/dφdη the flow coefficients v n , defined as [47] {v n } = e in(φ−Ψn) .
Here .. denotes the average over all particles for each heavy ion collision and over the entire statistical ensemble of events. The following discussion will not involve this later step. Rather, we shall discuss how different flow coefficients v n and event planes Ψ n can be identified in a CMB-like harmonic analysis of single heavy ion events.
To this end, we decompose Eq.(4) in spherical harmonics b l,m ≃ a l,m + n v n c l,m+n e −inΨn + c l,m−n e inΨn . (6) Here, b l,m are the coefficients of the spherical harmonic decomposition for S(θ, φ), and c l,m∓n = a l,m g(l, m, ∓n), where g(l, m, ∓n) = 2πN l,m∓n N l,m This type of equation is also encountered in the CMB data analysis [48] for cases where the statistical isotropy of the signal is broken by regular modulations. We now discuss this equation in more detail. The iterative solution of Eq. (6) has an especially simple form for even-n, n = 2k, k = 1, 2... In this case, the dominant component of the a l,m -coefficients is associated with m = 0 modes, and l = even. As we have pointed out in Section 3, the strongest component of the f (θ, φ) signal is the m = 0 component (blue dots in Fig.  4). This implies that, in Eq. (6), the last term in the brackets a l,m−n is maximal, if m.
For m = 0, Eq. (6) reads b l,0 = a l,0 + n v n c l,n e −inΨn + c l,−n e inΨn = a l,0 + 2 n v n |c l,n | cos(n(φ l,n − Ψ n )) , (7) where we have used the relation a l,−m = (−1) m a * l,m . For the scenario studied here, we can neglect in Eq. (7) the second term (∝ v n |c l,n |) relative to the first one, due to the inequality |a l,0 | ≫ |c l,n | which follows from the dominance of the m = 0 modes. We consequently obtain the following approximate solution: b l,0 ≃ a l,0 .
The coefficients b l,m vanish for l < |m|. As a consequence, in Eq. coupling, but this effect is significantly smaller (∼ v 2 ), in comparison with the discussed b 4,2 and b 4,4 components. These features of Eq. (6) motivate to use b n,n = a n,n + v n b n,0 g(n)e inΨn , as a basis for separating the contributions from different v n , where g(n) = 2πN n,0 N n,n 1 −1 dxP 0 n (x)P n n (x) (see Appendix A).
For a sufficiently large flow signal, such that |a n,n | ≪ v n |b n,0 |g(n), the iterative solutions of Eq. (8) for the amplitude of the flow v n and the event plane Ψ n for a single event are then given by v n ≃ |b n,n | g(n)|b n,0 | , nΨ n = φ n,n .
B. vn-modulation with n = odd The method presented above clearly illustrates how spherical harmonics with l = even give access to amplitudes of the v n -flow coefficients with n = 2, 4, 6... We now focus on the odd harmonics v n , n = 1, 3, 5, ... that are known to provide information about the fluid dynamic response to fluctuations present in the initial stage of heavy ion collisions. To this end, we study the b l,mbasis of Eq. (6) for the case l = even and n = odd. We start from the most powerful l = 2 and l = 4 components of b l,m under consideration (see Fig. 4). For the quadrupole component, l = 2, in the presence of even and odd flow harmonics, the "resonant" component b 2,m for v 1 is b 2,1 : (10) where the phase of the event plane Ψ 1 corresponds to v 1 -flow. In Eq. (10) the major contribution to b 2,1 is given by the first term, proportional to the most powerful component a 2,0 ≃ b 2,0 , and where φ 2,1 is the phase of the b 2,1 -component, and g 1 = 4π 1 0 dxP 0 2 (x)P 1 2 (x). It is obvious that for any odd n the corresponding solution for the amplitude and the phase is given by v n = |b n+1,n | g(n + 1)|b n+1,0 | , nΨ n = φ n+1,n .
where g(n + 1) = 4π 1 0 dxP 0 n+1 (x)P n n+1 (x). Similar to n = even case, the phase of the event plane for n = odd can be reconstructed also from φ n,n phase as follows: nΨ n = φ n,n .
Note that all methods of reconstructing the orientation of the reaction plane(s) in heavy ion collisions distinguish between the 'true' reaction plane Ψ n and the orientation of the event plane reconstructed from experimental data Ψ exp n . Uncertainties in extracting the true parameter from a statistically finite amount of data are characterized by quoting the finite resolution. That means that, in respect to the statistical ensemble of realizations, v n and Ψ n should be treated as random variables with probability density functions P vn and P Ψn . Performing the same estimation as in Eq. (9) and Eq. (12) for each event we can determine the number of realizations with v n and Ψ n within a given interval of uncertainty, and in fact, we can obtain the probability distribution functions P vn and P Ψn for a given ensemble of realizations. Then, one can define the first moments of the corresponding variables v n and Ψ n , for a statistically finite amount of data. We will illustrate this approach in the next section.

V. ELLIPTIC FLOW: THE CASE n = 2
Here, we further illustrate how the symmetries of the spherical harmonic representation of the particle distribution can be used to extract the elliptic flow signal following Eq. (9). As discussed in Section 3 and shown in Fig. 2, this signal is a modulation of the random background particle distribution by a factor ∝ v n cos(2φ−Ψ). In Fig. 6 we show the corresponding power spectrum using a flow of v 2 = 0.07. The power in the components D(l), l = 2, 4, 6, 8 is seen to be significantly higher than the asymptotically flat and small D(l) ≃ const values, obtained from a random realization without flow (see Fig.  4). It is in this way that the CMB-like harmonic analysis identifies -on the basis of a single event and without recourse to an event sample -the presence of a signal on top of random localized fluctuations. The enhancement of the even low-l multipoles originates from flow contributions to b 2,2 , b 4,2 , b 6,2 . Expressing these multipoles via Eq. (6), we find where (2l + 1)D(l) = m=1 |a l,m | 2 .
We discuss now in further detail the relation between the flow signal and the multipole components determined in a CMB-like harmonic analysis. In a Mollweide projection, one can visualize easily how the |b 2,2 | component changes with a flow signal varying between v 2 = 0 to v 2 = 0.5, see Fig. 7. In Fig. 8, we demonstrate this dependence beyond the level of single events by averaging over samples of 10 5 events with a specific flow amplitude v 2 . We see that there is a linear relation between v 2 and the ratio of the dominant quadrupole components |b 2,2 | and |b 2,0 |. The small error bars shown in Fig. 8 correspond to a 1σ standard deviation. They illustrate the degree to which in the present model study fluctuations at the single event are small compared to the event-averaged mean. An analogous conclusion can be drawn from Fig. 9 about the correlation between the true event plane (Ψ R ), which is known for each simulated event, and the phase φ 2,2 determined in the analysis (see Eq. (9)). The insert shows the distribution of the difference Ψ R − 0.5φ 2,2 which is a measure of the event plane resolution obtained in this method. The distribution is consistant with a Gaussian with mean value µ ≈ 0 and standard deviation σ = 0.0254rad. The tight correlation between both orientations is due to the fact that particles at all rapidities show on average the same azimuthal correlations and this contributes to constraining statistical fluctuations in φ 2,2 .
We finally indicate how a CMB-like harmonic analysis can be used to assign to an event a probability that the measured harmonic components result as fluctuations from a random background, rather than in response to a collective flow signal. To assess this question, we have generated the distribution of |a 2,2 | values obtained by decomposing 10 5 events without flow into spherical harmonics. The probability P of assigning a given event as a random background for, say v 2 = 0.01, can be defined as the number of realizations without flow above the limit of |b 2,2 | corresponding to v 2 = 0.01. We have found 9230 out of 10 5 random realizations, and the probability of correctly concluding an elliptic flow value of v 2 = 0.01 to an event with this |b 2,2 | analysis is therefore 1−P ≃ 90%. For a five times larger signal v 2 = 0.05, the corresponding level of detectability is better than 99, 999%. We parallel this analysis for the harmonic component |b 2,1 | and have confirmed the fact that the limits for v 2 = 0.01 and v 2 = 0.05 both lie around the average of the background distribution.
The exploratory analysis presented here involves information about the simulated events that is not directly accessible experimentally. For instance, experimental data do not allow to check directly the correlation between |b 2,2 | and the 'true flow signal v 2 ', or between φ 2,2 and the 'true reaction plane Ψ R , as done in Fig. 8 and Fig. 9. In this sense, our discussion serves mainly the purpose to illustrate that a CMB-like harmonic analysis can be applied to individual heavy ion events and that it has the sensitivity of associating flow values and event planes to single events in a meaningful way. We further note that while Figs. 8 and 9 are not directly measurable, closely related quantities are. For instance, one could plot the ratio of the quadrupole components |b 2,2 | and |b 2,0 | as a function of event centrality to characterize the nature of event-by-event fluctuations in an experimentally measured event sample. We also note that we have tested the sensitivity of the method to a coarser binning of the primary particle multiplicity information, as it could result from a finite binning of the detectors. We find no significant difference in the results when using 100 bins in η and 40 bins in φ. Also, methods have been developed for the CMB analysis to handle partial detector acceptance and efficiency etc. A more detailed discussion of these aspects will be left to further work.

VI. SUMMARY AND CONCLUSIONS
The present work is a bold exploratory step aimed at testing the use of CMB-like analysis tools for the study of high-multiplicity heavy ion collisions. It parallels other recent efforts, see e.g. [28,29]. We started from representing simplified simulations of heavy-ion collisions in a Mollweide projection commonly used in the CMB analysis. We then asked how pertinent features of heavy ion collisions, including collective flow and small-scale fluctuations manifest themselves visually in such a representation, and how they are characterized by standard CMB tools.
A CMB-like analysis of fluctuations is based on an expansion in terms of spherical harmonics and thus returns the two-parameter set of coefficients b l,m . For the characterization of collective flow phenomena, that show up as purely azimuthal dependencies, the azimuthal harmonics with flow amplitudes v n provide a complete set of functions. Therefore, the larger set of coefficients b l,m must arguably contain redundancies for the characterization of collective flow in heavy ion collisions. However, event-wise fluctuations in heavy ion collisions are not confined to the azimuthal φ-direction and may be expected to show locally no preferred orientation in the η-φ-plane. Our study indicates that the larger set b l,m , while showing redundancies for the characterization of φ-dependent phenomena, may offer novel opportunities towards characterizing small-scale fluctuations in heavy ion collisions and separating them from flow effects. For instance, we observed that due to the η → −η reflection symmetry of event-averaged identical nucleus-nucleus collisions, expansion coefficients with even integers l are inevitably dominated by event-wise fluctuations. The set of these coefficients may then be used as a baseline on top of which collective or global event properties can be established. We have discussed all possible generalizations of the method in Appendix A, focusing on θ-dependency of the amplitude of the flow .This baseline allowed, for instance, to check in Fig. 4 to what extent removal of the m = 0 mode subtracts efficiently the effects of a rapidity-dependent multiplicity distribution from a fluctuation analysis. And we have seen in the discussion of Fig. 6 how the larger set of coefficients b l,m can help visually and statistically to identify flow on top of a fluctuating background.
For the detection of flow-like azimuthal modulations, the method used here is limited by the level of statistical noise v n ≥ N − 1 2 , where N is the multiplicity. Thus, for v n ∼ 10 −2 the corresponding multiplicity required is in order of 10 4 particles, which is close to the parameters at the LHC. This is parametrically the same dependence as e.g. in a flow analysis in terms of second order cumulants. Higher order cumulant expansions can improve this statistical limit. In this sense, the CMB-like analysis of heavy ion collisions shown here provides an alternative for characterizing collective flow but does not offer obvious parametric advantages. However, it offers arguably a very well-suited setting for analyzing those event-wise fluctuations that remain once the dominant azimuthal modulations have been removed. In analogy to the questions asked in cosmology, one can then characterize what sets the scale of these remaining fluctuations. For instance, which physics underlies the fluctuations seen in a Mollweide projection of heavy ion collisions after removal of flow harmonics, such as in Fig. 3. How are these fluctuations sensitive to resonance decays, jet-like particle correlations or jets? On which scales could critical phenomena leave characteristic signatures? Which similarities or characteristic differences could be expected for the corresponding maps based on fluctuations in the distribution of electric charge, or strangeness content? We intend to follow up these and further questions in the analysis of real data, as well as in further model studies. Some of these questions are well known in CMB science, when we can use , for instance, optimal filter approach to amplify the point-like sources , incorporated into diffuse background. This method in combination with wavelet transform of the signal can potentially detect localized features with very small amplitude of heavy ion collisions, which can be associated with low amplitude jets. This work is in progress and will be published in a separate paper.
where s + (θ, φ) = f (θ, φ)e −inφ , s − (θ, φ) = f (θ, φ)e inφ . (A2) The distribution S(θ, φ) can have a significant θdependence. We weigh the left and right hand side of Eq.(A1) by the window function W (θ), In radioastronomy, this is referred to as the apodization approach and it serves to amplify some particular range of the (θ, φ) map, in order to reduce some very bright sources of contamination of the signal under investigation. Here, we use it to select bins in θ for a characterization of the dependence of the signal on pseudo-rapidity. The coefficients S l,m of the spherical harmonic decomposition can be expressed as S l,m = a l,m + D l,m , where (2l ′′ + 1)(2l ′ + 1)(2l + 1) 4π and a lm =F f (θ, φ). This equation is well known in CMB physics (see for instance [49]). The major differ-ence betwen modulation of the signal in form of Eq.(A16) and v n = const is related to redistribution of the power from even multipoles to odd ones, if V 2n+1,m has non-zero components. In the case when only V 2n,m -coefficients are non-vanishing, this quadrupole modulation will redistribute the power from the n-th mode to n − 1,n + 1 modes, where l = 2n and change the orientation of the quadrupole from very planar ( v 2 = const) to non-planar, depending on amplitude of 2, 1-component. We shall explore in future work to what extent Eqs. (A16) and (A17) allow to investigate more complex modulations of the particle distribution in high-multiplicity heavy ion collisions.