Mixed higher-order anisotropic flow and nonlinear response coefficients of charged particles in $\mathrm {PbPb}$ collisions at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76$ and 5.02$\,\text {TeV}$

Abstract

Anisotropies in the initial energy density distribution of the quark-gluon plasma created in high energy heavy ion collisions lead to anisotropies in the azimuthal distributions of the final-state particles known as collective anisotropic flow. Fourier harmonic decomposition is used to quantify these anisotropies. The higher-order harmonics can be induced by the same order anisotropies (linear response) or by the combined influence of several lower order anisotropies (nonlinear response) in the initial state. The mixed higher-order anisotropic flow and nonlinear response coefficients of charged particles are measured as functions of transverse momentum and centrality in $\mathrm {PbPb}$ collisions at nucleon-nucleon center-of-mass energies $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76$ and 5.02$\,\text {TeV}$ with the CMS detector. The results are compared with viscous hydrodynamic calculations using several different initial conditions, as well as microscopic transport model calculations. None of the models provides a simultaneous description of the mixed higher-order flow harmonics and nonlinear response coefficients.

1 Introduction

The azimuthal anisotropy of particle production in a heavy ion collision can be characterized by the Fourier expansion of the particle azimuthal angle distribution [1],

$$\begin{aligned} \frac{dN}{d\phi }=\frac{N}{2\pi }\sum _{n=-\infty }^{+\infty }V_n e^{-in\phi }, \end{aligned}$$

(1)

where $V_n=v_n\exp (in\varPsi _n)$ is the nth complex anisotropic flow coefficient [2]. The $v_n$ and $\varPsi _n$ are the magnitude and phase (also known as the nth order symmetry plane angle) of $V_n$, respectively. Anisotropic flow plays a major role in probing the properties of the produced medium in heavy ion collisions at the BNL RHIC [3,4,5,6] and CERN LHC [7,8,9]. Studies of flow harmonics higher than the second order [10,11,12], flow fluctuations [13,14,15,16], the correlation between the magnitude and phase of different harmonics [17,18,19,20,21,22,23,24], and the transverse momentum ($p_{\mathrm {T}}$) and pseudorapidity $(\eta )$ dependence of symmetry plane angles [25, 26], have led to a broader and deeper understanding of the initial conditions [3, 27] and the properties of the produced hot and dense matter. There are significant correlations between the symmetry plane angles of different orders [20], which indicate that higher-order mixed harmonics can be studied with respect to multiple lower-order symmetry plane angles.

In hydrodynamical models describing the quark-gluon plasma (QGP) created in relativistic heavy ion collisions, anisotropic flow arises from the evolution of the medium in the presence of an anisotropy in the initial-state energy density, as characterized by the eccentricities $\epsilon _n$ [10]. The magnitudes of the second- and third-order harmonic final state coefficients, $v_{2}$ and $v_{3}$, are to a good approximation linearly proportional to the initial-state anisotropies, $\epsilon _{2}$ and $\epsilon _{3}$, respectively [10, 17]. In contrast, $V_4$ and higher harmonics can arise from initial-state anisotropies in the same-order harmonic (linear response) or can be induced by lower-order harmonics (nonlinear response) [1, 28, 29]. More specifically, these harmonics can be decomposed into linear and nonlinear response contributions as follows [1, 28]:

$$\begin{aligned} \begin{aligned} V_4&= V_{4 L} + \chi _{422} V_2^2,\\ V_5&= V_{5 L} + \chi _{523} V_2 V_3,\\ V_6&= V_{6 L} + \chi _{624} V_2 V_{4L} + \chi _{633} V_3^2 + \chi _{6222} V_2^3,\\ V_7&= V_{7 L} + \chi _{725} V_2 V_{5L} + \chi _{734} V_3 V_{4L} + \chi _{7223} V_2^2 V_3, \end{aligned} \end{aligned}$$

(2)

where $V_{nL}$ denotes the part of $V_n$ that is not induced by lower-order harmonics [29,30,31], and the $\chi$ are the nonlinear response coefficients. Each nonlinear response coefficient has its associated mixed harmonic, which is $V_n$ measured with respect to the lower-order symmetry plane angle or angles. The strength of each nonlinear response coefficient determines the magnitude of its associated mixed harmonic. The $V_1$ terms are neglected in the decomposition in Eq. (2) because the correlation between $V_{n}$ and $V_{1}V_{n-1}$ was shown to be negligible after correcting $V_{1}$ for global momentum conservation [28]. This analysis focuses on the terms that only involve the two largest anisotropic flow coefficients $V_2$ and $V_3$ on the right-hand side of Eq. (2). The procedures used to extract both mixed-harmonic and nonlinear response coefficients are given in Sect. 4.

It is difficult to use measured $v_2$ and $v_3$ coefficients to evaluate hydrodynamic theories because these flow observables have a strong dependence on the initial anisotropies, which cannot be experimentally determined or tightly constrained. In contrast, most of the nonlinear response coefficients are not strongly sensitive to the initial anisotropies, which largely cancel in the dimensionless ratios used to determine these coefficients [1, 28, 31, 32]. As a result, their experimental values can serve as unique and robust probes of hydrodynamic behavior of the QGP [31].

Most previous flow measurements focused on $V_n$ (overall flow), i.e., $v_n$ with respect to $\varPsi _n$, which does not separate the linear and nonlinear parts of Eq. (2). Direct measurements of the mixed higher-order flow harmonics, $v_4$ and $v_6$ with respect to $\varPsi _2$, already exist at both RHIC [33] and LHC [11] energies, but were performed using the event plane method [34]. This method has been criticized for yielding an ambiguous measure lying somewhere between the event-averaged mean value $\left\langle v_n \right\rangle$ and the root-mean-square value $\sqrt{\left\langle {v_n^2} \right\rangle }$ of the $v_n$ distribution, depending on the resolution of the method [13, 16, 35]. This ambiguity can be removed by using the scalar-product method [35, 36], which always measures the root-mean-square values of $v_n$. The difference between the two methods is typically a few percent for $v_2$, $\sim$ 10% for $v_3$, and much larger for mixed harmonics [35].

This paper presents the mixed higher-order flow harmonics and nonlinear response coefficients for $n = 4$, 5, 6, and 7 using the scalar-product method. These variables are measured in $\mathrm {PbPb}$ collisions at nucleon-nucleon center-of-mass energies $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76$ and 5.02$\,\text {TeV}$, as functions of collision centrality and charged particle $p_{\mathrm {T}}$ in the region $|\eta |<0.8$. To compare the mixed flow harmonics with the overall flow coefficients, the higher-order flow harmonics with respect to the same-order symmetry plane, measured using the scalar-product method, are also presented.

2 The CMS detector

The central feature of the CMS apparatus is a superconducting solenoid of 6$\, \text {m}$ internal diameter, providing a nearly constant magnetic field of 3.8$\, \text {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. In this analysis, the tracker and the forward hadron (HF) calorimeter subsystems are of particular importance. The HF uses steel as an absorber and quartz fibers as the sensitive material. The two halves of the HF are located 11.2$\, \text {m}$ from the center of the interaction region, one on each end, and together they provide coverage in the range $3.0< |\eta | < 5.2$. These calorimeters are azimuthally subdivided into $20^{\circ }$ modular wedges and further segmented to form $0.175{\times }0.175$ $(\varDelta \eta {\times }\varDelta \phi )$ “towers”, where the angle $\phi$ is in radians. The silicon tracker measures charged particles within the range $|\eta | < 2.5$. It consists of 1440 silicon pixel and 15,148 silicon strip detector modules. For nonisolated particles of $1< p_{\mathrm {T}} < 10{\,\text {GeV/}c}$ and $|\eta | < 1.4$, the track resolutions are typically 1.5% in $p_{\mathrm {T}}$ and 25–90 (45–150)$\,\upmu \text {m}$ in the transverse (longitudinal) impact parameter [37]. The Beam Pick-up Timing for the eXperiments (BPTX) devices are located around the beam pipe at a distance of 175$\, \text {m}$ from the interaction region on both sides, and are designed to provide precise information on the LHC bunch structure and timing of the incoming beams. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [38]. The Monte Carlo simulation of the particle propagation and detector response is based on the Geant4 [39] program.

3 Event and track selections

This analysis is performed using minimum bias $\mathrm {PbPb}$ data collected with the CMS detector at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} =$ 5.02 and 2.76$\,\text {TeV}$ in 2015 and 2011, corresponding to integrated luminosities of 13$\,\mu \text {b}^{-1}$ and 3.9$\,\mu \text {b}^{-1}$, respectively. The minimum bias trigger [40] used in this analysis requires coincident signals in the HF calorimeters at both ends of the CMS detector with total energy deposits above a predefined energy threshold of approximately 1$\,\text {GeV}$ and the presence of both colliding bunches in the interaction region as determined using the BPTX. By requiring colliding bunches, events due to noise (e.g., cosmic rays and beam backgrounds) are largely suppressed. In the offline analysis, events are required to have at least one reconstructed primary vertex, which is chosen as the reconstructed vertex with the largest number of associated tracks. The primary vertex is formed by two or more associated tracks and is required to have a distance of less than 15$\, \text {cm}$ along the beam axis from the center of the nominal interaction region and less than 0.15$\, \text {cm}$ from the beam position in the transverse plane. An additional selection of hadronic collisions is applied by requiring at least three towers, each with total energy above 3$\,\text {GeV}$ in each of the two HF calorimeters. The average number of collisions per bunch crossing is less than 0.001 for the events used in this analysis, with a pileup fraction less than 0.05%, which has a negligible effect on the results. Events are classified using a centrality variable that is related to the degree of geometric overlap between the two colliding nuclei. Events with complete (no) overlap are denoted as centrality 0 (100)%, where the number is the fraction of events in a given class with respect to the total number of inelastic hadronic collisions. The centrality is determined offline via the sum of the HF energies in each event. Very central events (centrality approaching 0%) are characterized by a large energy deposit in the HF calorimeters. The results reported in this paper are presented up to 60% in centrality. The minimum bias trigger and event selections are fully efficient in this centrality range.

Track reconstruction [37, 41] is performed in two iterations to ease the computational load for high-multiplicity central $\mathrm {PbPb}$ collisions. The first iteration reconstructs tracks from signals (“hits”) in the silicon pixel and strip detectors compatible with a trajectory of $p_{\mathrm {T}} >0.9{\,\text {GeV/}c}$. The significance of the separation along the beam axis (z) between the track and the primary vertex, $d_z/\sigma (d_z)$, and the significance of the impact parameter relative to the primary vertex transverse to the beam, $d_{\mathrm {0}}/\sigma (d_{\mathrm {0}})$, must be less than 2. In addition, the relative uncertainty of the $p_{\mathrm {T}}$ measurement, $\sigma (p_{\mathrm {T}})/$ $p_{\mathrm {T}}$, must be less than 5%, and tracks are required to have at least 11 out of the possible 14 hits along their trajectories in the pixel and strip trackers. To reduce the number of misidentified tracks, the chi-squared per degree of freedom, $\chi ^2/\mathrm {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 having hits along the trajectory of the track. The second iteration reconstructs tracks compatible with a trajectory of $p_{\mathrm {T}} >0.2{\,\text {GeV/}c}$ using solely the pixel detector. These tracks are required to have $d_{z}/\sigma (d_{z}) < 6$ and a fit $\chi ^2/\mathrm {dof}$ value less than 9 times the number of layers with hits along the trajectory of the track. In the final analysis, first iteration tracks with $p_{\mathrm {T}} > 1.0{\,\text {GeV/}c}$ are combined with pixel-detector-only tracks that have $0.2< p_{\mathrm {T}} < 2.4{\,\text {GeV/}c}$. After removing duplicates [7], the merged track collection has a combined geometric acceptance and efficiency exceeding 60% for $p_{\mathrm {T}}$ $\approx 1.0{\,\text {GeV/}c}$ and $|\eta |<0.8$, as determined using the hydjet event generator [42]. When the track $p_{\mathrm {T}}$ is below 1${\,\text {GeV/}c}$, the acceptance and efficiency steadily drops, reaching approximately 40% at $p_{\mathrm {T}} \approx 0.3{\,\text {GeV/}c}$, which is the lower limit for $p_{\mathrm {T}}$ in this analysis.

4 Analysis technique

The analysis technique follows the method described in Refs. [1, 28] using detector information from both HF and the tracker. The notation $V_n=v_n\exp (in\varPsi _n)=\left\langle e^{in\phi } \right\rangle$ in Eq. (1) will be replaced by the measured complex flow vector $Q_n$ with real and imaginary parts defined as

$$\begin{aligned} {\text {Re}}(Q_n)= & {} \frac{1}{\sum {w_j}}\sum \limits _j^M {w_j}\cos \left( {n{\phi _j}} \right) \nonumber \\&- \left\langle \frac{1}{\sum {w_j}}\sum \limits _j^M {w_j}\cos \left( {n{\phi _j}} \right) \right\rangle , \end{aligned}$$

(3)

$$\begin{aligned} {\text {Im}}(Q_n)= & {} \frac{1}{\sum {w_j}}\sum \limits _j^M {w_j}\sin \left( {n{\phi _j}} \right) \nonumber \\&- \left\langle \frac{1}{\sum {w_j}}\sum \limits _j^M {w_j}\sin \left( {n{\phi _j}} \right) \right\rangle , \end{aligned}$$

(4)

where M represents the number of tracks or HF towers used for calculating the Q vector, $\phi _j$ is the azimuthal angle of the jth track or HF tower, and $w_j$ is a weighting factor equal to transverse energy for HF Q vectors. To correct for the tracking inefficiency, $w_j = 1/\varepsilon _j$ is the inverse of the tracking efficiency $\varepsilon _j(p_{\mathrm {T}}, \eta )$ of the jth track. Unlike the averages over particles in a single event in the definitions of $Q_n$, the angle brackets in Eqs. (3) and (4) denote an average over all the events within a given centrality range. Subtraction of the event-averaged quantity removes biases due to the detector acceptance.

The mixed higher-order harmonics in each $p_{\mathrm {T}}$ range are extracted using the scalar-product method as shown in Eqs. (5)–(9) [1], which describe the various harmonics measured with respect to symmetry plane angles of different orders. Equations (5)–(9) show $v_4$ with respect to the second-order, $v_5$ with respect to the second- and third-order, $v_6$ with respect to the second-order, $v_6$ with respect to the third-order, and $v_7$ with respect to the second- and third-order symmetry plane angles, respectively.

$$\begin{aligned} v_{4}\{\varPsi _{22}\}\equiv & {} \frac{{\mathrm {Re}}\langle Q_{4} Q_{2A}^{*}Q_{2B}^{*}\rangle }{\sqrt{{\mathrm {Re}}\langle Q_{2A}Q_{2A}Q_{2B}^{*}Q_{2B}^{*} \rangle }} \end{aligned}$$

(5)

$$\begin{aligned} v_{5}\{\varPsi _{23}\}\equiv & {} \frac{{\mathrm {Re}}\langle Q_{5} Q_{2A}^{*}Q_{3B}^{*}\rangle }{\sqrt{{\mathrm {Re}}\langle Q_{2A}Q_{3A}Q_{2B}^{*}Q_{3B}^{*} \rangle }} \end{aligned}$$

(6)

$$\begin{aligned} v_{6}\{\varPsi _{222}\}\equiv & {} \frac{{\mathrm {Re}}\langle Q_{6} Q_{2A}^{*}Q_{2B}^{*}Q_{2B}^{*}\rangle }{\sqrt{{\mathrm {Re}}\langle Q_{2A}Q_{2A}Q_{2A}Q_{2B}^{*}Q_{2B}^{*}Q_{2B}^{*} \rangle }} \end{aligned}$$

(7)

$$\begin{aligned} v_{6}\{\varPsi _{33}\}\equiv & {} \frac{{\mathrm {Re}}\langle Q_{6} Q_{3A}^{*}Q_{3B}^{*}\rangle }{\sqrt{{\mathrm {Re}}\langle Q_{3A}Q_{3A}Q_{3B}^{*}Q_{3B}^{*} \rangle }} \end{aligned}$$

(8)

$$\begin{aligned} v_{7}\{\varPsi _{223}\}\equiv & {} \frac{{\mathrm {Re}}\langle Q_{7} Q_{2A}^{*}Q_{2B}^{*}Q_{3B}^{*}\rangle }{\sqrt{{\mathrm {Re}}\langle Q_{2A}Q_{2A}Q_{3A}Q_{2B}^{*}Q_{2B}^{*}Q_{3B}^{*} \rangle }} \end{aligned}$$

(9)

Here, $Q_{nA}$ and $Q_{nB}$ are vectors from two different parts of the detector, specifically the positive and negative sides of HF, $Q_{n}$ is the vector from charged particles in each $p_{\mathrm {T}}$ range within $|\eta |<0.8$, and angle brackets denote the average (weighted by the number of particles) over all events within a given centrality range. The minimum $\eta$ gap between tracks used to find the charged-particle Q vector and towers used for the HF Q vectors is 2.2 units of $\eta$.

With the assumption that the linear and nonlinear terms in Eq. (2) are uncorrelated, the nonlinear response coefficients in each $p_{\mathrm {T}}$ range can be expressed as [1, 28],

$$\begin{aligned} \chi _{422}= & {} \frac{{\mathrm {Re}}\langle Q_{4} Q_{2A}^{*}Q_{2B}^{*}\rangle }{{{\mathrm {Re}}\langle Q_{2}Q_{2}Q_{2A}^{*}Q_{2B}^{*} \rangle }}, \end{aligned}$$

(10)

$$\begin{aligned} \chi _{523}= & {} \frac{{\mathrm {Re}}\langle Q_{5} Q_{2A}^{*}Q_{3B}^{*}\rangle }{{{\mathrm {Re}}\langle Q_{2}Q_{3}Q_{2A}^{*}Q_{3B}^{*} \rangle }}, \end{aligned}$$

(11)

$$\begin{aligned} \chi _{6222}= & {} \frac{{\mathrm {Re}}\langle Q_{6} Q_{2A}^{*}Q_{2B}^{*}Q_{2B}^{*}\rangle }{{{\mathrm {Re}}\langle Q_{2}Q_{2}Q_{2}Q_{2A}^{*}Q_{2B}^{*}Q_{2B}^{*} \rangle }}, \end{aligned}$$

(12)

$$\begin{aligned} \chi _{633}= & {} \frac{{\mathrm {Re}}\langle Q_{6} Q_{3A}^{*}Q_{3B}^{*}\rangle }{{{\mathrm {Re}}\langle Q_{3}Q_{3}Q_{3A}^{*}Q_{3B}^{*} \rangle }}, \end{aligned}$$

(13)

$$\begin{aligned} \chi _{7223}= & {} \frac{{\mathrm {Re}}\langle Q_{7} Q_{2A}^{*}Q_{2B}^{*}Q_{3B}^{*}\rangle }{{{\mathrm {Re}}\langle Q_{2}Q_{2}Q_{3}Q_{2A}^{*}Q_{2B}^{*}Q_{3B}^{*} \rangle }}, \end{aligned}$$

(14)

where the charged-particle $Q_{n}$ vector enters both the numerator and the denominator.

5 Systematic uncertainties

Six sources of systematic uncertainties are considered in this analysis. The systematic uncertainty due to vertex position selection is estimated by comparing the results with events from vertex position ranges $|v_z |<3$ $\, \text {cm}$ to $3<|v_z |<15$ $\, \text {cm}$. For both mixed harmonic and nonlinear response coefficients, this uncertainty is estimated to be 1–3%, with no dependence on $p_{\mathrm {T}}$ or centrality. Systematic uncertainty due to track quality requirements are examined by varying the track selections for $d_z/\sigma (d_z)$ and $d_{\mathrm {0}}/\sigma (d_{\mathrm {0}})$ from 1.5 to 5, the pixel track $d_{z}/\sigma (d_{z})$ from 5 to 10, and the fit $\chi ^2/\mathrm {dof}$ value from 7 to 18 times the number of layers with hits. The uncertainty is estimated to be 1–4% depending on $p_{\mathrm {T}}$ and centrality for both mixed harmonic and nonlinear response coefficients.

The charged-particle tracking efficiency depends on the efficiency of detecting different types of charged particles and the species composition of the set of particles. Two event generators (hydjet [42] and epos lhc [43]) with different particle composition are used to study the tracking efficiency, and the systematic uncertainty is obtained by comparing the results using efficiencies from the two generators mentioned above. The systematic uncertainty from this source is 3% for the mixed harmonics and less than 1% for the nonlinear response coefficients, with no dependence on $p_{\mathrm {T}}$ or centrality.

The sensitivity of the results to the centrality calibration is evaluated by varying the trigger and event selection efficiency by $\pm 2$%. The resulting uncertainty is estimated to be less than 1%. The minimum $\eta$ gap between the correlated charged particles and the Q vectors in the HF region is changed from 2.2 to 3.2 units of $\eta$ (achieved by changing the $\eta$ ranges of the HF Q vectors) to estimate the uncertainty due to short-range correlations from resonance decays and jets. This study results in a systematic uncertainty of 1–8%, depending on both $p_{\mathrm {T}}$ and centrality. This $\eta$ gap uncertainty also includes a possible physics effect from the $\eta$-dependent fluctuations of symmetry plane angles [26, 44], although a recent study from the ALICE experiment indicates that this effect is small for correlations between symmetry plane angles of different order [45].

When the same set of HF towers are used for different Q vectors in the equations of mixed harmonic and nonlinear response coefficients, the product of these Q vectors contains self-correlations. An algorithm for removing the duplicated terms when multiplying two or more Q vectors, the same as the approach of Ref. [46], is used. The algorithm only works perfectly when the detector has fine granularity and there is no merging of HF towers. Therefore, the difference before and after correcting for this effect is taken as the systematic uncertainty, yielding values which depend on centrality but are always less than 3%.

The different systematic sources described above are added in quadrature to obtain the overall systematic uncertainty, which is about 10% at low $p_{\mathrm {T}}$ and decreases to around 5% for $p_{\mathrm {T}}$ larger than 1${\,\text {GeV/}c}$. As a function of centrality, the overall systematic uncertainty ranges from 3 to 9% for different coefficients, with larger uncertainties for central events.

6 Results

The measurements in this paper are presented using tracks in the range of $|\eta | < 0.8$. Figure 1 shows the mixed higher-order flow harmonics, $v_4\{\varPsi _{22}\}$, $v_5\{\varPsi _{23}\}$, $v_6\{\varPsi _{222}\}$, $v_6\{\varPsi _{33}\}$, and $v_7\{\varPsi _{223}\}$ from the scalar-product method at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76$ and 5.02$\,\text {TeV}$ as a function of $p_{\mathrm {T}}$ in the 0–20% (upper row) and 20–60% (lower row) centrality ranges.

Fig. 1

Mixed higher-order flow harmonics, $v_4\{\varPsi _{22}\}$, $v_5\{\varPsi _{23}\}$, $v_6\{\varPsi _{222}\}$, $v_6\{\varPsi _{33}\}$, and $v_7\{\varPsi _{223}\}$ from the scalar-product method at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76$ and 5.02$\,\text {TeV}$ as a function of $p_{\mathrm {T}}$ in the 0–20% (upper row) and 20–60% (lower row) centrality ranges. Statistical (bars) and systematic (shaded boxes) uncertainties are shown

Fig. 2

Comparison of mixed higher-order flow harmonics, $v_4\{\varPsi _{22}\}$, $v_5\{\varPsi _{23}\}$, $v_6\{\varPsi _{222}\}$, $v_6\{\varPsi _{33}\}$ and $v_7\{\varPsi _{223}\}$ with the corresponding overall flow, $v_4\{\varPsi _{4}\}$, $v_5\{\varPsi _{5}\}$, $v_6\{\varPsi _{6}\}$, $v_6\{\varPsi _{6}\}$ and $v_7\{\varPsi _{7}\}$, respectively, at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 5.02\,\text {TeV}$ as a function $p_{\mathrm {T}}$ in the 0–20% (upper row) and 20–60% (lower row) centrality ranges. Statistical (bars) and systematic (shaded boxes) uncertainties are shown

It is observed that the shapes of the mixed higher-order flow harmonics as a function of $p_{\mathrm {T}}$ are qualitatively similar to the published overall flow harmonics with respect to $\varPsi _n$ [7, 11], first increasing at low $p_{\mathrm {T}}$, reaching a maximum at about 3–4${\,\text {GeV/}c}$, then decreasing at higher $p_{\mathrm {T}}$. This may indicate that, for each $p_{\mathrm {T}}$ region, the underlying physics processes that generate the flow harmonics are the same for the nonlinear and the linear parts. Similar to previous observation that the overall flow shows a weak energy dependence from RHIC to LHC energies [7, 8], the mixed harmonics are also found to be consistent between the two collision energies within the uncertainties, except for $v_4\{\varPsi _{22}\}$ and $v_5\{\varPsi _{23}\}$ at $p_{\mathrm {T}}$ larger than 3${\,\text {GeV/}c}$ in the mid-central collisions, with 5.02$\,\text {TeV}$ results slightly above 2.76$\,\text {TeV}$ results.

Fig. 3

Nonlinear response coefficients, $\chi _{422}$, $\chi _{523}$, $\chi _{6222}$, $\chi _{633}$, and $\chi _{7223}$ from the scalar-product method at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76$ and 5.02$\,\text {TeV}$ as a function of $p_{\mathrm {T}}$ in the 0–20% (upper row) and 20–60% (lower row) centrality ranges. Statistical (bars) and systematic (shaded boxes) uncertainties are shown. The results are compared with hydrodynamic predictions [30] at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76\,\text {TeV}$ with $\eta /s = 0.08$ and Glauber initial conditions in the 5–10% (blue lines) and 35–40% (dashed green lines) centrality ranges

A direct comparison of the mixed higher-order flow harmonics and overall flow at 5.02$\,\text {TeV}$ is presented in Fig. 2 as a function of $p_{\mathrm {T}}$ in the two centrality ranges. Hydrodynamic models predict that the contribution of the nonlinear response to the overall flow increases towards peripheral collisions for $v_4$ and $v_5$ [17, 29, 47]. From a comparison of the relative contribution in the two centrality ranges, the present results are consistent with these predictions, as well as an estimate by the ATLAS Collaboration using a two-component fit of the correlation between flow harmonics [21], and a recent study of the nonlinear mode by the ALICE Collaboration [45]. By comparing different harmonics, the contribution of the nonlinear response for $v_5$ is larger than those for the other harmonics in the centrality range 20–60%.

The nonlinear response coefficients, $\chi _{422}$, $\chi _{523}$, $\chi _{6222}$, $\chi _{633}$, and $\chi _{7223}$ are presented in Fig. 3 as a function of $p_{\mathrm {T}}$ in the two centrality ranges. It is observed that the odd harmonic coefficients $\chi _{523}$ and $\chi _{7223}$ are larger than those for the even harmonics for $p_{\mathrm {T}}$ less than 3${\,\text {GeV/}c}$ in the two explored centrality ranges. The values for the even harmonics first decrease slightly as $p_{\mathrm {T}}$ increases, reach a minimum at $p_{\mathrm {T}}$ about 2${\,\text {GeV/}c}$, and then slowly increase until appearing to plateau for $p_{\mathrm {T}}$ above 4${\,\text {GeV/}c}$. The results are compared with viscous hydrodynamic predictions [30] at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76\,\text {TeV}$ with $\eta /s = 0.08$ (where $\eta /s$ is the shear viscosity to entropy density ratio of the hydrodynamic medium, and here $\eta$ denotes shear viscosity rather than pseudorapidity) and Glauber initial conditions in two centrality ranges (5–10% and 35–40%) which roughly match those of the data (0–20% and 20–60%). In the model, as $p_{\mathrm {T}}$ increases from 0.3 to 1${\,\text {GeV/}c}$, the predicted coefficients increase for $n = 4$ and 5, but decrease and then increase for $n = 6$ and 7, with a much stronger $p_{\mathrm {T}}$ dependence than the data. The strong $p_{\mathrm {T}}$ dependence, attributed to the large variance of the flow angles $\varPsi _n$ at small $p_{\mathrm {T}}$ [30], is not observed in data for $n = 4$ and 5.

Fig. 4

Mixed higher-order flow harmonics, $v_4\{\varPsi _{22}\}$, $v_5\{\varPsi _{23}\}$, $v_6\{\varPsi _{222}\}$, $v_6\{\varPsi _{33}\}$, and $v_7\{\varPsi _{223}\}$ from the scalar-product method at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76$ and 5.02$\,\text {TeV}$, as a function of centrality. Statistical (bars) and systematic (shaded boxes) uncertainties are shown. Hydrodynamic predictions [1] with $\eta /s = 0.08$ (blue lines) at 2.76$\,\text {TeV}$ are shown in (b) and (e)

Figure 4 shows the mixed higher-order flow harmonics, $v_4\{\varPsi _{22}\}$, $v_5\{\varPsi _{23}\}$, $v_6\{\varPsi _{222}\}$, $v_6\{\varPsi _{33}\}$, and $v_7\{\varPsi _{223}\}$ from the scalar-product method, as a function of centrality in the $p_{\mathrm {T}}$ range from 0.3 to 3.0${\,\text {GeV/}c}$. Hydrodynamic predictions with a deformed symmetric Gaussian density profile as the initial conditions for $v_5\{\varPsi _{23}\}$ and $v_7\{\varPsi _{223}\}$ [1] at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76\,\text {TeV}$ are compared with the data. The model qualitatively describes $v_5\{\varPsi _{23}\}$ in the 0–40% centrality range but underestimates the result for more peripheral collisions. For $v_7\{\varPsi _{223}\}$, the predicted values are much smaller than the data, especially for centrality from 35 to 50%.

Fig. 5

Nonlinear response coefficients, $\chi _{422}$, $\chi _{523}$, $\chi _{6222}$, $\chi _{633}$, and $\chi _{7223}$ from the scalar-product method at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76$ and 5.02$\,\text {TeV}$, as a function of centrality. Statistical (bars) and systematic (shaded boxes) uncertainties are shown. The results are compared with predictions at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76\,\text {TeV}$ from AMPT [48] as well as hydrodynamics with a deformed symmetric Gaussian density profile as the initial conditions using $\eta /s = 0.08$ from Ref. [1], and from iEBE-VISHNU hydrodynamics with both Glauber and the KLN initial conditions using the same $\eta /s$ [28]

The nonlinear response coefficients, $\chi _{422}$, $\chi _{523}$, $\chi _{6222}$, $\chi _{633}$, and $\chi _{7223}$ are presented in Figs. 5 and 6, as a function of centrality in the $p_{\mathrm {T}}$ range from 0.3 to 3.0${\,\text {GeV/}c}$. The results are compared with predictions at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76\,\text {TeV}$ from the microscopic transport model AMPT [48, 49], a macroscopic hydrodynamic model using a deformed symmetric Gaussian density profile as the initial conditions with $\eta /s = 0.08$ [1], and from another hydrodynamic calculation (iEBE-VISHNU) with both Glauber and Kharzeev–Levin–Nardi (KLN) gluon saturation initial conditions using the same $\eta /s$ [28]. The model with Gaussian profile initial conditions gives a better description of the nonlinear response coefficients compared to other calculations, but it underestimates the values of $v_7\{\varPsi _{223}\}$ for centrality above 30%, as shown in Fig. 4. In Fig. 6, the same results are compared with the predictions from hydrodynamics $+$ hadronic cascade hybrid approach with the IP-Glasma initial conditions using $\eta /s = 0.095$ [50] at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 5.02\,\text {TeV}$ and from iEBE-VISHNU hydrodynamics with the KLN initial conditions using $\eta /s = 0$, 0.08 and 0.2 [28] at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76\,\text {TeV}$. All the calculations describe the $\chi _{422}$ well, but none of them are successful for $\chi _{523}$ and $\chi _{7223}$. The model calculations of $\chi _{7223}$ are quite different for various initial conditions and $\eta /s$, which suggests that the first-time measurement of $\chi _{7223}$ presented in this paper could provide strong constraints on models.

Fig. 6

The same results as in Fig. 5 but compared with predictions from a hydrodynamics $+$ hadronic cascade hybrid approach with the IP-Glasma initial conditions using $\eta /s = 0.095$ [50] at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 5.02\,\text {TeV}$ and from iEBE-VISHNU hydrodynamics with the KLN initial conditions using $\eta /s = 0$, 0.08 (the same curve as in Fig. 5) and 0.2 [28] at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76\,\text {TeV}$

7 Summary

The mixed higher-order flow harmonics and nonlinear response coefficients of charged particles have been studied as functions of transverse momentum $p_{\mathrm {T}}$ and centrality in $\mathrm {PbPb}$ collisions at $\sqrt{\smash [b]{s_{_{\mathrm {NN}}}}} = 2.76$ and 5.02$\,\text {TeV}$ using the CMS detector. The measurements use the scalar-product method, covering a $p_{\mathrm {T}}$ range from 0.3 to 8.0${\,\text {GeV/}c}$, pseudorapidity $|\eta |<0.8$, and a centrality range of 0–60%. The mixed higher-order flow harmonics, $v_4\{\varPsi _{22}\}$, $v_5\{\varPsi _{23}\}$, $v_6\{\varPsi _{222}\}$, $v_6\{\varPsi _{33}\}$, and $v_7\{\varPsi _{223}\}$ all have a qualitatively similar $p_{\mathrm {T}}$ dependence, first increasing at low $p_{\mathrm {T}}$, reaching a maximum at about 3–4${\,\text {GeV/}c}$, and then decreasing at higher $p_{\mathrm {T}}$. As a comparison, the overall $v_n$ harmonics ($n = 4$–7) with respect to their own symmetry planes are measured in the same $p_{\mathrm {T}}$, $\eta$, and centrality ranges. The relative contribution of the nonlinear part for $v_5$ is larger than for other harmonics in the centrality range 20–60%. In addition, the nonlinear response coefficients of the odd harmonics are observed to be larger than those of even harmonics for $p_{\mathrm {T}}$ less than 3${\,\text {GeV/}c}$. At $p_{\mathrm {T}}$ less than 1${\,\text {GeV/}c}$, a viscous hydrodynamic calculation with Glauber initial conditions and shear viscosity to entropy density ratio $\eta /s = 0.08$ predicts a much stronger $p_{\mathrm {T}}$ dependence for the nonlinear response coefficients. The coefficients, including the first-time measurement of $\chi _{7223}$, as a function of centrality, are compared with AMPT and hydrodynamic predictions using different $\eta /s$ and initial conditions. Compared to the data, none of the models provides a simultaneous description of the mixed higher-order flow harmonics and nonlinear response coefficients. Therefore, these results can constrain both initial conditions and transport properties of the produced medium.