The Angular Power Spectrum of Heavy Ion Collisions

TL;DR

This paper analyzes the angular power spectrum of particle multiplicities in heavy ion collisions to understand early state anisotropies and extract flow coefficients, providing new insights into Quark Gluon Plasma dynamics.

Contribution

It introduces the use of angular power spectrum analysis to study fluctuations and flow coefficients in heavy ion collision data, a novel approach in this context.

Findings

Odd modes follow a power-law with index close to one

Angular power spectrum accurately extracts flow coefficients

Results based on Pb-Pb collision data at 2.76 TeV

Abstract

Particles produced in heavy ion collisions carry information about anisotropies present already in the early state of the system and play a crucial role in understanding the Quark Gluon Plasma and its evolution. We explore the angular power spectrum of particle multiplicities in such heavy ion collisions to extract fluctuations in particle multiplicities on the surface a sphere. Results are presented for Pb-Pb data at $\sqrt{s_{NN}} = 2.76\mathrm{TeV}$, extracted from the ALICE open data portal. We find that odd modes of the power spectrum display a power-law behavior with corresponding index $\beta$, which is found to be close to unity. We also demonstrate that the angular power spectrum allows us to extract accurately the flow coefficients of non-central collisions.

Figures

Equations

$$f(\theta ,\varphi )=\sum _{l=0}^{l_{max}}\sum _{m=-l}^{m=l}{a}_{lm}{Y}_l^m(\theta ,\varphi ),$$

$$Y_l^m(\theta ,\varphi )=\sqrt{\frac{2l+1}{4\pi }\frac{(l-m)!}{(l+m)!}}{P}_l^m(\cos (\theta ))e^{im\varphi },$$

$$a_{lm}=\frac{4\pi }{N_{pix}}\sum _{j=0}^{N_{pix}-1}{Y}_{lm}^{\ast }({\theta }_j,{\varphi }_j)\bar{f}({\theta }_j,{\varphi }_j),$$

$$C_l=\frac{1}{2l+1}\sum _{m=-l}^{m=l}|a_{lm}{|}^2.$$

$$C_l^{m\ne 0}=\frac{1}{2l+1}\sum _{m=-l}^{m=l}|a_{lm}{|}^2-\frac{|a_{l0}{|}^2}{2l+1}.$$

$$\frac{dN}{d\varphi }\propto \frac{1}{2\pi }\left[1+2\sum _{n=1}^{\infty }{v}_n\cos (n(\varphi -{\psi }_n))\right],$$

$$a_{l0}$$

$$a_{lm}$$

$$b_{lm}=\sqrt{\frac{2l+1}{4\pi }\frac{(l-m)!}{(l+m)!}}{\int }_{{\theta }_i}^{{\theta }_f}\sin \theta P_{lm}(\cos \theta )d\theta$$

$$|v_n{|}^2=\frac{(2n+1)}{2}\cdot \frac{C_n^{m\ne 0}}{|b_{nn}{|}^2}\cdot \frac{|b_{00}{|}^2}{4\pi }.$$

Full text

††preprint: APS/123-QED

I Introduction

Heavy ion collisions at ultrarelativistic energies, such as currently studied at the top energy ranges of the Large Hadron Collider (LHC) at CERN, can typically generate more than 20,000 charged particles in each central collision Adam et al. (2017). In the first instants of the collision a strongly interacting Quark Gluon Plasma (QGP) is formed Arsene et al. (2005); Back et al. (2005); Adams et al. (2005); Adcox et al. (2005) which subsequently expands as an isolated system, cools, and hadronizes. A central current theme is the study of the collective expansion of the QGP, which can be described by viscous hydrodynamics and which provides information on the viscosity and other transport properties of this exotic state of matter. A powerful analysis method relies on the understanding of particle production and momentum resultant from colliding ions which hit each other with intermediate impact parameters. The azimuthal distribution of the charged particles created in the overlap (participant) zone can be characterized by flow coefficients $v_{n}$ of the harmonic modes contributing to the Fourier decomposition of the distribution Adam et al. (2016).

In ref. Naselsky et al. (2012) it was proposed to analyze large-multiplicity heavy ion collisions by means of statistical tools developed for the study of the Cosmic Microwave Background in cosmology. The idea is obvious: to view the QGP as the opaque very early universe, and to see the hadronization as emanating from the “last scattering surface” of the QGP. From that stage the color-singlet hadrons can escape and eventually, perhaps after residual fragmentations, reach the particle detector which we parametrize as in the Mollweide map of the sky. In the case of the Cosmic Microwave Background the event is unique and it is naturally analyzed by expanding the observed temperature distributions in terms of spherical harmonics, viewing the distribution as a function of two angular variables (polar angle $\theta$ and azimuthal angle $\phi$). Just as the associated angular power spectrum of the Cosmic Microwave Background contains information on the expansion of the Universe, we may hope glean information about the system of quark and gluon matter, and its expansion and eventual hadronization, from its angular power spectrum.

In this Letter we present such a study of the angular power spectrum based on the analysis of public heavy-ion data from Pb-Pb collisions at $\sqrt{s_{NN}}=2.76\mathrm{TeV}$ collected with the ALICE experiment and taken from the CERN open data portal ALICE Collaboration ; Bourjau (18). Earlier works in similar directions have been reported in  Naselsky et al. (2012); Sarwar et al. (2018) and  Llanes-Estrada and Muñoz Martinez (2018). The present study takes a different approach by basing the analysis entirely on the angular power spectrum of the collisions. Special attention needs to be paid to the fact that the heavy ion sky is cut by detector limitations in the direction of the polar caps, and we show how to overcome this by subtraction of the $m=0$ mode. For odd $l$-modes we observe a power-law behavior and we define the corresponding critical index $\beta$, a new and potentially interesting observable in heavy ion collisions. Finally, we compare our determinations of flow coefficients $v_{n}$ with results yielded from more conventional anisotropic flow analyses, finding good agreement and also some intriguing small discrepancies that suggest that a further focus on these methods may provide new insight into the flow patterns observed.

II Map and Spectrum

Charged particles produced in ultra relativistic heavy ion collisions at the LHC can be described in terms of their azimuthal angle $\phi$ and polar angle $\theta$ between the particle’s 3-momentum and the beam axis. Therefore, multiplicity distributions are expressed by a function $f(\theta,\phi)$, which can be expanded in spherical harmonics,

[TABLE]

with

[TABLE]

where $\theta\in[0,\pi]$, $\phi\in[0,2\pi)$, and $P_{l}^{m}(\cos(\theta))$ are the associated Legendre Polynomials. The $l$-modes are labeled as multipole moments, with $l=0$ being the monopole, $l=1$ the dipole, and so on. Here it is assumed that modes with $l>l_{max}$ hold insignificant signal power.

The software HEALPix (Hierarchical Equal Area isoLatitude Pixelation) Górski et al. (2005) partitions a spherical surface in pixels of equal area, providing a map of an event’s final particle distribution as a pixelation of Eq. (1). The number of divisions on the sphere depends on the chosen resolution. Each particle with coordinates ($\theta_{j},\phi_{j}$) maps to a pixel $j$ on the sphere, as shown in Fig. 1. The color coding follows particle density times map resolution (number of pixels).

The polar angle relates to pseudorapidity $\eta$ via $\theta=2\arctan(\exp(-\eta))$. The limitation in pseudorapidity to $|\eta|<0.9$ is imposed on all charged particles reconstructed with the Time Projection Chamber (TPC) of the ALICE detector  Selyuzhenkov and the ALICE Collaboration (2011). At the 0-5% centrality, the multiplicity of each event ranges from $\sim 2000$ to $\sim 3000$ and a little experimentation leads us to choose $l_{max}=23$. Through HEALPix, the sphere is then divided into 768 pixels of equal area and the region within $|\eta|<0.9$ ($44^{\mathrm{o}}\lesssim\theta\lesssim 136^{\mathrm{o}}$) for which our data are limited contains 544 pixels. To better visualize where the particles cluster, the maps in Fig. 1 have been smoothed by $5^{\mathrm{o}}$: particles were bundled in gaussian beams with full width at half maximum of $\sim 0.09$ radians ($5^{\mathrm{o}}$).

In order to correct for detector efficiency, we have chosen the following approach: first, all events corresponding to a certain centrality were added and thus mapped onto the same Mollweide projection (top of Fig. 1) regardless of the orientation of their reaction plane. This map is referred to as $F^{all}(\theta,\phi)$. Due to the random azimuthal orientation of each collision’s reaction plane, all remaining anisotropies along $\phi$ in $F^{all}(\theta,\phi)$ may then be attributed to detector defects. Second, each event map was divided by the all-event map within the $|\eta|<0.9$ range. Thus, given that $f(\theta,\phi)$ represents an event map, $f(\theta,\phi)\to\bar{f}(\theta,\phi)=f(\theta,\phi)/F^{all}(\theta,\phi)$, with $\bar{f}(\theta,\phi)$ the corrected map. The latter is depicted in Fig. 1 (bottom) and we use it to determine the coefficients $a_{lm}$:

[TABLE]

where $N_{pix}$ represents the total number of pixels. Additionally, it should be remarked that this is the zeroth order estimator, as pixelating $f(\theta,\phi)$ ($f(\theta,\phi)\to f(\theta_{j},\phi_{j})$) corresponds to taking its average within each pixel with surface area $\Omega_{pix}$ Górski et al. (2005). The higher order estimators are implemented in the facilities of HEALPix. From Eq. (3) the angular power spectrum $C_{l}$ of an individual heavy ion event is then defined by

[TABLE]

The multipole moments relate to the angular scale $\alpha$ of the distribution through $\alpha=\frac{180^{\mathrm{o}}}{l}$.

The fixed value $C_{0}=4\pi$ for all events is a consequence of the chosen normalization, obtained by dividing $\bar{f}(\theta,\phi)$ by the event multiplicity and multiplying it by $N_{pix}$. In the present case, $\bar{f}(\theta,\phi)$ is a piecewise function: while well defined within $44^{\mathrm{o}}\lesssim\theta\lesssim 136^{\mathrm{o}}$, it takes a null value otherwise. This causes a suppression of $l=4n$, $n=1,2,\ldots$ relative to the other even modes, as seen in Fig. 2, an effect first identified in ref. Llanes-Estrada and Muñoz Martinez (2018).

In order to test these important edge effects on the power spectrum, we have generated 8000 distributions that are isotropic on the surface of a unit sphere and have the same multiplicities as the 0-5% centrality. We then computed $C_{l}$ for each event after submitting $(\theta_{j},\phi_{j})$, $j=0,...,N_{pix}$ to the cut $44^{\mathrm{o}}\lesssim\theta\lesssim 136^{\mathrm{o}}$. The resulting averaged power spectrum is shown in Fig. 3 and compared to actual data. For a sphere with smooth uniform distribution $Y_{00}=1/\sqrt{4\pi}$ holds all the power. Consequently, $C_{0}=4\pi$ and $C_{l}=0$ for $l>0$. This remains approximately true for a discrete distribution that is isotropic on average, which we label ‘full sky’: the monopole yields a power spectrum value of $4\pi$, while the other moments have $C_{l}\sim 10^{-3}$ (see Fig. 4), only differing from zero due to the chosen finite multiplicity - the $\langle C_{l}\rangle$ in Fig. 4 come from events with multiplicity of $\sim 15000$. As the latter increases, $C_{l}$ approaches zero for $l>0$; this issue shall be addressed in more detail in ref. M. Machado (2019).

It is striking that the isotropic distributions shown in Fig. 3, which have the same multiplicities as the real events from 0-5% centrality, have almost identical averaged power spectra. The modes $l=4n$, $n=1,2,...$ are strongly suppressed relative to the other even $l$, demonstrating that such feature originates from the limited detector acceptance. Tiny differences can be observed: $C_{4,8,12}$ have slightly higher values for data than the isotropic distributions. This suggests that there are more fluctuations within a solid angle $\Omega=\pi/2$ for real data than for the isotropic case.

Since gross features of the angular power spectrum are so well reproduced by simulations, it is straightforward to conclude that the suppression of modes with $l=4n$ is an artifact of data being limited to the $\theta$-range shown. This is demonstrated in Fig. 4, where we show averaged power spectra of isotropic distributions with same multiplicity for different cuts in $\eta$. As the range in $\eta$ narrows, even $l$-modes become more enhanced.

For isotropic distributions, $\langle C_{l}\rangle$ for odd $l$-values are consistently low, as expected on account of parity symmetry and as seen in Fig. 4. Real data, even those associated with central collisions, show a markedly different behavior as shown in Fig. 5 for the case of 0-5% centrality. Because the data appear to follow a power law we have performed fits of averaged power spectra to the function $C_{l}=A\cdot l^{-\beta}+\mathcal{C}$. The exponent $\beta$ appears to be independent of centrality and it is very close to unity (a best fit to a constant yields a value $~{}1.068$, consistent with unity within $1\sigma$) as shown in Fig. 6. The standard deviation of each $\beta$ was calculated from the variance of the parameter estimate.

III Power Spectrum for $m\neq 0$

The limitation imposed by the acceptance of the TPC leads us to explore a new strategy that corrects the data from artificial suppressions of those even-$l$ modes that are just a consequence of the geometric limitation. Since we are free to define the angular power spectrum we need not tie ourselves to the standard definition employed in situations where there is an essentially uniform map of the full sky, with smaller fluctuations on top. We therefore define a modified angular power spectrum where the ’global’ $m=0$ mode is excluded:

[TABLE]

This definition is designed to remove detector acceptance effects while keeping essential physical information of all modes except those corresponding to $m=0$. Indeed, from Fig. 7 we see that the averaged power spectrum for the isotropic distributions is close to trivial, as expected. On simulated data there are no microscopic physical mechanisms to introduce a non-trivial signal and the power spectrum is essentially flat. On the other hand, real data display a clear peak in $l=2$, aside from also the enhanced $\langle C^{m\neq 0}_{1}\rangle$ and $\langle C^{m\neq 0}_{3}\rangle$ values. The same pattern occurs for all centralities (Fig. 8), suggesting the presence of real anisotropies in the considered events. The continuous enhancement in magnitude as the spectra become associated with more peripheral collisions is mainly due to decreasing multiplicities, since map sparsity leads to an increase in fluctuations, thus enhancing $C_{l}$ values.

IV Flow Extraction

Matter produced in heavy ion collisions exhibits strikingly what is known as collective flow Gustafsson et al. (1984); Appelshäuser et al. (1998), indicating that anisotropies present in the early stages of the collision result in increased particle emission in certain directions. In the presence of transverse flow (perpendicular to the collision axis), the final particle distribution along the azimuthal direction can be decomposed as Voloshin and Zhang (1996); Poskanzer and Voloshin (1998):

[TABLE]

where $v_{n}$ are the flow coefficients and $\psi_{n}$ the symmetry planes corresponding to the different $n$ modes. Coefficients $v_{1}$, $v_{2}$ and $v_{3}$ are denoted directed, elliptic and triangular flow, respectively.

Given our observed enhancement of low-$l$ modes in the averaged power spectra $\langle C^{m\neq 0}_{l}\rangle$ as compared to the isotropic case (see Fig. 7), it is natural to ask if this could be caused by flow and thus potentially provide an alternative method for determining it. To investigate this question, we associate $\bar{f}(\theta,\phi)$ with the left hand side of Eq. 6 for $44^{\mathrm{o}}\lesssim\theta\lesssim 136^{\mathrm{o}}$ ($|\eta|<0.9$) and zero otherwise, while constant in $\theta$. From the expansion in spherical harmonics (Eq. 1) and the analytical form of Eq. 3 we find

[TABLE]

where

[TABLE]

for all $l$ and $m$. Here $\theta_{i}$, $\theta_{f}$ define the initial and final values of the interval where $\bar{f}(\theta,\phi)$ is nonvanishing. Given that $\sin{\theta}\cdot P_{lm}(\cos{\theta})$ is symmetric around $\pi/2$, the only surviving coefficients will have $l$ and $m$ sharing the same parity. Additionally, flow coefficients $v_{n}$ as well as their angles of symmetry planes $\psi_{n}$ do not contribute to any of the $m=0$ modes. The event plane angle $\psi_{n}$ does not affect the angular power spectrum since $a_{lm}\propto e^{-in\psi_{n}}$ and only $|a_{lm}|^{2}$ enters there.

Solving for $|v_{n}|$ $n=1,2$, after combining Eq.(5) with Eqs.(7, 8), we find

[TABLE]

For the case at hand, $|b_{31}|^{2}$ is of $\mathcal{O}(10^{2})$ smaller than $|b_{31}|^{2}$, leading to the use of Eq.(9) in the calculation of $v_{3}$.

Increased sparsity of particles leads to higher $C_{l}$ values due to larger fluctuations and it is crucial to take this into account when studying their magnitude. We propose to correct for the trivial increase due to limited multiplicity by subtracting the spectrum of isotropic distributions with multiplicities corresponding to the chosen centralities: $\langle C^{m\neq 0}_{n}\rangle\to\langle C^{m\neq 0}_{n}\rangle-\langle C^{m\neq 0[iso]}_{n}\rangle$ in Eq. (9). To compare, we have generated $10^{6}$ isotropic Mollweide maps limited to $|\eta|<0.9$ with corresponding multiplicities. As in all previous cases, the angular power spectrum has been extracted from each map and averaged over all events.

The extraction of leading flow coefficients $v_{1}$, $v_{2}$ and $v_{3}$ from the power spectrum then proceeds as follows: For each centrality we compute $\langle C^{m\neq 0}_{n}\rangle-\langle C^{m\neq 0[iso]}_{n}\rangle$ and substitute it into Eq. (9). We have successfully tested the high accuracy of this method with Monte Carlo data that took specific $v_{i}$’s as input. Using this method we have calculated $v_{n}$ for the ALICE open data ALICE Collaboration and compared with the Q-cumulants method of flow analysis for two-particle correlations Bilandzic et al. (2011) applied to the same data set. In Fig. 9 these coefficients are labelled $v_{n}\{C_{l}\}$ and $v_{n}\{2,QC\}$, respectively. For directed and elliptic flows ($v_{1}$ and $v_{2}$), numbers agree better than for triangular flow coefficients $v_{3}$, which show small deviations, especially for peripheral centralities. We stress that when both methods are applied to Monte Carlo simulations, i.e., purely flow scenarios, they agree uniformly.

The angular power spectrum is simply the two-point correlation function in Fourier space. In that case, computing $C_{l}$ for events maps like in Fig. 1 (bottom) is akin to correlating pixel windows. The Q-cumulants method employed in this study takes the two-particle correlation function in the azimuthal direction only. Therefore, the first main difference between the two approaches is what they correlate: for the former, it is $(\theta_{i},\phi_{i})$ windows, while for the latter it is $\phi_{i}$. The Monte Carlo simulations considered $f(\theta,\phi)$ as both factorizable and having the same azimuthal distribution for all $\theta$, so it is not surprising that $v_{n}\{C_{l}\}$ and $v_{n}\{2,QC\}$ would agree. In light of this fact, a second set of MC events was generated with $v_{n}$ coefficients varying slightly with $\theta$ and the result persisted.

In regards to correlations arising from jets and resonance decays, it escapes the range of this paper to answer how the angular power spectrum responds to such effects. Having said that, it must be emphasized that the scale structure of a jet, for instance, is quite small in comparison to elliptic eccentricity. Therefore, jets are expected to mainly influence higher $l$-modes and have negligible to no contribution to lower $l$. It should also be mentioned that the Q-cumulants calculation in this study did not take into consideration pseudorapidity gaps. That means non-flow contributions have not been suppressed when computing $v_{n}\{2,QC\}$ either. They may have different responses to non-flow effects, as their two-point functions are computed in distinct dimensions, though that is an issue for future research M. Machado (2019).

V Discussion

We have explored some of the powerful methods of Comic Microwave Background analyses when applied to the study of heavy ion collisions with very large particle multiplicities. We have shown that it is crucial to take into account the limitations of detector coverage as compared with the (almost) full-sky coverage of the Cosmic Microwave Background. Detector limitations introduce artificial structures in the angular power spectrum $C_{l}$ that can swamp the physical information. This holds in particular for multipoles of even $l$ with $m=0$. For the ALICE detector and the publicly available data used in the present study the coverage is roughly between 45*∘* and 135*∘*, leading to a suppression of $l=4n$, $n=1,2,\cdots$ relative to the otherwise enhanced even $l$’s. Other ranges of cuts in the polar angle with otherwise isotropic Monte Carlo distributions demonstrate clearly the relative suppressions of the (shifting) even $l$-modes due simply to the geometry of the detector limitations. The odd $l$-modes, although suppressed on average as compared to the even $l$-modes due to the approximate parity symmetry between multiplicities in the forward and backward directions, provide intriguing new information about the events. Here, there are very pronounced differences between the angular power spectrum of real data as compared to simple distributions based on approximate isotropy (which have essentially vanishing odd components on account of parity). Actual data seem to obey quite accurately a power-law behavior $C_{l}=A\cdot l^{-\beta}+\mathcal{C}$ over a wide range of odd $l$-values. The exponent $\beta$ appears to be constant, independent of centralities. What could be the origin of such scaling law, restricted to the odd-$l$ sector? When looking in closer details we find that this scaling appears to be directly triggered by the spread in interaction points of the two colliding heavy ions. Although a miniscule effect at the detector level, collisions that occur slightly shifted with respect to the center of the detector lead to a small forward-backward asymmetry in the angular coverage of the event. This issue will be addressed in more details in a forthcoming paper M. Machado (2019).

Finally, we have demonstrated that the angular power spectrum can be used to compute flow coefficients $v_{n}$. Again care must be taken in order to compensate for the limited range of the TPC detector. We have noted that such effects are encoded strongly in the $a_{l0}$ coefficients and an efficient way to eliminate detector limitations in the $\theta$-direction is to compute the angular power spectrum $C^{m\neq 0}_{l}$ without contributions from $a_{l0}$. Flow coefficients were then extracted using the power spectrum for $m\neq 0$ and compared to the cumulant method Bilandzic et al. (2011) for two-particle correlations.

These results show that analyzing heavy ion collisions by means of the angular power spectrum is a promising new avenue. As a measure of two-point correlations between $(\theta_{i},\phi_{i})$ windows, it could show that particle distributions have $\eta$ and $\phi$ event-by-event dependencies not seen when taking only azimuthal two-point correlations, even with $\eta$ gaps. We have here focused on the simplest of all observables, particle multiplicity. It would be most interesting to extend this analysis to explore how jets can impact the angular power spectrum or how it looks like when considering particles within different $p_{T}$ intervals.

Acknowledgements.

M.M. is most grateful to Ante Bilandzic for crucial help in the first part of this project. We thank Pavel Naselsky, Hao Liu and You Zhou for illuminating discussions. This work was supported in part by the Danish National Research Foundation (DNRF). The research of M.M. is supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).