Centrality and transverse momentum dependence of $D^0$-meson production at mid-rapidity in Au+Au collisions at ${\sqrt{s_{\rm NN}} = \rm{200\,GeV}}$

TL;DR

This paper presents detailed measurements of $D^0$-meson production in Au+Au collisions at 200 GeV, revealing insights into charm quark energy loss, collective flow, and freeze-out properties in the quark-gluon plasma.

Contribution

It provides new high-precision $D^0$-meson spectra and flow measurements, using the Heavy Flavor Tracker, to study charm quark interactions and decoupling in heavy-ion collisions.

Findings

Charm quarks experience significant energy loss in the medium.

$D^0$ mesons exhibit smaller radial flow velocity than light hadrons.

Evidence of charm quarks gaining collective motion during medium evolution.

Abstract

We report a new measurement of $D^0$-meson production at mid-rapidity ($|y|$\,$<$\,1) in Au+Au collisions at ${\sqrt{s_{\rm NN}} = \rm{200\,GeV}}$ utilizing the Heavy Flavor Tracker, a high resolution silicon detector at the STAR experiment. Invariant yields of $D^0$-mesons with transverse momentum $p_{T}$ $\lesssim 9$\,GeV/$c$ are reported in various centrality bins (0--10\%, 10--20\%, 20--40\%, 40--60\% and 60--80\%). Blast-Wave thermal models are used to fit the $D^0$-meson $p_{T}$ spectra to study $D^0$ hadron kinetic freeze-out properties. The average radial flow velocity extracted from the fit is considerably smaller than that of light hadrons ($\pi,K$ and $p$), but comparable to that of hadrons containing multiple strange quarks ($\phi,\Xi^-$), indicating that $D^0$ mesons kinetically decouple from the system earlier than light hadrons. The calculated $D^0$ nuclear modification factors re-affirm that charm quarks suffer large amount of energy loss in the medium, similar to those of light quarks for $p_{T}$\,$>$\,4\,GeV/$c$ in central 0--10\% Au+Au collisions. At low $p_{T}$, the nuclear modification factors show a characteristic structure qualitatively consistent with the expectation from model predictions that charm quarks gain sizable collective motion during the medium evolution. The improved measurements are expected to offer new constraints to model calculations and help gain further insights into the hot and dense medium created in these collisions.

Tables

Table 1. Table 1: Estimated values of average number of binary collisions ( N bin subscript 𝑁 bin N_{\rm bin} ), average number of participants ( N part subscript 𝑁 part N_{\rm part} ) and trigger correction factors ( ε trg subscript 𝜀 trg \varepsilon_{\rm trg} , uncertainties negligible) for various centrality bins obtained from the MC Glauber model fit to the measured multiplicity distributions.

Centrality	$⟨N_{\mathrm{bin}}⟩$				$⟨N_{\mathrm{part}}⟩$			${\epsilon }_{\mathrm{trg}}$
0–10 %	938.8	$\pm$	26.3		319.4	$\pm$	3.4	1.0
10–20 %	579.9	$\pm$	28.8		227.6	$\pm$	7.9	1.0
20–40 %	288.3	$\pm$	30.4		137.6	$\pm$	10.4	1.0
40–60 %	91.3	$\pm$	21.0		60.5	$\pm$	10.1	0.92
60–80 %	21.3	$\pm$	8.9		20.4	$\pm$	6.6	0.65

Table 2. Table 2: Several key characteristic parameters for each subsystem of the STAR HFT detector.

Subsystem	Radius (cm)	Length (cm)	Thickness at $\eta$ $=$ 0 ($X_0$)	Pitch Size ($\mu {\text{m}}^2$)
PXL inner layer	2.8	20	0.52% (0.39%†)	20.7$\times$20.7
PXL outer layer	8.0	20	0.52%	20.7$\times$20.7
IST	14.0	50	1.0%	600$\times$6000
SSD††	22.0	106	1.0%	95$\times$40000

Table 3. Table 3: TPC and TOF selection cuts for K 𝐾 K and π 𝜋 \pi tracks.

Variable		$K^{\mp }$	${\pi }^{\pm }$
$p_T$ (GeV/$c$)	$>$	0.6	0.6
$\left|\eta \right|$		1.0	1.0
nHitsFit (TPC)	$>$	20	20
$\left|n{\sigma }_X\right|$		2.0	3.0
$\left|\Delta 1/\beta \right|$ (if TOF matched)		0.03	0.03

Table 4. Table 4: Topological cuts used for D 0 superscript 𝐷 0 D^{0} reconstruction in 0–10% most central collisions for separate p T subscript 𝑝 𝑇 p_{T} intervals.

$0-10\%|p_T$ (GeV/$c$)		(0,0.5)	(0.5,1)	(1,2)	(2,3)	(3,5)	(5,8)	(8,10)
Decay Length ($\mu$m)	$>$	100	199	227	232	236	255	255
${\mathrm{DCA}}_{12}$ ($\mu$m)		71	64	70	63	82	80	80
${\mathrm{DCA}}_{D^0}$ ($\mu$m)		62	55	40	40	40	44	44
${\mathrm{DCA}}_{\pi }$ ($\mu$m)	$>$	133	105	93	97	67	55	55
${\mathrm{DCA}}_K$ ($\mu$m)	$>$	138	109	82	94	76	54	54
$\cos \left(\theta \right)$	$>$	0.95	0.95	0.95	0.95	0.95	0.95	0.95

Table 5. Table 5: Summary of systematic uncertainties, in percentage, on the D 0 superscript 𝐷 0 D^{0} invariant yield in 0–10% and 60–80% collisions and on R CP subscript 𝑅 CP R_{\rm CP} (0–10%/60–80%).

Source	Systematic uncertainty [%]			Correlation in $p_T$
	0–10%	60–80%	$R_{\mathrm{CP}}$(0–10%/60–80%)
Signal extra.	1-6	1-12	2-13	uncorr.
Double mis-PID	1-7	1-7	negligible	uncorr.
${\epsilon }_{\mathrm{TPC}}$	5-7	5-8	3-7	largely corr.
${\epsilon }_{\mathrm{HFT}}$	3-15	3-20	3-20	largely corr.
${\epsilon }_{\mathrm{PID}}$	3	3	negligible	largely corr.
${\epsilon }_{\mathrm{vtx}}$	5	9-18	10-18	largely corr.
BR.	0.5		0	global
$N_{\mathrm{bin}}$	2.8	42	42	global

Table 6. Table 6: ⟨ β ⟩ delimited-⟨⟩ 𝛽 \langle\beta\rangle and ( q − 1 ) 𝑞 1 (q-1) from the Tsallis Blast-Wave fits to the D 0 superscript 𝐷 0 D^{0} data in different centralities.

Centrality	$⟨\beta ⟩$($c$)				$q-1$
0–10 %	0.263	$\pm$	0.018		0.066	$\pm$	0.008
10–20 %	0.255	$\pm$	0.022		0.068	$\pm$	0.010
20–40 %	0.264	$\pm$	0.015		0.070	$\pm$	0.007
40–60 %	0.251	$\pm$	0.023		0.074	$\pm$	0.011
60–80 %	0.217	$\pm$	0.037		0.075	$\pm$	0.010

Table 7. Table 7: D ¯ 0 superscript ¯ 𝐷 0 \overline{D}^{0} / D 0 superscript 𝐷 0 D^{0} ratio for various centrality bins obtained from the fit to data in Fig.  34 .

Centrality	${\overline{D}}^0$/$D^0$
0–10 %	1.104	$\pm$	0.021
10–20 %	1.071	$\pm$	0.019
20–40 %	1.060	$\pm$	0.015
40–60 %	1.073	$\pm$	0.022
60–80 %	0.943	$\pm$	0.039

Equations

$$R_{\mathrm{AA}}(p_T)=\frac{1}{⟨T_{\mathrm{AA}}⟩}\frac{d{N}_{\mathrm{AA}}/d{p}_T}{d{\sigma }_{pp}/d{p}_T}.$$

$$n{\sigma }_X=\frac{1}{R}\ln \frac{⟨dE/dx{⟩}_{mea.}}{⟨dE/dx{⟩}_X},$$

$$\frac{d^{2}N}{2\pi p_{T}d{p}_{T}dy}=\frac{1}{B.R.}\times \frac{N^{\mathrm{raw}}}{N_{\mathrm{evt}}2\pi p_T\Delta p_T\Delta y}$$

$$\times \frac{1}{{\epsilon }_{\mathrm{trg}}\times {\epsilon }_{\mathrm{TPC}}\times {\epsilon }_{\mathrm{HFT}}\times {\epsilon }_{\mathrm{PID}}\times {\epsilon }_{\mathrm{vtx}}},$$

$$R_{\mathrm{CP}}=\frac{d^{2}N/d{p}_{T}dy}{N_{\mathrm{bin}}}{|}_{\mathrm{cen}}\times \frac{N_{\mathrm{bin}}}{d^{2}N/d{p}_{T}dy}{|}_{\mathrm{peri}}.$$

$$\frac{d^{2}N}{2\pi p_{T}d{p}_{T}dy}=$$

$$\displaystyle\times\bigg{(}1+\frac{\sqrt{p_{T}^{2}+m_{0}^{2}}-m_{0}}{nT}\bigg{)}^{-n},$$

$$\frac{d^{2}N}{2\pi m_{T}d{m}_{T}dy}=\frac{dN/dy}{2\pi T_{\mathrm{eff}}(m_0+T_{\mathrm{eff}})}{e}^{-(m_T-m_0)/T_{\mathrm{eff}}}.$$

$$\displaystyle\frac{d^{2}N}{2\pi p_{T}dp_{T}dy}=\frac{dN}{dy}\frac{4(n-1)(n-2)}{2\pi(n-3)^{2}\langle p_{T}\rangle^{2}}\bigg{(}1+\frac{2p_{T}}{\langle p_{T}\rangle(n-3)}\bigg{)}^{-n},$$

$$T_{\mathrm{eff}}=T_{\mathrm{fo}}+m_0⟨{\beta }_T{⟩}^2.$$

$$\frac{dN}{p_{T}d{p}_T}=\frac{dN}{m_{T}d{m}_T}\propto$$

$$\displaystyle\int_{0}^{R}rdrm_{T}I_{0}\bigg{(}\frac{p_{T}\sinh\rho}{T_{\rm kin}}\bigg{)}K_{1}\bigg{(}\frac{m_{T}\cosh\rho}{T_{\rm kin}}\bigg{)},$$

$$\beta ={\beta }_s{\left(\frac{r}{R}\right)}^n,$$

$$\frac{dN}{m_{T}d{m}_T}\propto m_T{\int }_{-Y}^{+Y}\cosh (y)dy{\int }_{-\pi }^{+\pi }d\varphi {\int }_0^{R}rdr$$

$$\displaystyle\bigg{(}1+\frac{q-1}{T_{\rm kin}}\big{(}m_{T}\cosh(y)\cosh(\rho)-p_{T}\sinh(\rho)\cos(\phi)\big{)}\bigg{)}^{\textstyle-\frac{1}{q-1}}$$

Full text

††preprint: APS/123-QED

I Introduction

The heavy ion program at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) focuses on the study of strong interactions and Quantum Chromodynamics (QCD) at high temperature and density. Over the last few decades, experimental results from RHIC and LHC using light flavor probes have demonstrated that a strongly-coupled Quark-Gluon Plasma (sQGP) is created in these heavy-ion collisions. The most significant evidence comes from the strong collective flow and the large high transverse momentum ($p_{T}$) suppression in central collisions for various observed hadrons including multi-strange-quark hadrons $\phi$ and $\Omega$ Adams et al. (2005); Adcox et al. (2005); Muller et al. (2012); Adamczyk et al. (2016); Abelev et al. (2015).

Heavy quarks ($c$,$b$) are created predominantly through initial hard scatterings due to their large masses Lin and Gyulassy (1995); Cacciari et al. (2005). The modification to their production in transverse momentum due to energy loss and radial flow and in azimuth due to anisotropic flows is sensitive to heavy quark dynamics in the partonic sQGP phase Moore and Teaney (2005). Recent measurements of high-$p_{T}$ $D$-meson production at RHIC and LHC show a strong suppression in the central heavy-ion collisions Abelev et al. (2012); Adam et al. (2016); Sirunyan et al. (2018a); Adamczyk et al. (2014a). The suppression is often characterized by the nuclear modification factor $R_{\rm AA}$, defined as

[TABLE]

where $dN_{\rm AA}/dp_{T}$ and $d\sigma_{pp}/dp_{T}$ are particle production yield and cross section in A+A and $p$+$p$ collisions, respectively. The nuclear thickness function $T_{\rm AA}=\langle N_{\rm bin}\rangle/\sigma_{pp}^{\rm inel}$ is often calculated using a Monte-Carlo Glauber model, where $\langle N_{\rm bin}\rangle$ is the average number of binary collisions and $\sigma_{pp}^{\rm inel}$ is the total inelastic $p$+$p$ cross section. The $D$-meson $R_{\rm AA}$ is similar to that of light hadrons for p_{T}$$\,> 4 GeV/$c$, suggesting significant energy loss for charm quarks inside the sQGP medium. The measured $D$-meson anisotropic flow shows that $D$-mesons also exhibit significant elliptic and triangular flow at RHIC and LHC Abelev et al. (2013a, 2014); Sirunyan et al. (2018b); Adamczyk et al. (2017a). The flow magnitude when scaled with the transverse kinetic energy is similar to that of light and strange flavor hadrons. This indicates that charm quarks may have reached thermal equilibrium in these collisions at RHIC and LHC.

In this article, we report measurements of the centrality dependence of $D^{0}$-meson transverse momentum spectra at mid-rapidity ($|y|$ $<$ 1) in Au+Au collisions at $\sqrt{s_{{}_{\rm NN}}}=200$ GeV. The measurements are conducted at the Solenoidal Tracker At RHIC (STAR) experiment utilizing the high resolution silicon detector (the Heavy Flavor Tracker, HFT) Contin et al. (2018). The paper is organized in the following order: In Sec. II, we describe the detector setup and dataset used in this analysis. In Sec. III, we present the topological reconstruction of $D^{0}$ mesons in the Au+Au collision data, followed by Sec. IV and Sec. V for details on efficiency corrections and systematic uncertainties. We present our measurement results and physics discussions in Sec. VI. Finally, we end the paper with a summary in Sec. VII .

II Experimental setup and Dataset

The dataset used in this analysis consists of Au+Au collision events at ${\sqrt{s_{\rm NN}}=\rm{200\,GeV}}$ collected in the 2014 year run. The main detectors used in this analysis are the Time Projection Chamber (TPC), the HFT, the Time of Flight (TOF) detector and the Vertex Position Detector (VPD).

II.1 Tracking and Particle Identification Subsystems

Precision tracking for this analysis is achieved with the TPC and HFT detectors. Particle identification for stable hadrons is performed with a combination of the ionization energy loss ($dE/dx$) measurement with the TPC and the time-of-flight ($tof$) measurement with the TOF detector. The event start time is provided by the VPD. Both the TPC and TOF detectors have full azimuthal coverage with a pseudo-rapidity range of $|\eta|$ $<$ 1 Anderson et al. (2003); Llope (2012). The TPC and TOF subsystems have been extensively used in many prior STAR analyses, including $D$-meson measurements  Adamczyk et al. (2012, 2014a, 2016). The HFT detector provides measured space points with high precision that are used to extend track trajectories and offer high-pointing resolution to the vicinity of the event vertex.

II.2 Trigger and Dataset

The minimum bias trigger used in this analysis is defined as a coincidence between the east and west VPD detectors located at 4.4 $<$ $|\eta|$ $<$ 4.9 Llope et al. (2004). Each VPD detector is an assembly of nineteen small detectors, each consisting of a Pb converter followed by a fast plastic scintillator read out by a photomultiplier tube. To efficiently sample the collision events in the center of the HFT acceptance, an online cut on the collision vertex position along the beam line (calculated via the time difference between the east and west VPD detectors) $|V_{z}^{\rm VPD}|$ $<$ 6 cm is applied. The decrease in the coincidence probability in the VPD degrades the online VPD vertex resolution in peripheral low multiplicity events. These inefficiencies are corrected in the offline analysis with a method discussed in the next section.

Events used in this analysis are selected with the offline reconstructed collision vertex within 6 cm of the TPC and HFT centers along the beam direction to ensure uniform and large acceptance. The maximum total drift time of ionization electrons inside the TPC is about 40 $\upmu$s while the hadronic Au+Au collision rate is typically around 40 kHz for this dataset. There is a finite chance that more than one event is recorded in the TPC readout event frame. The VPD is a fast detector which can separate events from different bunch crossings (one bunch crossing interval at RHIC is 106 ns). In order to suppress the chance of selecting a wrong vertex from collisions happening in bunch crossings different from that of the trigger, the difference between the event vertex $z$ coordinate $V_{z}^{\rm TPC}$ and the $V_{z}^{\rm VPD}$ is required to be less than 3 cm. Approximately 9$\times 10^{8}$ minimum bias triggered events with 0–80% centrality pass the selection criteria and are used in this analysis.

II.3 Centrality Selection and Trigger Inefficiency

The centrality is selected using the measured charged global track multiplicity $N_{\rm ch}^{\rm raw}$ at mid-rapidity within $|\eta|\ \rm{<0.5}$ and corrected for the online VPD triggering inefficiency using a Monte Carlo (MC) Glauber simulation. 0–X% centrality is defined as the 0–X% most central in terms of the total hadronic cross section determined by the impact parameter between two colliding nuclei. In this analysis, the dependence of $N_{\rm ch}^{\rm raw}$ on the collision vertex position and the beam luminosity has been taken into account. The measured track multiplicity distribution from Au+Au 200 GeV from RHIC run 2014, corrected for the vertex and luminosity dependence, is shown in Fig. 1. The measured distribution is fit by the MC Glauber calculation in the high multiplicity region. One can observe that the fitted MC Glauber calculation matches the real data well for $N_{\rm ch}^{\rm raw}$ $>$ 100, while the discrepancy in the low multiplicity region shows the VPD trigger inefficiency. Figure 1 panel (b) shows the ratio between MC and data. Centrality is defined according to the MC Glauber model distribution shown in panel (a). Events in the low-multiplicity region are weighted with the ratio shown in panel (b) in all the following analysis as a correction for the inefficiency in trigger.

Table 1 lists the extracted values of average number of binary collisions ($N_{\rm bin}$), average number of participants ($N_{\rm part}$) and trigger inefficiency correction factors ($\varepsilon_{\rm trg}$) and their uncertainties in various centrality bins. The $\varepsilon_{\rm trg}$ correction factor is applied event-by-event in the analysis when combining centrality bins.

II.4 Heavy Flavor Tracker

The HFT Contin et al. (2018) is a high resolution silicon detector system that aims for topological reconstruction of secondary vertices from heavy flavor hadron decays. It consists of three silicon subsystems: the Silicon Strip Detector (SSD), the Intermediate Silicon Tracker (IST), and two layers of the PiXeL (PXL) detector. Table 2 lists the key characteristic parameters of each subsystem. The SSD detector was still under commissioning when the dataset was recorded, and therefore is not used in the offline data production and this analysis. The PXL detector uses the new Monolithic Active Pixel Sensors (MAPS) technology Contin et al. (2018). This is the first application of this technology in a collider experiment. It is specifically designed to measure heavy flavor hadron decays in the high multiplicity heavy-ion collision environment.

In the offline reconstruction, tracks are reconstructed in the TPC first and then extended to the HFT detector to find the best fit to the measured high resolution spatial points. A Kalman filter algorithm that considers various detector material effects is used in the track extension Kalman (1960). Considering the level of background hits in the PXL detector due to pileup hadronic and electromagnetic collisions, tracks are required to have at least one hit in each layer of the PXL and IST sub-detectors. Figure 2 shows the track pointing resolution to the primary vertex in the transverse plane ($\sigma_{\rm XY}$) in panel (a) and along the longitudinal direction ($\sigma_{\rm Z}$) in panel (b) as a function of total momentum ($p$) for identified particles in 0–80% centrality Au+Au collisions. The design goal for the HFT detector was to have a pointing resolution better than 55 $\upmu$m for 750 MeV charged kaon particles. Figure 2 demonstrates that the HFT detector system meets the design requirements. This performance enables precision measurement of $D$-meson production in high multiplicity heavy-ion collisions.

III $D^{0}$-meson reconstruction

$D^{0}$ and $\overline{D}^{0}$ mesons are reconstructed via the hadronic decay channel $D^{0}\rightarrow K^{-}+\pi^{+}$ and its charge conjugate channel with a branching ratio ($B.R.$) of 3.89%. In what follows, we imply $(D^{0}+\overline{D}^{0})/2$ when using the term $D^{0}$ unless otherwise specified. $D^{0}$ mesons decay with a proper decay length of $c\tau\sim 123\ \upmu$m after they are produced in Au+Au collisions. We utilize the high-pointing resolution capability enabled by the HFT detector to topologically reconstruct the $D^{0}$ decay vertices that are separated from the collision vertices, which drastically reduces the combinatorial background ($\sim$five orders of magnitude) and improves the measurement precision.

Charged pion and kaon tracks are reconstructed with the TPC and HFT. Tracks are required to have at least 20 measured TPC points out of maximum 45 to ensure a good momentum resolution. To enable high pointing precision, both daughter tracks are required to have at least one measured hit in each layer of the PXL and IST as described above. Particle identification is achieved via a combination of the ionization energy loss measurement in the TPC and the time-of-flight measurement in the TOF. The resolution-normalized $dE/dx$ deviation from the expected values is defined as:

[TABLE]

where $\langle{dE/dx}\rangle_{\rm{mea.}}$ and $\langle{dE/dx}\rangle_{X}$ represent measured and expected values with a hypothesis of particle $X$, and $R$ is the $dE/dx$ resolution (typically $\sim$ 8% Anderson et al. (2003)). The $n\sigma_{X}$ distribution should be close to a standard Gaussian for each corresponding particle species (mean $=$ 0, $\sigma=$ 1) with good $dE/dx$ calibration. Pion (kaon) candidates are selected by a requirement of the measured $dE/dx$ to be within three (two) standard deviations ($|n\sigma_{X}|$) from the expected value. When tracks have matched hits in the TOF detector, an additional requirement on the measured inverse particle velocity ($1/\beta$) to be within three standard deviations from the expected value ($|\Delta 1/\beta|$) is applied for either daughter track. Figures 3 and  4 show examples of the particle identification capability from the TPC and TOF. Tracks within the kinematic acceptance $p_{T}$ $>$ 0.6 GeV/$c$ and $|\eta|$ $<$ 1 are used to combine and make pairs. The choice of the $p_{T}$ $>$ 0.6 GeV/$c$ cut is an optimized consideration to balance the loss of signal acceptance when tightening the cut, and the increase in background due to the HFT fake matches when loosening this cut (see Sec. IV.2). The threshold has been varied for systematic uncertainty evaluation. See Sec. V for details. Table 3 lists the TPC and TOF selection cuts for daughter kaon and pion tracks used for $D^{0}$ reconstruction.

With a pair of two daughter tracks, pion and kaon, the $D^{0}$ decay vertex is reconstructed as the middle point on the distance of the closest approach between the two daughter trajectories. One of the dominant background sources is the random combination of $K\pi$ pairs directly from the collision point. With the selection of the following topological variables, the background level can be greatly reduced.

• •

Decay Length: the distance between the reconstructed decay vertex and the Primary Vertex (PV).

• •

Distance of Closest Approach (DCA) between the two daughter tracks ($\rm{DCA_{12}}$).

• •

DCA between the reconstructed $D^{0}$ and the PV ($\rm{DCA_{D^{0}}}$).

• •

DCA between the pion and the PV ($\rm{DCA_{\pi}}$).

• •

DCA between the kaon and the PV ($\rm{DCA_{K}}$).

• •

Angle between the $D^{0}$ momentum and the direction of the decay vertex with respect to the PV ($\theta$).

The schematic in Fig. 5 shows the topological variables used in the analysis, where $\vec{P}$ represents the $D^{0}$ momentum. The Decay Length and angle $\theta$ follow the formula: $\rm DCA_{D^{0}}$ = Decay Length $\times$ sin($\theta$). The cuts on the topological variables for this analysis are optimized using a Toolkit for Multivariate Data Analysis (TMVA) package integrated in the ROOT framework in order to obtain the greatest signal significance Voss et al. (2007). The Rectangular Cut optimization method from the TMVA package is chosen in this analysis, similar as in our previous publication Adamczyk et al. (2017a). The optimization is conducted for different $D^{0}$ $p_{T}$ bins and different centrality bins. Table 4 lists a set of topological cuts for 0–10% central Au+Au collisions.

Figure 6 shows the invariant mass distributions of $K\pi$ pairs in the $p_{T}$ region of 0–10 GeV/$c$ for 0–80% minimum bias and the 0–10% most central collisions, and 0–8 GeV/$c$ for 60–80% peripheral collisions, respectively. The reason of choosing a different $p_{T}$ range for the 60–80% centrality bin is because no signal is observed beyond the current statistics. The combinatorial background is estimated with the same-event (SE) like-sign (LS) pairs (blue histograms) and the mixed-event (ME) unlike-sign (US) (grey histograms) technique in which $K$ and $\pi$ from different events of similar characteristics ($V_{Z}$, centrality, event plane angle) are paired. The mixed-event spectra are normalized to the like-sign distributions in the mass range of 1.7–2.1 GeV/$c^{2}$. After the subtraction of the mixed-event unlike-sign combinatorial background from the same-event unlike-sign pairs (black open circles), the remainder distributions are shown as red solid circles in each panel. Compared to the previous $D^{0}$ measurement Adamczyk et al. (2014a), the $D^{0}$ signal significance is largely improved by a factor of $\sim$ 15 using the same amount of event statistics.

Figures 7 and  8 show the invariant mass distributions in the same centrality bins as Fig. 6 but for different $p_{T}$ ranges: 0 $<$ $p_{T}$ $<$ 0.5 GeV/$c$ in Fig. 7 and 6 $<$ $p_{T}$ $<$ 8 GeV/$c$ in Fig. 8.

After the combinatorial background is subtracted, the residual $K\pi$ invariant mass distributions are then fit to a Gaussian plus linear function. The linear function is used to represent remaining correlated background from either partial reconstruction of charm mesons or jet fragments. The $D^{0}$ raw yields are extracted from the Gaussian function fit results while different choices of fit ranges, background functional forms, histogram counting vs. fitting methods etc. have been used to estimate systematic uncertainties on the raw yield extraction. See Sec. V for details.

IV Efficiencies and Corrections

The reconstructed $D^{0}$ raw yields are calculated in each centrality, $p_{T}$ bin, and within the rapidity window $|y|$ $<$ 1. The fully corrected $D^{0}$ production invariant yields are calculated using the following formula:

[TABLE]

where B.R. is the $D^{0}\rightarrow K^{-}\pi^{+}$ decay branching ratio, (3.89$\pm$0.04)% Tanabashi et al. (2018), $N^{\rm raw}$ is the reconstructed $D^{0}$ raw counts, $N_{\rm evt}$ is the total numbers of events used in this analysis, $\varepsilon_{\rm trg}$ is the centrality bias correction factor described in Sec. II.2. The raw yields need to be corrected for the TPC acceptance and tracking efficiency - $\varepsilon_{\rm TPC}$, the HFT acceptance and tracking plus topological cut efficiency - $\varepsilon_{\rm HFT}$, the particle identification efficiency - $\varepsilon_{\rm PID}$, and the finite vertex resolution correction - $\varepsilon_{\rm vtx}$.

IV.1 TPC Acceptance and Tracking Efficiency - $\varepsilon_{\rm TPC}$

The TPC acceptance and tracking efficiency is obtained using the standard STAR TPC embedding technique, in which a small amount of MC tracks (typically 5% of the total multiplicity of the real event) are processed through the full GEANT simulation Brun et al. (1994), then mixed with the raw Data Acquisition (DAQ) data in real events and reconstructed through the same reconstruction chain as the real data production. The TPC efficiency is then calculated as the ratio of the number of reconstructed MC tracks with the same offline analysis cuts for geometric acceptance and other TPC requirements to that of the input MC tracks.

Figure 9 shows the TPC acceptance and tracking efficiency $\varepsilon_{\rm TPC}$ for $D^{0}$ mesons within $|y|$ $<$ 1 in various centrality classes in this analysis. The efficiencies include the TPC and analysis acceptance cuts $p_{T}$ $>$ 0.6 GeV/$c$ and $|\eta|$ $<$ 1 as well as the TPC tracking efficiency for both pion and kaon daughters. The lower efficiency observed in central collisions is due to the increased multiplicity resulting higher detector occupancy which leads to reduced tracking efficiency in these collisions.

IV.2 HFT Acceptance, Tracking and Topological Cut Efficiency - $\varepsilon_{\rm HFT}$

IV.2.1 Data-driven Simulation

In order to fully capture the real-time detector performance, the HFT-related efficiency is obtained using a data-driven simulation method in this analysis. The performance of inclusive HFT tracks is characterized by a TPC-to-HFT matching efficiency ($\varepsilon_{\rm HFT}^{\rm match}$) and the DCA distributions with respect to the primary vertex. The HFT matching efficiency $\varepsilon_{\rm HFT}^{\rm match}$ is defined as the fraction of reconstructed TPC tracks that satisfy the requirement on the number of HFT hits. In this analysis, the requirement is to have at least one hit in each PXL and IST layer. The $\varepsilon_{\rm HFT}^{\rm match}$ includes the HFT geometric acceptance and the tracking efficiency that associate HFT hits to the extended TPC tracks. It contains the true matches for which the reconstructed tracks pick up real hits induced by these charged tracks when passing through the HFT, as well as some random fake matches. The latter has a decreasing trend as a function of $p_{T}$ as the track pointing resolution gets better at high $p_{T}$, resulting in a smaller search window when associating HFT hits in the tracking algorithm. The DCA distributions are obtained for those tracks that satisfy the HFT hit requirement. Figure 10 shows an example of the HFT matching efficiency and the 1-D projection of the $\textup{DCA}_{\textup{XY}}$ distribution for single pions at 1.0$<$ $p_{T}$ $<$ 1.2 GeV/$c$ and 0–10% central collisions. Such distributions obtained from real data are fed into a MC decay generator for $D^{0}\rightarrow K^{-}\pi^{+}$, followed by the same reconstruction of $D^{0}$ secondary vertex as in the real data analysis. The same topological cuts are then applied and the HFT related efficiency for the $D^{0}$ reconstruction is calculated.

To best represent the real detector performance, we obtain the following distributions from real data in this Monte Carlo approach:

• •

Centrality-dependent $\textup{V}_{\textup{z}}$ distributions.

• •

HFT matching efficiency $\varepsilon_{\rm HFT}^{\rm match}$, including the dependence on particle species, centrality, $p_{T}$, $\eta$, $\phi$, and $\textup{V}_{\textup{z}}$.

• •

$\textup{DCA}_{\textup{XY}}$–$\textup{DCA}_{\textup{Z}}$ 2-dimensional (2D) distributions including the dependence on particle species, centrality, $p_{T}$, $\eta$, and $\textup{V}_{\textup{z}}$.

The $\textup{DCA}_{\textup{XY}}$–$\textup{DCA}_{\textup{Z}}$ 2D distributions are the key to represent not only the true matches, but also the fake matches when connecting the TPC tracks with HFT hits. The distributions are obtained in 2D to consider the correlation between the two quantities and this is necessary and essential to reproduce the 3D DCA position distributions observed in real data. The $\phi$ dependence of these distributions are integrated over due to computing resource limits. We have checked the $\phi$ dependence (by reducing other dependencies for the same reason) and it gives a consistent result compared to the $\phi$-integrated one.

In total, there are 11 ($\phi$) $\times$ 10 ($\eta$) $\times$ 6 ($\textup{V}_{\textup{z}}$) $\times$ 9 (centrality) $\times$ 2 (particles) 1D histograms (36 $p_{T}$ bins) used for the HFT matching efficiency distributions and 5 ($\eta$) $\times$ 4 ($\textup{V}_{\textup{z}}$) $\times$ 9 (centrality) $\times$ 2 (particles) $\times$ 19 ($p_{T}$) 2D histograms (144 $\textup{DCA}_{\textup{XY}}$ $\times$ 144 $\textup{DCA}_{\textup{Z}}$ bins) for 2D DCA distributions. The number of bins chosen is optimized to balance the need of computing resources as well as the stability of the final efficiency. All dimensions have been checked so that further increase in the number of bins (in balance we need to reduce the number of bins in other dimensions) will not change the final obtained efficiency.

The procedure for this data-driven simulation package for efficiency calculation is as follows:

• •

Sample $\textup{V}_{\textup{z}}$ distribution according to the distribution obtained from the real data.

• •

Generate $D^{0}$ at the event vertex position with desired $p_{T}$ (Levy function shape fitted to $D^{0}$ spectra Adamczyk et al. (2014a)) and rapidity (flat) distributions.

• •

Propagate $D^{0}$ and simulate its decay to $K^{-}\pi^{+}$ daughters following the decay probability.

• •

Smear daughter track momentum according to the values obtained from embedding.

• •

Smear daughter track starting position according to the $\textup{DCA}_{\textup{XY}}$–$\textup{DCA}_{\textup{Z}}$ 2D distributions from the reconstructed data.

• •

Apply HFT matching efficiency according to that extracted from the reconstructed data.

• •

Perform the topological reconstruction of $D^{0}$ decay vertices with the same cuts as applied in the data analysis and calculate the reconstruction efficiency.

The DCA and HFT matching efficiency distributions used as the input in this simulation tool can be obtained from the real data or the reconstructed data in MC simulation. The latter is used when we validate this approach using the MC GEANT simulation (see Sec. IV.2.2).

This approach assumes these distributions obtained from real data are good representations for tracks produced at or close to the primary vertices. The impact of the secondary particle contribution will be discussed in Sec. IV.2.4. The approach also neglects the finite event vertex resolution contribution which will be discussed in Sec. IV.3.

Lastly in this MC approach, we also fold in the TPC efficiency obtained from the MC embedding so the following presented efficiency will be the total efficiency of $\varepsilon_{\rm TPC}\times\varepsilon_{\rm HFT}$.

IV.2.2 Validation with GEANT Simulation

In this subsection, we will demonstrate that the data-driven MC approach has been validated with the GEANT simulation plus the offline tracking reconstruction with realistic HFT detector performance to reproduce the real $D^{0}$ reconstruction efficiency. We should point out that in this validation procedure, what we are after is the efficiency difference between two calculation methods: one from the MC simulation directly, and the other one from the data-driven simulation package using the reconstructed MC simulation data as the input.

The GEANT simulation uses the HIJING Gyulassy and Wang (1994) generator as its input with $D^{0}$ particles embedded to enrich the signal statistics. The full HFT detector materials (both active and inactive) have been included in the GEANT simulation as well as the offline track reconstruction. The pileup hits in the PXL detector due to finite electronic readout time have been added to realistically represent the HFT matching efficiency and DCA distributions. The overall agreement between the GEANT simulation and real data is fairly good, as can be seen in Fig. 10. The small deviations between real data and MC simulation are not considered in the systematic uncertainty estimation since the latter is not used to calculate the absolute efficiency directly, but to validate the data-driven simulation procedure as described below.

The increase in the HFT matching efficiency at low $p_{T}$ range is due to the increased fake matches (in contrast to true HFT matches) and the efficiency stays flat in the high $p_{T}$ range. The matching efficiency includes the tracking efficiency when associating the HFT hits as well as the HFT geometric acceptance. Therefore, the ratio has a strong dependence on the event $\textup{V}_{\textup{Z}}$ and the track $\eta$. The DCA distributions used in the package are 2-dimentional distributions, as $\textup{DCA}_{\textup{XY}}$ and $\textup{DCA}_{\textup{Z}}$ are strongly correlated.

With the tuned simulation setup, we use this sample to validate our data-driven simulation approach for $D^{0}$ efficiency calculation. We follow the same procedure as described in Sec. IV.2.1 to obtain the HFT matching efficiency as well as the 2D $\textup{DCA}_{\textup{XY}}$-$\textup{DCA}_{\textup{Z}}$ distributions for primary particles from the reconstructed data in this simulation sample. Then these distributions are fed into the data-driven simulation framework to calculate the $D^{0}$ reconstruction efficiency. The calculated $D^{0}$ efficiency from the data-driven simulation framework will be compared to the real $D^{0}$ reconstruction efficiency directly obtained from the GEANT simulation sample.

To validate the data-driven simulation tool, Fig. 11 shows the comparisons of several topological variables used in the $D^{0}$ reconstruction obtained from the GEANT simulation directly and from the data-driven simulation with reconstructed GEANT simulation data as the input in the most central (0–10%) centrality and in 2 $<$ $p_{T}$ $<$ 3 GeV/$c$. The topological variables shown here are $D^{0}$ decay length, DCA between two $D^{0}$ decay daughters, $D^{0}$ DCA with respect to the collision vertex, pion DCA and kaon DCA with respect to the collision vertex. As seen in this figure, the data-driven simulation tool reproduces all of these topological distributions quite well. The agreements for the other $p_{T}$ ranges are also decent.

Figure 12 shows the $D^{0}$ reconstruction efficiency $\varepsilon_{\rm TPC}$ $\times$ $\varepsilon_{\rm HFT}$ calculated with the following two methods in this GEANT simulation. The first method is the standard calculation by applying the tracking and topological cuts for reconstructed $D^{0}$ mesons in the simulation sample. In the second method, we employ the data-driven simulation method and take the reconstructed distributions from the simulation sample as the input and then calculate the $D^{0}$ reconstruction efficiency in the data-driven simulation framework. In panel (a) of Fig. 12, efficiencies from two calculation methods agree well in the whole $p_{T}$ region in central 0–10% Au+Au collisions, and the ratio between the two is shown in panel (b). This demonstrates that the data-driven simulation framework can accurately reproduce the real $D^{0}$ reconstruction efficiency in central Au+Au collisions.

IV.2.3 Efficiency for real data

We employ the validated data-driven simulation method for the real data analysis. Figure 13 shows comparisons of the same five topological variables between $D^{0}$ signals in real data and data-driven simulated distributions with real data as the input in central 0–10% collisions for $D^{0}$ mesons at 2 $<$ $p_{T}$ $<$ 3 GeV/$c$. The real data distributions are extracted by reconstructing $D^{0}$ signals with the same reconstruction cuts as in Sec. III except for the interested topological variable to be compared. The distributions for $D^{0}$ candidates are generated by statistically subtracting the background using the side-band method from the same-event unlike-sign distributions within the $D^{0}$ mass window. The cut on the interested topological variable is loosened, but one needs to place some pre-cuts to ensure reasonable $D^{0}$ signal reconstruction for the extraction of these topological variable distributions. These pre-cuts effectively reduce the low-end reach for several topological variables, e.g. the $D^{0}$ decay length. In the data-driven simulation method, charged pion and kaon HFT matching efficiencies and 2D DCA distributions are used as the input to calculate these topological variables for $D^{0}$ signals. Figure 13 shows that in the selected ranges, the data-driven simulation method reproduces topological variables distributions of $D^{0}$ signals, which supports that this method can be reliably used to calculate the topological cut efficiency.

Figure 14 shows the HFT tracking and topological cut efficiency $\varepsilon_{\rm HFT}$ as a function of $D^{0}$ $p_{T}$ for different centrality bins obtained using the data-driven simulation method described in this section with the input distributions taken from the real data. The smaller efficiency seen in central collisions is in part because the HFT tracking efficiency is lower in higher occupancy central collisions, and in addition because we choose tighter topological cuts in central collisions for background suppression.

IV.2.4 Secondary particle contribution

In the data-driven method for obtaining the efficiency correction, inclusive pion and kaon distributions are taken from real data as the input while the validation with GEANT simulation is performed with primary particles. There is a small amount of secondary particle contribution (e.g. weak decays from $K^{0}_{S}$ and $\Lambda$) to the measured inclusive charged pion tracks.

The impact of secondary particle contribution to the charged pions is studied using the HIJING events processed through the GEANT simulation and the same offline reconstruction. The fraction of secondary pions from weak decay of strange hadrons ($K^{0}_{S}$ and $\Lambda$) to the total inclusive charged pions within ${\rm DCA}<$ 1.5 cm cut is estimated to be around 5% at pion $p_{T}$ = 0.3 GeV/$c$ and decreases to be $<2\%$ above 2 GeV/$c$. This is consistent with what was observed before in measuring the prompt charged pion spectra Adams et al. (2004). There is another finite contribution of low momentum anti-protons and anti-neutrons annihilated in the detector material and producing secondary pions. The transverse momenta of these pions are mostly around 2–3 GeV/$c$ and the fraction of total inclusive pions is $\sim$ 10–12% at $p_{T}=2\textup{-}3$ GeV/$c$ based on this simulation and contribute $\sim$ 5–8% to the HFT matching efficiency. This is obtained using the GEANT simulation with GHEISHA hadronic package. With a different hadronic package, FLUKA, the secondary pion fraction in 2–3 GeV/$c$ region is significantly reduced to be negligible. The difference between the primary pions and the inclusive pions in the HFT matching efficiency has been considered as one additional correction factor in our data-driven simulation method when calculating the final efficiency. The maximum difference with respect to the result obtained using the GHEISHA hadronic package is used as the systematic uncertainty for this source. Figure 15 shows the secondary pion contribution in Au+Au collisions with FLUKA hadronic package. Panel (a) shows the fraction of different sources for secondary tracks including the weak decays, the scattering and the $\bar{p}/\bar{n}$ annihilation in the detector material. Panel (b) shows the HFT matching efficiencies for inclusive, prompt and secondary pions. Panel (c) is the ratio of the HFT matching efficiencies between the inclusive and the primary pions from panel (b). The effect of such secondary contribution to charged kaons is found to be negligible Adams et al. (2004).

IV.3 Vertex Resolution Correction - $\varepsilon_{\rm vtx}$

In the data-driven approach, $D^{0}$ mesons are injected at the event vertex. In the real data, the reconstructed vertex has a finite resolution with respect to the real collision vertex. This may have some effect on the reconstructed $D^{0}$ signal counts after applying the topological cuts in small multiplicity events where the event vertex resolution decreases. We carry out similar simulation studies as described in Sec. IV.2.1 for other centrality bins. Figure 16 shows the Full-Width-at-Half-Maximum (FWHM) of the difference in the vertex x-position of two randomly-divided sub-events in various centrality bins between data and MC simulation. We choose the FWHM variable here as the distributions are not particularly Gaussian. The MC simulation reproduces the vertex difference distributions seen in the real data reasonably well. This gives us confidence for using this MC simulation setup to evaluate the vertex resolution correction $\varepsilon_{\rm vtx}$.

To estimate the vertex resolution effect, we embed single PYTHIA $c\bar{c}$ event into a HIJING Au+Au event, and the whole event is passed through the STAR GEANT simulation followed by the same offline reconstruction as in the real data production. The PYTHIA $c\bar{c}$ events are pre-selected to have at least one $D^{0}\rightarrow K^{-}\pi^{+}$ decay or its charge conjugate to enhance the statistics. Figure 17 shows the comparison in the obtained $D^{0}$ reconstruction efficiency between MC simulation $(black)$ and data-driven simulation using reconstructed MC data as the input $(red)$ for 20–30% (left), 50–60% (middle) and 70–80% (right) centrality bins, respectively. The bottom panels show the ratios of the efficiencies obtained from the two calculation methods. In central and mid-central collisions, the data-driven simulation method can properly reproduce the $D^{0}$ real reconstruction efficiency. This is expected since the vertex resolution is small enough so that it has negligible impact on the obtained efficiency using the data-driven simulation method. However, in more peripheral collisions, the data-driven simulation method overestimates the $D^{0}$ reconstruction efficiency as shown in the middle and right panels. The vertex resolution correction factor $\varepsilon_{\rm vtx}$, denoted in Eq. 3, has a mild $p_{T}$ dependence but strong centrality dependence as shown in Fig. 18 for $p_{T}$ = 2 and 4 GeV/$c$. The brackets denote the systematic uncertainties in the obtained correction factor $\varepsilon_{\rm vtx}$. They are estimated by changing the multiplicity range in the HIJING + GEANT simulation so that the variation in the sub-event vertex difference distributions from the real data can be covered by distributions obtained from different simulation samples. The vertex resolution corrections are applied as a function of $p_{T}$ in each individual centrality class.

IV.4 PID Efficiency - $\varepsilon_{\rm PID}$ and Doubly-mis-PID Correction

The $D^{0}$ daughter particle identification (PID) cut efficiency includes contributions from the $dE/dx$ selection cut efficiency as well as the TOF matching and $1/\beta$ cut efficiency. To best estimate the selection cut efficiency, we select the enriched kaon and pion samples from $\phi,K_{S}^{0}$ decays following the same procedure as in Shao et al. (2006); Xu et al. (2010) and obtain the mean and width in the $dE/dx$ $n\sigma_{X}$ distributions. The $dE/dx$ cut efficiencies for pion and kaon daughter tracks are calculated correspondingly. The TOF $1/\beta$ cut efficiency is determined by studying the $1/\beta$ distributions for kaons and pions in the clean separation region, namely $p_{T}$ $<$ 1.5 GeV/$c$. There is a mild dependence for the offset and width of $\Delta 1/\beta$ distributions vs. particle momentum and our selection cuts are generally wide enough to capture nearly all tracks once they have valid $\beta$ measurements. The total PID efficiency of $D^{0}$ mesons is calculated by folding the individual track TPC and TOF PID efficiencies following the same hybrid PID algorithm as implemented in the data analysis. Figure 19 shows the total PID efficiencies for $D^{0}$ reconstruction in various centrality bins. The total PID efficiency is generally high and has nearly no centrality or $p_{T}$ dependence.

When the $D^{0}$ daughter kaon track is mis-identified as a pion track and the other daughter pion track is mis-identified as a kaon track, the pair invariant mass distribution will have a bump structure around the real $D^{0}$ signal peak, but the distribution is much broader in a wide mass region due to the mis-assigned daughter particle masses. Based on the PID performance study described above, we estimate the single kaon and pion candidate track purities. After folding the realistic particle momentum resolution, we calculate the reconstructed $D^{0}$ yield from doubly mis-identified pairs (double counting) underneath the real $D^{0}$ signal and the double counting fraction is shown in Fig. 20. The black markers show the fraction by taking all doubly mis-identified pairs in the $D^{0}$ mass window while the blue markers depict it with an additional side-band (SB) subtraction. The latter is used as a correction factor to the central values of reported $D^{0}$ yields while the difference between the black and blue symbols is considered as the systematic uncertainty in this source. The double counting fraction is below 10% in all $p_{T}$ bins, and there is little centrality dependence.

Figure 21 shows the total $D^{0}$ reconstruction efficiency from different centrality classes in Au+Au collisions including all of the individual components discussed above.

V Systematic Uncertainties

The systematic uncertainty on the final measured $D^{0}$ $p_{T}$ spectra can be categorized as the uncertainty of the raw $D^{0}$ yield extraction and the uncertainty of efficiencies and corrections.

The uncertainty of the raw yield extraction is estimated by a) changing the $D^{0}$ raw yield counting method from the Gaussian fit to histogram bin counting. b) varying invariant mass ranges for fit and for side bands and c) varying background estimation from mixed-event and like-sign methods. The maximum difference between these scenarios is then converted to the standard deviation and added to the systematic uncertainties. It is the smallest in the mid-$p_{T}$ bins due to the best signal significance and grows at both low and high $p_{T}$. The double counting contribution in the $D^{0}$ raw yield due to mis-PID is included as another contribution to the systematic uncertainty for the $D^{0}$ raw yield extraction as described in Sec. IV.4.

The uncertainty of the TPC acceptance and efficiency correction $\varepsilon_{\rm TPC}$ is estimated via the standard procedure in STAR by comparing the TPC track distributions between real data and the embedding data. It is estimated to be $\sim$5–7% for 0–10% collisions and $\sim$5–8% for 60–80% collisions, and is correlated for different centralities and $p_{T}$ regions.

The uncertainty of the PID efficiency correction is estimated by varying the PID selection cuts and then convoluting to the final corrected $D^{0}$ yield.

To estimate the uncertainty of the HFT tracking and topological cut efficiency correction $\varepsilon_{\rm HFT}$, we employ the following procedures: a) We vary the topological variable cuts such that the $D^{0}$ $\varepsilon_{\rm HFT}$ is changed to 50% and 150% from the nominal (default) efficiency and compare the efficiency–corrected final $D^{0}$ yields. The maximum difference between the two scenarios is then added to the systematic uncertainties. b) We also vary the lower threshold cut on the daughter $p_{T}$ between 0.3 to 0.6 GeV/$c$ and the maximum difference in the final corrected $D^{0}$ yield is also included in the systematic uncertainties. c) We add the systematic uncertainty due to limitation of the data-driven simulation approach, $\sim$5%, and the impact of the secondary particles, $\sim$2%, to the total $\varepsilon_{\rm HFT}$ systematic uncertainty.

With the corrected $D^{0}$ transverse momentum spectra, the nuclear modification factor $R_{\rm CP}$ is calculated as the ratio of $N_{\rm bin}$–normalized yields between central and peripheral collisions, as shown in the following formula:

[TABLE]

The systematic uncertainties in the raw signal extraction in central and peripheral collisions are propagated as they are uncorrelated, while the systematic uncertainties from the other sources are correlated or partially correlated in contributing to the measured $D^{0}$ yields. To best consider these correlations, we vary selection cuts simultaneously in central and peripheral collisions, and the difference in the final extracted $R_{\rm CP}$ value is then directly counted as systematic uncertainties in the measured $R_{\rm CP}$.

The nuclear modification factor $R_{\rm AA}$ is calculated as the ratio of $N_{\rm bin}$–normalized yields between Au+Au and $p$+$p$ collisions. The baseline for $p$+$p$ collisions is chosen the same as in Ref. Adamczyk et al. (2014a). The uncertainties from the $p$+$p$ reference dominates the systematic uncertainty for $R_{\rm AA}$. They include the 1$\sigma$ uncertainty from the Levy function fit to the measured spectrum and the difference between Levy and power-law function fits for extrapolation to low and high $p_{T}$, expressed as one standard deviation.

With the corrected $D^{0}$ and $\overline{D}^{0}$ transverse momentum spectra, the $\overline{D}^{0}/D^{0}$ ratio is calculated as a function of the transverse momentum. The systematic uncertainties in the raw signal extraction for $\overline{D}^{0}$ and $D^{0}$ are propagated as they are uncorrelated, while the systematic uncertainties from the other sources are correlated or partially correlated in contributing to the measured $\overline{D}^{0}/D^{0}$ ratio. As in the $R_{\rm CP}$ systematic uncertainty estimation, we vary selection cuts simultaneously for $D^{0}$ and $\overline{D}^{0}$, and the difference in the final extracted $\overline{D}^{0}/D^{0}$ value is then directly counted as systematic uncertainties for the measured $\overline{D}^{0}/D^{0}$ ratio.

Table 5 summarizes the systematic uncertainties and their contributions, in percentage, on the $D^{0}$ invariant yield in 0–10% and 60–80% collisions and $R_{\rm CP}$(0–10%/60–80%). In the last column we also comment on the correlation in $p_{T}$ for each individual source. Later when reporting $p_{T}$–integrated yields or $R_{\rm CP}$, systematic uncertainties are calculated under the following considerations: a) for $p_{T}$ uncorrelated sources, we take the quadratic sum of various $p_{T}$ bins; b) for sources that are largely correlated in $p_{T}$, we take the arithmetic sum as a conservative estimate.

VI Results and Discussion

VI.1 $p_{T}$ Spectra and Integrated Yields

Figure 22 shows the efficiency–corrected $D^{0}$ invariant yield at mid-rapidity ($|y|<1$) vs. $p_{T}$ in 0–10%, 10–20%, 20–40%, 40–60% and 60–80% Au+Au collisions. $D^{0}$ spectra in some centrality bins are scaled with arbitrary factors indicated on the figure for clarity. Dashed and solid lines depict fits to the spectra with the Levy function:

[TABLE]

where $m_{0}$ is the $D^{0}$ mass (1.864 GeV/$c^{2}$) and $dN/dy$, $T$ and $n$ are free parameters. The Levy function fit describes the $D^{0}$ spectra nicely in all centrality bins in our measured $p_{T}$ region.

We compare our new measurements with previous measurements using the STAR TPC only. The previous measurements are recently corrected after fixing errors in the TOF PID efficiency calculation Adamczyk et al. (2014a). Figure 23 shows the $p_{T}$ spectra comparison in 0–10%, 10-40% and 40–80% centrality bins in panel (a) and the ratios to the Levy fit functions in panels (b), (c), and (d), respectively. The new measurement with the HFT detector shows a nice agreement with the measurement without the HFT, but with significantly improved precision.

The measured $D^{0}$ spectra cover a wide $p_{T}$ region which allows us to extract the $p_{T}$–integrated total $D^{0}$ yield at mid-rapidity with good precision.   Figure  24 shows the $p_{T}$–integrated cross section for $D^{0}$ production per nucleon-nucleon collision $d\sigma^{NN}/dy|_{y=0}$ from different centrality bins for the full $p_{T}$ range shown in the top panel and for $p_{T}$ $>$ 4 GeV/$c$ shown in the bottom panel. The result from the previous $p$+$p$ measurement is also shown in both panels Adamczyk et al. (2012).

While $d\sigma^{NN}/dy|_{y=0}$ for $p_{T}$ $>$ 4 GeV/$c$ shows a clear decreasing trend from peripheral to mid-central and central collisions, that for the full $p_{T}$ range shows approximately a flat distribution as a function of $N_{\rm part}$, though the systematic uncertainty in the 60–80% centrality bin is a bit large. The values for the full $p_{T}$ range in mid-central to central Au+Au collisions are smaller than that in $p$+$p$ collisions with $\sim 1.5\sigma$ effect considering the large uncertainties from the $p$+$p$ measurements. The total charm quark yield in heavy-ion collisions is expected to follow the number-of-binary-collision scaling since charm quarks are believed to be predominately created at the initial hard scattering before the formation of the QGP at RHIC energies. However, the cold nuclear matter (CNM) effect including shadowing could also play an important role. In addition, hadronization through coalescence has been suggested to potentially modify the charm quark distribution in various charm hadron states which may lead to the reduction in the observed $D^{0}$ yields in Au+Au collisions Greco et al. (2004) (as seen in Fig. 24). For instance, hadronization through coalescence can lead to an enhancement of the charmed baryon $\Lambda_{c}^{+}$ yield relative to $D^{0}$ yield Oh et al. (2009); Zhao et al. (2018); Plumari et al. (2018), and together with the strangeness enhancement in the hot QCD medium, can also lead to an enhancement in the charmed strange meson $D_{s}^{+}$ yield relative to $D^{0}$ He et al. (2013); Zhao et al. (2018); Plumari et al. (2018). Therefore, determination of the total charm quark yield in heavy-ion collisions will require measurements of other charm hadron states over a broad momentum range.

VI.2 Collectivity

VI.2.1 $m_{T}$ Spectra

Transverse mass spectra can be used to study the collectivity of produced hadrons in heavy-ion collisions. Figure 25 shows the $D^{0}$ invariant yield at mid-rapidity ($|y|$ $<$ 1) vs. transverse kinetic energy ($m_{T}$-$m_{0}$) for different centrality classes, where $m_{T}=\sqrt{p_{T}^{2}+m_{0}^{2}}$ and $m_{0}$ is the $D^{0}$ meson mass at rest. Solid and dashed black lines depict exponential function fits inspired by thermal models to data in various centrality bins up to $m_{T}-m_{0}$ $<$ 3 GeV/$c^{2}$ using the fit function shown below:

[TABLE]

Such a method has been often used to analyze the particle spectra and to understand kinetic freeze-out properties from the data in heavy–ion collisions Kaneta (1999); Adams et al. (2005).

A power-law function (shown below) is also used to fit the spectrum in the 60–80% centrality bin:

[TABLE]

where $dN/dy$, $\langle p_{T}\rangle$, and $n$ are three free parameters.

The power-law function fit shows a good description of the 60–80% centrality data indicating that the $D^{0}$ meson production in this peripheral bin is close to the expected feature of perturbative QCD. The $D^{0}$ meson spectra in more central collisions can be well described by the exponential function fit at $m_{T}-m_{0}$ $<$ 3 GeV/$c^{2}$ suggesting the $D^{0}$ mesons have gained collective motion in the medium evolution in these collisions.

The obtained slope parameter $T_{\rm eff}$ for $D^{0}$ mesons is compared to other light and strange hadrons measured at RHIC. Figure 26 summarizes the slope parameter $T_{\rm eff}$ for various identified hadrons ($\pi^{\pm}$, $K^{\pm}$, $p$/$\bar{p}$, $\phi$, $\Lambda$, $\Xi^{-}$, $\Omega$, $D^{0}$ and $J/\psi$) in central Au+Au collisions at $\sqrt{s_{{}_{\rm NN}}}{\rm=200\,GeV}$ Adams et al. (2004); Abelev et al. (2007a); Adams et al. (2007); Adamczyk et al. (2014b). Point-by-point statistical and systematic uncertainties are added as a quadratic sum when performing these fits. All fits are performed up to $m_{T}-m_{0}<1\,\rm{GeV}/c^{2}$ for $\pi,\ K,\ p$, $<2$ GeV/$c^{2}$ for $\phi,\ \Lambda,\ \Xi$, and $<3$ GeV/$c^{2}$ for $\Omega,\ D^{0},\ J/\psi$, respectively.

The slope parameter $T_{\rm eff}$ in a thermalized medium can be characterized by the random (generally interpreted as a kinetic freeze-out temperature $T_{\rm fo}$) and collective (radial flow velocity $\langle\beta_{T}\rangle$) components with a simple relation Adams et al. (2005); Csorgo and Lorstad (1996); Kolb and Heinz (2003):

[TABLE]

Therefore, $T_{\rm eff}$ will show a linear dependence as a function of particle mass $m_{0}$ with a slope that can be used to characterize the radial flow collective velocity.

The data points clearly show two different systematic trends. $\pi,\ K,\ p$ data points follow one linear dependence while $\phi,\ \Lambda,\ \Xi^{-},\ \Omega^{-},\ D^{0}$ data points follow another linear dependence, as represented by the dashed lines shown in Fig. 26. Particles, such as, $\pi,\ K,\ p$ gain radial collectivity through the whole system evolution, therefore the linear dependence exhibits a larger slope. On the other hand the linear dependence of $\phi,\ \Lambda,\ \Xi^{-},\ \Omega^{-},\ D^{0}$ data points shows a smaller slope indicating these particles may freeze out from the system earlier, and therefore receive less radial collectivity.

VI.2.2 Blast-wave fit

The Blast-Wave (BW) model is extensively used to study the particle kinetic freeze-out properties Adams et al. (2004); Adamczyk et al. (2017b). Assuming a hard-sphere uniform particle source with a kinetic freeze-out temperature $T_{\rm kin}$ and a transverse radial flow velocity $\beta$, the particle transverse momentum spectral shape is given by Schnedermann et al. (1993):

[TABLE]

where $\rho=\tanh^{-1}\beta$, and $I_{0}$ and $K_{1}$ are the modified Bessel functions. The flow velocity profile is taken as:

[TABLE]

where $\beta_{s}$ is the maximum velocity at the surface and $r/R$ is the relative radial position in the thermal source. The choice of $R$ only affects the overall spectrum magnitude while the spectrum shape constrains the three free parameters $T_{\rm kin}$, $\langle\beta\rangle=2/(2+n)\beta_{s}$, and $n$.

In the modified Tsallis Blast-Wave (TBW) model, an additional parameter $q$ is introduced to account for the non-equilibrium feature in a system Tang et al. (2009). The particle transverse momentum spectral shape is then described by:

[TABLE]

When $q-1$ approaches zero, the TBW function returns to the regular BW function shown in Eq. 9.

Figure 27 shows the Blast-Wave and Tsallis Blast-Wave fits to the data in different centrality bins. The parameter $n$ in these fits is fixed to be 1 due to the limited number of data points and is inspired by the fit result for light-flavor hadrons ($\pi,K$ and $p$) Tang et al. (2009). The $p_{T}$ range in the BW fits is restricted to be less than 3$m_{0}$ where $m_{0}$ is the rest mass of $D^{0}$ mesons.

Figure 28 summarizes the fit parameters $T_{\rm kin}$ vs. $\langle\beta\rangle$ from the Blast-Wave model fits to different groups of particles: black markers for the simultaneous fit to $\pi,\ K,\ p$; red markers for the simultaneous fit to $\phi,\ \Xi^{-}$ and blue markers for the fit to $D^{0}$.   The data points for each group of particles represent the fit results from different centrality bins with the most central data point at the largest $\langle\beta\rangle$ value.   Similar as in the fit to the $m_{T}$ spectra, point-by-point statistical and systematic uncertainties are added in quadrature when performing the fit. The fit results for $\pi,\ K,\ p$ are consistent with previously published results Tang et al. (2009). The fit results for multi-strangeness particles $\phi,\ \Xi^{-}$, and for $D^{0}$ show much smaller mean transverse velocity $\langle\beta\rangle$ and larger kinetic freeze-out temperature, suggesting these particles decouple from the system earlier and gain less radial collectivity compared to light hadrons. The resulting $T_{\rm kin}$ parameters for $\phi,\ \Xi^{-}$ and for $D^{0}$ are close to the pseudocritical temperature $T_{c}$ calculated from a lattice QCD calculation at zero baryon chemical potential Bazavov et al. (2012), indicating negligible contribution from the hadronic stage to the observed radial flow of these particles.  Therefore, the collectivity that $D^{0}$ mesons obtain is mostly through the partonic stage re-scatterings in the QGP phase.

Table 6 lists the fitting parameters, $\langle\beta\rangle$ and $(q-1)$ for the $D^{0}$ data in different centralities. Results show a similar trend as the regular BW fit, i.e. the most central data point is located at the largest $\langle\beta\rangle$ value. The $(q-1)$ parameter in TBW, which characterizes the degree of non-equilibrium in a system, is found to be close to zero, indicating that the system is approaching thermalization in these collisions.

VI.3 Nuclear Modification Factors - $R_{\rm CP}$ and $R_{\rm AA}$

Figure 29 shows the calculated $D^{0}$ $R_{\rm CP}$ (see Eq. 4) with the 60–80% peripheral bin as the reference for different centrality bins 0–10%, 10–20%, 20–40% and 40–60%; the results are compared to other light and strange flavor mesons. The grey bands around unity depict the vertex resolution correction uncertainty on the measured $D^{0}$ data points, mostly originating from the 60–80% reference spectrum. The dark and light green boxes around unity on the right side indicate the global $N_{\rm bin}$ systematic uncertainties for the 60–80% centrality bin and for the corresponding centrality bin in each panel. The global $N_{\rm bin}$ systematic uncertainties should be applied to the data points of all particles in each panel.

The measured $D^{0}$ $R_{\rm CP}$ in central 0–10% collisions shows a significant suppression at $p_{T}>$ 5 GeV/$c$. The suppression level is similar to that of light and strange flavor mesons and the $R_{\rm CP}$ suppression gradually decreases when moving from central collisions to mid-central and peripheral collisions. The $D^{0}$ $R_{\rm CP}$ for $p_{T}$ $<$ 4 GeV/$c$ is consistent with no suppression, in contrast to light-flavor hadrons. Comparisons to dynamic model calculations for the $D^{0}$ $R_{\rm CP}$ will be discussed in Sec. VI.5.

The precision of the 60–80% centrality spectrum is limited due to the large systematic uncertainty in determining the $N_{\rm bin}$ based on the MC Glauber model. Figure 30 shows the $D^{0}$ $R_{\rm CP}$ for different centralities as a function of $p_{T}$ with the 40–60% centrality spectrum as the reference. The grey bands around unity in each panel represent the systematic uncertainties due to the vertex resolution contribution from the 40–60% centrality. The green boxes around unity depict the global $N_{\rm bin}$ systematic uncertainties for the 40–60% centrality bin and for each corresponding centrality bin. As a comparison, $R_{\rm CP}$ of charged pions, $K_{s}^{0}$ and $\phi$ in the corresponding centralities are also plotted in each panel. With much smaller systematic uncertainties, the observations seen before using the 60–80% centrality spectrum as the reference still hold.

Figure 31 shows the calculated $D^{0}$ $R_{\rm AA}$ (see Eq. 1) with the $p$+$p$ measurement Adamczyk et al. (2012) as the reference for different centrality bins 0–10% (a), 10–40% (b) and 40–80% (c), respectively. The new $R_{\rm AA}$ measurements are also compared to the previous Au+Au measurements using the STAR TPC after the recent correction Adamczyk et al. (2014a). The $p$+$p$ $D^{0}$ reference spectrum is updated using the latest global analysis of charm fragmentation ratios from Lisovyi et al. (2016) and also by taking into account the $p_{T}$ dependence of the fragmentation ratio between $D^{0}$ and $D^{*\pm}$ from PYTHIA. The new measurement with the HFT detector shows a nice agreement with the measurement without the HFT. The brackets on the data points depict the total systematic uncertainty dominated by the uncertainty in the $p$+$p$ reference spectrum. The first two and last two data points are empty circles indicating those are calculated with an extrapolated $p$+$p$ reference. The light and dark green boxes around unity on the right side indicate the global $N_{\rm bin}$ systematic uncertainties for the corresponding centrality bin in each panel and the total inelastic cross section uncertainty in $p$+$p$ collisions.

The measured $D^{0}$ $R_{\rm AA}$ in central (0–10%) and mid-central (10–40%) collisions show a significant suppression at the high $p_{T}$ range which reaffirms the strong interactions between charm quarks and the medium, while the new Au+Au data points from this analysis contain much improved precision. Figure 32 shows the $D^{0}$ $R_{\rm AA}$ in the 0–10% most central collisions compared to that of (a) average $D$-meson from ALICE and (b) charged hadrons from ALICE and $\pi^{\pm}$ from STAR Adam et al. (2016); Abelev et al. (2013b, 2007b). The $D^{0}$ $R_{\rm AA}$ from this measurement is comparable to that from the LHC measurements in Pb+Pb collisions at $\sqrt{s_{{}_{\rm NN}}}$ = 2.76 TeV despite the large energy difference between these measurements. The comparison to that of light hadrons shows a similar suppression at high $p_{T}$, while in the intermediate range, $D^{0}$ mesons seem to be less suppressed. From low to intermediate $p_{T}$ region, the $D^{0}$ $R_{\rm AA}$ in the central 0-10% collisions shows a characteristic bump structure that is consistent with the expectation from model predictions that charm quarks gain sizable collective motion during the medium evolution. The large uncertainty in the $p$+$p$ baseline need to be further reduced before making more quantitative conclusions.

VI.4 $\overline{D}^{0}$ and $D^{0}$ spectra and double ratio

Figure 33 shows the $p_{T}$ spectra of $\overline{D}^{0}$ and $D^{0}$ mesons separately in 0–10%, 10–20%, 20–40%, 40–60% and 60–80% centrality bins. Figure 34 shows the $\overline{D}^{0}$/$D^{0}$ ratio in various centrality bins. Dashed lines represent constant function fits to the $\overline{D}^{0}$/$D^{0}$ ratio in each centrality bin by combining the point-by-point statistical and systematic uncertainties. The $\overline{D}^{0}/D^{0}$ ratio has a small but significant deviation from unity in central and mid-central collisions. Table 7 lists the fitted results for the $\overline{D}^{0}$/$D^{0}$ ratio from various centralities. In the most central collisions, $\overline{D}^{0}$ yield is higher than the $D^{0}$ yield by $\sim$4.9$\sigma$. The total charm quark and anti-charm quark should be conserved since they are created in pairs. A thermal model calculation predicts that the $\Lambda_{c}^{-}$/$\Lambda_{c}^{+}$ ratio will be smaller than unity at RHIC due to the finite baryon density Andronic et al. (2003). This will then yield more $\overline{D}^{0}$ mesons formed than $D^{0}$ mesons in Au+Au collisions at RHIC. To verify the total charm quark conservation, one would need precise measurements of $D^{+}$/$D^{-}$, $D_{s}^{+}$/$D_{s}^{-}$ as well as $\Lambda_{c}^{+}$/$\Lambda_{c}^{-}$ ratios in the future.

VI.5 Comparison to Models

Over the past several years, there have been rapid developments in the theoretical calculations on the charm hadron production Rapp et al. (2018); Cao et al. (2018). Here we compare our measurements to several recent calculations based on the Duke model and the Linearized Boltzmann Transport (LBT) model Cao et al. (2016); Cao ; Xu et al. (2018).

The Duke model Cao et al. (2015); Xu et al. (2018) uses a Langevin stochastic simulation to trace the charm quark propagation inside the QGP medium. Both collisional and radiative energy losses are included in the calculation and charm quarks are hadronized via a hybrid approach combining both coalescence and fragmentation mechanisms. The bulk medium is simulated using a viscous hydrodynamic evolution followed by a hadronic cascade evolution using the UrQMD model Bleicher et al. (1999). The charm quark interaction with the medium is characterized using a temperature and momentum-dependent diffusion coefficient. The medium parameters have been constrained via a statistical Bayesian analysis by fitting the previous experimental data of $R_{\rm AA}$ and $v_{2}$ of light, strange and charm hadrons Xu et al. (2018). The extracted charm quark spatial diffusion coefficient at zero momentum $2\pi TD_{s}|_{p=0}$ is about 1–3 near $T_{c}$ and exhibits a positive slope for its temperature dependence above $T_{c}$.

The Linearized Boltzmann Transport (LBT) calculation Cao et al. (2016) extends the LBT approach developed before to include both light and heavy flavor parton evolution in the QGP medium. The transport calculation includes all $2\rightarrow 2$ elastic scattering processes for collisional energy loss and the higher-twist formalism for medium induced radiative energy loss. It uses the same hybrid approach as in the Duke model for charm quark hadronization. The heavy quark transport is coupled with a 3D viscous hydrodynamic evolution which is tuned for light flavor hadron data. The charm quark spatial diffusion coefficient is estimated via the equation $2\pi TD_{s}=8\pi/\hat{q}$ ($\hat{q}$ is the quark transport coefficient due to elastic scatterings) at parton momentum $p$ $=$ 10 GeV/$c$. The $2\pi TD_{s}$ is $\sim$3 at $T_{c}$ and increases to $\sim$6 at $T=500$ MeV Cao .

Figures 35 and 36 show the measured $D^{0}$ $R_{\rm CP}$ compared to the Duke and LBT model calculations with the 60–80% and 40–60% reference spectra respectively. The $R_{\rm CP}$ curves from these models are calculated based on the $D^{0}$ spectra provided by each group Cao et al. (2016); Cao ; Xu et al. (2018). The Duke model did not calculate the spectra in the 60–80% centrality bin due to a concern about the viscous hydrodynamic implementation. In Fig. 32 for the most central collisions, there are also calculations for the $D^{0}$ $R_{\rm AA}$ from the Duke and LBT groups, respectively. These two models also have the predictions for the $D^{0}$ $v_{2}$ measurements for Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{200\,GeV}}$ Adamczyk et al. (2017a). Both model calculations match our new measured $R_{\rm CP}$ data well. The much improved precision of these new measurements are expected to further constrain the theoretical model uncertainties in these calculations.

VII Summary

In summary, we report the improved measurement of $D^{0}$ production invariant yield at mid-rapidity ($|y|$ $<$ 1) in Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{200\,GeV}}$ with the STAR HFT detector.  $D^{0}$ invariant yields are presented as a function of $p_{T}$ in various centrality classes.  The $p_{T}$–integrated $D^{0}$ production cross section per nucleon-nucleon collisions in mid-central and central Au+Au collisions seem to be smaller than that measured in $p$+$p$ collisions by 1.5$\sigma$, indicating that CNM effects and/or hadronization through quark coalescence may play an important role in Au+Au collisions. This calls for precise measurements of $D^{0}$ production in $p/d$+A collisions to understand the CNM effects as well as other charm hadron states in heavy-ion collisions to better constrain the total charm quark yield.

The $\overline{D}^{0}$ yield is observed to be higher than the $D^{0}$ yield in the most central collisions, by $\sim$4.9$\sigma$ on average. This is potentially consistent with the expectation, due to the finite baryon density of the system at RHIC, that the $\Lambda_{c}^{-}$/$\Lambda_{c}^{+}$ ratio should be smaller than unity, which would result in larger $\overline{D}^{0}$ yield than $D^{0}$.

The $D^{0}$ spectra at low $p_{T}$ and $m_{T}$ regions are fit by the exponential function and the (Tsallis) Blast-Wave model to study the $D^{0}$ meson radial collectivity. The slope parameter extracted from the exponential function fit for $D^{0}$ mesons follows the same linearly increasing trend vs. particle mass as $\phi$, $\Lambda$, $\Xi^{-}$, $\Omega^{-}$ particles, but is different from the trend of $\pi,K$ and $p$ particles. The extracted kinetic freeze-out temperature and transverse velocity from the Blast-Wave model fit are comparable to the fit results of $\phi,\Xi^{-}$ multi-strange-quark hadrons, but different from those of $\pi,K$ and $p$.  These observations suggest that $D^{0}$ hadrons show a radial collective behavior with the medium, but freeze out from the system earlier and gain less radial collectivity compared to $\pi,K$ and $p$ particles. This observation is consistent with collective behavior observed in $v_{2}$ measurements. The fit results also suggest that $D^{0}$ mesons have similar kinetic freeze-out properties as multi-strange-quark hadrons $\phi,\Xi^{-}$.

The nuclear modification factors $R_{\rm CP}$ of $D^{0}$ mesons are presented with both 60–80% and 40–60% centrality spectra as the reference, respectively. The $D^{0}$ $R_{\rm CP}$ is significantly suppressed at high $p_{T}$ and the suppression level is comparable to that of light hadrons at $p_{T}>5\,\rm{GeV}/c$, re-affirming our previous observation Adamczyk et al. (2014a). This indicates that charm quarks lose significant energy when traversing through the hot QCD medium. The $D^{0}$ $R_{\rm CP}$ is above the light hadron $R_{\rm CP}$ at low $p_{T}$. We compare our $D^{0}$ $R_{\rm CP}$ measurements to two recent theoretical model calculations from LBT and the Duke groups. These two models have the $2\pi TD_{s}$ value around 1-3 near $T_{c}$ and agree with our new $R_{\rm CP}$ measurements. The nuclear modification factor $R_{\rm AA}$ of $D^{0}$ mesons in 0-10% central Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{200\,GeV}}$ is comparable to that from the ALICE measurement in Pb+Pb at ${\sqrt{s_{\rm NN}}=\rm{2.76\,TeV}}$. At $p_{T}$ $<$ 5 GeV/$c$, the $D^{0}$ $R_{\rm AA}$ shows a characteristic bump structure. Model calculations that predict sizable collective motion for charm quarks during the medium evolution can qualitatively describe our data. We expect the new data points with much improved precision can be used in the future to further constrain our understanding of the charm-medium interactions as well as to better determine the medium transport parameter.

VIII Acknowledgement

We thank the RHIC Operations Group and RCF at BNL, the NERSC Center at LBNL, and the Open Science Grid consortium for providing resources and support. This work is supported in part by the Office of Nuclear Physics within the U.S. DOE Office of Science, the U.S. National Science Foundation, the Ministry of Education and Science of the Russian Federation, National Natural Science Foundation of China, Chinese Academy of Science, the Ministry of Science and Technology of China and the Chinese Ministry of Education, the National Research Foundation of Korea, GA and MSMT of the Czech Republic, Department of Atomic Energy and Department of Science and Technology of the Government of India; the National Science Centre of Poland, National Research Foundation, the Ministry of Science, Education and Sports of the Republic of Croatia, RosAtom of Russia and German Bundesministerium fur Bildung, Wissenschaft, Forschung and Technologie (BMBF) and the Helmholtz Association.