Measurements of azimuthal anisotropy of nonprompt D$^0$ mesons in PbPb collisions at $\sqrt{s_\mathrm{NN}}$ = 5.02 TeV

Measurements of the elliptic ($v_2$) and triangular ($v_3$) azimuthal anisotropy coefficients are presented for D$^0$ mesons produced in b hadron decays (nonprompt D$^0$ mesons) in lead-lead collisions at $\sqrt{s_\mathrm{NN}}$ = 5.02 TeV. The results are compared with previously published charm meson anisotropies measured using prompt D$^0$ mesons. The data were collected with the CMS detector in 2018 with an integrated luminosity of 0.58 nb$^{-1}$. Azimuthal anisotropy is sensitive to the interactions of quarks with the hot and dense medium created in heavy ion collisions. Comparing results for prompt and nonprompt D$^0$ mesons can assist in understanding the mass dependence of these interactions. The nonprompt results show lower magnitudes of $v_2$ and $v_3$ and weaker dependences on the meson transverse momentum and collision centrality than those found for prompt D$^0$ mesons. The results are in agreement with theoretical predictions that include a mass dependence in the interactions of quarks with the medium.


Introduction
In ultra-relativistic heavy ion collisions, cold nuclear matter transforms into a state of strongly coupled matter where quarks and gluons are deconfined from hadrons, called the quark-gluon plasma (QGP) [1][2][3][4].One of the main features of the QGP is the collective motion of its constituents.A signature of this collective motion is the presence of a momentum anisotropy of the particles that emerge from the collision.In the plane perpendicular to the beam axis (transverse plane), this anisotropy results from the conversion to momentum space of the initial transverse spatial anisotropy of the colliding system, a phenomenon known as collective anisotropic flow.To quantify this effect, azimuthal (Fourier) coefficients of order n, v n = ⟨cos[n(ϕ − Ψ n )]⟩ can be used, where ϕ is the single-particle azimuthal angle and Ψ n is the first angle where the nth harmonic component has its maximum outgoing particle density [5,6].The measured v n values for low-momentum light hadrons can be described by relativistic hydrodynamics [7][8][9].The second Fourier coefficient, v 2 (referred to as elliptic flow), arises from the combined effects of the elongated shape and the event-by-event fluctuations of the overlap region of the colliding nuclei [10].The third coefficient, v 3 (referred to as triangular flow), is predominantly due to the fluctuations.
Because of their large mass, bottom (b) and charm (c) quarks are produced at the earliest stage of the collision [11], but their azimuthal distributions can be affected by their interaction with the medium as they travel through the QGP [12].At low transverse momentum (p T ), it is expected that charm quarks thermalize and, therefore, they should follow the motion of the gluons and light flavor quarks [13].In the high p T region, the azimuthal asymmetry could be caused by a path length dependence of the parton energy loss [14,15].Hence, studying the flow coefficients of heavy quarks in heavy ion collisions can provide important inputs for understanding the properties of the QGP.Significant v 2 and v 3 coefficients of charm hadrons in lead-lead (PbPb) collisions have already been measured at the CERN LHC, proving that charm quarks exhibit collective flow [16][17][18][19][20][21][22].An early measurement by the CMS Collaboration of the elliptic flow of J/ψ from b hadron decays was consistent with zero, albeit with large uncertainties [23].Recently, the ALICE Collaboration measured the v 2 of electrons from b hadron decays [24], and the ATLAS Collaboration measured both v 2 and v 3 of muons from b hadron decays [19].These results indicate that b hadron elliptic flow is nonzero, and also show a clear mass ordering, with smaller measured values of Fourier coefficients than those found for particles containing lighter flavor quarks.Compared to these J/ψ and lepton measurements, the study of the b → D 0 channel is a promising opportunity due to the large branching fraction of this decay, approximately 70% [25].This enables lower p T coverage, thereby increasing sensitivity to the low-p T flow of b hadrons.In addition, because of its large mass, the D 0 momentum is more closely correlated with the b quark momentum than is the case for leptons.
In this Letter, the first v 2 and v 3 measurements of D 0 mesons from b hadron decay (nonprompt D 0 mesons) in large collision systems are reported, using PbPb collisions at a center-of-mass energy per nucleon pair of √ s NN = 5.02 TeV.Results for D 0 mesons in the rapidity range |y| < 1 are shown as a function of their p T , spanning 1-30 GeV/c, and in three classes of PbPb centrality: 0-10%, 10-30%, and 30-50%.These centrality percentages represent the fraction of the total inelastic hadronic cross section, with 0% corresponding to full overlap of the two colliding nuclei.The broad range in p T enables a comprehensive study of different flow generation mechanisms.The results are compared with azimuthal anisotropy measurements of D 0 mesons that do not come from b hadron decays (prompt D 0 mesons) [18], as well as with theoretical predictions.Tabulated results are provided in the HEPData record for this analysis [26].

Experimental setup and data sample
The CMS apparatus [27] is a multipurpose, nearly hermetic detector, designed to trigger on [28,29] and identify electrons, muons, photons, and (charged and neutral) hadrons [30][31][32].A global algorithm [33] aims to reconstruct all individual particles in an event, combining information provided by the all-silicon inner tracker and by the crystal electromagnetic and brassscintillator hadron calorimeters, operating inside a 3.8 T superconducting solenoid, with data from the gas-ionization muon detectors embedded in the flux-return yoke outside the solenoid.The forward hadron (HF) calorimeter uses steel as an absorber and quartz fibers as the sensitive material.The two halves of the HF are located 11.2 m from the interaction region, one on each end, and together they provide coverage in the pseudorapidity range 3.0 < |η| < 5.2.The HF calorimeters are subdivided into "towers" with ∆η×∆ϕ = 0.175×0.175,and energy deposited in a tower is treated as a detected hadron in this analysis.They serve as luminosity monitors and the total energy deposited in both HF calorimeters is used for centrality determination [34].Events are filtered using a two-tiered trigger system.The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors [28].The second level, known as the high-level trigger [29], consists of a farm of processors running a version of the full event reconstruction software.
The data analyzed in this Letter consist of minimum bias (MB) PbPb events corresponding to an integrated luminosity of 0.58 nb −1 [35,36].The MB selection requires signals above readout thresholds in the range of ∼6-12 GeV on both sides of the HF calorimeters [29].The events are further filtered to remove background events (beam-gas interactions and nonhadronic collisions) by applying the procedure described in Ref. [37].The final results include only events with at least one vertex associated with two or more tracks (primary vertex) within 15 cm from the nominal interaction point along the beam direction, and at least two towers with energy larger or equal to 4 GeV in each of the HF detectors.Finally, in order to suppress the contamination from events with multiple ion collisions, an event is rejected if the shapes of the clusters in the pixel part of the inner tracker are not compatible with those expected in single PbPb collisions.The total event selection efficiency is (97 ± 2)%.
Simulated events from Monte Carlo (MC) generators are used to study the kinematics of D 0 mesons in PbPb collisions.The D 0 mesons are generated with PYTHIA 8.212 [38], tune CP5 [39], while their decays are modeled with EVTGEN 1.3 [40].Both prompt and nonprompt D 0 samples are produced.The decay products of the D 0 mesons are then embedded into MB events generated using HYDJET 1.9 [41], tuned to reproduce the charged particle multiplicity and spectra.The response of the CMS detector to these combined events is simulated with GEANT4 [42].

Analysis procedure
Inclusive (both prompt and nonprompt) D 0 and D 0 meson candidates are reconstructed via their decay channels: D 0 → π + + K − and D 0 → π − + K + .Both D 0 and D 0 are referred to as D 0 in this Letter.All tracks used in the analysis have p T > 1 GeV/c and |η| < 1.2.To ensure the best track quality, only tracks that satisfy high purity [32] criteria are considered.In addition, the relative uncertainty in the track p T measurement is required to be less than 10%.The number of hits along the track trajectory in the tracker must be greater or equal to 11.The χ 2 of the track fit divided by the total number of degrees of freedom of the fit, and also divided by the total number of tracker layers having a hit along the track path, must be less than 0.18.Candidates are formed by combining pairs of tracks from oppositely-charged particles and requiring an invariant mass (m inv ) within a ±200 MeV/c 2 window of the world-average D 0 meson mass of 1865 MeV/c 2 [25].Since no particle identification is available, both possible particle assignments (pion and kaon) are considered for each pair of selected candidate decay tracks.
Further selection is done by applying a boosted decision tree (BDT) algorithm implemented using the TMVA package [43].For the BDT training, the signal candidates are taken from the simulated events in which reconstructed D 0 mesons are required to match the generated nonprompt D 0 particles.The background training sample consists of two components combined together: combinatorial background and prompt D 0 production.The ratio of these two components in the background training sample is the same as found in the data.The combinatorial background sample is composed of D 0 candidates reconstructed from data whose mass is three to six standard deviations away from the nominal D 0 meson mass (∼1.795-1.830and ∼1.900-1.935GeV/c 2 ).The prompt D 0 meson component of the background sample consists of D 0 candidates that correspond to simulated prompt D 0 mesons in the MC events.In this way, the training is optimized to favor the D 0 mesons produced in b hadron decays.Training is performed separately for each p T range and centrality class.The variables related to D 0 mesons used to discriminate the signal from the background are: χ 2 probability for the D 0 vertex fit, the distance between the D 0 and primary vertices and its significance, as well as the angle between the momentum of the D 0 meson candidate and the line connecting the primary and D 0 vertices (pointing angle).In addition, related to the decay products of the D 0 meson candidate, the variables used are: the significance of the distances of closest approach (DCA) to the primary vertex (both along and perpendicular to the beam direction), and the number of hits assigned to a particular track.These variables are chosen by analyzing their importance in the decision tree and the correlation matrix among all variables.The BDT cut is optimized to provide the maximal nonprompt D 0 meson signal significance for each analysis bin.Each D 0 candidate is weighted by a correction factor, defined as the inverse of the product of the acceptance and selection efficiencies as functions of p T and centrality.This factor is calculated from the MC sample using the ratio of the number of reconstructed D 0 mesons (that satisfy selection criteria) to the total number of generated nonprompt D 0 mesons in a certain kinematic range.Potential discrepancies between data and simulated p T spectra could lead to an inaccuracy in the correction factor.Therefore, before performing this correction, the p T spectrum of D 0 mesons in the simulation is weighted to match that in the data.The weighting factor is calculated in each centrality range from the ratio of the uncorrected p T distribution of nonprompt D 0 mesons, extracted from the invariant mass and DCA template fits to data, and the model distribution, obtained with the same reconstruction and selection criteria applied to data.
The measurement of inclusive D 0 meson v n coefficients employs the scalar-product (SP) method [44] used in previous CMS publications for prompt D 0 mesons [17,18].
In this method, the v n coefficients of all D 0 candidates (v sig+bkg n ) are determined using where the Q-vectors are defined as Q n ≡ ∑ j=1 w j e inϕ j , the asterisk symbol (*) means complex conjugation, and the subscripts A, B, and C refer to different subevents.In each event, the sum for Q nA and Q nB is over all hadrons detected in HF which are above a threshold energy, while the sum for Q nC includes all reconstructed tracks above a p T threshold.The Q D 0 n signifies the flow vector of a D 0 meson candidate within a particular kinematic range.The weight w j is a dimensionless quantity and corresponds to the energy deposited in the HF tower in GeV for Q nA and Q nB , or to the track p T in GeV/c, at azimuthal angle ϕ j for Q nC , or w j = 1 in the case of D 0 meson candidates.The subevent A (B) uses the negative (positive) side of the HF when the D 0 meson candidate is at positive pseudorapidity, and vice versa.This choice of subevents avoids autocorrelations and results in an η gap of at least three units between the D 0 meson daughters and particles from the underlying subevent, thereby suppressing short-range correlations.Flattening and recentering procedures are applied to the Q-vectors related to HF and the tracker, for removing detector acceptance effects [45,46].The averages ⟨Q nA Q * nB ⟩, ⟨Q nA Q * nC ⟩, and ⟨Q nB Q * nC ⟩ are found by considering all selected events, while the average ⟨Q D 0 n Q * nA ⟩ includes all D 0 meson candidates in all selected events.To extract the inclusive D 0 meson flow harmonics, the individual D 0 candidates are divided into bins of their flow vector SPs (i.e., Eq. ( 1) but not averaged over all D 0 candidates).A separate mass spectrum is generated in each of these SP bins and an invariant mass fit is performed to get the SP distribution of the D 0 signal.The inclusive D 0 v n values can be obtained by: where v i n is the center of the ith SP bin, and the Y i is the D 0 yield in the same bin.The former fit is used for determining the total D 0 yields and the latter for determining the fraction of nonprompt D 0 mesons.
The m inv distributions are fitted with five components, as shown in Fig. 1: the sum of two Gaussian functions with the same mean but different widths for the D 0 signal, S(m inv ); an additional Gaussian function to describe the invariant mass shape of those D 0 candidates that were given an incorrect mass assignment when swapping the pion and kaon designations, SW(m inv ); two Crystal Ball functions [47] to describe the processes D 0 → π + π − (S(m π + π − )) and D 0 → K + K − (S(m K + K − )); and a third-order polynomial to model the combinatorial background, Bkg(m inv ).
The contributions from the processes D 0 → π + π − and D 0 → K + K − are also resulting from using an incorrect K or π assignment.The ratio between the S(m inv ) and SW(m inv ) function yields is fixed according to values obtained from simulated events.The widths and means of the S(m inv ), SW(m inv ), and Crystal Ball functions are initialized using results obtained from simulation studies.To account for potential differences in momentum resolution between the model and the data, the widths are allowed to vary during the data fitting process.Specifically, a single scale factor, common to all three functions, is employed to adjust the component of the widths arising from momentum resolution.In addition, the means are also allowed to be different between the simulation and data.
The fit of the invariant mass spectrum gives the yield for inclusive D 0 mesons.For the extraction of the two individual components (prompt and nonprompt), a similar procedure as in Refs.[48,49] is followed.Distributions of the DCA between the D 0 meson trajectory and the primary vertex are fitted with a linear combination of prompt and nonprompt D 0 DCA templates obtained from MC in each bin of p T , centrality, and v i n .The data distribution was obtained by extracting the D 0 signal from fits to the invariant mass distribution in each DCA interval.In generating the simulated templates, the difference of the DCA values with and without detector effects for each D 0 meson is scaled.The same scale factor is chosen for both prompt and nonprompt D 0 mesons, with its value found by minimizing the χ 2 of the combined template fit.
The left panel of Fig. 1 shows an example of a fit to the mass spectrum for D 0 candidates in the p T interval 6-8 GeV/c for the centrality class 30-50% and the SP range 0.05 < v i 2 < 0.10.The right panel shows an example of the template fit of the inclusive D 0 meson yields, extracted as a function of DCA in the same D 0 kinematic region, but averaged over all v n .

Systematic uncertainties
Systematic uncertainties are estimated by varying the conditions and procedures of the v n measurements.The absolute differences between the nominal and the alternative analyses are considered as one root-mean-square deviations.In the fit to m inv , the functional form of the combinatorial background was varied using a second-order polynomial and an exponential function.In the 0-10% centrality bin, the uncertainty decreases with p T up to 5 GeV/c, remains constant from 5 to 10 GeV/c, and increases above 10 GeV/c.In more peripheral collisions, the uncertainty from this source increases monotonically with p T over the entire p T range.The mass shape for the signal region is varied by using a triple-Gaussian function and negligible differences are observed.The effect of the efficiency is studied by changing the efficiency correction values.In this alternative approach, the efficiency is based solely on simulation and does not consider any potential differences between the MC and data spectra.The largest effect is observed at low p T .The systematic uncertainty from the determination of the nonprompt D 0 meson fraction is evaluated by performing the template fit using an alternative variable, the DCA significance (DCA/(DCA uncertainty)).The differences observed between data and simulation for the DCA measurement and uncertainty are correlated, as they come from the momentum resolution of the same particles.Therefore, this effect is partially canceled in the DCA significance, obtained by dividing the DCA by its uncertainty.However, as DCA significance arises from the D 0 reconstruction, the widths of the distributions cannot be scaled to further match the data, in the same way it was done for DCA.The difference in results between the DCA and DCA significance template fits is quoted as a systematic uncertainty.This variation of the template fit gives the largest contribution to the systematic uncertainty, especially at the highest p T in central collisions.The effect of the BDT selection is estimated from the MC samples by comparing the results with and without applying it, in bins of v i n , and taking the difference as a systematic uncertainty.No p T dependence is observed for this source.Table 1 summarizes the estimates of the systematic uncertainties for v 2 and v 3 from each source.The total systematic uncertainty is obtained by adding individual uncertainties in quadrature.The nonprompt D 0 meson flow harmonics are presented in Fig. 2 together with values for prompt D 0 mesons from Ref. [18].The results, averaged in the range 2 < p T < 10 GeV/c, of v 2 = 0.0289 ± 0.0069, 0.0477 ± 0.0054, and 0.060 ± 0.010 in the centrality intervals 0-10%, 10-30%, and 30-50%, respectively, indicate nonzero values of elliptic flow for B mesons.The v 2 coefficient in the case of b hadron daughters has its maximum value at p T of ∼5 GeV/c and is significantly lower than in the prompt D 0 meson case.This difference becomes more pronounced going from central to peripheral collisions.These results for the mass ordering of collective flow extend the previously observed relationship between charm and light flavor hadrons, where the flow of charm hadrons has a lower magnitude than that for light flavor particles.Measurements also suggest an increase of v 2 towards peripheral collisions, as in the case of light hadrons.This observation agrees with the paradigm of flow as a consequence of initial space anisotropy [10].The elliptic flow of nonprompt D 0 mesons is found to be of the same magnitude as what the ATLAS Collaboration observed for nonprompt muons [19].
The nonprompt v 3 coefficient results have large statistical uncertainties so that neither the p T nor the centrality dependence can be determined.As was the case for v 2 , the triangular flow signal is weaker for nonprompt than for prompt D 0 mesons, but the difference is not as large as that for the elliptic flow.Measurements, averaged over all three centrality intervals, suggest positive v 3 for nonprompt D 0 mesons in the 4 < p T < 6 GeV/c range with a significance of 3.0 standard deviations.The ATLAS Collaboration v 3 results for b hadrons decaying into muons with p T > 4 GeV/c were found to be consistent with zero [19].Figure 3 shows the measured nonprompt D 0 meson v n coefficients compared with theoretical calculations that have different modeling of the b quark flow.The PHSD model is a microscopic off-shell transport model based on a Boltzmann transport equation approach, and includes only collisional energy loss [50].The TAMU model computes the space-time evolution of the heavy-quark phase space distribution in the QGP using the Fokker-Planck equation, implemented via Langevin dynamics, and also has no radiative energy loss [53].The LGR model uses the Langevin approach with gluon radiation, emphasizing the in-medium diffusion dynamics.This model aims to describe heavy-quark evolutions from low to intermediate p T range, p T < 15 GeV/c [54,55].The LBT model is based on a linearized Boltzmann approach coupled to a hydrodynamic background, and employs both collisional and radiative energy loss [51,52].The CUJET3.0 is essentially a jet energy loss framework based on a nonperturbative QGP medium [56,57].The LBT and CUJET3.0approaches are applicable only for the high-p T range and do not have predictions for lower momenta.The CUJET3.0 model implements hadronization through a parton fragmentation function, while all other models combine coalescence and fragmentation.All models can qualitatively describe the p T dependence of the data.The LGR model reproduces the measurements for the centrality 30-50% and in the low-and the intermediate-p T ranges, where b quarks are expected to follow Brownian motion, i.e., small random momentum fluctuations caused by collisions with thermal particles.For the central events, 0-10%, in the range 4 < p T < 10 GeV/c, the results from this model lie systematically below the data.In the same p T range, the PHSD model provides a good description of the data, for centralities 0-10% and 10-30%, but it cannot describe the data at higher centralities.In contrast, this same PHSD model underpredicts the v 2 Fourier coefficients for D 0 mesons that are not the product of b hadron decays [18].The TAMU model has predictions in the centrality range 20-40%, therefore only an indirect comparison with data is possible, but the predictions describe well the measurement performed in the centrality range 10-30%.At higher p T , where the anisotropy is driven by the path-length dependence of parton energy loss, the LBT model gives lower predictions than CUJET3.0,but both models match data within uncertainties.The PHSD, LBT, and LGR models account for event-by-event fluctuations of the initial geometry, but PHSD is the only model that has predictions for v 3 coefficients.While precision is limited for both the model and the data, they both reach similar maximum values in all centralities.

Summary
In summary, the elliptic (v 2 ) and triangular (v 3 ) flow harmonics of D 0 mesons that originate in b hadron decays (nonprompt D 0 mesons) are measured in lead-lead collisions at √ s NN = 5.02 TeV.The v 2 results show a weak transverse momentum (p T ) dependence and suggest a slight increase for more peripheral collisions.An indication of a nonzero v 3 coefficient is found for nonprompt D 0 mesons with 4 < p T < 6 GeV/c.The magnitudes of the flow coefficients are lower for nonprompt D 0 than for prompt D 0 mesons.This magnitude difference is more pronounced in the case of v 2 .The qualitative agreement between the results and theoretical models, which include a mass hierarchy in quark interactions with the quark-gluon plasma, extends our understanding of heavy quark interactions with the medium.

Figure 1 :
Figure 1: An example of the fit to the invariant mass spectrum (left panel) and an example of the template fit of the inclusive D 0 meson yields, extracted as a function of DCA (right panel).The former fit is used for determining the total D 0 yields and the latter for determining the fraction of nonprompt D 0 mesons.

Figure 2 :
Figure 2: The elliptic, v 2 (upper panels), and the triangular, v 3 (lower panels), flow coefficients of nonprompt and prompt (from Ref. [18]) D 0 mesons as functions of their p T and in three bins of centrality.The bars and the boxes represent statistical and systematic uncertainties, respectively.

Figure 3 :
Figure 3: The elliptic, v 2 (upper panel), and the triangular, v 3 (lower panel), flow coefficients of nonprompt D 0 mesons as functions of their p T and in three bins of centrality.The bars and the boxes represent statistical and systematic uncertainties, respectively.The colored bands show theoretical predictions [50-57].

Table 1 :
Summary of systematic uncertainties in the absolute differences between the nominal and the alternative analyses for v 2 and v 3 .