Measurement of D$^0$-meson + hadron two-dimensional angular correlations in Au+Au collisions at $\sqrt{s_{\rm NN}} = $ 200 GeV

Open heavy flavor hadrons provide unique probes of the medium produced in ultra-relativistic heavy-ion collisions. Due to their increased mass relative to light-flavor hadrons, long lifetime, and early production in hard-scattering interactions, they provide access to the full evolution of the partonic medium formed in heavy-ion collisions. This paper reports two-dimensional (2D) angular correlations between neutral $D$-mesons and unidentified charged particles produced in minimum-bias Au+Au collisions at $\sqrt{s_{\rm NN}}$ = 200 GeV. $D^0$ and $\bar{D}^0$ mesons are reconstructed via their weak decay to $K^{\mp} \pi^{\pm}$ using the Heavy Flavor Tracker (HFT) in the Solenoidal Tracker at RHIC (STAR) experiment. Correlations on relative pseudorapidity and azimuth $(\Delta\eta,\Delta\phi)$ are presented for peripheral, mid-central and central collisions with $D^0$ transverse momentum from 2 to 10 GeV/$c$. Attention is focused on the 2D peaked correlation structure near the triggered $D^0$-meson, the {\em near-side} (NS) peak, which serves as a proxy for a charm-quark containing jet. The correlated NS yield of charged particles per $D^0$-meson and the 2D widths of the NS peak increase significantly from peripheral to central collisions. These results are compared with similar correlations using unidentified charged particles, consisting primarily of light-flavor hadrons, at similar trigger particle momenta. Similar per-trigger yields and widths of the NS correlation peak are observed. The present results provide additional evidence that $D^0$-mesons undergo significant interactions with the medium formed in heavy-ion collision and show, for the first time, significant centrality evolution of the NS 2D peak in the correlations of particles associated with a heavy-flavor hadron produced in these collisions.


Au+Au collisions at
√ sNN = 200 GeV. D 0 andD 0 mesons are reconstructed via their weak decay to K ∓ π ± using the Heavy Flavor Tracker (HFT) in the Solenoidal Tracker at RHIC (STAR) experiment. Correlations on relative pseudorapidity and azimuth (∆η, ∆φ) are presented for peripheral, mid-central and central collisions with D 0 transverse momentum from 2 to 10 GeV/c. Attention is focused on the 2D peaked correlation structure near the triggered D 0 -meson, the near-side (NS) peak, which serves as a proxy for a charm-quark containing jet. The correlated NS yield of charged particles per D 0 -meson and the 2D widths of the NS peak increase significantly from peripheral to central collisions. These results are compared with similar correlations using unidentified charged particles, consisting primarily of light-flavor hadrons, at similar trigger particle momenta. Similar per-trigger yields and widths of the NS correlation peak are observed. The present results provide additional evidence that D 0 -mesons undergo significant interactions with the medium formed in heavy-ion collision and show, for the first time, significant centrality evolution of the NS 2D peak in the correlations of particles associated with a heavy-flavor hadron produced in these collisions.

I. INTRODUCTION
Heavy-flavor (HF) quark (charm and beauty) production in relativistic heavy-ion collisions provides a unique probe of the produced deconfined partonic matter. This is because the energy scale (3 GeV for charm production) is sufficiently large such that the production mechanisms can be calculated from perturbative QCD (e.g. FONLL [1][2][3]). The charm quark contained in the finalstate particle is very likely produced in the initial collision stages [1,4]. The charm quark and/or charmed hadron can therefore access the many-body QCD dynamics in the very early collision stage when the partonic density is largest. This enables experimental studies of (1) heavy flavor quark or hadron interactions in the medium, (2) medium modifications of heavy flavor quark fragmentation [5] and hadronization, and (3) dissociation mechanisms of HF mesons with hidden flavor (e.g. J/ψ) [3]. Experimentally, open HF mesons can be indirectly observed via semi-leptonic decay modes to single electrons, or directly via weak decay channels, e.g. D 0 → K − π + andD 0 → K + π − . The latter two decay channels are used in the present analysis, and D 0 will be used to represent both D 0 andD 0 throughout this paper.
In recent years, HF yields in the form of transverse momentum (p T ) spectra and nuclear modification factor R AA [6][7][8], and azimuthal anisotropy amplitude v 2 [9][10][11] have been reported from relativistic heavy-ion collision experiments at the Relativistic Heavy Ion Collider (RHIC) and at the Large Hadron Collider (LHC). Additionally, HF correlations using D-mesons have been studied by the ALICE collaboration in p+p and p+Pb collisions, and the results have been shown to be consistent with PYTHIA [12]. While overall charm quark production (cc) follows binary nucleon + nucleon collision scaling, the yields at p T > 2 GeV/c are suppressed in heavy-ion collisions, relative to binary scaling expectations as seen in the measurements of open-charm hadron R AA [6][7][8]. Similar amounts of suppression are also observed for light flavor (LF) meson production. In addition, D 0 -meson v 2 , as a function of transverse kinetic energy, is also similar to LF results and is consistent with the number of constituent quark scaling observed for LF hadrons [9][10][11] . Both results suggest that the charm quark or meson interacts with the medium similarly to LF quarks or mesons. For HF quarks, QCD predicts less radiative energy loss than for low mass quarks ("dead cone effect") [13] due to the suppression of gluon radiation at forward angles, below M quark /E quark [14]. Collisional energy loss scales with inverse mass [14], suggesting further reduction in interaction effects relative to LF mesons. The surprisingly strong charm quark interaction effects implied by measurements of open-charm R AA and v 2 therefore motivate additional measurements and new studies to better understand these interactions, e.g. measurements of possible long-range correlations on pseudorapidity.
Transverse momentum integrated, two-dimensional (2D) angular correlations of unidentified chargedparticles from minimum-bias Au+Au collisions at RHIC energies have been measured by the STAR experiment [15]. These correlations exhibit a sudden onset, starting near mid-central collisions, of an increase in the per-trigger amplitude and width along relative pseudorapidity of the near-side (NS), 2D correlation peak. Similarly, this elongated correlation structure, commonly referred to as the "ridge," 1 has been reported in p Tselected, trigger-associated correlations in Au+Au collisions [16,17] and in Pb+Pb collisions at the LHC [18][19][20]. The primary goals of the present analysis are to study the centrality dependence of 2D angular correlations of D 0 -mesons plus associated charged hadrons, and to determine if the ridge phenomenon also occurs for HF mesons.
In general, the HF R AA and v 2 data from RHIC and the LHC have been described by a variety of models. The principle physics issues considered in the models include: (1) description of the initial-state including shadowing and saturation [21]; (2) the rapid HF formation time restricting the HF parton cascade [22]; (3) radiative, collisional (diffusion), and dissociative interac-tions [3,14]; (4) longitudinal color-field (glasma) effects on HF fragmentation [23]; (5) HF hadronization based on fragmentation, recombination or a mixture of both [24][25][26][27][28]. Transport models [2,26,27,29] or stochastic transport of the HF quark or hadron within a hydrodynamic medium [25,28,30] are typically assumed. By reporting an experimental quantity, other than R AA and v 2 , which gives access to charm-jet and flow-related physics simultaneously, more sensitivity to the many-body, nonperturbative QCD interactions is possible. Future theoretical analyses of the data presented here may lead to a better understanding of those interactions and of the medium itself.
In the present work we report 2D angular correlations on relative pseudorapidity ∆η = η D 0 − η h ± and relative azimuthal angle ∆φ = φ D 0 −φ h ± of charged-hadrons with D 0 andD 0 mesons produced in minimum-bias Au+Au collisions at √ s NN = 200 GeV. This type of analysis permits the (∆η, ∆φ) dependences of the correlations to be separated, allowing any possible ∆η-dependent correlations to be distinguished from the ∆η-independent azimuthal harmonics, such as elliptic flow. In this paper, we focus on the centrality evolution of the 2D angular distribution and number of associated charged hadrons on the NS within |∆φ| ≤ π/2 which have become correlated with the triggered D 0 -meson as the charm quark fragments or the D 0 -meson propagates through the medium. This paper is organized as follows. The analysis method is described in Sec. II and the data processing steps and other details are discussed in Sec. III and in two appendices. The correlation data and fitting results are presented in Sec. IV. Systematic uncertainties are discussed in Sec. V and our results are further discussed and compared with predictions from PYTHIA and LF correlations measured by STAR in Sec. VI. A summary and conclusions from this work are given in Sec. VII.

II. ANALYSIS METHOD
In conventional two-particle correlation analyses with charged particles, the number of particle pairs from the same-event (SE) in each (∆η, ∆φ) bin is summed for a collection of collision events (e.g. in a centrality class). The particle pairs are primarily uncorrelated, but do include correlation effects. The uncorrelated pair background can be estimated by similarly counting pairs of particles where the two particles in each pair are sampled from different events (mixed-events, ME). The mixed-events are required to have similar multiplicities and primary collision vertex positions along the beam axis. These quantities are defined as ρ SE (∆η, ∆φ) and ρ ME (∆η, ∆φ), respectively, which are pair densities (number of pairs per (∆η, ∆φ) bin area). The normalized difference, ∆ρ ≡ ρ SE − αρ ME , approximates the two-particle correlation in each (∆η, ∆φ) bin, where α is the ratio of the total number of SE pairs to ME pairs, α ≡ N pairs,SE /N pairs,ME , which normalizes ρ ME to the same overall scale as ρ SE . Both pair counts are affected by detector acceptance and particle reconstruction inefficiency. These effects can be corrected in each (∆η, ∆φ) bin by dividing ∆ρ by the normalized ME distribution, αρ ME . The ratio ∆ρ/αρ ME (∆η, ∆φ) [see Eq. (1) in Ref. [15]] is the underlying acceptance and efficiency corrected correlation quantity from which other quantities, e.g. per-trigger correlations, can be derived.
For D 0 -meson plus charged-particle (D 0 + h ± ) correlations, where the short-lived D 0 cannot be detected directly, only statistical reconstruction is possible because the number of produced D 0 -mesons must be inferred from the invariant mass distribution constructed from the daughter particle momentum vectors. The cleanest decay channel for D 0 reconstruction is the weak decay to unlike-sign K ∓ π ± pairs (B.R = 3.93%) [31]. Random, combinatoric K ∓ π ± pairs which pass all the track and vertex reconstruction cuts can be drastically reduced based on the optimized secondary decay vertex parameters. These parameters are reconstructed using particle trajectories (tracks) measured in the STAR Time Projection Chamber (TPC) [32] and Heavy Flavor Tracker (HFT) [33].
Even so, background pairs remain and it is impossible to distinguish random K ∓ π ± pairs from true D 0 decay daughters. Correlations between those random Kπ pairs and other charged hadrons must be accounted for and removed in the analysis. These background correlations can be estimated by correlating random Kπ pairs from side-bands (SB) in the Kπ invariant mass distribution, with other charged particles in the event, as further explained in Sec. III. Unless stated explicitly, both K ∓ π ± pairs are included and denoted simply as Kπ.
In addition, a significant fraction of the D 0 -mesons (∼17% [31]) are produced from decays of the D ⋆± resonance. Correlations between these daughter D 0 -mesons and other charged particles are of physical interest because the daughter D 0 -meson contains the original cquark and most of the parent D ⋆± momentum. However, the decay length of the D ⋆± (∼0.12 nm) allows the decay channel D ⋆± → D 0 + π soft (B.R. = 67.7%, π soft are lowmomentum soft pions) to occur well outside the medium. The resulting D 0 +π soft correlation reflects only the decay kinematics and is considered a background correlation which must also be removed from the measurements.
The true D 0 + h ± correlations are calculated by subtracting the random Kπ + h ± background correlations and the preceding correlated D 0 + π soft pairs from the measured quantity ∆ρ/αρ ME , where the latter uses all Kπ pairs in the D 0 signal range of the Kπ invariant mass distribution. The basic correlation quantity, derived in Appendix A, is given by where S and B are the deducd D 0 signal and background yields from the Kπ invariant mass distribution near the D 0 mass, as explained in the next section. Subscripts D 0 + h, sig, SB, and D 0 + π soft indicate the pair or Kπ invariant mass region used. The quantity on the left-hand side of Eq. (1) is symbolic only and must be determined by the measured quantities on the right-hand side. Other technical details for these calculations are explained in the next section.
III. DATA AND TECHNICAL DETAILS A. Event and particle selection The data for this analysis were collected by the STAR experiment [34] at RHIC during the 2014 run period. The dataset includes approximately 900 million minimumbias Au+Au √ s N N = 200 GeV collision events for which coincident signals between the two, symmetrically positioned Vertex Position Detectors (VPD) [35] were required. In addition, the primary collision vertex (PV) for each accepted event was required to be within the tracking fiducial region of the HFT [33], ±6 cm along the beam axis (z-axis), to ensure uniform HFT acceptance. Charged particle trajectories were initially reconstructed using the TPC in the presence of a 0.5 T uniform magnetic field parallel with the beam axis. All used tracks were required to have at least 20 reconstructed space points ("hits") (out of a possible 45) in the TPC, a ratio of the number of found hits to the maximum number expected > 0.52 to remove split tracks, a Kalman filter least-squares fitted χ 2 /NDF < 3, and a distance of closest approach (DCA) to the PV < 3 cm. All tracks used in the analysis must also fall within the acceptance of the STAR TPC: p T > 0.15 GeV/c, |η| ≤ 1 (pseudorapidity), and full 2π in azimuth. This sub-set of reconstructed charged particles in the TPC are referred to as TPC tracks and were used to determine collision centrality (see next sub-section). In addition, all tracks used to construct the correlations were required to match to at least one hit in each of the inner three layers of the HFT [33], including two in the silicon pixel detector (PXL) layers, and one in the Intermediate Silicon Tracker (IST) [36]. The spatial resolution of projected tracks near the PV is greatly improved, e.g. from ∼ 1 cm for TPC tracks to ∼ 30 µm, for p T > 1.5 GeV/c, when HFT hits are included. Furthermore, the fast-timing of the IST suppresses track pileup contamination from collisions occurring before or after the triggered collision.
The trigger D 0 andD 0 were reconstructed via the hadronic decay channel D 0 → K + π (cτ = 123 µm for D 0 ). The secondary decay vertices are reconstructed using the above TPC+HFT tracks and an optimized set of limits on the accepted decay topology parameters (see Fig. 1). The parameter limits were taken from a previous STAR D 0 analysis [9] after adjusting for the 2 < p T,D 0 < 10 GeV/c range used in the present analysis. The limits for the five geometrical cuts in Fig. 1 − the decay length, DCA between daughters, DCA of the reconstructed D 0 to the PV, DCA of the pion daughter to the PV, and DCA of kaon daughter to the PV are, respectively, > 212 µm, < 57 µm, < 38 µm, > 86 µm, and > 95 µm. In addition, pion and kaon identification requirements based on measured ionization energy loss (dE/dx) in the TPC were imposed on the D 0 decay daughter candidates. Those cuts required that the fitted dE/dx for the assigned space-points be < 2σ from the expected mean.
The invariant mass distribution was constructed for all unlike-sign (US) and like-sign (LS) Kπ pairs, where the LS approximates the invariant mass background. The LS distribution was normalized to the US data in the invariant mass range from 2.0 < M Kπ < 2.1 GeV/c 2 and subtracted from the US distribution. The remaining background was fit with a straight-line and then subtracted. Other functional shapes (e.g. exponential, polynomial) were tested and the resulting variations in the signal (S) and background (B) yields were found to be small (< 1%). The systematic effect on the correlations due to the extraction of the S and B yields are discussed in Sec.V. The results of this procedure are summarized in Fig. 2. The signal and background yields were calculated using bin counting in the invariant mass distribution in the range 1.82 < M Kπ < 1.90 GeV/c 2 , where S+B was calculated using the raw, unlike-sign distribution, and S was calculated from the fully-subtracted distribution. All trigger D 0 s used in the present analysis were restricted within 2 < p T,D 0 < 10 GeV/c to maximize statistical significance. Correlations constructed with p T,D 0 < 2 GeV/c exhibited large fluctuations in the correlation structures for small changes in the topological cuts, far beyond statistical uncertainty, and were therefore excluded from the analysis. Residual structure in the invariant mass background below 1.7 GeV/c 2 was found to be predominately from other D 0 -meson decays [9].

B. Collision centrality determination
The minimum-bias event sample was divided into three centrality classes, using the observed event-wise number of TPC tracks with |η| ≤ 1 and p T ≥ 0.15 GeV/c according to the method in Ref. [15]. The measured multiplicity frequency distribution was corrected for TPC tracking efficiency, thus determining TPC track multiplicity limits corresponding to centrality fractions 0-20%, 20-50% and 50-80% of the total reaction cross section. Centrality was based on multiplicities within |η| ≤ 1, instead of |η| ≤ 0.5, in order to avoid significant artifacts in the angular correlations along ∆η as explained in Ref. [15]. Additional corrections due to small (few percent) variations in TPC tracking efficiency as functions of PV position and beam + beam run-time luminosity were negligible. C. Construction of pair histograms D 0 -candidate + associated charged-hadron pairs from the same event are formed on binned coordinates (∆η, ∆φ) and summed over all events in each centrality class. Particles used as D 0 daughter-candidates are excluded from the associated hadrons. The D 0 trigger is defined to be any reconstructed Kπ pair passing all the above cuts and falling within the invariant mass range 1.82 < M Kπ < 1.90 GeV/c 2 . This selection includes both real D 0 s as well as combinatorial background Kπ pairs. To estimate the correlations from this background, pair histograms on (∆η, ∆φ) are also constructed using Kπ pairs from two side-band regions in the invariant mass spectrum defined by: left side-band, 1.7 < M Kπ < 1.8 GeV/c 2 ; right side-band, 1.92 < M Kπ < 2.10 GeV/c 2 . The different widths of the left and right side-bands are chosen to use approximately the same yield of background Kπ pairs from each side-band. An efficiency correction (weight) is applied to each individual Kπ + associated hadron pair. The pair-weight is defined as where ǫ D 0 , ǫ K , ǫ π and ǫ h are the individual reconstruction efficiencies for the D 0 , K, π and charged-hadron, respectively. Overbars in Eq.(2) denote averages. Ratios S/(S + B) and B/(S + B) are the probabilities that the candidate Kπ pair is actually from a D 0 decay or is random, respectively. In the side-bands all Kπ pairs are considered random. The K, π and charged-hadron TPC tracking efficiencies are taken from the analysis in Ref. [37]. Those efficiencies are then multiplied by the additional p T -dependent, TPC+HFT track matching ratio to get the quantities ǫ K , ǫ π and ǫ h . The D 0 efficiency is the ratio of the raw yield of D 0 mesons as a function of p T , using the above cuts, to the published invariant yield [6]. The overall shape of the above Kπ + hadron pair distribution is dominated by the finite pseudorapidity acceptance which introduces an approximate triangular shape on ∆η. This overall shape plus any other acceptance artifacts caused by the TPC sector edges, electronics outages, etc., can be corrected by dividing, bin-bybin, a similarly constructed mixed-event distribution as explained in Sec. II. Accurate acceptance corrections require that the PV location of each pair of mixed-events are sampled within sufficiently narrow sub-bins along the beam-axis, where the 12 cm range was divided into 10 uniform sub-bins. Correlation artifacts can also occur when the centralities of the mixed-events differ too much. Restricting the range of multiplicties in the eventmixing sub-bins to < 50, assuming 2 units in η, was previously shown to be sufficient [15,38,39]. The latter resulted in 16 multiplicity sub-bins, for a total of 160 event-mixing sub-bins in this analysis. Mixed-event distributions were constructed for each of the three Kπ invariant mass ranges discussed above. Efficiency corrections were also applied to each mixed-event Kπ + hadron pair as in Eq. (2).

D. Symmetrization on ∆η and ∆φ
For the present analysis with identical, unpolarized colliding ions, where particles within a symmetric pseudorapidity range centered at η = 0 are used, we may project the above pair histograms onto absolute value binned coordinates (|∆η|, |∆φ|) without loss of information. Pair counts and statistical errors can then be copied to corresponding bins in the other quadrants in the full (∆η, ∆φ) space for visual display. An odd number (13) of uniform ∆η bins within −2 ≤ ∆η ≤ 2 were used and multiples of four ∆φ bins (12) were assumed within full 2π, where ∆φ bins are centered at ∆φ = 0 and π. Bins centered at (∆η, ∆φ) = (0, 0) and (0, π) therefore have approximately 1/4 the number of pairs as nearby bins centered at non-zero (∆η, ∆φ). Other bins centered at either ∆η = 0, ∆φ = 0, or ∆φ = π similarly have ∼ 1/2 the number of pairs. Statistical errors in these bins are approximately factors of 2 and √ 2 larger, respectively, than that in neighboring non-zero (∆η, ∆φ) bins. Statistical errors in the final correlations in Eq. (1) are determined by the SE and ME pair counts in each (|∆η|, |∆φ|) bin, including the error contributions from the D 0 signal region and the two side-bands of the Kπ invariant mass distribution, assuming uncorrelated uncertainties. Absolute statistical uncertainties in the correlations for each centrality class are approximately ±0.0095 for 50-80%, ±0.004 for 20-50%, and ±0.002 for the 0-20% centrality class. Errors generally increase by almost a factor of two at the outermost bins on ∆η.

E. D ⋆± correction
In Sec. II, the background contribution of D ⋆± decays occurring outside the medium to D 0 +π soft pairs was discussed. The charm-quark which forms the D ⋆± is created in an initial hard-scattering interaction and is therefore sensitive to the evolution of the medium. However, the D ⋆± → D 0 + π soft decay occurs outside the medium and the daughter D 0 +soft-pion angular correlation is a result of vacuum decay. This contributes to the measured D 0 + h ± correlation mainly in the (∆η, ∆φ) = (0, 0) bin. These decay daughter-pairs are treated as background.
The number of such D ⋆± decays can be measured via a three-body invariant mass distribution constructed as M Kππ soft − M Kπ where the D ⋆± appears as a peak in the range [0.143,0.147] GeV/c 2 [40]. The D ⋆± yield and its background reference are normalized and used to correct the final correlations as described in Sec. II and in Eq. (1). The peak amplitudes of the final correlations in the (0,0) bin are reduced by approximately 0.037 ± 0.012(stat.) ± 0.004(syst.), 0.046 ± 0.002(stat.) ± 0.018(syst.), and 0.015 ± 0.001(stat.) ± 0.013(syst.) in the 50-80%, 20-50% and 0-20% centrality classes, respectively. The quoted statistical errors are the statistical uncertainties in the correlation amplitude for D 0 +hadron correlations coming from a D * decay. The systematic uncertainties are estimated from the variation in the deduced D * yields associated with 20% changes in the magnitudes of the combinatoric background in the invariant mass distributions used to identify D * decays. The total systematic effects from the D * ± contamination and its correction are further discussed in Sec. V.

IV. RESULTS
The per-pair normalized D 0 + h ± symmetrized correlations for centralities 50-80%, 20-50% and 0-20% with 2 < p T,D 0 < 10 GeV/c, defined in Eq. (1), are shown in the left-hand column of panels in Fig. 3. Significant structures are visible in the correlations which exceed the statistical fluctuations, including ∆η-independent near-side and away-side (AS) structures on ∆φ, and a near-side 2D peak which increases in width on ∆η with centrality. The structures are similar to the dominant features reported previously in unidentified and LF identified 2D dihadron angular correlations [15][16][17] from Au+Au collisions at 200 GeV. The broadening on ∆η of the NS 2D peak with increasing centrality, observed in these D 0 +h ± correlations, is similar to that reported for unidentified dihadron correlations [15].

A. Model Fitting
A quantitative representation of the D 0 + h ± correlations and centrality trends is facilitated by fitting the data with a model with a minimum number of elements which are chosen to describe the visible features in the data. The D 0 + h ± correlations are visually similar to previously reported, unidentified dihadron correlations; we therefore adopted the fitting model in Ref. [15]. Additional and/or alternate model elements were included in the study of systematic uncertainties, discussed in the next section. We assume a NS 2D Gaussian centered at (∆η, ∆φ) = (0, 0), an AS 2D Gaussian centered at (0, π), a ∆η-independent quadrupole, and an overall constant offset. Both 2D Gaussians are required to be periodic on ∆φ. The model is given by where near-side Gaussian terms at ∆φ = ±2π, etc. , and away-side Gaussians at ∆φ = −π, ±3π, etc., are not listed but are included in the model.
The eight fitting parameters A 0 , A Q , A NS , σ ∆η,NS , σ ∆φ,NS , A AS , σ ∆η,AS , and σ ∆φ,AS were, in general, allowed to freely vary to achieve the best description of the correlations based on minimum χ 2 . A few physically motivated restrictions were imposed, however. The quadrupole amplitude is equal to the product of the single-particle azimuthal anisotropy amplitudes v D 0 2 v h ± 2 , assuming factorization. Because both v D 0 2 > 0 and v h ± 2 > 0 in this collision system, the quadrupole amplitude was required to be non-negative [9,15,41]. For the 20-50% and 0-20% centrality classes, the AS Gaussian width on ∆φ increased sufficiently that the periodic Gaussian distribution reached the dipole limit 2 . For these two centrality classes the AS 2D Gaussian was replaced with A D cos(∆φ−π) exp(−∆η 2 /2σ 2 ∆η,AS ). For the cases with an undefined σ ∆η,AS , the fits were consistent with no AS dependence on ∆η, and the term was therefore dropped from the fit function. Statistical fluctuations exacerbated the appearance of false, local minima in the χ 2 space, resulting in false fitting solutions with unphysically narrow structures. The multi-dimensional χ 2 space was mapped and necessary search limits were imposed to avoid these false solutions.
The model fits and the residuals in terms of their statistical significance, |nσ resid | = |(data -model)/(statistical error)|, are shown in the middle and right-hand columns of panels in Fig. 3. The residuals are generally consistent with statistical errors except in the outermost ∆η bins which were omitted from the fitting procedure. The fitting model in Eq. 3 exhausts the statistically significant information in the measurements. The centrality dependences of A Q , σ ∆η,NS , and σ ∆φ,NS , determined within the ∆η acceptance from −2 to +2 units, are shown in Fig. 4. Statistical errors and systematic uncertainties, discussed in the next section, are shown by the vertical error bars and cross-bars, respectively. The optimum fit parameters and errors are listed in Table I. The azimuthal width of the AS 2D Gaussian increases with centrality, reaching the dipole limit in the mid-central and most-central bins, providing evidence of rescattering in the medium.   The efficiency and acceptance corrected average number of associated charged particles correlated with each D 0 trigger, the NS associated yield, is approximately where The details of the derivation and calculation of the associated yield from the correlations are found in Appendix B. The first term on the right hand side of Eq. 4 is the efficiency corrected, charged particle multiplicity in the centrality class, obtained by interpolating dN ch /dη from Table III in Ref. [15] for 200 GeV Au+Au collisions. The volume of the NS correlation peak, V NS−peak in Appendix B, is given by the integral on the right hand side where F NS−peak in the present analysis is assumed to be the NS 2D Gaussian in Eq. (3). Here, we are including all NS correlations, other than the quadrupole, in the NS yield per D 0 . The ∆η acceptance correction factor in Appendix B, not included above, is approximately one. The yield per D 0 -trigger in the NS 2D Gaussian peak gives the average number of hadrons correlated with each D 0 within the acceptance. This number, shown in Fig. 5 for the assumed model description, increases significantly with centrality as do the widths, especially the width on ∆η, shown in Fig. 4. The null hypothesis, which is that the NS correlations per D 0 -trigger are not affected by the increasing size and density of the medium, i.e. a linear superposition hypothesis [15], is strongly violated by these results.

V. SYSTEMATIC UNCERTAINTIES
We estimated systematic uncertainties in the correlation data, in the fitting model parameters, and yields per trigger. Systematic uncertainties in the 2D D 0 + h ± correlation data come from multiple sources in the analysis including: (1)  widths and positions, (3) efficiency correction method, (4) B-meson feed-down contribution to the D 0 + h ± correlations, (5) non-primary (secondary) particle contamination, (6) uncertainty in the correction for D ⋆± → D 0 + π sof t contamination, (7) particle identification of the D 0 daughters, (8) D 0 signal and background yield estimates, (9) event-mixing multiplicity and z-vertex subbin widths, (10) PV position and beam+beam luminos-ity effects on event-wise multiplicity determination, and (11) pileup from untriggered, out-of-time collision events in the TPC. Uncertainties in the correlations were estimated for each of these sources and most were found to be either negligible or indistinguishable from statistical noise. Sources resulting in non-negligible uncertainties are discussed below along with a few others. Systematic uncertainties (1)-(3) were estimated by varying each cut or correction individually and examining the bin-wise changes in the correlations. The vast majority of the effects were less than 2σ of the statistical fluctuations, and were indistinguishable from noise. Variations in the topological cuts for the D 0 -candidate daughter kaon and pion DCA to the PV produced larger changes in the correlations, up to 4σ in the most-central data, which were modeled and included as a systematic uncertainty. All other systematic uncertainties from sources (1)-(3) were negligible, and therefore not included in the final systematic quadrature sum.
Systematic uncertainties from B-meson feed-down to D 0 -mesons were studied extensively in Ref. [9] where it was estimated that about 4% of the D 0 sample are from feed-down. Contributions to the true, primary D 0 + h ± correlations could be as much as 4% in the total correlation amplitude, affecting the overall normalization of the correlations. Uncertainties arising from non-primary (secondary) particle contamination were estimated in Ref. [15] for unidentified charged-particle correlations. In the present analysis secondary particle contamination is suppressed for the D 0 candidates due to the PV resolution afforded by the HFT. The remaining contamination in the associated particle sample contributes about ±1.5% overall uncertainty in the correlation amplitudes.
Systematic uncertainty in the magnitude of the D ⋆ → D 0 π sof t contamination in the (∆η, ∆φ) = (0, 0) bin was estimated by adjusting the scale of the combinatoric background in the M Kππ soft − M Kπ invariant mass distribution by 20% based on background fluctuations, causing the extracted D * ± yield to vary (see Sec. III E). The reduction of pairs in the (∆η, ∆φ) = (0, 0) bin from the D * ± correction altered the correlations by changing the shape of the NS jet-like peak, causing a subsequent alteration of the model fit parameters. These variations are included as additional systematic uncertainties on the extracted fit parameters.
Misidentified D 0 decay daughter particles (K ↔ π) are broadly dispersed in the M Kπ distribution and have negligible contributions as previously reported [9]. Variations in the estimate of D 0 signal and background (factors S and B) from 10% to 20%, as a result of changes in the background subtraction (Fig. 2), had negligible effect on the final correlations. Requiring all tracks used in the analysis to include one hit in the IST, which resolves particles from separate beam bunch crossings [36], essentially eliminates pileup contamination. Event-mixing sub-bin widths were sufficiently narrow to eliminate artifacts, and the small variations in track reconstruction efficiency with PV position and luminosity negligibly af-fect event-wise multiplicity determination.
Typical systematic uncertainties in the measured correlations in Fig. 3 are between 3-10%. The above, nonnegligible uncertainties affect the model parameters as indicated in the discussion.
The largest source of systematic uncertainty in the model parameters and NS peak yield is due to the choice of fitting model. The nominal fit-model was introduced in Sec. IV. Systematic uncertainties were estimated by including a cos(3∆φ) (sextupole), or by replacing the ∆ηdependent part of the NS 2D Gaussian in Eq. (3) with either a Lorentzian function 3 (leptokurtic) or a raisedcosine function 4 (platykurtic). To be included in the systematic uncertainty estimates, alternate fit models were required to have similar χ 2 and residuals as the nominal fits, with one unique χ 2 -minimum corresponding to physically reasonable parameters. Finally, the nominal fitting model was applied to a rebinned version of the correlations assuming 11 ∆η bins and 16 ∆φ bins, which resulted in small changes in the fit parameters.
Each of the above positive and negative systematic uncertainties in the correlation quantities resulting from the non-negligible sources of systematic uncertainty discussed in this section were added in quadrature, where positive and negative errors were combined separately. The nominal fitting model results, statistical errors, and the combined systematic uncertainties are listed in Table I.

VI. DISCUSSION
It is interesting to compare these results with expectations based on HF production models and on other correlation results. Angular correlations for D 0 +h ± predicted by perturbative QCD with conventional fragmentation, as in pythia [27,42,43] [27,43]. The decay-daughters from K 0 S and Λ were excluded from the associated track sample. 6 The modified NS 2D Gaussian used to fit the pythia correlations is is approximately twice the measured yield, but is within ∼ 2σ of the measurement. The NS yield per trigger and the NS ∆η width increase in the 20-50% and 0-20% centralities, relative to the perturbative QCD predictions, is similar to that reported in Ref. [15] for unidentified dihadron correlations. The present centrality trends for the assumed fitting model are consistent with the onset of significant increases in the NS peak amplitude and ∆η width at approximately 40-50% centrality [15], and are consistent with the appearance of a near-side ridge in these D 0 + hadron correlations. In addition, the increase in per-trigger yield and near-side peak widths in the 20-50% and 0-20% centrality classes relative to the pythia predictions, coincides with the suppression observed in the nuclear modification factor, R AA , for D 0 yields over the same 2-10 GeV/c p T range [6]. Both the present and the previous R AA observations imply increased D 0 + medium interactions in more central collisions.
The reduction of the D 0 yield in more central collisions in the present p T range (R AA [6,7]) could imply that the relative yield of D 0 triggers emitted from the interior of the collision medium is suppressed. This resulting surface bias in the observed D 0 sample makes the observed increases in amplitude and width of the NS D 0 +h ± correlations, in more central collisions, even more remarkable.
We also compare the D 0 + h ± correlation structures with similar, unidentified charged-particle correlations for 200 GeV Au+Au collisions [44] in which p T of the "trigger" particle is binned while the associated particle's p T ≥ 0.15 GeV/c. In Ref. [44] the minimumbias 200 GeV Au+Au data were divided into centrality classes 0-5%, 5-10%, 10-20%, · · · , 60-70%, plus several trigger p T bins, using the same |η| ≤ 1 and full 2π azimuth acceptance as in the present analysis. A similar fitting model to that used here, consisting of an offset, dipole, quadrupole and NS 2D Gaussian, was assumed in Ref. [44]. The D 0 + h ± and unidentified dihadron correlations are compared using a common trigger-particle p T as was assumed in Ref. [9] for the analysis of D 0 v 2 . The light-and heavy flavor results could also be compared via a common trigger-particle velocity [45] assuming that diffusion in a dispersive medium is the dominant process. The highest three trigger p T bins in Ref. [44], [2.1,3.1] GeV/c, [3.1,4.7] GeV/c and [4.7,7.0] GeV/c, offer the best overlap with the D 0 p T range, [2,10] GeV/c, used here. Application of Eq. (4) to the dihadron analysis gives the approximate number of correlated charged particles per trigger-particle in the NS 2D peak as (dN ch /2πdη)V NS−peak , whereN ch is the efficiency corrected, average number of charged-particles with p T ≥ 0.15 GeV/c, and V NS−peak is the volume integral of the NS 2D Gaussian within the acceptance.
Centrality fraction weighted results from Ref. [44] for the 50-70%, 20-50% and 0-20% centralities are shown in Figs. 4 and 5 for the per-trigger NS 2D peak yields and for the Gaussian widths on ∆φ and ∆η for trigger p T bins [2.1,3.1] GeV/c ( p T = 2.56 GeV/c) and [4.7,7.0] GeV/c ( p T = 5.7 GeV/c) by the upside down triangles and dot symbols, respectively. The centrality trend of the pertrigger yields for the lower-p T LF dihadron correlations follows the D 0 +h ± trend fairly well; the higher p T results do not increase as rapidly with centrality. The dihadron azimuthal widths for p T = 2.56 GeV/c are similar in magnitude to the D 0 + h ± widths, and follow a similar trend with centrality. The σ ∆η,NS width for p T = 2.56 GeV/c increases by a factor of 4 in the most-central bin relative to the most-peripheral bin, while the D 0 +h ± results increase from peripheral to mid-centrality, but are constant within errors from mid-central to most-central.
The dihadron and D 0 + h ± quadrupole amplitudes are compared in Fig. 4 where the 20-50% and 0-20% centrality results are consistent within errors. However, the D 0 + h ± quadrupole amplitude is smaller than the dihadron amplitude for peripheral collisions.
Finally, we compare the D 0 -meson azimuthal anisotropy parameter v D 0 2 inferred from the present analysis with a previous STAR Collaboration measurement [9] which used both event plane and two-particle correlation methods. The quadrupole amplitude is (see Table I compares well with the previous measurement which included the 10-40% centrality. In conventional v 2 analyses [9, 10] η-gaps are used to reduce the contribution of "non-flow" corre-lations, e.g. jet fragmentation. The extended near-side correlation structure on ∆η in the two more central bins in Fig. 3 indicates that either large η-gaps or higher-order cumulant methods [41] are required for conventional v 2 analyses in the more central collisions.

VII. SUMMARY AND CONCLUSIONS
In this paper we report the first measurement of two-dimensional angular correlations between D 0 -mesons and unidentified charged hadrons produced in relativistic heavy-ion collisions. Attention was focused on the centrality evolution of the near-side, 2D correlation peak widths and associated yields. Results for the associated hadron yield per D 0 -meson trigger and the 2D widths of the correlated angular distribution, obtained by fitting the data with a 2D Gaussian model, were used to characterize the centrality dependence. We find that the associated hadron per D 0 -meson yield and 2D widths increase significantly from peripheral to central collisions. With the D 0 -meson serving as a proxy for a charm-quark jet, this measurement is a first attempt to understand heavy flavor jets in heavy-ion collisions at RHIC energies.
The increase in the near-side correlation yield and width coincides in both the centrality and D 0 -meson p T ranges where the D 0 -meson nuclear modification factor R AA is suppressed. Both results are consistent with the expectation that the interactions between the charm quark and the medium increase with centrality. The present results complement previous studies of D 0 -meson spectra, R AA and v 2 . The centrality trends and magnitudes of the NS 2D Gaussian fit parameters, qualitatively agree with a similar analysis of dihadron 2D correlations for 200 GeV Au+Au minimum-bias collisions for similar p T and centrality ranges. These results imply that the effective strength and centrality dependence of heavy flavor particle interactions with the medium are similar to that observed for light flavor particles, as seen in previous, complementary studies.
In conclusion, the near-side, non-quadrupole correlated hadrons, which are associated with D 0 -mesons, display a large increase in per-trigger yield and 2D widths, especially the width along relative pseudorapidity, for collisions more central than about 50%. This ridge formation phenomenon has been observed in light flavor dihadron correlations at both the RHIC and the LHC and is now observed in D 0 -meson + hadron correlations in Au+Au collisions at 200 GeV. ∆ρ D 0 +h αρ ME,D 0 +h = ∆ρ sig − ∆ρ D 0 +π soft − ∆ρ BG+h αρ ME,D 0 +h .
Substituting these results into Eq. (7) and using the average of the left and right side-bands to estimate the background correlations, gives the final expression in Eq. (1) of the main text: ∆ρ D 0 +h αρ ME,D 0 +h = S + B S ∆ρ sig − ∆ρ D 0 +π soft αρ ME,sig − B S ∆ρ SB α SB ρ ME,SB .
Appendix B: Calculation of the near-side correlated yield per D 0 trigger.
The correlated pair yield per D 0 trigger in the NS 2D peaked correlation structure, Y NS−peak /N D 0 , is estimated by summing that portion of the correlation fitting model in Eq. (3) over the (∆η, ∆φ) acceptance, including efficiency and acceptance corrections, and dividing by the efficiency corrected number of D 0 mesons, N D 0 , used in the analysis. This estimate is given by [12] Y NS−peak /N D 0 = 1 N D 0 × ∆η,∆φ   ∆n D 0 +h αρ ME,D 0 +h αρ max ME,D 0 +h where ∆n D 0 +h = δ ∆η δ ∆φ ∆ρ D 0 +h in each (∆η, ∆φ) bin, and δ ∆η , δ ∆φ are the bin widths on ∆η and ∆φ. Also in Eq. (11), αρ max ME,D 0 +h is the maximum value of the normalized, mixed-event pair distribution, evaluated by averaging over the ∆φ bins for ∆η = 0. The ratio in the denominator represents the detector acceptance distribution normalized to 1.0 at the maximum. Rearranging Eq. (11) gives Y NS−peak /N D 0 = αρ max ME,D 0 +h N D 0 × ∆η,∆φ δ ∆η δ ∆φ ∆ρ D 0 +h αρ ME,D 0 +h NS−peak = αρ max ME,D 0 +h where the summation in the second line of Eq. (12) is defined as V NS−peak , the volume of the NS peak correlation structure. The ratio on the RHS of the third line of Eq. (12) can be estimated from the measured numbers of D 0 and D 0 + h ± ME pairs, provided both numerator and denominator are corrected for inefficiencies. A simpler form is given in the following in which the required efficiency corrected quantities are more readily obtained.
The maximum value of the efficiency corrected, normalized mixed-event density equals the fraction of the total number of D 0 + h ± pairs in a ∆η = 0, ∆φ bin per bin area. This is given by where ε is the number of events in the centrality class, n D 0 andn h are the efficiency corrected, event-averaged number of D 0 mesons and associated h ± particles in the acceptance, N ∆η and N ∆φ are the numbers of ∆η and ∆φ bins, where N ∆η is odd and N ∆φ is a multiple of four. The second ratio on the RHS of Eq. (13) is the fraction of D 0 + h ± pairs in an average ∆η = 0, ∆φ bin. The efficiency corrected number of D 0 mesons is N D 0 = εn D 0 . The ratio in Eq. (12) simplifies to αρ max ME,D 0 +h (14) where N ∆η N ∆φ δ ∆η δ ∆φ = 4πΩ η , and Ω η is the single particle pseudorapidity acceptance which equals 2 units for the STAR TPC [34]. In the last step we assumed that the number of K, π daughters is much less than the event multiplicity, such thatn h is well approximated by event multiplicity N ch . The final NS-peak correlated yield per D 0 trigger is given by (15) where dN ch /2πdη is efficiency corrected [15].