$p_T$-Dependent Particle Number Fluctuations From Principal Component Analyses in Hydrodynamic Simulations of Heavy-Ion Collisions

TL;DR

This study uses principal component analysis on hydrodynamic simulations of heavy-ion collisions to compare flow and multiplicity fluctuations with experimental data, revealing key differences in $p_T$ dependence.

Contribution

It introduces a toy model to analytically understand PCA modes and their dependence on fluctuations, highlighting discrepancies in $p_T$ trends between simulations and experiments.

Findings

Flow PCA components match experimental trends.

Multiplicity PCA shows simulation-experiment $p_T$ trend differences.

Hydrodynamic models exhibit larger transverse momentum fluctuations.

Abstract

We carry out a principal component analysis of fluctuations in a hydrodynamic simulation of heavy-ion collisions, and compare with experimental data from the CMS collaboration. The leading and subleading principal components of elliptic and triangular flow reproduce the trends seen in data. By contrast, the principal components of multiplicity fluctuations show an interesting difference in their $p_T$ dependence for simulations compared to experimental data. Specifically, the leading component increases with $p_T$ in hydrodynamics, while it is constant in experiment. In order to understand how the leading and subleading modes arise, we construct a toy model where the principal components have a simple analytic form. We show how the PCA components depend on fluctuations of the average transverse momentum and of the total multiplicity, as well as correlations between the two, and we verify that hydrodynamic simulations agree with the predictions of the toy model. The difference in the momentum trend is likely due to the fact that hydrodynamic models typically have transverse momentum fluctuations that are larger than seen experimentally.

Tables

Table 1. Table 1: Values of the variances and covariance of N 𝑁 N and p ¯ T subscript ¯ 𝑝 𝑇 \bar{p}_{T} at LHC and RHIC in our hydrodynamical calculation using NeXSPheRIO. The number in italics are approximate values extracted from CMS data through a rough fit of the principal components v 0 ( 1 ) ​ ( p T ) superscript subscript 𝑣 0 1 subscript 𝑝 𝑇 v_{0}^{(1)}(p_{T}) and v 0 ( 2 ) ​ ( p T ) superscript subscript 𝑣 0 2 subscript 𝑝 𝑇 v_{0}^{(2)}(p_{T}) , shown in the top panels of Fig.  2 , using Eq. ( 18 ).

energy	centrality		$\frac{{\sigma }_N}{⟨N⟩}$	$\frac{{\sigma }_{p_T}}{⟨{\overline{p}}_T⟩}$	$\sqrt{\frac{⟨\delta N\delta {\overline{p}}_T⟩}{⟨N⟩⟨{\overline{p}}_T⟩}}$
2.76 TeV	0-5 %	hydro	0.12	0.026	0.041
		CMS	0.09	0.010	0.000
	20-30 %	hydro	0.16	0.041	0.070
		CMS	0.13	0.019	0.020
200 GeV	0-10 %	hydro	0.11	0.017	0.017
	20-30%	hydro	0.12	0.025	0.031

Equations

$$\frac{dN}{d\vec{p}}$$

$$=\sum _{n=-\infty }^{+\infty }{\mathcal{V}}_n(p_T)e^{-in\varphi }$$

$${\mathcal{V}}_{n\Delta }(p_T^a,p_T^b)\equiv ⟨{\mathcal{V}}_n(p_T^a){\mathcal{V}}_n^{\ast }(p_T^b)⟩-⟨{\mathcal{V}}_n(p_T^a)⟩⟨{\mathcal{V}}_n^{\ast }(p_T^b)⟩.$$

$${\mathcal{V}}_{n\Delta }(p_T^a,p_T^b)=\sum _{\alpha }{\mathcal{V}}_n^{(\alpha )}(p_T^a){\mathcal{V}}_n^{(\alpha )}(p_T^b),$$

$${\mathcal{V}}_n(p_T)=\sum _{\alpha }{\xi }_n^{(\alpha )}{\mathcal{V}}_n^{(\alpha )}(p_T),$$

$$v_n^{(\alpha )}(p_T)=\frac{{\mathcal{V}}_n^{(\alpha )}(p_T)}{⟨{\mathcal{V}}_0(p_T)⟩}.$$

$$\frac{1}{2\pi }\frac{dN}{dyd{p}_T}={\mathcal{V}}_0(p_T)=\frac{2p_{T}N}{\pi {\bar{p}}_T^2}{e}^{-\frac{2p_T}{{\bar{p}}_T}}$$

$$N$$

$${\bar{p}}_T$$

$$\frac{\delta {\mathcal{V}}_0(p_T)}{⟨{\mathcal{V}}_0(p_T)⟩}=\frac{\delta N}{⟨N⟩}-2\frac{\delta {\bar{p}}_T}{⟨{\bar{p}}_T⟩}+2\frac{p_T\delta {\bar{p}}_T}{⟨{\bar{p}}_T{⟩}^2}.$$

$${\mathcal{V}}_{0\Delta }(p_T^a,p_T^b)\equiv ⟨\delta {\mathcal{V}}_0(p_T^a)\delta {\mathcal{V}}_0(p_T^b)⟩,$$

$$\frac{{\mathcal{V}}_{0\Delta }(p_T^a,p_T^b)}{⟨{\mathcal{V}}_0(p_T^a)⟩⟨{\mathcal{V}}_0(p_T^b)⟩}$$

$$v_0^{(\alpha )}(p_T)=a^{(\alpha )}+b^{(\alpha )}\frac{p_T}{⟨{\bar{p}}_T⟩},$$

$$v_0^{(1)}(p_T)\simeq \frac{{\sigma }_N}{⟨N⟩}.$$

$$v_0^{(1)}(p_T)$$

$$v_0^{(2)}(p_T)$$

$$\frac{{\mathcal{V}}_{0\Delta }(p_T^a,p_T^b)}{⟨{\mathcal{V}}_0(p_T^a)⟩⟨{\mathcal{V}}_0(p_T^b)⟩}=v_0^{(1)}(p_T^a)v_0^{(1)}(p_T^b)+v_0^{(2)}(p_T^a)v_0^{(2)}(p_T^b).$$

Full text

I Introduction

The expansion of the matter formed in nucleus-nucleus collisions at relativistic energies produces a collective transverse flow. This flow is the response to the density gradients in the initial fireball. It is azimuthally asymmetric because the initial fireball is anisotropic and contains hot spots. These inhomogeneities are of interest: they reflect the poorly known mechanism of energy deposition, via the strong interaction, when two nuclei collide, and their influence on the final flow depends on fluid properties, which are also poorly known (e.g. shear and bulk viscosities). A lot of work has been done to relate initial inhomogeneities and final flow of produced particles. In particular the mapping between initial conditions and anisotropic flow has been studied globally and event-by-event Gardim:2011xv ; Gardim:2014tya ; Niemi:2012aj ; Niemi:2015qia ; Fu:2015wba ; Noronha-Hostler:2015dbi . To get more detailed information on fluctuations in the initial state, a useful observable is the factorization breaking ratio Gardim:2012im ; Kozlov:2014fqa ; Heinz:2013bua ; Khachatryan:2015oea ; Gardim:2017ruc ; Zhao:2017yhj ; McDonald:2016vlt , which encodes the correlations of flow harmonics at different transverse momenta or pseudorapidities. More recently a new more precise tool was proposed, the Principal Component Analysis (PCA) for event-by-event fluctuations Bhalerao:2014mua ; Mazeliauskas:2015vea ; Mazeliauskas:2015efa and first experimental results for such an analysis have been presented by the CMS collaboration Sirunyan:2017gyb . The aim of this paper is to present a hydrodynamical study of these observables and point out an interesting difference between data and some hydrodynamic simulations for the $n=0$ leading and subleading components, corresponding to multiplicity fluctuations. These components are sensitive to physics not explored by anisotropic flow and can put new constraints on initial conditions models, in particular on the transverse size of the fireball and its fluctuations.

II Principal Component Analysis

Principal Component Analysis is a common technique for finding patterns in data of high dimension. One tries to find new variables that incorporate as much as possible of the variations. This amounts to diagonalizing the covariance matrix (e.g. pcaref1 ). It was first suggested to use it to study event-by-event fluctuations in relativistic nuclear collisions in Bhalerao:2014mua . Consider a set of collisions or events. For each event, the single particle distribution can be expanded as

[TABLE]

where $d\vec{p}=dydp_{T}d\phi$, $\phi$ is the azimuthal angle of the particle momentum. $\mathcal{V}_{n}(p_{T})$ is a Fourier coefficient (without the usual normalization by multiplicity) which is complex for $n\not=0$. Its magnitude and orientation vary for each event.

For each transverse momentum bin, the variance can be computed $\langle|\mathcal{V}_{n}(p_{T}^{a})|^{2}\rangle-|\langle\mathcal{V}_{n}(p_{T}^{a})\rangle|^{2}$ (the average is performed over events) but brings no information about possible relationship between different bins. To investigate how different bins are correlated, one constructs the covariance matrix:

[TABLE]

The terms $\langle\mathcal{V}_{n}(p_{T}^{a})\rangle$ are zero by azimuthal symmetry, except for $n=0$.

This covariance matrix is real, symmetric, positive-semidefinite. It can be diagonalized and re-written in term of its real orthogonal eigenvectors $\mathcal{V}_{n}^{(\alpha)}(p_{T})$

[TABLE]

from where one can express the flow vector in a given event as

[TABLE]

where $\xi_{n}^{(\alpha)}$ are coefficients that vary from event to event (specifically, uncorrelated, random complex numbers with zero mean and unit variance). Terms in the right-hand side of Eq. (4) are ordered according to the magnitude of the eigenvalues. Even by truncating the sum to the first two or three terms, one typically obtains a very good approximation to the left-hand side. The largest component $\mathcal{V}_{n}^{(1)}(p_{T})$ is called the leading mode, $\mathcal{V}_{n}^{(2)}(p_{T})$ the subleading mode, etc. For comparison with standard flow, it is useful to introduce the following scaled principal components

[TABLE]

Once the dominant terms in Eq. (4) are determined (i.e. patterns are found in our high dimension data), the physical meaning of these terms must be investigated. This was done in Bhalerao:2014mua ; Mazeliauskas:2015vea ; Mazeliauskas:2015efa ; Bozek:2017thv ; Liu:2019jxg and is discussed in section III ($n=2,3$) and IV ($n=0$).

III Results for anisotropic flow

In this section and the next, we present results obtained from a hydrodynamic simulation for a perfect fluid expanding in 3+1 dimensions starting from NeXus initial conditions Drescher:2000ha . The code used, NeXSPheRIO, has been shown to lead to a consistent description of many flow data at top RHIC energies Qian:2007ff ; Andrade:2006yh ; Andrade:2008xh ; Andrade:2008fa ; Gardim:2011qn ; Gardim:2012yp ; Takahashi:2009na ; Qian:2012qn .

We also have some data accumulated for two centrality windows (0-5 and 20-30%, where centrality is defined according to the number of participant nucleons) at $\sqrt{s}=2.76$ TeV and their compatibility with flow observables more subtly related to fluctuations (scaled harmonic flow distributions, factorization breaking ratio) has been tested RHICdist . This code is therefore an interesting tool for a first investigation of the PCA results obtained recently by CMS at the LHC Sirunyan:2017gyb .

For $n$=2–3, we show the first two scaled principal components and comparison with CMS data in Fig. 1. Our cuts are $|\eta|<2.5$ (equivalent to CMS) but $p_{T}>0.5$ GeV, slightly higher than CMS $p_{T}>0.3$ GeV. We used similar $p_{T}$ bins as experimentally.

The leading component is straightforward to interpret Bhalerao:2014mua ; Mazeliauskas:2015vea ; Mazeliauskas:2015efa . If it dominates, Eq. (4) yields $\mathcal{V}_{n\Delta}(p_{T}^{a},p_{T}^{b})\sim\mathcal{V}_{n}^{(1)}(p_{T}^{a})\mathcal{V}_{n}^{(1)}(p_{T}^{b})$ i.e. there is flow factorization. The event flow defined by Eq. (5) reduces to $\mathcal{V}_{n}(p_{T})\sim\xi_{n}^{(1)}(p_{T})\mathcal{V}_{n}^{(1)}(p_{T})$, i.e., the leading component corresponds to usual anisotropic flow. Concentrating on the region from 0 to 2 GeV, we see that our hydro simulation slightly overestimates the leading components. Inclusion of viscosity would damp them and improve agreement with data, as explicitly shown for the $p_{T}$-integrated $n=3$ leading component in Ref. Mazeliauskas:2015vea .

Higher-order principal components encode the information about the momentum dependence of flow fluctuations. They are in particular responsible for the breaking of factorization of two-particle correlations Gardim:2012im . This effect is often quantified using the factorization breaking ratio $r_{n}$ Khachatryan:2015oea , which is a function of two variables $p_{T}^{a}$ and $p_{T}^{b}$. Higher-order principal components express the same information in a simpler way, because they are functions of a single variable $p_{T}^{a}$. We only show the subleading component. In the range 0 to 2 GeV, our simulations capture the main features of the data. The subleading component changes sign as a function of $p_{T}$, which is imposed by orthogonality with the leading mode. The fact that our $p_{T}$ cut is slightly higher than in data shifts this crossing point to the right. Note that inclusion of viscosity is not expected to change significantly the magnitude of the subleading mode, as was shown for $n=3$ in Ref. Mazeliauskas:2015vea .

IV Results for multiplicities

We now discuss multiplicity fluctuations, corresponding to $n=0$ principal components. The comparison between our results and CMS data is displayed in the top panels of Fig. 2. There is rough overall agreement, but not as good as in Fig. 1. The leading component is rather independent of $p_{T}$ in experiment, while it increases with $p_{T}$ in our hydrodynamic calculation. The increase is less strong at RHIC energies (bottom panel of Fig. 2). The increase at LHC energies is not specific to our implementation, as it has been seen by other groups Mazeliauskas:2015efa ; Matt . Such qualitative disagreement between hydrodynamics and experimental data is rare, therefore, we investigate its origin in detail.111Note that the transport model AMPT without hydrodynamics predicts a flat leading component, as seen in data.

In order to understand the principal components for $n=0$, we introduce a toy model where the fluctuation of the multiplicity in a $p_{T}$ bin originates from two sources:

1. fluctuations of the total multiplicity $N$.
2. fluctuations of the mean transverse momentum $\bar{p}_{T}$. We assume that the $p_{T}$ spectrum is exponential:

[TABLE]

where $N$ is the total multiplicity per unit rapidity and $\bar{p}_{T}$ is the mean transverse momentum in one event. Next, we allow $N$ and $\bar{p}_{T}$ in a given event to deviate from the event-averaged total multiplicity $\langle N\rangle$, and the event-averaged mean transverse momentum $\langle\bar{p}_{T}\rangle$ in a centrality bin, respectively:

[TABLE]

Expanding Eq. (7) to first order in $\delta N$ and $\delta\bar{p}_{T}$, one obtains:

[TABLE]

The covariance (3) is then given by

[TABLE]

where angular brackets denote an average over events in a centrality bin. Inserting Eq. (10) into Eq. (11), one obtains:

[TABLE]

where $\sigma_{N}^{2}\equiv\langle\delta N^{2}\rangle$ and $\sigma_{p_{T}}^{2}\equiv\langle\delta\bar{p}_{T}^{2}\rangle$ denote the variance of the multiplicity and mean $p_{T}$, respectively. Inspection of the dependence on $p_{T}^{a}$ and $p_{T}^{b}$ shows that the scaled principal components defined by Eqs. (4) and (6) can only be of the form

[TABLE]

i.e., they are linear in $p_{T}$. Since they span a two-dimensional space, this in turn implies that there are at most two principal components (remember that principal components are mutually orthogonal). The full analytic expressions of these principal components are cumbersome. Therefore, we make further simplifying assumptions, by identifying the leading terms in Eq. (14).

Table 1 gives the values of the relative fluctuations of $N$ and $\bar{p}_{T}$ in our hydrodynamic calculation, as well as their covariance. The relative fluctuations of $N$ are larger by an order of magnitude, which is explained by the large width of the centrality bin. In the limit where $\sigma_{p_{T}}$ and $\langle\delta N\delta\bar{p}_{T}\rangle$ can be neglected, only the first term remains in the right-hand side of Eq. (14). The covariance matrix trivially factorizes, i.e., there is only one principal component. The scaled principal component, defined by Eq. (6), is:

[TABLE]

It is independent of $p_{T}$. Thus, the fact that our hydrodynamic calculation reproduces the magnitude of $v_{0}^{(1)}(p_{T})$ at low $p_{T}$ (i.e., for the bulk of produced particles) simply means that it has the correct multiplicity fluctuations. These are largely dominated by the width of the centrality bin used for the analysis, or, equivalently, by impact parameter fluctuations.

We now consider the more general case where $\sigma_{p_{T}}/\langle\bar{p}_{T}\rangle$ and $\langle\delta N\delta\bar{p}_{T}\rangle/\langle N\rangle\langle\bar{p}_{T}\rangle$ are not zero, but can still be treated as small quantities. Then, to leading order in these quantities, the scaled principal components are:

[TABLE]

One can check that with these expressions, the decomposition (4) is satisfied. In terms of the scaled components, this equation can be written:

[TABLE]

Inserting Eq. (18) into Eq. (19), and expanding to first order in $\langle\delta N\delta\bar{p}_{T}\rangle$ and $\sigma_{p_{T}}^{2}$, one recovers Eq. (14) except for the second and third terms of the first line, which are subleading corrections to the first term.

Equation (18) is a refinement of the zeroth-order result, Eq. (16). A subleading mode $v_{0}^{(2)}(p_{T})$ appears, which is directly proportional to $\sigma_{p_{T}}/\langle\bar{p}_{T}\rangle$. The connection between the subleading mode and $p_{T}$ fluctuations was already made in Ref. Mazeliauskas:2015efa . The change of sign of the subleading mode occurs at $p_{T}=(3/4)\langle\bar{p}_{T}\rangle$, which depends little on centrality. Figure 2 displays a comparison between Eq. (18) and the result from the full hydrodynamic calculation. Agreement is very good at RHIC and a little worse at LHC (presumably due to the different lower $p_{T}$ cuts). We therefore conclude that Eq. (18) captures the physics of the first two $n=0$ modes.

Using CMS data on principal components, shown in Fig. 2, one can estimate the quantities appearing in the right-hand side of Eq. (18). The corresponding numbers are reported in Table 1, and should be considered rough figures. As explained above, the value of $\sigma_{N}/\langle N\rangle$ is given by the value of $v_{0}^{(1)}(p_{T})$ at low $p_{T}$. The value of $\langle p_{T}\rangle$ is inferred from the value of $p_{T}$ for which $v_{0}^{(2)}(p_{T})$ crosses the horizontal line. This gives $\langle p_{T}\rangle\sim 0.75~{}GeV$, in reasonable agreement with the value $0.81~{}GeV$ obtained by direct integration of $p_{T}$-spectra in the same range (0.3 to 3 GeV) Acharya:2018qsh . The value of $\sigma_{p_{T}}/\langle\bar{p}_{T}\rangle$ is then estimated by fitting $v_{0}^{(2)}(p_{T})$ at low $p_{T}$, and the resulting values agree with those from a dedicated analysis Abelev:2014ckr . The covariance $\langle\delta N\delta\bar{p}_{T}\rangle/\langle N\langle\langle\bar{p}_{T}\rangle$ is finally inferred from the $p_{T}$ dependence of $v_{0}^{(1)}(p_{T})$. While the values of $\sigma_{N}/\langle N\rangle$ from the hydrodynamic calculation are in reasonable agreement with data, values of $\sigma_{p_{T}}/\langle\bar{p}_{T}\rangle$ are too large by a factor $\sim 2$, and the discrepancy is even worse for the covariance. We come back to this point below.

The motivation for building the toy model was to understand under which condition the leading mode is independent of $p_{T}$, or rises with $p_{T}$. The first line of Eq. (18) shows that a rise with $p_{T}$ can be ascribed to a positive correlation between the mean transverse momentum and the multiplicity, represented by the quantity $\langle\delta N\delta\bar{p}_{T}\rangle$. The fact that this rise is seen in hydrodynamic calculations, not in data, implies that hydrodynamic calculations overestimate $\langle\delta N\delta\bar{p}_{T}\rangle$, as illustrated by the numbers in Table 1 (the covariance extracted from CMS data for central collisions is compatible with 0). This can be related to the fact that hydrodynamic models yield too large $\delta\bar{p}_{T}$ in general, as pointed out by a study of transverse momentum fluctuations Bozek:2012fw . Since transverse momentum fluctuations in hydrodynamics originate from fluctuations in the transverse size of the interaction region Broniowski:2009fm ,222At a given centrality, a smaller size implies a larger density and temperature, hence a larger mean transverse momentum., this in turn implies that existing models of initial fluctuations tend to overestimate the size fluctuations.

The conclusion of this study is that a model which predicts the right multiplicity and $p_{T}$ fluctuations should capture the first two principal components for $n=0$. The reason why our hydrodynamical model predicts a rise of the leading mode with $p_{T}$, which is not seen in data, can be related to the fact that that $\sigma_{p_{T}}$ is larger in our model than in data. Despite the fact that our hydrodynamic calculation at 2.76 TeV does not reproduce CMS data, we expect our predictions for 200 GeV collisions, shown in the bottom panels of Fig. 2, should correctly predict the first two modes of multiplicity fluctuations at RHIC. The reason is that the values of $\sigma_{N}/\langle N\rangle$ from Table 1 are comparable with experimental values from PHENIX Adare:2008ns , and the values of $\sigma_{p_{T}}/\langle\bar{p}_{T}\rangle$ are slightly too large compared to STAR data Adams:2005ka , but in fair agreement.

V Conclusion

We have compared results from a hydrodynamic simulation using the code NeXSPheRIO with recent experimental data by CMS, on the Principal Component Analysis. The trends for the leading and subleading components of elliptic and triangular flow are in fair agreement with data. In contrast, for multiplicity fluctuations, we have pointed out a qualitative disagreement: The leading component increases with $p_{T}$ in hydrodynamics (here as well as in Mazeliauskas:2015efa ; Matt ) while it is constant in data at LHC energies. We have constructed a toy model which gives result in good agreement with the full hydrodynamic calculation. In this toy model, the subleading component is proportional to the standard deviation of the mean $p_{T}$, $\sigma_{p_{T}}$. The leading component is close to $\sigma_{N}/\langle N\rangle$ at low $p_{T}$, but increases with $p_{T}$ if the fluctuations of $p_{T}$ are large and correlated with the fluctuations of the multiplicity.

We have thus related $n=0$ results from the principal component analysis to multiplicity and transverse momentum fluctuations. Fluctuations in $N$ and $\bar{p}_{T}$ have been attracting attention for a long time because they may probe the QCD phase transition (see e.g. Stephanov:2005iu ), as well as initial inhomogeneities (see for example Gavin:2011gr ; Broniowski:2009fm ). The principal components are sensitive not only to the width of multiplicity and transverse momentum fluctuations, but also to their mutual covariance. They open a new window on initial fluctuations, which can be used to rule out initial condition models.

Acknowledgments

F.G. acknowledges support from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP grant 2018/18075-0) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq grant 310141/2016-8). F.G.G. was supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq grant 205369/2018-9 and 312932/2018-9) and FAPEMIG (grant APQ-02107-16). P.I. thanks support from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). M.L. acknowledges support from FAPESP projects 2016/24029-6 and 2017/05685-2. F.G., F.G.G. and M.L. acknowledge support from project INCT-FNA Proc. No. 464898/2014-5 and F.G., M.L. and J.-Y.O. from USP-COFECUB (grant Uc Ph 160-16, 2015/13).