Probing charm quark dynamics via multiparticle correlations in PbPb collisions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV

Multiparticle azimuthal correlations of prompt D$^0$ mesons are measured in PbPb collisions at a nucleon-nucleon center-of-mass energy of $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV. For the first time, a four-particle cumulant method is used to extract the second Fourier coefficient of the azimuthal distribution ($v_2$) of D$^0$ mesons as a function of event centrality and the D$^0$ transverse momentum. The ratios of the four-particle $v_2$ values to previously measured two-particle cumulant results provide direct experimental access to event-by-event fluctuations of charm quark azimuthal anisotropies. These ratios are also found to be comparable to those of inclusive charged particles in the event. However, hints of deviations are seen in the most central and peripheral collisions. To investigate the origin of flow fluctuations in the charm sector, these measurements are compared with models implementing fluctuations of charm quark energy loss via collisional or radiative processes in the quark-gluon plasma. These models cannot quantitatively describe the data over the full transverse momentum and centrality ranges, although the calculations with collisional energy loss provide a better description of the data.


1
A strongly coupled quark-gluon plasma (QGP), has been studied in nucleus-nucleus collisions at the BNL RHIC [1][2][3][4] and CERN LHC [5,6]. This medium exhibits the behavior of a nearlyperfect liquid [7,8]. The azimuthal anisotropy of produced hadrons, resulting from pressuredriven expansion, is a powerful tool to study QGP dynamics and can be characterized by the Fourier coefficients (v n ) of the hadrons' azimuthal angle (φ) distribution [9]. The second-order Fourier coefficient (v 2 ), known as elliptic flow, of low transverse momentum (p T ) particles reflects the QGP response to the average initial collision geometry and its event-by-event fluctuations [9]. The v 2 coefficient is also influenced by the path length dependence of parton energy loss at high p T [10][11][12].
The magnitude of event-by-event fluctuations [29] of azimuthal anisotropy harmonics from heavy-flavor quarks has not been experimentally measured. Multiparticle correlation techniques involving four or more particles, v 2 {n}, with n ≥ 4 [30], allow direct access to cumulants of the v 2 probability density distribution. The technique has been widely applied in the light-flavor sector to extract the magnitude of v 2 fluctuations, which is then used to constrain fluctuations of the initial-state geometry. It has been recently proposed that for hard probes (such as high-p T jets, and heavy-flavor hadrons), fluctuations of anisotropy harmonics are not only influenced by the initial-state geometry, but are also sensitive to final-state fluctuations of energy loss when these hard probes propagate in the QGP medium [31]. Therefore, measurements of v 2 {4} and its ratio to v 2 {2} for heavy-flavor hadrons have the potential to set constraints on the mechanism of heavy-quark energy loss, especially how it fluctuates on an event-by-event basis in QGP.
In this letter, the prompt D 0 meson v 2 coefficient is measured for the first time using fourparticle correlations, and the ratio v 2 {4}/v 2 {2} is presented. These measurements use data from lead-lead (PbPb) collisions at a nucleon-nucleon center-of-mass energy of √ s NN = 5.02 TeV, collected by the CMS detector at the LHC in 2018. The behavior of v 2 is examined in the rapidity (y) range |y| < 1 over the p T range of 2-15 GeV, and in the event centrality classes (i.e. the percentage ranges of the total inelastic hadronic cross section) of 10-30% and 30-50%. A 0% centrality corresponds to the largest overlap of the two nuclei. The centrality dependence of v 2 is also measured over the broader range of 5-60% for 2 < p T < 8 GeV. Tabulated results are provided in the HEPData record for this analysis [32].
The CMS apparatus [33] is a multipurpose, nearly hermetic detector, designed to trigger on [34, 35] and identify electrons, muons, photons, and hadrons [36][37][38][39]. In this analysis, the information from two subdetectors were used: the silicon inner tracker, which measures charged particles within the range of pseudorapidity |η| < 3.0; and the hadronic forward (HF) calorimeters, made of steel and quartz fibres, which extend the pseudorapidity coverage provided by the barrel and endcap detectors to about |η| < 5.0, and are segmented to form 0.175×0.175 (∆η×∆φ) towers.
The data analyzed consist of 4.27×10 9 minimum bias events, corresponding to an integrated luminosity of 0.58 nb −1 . The events are triggered by requiring signals above thresholds in the range of ∼6-12 GeV in both sides of the HF calorimeters [35]. Events must also have at least one reconstructed primary vertex within 15 cm of the interaction point along the beam axis. The primary vertex is selected as the one with the highest track multiplicity in the event. The effects from concurrent interactions in the same bunch crossing were shown to be negligible. The centrality is calculated using the HF calorimeters [40].
Monte Carlo (MC) event samples are simulated containing either prompt or nonprompt D 0 mesons; the latter originate from beauty hadron decays. The simulated events are generated using PYTHIA 8.212 [41], tune CP5 [42], and embedded into MC PbPb events from HYDJET 1.9 [43]. The prompt D 0 meson event sample is employed to define signal selections and efficiency corrections, while the other sample is used to estimate systematic uncertainties from nonprompt D 0 contamination.
Both D 0 and D 0 mesons are reconstructed via the process D 0 → π + + K − (D 0 → π − + K + ), with a branching fraction of (3.95  This analysis shares the same data sets and uses a similar procedure to that described in Ref. [18], in which the D 0 meson v 2 is measured using the two-particle correlation (or cumulant) method, v 2 {2}(D 0 ), where the D 0 meson v 2 signal is extracted by correlating a D 0 meson with reference particles measured in the HF detectors. To measure the differential second-order (elliptic) harmonic from the four-particle cumulant, v 2 {4}(D 0 ) [30], a first step involves either two-or four-particle correlations calculated using energy deposits in the HF towers to obtain elliptic harmonics of reference particles. Here, each HF tower is used to represent one or more particles with a weight applied corresponding to its deposited transverse energy in the calculation of cumulants when averaging over all HF towers, as detailed below. The two-and four-particle azimuthal correlations for the nth harmonic are defined as Here, φ j (j = 1, . . . , 4) are the azimuthal angles of one unique combination of multiple particles in an event and the double average symbol · · · indicates that the average is taken over all unique particle combinations and for all events. In addition, the superscripts a and b indicate towers chosen from two different HF calorimeters, HF− (−5 < η < −3) or HF+ (3 < η < 5).
In a second step, the four-particle cumulant of reference particle azimuthal correlations, To measure the prompt D 0 meson v 2 coefficient, the φ a 1 from an HF tower in Eq.(1) is replaced with a D 0 candidate's azimuthal angle selected within the tracker acceptance |η| < 2.4. To suppress the nonflow effects from sources such as resonance decays or jets, in the two-particle cumulant method, a tower with φ 2 is selected from the HF calorimeter (HF+ or HF−) having the opposite η sign as that of the D 0 candidate. For the four-particle correlations method, φ a is picked from the HF detector having the same η sign as the D 0 candidate, but φ b 3 and φ b 4 are chosen from the other HF detector. Studies performed with simulated events indicate that nonflow effects are negligible when measuring v 2 {2} with these η gaps [48, 50]. These modified particle correlators involving a D 0 meson are denoted by 2 and 4 . The differential four-particle cumulant of D 0 mesons is then defined as [30, 49], Finally, with respect to the reference four-particle cumulants, the differential four-particle v n (D 0 ) coefficients are extracted as in Refs. [30,49], [18] is performed. First, the invariant mass spectrum of all D 0 candidates is fit using a formula containing five components: (i) A sum of two Gaussian functions having the same mean but different widths are used for the D 0 signal; (ii) a single Gaussian function describes the invariant mass spectrum of D 0 candidates with an incorrect mass assignment resulting from the exchange of the kaon and pion designations; (iii) a Crystal Ball function [51] is used for the processes D 0 → K + K − [52]; (iv) another Crystal Ball function to describe D 0 → π + π − [52]; (v) a third-order polynomial is used to model the combinatorial background. The first four components are initialized by values calculated using simulated events, and their widths are allowed to vary with a common scale factor during the fit to data. Using the signal and background D 0 candidate yield fraction extracted from the invariant mass fit, the measured v 2 data of all D 0 candidates, v   [11,12]) with the option of turning on energy loss by gluon radiation or alternatively by elastic collisions described by Langevin dynamics during the heavy-quark propagation, are also shown in the upper panel of Fig. 2. The radiative energy loss process is expected to be the dominant phenomenon in the high-p T region. Langevin dynamics, which describe the propagation of heavy quarks in the medium as a Brownian motion, can account for collisional processes using Langevin-like equations [11] in the low-and intermediate-p T region. Both models seem to capture the general trends of the data, without reproducing them quantitatively.
To further investigate the underlying physics processes behind elliptic flow fluctuations of charm quarks, the ratios v 2 {4}/v 2 {2} are presented as a function of p T , up to 15 GeV, in the lower panel of Fig. 2. Generally speaking, a larger deviation of v 2 {4}/v 2 {2} ratios from unity indicates a larger magnitude of flow fluctuations. The same ratios for charged particles (dominated by light-flavor hadrons) are shown. The ratios for prompt D 0 mesons are consistent with those for charged particles. The roughly flat behavior of the ratios at low p T suggests that initial-state geometry fluctuations are likely the dominant source of flow fluctuations there [11]. The ratios based on the DABMod model for D 0 mesons [11,12], also shown in Fig. 2 (bottom), lie systematically above the data, suggesting an underestimation of the magnitude of flow fluctuations in the data.
The p T -integrated results of v 2 {4} for 2 < p T < 8 GeV and |y| < 1 are shown as a function of centrality from 5 to 60% in Fig. 3  are plotted for comparison. The v 2 {4}/v 2 {2} ratios are shown in the lower panel of Fig. 3. The prompt D 0 data are also compared to those of inclusive charged particles within the range |η| < 1 and 2 < p T < 8 GeV.
Similar to the D 0 meson v 2 {2} coefficient, the D 0 meson v 2 {4} value increases with centrality in the 5-40% range, and then decreases for more peripheral collisions. This trend is qualitatively reproduced by calculations incorporating an interplay of initial-state geometry and parton energy loss in QGP. Within the 10-40% centrality range, the v 2 {4}/v 2 {2} ratios are almost identical between prompt D 0 mesons and inclusive charged particles within uncertainties. This indicates that, within this centrality range, the dominant source of flow fluctuations for heavyflavor is similar to that for soft light-flavor particles, namely initial-state geometry fluctuations, and therefore the contribution from final-state fluctuations is small. The hint of different trends in v 2 {4}/v 2 {2} between D 0 mesons and charged particles seen in the most central and most peripheral events could indicate that fluctuations from final-state effects, such as parton energy loss, in hard processes become visible for charm mesons [11]. For example, as the system size becomes smaller for peripheral events, the number of scatterings a hard probe experiences with QGP will decrease, leading to larger fluctuations in the energy loss on an event-by-event basis. However, the experimental uncertainties are still large, with the difference of ∼2 standard deviations between the values. Calculations based on the DABMod model [11,12] assuming collisional (or Langevin dynamics) and radiative energy loss processes are also shown in Fig. 3. A better description of the experimental data is obtained using the Langevin dynamics, although no increase or decrease for the most central or peripheral events, respectively, is predicted.