Studies of charm and beauty hadron long-range correlations in pp and pPb collisions at LHC energies

Measurements of the second Fourier harmonic coefficient ($v_2$) of the azimuthal distributions of prompt and nonprompt D$^0$ mesons (the latter coming from beauty hadron decays) produced in pp and pPb collisions are presented. The data samples are collected by the CMS experiment at nucleon-nucleon center-of-mass energies of 13 and 8.16 TeV, respectively. In high multiplicity pp collisions, $v_2$ signals for prompt charm hadrons are reported for the first time, and are found to be comparable to those for light-flavor hadron species over a transverse momentum ($p_\mathrm{T}$) range of 2-6 GeV. Compared at similar event multiplicities, the prompt D$^0$ meson $v_2$ values in pp and pPb collisions are similar in magnitude. The $v_2$ values for open beauty hadrons are extracted for the first time via nonprompt D$^0$ mesons in pPb collisions. For $p_\mathrm{T}$ in the range of 2-5 GeV, the results suggest that $v_2$ for nonprompt D$^0$ mesons are smaller than those for prompt D$^0$ mesons. These new measurements indicate a positive charm hadron $v_2$ in pp collisions and suggest a mass dependence in $v_2$ between charm and beauty hadrons in the pPb system. These results provide insights into the origin of heavy-flavor quark collectivity in small systems.


Introduction
Strong collectivity in high-energy nucleus-nucleus (AA) collisions at the BNL RHIC [1-4] and at the CERN LHC [5, 6], has indicated the formation of a hot, strongly interacting quark gluon plasma (QGP), which exhibits nearly ideal hydrodynamic behavior [7][8][9]. The collective phenomena manifests itself in long-range (large pseudorapidity gap) particle correlations [10][11][12][13][14][15]. Although not originally expected, similar long-range collective azimuthal correlations are also being observed in small colliding systems with high final-state particle multiplicity, such as proton-proton (pp) [16][17][18][19][20], proton-nucleus (pA) [21][22][23][24][25][26][27][28][29][30][31], and lighter nucleus-nucleus systems [31][32][33][34]. This observation raised the question of whether a fluid-like QGP medium with a size significantly smaller than in AA collisions is created in these other systems [35][36][37]. At the same time, there is no observation of long-range correlations in e + e − and ep collisions, which correspond to an even smaller system size compared to small hadronic collisions [38,39]. In the context of hydrodynamic models, the observed azimuthal correlation structure of emitted particles is typically characterized by its Fourier components [40]. The second and third Fourier anisotropy coefficients are known as elliptic (v 2 ) and triangular (v 3 ) flow, which most directly reflect the QGP medium response to the initial collision geometry and its fluctuations, respectively [41][42][43][44]. The experimental measurements in the small systems are consistent with the dominance of strong final-state interactions [35,37,[45][46][47], such as a hydrodynamic expansion of a tiny QGP droplet [35,37]. Alternative scenarios based on gluon saturation in the initial state can also capture the main features of the correlation data, and are conjectured to play a dominant role as the event multiplicity decreases [35,36].
Heavy-flavor quarks (charm and bottom) are produced via hard scatterings in the very early stages of the high energy collisions. These quarks are available to probe both initial-and finalstate effects of the collision dynamics [48,49]. In small colliding systems, the study of heavy-flavor hadron collectivity has the potential to disentangle possible contributions from both initial-and final-state effects. In particular, heavy flavor hadrons may be more sensitive to possible initial-state gluon saturation effects. Recent observation of a significant elliptic flow signal for prompt D 0 [57] and prompt J/ψ [58, 59] mesons in pPb collisions provided the first evidence for charm quark collectivity in small systems. Surprisingly, despite the mass differences, the observed v 2 signal for prompt J/ψ mesons in pPb collisions is found to be comparable to that of prompt D 0 mesons and light-flavor hadrons at a given particle transverse momentum (p T ). This behavior cannot be explained by the finalstate effects of a QGP medium, as the contribution from recombinations to J/ψ production is not expected to be significant in small systems [60]. This finding may imply the existence of initial-state correlation effects [61]. Further detailed investigations are important to address many open questions for understanding the origin of heavy-flavor quark collectivity in small systems. These include the multiplicity dependence of charm quark collectivity in both pPb and pp systems and the details of collective behavior of beauty quarks.
This Letter presents the first measurement of the elliptic flow (v 2 ) for prompt D 0 mesons in pp collisions at center-of-mass energy √ s = 13 TeV and for nonprompt D 0 mesons (from decays of beauty hadrons) in pPb collisions at nucleon-nucleon center-of-mass energy √ s NN = 8.16 TeV, using long-range two-particle angular correlations. The v 2 harmonic coefficient is determined over the 2-8 GeV p T range for prompt D 0 mesons as a function of multiplicity with results for the pp and pPb collisions. The nonprompt D 0 meson v 2 values are extracted in highmultiplicity pPb collisions for two transverse momentum ranges 2-5 and 5-8 GeV, and are compared to previous measurements of prompt D 0 mesons and light flavor hadrons.

Experimental apparatus and data sample
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume, there are four primary subdetectors including a silicon pixel and strip tracker detector, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Iron and quartz-fiber Cherenkov hadron forward calorimeters cover the pseudorapidity (η lab ) range 2.9 < |η lab | < 5.2 in laboratory frame. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. The silicon tracker measures charged particles within the range |η lab | < 2.5. For charged particles with 1 < p T < 10 GeV and |η lab | < 1.4, the track resolutions are typically 1.5% in p T and 25-90  µm in the transverse (longitudinal) impact parameter [62]. A 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. [63].
The event samples were collected by the CMS experiment with a two-level trigger system [64]: at level-1 events are selected by custom hardware processors while the high-level trigger uses fast versions of the offline software. The pPb data at √ s NN = 8.16 TeV used in this analysis were collected in 2016, and correspond to an integrated luminosity of 186.0 nb −1 [65]. The beam energies are 6.5 TeV for the protons and 2.56 TeV per nucleon for the lead nuclei. Because of the asymmetric beam conditions, particles selected in this analysis from midrapidity in the laboratory frame (|y lab | < 1) correspond to rapidity in the nucleon-nucleon center-of-mass frame of −1.46 < y cm < 0.54, with positive rapidity corresponding to the proton beam direction. The pp data at √ s = 13 TeV were collected in 2017 and 2018 with integrated luminosities of 1.27 pb −1 and 10.22 pb −1 during special runs with low beam intensity, resulting in an average number of concurrent pp collisions of about 1 per bunch crossing. The event reconstruction, event selections, and triggers (minimum bias and high multiplicity) are identical to those described in Refs. [19,66,67]. Similar to previous CMS correlation measurements, the pPb and pp data are analyzed for several multiplicity (N offline trk ) classes, where N offline trk is the number of offline selected tracks [19,62] with |η lab | < 2.4 and p T > 0.4 GeV.

Prompt and nonprompt D 0 meson reconstruction and selection
The D 0 (and its charge conjugate state D 0 ) mesons are reconstructed through the hadronic decay channel D 0 → K − π + (D 0 → K + π − ). In order to suppress the combinatorial background and improve the momentum and mass resolution, high-purity [62] tracks with p T > 0.7 GeV, |η lab | < 2.4, and relative uncertainty in p T < 10% are used. For each pair of selected tracks, two D 0 candidates are considered by assuming that one of the tracks has the pion mass while the other track has the kaon mass, and vice versa.
The D 0 candidates are selected based upon the following variables: their daughter charged particle track kinematics; the number of valid hits and relative p T uncertainties; the χ 2 probability of daughter tracks to originate from a common decay vertex; the three-dimensional distance (normalized or not by its uncertainty) between the primary and decay vertices; and the pointing angle (defined as the angle between the line segment connecting the primary and decay vertices and the momentum vector of the reconstructed particle candidates in the plane transverse to the beam direction). The selection is optimized separately for pp and pPb collisions, and for all p T ranges, using a multivariate technique that employs the boosted decision tree (BDT) algorithm [68], in order to maximize the statistical significance of the prompt or nonprompt D 0 meson signals. The Monte Carlo (MC) signal simulated samples are produced with PYTHIA 8.209 [69] tune CUETP8M1 [70] (embedded into EPOS LHC [71] for the case of pPb analysis) for both prompt and nonprompt D 0 events. The background samples for the multivariate training are taken from data. The training variables related to D 0 mesons include: the χ 2 probability for D 0 vertex fitting; the three-dimensional distance (with and without being normalized by its uncertainty) between the primary and decay vertices; and the three-dimensional pointing angle. The training variables related to the decay products are: p T ; pseudorapidity and the longitudinal and transverse track impact parameter significance. In the BDT training for prompt D 0 signals, same-sign (SS) π ± K ± candidates are used, which contain predominantly combinatorial background. For optimizing nonprompt D 0 signals, both prompt D 0 signals and combinatorial candidates are considered as dominant background to be suppressed. For this reason, opposite-sign (OS) candidates (although including fractions <5% of nonprompt D 0 signals) are used for the background training sample. This approach is found to give better performance for achieving higher nonprompt D 0 fractions than using SS background candidates, especially at higher p T .
The optimal selection criterion is the working point with the highest signal significance of prompt and nonprompt D 0 signals. For extracting the nonprompt D 0 yield, the distributions of distance of closest approach (DCA) of the D 0 meson momentum vector, relative to the primary vertex, are fitted using the template probability distribution functions (PDs) for prompt and nonprompt D 0 signals derived from MC simulation. The residual nonprompt fraction in the BDT prompt-trained sample is found to be no more than 7%, while in the BDT nonprompttrained sample, the optimal selection yields a nonprompt fraction up to 20%. This procedure is further outlined in Section 4.

Data analysis
The azimuthal anisotropies of D 0 mesons and strange hadrons are extracted from their longrange (|∆η| > 1) two-particle azimuthal correlations of D 0 candidates with charged particles, as described in Refs. [19,26]. Taking the D 0 meson as an example, the two-dimensional (2D) correlation function is constructed by pairing each D 0 candidate with reference primary chargedparticle tracks with 0.3 < p T < 3.0 GeV (denoted "ref" particles), and calculating where ∆η and ∆φ are the differences in pseudorapidity η lab and azimuthal angle φ of each pair. The same-event pair distribution, S(∆η, ∆φ), represents the yield of particle pairs normalized by the number of D 0 candidates (N D 0 ) from the same event. The mixed-event pair yield distribution, B(∆η, ∆φ), is constructed by pairing D 0 candidates in each event with the reference primary charged-particle tracks from 10 different randomly selected events, from the same N offline trk range, and with a primary vertex falling in the same 2 cm wide range of reconstructed z coordinates. The B(0, 0) represents the value of B(∆η, ∆φ) at ∆η = 0 and ∆φ = 0. It is evaluated by interpolating the four nearest bins with a bin width of 0.3 in ∆η and 1/16π in ∆φ bilinearly. The interpolation shows a negligible effect on the measurements. The analysis procedure is performed in each D 0 candidate p T range by dividing it into 14 intervals of invariant mass. The correction for the acceptance and efficiency (derived from simulations using PYTHIA for pp and PYTHIA+EPOS for pPb) of the D 0 meson yield is found to have a negligible effect on the measurements, and is not applied. The corresponding effects are discussed in Section 5. The ∆φ correlation functions averaged over |∆η| > 1 (to remove short-range correlations, such as jet fragmentation) are obtained from the projection of 2D correlation functions and fitted by the first three terms of a Fourier series: 2V n∆ cos(n∆φ) . (2) Here, V n∆ are the Fourier coefficients and N assoc represents the total number of pairs per D 0 candidate. The inclusion of additional Fourier terms to the fit has negligible effect. By assuming V n∆ to be the product of single-particle anisotropies [72], , the v n anisotropy harmonics for D 0 candidates can be extracted from the equation: Because of the limited statistical precision of the available data, only the elliptic anisotropy harmonic results are reported in this analysis.
To extract the V 2∆ values of the inclusive D 0 meson signal (V S 2∆ ), a two-step fit to the invariant mass spectrum of D 0 candidates and their V 2∆ as a function of the invariant mass, The mass spectrum fit function is composed of five components: 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 D 0 candidates with an incorrect mass assignment from the exchange of the pion and kaon designations, SW(m inv ); Crystal Ball (CB) functions [73] to describe processes D 0 → π + π − (S(m π + π − )) and D 0 → K + K − (S(m K + K − )); and a third-order polynomial to model the combinatorial background, B(m inv ). The contributions from the processes D 0 → π + π − and D 0 → K + K − are the results of mislabelling K as π, or vice versa. These two components are emulated by two CB functions at two sides away from the peak region. The width and the ratio of the yields of SW(m inv ) and S(m inv ) and the CB function shape are fixed according to results obtained from simulation studies using PYTHIA for pp collisions and PYTHIA+EPOS for pPb collisions. where Here V B 2∆ (m inv ) for the background D 0 candidates is modeled as a linear function of the invariant mass, and α(m inv ) is the D 0 signal fraction. The K-π swapped, D 0 → π + π − and D 0 → K + K − components are included in the signal fraction because these candidates are from genuine D 0 mesons and should have the same v 2 value as that of the D 0 signal. Figure 1 shows an example of fits to the mass spectrum and V S+B 2∆ (m inv ), for the BDT prompttrained sample in the p T interval 4-6 GeV for the multiplicity range N offline trk ≥ 100 in pp collisions. Similar fits in pPb data can be found in Ref.
[57], which are not repeated here. For extracting the V 2∆ values of nonprompt D 0 mesons, the measurement and fitting procedure described above are repeated in three separate DCA ranges, containing very different nonprompt D 0 fractions. A linear fit by the functional form, to the measured D 0 V 2∆ values as a function of nonprompt D 0 fraction is performed to extrapolate to the V 2∆ value at a nonprompt fraction of 100%. The f b→D represents the nonprompt D 0 fraction. The v 2 values of nonprompt D 0 are evaluated by using Eq. (3). Figure 2 shows an example of fits to the mass spectrum and V S+B 2∆ (m inv ) for the BDT nonprompt-trained sample for DCA < 0.008 cm and 0.008 < DCA < 0.014 cm, in the p T interval 2-5 GeV, for the multiplicity  Figure 2: Example of fits to the invariant mass spectrum and V S+B 2∆ (m inv ), for the BDT nonprompt-trained sample in pPb collisions. The left plot shows the fit for DCA < 0.008 cm and the right plot is for 0.008 < DCA < 0.014 cm.
range 185 ≤ N offline trk < 250 in pPb collisions. The resulting D 0 signal V 2∆ distributions contain contributions from both prompt and nonprompt D 0 mesons.
Inclusive D 0 meson yields, extracted as a function of DCA, by fitting the invariant mass distribution in each DCA bin, are shown in Fig. 3 (left). A template fit to the DCA distribution is performed using template distributions of prompt and nonprompt D 0 mesons obtained from MC simulation to estimate the nonprompt D 0 fractions in each of the three DCA regions used to extract inclusive D 0 V 2∆ , as described above. The inclusive D 0 V 2∆ values from the three DCA regions are then plotted as a function of the corresponding nonprompt D 0 fraction, shown in Fig. 3 (right), for 2 < p T < 5 GeV and 5 < p T < 8 GeV, respectively. The measurements are well described by a linear-function fit, which is shown as a red line in Fig. 3.
The residual contribution of back-to-back dijets to the measured v 2 results is corrected by subtracting correlations from low-multiplicity events, following an identical procedure established in Refs. [19,72]. The Fourier coefficients, V n∆ , extracted from Eq. (2) for N offline trk < 35(20), in pPb (pp) collisions, are subtracted from the V n∆ coefficients obtained in the high-multiplicity region, with Here, Y jet represents the jet yield. It is the difference between integrals of the short-range (|∆η| < 1) and long-range (|∆η| > 2) event-normalized associated yields for each multiplicity class. The ratio Y jet /Y jet (N offline trk < 35) is introduced to account for the enhanced jet correlations resulting from the selection of higher-multiplicity events. It is observed that the values of jet yield ratio show little dependence on p T over the full p T range. For the measurement of nonprompt D 0 mesons, all quantities in Eq. (7) are first extrapolated to values at a nonprompt D 0 fraction of 100%, following the same approach as in Fig. 3, before applying the subtraction procedure. Elliptic flow (v sub 2 ), corrected for residual jet correlations, is obtained from V sub 2∆ using Eq. (3).  Other sources of systematic uncertainty include the background mass PD, the D 0 meson yield correction (acceptance and efficiency correction), the background V 2∆ PD, and the jet subtraction method. Changing the background mass PD to a second-order polynomial or an exponential function shows negligible systematic effects. To evaluate the uncertainties arising from the p T -dependent D 0 meson yield correction, the v 2 values are extracted from the corrected signal D 0 distributions and compared to the uncorrected v 2 values as a conservative estimate. This yields an uncertainty of less than 0.013. The systematic uncertainties from the background v 2 PD are evaluated by changing v B 2 (m inv ) to a second-order polynomial function of the invariant mass, yielding an uncertainty of less than 0.005. To study potential trigger biases, a comparison to high-multiplicity pPb data for a given multiplicity range that were collected using a lower threshold trigger with 100% efficiency is performed. The uncertainty from trigger bias is quoted as 0.001. Though data collected with low beam intensity are used in this analysis, there are still additional collisions besides the one of interest per bunch crossing, which are known as For the measurement of prompt D 0 mesons, the contribution from nonprompt D 0 mesons is significantly suppressed. No explicit correction is applied and a systematic uncertainty is quoted instead. Based on the prediction for AA collisions that B mesons have a smaller v 2 than lightflavor particles because of the larger mass of the b quark [74][75][76], the nonprompt D 0 v 2 values are assumed to lie between 0 and those of strange hadrons. The v 2 for prompt D 0 is thus reestimated with the bounds of nonprompt D 0 v 2 and the extracted nonprompt D 0 fractions and the change in v 2 signal is found to be smaller than 0.008. For the measurement of nonprompt D 0 mesons, a major systematic uncertainty comes from the determination of nonprompt D 0 fraction in different DCA regions. The DCA template distributions of prompt and nonprompt D 0 mesons from MC simulation are smeared via scaling the width of these distributions. The variation of DCA width is 2-8%, based on the best χ 2 fit to data. The resulting variation in the extracted nonprompt D 0 v 2 are quoted as a systematic uncertainty of 0.007.

Systematic uncertainties
All sources of systematic uncertainties are added in quadrature to obtain the total systematic uncertainty. The total systematic uncertainties for prompt and nonprompt D 0 mesons in pPb collisions yield 0.005-0.018 and 0.016-0.017, respectively. For prompt D 0 mesons in pp collisions, the total systematics uncertainties are quoted as 0.013-0.052.

Results
The v sub 2 results of prompt D 0 mesons in pp collisions at √ s = 13 TeV, are presented in Fig. 4 as a function of p T for |y lab | < 1, with N offline trk ≥ 100. Published data for light-flavor hadrons including inclusive charged particles (dominated by pions), K 0 S mesons and Λ baryons are also shown for comparison [19]. The positive v 2 signal (0.061 ± 0.018 (stat) ± 0.013 (syst)) over a p T range of ∼2-4 GeV for prompt charm hadrons provides indications of the collectivity of charm quarks in pp collisions, with a declining trend toward higher p T . The v 2 magnitude for prompt D 0 mesons is found to be compatible with light-flavor hadron species, though slightly smaller by about one standard deviation. The results suggest that collectivity is being developed for charm hadrons in pp collisions, comparable (or slightly weaker) than that for light-flavor hadrons. This finding is similar to the observation made in pPb collisions at To further investigate possible system size dependence of collectivity for charm hadrons in small colliding systems, v 2 for prompt D 0 mesons in pPb and pp collisions are both measured in different multiplicity classes. The prompt D 0 v 2 as a function of event multiplicity for three different p T ranges: 2 < p T < 4 GeV, 4 < p T < 6 GeV, and 6 < p T < 8 GeV are presented in Fig. 5. At similar multiplicities of N offline trk ∼ 100, the prompt D 0 v 2 values are found to be comparable within uncertainties in pp and pPb systems. For 2 < p T < 4 GeV, the measured results of prompt D 0 provide indications of positive v 2 down to N offline trk ∼ 50 with a significance of more than 2.4 standard deviations in pPb collisions, while for 6 < p T < 8 GeV the prompt D 0 v 2 signal tends to diminish in the low multiplicity regions. No clear multiplicity dependence can be determined for pp data, because of large statistical uncertainties at low multiplicities.
The v sub 2 results for nonprompt D 0 mesons from beauty hadron decays are shown in Fig. 6 as a function of p T for pPb collisions at 8.16 TeV with 185 ≤ N offline trk < 250. The extracted v sub 2 values are −0.008 ± 0.028 (stat) ± 0.016 (syst) for 2 < p T < 5 GeV and 0.057 ± 0.029 (stat) ± 0.017 (syst) for 5 < p T < 8 GeV. At low p T , the nonprompt D 0 v 2 is consistent with zero, while at high p T , a hint of a positive v 2 value for beauty mesons is suggested but not significant within statistical and systematic uncertainties. Previously published v 2 data for prompt D 0 mesons and strange hadrons are also shown [57].
At p T ∼ 2-5 GeV, the nonprompt D 0 meson v 2 from beauty hadron decays is observed to be smaller than that for prompt D 0 mesons with a significance of 2.7 standard deviations. And nonprompt D 0 mesons carry >50% of B transverse momenta based on MC simulations using PYTHIA 8.209 [69] tune CUETP8M1 [70]. These studies suggest a flavor hierarchy of the collectivity signal that tends to diminish for the heavier beauty hadrons. This indication is qual- itatively consistent with the scenario of v 2 being generated via final-state rescatterings, where heavier quarks tend to develop a weaker collective v 2 signal [49].
Correlations at the initial stage of the collision between partons originating from projectile protons and dense gluons in the lead nucleus are able to generate sizable elliptic flow in the color glass condensate (CGC) framework [35, 61,77]. These CGC calculations of v 2 signals for prompt J/ψ mesons, as well as prompt and nonprompt (from B meson decay) D 0 mesons, are compared with data in Fig. 6. The qualitative agreement between data and theory suggest that initial-state effects may play an important role in the generation of collectivity for these particles in pPb collisions. The CGC framework also predicts a flavor hierarchy between prompt and nonprompt D 0 for p T ∼ 2-5 GeV, again consistent with the data within uncertainties.

Summary
The first measurements of elliptic azimuthal anisotropies for prompt D 0 mesons in protonproton (pp) collisions at center-of-mass energy √ s = 13 TeV, and for nonprompt D 0 mesons from beauty hadron decays in proton-lead (pPb) collisions at nucleon-nucleon center-of-mass energy √ s NN = 8.16 TeV are presented. In pp collisions with multiplicities of N offline trk ≥ 100, the second Fourier harmonic coefficient (v 2 ) of the azimuthal distributions for prompt D 0 mesons are measured over the transverse momentum (p T ) range of 2-8 GeV, with indications of positive v 2 signals over the p T range of 2-4 GeV. These values are found to be comparable (or slightly smaller) to those of light-flavor hadron species. At similar event multiplicities, the prompt D 0 meson v 2 signals in pp and pPb collisions are found to be comparable in magnitude. The v 2 values of open beauty hadrons are extracted for the first time via non-prompt D 0 mesons in pPb collisions, with magnitudes smaller than those for prompt D 0 mesons for p T ∼ 2-5 GeV. The new measurements of charm hadron v 2 in the pp system and the indications of mass dependence of heavy-flavor hadron v 2 in the pPb system provide insights into the origin of heavy-flavor quark collectivity in small colliding systems.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses.     [11] STAR Collaboration, "Long range rapidity correlations and jet production in high energy nuclear collisions", Phys.  [60] X. Du and R. Rapp, "In-medium charmonium production in proton-nucleus collisions", JHEP 03 (2019) 015, doi:10.1007/JHEP03(2019)015, arXiv:1808.10014.