Measurement of higher-order harmonic azimuthal anisotropy in PbPb collisions at a nucleon-nucleon center-of-mass energy of 2.76 TeV

Measurements are presented by the CMS Collaboration at the Large Hadron Collider (LHC) of the higher-order harmonic coefficients that describe the azimuthal anisotropy of charged particles emitted in sqrt(s[NN]) = 2.76 TeV PbPb collisions. Expressed in terms of the Fourier components of the azimuthal distribution, the n = 3-6 harmonic coefficients are presented for charged particles as a function of their transverse momentum (0.3<pt<8.0 GeV), collision centrality (0-70%), and pseudorapidity (abs(eta)<2.0). The data are analyzed using the event plane, multiparticle cumulant, and Lee-Yang zeros methods, which provide different sensitivities to initial-state fluctuations. Taken together with earlier LHC measurements of elliptic flow (n = 2), the results on higher-order harmonic coefficients develop a more complete picture of the collective motion in high-energy heavy-ion collisions and shed light on the properties of the produced medium.


Introduction
In the collision of two heavy ions moving relativistically, a high-density energetic state of matter is created in the overlap region of the two Lorentz-contracted nuclei. Earlier studies at the Relativistic Heavy-Ion Collider (RHIC), where gold nuclei were collided at nucleon-nucleon center-of-mass energies up to √ s NN = 200 GeV [1][2][3][4], found that the particles produced in rare, high-momentum-transfer scatterings encounter a dense medium with high stopping power for colored probes. The low-momentum particles that comprise the bulk of the medium exhibit strong azimuthal anisotropies that indicate a collective fluid expansion. These findings have been interpreted as manifestations of a strongly interacting quark-gluon plasma. The created medium is found to behave as a nearly perfect fluid [5][6][7][8] with a ratio of shear viscosity η to entropy density s approaching the conjectured lower limit of η/s ≥h/(4πk B ), found by considering the implications of the Heisenberg uncertainty principle for a viscous plasma [9]. Pressure gradients that develop in the fluid during the collision result in an anisotropic momentum distribution of the outflowing matter, which, in turn, leads to a preferential emission of particles in the short direction of the lenticular-shaped overlap region [10][11][12]. The hydrodynamic behavior suggests that local thermal equilibrium may be achieved very rapidly in the hot medium, with the observed anisotropy in particle emission therefore being sensitive to the basic properties of the created system, such as the equation of state, the η/s value, and the speed of sound in the medium. The anisotropy also depends on the initial conditions, allowing the investigation of whether a Glauber-like picture of individual nucleon collisions [13] prevails or if gluon saturation effects, as found in the color glass condensate (CGC) model [14,15], play an important role.
More recently, at the Large Hadron Collider, the azimuthal anisotropy measurements have been extended to a much higher energy, with PbPb collisions at √ s NN = 2.76 TeV [16][17][18]. Moreover, the azimuthal distributions are being investigated with greater precision with the exploration of higher-order anisotropies at both the RHIC [19][20][21][22][23] and LHC [24][25][26] facilities. The azimuthal dependence of the particle yield N can be written in terms of harmonic expansion coefficients v n ; with [27] E d 3 N dp 3 = 1 2π d 2 N p T dp T dy 1 + where φ, E, y, and p T are the azimuthal angle, energy, rapidity, and transverse momentum of the particles, respectively. In analyzing an experimental distribution, the reference angle Ψ needs to be determined empirically and typically corresponds to the direction of the greatest azimuthal density of outgoing particles. The higher-order harmonics are expected to be particularly sensitive to fluctuations in the initial conditions [28][29][30][31][32][33][34][35][36][37][38][39][40] and to the shear viscosity of the created medium [40][41][42][43][44].
This paper presents results from the Compact Muon Solenoid (CMS) Collaboration on higherorder harmonic anisotropy components for PbPb collisions at √ s NN = 2.76 TeV using the event plane [27], multiparticle cumulant [45], and Lee-Yang zeros [46,47] methods to exploit the different sensitivities of these methods to initial state fluctuations. In addition to correlations resulting from flow and initial state fluctuations, there exist other sources of azimuthal correlations, such as those arising from resonance decays and jets. These "nonflow" correlations either do not or only weakly depend on the bulk motion of the medium and need to be considered when determining the "true," global collective flow behavior. The methods discussed in this paper are affected differently by the nonflow effects and steps are taken to minimize the influence of these processes when possible. This work extends the previously published CMS results on elliptic flow (the n = 2 harmonic) [18]. The data and event selection used here are identical to those employed in the elliptic-flow analysis, and the current discussion of the experimental methods summarizes a more extensive discussion found in the earlier paper. New results are presented for harmonics n = 3, 4, 5, and 6 of charged particle distributions as a function of their transverse momentum (0.3 < p T < 8.0 GeV/c), collision centrality (0-70%), and pseudorapidity (|η| < 2.0). The pseudorapidity is defined in terms of the polar angle θ with η = − ln[tan(θ/2)]. The collision centrality reflects the degree of overlap of the two colliding nuclei, with 0% central events corresponding to impact parameter (i.e., the distance between the centers of the two colliding nuclei at closest approach) b = 0. Some of the earlier CMS elliptic-flow results are included in order to present a more complete description of the anisotropy behavior. The CMS Collaboration results on higher-order harmonic anisotropies obtained using the two particle correlation technique [48] are also included here for comparison.
The paper is organized as follows. Section 2 presents an overview of the experimental procedures and different methods that are used in the analysis. Various sources of systematic uncertainties are discussed. This section also develops the Glauber model eccentricities used in discussing the experimental results. Section 3 presents the "differential" and yield-weighted average harmonic coefficients ("integral" flow) for the different methods. The pseudorapidity dependence is presented for azimuthal anisotropy coefficients evaluated by the event plane method. Section 3 also contains a comparison of the new CMS results to previously published results of the ALICE and ATLAS Collaborations, as well as lower-energy results obtained by the PHENIX Collaboration at RHIC. Section 4 presents a discussion of the results and Section 5 gives an overall summary.

Experimental Details
The measurements were done with the CMS detector using data obtained in √ s NN = 2.76 GeV PbPb collisions during the 2010 heavy-ion run at the LHC. The analysis uses the same data and techniques as for the elliptic-flow study of Ref. [18], allowing for a direct comparison with the results of that study.

Experimental setup
The CMS detector consists of a silicon tracker, a crystal electromagnetic calorimeter, and a brass/scintillator hadronic calorimeter housed within a superconducting solenoid 6 m in diameter that provides a 3.8 T magnetic field. Muons are identified in gas-ionization chambers that are embedded in a steel flux return yoke. The muon information is not used in the current analysis. The CMS detector includes extensive charged particle tracking and forward calorimetry. The inner tracker, consisting of silicon pixel and strip detector modules, reconstructs the trajectories of charged particles within the pseudorapidity range |η| < 2.4. In the forward region, two steel/quartz-fiber Cherenkov hadron forward (HF) calorimeters cover a pseudorapidity range of 2.9 < |η| < 5.2. These calorimeters are azimuthally subdivided into 20 • modular wedges and further segmented to form 0.175 × 0.175 (∆η × ∆φ) "towers," where the angle φ is in radians. A set of scintillator tiles, the beam scintillator counters (BSC), are mounted on the inner side of the HF calorimeters and are used for triggering and beam-halo rejection. The BSCs cover the range 3.23 < |η| < 4.65. CMS uses a right-handed coordinate system where the x, y, and z axes are aligned with the radius of the LHC ring, the vertical direction, and the counterclockwise-beam direction, respectively, with the origin located at the center of the nominal interaction region. A more detailed description of the CMS detector can be found elsewhere [49]. In this analysis, the azimuthal anisotropy correlations are determined based on the charged particles reconstructed from the tracks in the silicon tracker. The event plane analysis also uses information from the HF calorimeters to establish event planes that have a large separation in pseudorapidity from the tracks used to determine the anisotropy harmonics.
Track reconstruction is accomplished by starting with a set of three reconstructed clusters in the inner layers of the silicon pixel detector that are compatible with a helical trajectory, have a p T above a minimum value, and are within a selected region around the reconstructed primary collision vertex, as determined through an iterative process. In determining tracks, trajectories with a minimum p T value of 0.9 GeV/c are first propagated outward through sequential silicon strip layers using a combinatorial Kalman filter algorithm [50]. Then, trajectories in the range of 0.2 < p T < 1.8 GeV/c are determined using only signals in the pixel detector, without requiring continuation of the trajectories to the silicon strip detector layers. The reconstructed tracks from both procedures are then merged, removing duplicate tracks by giving preference to the tracks extending into the silicon strip layers. A track will be "misreconstructed" if it is generated by the passage of more than one charged particle or if other spurious signals in the pixel or strip detectors contribute to its determination. The effect of misreconstructructed tracks is considered in determining the systematic uncertainties for the different methods, as discussed below.

Event and track selection
Minimum bias PbPb events are triggered by coincident signals from both ends of the CMS detector in either the BSC or HF detectors. This trigger is required to be in coincidence with the presence of both colliding ion bunches in the interaction region. Additional offline event selections are applied in order to obtain a pure sample of inelastic hadronic collisions, thus removing contamination from noncollision beam backgrounds and from ultra-peripheral collisions (UPC) where an electromagnetic interaction leads to the breakup of one or both Pb nuclei [51]. These offline selections include the requirement of proper timing of the BSC signals from both sides of the detector, a coincidence of at least three HF towers on each side of the interaction point, with at least 3 GeV energy deposited in each tower, a reconstructed vertex compatible with the expected collision region, and the shape of reconstructed clusters from the pixel detector being consistent with having been produced by particles originating from the primary collision vertex. These selections are described in more detail in Ref. [52].
Events used in this analysis are required to have a longitudinal vertex position within 10 cm of the nominal interaction point of the detector in order to consistently measure charged particle distributions over the tracker rapidity range. After all selections, 22.6 million events remain in the final sample, corresponding to an integrated luminosity of approximately 3 µb −1 . This final sample is the same data set used in the elliptic-flow study and is described in further detail in Ref. [18].

Centrality
The centrality of a collision is a measure of the degree of overlap of the colliding ions. Several quantities depend on the centrality and can be used for its determination. In this analysis, the centrality is based on the total energy deposited in both HF calorimeters, with the distribution of the total energy for all events divided into 40 centrality bins, each representing 2.5% of the total PbPb interaction cross section. More central events will have a larger particle multiplicity and therefore greater total energy in the HF calorimeters. Events falling into adjacent centrality bins are then combined to form the 5% and 10% bins used to present the final results. The measured charged particle multiplicity distribution does not represent the full interaction cross sec-tion because of inefficiencies in the minimum bias trigger and the event selection. Monte Carlo simulations (MC) are used to estimate the multiplicity distribution of charged particles in the regions where events are lost because of low trigger efficiency to correct for the inefficiencies. Comparing the simulated to the measured distribution, we determine the combined efficiency for the minimum bias trigger and the event selection to be (97 ± 3)%. Further discussion of the CMS centrality determination can be found in Ref. [53].

Analysis methods
To extract the v n coefficients, the event plane [27], the generating-function-based multiparticle cumulant [45], and the Lee-Yang zeros [46,47] methods are used, as described in the next three subsections. An earlier CMS paper [48] studied the higher harmonic coefficients using the two particle correlation method. For completeness, these earlier results are also presented. The two particle correlation method is similar to a two particle cumulant analysis although, in this case, with the requirement of a gap in pseudorapdity between the particles in a pair.

Event Plane method
The event plane method [27] measures the azimuthal anisotropy with respect to an event plane angle of a given order m that is determined in a different pseudorapidity region from that where the anisotropy coefficient v n is being measured. Expressed in terms of the event plane angle Ψ m corresponding to harmonic m, the azimuthal distribution can be written as [27] dN where dN is the number of charged particles emitted into a differential azimuthal angular range d(φ − Ψ m ), n = km, and the dependence of v obs n on the event plane harmonic m is explicitly noted. The event plane harmonic m can take any integer value greater than 0. Generally, when higher-order event planes are considered (m > 2), only the k = 1, n = m term is found to be needed to describe the corresponding azimuthal distribution. For m > 1, it is not possible to describe the correlations corresponding to n < m. It has recently been noted [54] that participant fluctuations can lead to terms where n is not an integer multiple of m. These mixed harmonic terms will not be studied in this paper.
The method assumes that the event plane angle is a pseudorapidity independent global observable. For this analysis, event plane angles are calculated using the measurements of transverse energy obtained with elements of the azimuthally symmetric HF calorimeters, using where φ i is the azimuthal position of the ith-tower element. The bracket indicates a sum over tower elements. The transverse energy in each tower i is used as the weight w i . For each fundamental harmonic m we define two event planes Ψ m (HF−) (−5 < η < −3) and Ψ m (HF+) (3 < η < 5), corresponding to the HF calorimeters on either side of the CMS detector. A standard event plane flattening procedure is employed to avoid having an azimuthal bias introduced by detector effects [18,27]. The "flattened" event plane distributions show no significant azimuthal anisotropy in any harmonic order m when plotted with respect to the laboratory frame.
The differential-anisotropy parameters v n (p T , η) are then determined with v obs n (p T , and v obs n (p T , where In Eqs. (4) and (5) the first sum is over all particles i with azimuthal angles φ i,j within a given pseudorapidity range in an event j with a given Ψ m and then a sum is taken over all events. Particles with η < 0 are correlated with HF+, and those with η > 0 are correlated with HF−. In this manner, a minimum gap of 3 units in pseudorapidity is maintained between the particles used in the event plane angle determination and those used to determine the azimuthal anisotropy harmonic, helping in the suppression of short-range nonflow effects.
The observed harmonic coefficient v obs n depends on the resolution of the event plane angles and is therefore sensitive to both the particle multiplicity and the magnitude of the azimuthal asymmetry in the pseudorapidity range used to determine the event plane angle. The final anisotropy values are obtained by correcting for the event plane angle resolution using a cor- To determine R n {Ψ m } we use the three-subevent method [27]. A subevent can be defined by restricting the particles used in an event plane determination to those found in a particular pseudorapidity range. The corresponding subevent angle is calculated only using particles in this limited region of pseudorapidity. The resolution of Ψ a m associated with subevent a (e.g., HF−) is determined using additional separate subevent angles Ψ b m (e.g., HF+) and Ψ c m , with For this analysis Ψ c m is determined using the silicon tracker detector. Charged particles with |η| < 0.75 are used to determine the Ψ c m subevent angle based on Eq. (3). In this case the weight w i of particle i is taken as the transverse momentum of the particle. The resolution correction values used in the analysis are shown in Fig. 1. For the second-order event plane (Ψ 2 ; m = 2) the resolution corrections are shown for the n = 2, 4, and 6 harmonic terms of Eq. (2). In general, the symmetry of the HF+ and HF− detectors leads to very similar R n {Ψ m } values for the event planes corresponding to the two detectors. However, when the resolution correction factors become sufficiently small, statistical fluctuations start to have a significant influence on the observed values, resulting in differences between the HF+ and HF− values. In Fig. 1 this is evident for R n {Ψ m } values less than 0.1. The contribution of the resolution correction to the overall v n systematic uncertainty is based on the observed difference in the HF+ and HF− results.
At low p T (< 0.8 GeV/c) the fraction of misreconstructed tracks is significant, reaching levels of up to 5% at mid-rapidity and up to 25% at forward rapidity (|η| ≈ 2) for the most central events. In this kinematic region, the v n signal is small. Full CMS Monte Carlo simulations based on the HYDJET [55] and AMPT [56] event generators indicate that the v n signal from misreconstructed tracks is approximately constant in this low-p T region with a value that can be larger than that of the properly reconstructed tracks. These studies suggest that the v n values from misreconstructed tracks can be characterized by a scaling factor κ, with v mis n = κ v n , and where the efficiency-and yield-weighted average v n is performed over the range 0.3 < p T < 3 GeV/c. Letting f represent the proportion of misreconstructed tracks, the observed v obs n value is related to the "true" flow signal v real n and that of the misreconstructed tracks v mis n by v obs Centrality (%)  The scaling factor κ is set by finding the value that leads to the least sensitivity of the final results when varying the track criteria and hence the number of misreconstructed tracks, as determined by the Monte Carlo simulations. The value is dependent on the harmonic and found to be κ = 1.3 ± 0.1 for v 2 , 1.0 ± 0.4 for v 3 , and 0.8 ± 0.6 for all higher harmonics, where the uncertainties are based on the observed sensitivity of the final results on the track criteria. The correction for misreconstructed tracks is negligible for p T values above ≈0.8 GeV/c.

Cumulant method
The cumulant method measures flow based on a cumulant expansion of multiparticle azimuthal correlations. An "integral" flow, or reference flow, of order m is first determined using an integral generating function of the multiparticle correlations in a complex (imaginary) plane [45]. The reference flow plays a similar role to that of the event plane determination in the event plane method. The integral generating function is constructed using all particles in a broad (p T , η) window, averaging over the events in a given centrality class. Then, the differential flow, i.e., the flow in a narrower phase space window, is measured with respect to the reference flow. A particle in the differential bin is correlated to the particles used for the reference flow through a differential generating function. Expressing the differential flow in terms of p T bins, a reference flow of order m can be used to determine differential flow of order n, where n is an integral multiple of m. To avoid autocorrelations, if a given particle is used in determining the differential flow, the particle will be excluded in the calculation of the reference flow. The generating functions for the reference and for differential flows are calculated at several different points in the complex plane and then interpolated. The interpolation points are defined by z j,k = r 0 j exp(i2kπ/8), where r 0 j is the radius of a point and 2kπ/8 its polar angle. In this analysis, the counting indices are taken as j = 1, 2, 3 and k = 0-7. Two (three) values for the radius parameter r 0 are used for the differential (reference) flow. The radius parameters are determined according to the detected charged particle multiplicity and the number of events analyzed in each centrality class. To reduce the effect of fluctuations in the event-by-event multiplicity on the flow measurements, for each given centrality interval, 80% of the mean multiplicity M of particles for that interval are selected at random in each event to determine the reference flow. In addition, the transverse momentum is restricted to p T < 3 GeV/c in the reference flow to reduce the nonflow contribution arising mainly from jets. The cumulant v 3 coefficient is measured with four particle correlations and denoted v 3 {4}, with the differential and integral flow both of order m = n = 3. The use of a four particle correlation strongly reduces the influence of the nonflow effects that are evident in two particle cumulant analyses [18]. The cumulant v 4 coefficient is calculated relative to the integral v 2 behavior using a five particle correlation and is denoted v 4 {5}. It is not possible to extract flow coefficients with n > 4 using this method since the signals become too small. However, v 6 coefficients calculated relative to the integral v 2 behavior are obtained using the Lee-Yang zeros method, as discussed below.
If the detector had uniform acceptance and full efficiency, a reference flow value calculated for a given pseudorapidity range would be equivalent to the yield-weighted flow coefficient for the same phase space. The need for an efficiency correction for the differential-flow coefficients is avoided by choosing sufficiently small p T bins such that the efficiency does not change significantly across the bin. However, this is not true for the reference flow and, to account for efficiency and acceptance effects, the yield-weighted cumulant v n values presented in this paper are based on the corresponding differential flow coefficient v n (p T ) weighted by the efficiency and acceptance corrected yields, similar to what is done for the event plane analysis. For both the multiparticle cumulant method and the Lee-Yang zeros method discussed in the next subsection, the influence of misreconstructed tracks at low-p T values is found to have only a small effect on the final results and is included as part of the systematic uncertainty.

Lee-Yang zeros method
The Lee-Yang zeros (LYZ) method [46,47] studies directly the large-order behavior of the cumulant expansion of the azimuthal correlations. As done for the cumulant analysis, a reference flow of order m is first determined based on a complex generating function that is calculated over a large range of momentum and pseudorapidity for a given centrality range. The generating function can be expressed in terms of either a sum or a product of individual particle terms. The product form is used in this analysis since it is expected to be less sensitive to nonflow and autocorrelation effects [47]. An estimate of the reference flow is found in terms of the location of the first minimum of the generating function calculated for a fixed projection angle. The analysis uses five different projection angles and averages the results to reduce statistical uncertainties. Charged particles with 0.3 < p T < 12 GeV/c and |η| < 2.4 are used to calculate the reference flow, which corresponds to the yield-weighted average flow in the indicated phase space, but neglects efficiency and acceptance effects.
Once the reference flow has been established, the differential flow v n {LYZ}(p T ) in a limited p T and pseudorapidity range can be determined with respect to the generating function of the reference flow. Again, the differential flow harmonic n can be any integral multiple of the reference flow harmonic m. As with the cumulant analysis, the integral v n {LYZ} coefficients presented in this paper are obtained by taking a yield-weighted average of the differential flow results, after correcting the yield for efficiency and acceptance effects. The new CMS Lee-Yang zeros results are for the n = 4 and 6 harmonics based on the m = 2 order reference flow. Measurements of the n = 3 and 5 harmonics are not possible because of their small magnitudes. Details of this part of the analysis can be found in Ref. [18].

Systematic uncertainties
The systematic uncertainties include those common to all methods, as well as method specific ones. Tables of uncertainties for the results presented in this paper are given in the Appendix. Common uncertainties include those resulting from the tracking efficiency and the centrality determination. Protons, pions, and kaons can have different v n values, and therefore species dependent differences in their tracking efficiency will affect the unidentified, charged particle results. The influence of particle composition is studied by determining the efficiency for identified particle detection using simulations of the CMS detector and then making different assumptions of the p T dependence of the particle mix based on the HYDJET event generator [55] and on assuming a similar behavior as found at RHIC [57]. As previously done to account for this effect in the v 2 measurement [18], a conservative 0.5% uncertainty, independent of p T , η, and centrality, is assumed. The sensitivity of the harmonic-flow coefficients to the centrality calibration is evaluated by varying the trigger efficiency by ±3%. The resulting uncertainty in v n is of the order of 1% and this value is adopted independent of p T and centrality. The uncertainty in the overall charged particle efficiency corrections, which only affects the yieldweighted average v n values, is evaluated by determining the efficiency based on the HYDJET model and, separately, by embedding simulated pions into recorded PbPb collision event data. The difference between the two resulting efficiencies gives at most a 0.5% change in the yieldweighted average v n values, which is taken as a systematic uncertainty.
Misreconstructed tracks affect both the differential v n (p T ) results and the yield-weighted average v n measurements of all methods, although not necessarily to the same extent. Therefore, separate studies are performed for each method. The effect of misreconstructed tracks, evaluated by varying the track quality criteria and labeled as the "Track Quality Requirements" uncertainty in the Tables 3 to 12, generally accounts for the largest single contribution to the systematic uncertainty, especially at low p T and for the most central events. Different sets of track quality requirements on the pixel tracks are used, and the ratio of the results provide an estimate of the systematic uncertainty from this source in different p T ranges. Track quality requirements include having the track pointing back to within a specified range of the reconstructed vertex location and having a specified goodness-of-fit for the track.
For the event plane method, the uncertainty in the resolution correction value is primarily a consequence of its statistical uncertainty and is generally small compared to the track quality requirement uncertainty. This is seen in Tables 3 to 8 where the systematic uncertainties for the v n (p T ) values obtained using the event plane method are presented. However, for the v 5 and v 6 results the resolution correction uncertainty becomes comparable to that for the track quality requirement uncertainty for mid-central events. The various systematic uncertainties are taken to be uncorrelated and added in quadrature in the measurements of the v n coefficients.
In addition to the systematic uncertainty terms common to all methods, the cumulant analyses are also influenced by the choice of the r 0 radius parameter and by the effect of fluctuations in the event-by-event multiplicity on the reference flow. These uncertainties are estimated by varying the r 0 parameter and comparing the flow results with and without the selection of 80% of the mean multiplicity. Tables 9 and 10 show the systematic uncertainties associated with the v 3 {4} and v 4 {5} results. The p T dependence of v 3 {4} could only be measured up to 4 GeV/c because of the small amplitude of this coefficient.
The effect of fluctuations in the event-by-event multiplicity is studied for the LYZ method by analyzing the events in finer 2.5% centrality bins. Tables 11 and 12 show the systematic uncertainties in the v 4 {LYZ} and v 6 {LYZ} results, respectively. In these cases, the total uncertainties, again found by adding the component uncertainties in quadrature, are dominated by the track quality requirement uncertainties.
The systematic uncertainties for the yield-weighted, average v n values with 0.3 < p T < 3.0 GeV/c are also shown in Tables 3 to 12. This p T range is the same as used in the previous CMS ellipticflow analysis [18] and covers the p T range dominated by hydrodynamic flow.

Glauber model eccentricity calculations
The Glauber model treats a nucleus-nucleus collision as a sequence of independent nucleonnucleon collisions (see Ref. [13] and references therein). The model can be used to obtain anisotropy parameters based on the transverse location of the participant nucleons [58]. These, in turn, are expected to be reflected in the observed particle anisotropies. The Glauber model assumes that the nucleons in a nucleus are distributed according to a Woods-Saxon density where ρ 0 is the nucleon density in the center of the nucleus, R is the nuclear radius, a is the skin depth, and w characterizes deviations from a spherical shape. For 208 Pb, the parameters R = 6.62 ± 0.13 fm, a = 0.546 ± 0.055 fm, and w = 0 are used [59]. The model assumes that nucleons in each nucleus travel on straight-line trajectories through the colliding system and interact according to the inelastic nucleon-nucleon cross section, σ NN inel , as measured in pp collisions. A value of σ NN inel = 64 ± 5 mb, which is found through an interpolation of values obtained at different center-of-mass energies [60], is used in the calculations at √ s NN = 2.76 TeV.
The spatial anisotropies of order n based on a participant plane angle of order m are calculated using the transverse location of each participant [37], with where, for a participant located at coordinates {x, y} in the transverse plane, r n ⊥ = x 2 + y 2 n , and φ = arctan (y/x). The "participant plane" angle Φ m is then found by summing over all participant particles with For n = m we define n = n,n . With this definition, n (or n,m ) can only take positive values and represents the maximum asymmetry for each collision. It has been demonstrated that a common behavior is achieved for the elliptic-flow coefficient v 2 in AuAu and CuCu collisions at √ s NN = 200 GeV when scaled by the eccentricity n and plotted as a function of the transverse charged particle areal density [29]. This scaling persist at LHC energies [18]. Table 1 lists the results for the average number of participants, N part , and the root-meansquared evaluation of the eccentricities, 2 n or 2 n,m , and the corresponding systematic uncertainties for the centrality bins used in this analysis. The method used to convert from impact parameter to centrality is discussed in Ref. [52]. The uncertainties in the parameters involved in the Glauber model calculations contribute to the systematic uncertainty in N part and n for a given centrality bin.

Results
This section presents the results for the higher-order harmonic coefficients. Previously published two particle correlation results [48] from the CMS Collaboration are also presented for completeness. The p T dependence of the coefficients at mid-rapidity is presented first, comparing the values obtained using the different analysis methods. This is followed by the measurements of the yield-weighted average v n values, which are given in terms of both their centrality and pseudorapidity dependencies. We conclude the section with comparisons of the CMS measurements to published results of the ALICE [24,25] and ATLAS [26] Collaborations. Figure 2 shows the v 3 coefficient results for |η| < 0.8 based on the event plane v 3 {Ψ 3 } and four particle cumulant v 3 {4} methods. The analyses employ the same event selection as used in the previous elliptic-flow (n = 2) study of Ref. [18]. Also shown are the two particle correlation results v 3 {2, |∆η| > 2} of Ref. [48]. The two particle correlation method is similar to a two particle cumulant analysis, as used in Ref. [18], but requires a pseudorapidity gap between the two particles, which was not the case for the two particle cumulant analysis. For the two particle correlation method, charged particles with |η| < 2.5 are first selected. To reduce the effect of nonflow processes, particle pairs are chosen with the requirement of a pseudorapidity gap between the particles in each pair of 2 < |∆η| < 4. This method, as applied to LHC data, is described in detail in Refs. [24,48].

The p T dependence of v n at mid-rapidity
The event plane and two particle correlation results are found to be very similar, although the results from the two particle correlations have systematically smaller values. This suggests similar sensitivity to initial state fluctuations and nonflow effects for the current implementations of the two methods when a large pseudorapidity gap is required for both analyses. We also note that the values of the harmonic coefficients determined from the two particle correlation method correspond to v 2 n . For near-perfect event plane resolution with R n ≈ 1 (Eq. (6)), the event plane results are expected to approach v n , whereas for lower values of R n , the event plane method gives values closer to v 2 n [29]. We discuss this later in the paper. With the R n>2 values shown in Fig. 1, the event plane method is expected to produce values very close to v 2 n in the absence of fluctuations. The smaller values found for the two particle correlation method is likely a consequence of sampling a larger pseudorapidity range than the event plane analysis, which only considered particles with |η| < 0.8. As shown later in this paper, the spectrum-weighted harmonic coefficients are found to have their maximum values near η ≈ 0 (see Section 3.2). The pseudorapidity gap used in the two particle correlations prevents the selection of both particles in a pair from the mid-rapidity, maximum v 3 region, thus assuring a somewhat smaller v 3 result as compared to the event plane analysis. The significantly smaller v 3 {4} values for the four particle cumulant method can be attributed to the difference in how fluctuations affect two particle and higher-order correlations [36]. Assuming a smooth overlap region and in the absence of fluctuations, the v 3 harmonic is expected to vanish based on the symmetry of the overlap region. For the two particle correlation and event plane results the fluctuations are expected to increase the harmonic coefficients with respect to those expected without initial state fluctuations, whereas for fourth-(and higher-) order particle correlations, the fluctuations can lower the values [36]. It was shown in Ref. [18] that effects of fluctuations can largely account for the differences observed in the v 2 values for the different methods. Fluctuations are expected to dominate the initial state geometry of the odd-harmonic, n = 3 asymmetry, as discussed in Refs. [30,37,38,61].   The event plane analysis (filled circles) is with |η| < 0.8. The two particle correlation results (open circles) are from a previous CMS measurement [48]. Statistical (error bars) and systematic (light gray boxes) uncertainties are shown. very similar results, with only a small dependence on centrality for most of the range.
Finally, Fig. 5 shows the event plane results for v 6 based on both the m = 2 and 6 event planes, as well as the LYZ results based on the m = 2 integral reference flow. In this case, the event plane results based on the second-order reference distribution are consistently smaller than those found for either the LYZ method, which are also based on a second-order integral-flow behavior, or the event plane method using a sixth-order reference distribution. The higher values and relatively weak centrality dependence found for the v 6 {Ψ 6 } results are consistent with these values being strongly influenced by fluctuations.
Summarizing the results for this section, the differential azimuthal harmonics are found to have their strongest dependence on centrality when the "reference" particles are based on the second-order participant planes, as is the case for the A weaker centrality dependence is observed in the other cases where the higher-order (m > 2) reference plane is of the same order as the harmonic being studied. This weak centrality dependence suggests a reduced influence of the average geometry of the overlap region, as might be expected if fluctuations in the participant locations dominate the results.

Yield-weighted average anisotropies
The centrality dependence of the yield-weighted average v n values with 0.3 < p T < 3.0 GeV/c is shown in Fig. 6 for the different methods. This is the p T range for which a significant hydrodynamic-flow contribution is expected. For completeness, the earlier n = 2 results from Ref. [18] are also shown. As noted for the p T dependent results, a stronger centrality dependence is found for analyses based on the m = 2 reference plane. This is more clearly seen in Fig. 7 where the analyses based on a second order, m = 2 reference plane are shown in Fig. 7a and those based on higher order, m > 2 reference planes are shown in Fig. 7b. The v 6 {Ψ 6 } results do not show a centrality dependence, although the large statistical uncertainties may mask a weak dependence. Figure 8 shows the pseudorapidity dependence for the yield-weighted average event plane v n values with n = 2, 3, and 4, and with 0.3 < p T < 3.0 GeV/c. The values for the n = 5 and 6 harmonics are found to be too small to establish their dependence on pseudorapidity. The data are sorted into ten pseudorapidity bins of ∆η = 0.4 spanning the range −2.0 < η < 2.0. The distributions have their maximum values near mid-rapidity, with the fractional decrease from |η| = 0 to |η| = 2 being similar for the different centrality ranges in a given harmonic. The v 2 results are from Ref. [18] and included for completeness. Statistical (error bars) and systematic (light gray boxes) uncertainties are shown. The different results found for a given v n reflect the role of participant fluctuations and the variable sensitivity to them in each method, as discussed in the text.

Comparison with other results
The current results extend and largely confirm previous results published by the ALICE [24,25] and ATLAS [26] Collaborations on higher-order harmonic correlations. Representative comparisons of the CMS results with those of these other two collaborations are shown in Figs. 9 to 12.
Corresponding results by the PHENIX Collaboration for AuAu collisions at √ s NN = 200 GeV are also shown [23]. Differences in the centrality and pseudorapidity ranges chosen by the different collaborations need to be considered in comparing the results. Table 2 summarizes the experimental conditions for the different measurements. Figure 9 compares results of CMS, ATLAS, and ALICE for the p T dependent v 3 coefficient. The ATLAS results for v 3 {Ψ 3 } are consistently lower than the CMS results for all but the most peripheral centrality bin. This is expected based on the larger pseudorapidity range being used for the ATLAS measurement. Good agreement is seen between the two particle correlation results of the CMS and ALICE Collaborations. This suggests that the pseudorapidity gap of     |∆η| > 0.8 employed by ALICE is already sufficient to remove most of the dijet contribution to these correlations.

Method(s)
The comparison of v 4 and v 5 values found by the three experiments lead to similarly consistent results. The CMS and ATLAS results are also very similar for v 6 {Ψ 6 }(p T ) within the respective uncertainties, as seen in Fig. 12. Figures 9-11 also show the predictions of the IP-Glasma+MUSIC model of Ref. [62], as discussed in the next section. The model calculations cover the hydrodynamic-dominated region of the p T distributions.
The lower energy v 3 {Ψ 3 } and v 4 {Ψ 4 } results of the PHENIX Collaboration for AuAu collisions at √ s NN = 200 GeV are also shown in Figs. 9 and 10, respectively. The n = 3 AuAu results are systematically lower than those obtained by the higher-energy LHC measurements, consistent with what was previously observed for the elliptic-flow, n = 2 harmonic [18]. A different picture is suggested by the n = 4 distributions, where now the RHIC results are systematically higher than those observed at the LHC, although with large systematic uncertainties.   Table 2. The predictions of the IP-Glasma+MUSIC model [62] for PbPb collisions at √ s NN = 2.76 TeV are shown by the solid lines in the 0-5%, 10-20%, 20-30%, 30-40%, and 40-50% panels for 0 < p T < 2 GeV/c.

Discussion
There is considerable interest in how the spatial anisotropies, as characterized by spatial anisotropy parameters n , created early in the collision of two ultra-relativistic heavy ions get transformed  Table 2. The predictions of the IP-Glasma+MUSIC model [62] for PbPb collisions at √ s NN = 2.76 TeV are shown by the solid lines in the 0-5%, 10-20%, 20-30%, 30-40%, and 40-50% panels for 0 < p T < 2 GeV/c.  Table 2. The predictions of the IP-Glasma+MUSIC model [62] are shown by the solid lines in the 0-5%, 10-20%, 20-30%, 30-40%, and 40-50% panels for 0 < p T < 2 GeV/c.  Table 2. into the experimentally observed azimuthal anisotropy of emitted particles [30, 33-35, 37, 38, 61, 63-65]. The higher-order harmonics are expected to be more sensitive to the details of the collision geometry and its event-by-event fluctuations. This section develops the scaling behavior of the experimental v n coefficients in terms of the Glauber model n values and also explores the effect of fluctuations on the different analysis methods.
It is now recognized that the different experimental methods used in determining the v n coefficients are related differently to the underlying n values. For example, v n {Ψ n } coefficients obtained with near-unity values for the event plane resolution factor R are expected to scale with n , whereas these coefficients scale with The details of the eccentricity scaling are model dependent and beyond the scope of this paper. However, to achieve an overview of the geometry scaled behavior, we present in Fig. 13 the yield-weighted average v n results of Fig. 6 as a function of the Glauber model 2 n,m azimuthal asymmetries discussed in Section 2.6. In general, the v n coefficients are found to increase monotonically with the Glauber model eccentricities for the most central events, up to the maxima in the distributions shown in Fig. 6, although large uncertainties are affecting the n > 4 event plane results and some method differences are observed. For n = 2, both the event plane and four particle cumulant results show similar behavior for the most central events, with the overall magnitude of the v 2 {4} coefficients being smaller. The observed difference is consistent with the fourth-order cumulant results scaling with the four particle cumulant eccentricity 2 {4}, as shown in Ref. [18]. A much larger difference is observed between the event plane and cumulant results for n = 3, as would be expected if the odd harmonics are dominated by fluctuation effects, which are strongly suppressed in the multiparticle cumulant analysis. The higher-order harmonic event plane results with n = m show relatively weak scaling with the Glauber geometry, also suggesting significant fluctuation components. For n = 4 the harmonic components based on a second-order reference plane, as is the case for v 4 {5}, v 4 {Ψ 2 }, and v 4 {LYZ}, are found to have a much stronger dependence on the Glauber eccentricity for more central events than is evident for the analysis based on the fourth-order event plane. Figure 14 shows the eccentricity scaled v n {Ψ n } values as a function of the harmonic order n for five different p T ranges and for four different centrality ranges. For all but the most central events, the v n values are found to decrease with increasing harmonic number. The rate of this decrease is expected to be sensitive to the shear viscosity of the medium, which leads to greater damping of the higher-order harmonic anisotropies [42,43]. For the most central events, the scaled v 3 coefficient is found to become larger than v 2 for the highest p T bin of 3.5 < p T < 4.0 GeV/c. Overall, the observed fall-off with harmonic order is very regular. The more central events demonstrate a fall-off that is steeper than an exponential in n. For midcentral events, however, the fall-off appears to scale as an exponential in n. Recent papers have suggested that the higher-order harmonic components may reflect a strong nonlinear response, particularly for noncentral events, with the higher-order harmonics dependent on mixtures of lower-order eccentricities [54,66].
Event-by-event fluctuations in the location of the participant nucleons can have a large and method dependent influence on the harmonic coefficients [29,36]. Expressing the fluctuations in terms of the azimuthal anisotropy in the participant plane v, where the harmonic number is suppressed, the magnitude of the fluctuations is given by It can then be  Figure 13: (Color online) Yield-weighted average azimuthal asymmetry parameters v n , for n = 2-6 with 0.3 < p T < 3.0 GeV/c, as a function of the corresponding Glauber model rms anisotropy parameters < 2 n >. The CMS v 2 results are from Ref. [18] and included for completeness. The v 4 {Ψ 2 } and v 6 {Ψ 2 } results are plotted against < 2 4,2 > and < 2 6,2 >, respectively, as given in Table 1. Statistical (error bars) and systematic (light gray boxes) uncertainties are indicated.
shown [36] that to leading order in σ v , two and four particle correlations are affected differently, with and The event plane method leads to an intermediate value, with where α is a parameter between 1 and 2 that is determined empirically in terms of the event plane resolution factor R ( Fig. 1) [29]. Multiparticle correlations with greater than four particles are expected to give results similar to those of four particle correlations. For harmonics with n > 2, the event plane resolutions lead to v {EP} ≈ v {2} using the parameterization of α in terms of the event plane resolution R given in Ref. [29]. Based on the larger R values observed for the CMS mid-central events with n = 2, the ratio of v 2 {2} to v 2 {Ψ 2 } is expected to be about 1.02 [36]. Motivated by Eqs. (11) and (12), as well as the approximate equality of v n {2} and v n {Ψ n }, Fig. 15 shows the ratio for n = 2 and 3. The cumulant multiparticle correlation harmonics are indicated by v n {Cum}, using four particle correlations for both n = 2 and 3. If the magnitude of the fluctuations is relatively small compared to the corresponding harmonic anisotropy term, this ratio should approach σ v / v . The points shown as solid squares in Fig. 15 are obtained by scaling the observed n = 2 event plane results to the limiting value associated with poor event plane resolution (i.e., v 2 n ). The resulting elliptic-flow (n = 2) fluctuation fraction near 0.4 is similar to that observed at RHIC for AuAu collisions at √ s NN = 200 GeV [67,68], although the more recent STAR results [68] are systematically higher than those observed at the LHC. The relative fluctuations are much larger for the n = 3 harmonic than for the elliptic flow, n = 2 harmonic.
In a recent event-by-event fluctuation analysis by the ATLAS Collaboration [69], the σ v / v ratio was measured directly by unfolding the event-by-event v n distributions with the multiplicity dependence of the measurements. These results are shown in Fig. 15 for the n = 2 and 3 harmonics. The CMS and ATLAS n = 2 results are in good agreement except for the most peripheral bins, where the CMS results are higher. The ATLAS analysis for the n = 3 harmonic leads to a relatively constant value of σ v / v with N part of approximately 0.53. For this higher harmonic the leading-order assumption made for Eqs. (11) and (12)   Another ratio that has received considerable attention in characterizing the azimuthal asymmetry is v 4 /v 2 2 [20,22,34,43,63,[70][71][72], where the n=2 and 4 harmonics are determined with respect to the elliptic-flow event plane. It is now recognized that this ratio, with a value of 0.5 obtained through ideal hydrodynamics [71], is strongly affected by flow fluctuations and nonflow correlations. The comparisons of theory to the experimental results need to account for how the results of the different analysis methods relate to the event-by-event v n asymmetry [34]. Figure 16 shows the ratio v 4 {Ψ 2 }/v 2 2 {Ψ 2 } for two different p T ranges as a function of centrality. In both cases the ratio initially decreases, but then remains relatively constant for centralities greater than ≈ 20%. The ratio using the yield-weighted average v n values over the larger p T range of 0.3 < p T < 3.0 GeV/c is systematically larger than that found in the p T range of 1.2 < p T < 1.6 GeV/c. The AuAu results obtained at √ s NN = 200 GeV by the PHENIX Collaboration for the range 1.2 < p T < 1.6 GeV/c are also shown. The CMS results are systematically higher by about 10%.
It has been shown that if the harmonic coefficients reflect ideal hydrodynamics with additional participant fluctuation effects, then the expected ratio is given by [34] v where β depends on the event plane resolution parameter R. The dashed line in Fig. 16 shows this result using a smooth fit to the σ v / v behavior found in Fig. 15 and the subevent results shown in Fig. 6 of Ref. [34] for β. The ideal hydrodynamic picture underestimates the observed  The PHENIX results (stars) for 1.2 < p T < 1.6 GeV/c are also shown for AuAu collisions at √ s NN = 200 GeV [22]. The dotted line is based on the method presented in Ref. [34] that allows the expected v 4 /v 2 2 ratio for an event plane analysis to be calculated based on the ideal hydrodynamics limit of 0.5 and the observed relative fluctuation ratio σ/ v n , as discussed in the text. Statistical (error bars) and systematic (light gray boxes) uncertainties are indicated.
In general, the flow harmonics are related to the initial state eccentricities through proportionality constants that depend on medium properties. It has been suggested that a greater sensitivity to initial state conditions might be achieved by studying the ratios of azimuthal anisotropy coefficients based on correlations of different numbers of particles or of mixed order [61,73]. One such ratio based on correlations with different numbers of particles is given by 2v 4 [73]. Figure 17 shows the quantity 2 − v 4 3 {4}/v 4 3 {Ψ 3 } as a function of centrality, together with the corresponding CGC and Glauber model predictions [73] . The yield-weighted average flow coefficients with 0.3 < p T < 3.0 GeV/c are used in determining the results from this analysis. The event plane results for v 3 {Ψ 3 } correspond to relatively low values of the resolution correction factor (see Fig. 1) and, consequently, are expected to give similar results to v 3 {2}. Within the current experimental uncertainties, it is not possible to clearly state which model works best. However, the results do suggest that a better determination of this ratio could help establish whether the initial state geometry is better described in a Glauber or CGC picture.
The STAR Collaboration at RHIC has recently released results for third-harmonic flow of charged particles from AuAu collisions at √ s NN = 200 GeV [74] and has presented them in terms of the ratio plotted in Fig. 17. Their results are consistent with having a constant value near 2 for the full centrality range, although with relatively large uncertainties. The four particle cumulant suppresses Gaussian fluctuations, and it is suggested that the larger multiplicities achieved at LHC energies may make these higher-energy results more sensitive to such fluctuations [74].
In a recent Letter [62] the higher-order harmonic coefficients have been predicted based on a calculation that uses the impact parameter dependent Glasma (IP-Glasma) model [15,65] to determine the early time evolution, and then switches to a relativistic hydrodynamic description using a (3+1)-dimensional relativistic viscous hydrodynamic simulation (MUSIC) [41]. The IP-Glasma model includes not only the quantum fluctuations associated with the distribution of nucleons, as reflected in Glauber model calculations, but also fluctuations in the color charge distributions inside a nucleon. These color charge fluctuations result in smaller-scale structure in the initial energy density profile than would be present if only sources of nucleon dimensions are considered. The results are shown by the curves in Figs. 9-11. Good agreement is found with the observed v n {p T } behavior in the lower-p T ranges.

Summary
Results from the CMS Collaboration have been presented on higher-order harmonic anisotropies of charged particles for PbPb collisions at √ s NN = 2.76 TeV. The harmonic coefficients v n have been studied as a function of transverse momentum (0.3 < p T < 8.0 GeV/c), centrality (0-70%), and pseudorapidity (|η| < 2.0) using the event plane, cumulant, and Lee-Yang zeros methods. The event plane method results are obtained with a pseudorapidity gap of at least three units, with the event plane determined in the range 3 < |η| < 5, suppressing the contribution of nonflow effects.
Comparisons of the event plane results with those of the cumulant and Lee-Yang zeros analyses suggest a strong influence of initial state fluctuations on the azimuthal anisotropies. The weak centrality dependence found for the event plane results based on event planes of harmonic order greater than two is also consistent with the presence of a strong fluctuation component. The pseudorapidity dependence of the higher-order azimuthal anisotropy parameters based on the event plane method is similar to that observed for elliptic flow, with only a modest decrease from the mid-rapidity values out to the limits of the measurement at |η| = 2.0. The mid-rapidity values are compared to those obtained by the ALICE [24,25] and ATLAS [26] Collaborations, and found to be in good agreement. A comparison is also done with lower-energy AuAu measurements by the PHENIX Collaboration at √ s NN = 200 GeV, with only small differences found with the much higher energy LHC data.
The results obtained for v 3 are compared to predictions from both the CGC and Glauber models. Both of these initial state models are found to be consistent with the data. It is noted that a calculation that employs IP-Glasma-model initial conditions for the early time evolution, followed by a viscous hydrodynamic development of the plasma, is quite successful in reproducing the observed v n (p T ) results in the low-p T , flow-dominated region [62].
The measurements presented in this paper help to further establish the pattern of azimuthal particle emission at LHC energies. Recent theoretical investigations have significantly increased our understanding of the initial conditions and hydrodynamics that lead to the experimentally observed asymmetry patterns. However, further calculations are needed to fully explain the method dependent differences seen in the data for the anisotropy harmonics. These differences can be attributed to the role of fluctuations in the participant geometry. Understanding the role of these fluctuations is necessary in order to establish the initial state of the created medium, thereby allowing for an improved determination of its properties. The current results are directly applicable to the study of the initial spatial anisotropy, time development, and shear viscosity of the medium formed in ultra-relativistic heavy-ion collisions.

A Systematic uncertainties
The systematic uncertainties for the results presented in this paper are given in Tables 3 to 12.   Table 3: Systematic uncertainties in the v 3 {Ψ 3 } values as a function of centrality in percent. Common uncertainties are shown at the top of the table, followed by those specific to the differential (p T dependent) and integral (|η| dependent) measurements.

Source
Centrality 0-10% 10-50% 50-70% Particle 0.   Table 6: Systematic uncertainties in the v 5 {Ψ 5 } values as a function of centrality in percent. Common uncertainties are shown at the top of the table, followed by those specific to the differential (p T dependent) and integral (|η| dependent) measurements.    Table 9: Systematic uncertainties in the v 3 {4} values as a function of centrality in percent. Common uncertainties are shown at the top of the table, followed by those specific to the differential (p T dependent) and integral (|η| dependent) measurements.  Common uncertainties are shown at the top of the table, followed by those specific to the differential (p T dependent) and integral (|η| dependent) measurements.  Table 11: Systematic uncertainties in the v 4 {LYZ} values as a function of centrality in percent. Common uncertainties are shown at the top of the table, followed by those specific to the differential (p T dependent) and integral (|η| dependent) measurements.