Non-Gaussian elliptic-flow fluctuations in PbPb collisions at $\sqrt{\smash[b]{s_{_\text{NN}}}} =$ 5.02 TeV

Event-by-event fluctuations in the elliptic-flow coefficient $v_2$ are studied in PbPb collisions at $\sqrt{s_{_\text{NN}}} =$ 5.02 TeV using the CMS detector at the CERN LHC. Elliptic-flow probability distributions ${p}(v_2)$ for charged particles with transverse momentum 0.3 $<p_\mathrm{T}<$ 3.0 GeV and pseudorapidity $| \eta |<$ 1.0 are determined for different collision centrality classes. The moments of the ${p}(v_2)$ distributions are used to calculate the $v_{2}$ coefficients based on cumulant orders 2, 4, 6, and 8. A rank ordering of the higher-order cumulant results and nonzero standardized skewness values obtained for the ${p}(v_2)$ distributions indicate non-Gaussian initial-state fluctuation behavior. Bessel-Gaussian and elliptic power fits to the flow distributions are studied to characterize the initial-state spatial anisotropy.


Introduction
Ultrarelativistic heavy ion collisions at both the BNL Relativistic Heavy Ion Collider (RHIC) and the CERN Large Hadron Collider (LHC) create a hot and dense state of matter that consists of strongly interacting quarks and gluons, the "quark-gluon plasma" (QGP) [1][2][3][4][5][6][7]. Measurements of azimuthal particle correlations resulting from these collisions reveal properties of the QGP, but also of the initial state of a heavy-ion collision. In particular, the overall shape and fluctuations in the initial-state transverse energy density transformed by the hydrodynamic evolution of the medium into anisotropies in the final-state momentum space for the emitted particles [8][9][10], as reflected in the azimuthal chargedparticle density. The early RHIC measurements of the azimuthal correlations showed that the QGP could be described well by hydrodynamic models [11], with a shear viscosity to entropy density ratio (η/s) that is of the order of the lowest possible value for a quantum fluid [12,13]. The azimuthal charged-particle density can be characterized by a Fourier expansion, with (1) Here, the nth-order flow vector for a given event is v n ≡ (v n cos n , v n sin n ), where n is the angle of the intrinsic nth-order flow E-mail address: cms -publication -committee -chair @cern .ch. symmetry plane, as determined by the geometry of the participant nucleons. The experimentally accessible "event plane" angle, obs n , is based on the direction of maximum outgoing particle density and is, on average, in the same direction as n , but fluctuates about n because of resolution effects due to finite particle multiplicities.
By calculating the flow coefficients over a large number of events, the underlying probability distribution functions of individual Fourier coefficients can be determined. While the mean values of the v n distributions can be related to the overall shape of the interaction region, the higher order moments can be used to constrain the origin and the nature of the initial-state fluctuations and help disentangle the initial-state effects from the subsequent evolution of the medium [14,15]. Here, an event-by-event analysis is performed where it is possible to reduce the sensitivity of the results to nonflow correlations [16] and to clearly establish higher-order moments of the n = 2 (elliptic) distribution function.
The mean of this distribution, v 2 , is largely determined by the lenticular shape of the collision overlap region.
While the final-state particle distribution is characterized by the v n coefficients, the initial-state spatial anisotropy can be characterized by a harmonic expansion in terms of eccentricity vectors ε n [17][18][19][20]. For a given impact parameter, fluctuations in the initialstate transverse energy density lead to event-by-event differences in the orientation and magnitude of the ε n vectors with respect to the experimentally inaccessible "reaction plane," defined by the collision impact parameter and beam directions. The presence of a nonzero viscosity will degrade the correspondence between initial- and final-state anisotropies [11,21]. Still, an almost linear dependence is expected for the lowest order n = 2 [22-26] and n = 3 [9,18] harmonics, with v n = k n ε n [19]. Here, v n ≡ | v n |, ε n ≡ | ε n |, and k n is the flow response coefficient. The probability distribution functions of the magnitudes of the ε n vectors, p(ε n ), can be related to the corresponding p(v n ) distribution assuming a linear response, according to: where the k n term is expected to depend on the hydrodynamic evolution of the medium [27,28].
The elliptic-flow p(v 2 ) distribution can be characterized using the experimentally determined multiparticle cumulant flow har- The unitless standardized skewness of a probability distribution is a measure of the asymmetry about its mean. For the case of elliptic flow, the standardized skewness with respect to the reaction plane can be estimated using the cumulant flow harmonics as in Ref. [33]: Hydrodynamic calculations find this estimate to be in good agreement with the actual skewness except for the most peripheral events [33]. The standardized skewness estimate vanishes for fluctuations that arise from an isotropic Gaussian transverse initial-state energy density profile. In this case, the p(v 2 ) distribution is found by taking an integral over the azimuthal dependence of the two-dimensional Gaussian function [31,34]. The resultant, onedimensional distribution has a Bessel-Gaussian shape, where the [35][36][37], where the approximate equalities are within a few percent, suggests that the v 2 fluctuations can be well described by a two-dimensional Gaussian function [31].
Still, non-Gaussian fluctuations are expected in the initial-state energy density [33], which should lead to differences in the higher order cumulant coefficients. Such differences have been reported by the ATLAS Collaboration [16] in a similar measurement of peripheral PbPb collisions to that reported here. The precision of the LHC measurements allows for these differences to be explored in detail, giving a new method to investigate the initial-state behavior. The elliptic power function has been suggested to describe the asymmetric behavior of the p(ε n ) distributions [14, 15,38], noting that the Bessel-Gaussian distribution reproduces neither Glauber Monte Carlo nor IP-Glasma results other than for very central events [14]. This function is based on the assumption that the initial energy density profile of the collision is a superposition of N point-like, independent sources. In terms of the harmonic-flow coefficients and assuming a linear response, where ε 0 is approximately equal to the mean eccentricity in the reaction plane and α, which is approximately proportional to N, describes the size of the eccentricity fluctuations. The elliptic power distribution reduces to a Gaussian, Bessel-Gaussian, or power distribution form with the appropriate choice of parameters [39] and has the advantage of naturally incorporating the unit constraint on eccentricity, where | n | < 1. In this Letter, the p(v 2 ) distributions for charged particles in the pseudorapidity range |η| < 1.0 and with transverse momenta distributions, with these results used to estimate the standardized skewness of the flow distribution. Elliptic power and Bessel-Gaussian fits to the flow distributions are presented to gain further insight into the initial-state and its fluctuations.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.

Event and track selection
This analysis is based on a PbPb minimum bias data set with √ s NN = 5.02 TeV and corresponding to an integrated luminosity of 26 μb −1 , collected in 2015. The minimum-bias trigger used requires coincident signals in the HF calorimeters at both ends of the CMS detector with energy deposits above a predefined energy threshold of approximately 1 GeV and the presence of both colliding bunches at the interaction point as determined using beam pickup timing monitors. By requiring colliding bunches, events due to noise (e.g., cosmic rays and beam backgrounds) are largely suppressed. Events are further selected offline by requiring at least three towers with an energy above 3 GeV in each of the two HF calorimeters. The primary vertex for each event is chosen as the reconstructed vertex with the largest number of associated tracks. Primary vertices are required to have at least two associated tracks and to be located within 15 (0.2) cm of the nominal collision point along the longitudinal (transverse) direction. To suppress contamination from events with multiple collisions in the same bunch crossing (pileup), the procedure outlined in Ref.
Here, compatibility scores based on the number of pixel clusters with widths compatible with particles originating from each primary vertex are determined and events with primary vertices with compatibility scores below a predefined threshold are rejected as pileup. After applying the selection criteria, the average number of collisions per bunch crossing is less than ≈0.001 for the events used in this analysis, with a pileup fraction <0.05%.
Track reconstruction [40,42] is performed in two iterations to ease the computational load for high-multiplicity central PbPb collisions. The first iteration reconstructs tracks from signals ("hits") in the silicon pixel and strip detectors compatible with a trajectory of p T > 0.9 GeV/c. These tracks are required to have consistency with originating from the primary vertex, having a longitudinal association significance (d z /σ d z ) and a distance of closest approach significance (d 0 /σ d 0 ) each less than 3. In addition, the p T resolution [40,42] for each track, σ p T /p T , is required to be less than 10% and tracks are required to have at least 11 out of the 14 possible hits along their trajectory in the pixel and strip trackers. To reduce the number of misidentified tracks, which can occur when the hit pattern is consistent with more than one possible track solution, the chi-squared per degree of freedom, χ 2 /dof, associated with fitting the track trajectory through the different pixel and strip layers must be less than 0.15 times the total number of layers with hits along the trajectory of the track. The second iteration reconstructs tracks compatible with a trajectory of p T > 0.2 GeV/c using solely the pixel detector. These tracks are required to have longitudinal association significance d z /σ d z < 8 and a fit χ 2 /dof value less than 12 times the number of layers with hits along the trajectory of the track. In the final analysis, first iteration tracks with p T > 1.0 GeV/c are used together with pixel-detector-only tracks with p T < 2.4 GeV/c after removing duplicates. Track reconstruction for the merged iterations has a combined geometric acceptance and efficiency exceeding 60% for p T ≈ 1.0 GeV/c and |η| < 1.0. When the track p T is below 1 GeV/c, the acceptance and efficiency steadily drops, reaching approximately 40% at

Analysis technique
Analyses of flow harmonics using multiparticle cumulants were initially introduced as a way to minimize nonflow effects [30]. These analyses have been based on either the generating function formalism [30] or, more recently, through direct calculation [43]. The unfolding procedure employed here, as introduced by the ATLAS collaboration [16], is expected to give similar results to a multiparticle cumulant analysis, but with reduced sensitivity to multiplicity fluctuations and nonflow effects [44].
The event-by-event v 2 coefficients and phases in Eq.
(1) can be estimated with v obs where φ i is the azimuthal angle of the track, obs 2 is the event plane angle for the 2nd harmonic, the angular brackets denote an efficiency weighted average over all particles in a given range of phase space for an event, and w i = 1/ε i is the inverse of the tracking efficiency ε i (p T , η) of the ith track. The analysis does not require the explicit calculation of the event plane angle for each event. In the absence of particle correlations unrelated to the hydrodynamic flow behavior ("nonflow"), the observed event-byevent flow vectors of Eq. (6) will approach the true underlying flow vectors as the particle multiplicity becomes large. In addition to the efficiency weighting, a standard recentering procedure [45], where the event average x-and y-components of the flow vector are required to equal zero, is applied to further suppress acceptance biases.
Events are sorted into different centrality classes, as determined by the transverse energy deposited in the HF calorimeters [6], and the magnitudes of the estimated flow vectors are used to construct the "observed" p(v obs 2 ) distributions for each class. Finite particle multiplicities result in a statistical fluctuation of the v obs 2 estimate for a given event about the true underlying v 2 value by a response function p(v obs 2 |v 2 ). This, in turn, results in a p(v obs 2 ) distribution that is broader than the underlying p(v 2 ) behavior. The observed distribution can be expressed as a convolution of the underlying flow behavior and the response function A data-based technique, first introduced by the ATLAS Collaboration [16], was used to build the response function in Eq. (7). This technique divides the full event sample into two symmetric subevents (a and b) based on pseudorapidity. Given that v 2 (η) is symmetric about η = 0 on average for the symmetric PbPb system, the physical flow signal cancels in the distribution of flow vector differences from each subevent p( v a n − v b n ). The resulting distribution contains residual effects from multiplicity-related fluctuations and nonflow effects [44] and provides a basis for building the response function. The ability of the analysis procedure to suppress nonflow effects was studied by introducing a v 2 signal on top of hijing 1.383 [46] simulated events, which contain nonflow. The EbyE analysis is found to recover the "truth" to within 0.1%. To unfold the effects of multiplicity-related fluctuations, the D'Agostini iterative method with early stopping (regularization) [47][48][49] was used to obtain a maximum likelihood estimate of the underlying p(v 2 ) behavior. The analysis was done using the RooUnfold [50] package of the root data analysis framework [51]. The unfolding procedure becomes increasingly sensitive to statistical fluctuations when the number of iterations is allowed to run to large values, resulting in unphysical oscillations in the low event count tails of the unfolded distribution. The regularization criterion used to suppress these oscillations is to apply the response function to each unfolding iteration ("refolding") and compare the resulting distribution to the observed one. Iterations are stopped when the χ 2 /dof between the refolded and observed distribution is approximately equal to one. After this final unfolding iteration is reached, the resulting distribution is truncated above v 2 + 4σ v 2 to further suppress any residual artifacts in the tails that result from the unfolding procedure. Representative final unfolded distributions are shown in Fig. 1. In addition, p(v obs are plotted for each centrality to illustrate the statistical resolution effects present prior to unfolding. The fits shown in Fig. 1 are discussed in Section 6.

Systematic uncertainties
A number of potential sources of systematic uncertainties for the v 2 {m} values extracted from the unfolded p(v 2 ) distributions were considered. The systematic uncertainties that arise from the vertex z position were investigated by splitting the default vertex range into two windows of |z vtx | < 3.0 cm and 3.0 < |z vtx | < 15.0 cm and comparing the results from the two ranges. The resulting uncertainties range from 5% for central events, decreasing to 0.5% for mid-central events. To estimate the bias from misidentified tracks, the track quality criteria described in Section 3 were varied. Two scenarios were considered, with one increasing and the other decreasing the probability of misidentifying a track. The results of these two scenarios were compared to the values obtained in the default analysis. The resulting uncertainties range from 2% for central events to 1% for mid-central events. To estimate the systematic uncertainty in the choice of response function, the unfolding procedure was repeated using an analytic response function obtained from a Gaussian fit to the data-driven statistical resolution distribution [16]. The resulting uncertainties are 3% for central events and decrease to 1% for mid-central events. Other sources of potential systematic bias were explored and found to be negligible. To assess the potential bias from residual pileup events, the threshold for determining pileup events was raised to decrease the probability of including events with multiple collisions in the analysis. The bias from unfolding regularization was studied by modifying the χ 2 /dof goodness-of-fit regularization criteria and comparing the cases when the refolding χ 2 /dof cutoff is 2.0 relative to when it is 1.0. To test the potential bias that might result from the 4σ truncation of the final unfolded distributions, the truncation point was varied between 3.5σ and 4.5σ . To assess the uncertainty on the choice of the prior, the unfolding was repeated using priors that were systematically transformed to have 10% larger and smaller means than the default prior. No significant bias was found with these variations of the prior. The total systematic uncertainties were obtained by adding the contribution from each source in quadrature. The v 2 values calculated for the different cumulant orders have a total systematic uncertainty of the order of 5% for central collisions, which decreases to 1% in midcentral collisions. For most centralities, the uncertainties are smaller than the symbol size.
As all of the systematic uncertainties are expected to be correlated between the different cumulant orders, with the same data used in the calculation of each order, all of the above studies were also performed for the ratios of different orders and for the skewness estimate given by Eq. (4). For the ratios, the total systematic uncertainty is found as 1% for central collisions, decreasing to 0.1% for mid-central collisions. The standardized skewness is very sensitive to small fluctuations in the cumulant flow harmonics, resulting in a systematic uncertainty of 100% for central collisions that reduces to 20% for mid-central collisions.

Results
The cumulant elliptic-flow harmonics obtained from the moments of the unfolded p(v 2 ) distributions using Eq. (3) are shown in Fig. 2 for cumulant orders 2, 4, 6, and 8. It was not possible to obtain 0-5% central results for v 2 {4} and v 2 {6} because the righthand side of Eq. (3) was found to be negative for these values. This behavior might be a consequence of volume fluctuations dominating the cumulant behavior for these central events, as discussed in Ref. [52]. The cumulant results exhibit the previously observed The centrality-dependent ratios for the elliptic-flow coefficients obtained for different cumulant orders are shown in Fig. 3. For most centrality ranges, the ratios indicate a rank ordering of the cumulants, with differences on the order of a few percent and with Glauber initial conditions [53] and an η/s value of 0.08 is shown by the shaded band. This simulation is for pions with 0.2 < p T < 3.0 GeV/c in PbPb collisions at √ s NN = 2.76 TeV [33]. Also shown are results from the ATLAS Collaboration [37] for PbPb collisions at 2.76 TeV and for charged particles with 0.5 < p T < 20.0 GeV/c and |η| < 2.5. The calculation is consistent with the experimental results found at both beam energies. The similarity between experimental results with 2.76 and 5.02 TeV is consistent with the small changes in the initial-state eccentricities expected between these energies [54] and the expectation that the cumulant flow harmonic ratios follow those of the corresponding eccentricity ratios [33]. values for PbPb collisions at 2.76 TeV from Ref. [33] are also shown and found to be consistent with the current measurements. Within the hydrodynamic model and allowing for a finite skewness of the event-by-event v 2 distribution, the small splitting between the cumulant orders is ex- [33]. Experimentally, we find a value for this splitting ratio of 0.143 ± 0.008 (stat) ± 0.014 (syst) for 20-25% central events, with the ratio increasing to 0.185 ± 0.005(stat) ± 0.012(syst) as the centrality increases to 55-60%. The observed values might suggest higher order terms in a cumulant expansion of the v 2 distribution are required to account for the skewness. This relationship was recently examined by the ALICE collaboration in Ref. [55] using a q-cumulant analysis, with results comparable to the findings in this paper when considering systematic uncertainties and a different kinematic range for the ALICE measurement.
Both elliptic power and Bessel-Gaussian parametrizations used for fits such as shown in Fig. 1 assume a linear response between eccentricity and flow, but only the elliptic power law allows for a finite skewness. For a Bessel-Gaussian distribution, the skewness is equal to zero. This feature results in the elliptic power function being in better agreement with the observed fluctuation behavior The fit parameters for the elliptic power function are shown in Fig. 5 for the different centrality bins. As also found in Ref. [15], the fits do not converge for central collisions where the distributions become very close to a Bessel-Gaussian form. Consequently, the parameters are shown for centralities >15%. The experimental k 2 values show only a weak centrality dependence. Viscous hydrodynamic calculations indicate that deviations from thermal equilibrium should lead to a reduced correspondence between the initial-state geometry and the flow signal in peripheral collisions [27,28]. This effect is suggested in Fig. 5 by the decrease in the k 2 value with increasing centrality, although the systematic uncer-  [15] using viscous hydrodynamics with Glauber initial conditions and an η/s value of 0.19 to determine the response coefficient k 2 . Glauber (blue shaded band) and IP-Glasma (red shaded band) model calculations from Ref. [15] are shown for the α and ε 0 parameters. The systematic uncertainties account for the highly correlated parameters of the elliptic power function fit and for the bin-to-bin correlations in the unfolded distributions introduced by the unfolding procedure.
tainties are too large for this to be a definitive observation. The calculated decrease is greater than observed, although within the systematic uncertainties of the measurement. The eccentricity parameter of the power law fit, ε 0 , is found to first increase, and then level off with increasing centrality. The leveling occurs for centralities > 40%, which is also where the v 2 values start to level off and then decrease. The α parameter, which reflects the number of sources in the power-law fit, is found to steadily decrease with increasing centrality, as expected.
Theoretical predictions at 2.76 TeV from Ref. [15] are compared to the current analysis in Fig. 5. A viscous hydrodynamic calculation with Glauber initial conditions and an η/s value of 0.19 is in agreement with the experimental k 2 values. This coefficient is expected to have only a weak dependence on the initial state, with its centrality dependence largely determined by the viscosity of the medium [15]. Predictions obtained using Glauber and IP-Glasma [56,57] initial conditions, where the IP-Glasma model includes gluon saturation effects, are shown for the ε 0 and α parameters. These latter two calculations qualitatively capture the observed behavior for the α-parameter, but a significant difference is found in comparing the theoretical ε 0 values with experiment.
This difference might reflect a nonlinear response term, which will alter the magnitude of the flow response coefficient and consequently the ε 0 and α parameters, as suggested in Ref. [15].

Summary
In summary, a non-Gaussian behavior is observed in the eventby-event fluctuations of the elliptic flow v 2 coefficients in PbPb collisions recorded by the CMS detector at √ s NN = 5.02 TeV. The probability distributions p(v 2 ) for 5%-centrality bins between 5% and 60% centrality are found by unfolding statistical resolution effects from measured flow distributions. The v 2 coefficients corresponding to different cumulant orders are calculated from the Based on the elliptic power function fits, the centrality dependence of the flow response coefficient, which relates the final state geometry to the initial state energy density distribution, is found to be consistent with model calculations. However, the observed eccentricities are smaller than predictions based on either the Glauber model or the IP-Glasma model initial conditions with an assumed linear flow response. This difference might indicate the need for a nonlinear response term. The current results illustrate that LHC experiments now have the precision to explore the details of the initial-state fluctuations.

Acknowledgements
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses.