# Modeling of heavy-flavor pair correlations in Au-Au collisions at 200 A GeV at the BNL Relativistic Heavy Ion Collider

## Abstract

We study the nuclear modification of angular and momentum correlations between heavy quark pairs in ultra-relativistic heavy-ion collisions. The evolution of heavy quarks inside the thermalized medium is described via a modified Langevin approach that incorporates both elastic and inelastic interactions with the medium constituents. The spacetime evolution of the fireball is obtained from a (2+1)-dimensional viscous hydrodynamics simulation. The hadronization of heavy quarks is performed utilizing a hybrid model of fragmentation and coalescence. Our results show that the nuclear modification of the transverse momentum imbalance of pairs reflects the total energy loss experienced by the heavy quarks and may help us probe specific regions of the medium. The angular correlation of heavy flavor pairs, especially in the low to intermediate transverse momentum regime, is sensitive to the detailed energy loss mechanism of heavy quarks inside the QGP.

## I Introduction

It has been shown that a new state of matter, commonly referred to as the strongly interacting quark-gluon plasma (sQGP), can be produced in ultra-relativistic nucleus-nucleus collisions, such as those performed at the CERN Large Hadron Collider (LHC) and the BNL Relativistic Heavy-ion Collider (RHIC). This highly-excited nuclear matter exhibits many remarkable properties (1); (2); (3); (4); (5); (6); (7); (8), among most important observations is the strong anisotropic collective flow exhibited by the hot and dense fireball and the final emitted hadrons (9); (10); (11), which can be nicely described by relativistic hydrodynamics simulations (12); (13); (14); (15); (16). Another result is the strong modification of large transverse momentum () hadrons and jets as compared to the reference of proton-proton collisions (17); (18); (19). Many high- phenomena and jet-related observables can be well understood from parton energy loss which originates from the interaction between the propagating hard partons and the hot and dense QCD matter (20); (21).

Heavy quarks and heavy flavor mesons are of particular importance. Because of their large masses, heavy quarks are expected to lose less energy than light flavor partons due to the so-called “dead cone” effect, i.e., the phase space of collinear radiation is suppressed by the kinematics compared to light partons (22). In contrast, experimental measurements have shown that the nuclear modification factors and the anisotropy coefficients for heavy flavor mesons are comparable to light flavor hadrons (23); (24); (25); (26); (27); (28).

Many phenomenological models have been developed to investigate the propagation and energy loss of heavy quarks inside the hot and dense nuclear medium (29); (30); (31); (32); (33); (34); (35); (36). In our earlier work (37); (38); (39); (35); (36), we have developed a transport model to simulate the evolution dynamics of heavy flavors in relativistic heavy-ion collisions. Both collisional and radiative energy losses of heavy quarks are included in our model via an improved Langevin approach, where the gluon radiation is treated as a recoil force exerted on the propagating heavy quarks. The recoil force is obtained from the medium-induced gluon radiation spectrum as calculated from perturbative QCD energy loss formalisms. The hadronization process of heavy quarks after passing through the QGP is simulated utilizing a hybrid model of fragmentation and coalescence, where the momentum dependence of the relative probability of the two hadronization mechanisms is calculated according to the Wigner functions of the heavy hadrons. The interactions of heavy flavor mesons with the hadron gas were also incorporated in Ref. (36) by utilizing the Ultra-relativistic Quantum Molecular Dynamic Model (UrQMD) (40). With the above setup, we can provide a reasonable description of the nuclear modification of heavy flavor mesons at the LHC and RHIC (36).

In addition to single inclusive observables, correlations of final state particles have played a very important role in the study of jet quenching and parton energy loss. In the light flavor sector, there have been a wealth of experimental measurements and theoretical studies on final state correlations, such as back-to-back dihadron and photon-hadron correlations (41); (42); (43); (44); (45); (46), as well as the momentum imbalance of back-to-back jet pairs and photon-jet pairs (47); (48); (49); (50); (51); (52); (53). In the heavy flavor sector, the angular correlations of heavy flavor pairs have been studied in Ref. (54); (55); (56); (57); (58); (34); (59); (60); (61), and the momentum imbalance of heavy flavor pairs has been investigated in Ref. (62); (63). In this article, we present a detailed study of both angular and momentum correlations for heavy flavor pairs using our heavy flavor transport model built in Ref. (37); (38); (39); (35). The dependence on the collisions centrality and the transverse momentum cuts are investigated for both angular and momentum correlations. We will show that while the momentum imbalance of back-to-back heavy flavor meson pairs is mostly sensitive to the total energy loss of heavy quarks in the dense QGP medium, the angular correlations probe the transverse momentum broadening of heavy quark propagation and are sensitive to different parton energy loss mechanisms, especially in the low transverse momentum region.

This article is organized as follows. In Sec.II, we present our theoretical framework for simulating the in-medium evolution and hadronization of heavy quarks in relativistic heavy-ion collisions. In Sec.III, we provide a detailed investigation on the angular correlations of heavy flavor pairs and show their sensitivity to different energy loss mechanisms. The transverse momentum imbalance between heavy flavor pairs is presented in Sec.IV. Sec.V contains our conclusions.

## Ii Heavy flavor dynamics

In the limit of multiple scatterings, the momentum exchange of heavy quarks in each scattering with the medium constituents is small. In this case, the evolution of heavy quarks inside a thermalized medium via multiple quasi-elastic scatterings can be described by the Langevin equation. However, as one moves to the high region, quasi-elastic scatterings alone are no longer sufficient to describe the heavy quark evolution and energy loss in hot and dense nuclear medium, and the contribution from medium-induced gluon radiation has to be included. To incorporate both energy loss processes, we follow our previous framework (35); (36) and modify the classical Langevin equation as follows:

(1) |

In the above equation, the first two terms on the right hand side are the drag force and the thermal random force experienced by heavy quarks when they scatter with the medium constituents; they are inherited from the classical Langevin equation. In this work, we assume the fluctuation-dissipation theorem still holds, , where is the momentum space diffusion coefficient of heavy quark defined as . On the right hand of the above equation, a third term is introduced to describe the recoil force experienced by heavy quarks when they radiate gluons. The probability for a heavy quark to emit a gluon within the time interval is evaluated according to the average number of radiated gluons:

(2) |

where is chosen to be sufficiently small to ensure . is the gluon radiation spectrum and in our work is obtained from the higher-twist energy loss calculation (64); (65); (66):

(3) |

In the equation, is the fractional energy carried by the radiated gluon, is its transverse momentum, is its formation time and is the vacuum splitting function. is known as the gluon transport coefficient, which is related to via . With this setup, there is only one free parameter in our model. In the literature, one often quotes the spatial diffusion coefficient of heavy quark to characterize its interaction strength with the medium; it is related to the momentum drag and diffusion as: . We will follow such practice and quote the values of in the following discussions.

Using our transport model, we can simulate the propagation of heavy quarks inside the QGP matter created in relativistic heavy-ion collisions. The spacetime evolution of the local temperature and fluid velocity of the dense nuclear matter is simulated with a (2+1)-dimensional viscous hydrodynamic model (67). In each time step, we boost the heavy quark into the local rest frame of the fluid cell in which the energy and momentum of a given heavy quark is updated according to our modified Langevin equation before it is boosted back to the global (laboratory) frame.

Before the start of the hydrodynamical evolution (at fm/c), the initial entropy density is obtained via a Monte-Carlo (MC) Glauber model. The same MC Glauber model is utilized to initialize the spatial distribution of initial produced heavy quarks. The initial momentum space distribution of heavy quarks is calculated in leading order perturbative QCD when calculating the production and nuclear modification of single inclusive heavy flavor mesons. The nuclear shadowing effect in nucleus-nucleus collisions has been included by utilizing the EPS09 parametrization (68). In the study of heavy flavor correlations, the initial angular and transverse momentum distributions of heavy flavor pairs are obtained from Pythia 8 simulation. The evolution of heavy quarks inside the QGP ends when they reach fluid cells with local temperature below (165 MeV). On the freeze-out hypersurface, heavy quarks hadronize into color neutral bound states using our hybrid model of fragmentation and coalescence developed in Ref. (35); (36). In such a hybrid model, high momentum heavy quarks tend to fragment directly into hadrons, but low momentum heavy quarks are more probable to combine with thermal partons from the medium to form heavy flavor mesons. The relative probability between these two processes is determined by the Wigner functions based on an instantaneous coalescence model (69). Within our coalescence model, the momentum distribution of the produced hadron can be calculated once a heavy quark is chosen to recombine with a thermal parton. If the heavy quark is selected to fragment, Pythia 6 (70) is used to simulate the fragmentation process using Peterson fragmentation function. The rescattering of the produced heavy mesons in a hadron gas will not be included in this work, since our previous study (36) has shown that its effect can be reasonably approximated by slightly increasing the momentum space transport coefficient of heavy quarks in the QGP phase.

In our model, the transport coefficient of heavy quarks is determined by comparing to the experimental data on heavy meson nuclear modification at high and then applied to calculations at lower and other heavy flavor observables. In Fig. 1, we provide calculations of the meson in 2.76 TeV central Pb-Pb collisions and 200 GeV Au-Au collisions. Different energy loss mechanisms are utilized and compared to the experimental data. We can see that with the simultaneous inclusion of collisional and radiative processes, the use of provides the best description for the dependence of the meson . If only quasi-elastic scattering or only gluon radiation are included, one can also obtain a reasonable description of the experimental data if the diffusion coefficient is adjusted to or , respectively. The uncertainty introduced by a different velocity dependence of the longitudinal and transverse transport coefficients ( and ) for the collisional energy loss (denoted by the star notation) is investigated as well: the velocity dependence calculated to the leading logarithm is taken from Ref. (29). In Fig. 1 we show that an adjustment of the diffusion constant to is able to reproduce our earlier results obtained with velocity-independent transport coefficients. Fig. 1 implies that the single particle spectrum appears not sufficient to distinguish between the contributions from different energy loss mechanisms. In the next section, we will show that this large uncertainty can be reduced by investigating the correlations of heavy flavor pairs. In this work we perform the calculations for 200 GeV Au-Au collisions at RHIC.

## Iii Angular correlations of heavy flavor pairs

In this section, we study the angular correlations of heavy flavor pairs and their nuclear modifications in relativistic heavy-ion collisions. The dependence of the nuclear modification effect on the centrality cut and the transverse momentum cuts will be investigated in some details.

As a first academic study, we initialize heavy quark pair production using LO pQCD calculation, where and pairs are assumed to be back-to-back and have the same magnitudes of transvese momentum. We calculate the angular correlation functions of pairs after they traverse a realistic QGP medium created in central Au-Au collisions at RHIC. The result is shown in Fig. 2, in which denotes the number of pairs (later for pairs within the selected trigger class) so that our correlation function is normalized to unity. The same values of transport coefficients are used here as in Fig. 1. We can see that while different energy loss mechanisms provide similar amount of suppression of single mesons, they result in very different angular correlation functions of pairs in the final state. Since the medium-induced gluon radiation is mostly dominated by the transverse momentum kicks exerted on the radiated gluons, the effect on the propagation directions of the heavy quarks is very mild. One can see that the angular distribution of the final state pairs still peaks around if only radiative energy loss is included. However, quasi-elastic scatterings are very effective in changing the heavy quark propagation directions, and the final angular distribution is widely spread out. Such result implies that the angular correlations of heavy flavor pairs provide a possibility to distinguish different energy loss mechanisms.

One interesting observation in Fig. 2 is that purely collisional energy loss leads to a peak around in the angular distribution. This indicates that a lot of low momentum pairs move collinearly in the final state even though they are initialized as back-to-back pairs. This is caused by the strong radial flow of the hydrodynamic background which boosts the pairs into its own expanding direction. Such effect was historically referred to as the “partonic wind effect” in Ref. (55); it might lead to some enhancement of quarkonium regeneration in a QGP matter. We verify this effect by solving the Langevin equation in the global (laboratory) frame, i.e., the medium flow effect is shut off manually and only the temperature profile affects the heavy quark evolution. With the flow effect turned off, the peak around disappears and the angular distribution of pairs becomes flat. The flat distribution also means that the final charm quark propagation directions can be entirely randomized by the strong quasielastic scatterings inside the QGP. In Fig. 2, we also observe that collisional energy loss with both velocity independent and velocity dependent and lead to similar angular correlation functions as long as the diffusion coefficients are properly adjusted to describe the observed meson as shown in Fig. 1. In the rest of this work, we will adopt the velocity independent as assumed in our earlier study (36).

The above study of the correlations between pairs is only for academic purposes. For the final state correlations of pairs, we first note that there might be more than one pair produced in one single event of heavy-ion collisions (which depends on the kinematic cuts and collision centrality). To take this effect into account, for each event we pair every meson with all mesons, and analyze all final state pairs. It is clear that there are uncorrelated pairs (background) contributing to the final correlations, especially with low cuts and in more central collisions.

For a more realistic simulation, Pythia 8 is utilized to simulate the initial pair production. We adopt the tunings of the Pythia 8 generator (71) that reproduce the spectrum of single charm hadrons in proton-proton collisions at RHIC as measured by the STAR Collaboration (72) [see Fig. 3], and then use the same parameter setting to obtain the angular correlation functions of pairs in proton-proton collisions as shown in Fig. 3. We can see that the angular correlations are pretty flat if no transverse momentum cut is applied. If one imposes a 2 GeV cut for the trigger or meson, a peak appears around . With higher and higher cuts applied for the trigger and associated heavy mesons, the away side peak becomes more and more prominent. Another interesting observation is that with a high enough cut, a near-side peak begin to appear; this is the contribution from the initial gluon splitting process. The above obtained results serve as a baseline for the following study of correlations in Au-Au collisions at RHIC.

In Fig. 4, we show the angular correlation functions of pairs for 200 GeV central Au-Au collisions at RHIC. Fig. 4 shows the result when a 2 GeV cut is applied to the trigger meson. Similar to pair correlation shown in Fig. 2, the purely radiative energy loss does not strongly affect the angular correlation very much [comparing the dashed line in Fig. 4 and the solid line in Fig. 3]. However, purely collisional energy loss is very effective in suppressing the away-side peak. Elastic scatterings also produce a near side peak around , this effect diminishes when the coupling between the heavy quarks and the collective flow of the QGP is switched off. The realistic scenario in which a combination of both energy loss mechanisms is enabled falls between the above two curves. While these findings may help distinguish between different energy loss mechanisms, note that systematic uncertainties exist for the radiative energy loss at low since currently the detailed balance between gluon emission and absorption is not implemented and we employ a low momentum cut-off below which heavy quarks can only interact via elastic scattering to retain detailed balance for heavy quarks that have thermalized with the bulk QCD medium. Another feature of correlations is the large background contributed from uncorrelated pairs, as compared to angular correlations. In our calculation, there are on average around 10 pairs of ’s produced in each central Au-Au collision; the background contribution diminishes when one moves to peripheral collisions or larger cuts are applied. For example, drops below 1 per trigger event in Fig. 4 when the centrality is greater than about 50%. Fig. 4 shows the same results, but with the application of a 4 GeV cut for the trigger or meson and a 2 GeV cut for the associated meson. One can see that the difference between collisional and radiative energy losses becomes smaller than that in Fig. 4. The difference almost disappears with higher enough cuts applied, see e.g., Fig. 4. As shown by the initial spectrum of charm quarks at RHIC energy [Fig. 3], with the cuts applied in Fig. 4 and 4, on average there should be less than one pair of produced in each event, and thus very few uncorrelated pairs contribute to these two scenarios.

To further quantify the medium modification of the angular correlations of heavy flavor meson pairs, we calculate the variances of the angular distributions for both near side and away side . The result is shown in Fig. 5 where the variances are plotted as functions of the participant numbers in Au-Au collisions at RHIC. The points at correspond to the baseline proton-proton collisions. Note that the variance of is defined as . And for the purpose of a better illustration, we have drawn a horizontal dashed line at which is the variance for a uniform distribution in or in ; a smaller variance than implies a peak around or , and a larger variance than implies a dip on the near or away side.

Figure 5 shows the variances for the near side correlations. When a 2 GeV cut is applied to the trigger (the upper panel), there is a dip structure in the correlations in proton-proton collisions and peripheral Au-Au collisions (therefore, the variance is larger than the value of a uniform distribution). With purely radiative energy loss, such dip gradually approaches a uniform distribution as one moves to more central collisions. However, with purely collisional energy loss, the variance first drops below when is not very large, and then increases and gradually approaches to the uniform distribution value as one moves to more central collisions. The smaller variance than in central and mid-central collisions indicates there is a peak formed around due to the “partonic wind effect” as discussed earlier. The non-monotonic behavior of the variance as a function of originates from the increasing contribution from the uncorrelated pairs in more central collisions which smoothes out the angular distribution and make the variance approach the uniform distribution value.

With the application of larger cuts [the middle and lower panels in Fig. 5], the variance at the near side starts from a value below for proton-proton collisions and decreases with the increase of in Au-Au collisions. This behavior can be understood from the in-medium energy loss of heavy quarks as follows. The final mesons within certain cuts originate from initial quark which have much higher momentum and thus smaller variances for their angular distributions [see Fig. 3]. Since the angular correlations do not change very much at high after in-medium evolution [see e.g., Fig. 4], one obtains more peaked angular distribution for heavy flavor pairs for more central collisions. We further verify this by applying the transverse momentum cuts to the initial heavy quarks instead of the final observed heavy mesons. For such fixed sample of heavy quark/meson pairs, we do find that the variance of the angular distribution increases monotonically towards when increasing . Again, we find that the difference among different energy loss mechanisms is larger for low meson pairs and becomes smaller for larger meson pairs.

The variances for the away-side angular distribution are shown in Fig. 5. With small cut (2 GeV), purely radiative energy loss smoothes the away-side peak in correlation distribution and the variance monotonically approaches towards the uniform distribution value. However, with purely collisional energy loss, the variance first increases above the uniform distribution limit which indicates a depletion of pairs around , and then decreases towards due to the smoothing of the background contribution. With higher cuts, there is the interplay between two combinational effects: the transverse momentum broadening and the energy loss of heavy quarks inside the QGP. The former increases the variance towards the uniform distribution limit and is more efficient for relatively lower , while the latter enhances the contribution from the charm quarks with higher initial which reduces the variance of the final angular distribution. As a consequence, a non-monotonic behavior of the variance is observed as a function of collision centrality. In the next section, we will disentangle these two contributions – momentum broadening and energy loss – to the value of the variance by utilizing the momentum imbalance of pairs. An alternative way to distinguish these two contributions – using the half width at the half maximum – was discussed in Ref. (34).

## Iv Momentum imbalance of pairs

In addition to the nuclear modification of heavy flavor pair angular correlations, another interesting observable for studying the correlations of heavy mesons is their momentum imbalance. For this study, we may define the transverse momentum ratio: for pairs, where and are the transverse momenta of the leading and subleading mesons, whose azimuthal angles are denoted as and . For each event, we select the or meson that has the highest transverse momentum as the leading (trigger) meson. On the back side, we look for the anti-particle ( or meson) with the highest transverse momentum and select it as the subleading (associated) meson. An angular cut is also applied for the away-side subleading meson: ; this should be helpful to reduce the background contribution and select more truly correlated pairs.

In Fig. 6, we show the event distribution (normalized to one or trigger) of the ratio () of pairs for different centrality bins and kinematic cuts. The results are shown for two different cuts: GeV, GeV in Fig. 6(a), and GeV, GeV in Fig. 6(b). Two observations can be found in Fig. 6: as one moves from proton-proton collisions to peripheral and subsequently to central Au-Au collisions, (1) a smaller number of pairs per triggered event can be seen, and (2) the distributions for both sets of selections tend to shift to the region of smaller (i.e., larger momentum imbalance). These are both due to stronger energy loss of heavy quarks inside the dense nuclear medium in more central collisions. As discussed earlier, there should be very few uncorrelated pairs contributing to the correlation functions here with the chosen cuts. In the subfigures of Fig. 6, we also present the ratios between nucleus-nucleus collisions and proton-proton collisions, i.e., as defined by

(4) |

One can see that the values of are above unity at small and below at large . At large , it decreases with the centrality of the collisions.

To have a cleaner description of the transverse momentum imbalance of pairs as a function of collision centrality, we show in Fig. 7 the average value of the ratio (or ) between the subleading and the leading mesons as a function of for the two sets of kinematic cuts, where corresponds to the proton-proton reference. One can see that as one moves from proton-proton collisions to central Au-Au collisions, the transverse momentum imbalance of pairs increases ( value changes about 6%). It might be easier to observe a larger medium modification to the momentum imbalance of pairs at the LHC with even higher cuts, as shown in Ref. (63).

With the knowledge of one may further classify the events within the same centrality bin. It has been shown in Ref. (52) that for -jet larger corresponds to events in which hard partons are produced at the edges of the QGP fireballs while smaller values are contributed by those pairs produced at the centers. In Fig. 8 we implement a similar study to investigate the correlation between the values of pairs and the production positions of the initial pairs. We observe that for heavy flavor pairs, smaller corresponds to events initially produced at the edge of the QGP fireballs in which one heavy quark travels outside the medium without much interaction while its partner traverse the whole QGP fireball and loses significantly amount of energy. For larger values, initial pairs are more likely to spread out smoothly over the QGP. Our study thus confirms the pattern found in (62) (Fig. 9). This allows us to utilize this momentum imbalance to probe different regions of the dense nuclear matter and gain knowledge on the path length dependence of parton energy loss as well. Note that in Fig. 6-8, we present our calculations with both collisional and radiative energy loss switched on. We have verified that since the momentum imbalance is directly related to the total amount of energy loss of heavy quarks, its value does not strongly depend on the detailed energy loss mechanism as long as the transport coefficient is properly adjusted to describe the meson .

The momentum imbalance allows us to refine the previous study on the angular correlation function in Sec. III. In order investigate the competition between the two ingredients that affect the variance of the angular correlation function – momentum broadening and energy loss – as shown in Fig. 5, the values of the variance are now calculated for different bins. Here we utilize the away-side variance obtained with the intermediate triggers as an example, which has the best potential to constrain the different energy loss models. In Fig. 9 we can observe when is small (), it is easier for the large energy loss to overwhelm the effect of the momentum broadening. With the pure radiative energy loss, the away side variance decreases monotonically; to the contrary, the pure collisional energy loss would increase this variance in more central collisions since it smears the direction of heavy quarks more effectively inside the QGP; and the combined energy loss mechanism leads to a relatively flat variance when the participant number is large. However, in higher bins (or with smaller energy loss), the effect of momentum broadening starts to dominate when is large and thus increases the away side variance towards the value of a uniform distribution. And the turning point between the decrease and increase of the variance starts earlier when a higher the cut is applied. Thus, the combination of the momentum imbalance and the angular correlation function may provide us a quantitative tool to study the detailed energy loss dynamics of heavy quarks inside the QGP.

## V Conclusions

In this work, we have studied the correlations of heavy meson pairs in relativistic heavy-ion collisions. The transport of open heavy quarks inside a dense nuclear matter is simulated utilizing our modified Langevin equation that incorporates both quasi-elastic scattering and medium-induced gluon radiation. The spacetime evolution of QGP fireball produced in Au-Au collisions at RHIC is simulated with a (2+1)-dimensional viscous hydrodynamic model. A hybrid model of fragmentation and heavy-light quark coalescence is applied to simulate the hadronization of heavy quarks into their mesonic states.

We have explored two relatively new observables for the nuclear modification of heavy flavors: the angular correlations and the transverse momentum imbalance of heavy meson pairs. Our results have shown that the transverse momentum imbalance of heavy meson pairs provides another measure of the energy loss of heavy quarks inside the QGP. Compared to the nuclear modification of single inclusive meson spectrum, one of the advantages of this observable is that it is not only able to quantify the total energy loss of the heavy quark, but also provides an opportunity to help probe different areas of the QGP fireball. It is found that with the increase of collision centrality, the average ratio between the subleading and the leading heavy flavor mesons decreases, which indicates a larger imbalance caused by the larger energy loss of heavy quarks. One the other hand, the angular correlation of pairs provides a good candidate to characterize the momentum broadening of heavy quarks inside the QGP, and is also sensitive to different energy loss mechanisms. We find that while purely radiative energy loss does not affect very much the angular distribution between pairs, the collisional energy loss is very effective in randomizing the propagation direction of heavy quarks, especially at low to intermediate . To better quantify the medium modification of the heavy flavor pair correlations, we have also calculated the variances for both near-side and away-side angular distributions and investigate them as a function of collision centrality and transverse momentum cuts. By utilizing different bins of the momentum imbalance, the competition between the two effects – momentum broadening and energy loss – on the medium modification of the variance can be clearly seen on the away side when intermediate cuts are implemented.

Our study constitutes an important contribution towards a more quantitative understanding of the energy loss and momentum broadening of heavy quarks inside a hot and dense nuclear matter. If heavy flavor meson pair correlations are measured in the future, one can not only infer how much energy is lost from heavy quarks, but also gain some insights on the detailed interaction mechanism between heavy quarks and the QCD medium; the later seems quite difficult by only studying the nuclear modification of single inclusive heavy flavor meson spectra. Nevertheless, it is worth noticing that in order to draw further quantitative conclusion on the heavy flavor dynamics from these correlation functions, it is crucial to develop a sound method for background subtraction, especially in the low regime and in central collisions.

## Acknowledgments

We are grateful to valuable discussions with X.-N. Wang and M. Nahrgang, and computational resources provided by the Open Science Grid (OSG). This work is funded by the Director, Office of Energy Research, Office of High Energy and Nuclear Physics, Division of Nuclear Physics, of the U.S. Department of Energy under Contract Nos. DE-AC02-05CH11231 and DE-FG02-05ER41367, and within the framework of the JET Collaboration, and by the Natural Science Foundation of China (NSFC) under Grant No. 11375072.

### References

- M. Gyulassy and L. McLerran, Nucl. Phys. A750, 30 (2005), nucl-th/0405013.
- PHENIX, K. Adcox et al., Nucl. Phys. A757, 184 (2005), nucl-ex/0410003.
- BRAHMS, I. Arsene et al., Nucl. Phys. A757, 1 (2005), nucl-ex/0410020.
- B. B. Back et al., Nucl. Phys. A757, 28 (2005), nucl-ex/0410022.
- STAR, J. Adams et al., Nucl. Phys. A757, 102 (2005), nucl-ex/0501009.
- B. Muller and J. L. Nagle, Ann. Rev. Nucl. Part. Sci. 56, 93 (2006), arXiv:nucl-th/0602029.
- B. Muller, J. Schukraft, and B. Wyslouch, Ann. Rev. Nucl. Part. Sci. 62, 361 (2012), arXiv:1202.3233.
- G.-Y. Qin, Int. J. Mod. Phys. E24, 1530001 (2015), arXiv:1502.02554.
- PHENIX Collaboration, A. Adare et al., Phys. Rev. Lett. 107, 252301 (2011), arXiv:1105.3928.
- STAR Collaboration, L. Adamczyk et al., Phys. Rev. C88, 014904 (2013), arXiv:1301.2187.
- ATLAS Collaboration, G. Aad et al., Phys. Rev. C86, 014907 (2012), arXiv:1203.3087.
- U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826.
- C. Gale, S. Jeon, and B. Schenke, Int. J. Mod. Phys. A28, 1340011 (2013), arXiv:1301.5893.
- H. Song, Pramana 84, 703 (2015), arXiv:1401.0079.
- P. Romatschke, Int. J. Mod. Phys. E19, 1 (2010), arXiv:0902.3663.
- P. Huovinen, Int. J. Mod. Phys. E22, 1330029 (2013), arXiv:1311.1849.
- PHENIX, K. Adcox et al., Phys. Rev. Lett. 88, 022301 (2002), nucl-ex/0109003.
- STAR, C. Adler et al., Phys. Rev. Lett. 89, 202301 (2002), nucl-ex/0206011.
- ALICE Collaboration, K. Aamodt et al., Phys. Lett. B696, 30 (2011), arXiv:1012.1004.
- S. A. Bass et al., Phys. Rev. C79, 024901 (2009), arXiv:0808.0908.
- JET, K. M. Burke et al., Phys. Rev. C90, 014909 (2014), arXiv:1312.5003.
- Y. L. Dokshitzer and D. E. Kharzeev, Phys. Lett. B519, 199 (2001), arXiv:hep-ph/0106202.
- PHENIX Collaboration, A. Adare et al., Phys. Rev. C84, 044905 (2011), arXiv:1005.1627.
- PHENIX, A. Adare et al., Phys. Rev. C91, 044907 (2015), arXiv:1405.3301.
- STAR collaboration, D. Tlusty, Nucl. Phys. A904-905, 639c (2013), arXiv:1211.5995.
- STAR Collaboration, L. Adamczyk et al., Phys. Rev. Lett. 113, 142301 (2014), arXiv:1404.6185.
- ALICE Collaboration, A. Grelli, Nucl. Phys. A904-905, 635c (2013), arXiv:1210.7332.
- ALICE Collaboration, B. Abelev et al., Phys. Rev. Lett. 111, 102301 (2013), arXiv:1305.2707.
- G. D. Moore and D. Teaney, Phys. Rev. C71, 064904 (2005), hep-ph/0412346.
- P. Gossiaux, J. Aichelin, T. Gousset, and V. Guiho, J. Phys. G37, 094019 (2010), arXiv:1001.4166.
- M. He, R. J. Fries, and R. Rapp, Phys. Rev. C86, 014903 (2012), arXiv:1106.6006.
- C. Young, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C86, 034905 (2012), arXiv:1111.0647.
- J. Uphoff, O. Fochler, Z. Xu, and C. Greiner, Phys. Lett. B717, 430 (2012), arXiv:1205.4945.
- M. Nahrgang, J. Aichelin, P. B. Gossiaux, and K. Werner, Phys. Rev. C90, 024907 (2014), arXiv:1305.3823.
- S. Cao, G.-Y. Qin, and S. A. Bass, Phys. Rev. C88, 044907 (2013), arXiv:1308.0617.
- S. Cao, G.-Y. Qin, and S. A. Bass, Phys. Rev. C92, 024907 (2015), arXiv:1505.01413.
- S. Cao and S. A. Bass, Phys. Rev. C84, 064902 (2011), arXiv:1108.5101.
- S. Cao, G.-Y. Qin, and S. A. Bass, J. Phys. G40, 085103 (2013), arXiv:1205.2396.
- S. Cao, G.-Y. Qin, S. A. Bass, and B. Muller, Nucl. Phys. A904-905, 653c (2013), arXiv:1209.5410.
- S. A. Bass et al., Prog. Part. Nucl. Phys. 41, 225 (1998), nucl-th/9803035.
- STAR, C. Adler et al., Phys. Rev. Lett. 90, 082302 (2003), nucl-ex/0210033.
- STAR, J. Adams et al., Phys. Rev. Lett. 97, 162301 (2006), arXiv:nucl-ex/0604018.
- PHENIX, A. Adare et al., Phys. Rev. C80, 024908 (2009), arXiv:0903.3399.
- STAR, B. I. Abelev et al., Phys. Rev. C82, 034909 (2010), arXiv:0912.1871.
- H. Zhang, J. F. Owens, E. Wang, and X.-N. Wang, Phys. Rev. Lett. 98, 212301 (2007), arXiv:nucl-th/0701045.
- G.-Y. Qin, J. Ruppert, C. Gale, S. Jeon, and G. D. Moore, Phys. Rev. C80, 054909 (2009), arXiv:0906.3280.
- ATLAS, G. Aad et al., Phys. Rev. Lett. 105, 252303 (2010), arXiv:1011.6182.
- CMS, S. Chatrchyan et al., Phys. Rev. C84, 024906 (2011), arXiv:1102.1957.
- ATLAS, G. Aad et al., Phys. Rev. Lett. 114, 072302 (2015), arXiv:1411.2357.
- G.-Y. Qin and B. Muller, Phys. Rev. Lett. 106, 162302 (2011), arXiv:1012.5280.
- W. Dai, I. Vitev, and B.-W. Zhang, Phys. Rev. Lett. 110, 142001 (2013), arXiv:1207.5177.
- G.-Y. Qin, Eur. Phys. J. C74, 2959 (2014), arXiv:1210.6610.
- X.-N. Wang and Y. Zhu, Phys. Rev. Lett. 111, 062301 (2013), arXiv:1302.5874.
- X. Zhu et al., Phys. Lett. B647, 366 (2007), arXiv:hep-ph/0604178.
- X. Zhu, N. Xu, and P. Zhuang, Phys. Rev. Lett. 100, 152301 (2008), arXiv:0709.0157.
- P. Gossiaux, V. Guiho, and J. Aichelin, J. Phys. G32, S359 (2006).
- Y. Akamatsu, T. Hatsuda, and T. Hirano, Phys. Rev. C80, 031901 (2009), arXiv:0907.0588.
- M. Younus and D. K. Srivastava, J. Phys. G40, 065004 (2013), arXiv:1301.5762.
- S. Cao, G.-Y. Qin, S. A. Bass, and B. Muller, J. Phys. Conf. Ser. 446, 012035 (2013).
- S. Cao, G.-Y. Qin, and S. A. Bass, Nucl.Phys. A932, 38 (2014), arXiv:1404.1081.
- S. Cao, G.-Y. Qin, and S. A. Bass, Nucl. Phys. A931, 569 (2014), arXiv:1408.0503.
- P. Gossiaux, R. Bierkandt, and J. Aichelin, Phys. Rev. C79, 044906 (2009), arXiv:0901.0946.
- J. Uphoff, F. Senzel, Z. Xu, and C. Greiner, Phys. Rev. C89, 064906 (2014), arXiv:1310.1340.
- X.-F. Guo and X.-N. Wang, Phys. Rev. Lett. 85, 3591 (2000), arXiv:hep-ph/0005044.
- A. Majumder, Phys. Rev. D85, 014023 (2012), arXiv:0912.2987.
- B.-W. Zhang, E. Wang, and X.-N. Wang, Phys. Rev. Lett. 93, 072301 (2004), arXiv:nucl-th/0309040.
- Z. Qiu, C. Shen, and U. Heinz, Phys. Lett. B707, 151 (2012), arXiv:1110.3033.
- K. Eskola, H. Paukkunen, and C. Salgado, JHEP 0904, 065 (2009), arXiv:0902.4154.
- Y. Oh, C. M. Ko, S. H. Lee, and S. Yasui, Phys. Rev. C79, 044905 (2009), arXiv:0901.1382.
- T. Sjostrand, S. Mrenna, and P. Z. Skands, JHEP 0605, 026 (2006), arXiv:hep-ph/0603175.
- S. Shi, X. Dong, and M. Mustafa, (2015), arXiv:1507.00614.
- STAR Collaboration, L. Adamczyk et al., Phys. Rev. D86, 072013 (2012), arXiv:1204.4244.