DAB-MOD sensitivity study of heavy flavor $R_{AA}$ and azimuthal anisotropies based on beam energy, initial conditions, hadronization, and suppression mechanisms

TL;DR

This study uses the DAB-MOD simulation to analyze how initial conditions, hadronization, and suppression mechanisms affect heavy flavor meson suppression and anisotropies in heavy-ion collisions, comparing models and experimental data.

Contribution

It provides a systematic comparison of transport models and initial condition effects on heavy flavor observables using state-of-the-art simulations.

Findings

Different transport models yield consistent $R_{AA}$ and flow results.

Initial conditions significantly influence heavy flavor anisotropies.

Heavy-light quark coalescence impacts the coupling of heavy quarks to the medium.

Abstract

Heavy flavor probes provide important information about the in-medium properties of the quark-gluon plasma produced in heavy-ion collisions. In this work, we investigate the effects of 2D+1 event-by-event fluctuating hydrodynamic backgrounds on the nuclear suppression factor and momentum anisotropies of heavy flavor mesons and non-photonic electrons. Using the state-of-the-art D and B mesons modular simulation code (called "DAB-MOD"), we perform a systematic comparison of different transport equations in the same background, including a few energy loss models --- with and without energy loss fluctuations --- and a relativistic Langevin model with different drag parametrizations. We present the resulting D and B mesons $R_{AA}$, $v_2$, $v_3$, and $v_4$ as well as multi-particle cumulants, in AuAu collisions at $\sqrt{s_\text{NN}}=200\, \text{GeV}$ and PbPb collisions at $\sqrt{s_\text{NN}}=2.76 \, \text{TeV}$ and $\sqrt{s_\text{NN}}=5.02 \, \text{TeV}$, and compare them to the available experimental data. The $v_2\{4\}/v_2\{2\}$ ratio, which is known to be a powerful probe of the initial conditions and flow fluctuations in the soft sector, is also studied in the context of heavy flavor. We also investigate the correlations between the transverse anisotropies of heavy mesons and all charged particles to better understand how heavy quarks couple to the hydrodynamically expanding quark-gluon plasma. We study the influence that different initial conditions and the implementation of heavy-light quark coalescence has on our results.

Figures

Tables

Table 1. Table 1: Values of the coupling factors for charm quarks determined for each transport model, collision energy, and decoupling temperature. These values are obtained using MCKLN initial conditions.

Coupling factors for charm quarks at $T_{\text{d}}=120$ $\backslash$ 160 MeV	RHIC AuAu $\sqrt{s_{\text{NN}}}=200\mathrm{GeV}$	LHC PbPb $\sqrt{s_{\text{NN}}}=2.76\mathrm{TeV}$	LHC PbPb $\sqrt{s_{\text{NN}}}=5.02\mathrm{TeV}$
$\alpha$ without fluctuations	0.393 $\backslash$ 0.623	1.0 $\backslash$ 1.624	0.708 $\backslash$ 1.011
$\alpha$ with uniform fluct.	0.649 $\backslash$ none	1.7 $\backslash$ none	0.993 $\backslash$ none
$\alpha$ with linear fluct.	0.77 $\backslash$ none	2.024 $\backslash$ none	1.130 $\backslash$ none
$\alpha$ with gaussian fluct.	0.43 $\backslash$ none	1.1 $\backslash$ none	0.751 $\backslash$ none
$\xi$	11.57 $\backslash$ 15.16	30.28 $\backslash$ 40.05	14.76 $\backslash$ 17.16
	0.48 $\backslash$ 0.34	0.227 $\backslash$ 0.169	0.5 $\backslash$ 0.41
	0.639 $\backslash$ 0.921	1.039 $\backslash$ 1.577	0.622 $\backslash$ 0.828

Table 2. Table 2: Values of the coupling factors for bottom quarks determined for each transport model, collision energy, and decoupling temperature. These values are obtained using MCKLN initial conditions.

Coupling factors for bottom quarks at $T_{\text{d}}=120$ $\backslash$ 160 MeV	RHIC AuAu $\sqrt{s_{\text{NN}}}=200\mathrm{GeV}$	LHC PbPb $\sqrt{s_{\text{NN}}}=2.76\mathrm{TeV}$	LHC PbPb $\sqrt{s_{\text{NN}}}=5.02\mathrm{TeV}$
$\alpha$ without fluctuations	0.264 $\backslash$ 0.4	0.72 $\backslash$ 1.12	0.667 $\backslash$ 0.823
$\alpha$ with uniform fluct.	0.316 $\backslash$ none	0.857 $\backslash$ none	0.824 $\backslash$ none
$\alpha$ with linear fluct.	0.339 $\backslash$ none	0.921 $\backslash$ none	0.913 $\backslash$ none
$\alpha$ with gaussian fluct.	0.265 $\backslash$ none	0.76 $\backslash$ none	0.624 $\backslash$ none
$\xi$	7.6 $\backslash$ 10	21.52 $\backslash$ 27.06	none $\backslash$ none
	0.648 $\backslash$ 0.486	0.32 $\backslash$ 0.226	0.516 $\backslash$ 0.411
	0.606 $\backslash$ 0.808	3.21 $\backslash$ 2.26	0.681 $\backslash$ 0.884

Table 3. Table 3: The approximate low p T subscript 𝑝 T p_{\text{T}} branching ratios or prompt hadronic ratios obtained with different configurations of the coalescence model (using a typical | 𝐮 | ≈ 0.6 𝐮 0.6 |{\bf u}|\approx 0.6 ) compared to the low p T subscript 𝑝 T p_{\text{T}} experimental data in AA collisions (or pp collisions if specified).

Approximate values at low $p_{\text{T}}$	Experimental Acharya et al. (2018d); Radhakrishnan (2018); Vermunt (2019); Acharya et al. (2018e); Gladilin (2015); Cacciari and Nason (2003)	“basic” coalescence	with extra thermal factors
direct $\text{c}\to {\text{D}}^0$	in pp: 0.17	0.06	0.10
prompt $\text{c}\to {\text{D}}^0$	in pp: 0.55	0.36	0.32
${\text{D}}^{+}/{\text{D}}^0$	0.47	0.32	0.45
${\text{D}}^{\ast +}/{\text{D}}^0$	0.45	0.5	0.4
${\text{D}}_s^{+}/{\text{D}}^0$	0.35	0.31	0.34
${\text{D}}_s^{+}/{\text{D}}^{+}$	0.75	0.99	0.76
${\Lambda }_c^{+}/{\text{D}}^0$	@ RHIC: 1.5 - @ LHC: 0.7	0.74	0.94

Equations

$$\frac{\text{d}E}{\text{d}x}(T,v_{\text{Q}})=-f(T,v_{\text{Q}})\zeta {\Gamma }_{\text{flow}},$$

$$\varGamma_{\text{flow}}=\gamma\big{[}1-v_{\text{flow}}\cos(\varphi_{\text{Q}}-\varphi_{\text{flow}})\big{]},$$

$${\Gamma }_{\text{flow}}^{\mathrm{exact}}=\gamma \sqrt{1-2\frac{v_{\text{flow}}}{v_{\text{Q}}}\cos ({\phi }_{\text{Q}}-{\phi }_{\text{flow}})+\frac{v_{\text{flow}}^2}{v_{\text{Q}}^2}-v_{\text{flow}}^2{\sin }^2({\phi }_{\text{Q}}-{\phi }_{\text{flow}})},$$

$$\text{d}{x}_i=\frac{p_i}{E}\text{d}t,$$

$$\text{d}{p}_i=-\Gamma (\mathbf{p})p_i\text{d}t+\sqrt{\text{d}t}\sqrt{\kappa }{\rho }_i,$$

$$\kappa =2E\Gamma T=2T^2/D,$$

$$p_i^{\prime }=p_i^{\prime }-\Gamma ({\mathbf{p}}^{\prime })p_i^{\prime }\Delta t^{\prime }+\sqrt{\Delta t^{\prime }}\sqrt{\kappa }{\rho }_i,$$

$$\Delta t^{\prime }=\Delta t\ast \frac{p^{\mu }{v}_{\mu }^{\mathrm{flow}}}{E}=\Delta t\ast \gamma \left[1-\frac{\Vert \mathbf{p}\Vert \Vert {\mathbf{v}}_{\mathrm{flow}}\Vert }{E}\cos ({\phi }_{\mathrm{Q}}-{\phi }_{\mathrm{flow}})\right],$$

$$x_i(t+\Delta t)=x_i(t)+\frac{p_i}{E}\Delta t.$$

$$D_{M\&T}=k_{M\&T}/(2\pi T),$$

$${\Gamma }_{G\&A}=k_{G\&A}\ast A_{G\&A}.$$

$$f(z)\propto \frac{1}{z(1-1/z-{ϵ}_{\text{Q}}/(1-z){)}^2}$$

$$R_{\text{AA}}(p_{\text{T}},\varphi )=\frac{1}{\mathcal{N}}\frac{\text{d}{N}_{\text{AA}}/\text{d}{p}_{\text{T}}\text{d}\varphi }{\text{d}{N}_{pp}/\text{d}{p}_{\text{T}}},$$

$$V_n=v_n{\text{e}}^{\text{i}n{\psi }_n}.$$

$$V_n(p_{\text{T}})=v_n(p_{\text{T}}){\text{e}}^{\text{i}n{\psi }_n(p_{\text{T}})}.$$

$$\operatorname{Re}\{V_{n}V_{n}^{*}(p_{\text{T}})\}=v_{n}v_{n}(p_{\text{T}})\,\cos\big{[}n\big{(}\psi_{n}-\psi_{n}(p_{\text{T}})\big{)}\big{]}.$$

$$v_n\{2\}(p_{\text{T}})$$

$$v_n\{4\}(p_{\text{T}})$$

$$d_n\{2\}(p_{\text{T}})$$

$$c_{n,j}\{2\}$$

$$d_n\{4\}(p_{\text{T}})$$

$$\displaystyle=\big{\langle}2\,c_{n,j}\{2\}\,d_{n,j}\{2\}(p_{\text{T}})-\langle v_{n}^{2}V_{n}V_{n}^{*}(p_{\text{T}})\rangle_{j}\big{\rangle}$$

$$-c_{n,j}\{4\}$$

$$\frac{v_n\{4\}(p_{\text{T}})}{v_n\{2\}(p_{\text{T}})}=\frac{v_n\{4\}}{v_n\{2\}}\left[1+{\left(\frac{v_n\{2\}}{v_n\{4\}}\right)}^4\left(\underset{\text{soft fluctuations}}{\underbrace{\frac{⟨v_n^4⟩}{⟨v_n^2{⟩}^2}}}-\underset{\text{hard fluctuations}}{\underbrace{\frac{⟨v_n^2{V}_n{V}_n^{\ast }(p_{\text{T}})⟩}{⟨v_n^2⟩⟨V_n{V}_n^{\ast }(p_{\text{T}})⟩}}}\right)\right],$$

$$\frac{⟨v_n^2{V}_n{V}_n^{\ast }(p_{\text{T}})⟩}{⟨v_n^2⟩⟨V_n{V}_n^{\ast }(p_{\text{T}})⟩}$$

$$\frac{v_n\{4\}(p_{\text{T}})}{v_n\{2\}(p_{\text{T}})}$$

$$\frac{⟨v_n^2{V}_n{V}_n^{\ast }(p_{\text{T}})⟩}{⟨v_n^2⟩⟨V_n{V}_n^{\ast }(p_{\text{T}})⟩}$$

$${\mathrm{P}}_{\mathrm{coal}}[q,Q\to M]({\mathbf{p}}_Q,\mathbf{u})=N\int {\text{d}}^3{\mathbf{p}}_q{f}_M({\mathbf{p}}_q,{\mathbf{p}}_Q)n_q({\mathbf{p}}_q,\mathbf{u},T_{\text{d}})\delta ({\mathbf{p}}_M-{\mathbf{p}}_q-{\mathbf{p}}_Q),$$

$$f_M({\mathbf{p}}_q,{\mathbf{p}}_Q)=g_M{h}_M\frac{(2\sqrt{\pi }\sigma {)}^3}{(2\pi {)}^3}{e}^{-{\mathbf{p}}_{\mathrm{rel}}^2{\sigma }^2},$$

$${\mathbf{p}}_{\mathrm{rel}}=\frac{m_q{\mathbf{p}}_Q^{\prime }-m_Q{\mathbf{p}}_q^{\prime }}{m_q+m_Q},$$

$$n_q({\mathbf{p}}_q,\mathbf{u},T_{\text{d}})=\frac{g_q}{e^{\sqrt{{{\mathbf{p}}_q^{\mathrm{cell}}}^2+m_q^2}/T_{\text{d}}}+1}=\frac{g_q}{e^{p_q.u/T_{\text{d}}}+1},$$

$${\mathrm{P}}_{\mathrm{coal}}[q_1,q_2,Q\to B]({\mathbf{p}}_Q,\mathbf{u})=N\int {\text{d}}^3{\mathbf{p}}_{q_1}{\text{d}}^3{\mathbf{p}}_{q_2}{f}_B({\mathbf{p}}_{q_1},{\mathbf{p}}_{q_2},{\mathbf{p}}_Q)n_{q_1}({\mathbf{p}}_{q_1},\mathbf{u},T_{\text{d}})\times$$

$$\times n_{q_2}({\mathbf{p}}_{q_2},\mathbf{u},T_{\text{d}})\delta ({\mathbf{p}}_B-{\mathbf{p}}_{q_1}-{\mathbf{p}}_{q_2}-{\mathbf{p}}_Q),$$

$$f_B({\mathbf{p}}_{q_1},{\mathbf{p}}_{q_2},{\mathbf{p}}_Q)=g_B{h}_B\frac{(2\sqrt{\pi }{)}^6({\sigma }_1{\sigma }_2{)}^3}{(2\pi {)}^6}{e}^{-{\mathbf{p}}_{rel1}^2{\sigma }_1^2-{\mathbf{p}}_{rel2}^2{\sigma }_2^2},$$

$${\mathbf{p}}_{rel1}=\frac{m_{q_2}{\mathbf{p}}_{q_1}^{\prime }-m_{q_1}{\mathbf{p}}_{q_2}^{\prime }}{m_{q_1}+m_{q_2}}$$

Full text

I Introduction

The quark-gluon plasma (QGP) — a deconfined state of mater composed of strongly interacting quarks and gluons — has been successfully reproduced in the laboratory both at the Large Hadron Collider (LHC) and at the Relativistic Heavy Ion Collider (RHIC). While it remains unknown how the strongly fluctuating color fields present in the initial stages of heavy-ion collisions Schenke et al. (2012) may later evolve to a state that displays hydrodynamic behavior (for a review, see Romatschke and Romatschke (2019)), it is well established by now that the soft, low transverse momentum $p_{\text{T}}$ sector of the QGP determined by the particle spectra and flow harmonics of light hadrons can be generally well-described by event-by-event relativistic viscous hydrodynamic simulations Niemi et al. (2016); Noronha-Hostler et al. (2016a); Eskola et al. (2018); Giacalone et al. (2018); Gale et al. (2013); Bernhard et al. (2016). Some of the properties of the early stages of heavy-ion collisions can be investigated using hard probes, such as light flavor jets and heavy flavor hadrons, which are created at very early times via hard scattering processes and then propagate through the evolving medium. Because of the large mass of the heavy quarks, many times larger than the cross-over temperature from the QGP to the hadron resonance gas Borsanyi et al. (2010), they are not likely to be created within the hydrodynamic evolution of the QGP and, thus, heavy flavor content is preserved from the initial stages well into hadronization Moore and Teaney (2005). Therefore, heavy flavor observables can provide key information about the early time dynamics of heavy-ion collisions.

Furthermore, heavy mesons in the intermediate/high $p_{\text{T}}$ region lose energy mostly from radiative processes whereas the low $p_{\text{T}}$ regime is dominated by collisional processes that may be described via Langevin-like equations Andronic et al. (2016); Moore and Teaney (2005). The two main experimental measurements involving open heavy flavor D and B mesons that we are interested in are the nuclear modification factor, $R_{\text{AA}}$, and the azimuthal anisotropies, $v_{n}(p_{\text{T}})$. $R_{\text{AA}}$ mainly encodes how much energy is lost in the medium compared to elementary pp collisions and, since the farther the heavy quark travels through the QGP medium the more energy it loses, one expects that heavy quarks lose more energy in media formed by heavy-ion collisions at high energies. Additionally, because each event has on average certain geometrical features (an almond-like shape especially for mid-central collisions), one expects that heavy quarks lose more energy along the long axis of the almond shape (in-plane) vs. the short axis (out-of-plane), which leads to a final elliptical azimuthal anisotropy $v_{2}$ defined by the second Fourier coefficient of the particle spectra  Wang (2001); Gyulassy et al. (2001).

However, quantum mechanical fluctuations involving the position of the incoming nucleons, and likely also the fluctuations of quark and gluon fields, lead to other initial geometrical patterns such as, e.g., a triangular geometry, producing $v_{3}$ and other even higher order harmonics Nahrgang et al. (2015); Prado et al. (2017), which have been measured in the heavy flavor sector at RHIC Lomnitz (2016) and the LHC Aaboud et al. (2018). Furthermore, we note that the overall flow of the medium has been shown to influence energy loss mechanisms Brewer et al. (2018) and hard probes can provide important information about the properties of the initial state Cao et al. (2015a); Andres et al. (2019). More specifically, it was recently shown in Andres et al. (2019) that the time $\tau_{0}$ after which hard probes begin to lose energy appears to be larger than was initially expected. Further proof of the influence of event-by-event fluctuations of the initial conditions on the heavy flavor sector can be seen in event-shape engineering calculations Prado et al. (2017) and measurements Acharya et al. (2018a) of heavy flavor flow harmonics.

Theoretical comparisons to the measured nuclear modification factor demonstrate that energy loss models alone significantly underpredict $R_{\text{AA}}$ in the low $p_{\text{T}}$ sector (azimuthal anisotropies experience the same issue) though they can describe the high $p_{\text{T}}$ regime quite well Prado et al. (2017). On the other hand, Langevin or Boltzmann based models perform quite well in the low $p_{\text{T}}$ sector but the addition of coalescence is needed to reproduce experimental data Greco et al. (2003a, b); Fries et al. (2003a); Greco et al. (2004); Oh et al. (2009); Nahrgang et al. (2015); Cao et al. (2016). In this manner, the $p_{\text{T}}$ dependence of $R_{\text{AA}}$ can be a useful tool to understand different regimes of heavy flavor energy loss. The heavy flavor $v_{n}$’s depend on the assumptions regarding the heavy flavor model as well but they are also sensitive to the choice of initial conditions, which are strongly related to final result for the flow harmonics (proven by the strong correlation quantified by the Pearson coefficient between the initial eccentricity vector and the final $V_{n}$ vector). There is likely also a sensitivity to the choice of shear viscosity to entropy density ratio, $\eta/s$, see Betz et al. (2017) and Esha et al. (2017). Thus, special care must be taken to first fully constrain the initial conditions and medium viscosity using the soft sector before calculating heavy flavor flow harmonic observables. This has been done in recent years in Prado et al. (2017) and also by DUKE/SUBATECH Cao et al. (2015b); Nahrgang et al. (2015), which has led to a better theoretical understanding of heavy flavor momentum anisotropies.

It is important to note that event-by-event initial state fluctuations are not the only source of fluctuations that can affect heavy flavor $R_{\text{AA}}$ and $v_{n}$. In fact, energy loss fluctuations can also occur in a given event and that has already been shown to affect light flavor high $p_{\text{T}}$ flow harmonics Zapp (2014); Betz and Gyulassy (2014); Betz et al. (2017). Thus, a systematic study of the effect of energy loss fluctuations should also be performed in the heavy flavor sector. Furthermore, the advent of multi-particle flow cumulants involving hard probes (where the hard probe particle of interest is correlated with other reference soft particles Sirunyan et al. (2018a); Betz et al. (2017)) provides a unique opportunity to study heavy flavor from an entirely new angle.

In order to confront this complicated emerging picture of the heavy flavor sector where experimental observables are influenced by multiple competing factors, in this paper we systematically study the effects of $2\text{D}+1$ (i.e., longitudinally boost invariant) event-by-event fluctuating hydrodynamic backgrounds on the nuclear suppression factor and momentum anisotropies of heavy flavor mesons and non-photonic electrons. The bulk dynamics of the medium is simulated using the event-by-event relativistic viscous hydrodynamic model, v-USPhydro Noronha-Hostler et al. (2013, 2014), coupled to either MCKLN initial conditions Drescher et al. (2006); Drescher and Nara (2007a, b) or Trento initial conditions (tuned to IP-Glasma Gale et al. (2013)) Moreland et al. (2015), with hydrodynamic parameters constrained to describe the relevant properties of the soft sector. The heavy flavor sector is described using the state-of-the-art D and B mesons modular simulation code (called “DAB-MOD”)111A preliminary version of DAB-MOD, which did not yet include coalescence and Langevin dynamics, was used in Prado et al. (2017)., with which we perform a systematic study of different transport equations, including a few energy loss models (with and without energy loss fluctuations) and a relativistic Langevin model with different drag parametrizations, and investigate their effect on $R_{\text{AA}}$ and $v_{n}$. We present calculations of D and B meson $R_{\text{AA}}$, $v_{2}$, $v_{3}$, and $v_{4}$ as well as multi-particle flow cumulants in AuAu collisions at $\sqrt{s_{\text{NN}}}=$200\text{,}\mathrm{GeV} and PbPb collisions at $\sqrt{s_{\text{NN}}}=$2.76\text{\,}\mathrm{TeV} and $\sqrt{s_{\text{NN}}}=$5.02\text{,}\mathrm{TeV}$$, and compare them to the available experimental data. Heavy flavor multi-particle flow cumulants are investigated via the $v_{2}\{4\}/v_{2}\{2\}$ ratio, which plays an important role in determining the initial conditions and flow fluctuations in the soft sector. We also study the correlations between the momentum anisotropies of heavy mesons and the flow of all charged particles to better understand how heavy quarks couple to the expanding medium.

This paper is organized as follows. In the next section we explain the different assumptions and details regarding the modeling of heavy flavor dynamics included in DAB-MOD. Sec. III investigates how the transport model assumptions and variations in the decoupling temperature affect the observables considering the case of MCKLN initial conditions without coalescence. We study how the results change when one replaces MCKLN by Trento initial conditions in Sec. IV. Sec. V describes how coalescence is now implemented in DAB-MOD and how it affects our results. Our final remarks are presented in Sec. VI.

II Details of the DAB-MOD simulation

To study the propagation of charm and bottom quarks inside the medium created in high energy heavy-ion collisions, we developed a modular Monte Carlo simulation, DAB-MOD, that allows for a variety of backgrounds and transport models to be implemented.

II.1 The background hydrodynamically expanding medium

Heavy quark transport models require background medium profiles that provide the temperature and flow velocity of the fluid cells at each time step along the heavy quark trajectories. In this work, we use either the event-by-event Monte Carlo Kharzeev-Levin-Nardi (MCKLN ) initial conditions — an implementation of a Color Glass Condensate $k_{\text{T}}$-factorization model Drescher et al. (2006); Drescher and Nara (2007a, b) — or the event-by-event Trento initial conditions tuned to IP-Glasma Moreland et al. (2015).

These initial profiles are evolved from an initial time $\tau_{0}=$0.6\text{,}\mathrm{fm} using the v-USPhydro code Noronha-Hostler *et al.* ([2013](#bib.bib34), [2014](#bib.bib35), [2016b](#bib.bib46)), a $2\text{D}+1$ relativistic viscous hydrodynamical model (with a Cooper-Frye freeze-out prescription) that has passed the standard accuracy tests of the field Marrochio *et al.* ([2015](#bib.bib47)). The viscous hydrodynamic evolution is encoded in a shear viscosity to entropy density ratio that we set for MCKLN (Trento) initial conditions to be $\eta/s=0.08$ ($\eta/s=0.05$) for AuAu collisions at $\sqrt{s_{\text{NN}}}=$200\text{\,}\mathrm{GeV} and to $\eta/s=0.11$ ($\eta/s=0.05$) and $\eta/s=0.05$ ($\eta/s=0.047$) for PbPb collisions at $\sqrt{s_{\text{NN}}}=$2.76\text{,}\mathrm{TeV} and $\sqrt{s_{\text{NN}}}=$5.02\text{\,}\mathrm{TeV} collisions, respectively, following Refs. Noronha-Hostler et al. (2016b, c). This model describes experimental data in the soft sector reasonably well and, thus, all the hydrodynamic parameters that involve the medium description are fixed in the present study and are not seen as free parameters of our heavy quark analysis. To obtain sufficient statistics for the heavy quark observables, we use $\sim$ 1000–2000 hydrodynamic events per 10% centrality range. Finally, we neglect the possible effects of heavy quark energy loss on the evolution of the medium Andrade et al. (2014).

In Fig. 1 we show the results for the $p_{\text{T}}$-integrated two particle flow cumulant, $v_{n}\{2\}$, computed using either MCKLN or Trento initial conditions coupled to v-USPhydro for PbPb 5.02 TeV, PbPb 2.76 TeV, and AuAu 200 GeV collisions. We find that MCKLN initial conditions work best at RHIC energies but they can work reasonably well even at LHC if one is interested only in two particle correlations. Additionally, at RHIC MCKLN leads to a slightly smaller $v_{3}\{2\}$ compared to Trento, which in turn is slightly closer to the experimental data from STAR Adams et al. (2005). This small difference may be alleviated with a temperature dependent $\eta/s$. Considering the results for PbPb 2.76 TeV collisions compared to ATLAS data, we find a reasonably good agreement to experimental data for both MCKLN and Trento initial conditions. However, at the top LHC energy more deviations from experimental data are seen for MCKLN initial conditions while Trento’s description improves. Fig. 1 shows for the first time a comparison between these results and LHC run 2 PbPb data from ALICE at 5.02 TeV Acharya et al. (2018b, c); Adam et al. (2016). Trento initial conditions provide the best fit to ALICE/ATLAS data, although, again, MCKLN is relatively near the data for $v_{2}\{2\}$ and $v_{3}\{2\}$.

The main caveat is that 4-particle correlations and, more specifically, the ratio between $v_{n}\{4\}$ and $v_{n}\{2\}$ which is known to be a good constraint on initial conditions Giacalone et al. (2017) (and is nearly medium independent), paint a different story. In Fig. 2 the ratio $v_{n}\{4\}/v_{n}\{2\}$ for $n=2,3$ is plotted for different beam energies across centrality for Trento vs MCKLN initial conditions Alba et al. (2018). The ratio $v_{2}\{4\}/v_{2}\{2\}$ approaches 1 when there are less $v_{n}$ fluctuations (i.e. a narrower distribution) and it decreases with increasing $v_{n}$ fluctuations. Thus, we find that MCKLN generally predicts less fluctuations than Trento tuned to IP-Glasma. At RHIC this appears to be a good predictor for $v_{2}\{4\}/v_{2}\{2\}$ and MCKLN captures the fluctuations data from STAR well, whereas Trento predicts too much fluctuations. In contrast, at LHC energies we find that MCKLN predicts too few fluctuations whereas Trento provides a nearly perfect fit to experimental data222Note that we are avoiding issues found in ultracentral collisions since we do not expect heavy flavor observables to have sufficient statistics to be considered in that case.. Thus, we conclude that MCKLN works best at RHIC and reasonably well at LHC run 1 but it misses LHC run 2 data whereas Trento manages to do the best job on a global analysis level. However, we note that its predictions for the $v_{2}\{4\}/v_{2}\{2\}$ fluctuations at RHIC is on the low side.

Finally, we compare our results for $v_{3}\{4\}/v_{3}\{2\}$ as well. Even with $\sim 35,000+$ events, statistical error bars (computed using jackknife resampling) are an issue on our side. That being said, we find that both MCKLN and Trento are within error bars of each other. It does appear that MCKLN may have less $v_{3}$ fluctuations compared to Trento but it is difficult to say from hydrodynamics at this level. In Giacalone et al. (2017) it was found that $\varepsilon_{3}\{4\}/\varepsilon_{3}\{2\}$ was consistently larger for MCKLN as well.

II.2 Heavy quark initial conditions

Because of their large masses, heavy quarks are produced at the very beginning of the collisions in parton scatterings that can be described by perturbative QCD (pQCD). Neglecting the possible effects of shadowing at mid-rapidities, which are observed to be non negligible for $p_{\text{T}}<10$ GeV Jena (2014); Song et al. (2016); Kusina et al. (2018), we assume that the amount of heavy quarks produced per binary nucleon-nucleon collisions and their initial momentum distributions to be the same as in the reference proton-proton collisions. Using Monte Carlo, we then sample the heavy quark initial transverse momenta by using the distributions coming from pQCD FONLL calculations Cacciari et al. (1998, 2001) in proton-proton collisions, choosing their central predictions in the renormalization scale range. The initial azimuthal directions of the heavy quark momenta are chosen randomly. We also sample the spatial distribution of the heavy quarks following the medium initial energy density of each hydrodynamic event, which reflects the initial parton scatterings of the considered collision. Even though in reality the number of heavy quarks at mid-rapidities per heavy-ion collision is on the order of ten, we largely oversample the number of heavy quarks to $\sim 10^{7}$ for each hydrodynamic event. This large oversampling is necessary for the statistics of some of the observables under study (e.g. the higher order particle cumulants) and it can be justified by the limited number of hydrodynamic events used in the simulation as compared to the large number of collisions in the experiments. Each of our hydrodynamic events can then be seen as a typical event with geometrical properties that correspond to a large number of collisions in the experiment.

II.3 Evolution of the heavy quarks in the medium

During the first stage of the collision ($\tau\lesssim 1$ fm/c) we neglect the effects of heavy quark energy loss in the cold nuclear matter of the colliding ions Arleo and Peigné (2014). During the deconfined stage of the collision ($1\lesssim\tau\lesssim 10$ fm/c), because the heavy quark typical energy scale (such as its mass $m>1$ GeV) is much larger than the medium scale (temperature $T\sim 100$–400 MeV), we assume each heavy quark to propagate and interact inside the expanding medium either through an energy loss process along a straight line or via Brownian motion described by relativistic Langevin dynamics.

II.3.1 Energy loss models

To easily study the dependence of the observables on some common variables, we use a simple parametric model for the heavy quark energy loss per unit length previously introduced in Refs. Betz and Gyulassy (2014); Horowitz and Gyulassy (2011), which is given by

[TABLE]

where $T$ is the local medium temperature experienced by the heavy quark, $v_{\text{Q}}$ is the heavy quark velocity in the global laboratory frame, $f(T,v_{\text{Q}})$ is a function encoding the energy loss parametrization, $\zeta$ is a random variable related to the energy loss fluctuations, and $\varGamma_{\text{flow}}$ takes into account the boost from the rest frame of the moving medium cell to the global laboratory frame Baier et al. (2007). The $\varGamma_{\text{flow}}$ factor is given by

[TABLE]

where \gamma=1\Big{/}\sqrt{1-v_{\text{flow}}^{2}}, $v_{\text{flow}}$ and $\varphi_{\text{flow}}$ are respectively the local medium cell velocity and azimuthal angle, and $\varphi_{\text{Q}}$ is the azimuthal angle defined by the propagating heavy quark in the transverse plane. Note that the relation (2) is a light quark jet ($p\gg m$) approximation to the general formula,

[TABLE]

which can be derived following a procedure similar to the one used in Baier et al. (2007) but without assuming $v_{\text{Q}}\to 1$. Discrepancies between the two expressions appear when $p_{\text{Q}}\lesssim 2m_{\text{Q}}$, i.e. when $p_{\text{c}}<3$ GeV for charm and $p_{\text{b}}<10$ GeV for bottom quarks. The consequences of this generalized formula on the observables will be investigated in a future work.

As shown in Fig. 3, we have tested a few parameterizations of the $f(T,v_{\text{Q}})$ function in order to select the ones that can roughly reproduce $R_{\text{AA}}$ data in $\sqrt{s_{\text{NN}}}=$2.76\text{,}\mathrm{TeV}$$ PbPb collisions in the 0–10% centrality range. The parametrizations $f=\xi T^{2}$ and $f=\delta\gamma_{\text{Q}}v_{\text{Q}}T^{2}$ (with $\gamma_{\text{Q}}=1/\sqrt{1-v_{\text{Q}}^{2}}$) are inspired by conformal AdS/CFT calculations Gubser (2006), whereas $f=\alpha$ and $f=\beta\gamma_{\text{Q}}v_{\text{Q}}$ are inspired by Ref. Das et al. (2015), which showed that a non-decreasing drag coefficient near the cross-over transition is favored for a simultaneous description of heavy flavor $R_{\text{AA}}(p_{\text{T}})$ and $v_{2}(p_{\text{T}})$ (this is also supported by $T$-matrix calculations van Hees et al. (2008); Riek and Rapp (2010)). Finally, we also consider the temperature dependent non-conformal drag force dependent model $f=\lambda T^{2}F_{\rm drag}$, where $F_{\rm drag}$ has been evaluated using holographic models that describe lattice QCD thermodynamics in Rougemont et al. (2016). The $\alpha$, $\beta$, $\delta$, $\lambda$ and $\xi$ factors are proportionality coefficients. Fig. 3 shows that the two energy loss models which are independent of the heavy quark velocity lead to very similar $R_{\text{AA}}$, which increases with increasing $p_{\text{T}}$. Their temperature dependence does not seem to play any important role for $R_{\text{AA}}$. The two energy loss models which are velocity dependent give also similar results but lead to $R_{\text{AA}}$ continuously decreasing with $p_{\text{T}}$. The non-conformal drag force model has a strong dependence on the parton masses, leading to different trends for D and B mesons. Thus, the non-conformal drag force and velocity dependent models are favored by the low-$p_{\text{T}}$ $R_{\text{AA}}$ data. Nevertheless, as the high $p_{\text{T}}$ data favors energy loss models which lead to increasing $R_{\text{AA}}$ with increasing $p_{\text{T}}$, we limit our study in this paper to the two velocity independent energy loss models defined above. As further explained in Section II.5, their two factors $\alpha$ and $\xi$ are tuned in this work via a single experimental observable.

To test the effects of energy loss fluctuations on the observables, we adopt a description put forward in Betz and Gyulassy (2014) where energy loss fluctuations in Eq. (2) are implemented via a random variable $\zeta$ for each heavy quark propagating in the medium. Three different probability distributions for $\zeta$ have been implemented for the constant energy loss model ($f=\alpha$) in this work: a Gaussian distribution $f_{\zeta}^{\rm gauss}(\zeta)=1/(\sqrt{2\pi}\sigma)\,{\rm exp}[-(\zeta-1)^{2}/2\sigma^{2}]$ with the $\sigma$ parameter taken to be $0.3$, a uniform distribution $f_{\zeta}^{\rm uniform}(\zeta)=1/2$ with $0\leq\zeta\leq 2$, and finally a linear distribution $f_{\zeta}^{\rm linear}(\zeta)=2/3-2/9\,\zeta$ with $0\leq\zeta\leq 3$.

II.3.2 Relativistic Langevin dynamics

Due to the large separation of scales, $T\ll m_{\text{Q}}$, one can also describe the propagation of heavy quarks in the deconfined medium as a stochastic Brownian motion. The Brownian approximation is particularly valid for bottom quarks but only partially valid for charm quarks when the temperature exceeds 250 MeV Das et al. (2014). Brownian motion is commonly described via the Fokker-Planck equation and simulated in practice for each individual heavy quark through a relativistic Langevin equation, which naturally encodes both some energy gain through a fluctuating force and some energy loss through a drag term Moore and Teaney (2005); van Hees et al. (2006); He et al. (2012, 2013); Akamatsu et al. (2009); Horowitz (2015); Cao et al. (2015b); Cao and Bass (2011); Nahrgang et al. (2018); Xu et al. (2017); Nahrgang (2015). Assuming the diffusion coefficients to be isotropic and the momentum space diffusion coefficient $\kappa$ to be independent of the heavy quark momentum $\bm{p}$, one can write the relativistic Langevin equation as

[TABLE]

where the index $i=x,y$ corresponds to transverse plane coordinates, $\Gamma$ is the drag coefficient, and $\rho_{i}$ is the fluctuating force described classically by white noise in a Markovian process. For the heavy quarks to be able to reach the correct thermal equilibrium He et al. (2013), the relativistic Einstein fluctuation-dissipation relation between the diffusion coefficients sets

[TABLE]

where $D$ is the spatial diffusion coefficient. As the deconfined medium is rapidly expanding, one needs to perform the necessary Lorentz boosts between the local rest frame of the moving medium (in which the heavy quark interacts) and the global laboratory frame (where it propagates). To do so, we first boost the heavy quark 4-momentum from the global frame $p=(E,\bm{p})$ to the local rest frame of the medium cell $p^{\prime}=(E^{\prime},\bm{p^{\prime}})$, where we calculate the momentum update via the Langevin equation with a pre-point Ito implementation of the stochastic scheme,

[TABLE]

for a time step $\Delta t^{\prime}$ of the implementation in the local rest frame of the medium, given by

[TABLE]

where $\Delta t$ is the related time step in the global frame, $v^{\rm\,flow}=(\gamma,\gamma\bm{v}_{\rm flow})$ is the local 4-velocity of the medium cell, and $\gamma=1/\sqrt{1-\bm{v}_{\rm flow}^{2}}$ is the corresponding Lorentz factor. We then boost the heavy quark 4-momentum back to the global frame, where we compute its propagation in the transverse plane via

[TABLE]

Finally, we repeat these operations until reaching the hadronization temperature.

In this work, two parametrizations of the diffusion coefficients have been chosen. The first one, denoted by “M&T” in the following, is inspired by Moore and Teaney’s leading order QCD description of the scattering process and a Debye mass correction of the gluon propagator Moore and Teaney (2005). In this model, the spatial diffusion coefficient reduces to

[TABLE]

where $k_{\rm M\&T}$ is a factor — estimated to be around 6 in their work — that we use in this paper as a tunable parameter. The second parametrization adopted in this work, denoted by “G&A” in the following, comes from Gossiaux and Aichelin’s collisional model based on a running coupling constant and an optimized hard thermal loop correction of the gluon propagator Gossiaux and Aichelin (2009). To obtain a tunable parameter as in the other transport models explored in this work, we multiply the drag $A_{\rm G\&A}(T,p)$ [c/fm] directly obtained from this model by a factor $k_{\rm G\&A}$, such that we obtain

[TABLE]

A comparison between the two corresponding drag coefficients can be found in Fig. 5 (right). At fixed temperature $T=0.3$ GeV, the M&T parametrization gives a significantly larger drag at low momentum whereas it saturates at larger momentum where the G&A becomes larger and increases linearly. The Fig. 5 (left) illustrates the temperature dependence of the dimensionless quantity $D(2\pi T)$ computed at zero momentum in the two models. While this quantity is constant by definition in the M&T parametrization, for G&A it acquires a temperature dependence approaching the results found using M&T from below. As compared to other models Andronic et al. (2016); Xu et al. (2018); Cao et al. (2019), the present drag and diffusion coefficients are rather small.

II.4 Fragmentation and decay

We assume each heavy quark to propagate until it reaches a cell where the temperature of the medium is lower or equal to a decoupling temperature $T_{\text{d}}$. Inspired by other works Betz et al. (2017); Nahrgang et al. (2016), and by lattice QCD results on the cross-over transition Aoki et al. (2006, 2009); Borsanyi et al. (2010); Mukherjee et al. (2016), we choose a range of $T_{\text{d}}$ from 120 MeV to 160 MeV in order to assess part of the uncertainties related to the complicated process of hadronization. In this work, fragmentation is assumed to be the same as in pp collisions (for medium modified fragmentation functions see Arleo (2009); Cacciari et al. (2013)). To perform the fragmentation, we use the Peterson fragmentation function

[TABLE]

as a probability distribution to obtain the fraction $z$ of the fragmenting heavy quark $E_{\text{Q}}+p_{\text{Q}}$ taken by the daughter hadron $E_{\text{H}}+p_{\text{H}}=z(E_{\text{Q}}+p_{\text{Q}})$. The Peterson function parameters $\epsilon_{\text{c}}$ and $\epsilon_{\text{b}}$ are chosen such as to reproduce ${{\text{D}}^{0}}$ and ${{\text{B}}^{0}}$ meson FONLL spectra in pp collisions Cacciari et al. (1998, 2001). Following the fragmentation, we do not consider the possible interactions between the produced heavy mesons and the hadronic gas Torres-Rincon et al. (2016). Finally, we perform the decays of the ${{\text{D}}^{0}}$ and ${{\text{B}}^{0}}$ mesons using Pythia 8 Sjostrand et al. (2008), focusing only on their semileptonic channels.

II.5 Determination of the tunable coupling factor of each transport model

In each of the transport models described above there is a free parameter, i.e. the coupling factors $\alpha$ and $\xi$ for the chosen energy loss models and the factors $k_{\rm M\&T}$ and $k_{\rm G\&A}$ for the Langevin models, which must be fixed. This is done by finding the value that gives the best description of differential $R_{\text{AA}}$ for $p_{\text{T}}\gtrsim 10$ GeV at LHC energies or $p_{\text{T}}\gtrsim 5$ GeV at RHIC obtained from the simulation compared to available experimental data for each model, heavy quark type, limiting values of the decoupling temperature considered ($T_{\text{d}}=120$ and $160$ MeV), and collision energy. To set the factors for the charm quark simulations, we use ${{\text{D}}^{0}}$ meson differential $R_{\text{AA}}$ data in one centrality range (0–10% unless specified otherwise). The obtained values of these constants are summarized in table 1. As almost no data for ${{\text{B}}^{0}}$ mesons is available, to obtain the constants for bottom quarks we use heavy flavor electron differential $R_{\text{AA}}$ data in most central collisions. The idea is to use the already fixed coupling factor for charm quarks and vary the bottom factor to find its value such that the total contribution to the electron $R_{\text{AA}}$ matches the data. Note that here we use more recent experimental data in PbPb collisions at $\sqrt{s_{\text{NN}}}=$5.02\text{,}\mathrm{TeV}$$ than in our previous paper Prado et al. (2017). The values of the constants for bottom quarks can be found in table 2. The values for the overall coupling factors for charm and bottom quarks are of the same order of magnitude and their ratios are not related to the quark mass ratio. For the energy loss models, bottom quark constants are generally smaller than charm quark factors, i.e. bottom quarks require a smaller coupling to the medium to fit the most central $R_{\text{AA}}$ data than charm quarks.

II.6 Examples of heavy quark spatial evolution

As explained above, in DAB-MOD we oversample each initial condition with a large number of heavy quarks. An example of the heavy quarks initial coordinates in the transverse plane for a central collision is shown at the top-left plot of Fig. 6. The distribution of heavy quarks follows the spatial fluctuations in energy density of the underlying medium. Each heavy quark then evolves and interacts with the bulk following a given transport model until it reaches a cell where the temperature of the medium is lower or equal to the decoupling temperature $T_{\text{d}}$. The corresponding final distributions for the decoupling temperatures $T_{\text{d}}=160$ MeV and $T_{\text{d}}=120$ MeV are shown, respectively, on the center-top and right-top of Fig. 6. At $T_{\text{d}}=160$ MeV, the heavy quark spatial distribution has a size similar to the initial distribution and is quite homogeneous. One can also note the presence of a relatively thick peripheral area including some high density “filaments” which are mainly composed of high $p_{\text{T}}$ heavy quarks and which density patterns are clearly correlated to the initial density patterns. At $T_{\text{d}}=120$ MeV, the distribution is wider and mainly homogeneous. The peripheral area is now thinner and more ring-shaped with the filaments being closer to each other as compared to $T_{\text{d}}=160$ MeV. Going to smaller $T_{\text{d}}$ seems therefore to partially wash out the initial density fluctuations. The same observations can be made in the case of a more peripheral event as shown at the bottom of Fig. 6 and with other transport models.

III Effect of transport model and decoupling temperature on observables in the case of MCKLN initial conditions and no coalescence

In this section we use MCKLN initial conditions and only fragmentation as the hadronization mechanism. A systematic study is performed where we:

• •

compare different energy loss parameterizations vs. different diffusion coefficients within the Langevin model

• •

vary the distributions of energy loss fluctuations

• •

determine the effect of the decoupling temperature $T_{\text{d}}$

• •

investigate the differences between D and B mesons, ${{\text{e}}^{\pm}}$, and muons

• •

study the centrality and beam energy dependence of our calculations.

III.1 Nuclear modification factor

The nuclear modification factor is essentially the ratio between the particle spectrum in AA collisions, $\mathop{}\!\text{d}N_{\text{AA}}/\mathop{}\!\text{d}p_{\text{T}}$, and the spectrum in pp collisions, $\mathop{}\!\text{d}N_{pp}/\mathop{}\!\text{d}p_{\text{T}}$, with a normalization factor $\mathcal{N}$ defined in terms of the number of binary collisions Miller et al. (2007). This gives:

[TABLE]

where we leave our calculations dependent on $\phi$ (the azimuthal angle in the plane transverse to the beam direction) in order to calculate the Fourier harmonics. Here, in practice, the $N\times\mathop{}\!\text{d}N_{pp}/\mathop{}\!\text{d}p_{\text{T}}$ spectrum is obtained using the same ingredients as described in Sec. II but in this case one turns off the heavy quark transport equations while hadronization is still done via fragmentation. In this work, the calculations are assumed to be boost invariant so comparisons are made only to mid-rapidity experimental data. Integrating over $\phi$ then reproduces the typical $R_{\text{AA}}(p_{\text{T}})$ that can be compared to experiment. Because we oversample the number of heavy quarks for each hydrodynamic background, we are able to reconstruct the entire $R_{\text{AA}}(p_{\text{T}},\phi)$ for each event.

As explained in Sec. II.5, the values of the free parameters of the transport models are determined by the best fit to most central high $p_{\text{T}}$ differential $R_{\text{AA}}$ experimental data and, thus, the magnitude of $R_{\text{AA}}(p_{\text{T}})$ in most central collisions obtained with DAB-MOD is not a prediction. However, its centrality and $p_{\text{T}}$ dependence are legitimate predictions of the model.

III.1.1 ${{\text{D}}^{0}}$ meson

In Fig. 7 we compare the results from different transport models for ${{\text{D}}^{0}}$ mesons in the 0–10% centrality of AuAu 200 GeV (left) and PbPb 5.02 TeV (right) collisions. At first sight, the four chosen transport models exhibit correct trends at intermediate/high $p_{\text{T}}$ when compared to data. However at low $p_{\text{T}}$, the energy loss and Langevin models lead to very different behaviors, clearly favoring Langevin models in the comparison to data. This difference originates in the fluctuating force within the Langevin approach: the balance between this force and the energy loss brought by the drag term maintains a certain amount of heavy quarks within the typical momentum range of the plasma particles. On the other hand, the energy loss models, where there is no energy gain, lead to a continuous decrease of $R_{\text{AA}}$ towards lower $p_{\text{T}}$ until a certain value (much) smaller than 1 GeV where a large amount of heavy quarks “accumulate” (not shown in the present figures).

Looking more carefully at Fig. 7 one can see that in Langevin models the choice of the diffusion coefficient leads to different results at high $p_{\text{T}}$. M$\&$T is able to reasonably capture the high $p_{\text{T}}$ behavior of $R_{\text{AA}}(p_{\text{T}})$ data whereas G$\&$A is found to be below the data333If we would have compared these predictions to the ALICE collaboration $R_{\text{AA}}(p_{\text{T}})$ data Acharya et al. (2018d), the conclusion might have been different since that data is significantly lower at high $p_{\text{T}}$.. At low $p_{\text{T}}$ they exhibit nearly identical behavior. While the results are somewhat lower than experimental data, we find (see Sec. V) that the addition of coalescence shifts the $R_{\text{AA}}(p_{\text{T}})$ curve upwards at low $p_{\text{T}}$ matching experimental data. The results for the transport models for different beam energies are also shown in Fig. 7. As explained in Sec. II.3.1, for the energy loss models we consider either a constant energy loss function, $f(T,v_{\text{Q}})=\alpha$, or a temperature dependent one where $f(T,v_{\text{Q}})=\zeta T^{2}$. Despite this difference, the ${{\text{D}}^{0}}$ $R_{\text{AA}}(p_{\text{T}})$ results are almost identical, which shows that this temperature dependence is not actually realized in practice in our calculations. Additionally, we do not find a strong dependence with beam energy when it comes to the qualitative differences between different types of heavy quark evolution. We consistently find that the Langevin model produces the largest $R_{\text{AA}}(p_{\text{T}})$ at low $p_{\text{T}}$ and that all the transport models behaves correctly at high $p_{\text{T}}$.

The point where a heavy quark decouples from the medium is still a source of uncertainty in heavy flavor modeling. In this work we investigate this issue by considering a range of decoupling temperatures between 120 and 160 MeV, which creates the bands in our theoretical calculations in following sections. In Fig.  7 $R_{\text{AA}}(p_{\text{T}})$ is plotted444We have checked that all experimental observables show either a monotonic increase or decrease with $T_{\text{d}}$ so it is sufficient to plot just the extrema to obtain the bands. for $T_{\text{d}}=120$ MeV and $T_{\text{d}}=160$ MeV. $R_{\text{AA}}(p_{\text{T}})$ for both beam energies does not exhibit a strong dependence on the decoupling temperature, which can be explained within our model by the magnitude of the high $p_{\text{T}}$ $R_{\text{AA}}$ being calibrated with the data for each decoupling temperature. Nevertheless, a significant difference can be seen between the decoupling temperatures at low $p_{\text{T}}$ within the Langevin model: $R_{\text{AA}}$ at $T_{\text{d}}=120$ MeV is larger than at $T_{\text{d}}=160$ MeV. In other words, in this situation the low $p_{\text{T}}$ heavy quarks are less suppressed if they are coupled to the plasma for longer. This ordering can be explained by the fluctuating force having a longer time to act when $T_{\text{d}}=120$ MeV to increase the momenta of (the larger amount of) heavy quarks initially at lower $p_{\text{T}}$. This effect outclasses an expected countereffect: at lower temperatures the thermal distribution of the heavy quarks shifts towards lower momenta, such that the $R_{\text{AA}}$ at $T_{\text{d}}=120$ MeV should be smaller than at $T_{\text{d}}=160$ MeV. This shows that quasi-equilibrium is far from being reached for the heavy quarks when $T_{\text{d}}=160$ MeV. As we will see below, the flow harmonics exhibit a different behavior and, in general, one finds that lower decoupling temperatures produce larger azimuthal anisotropies.

We also include different types of energy loss fluctuations in the case of the constant energy loss model (see Sec. II.3.1 for details). As shown in Fig. 8, Gaussian fluctuations do not make a significant difference. However, linear or uniform fluctuations suppress $R_{\text{AA}}$ at high $p_{\text{T}}$, showing a disagreement with the experimental data, although they increase $R_{\text{AA}}$ at low/intermediate $p_{\text{T}}$ improving the overall agreement with the data. However, all the energy loss results are within the error bars of the experimental data points so, unfortunately, no conclusion can be drawn for now regarding the specific form of the distribution of energy loss fluctuations within our model.

To test if the calibration made in the most central collisions gives reasonable results for the different models in other centralities, some predictions for more peripheral collisions are shown in Fig. 9. The different models show similar features found in 0–10% central collisions and fit correctly the less suppressed $R_{\text{AA}}$ data at intermediate/high $p_{\text{T}}$.

III.1.2 ${{\text{B}}^{0}}$ meson

Because ${{\text{B}}^{0}}$ production is less subject to initial state effects and final hadronic rescaterings compared to ${{\text{D}}^{0}}$, its observables are cleaner probes of the deconfined medium making it possible to study how changing the heavy quark mass affects its dynamics and hadronization. As explained in Sec. II.5, we use heavy flavor electron $R_{\text{AA}}$ data in the most central collisions to calibrate the bottom quark transport model coupling factors, such that all the ${{\text{B}}^{0}}$ meson results are predictions of the model. As shown in Fig. 10, the (limited) comparison between ${{\text{B}}^{+}}$ 0–100% experimental data and the ${{\text{B}}^{0}}$ meson $R_{\text{AA}}(p_{\text{T}})$ obtained with DAB-MOD in 0–50% centrality range illustrates the consistency of our calibration.

In Fig. 11 we compare different energy loss parametrizations to Langevin results for ${{\text{B}}^{0}}$ mesons in the 0–10% centrality class of AuAu 200 GeV (left) and PbPb 5.02 TeV (right) collisions. Comparing to Fig. 7, one can see that the ${{\text{B}}^{0}}$ and the ${{\text{D}}^{0}}$ meson $R_{\text{AA}}(p_{\text{T}})$ exhibit similar qualitative behavior. Nevertheless, at low and intermediate $p_{\text{T}}$, the ${{\text{B}}^{0}}$ $R_{\text{AA}}$ is generally less suppressed by a factor $\sim 2$ while its value at high $p_{\text{T}}$ is larger than that found for the ${{\text{D}}^{0}}$ meson. We note that in this $p_{\text{T}}$ regime the discrepancy between energy loss and Langevin approaches is much larger for ${{\text{B}}^{0}}$ than for ${{\text{D}}^{0}}$ mesons (there is also a larger difference between different Langevin parametrizations). At very high $p_{\text{T}}$, the values for ${{\text{B}}^{0}}$ and ${{\text{D}}^{0}}$ $R_{\text{AA}}$ are very similar since their mass difference becomes negligible when compared to the momentum. Finally, since the heavy quark mass has no impact on energy loss fluctuations, $R_{\text{AA}}(p_{\text{T}})$ of both ${{\text{B}}^{0}}$ and ${{\text{D}}^{0}}$ mesons computed within the constant energy loss model are similarly modified by energy loss fluctuations.

III.1.3 ${{\text{e}}^{\pm}}$ and ${{\mu}^{\pm}}$ from heavy flavor

Since very limited B meson data is available, one can instead study electrons (or muons) from heavy flavor decays since those are more commonly measured. In Fig. 12 the results for the heavy flavor electron $R_{\text{AA}}$ from semileptonic decays, ${\text{c}},{\text{b}}\rightarrow{{\text{e}}^{\pm}}$, are shown for two beam energies. As with the D and B meson $R_{\text{AA}}(p_{\text{T}})$, the same hierarchy involving the heavy quark evolution is seen in the electrons and the results fit reasonably well the data at intermediate and high $p_{\text{T}}$ for different beam energies and centralities (see Fig. 13). An exception can be observed in AuAu collisions at $\sqrt{s_{\text{NN}}}=$200\text{,}\mathrm{GeV}$$ where the models have trouble to reproduce the sudden increase of $R_{\text{AA}}$ in the data found around $p_{\text{T}}=6$ GeV. However, we note that those points are located at $p_{\text{T}}<10$ GeV, which is a regime that can be affected by coalescence and other initial/final hadronic effects. Additionally, the error bars in that region are quite large.

The ATLAS collaboration recently compared our results to muon data in PbPb $\sqrt{s_{\text{NN}}}=$2.76\text{,}\mathrm{TeV}$$ collisions Aaboud et al. (2018) and found that our constant energy loss results were able to reproduce $R_{\text{AA}}$ for $p_{\text{T}}\gtrsim 10$ GeV (at that point we did not have Langevin evolution in DAB-MOD). Note that the results obtained in the electron and muon semileptonic decay channels are almost identical down to 1 GeV as the mass of the leptons that stem from heavy flavor decays is negligible compared to their momentum (see, for instance, Figs. 14 and 32).

Fig. 15 shows a comparison of the separated charm and bottom hadron decay electrons ${\text{c}}\rightarrow{{\text{e}}^{\pm}}$ and ${\text{b}}\rightarrow{{\text{e}}^{\pm}}$, respectively, for two different centrality ranges. Similarly to RHIC data, we observe that the electrons from bottom hadrons are less suppressed than those from charmed hadrons in the low/intermediate $p_{\text{T}}$ regime. The electron $R_{\text{AA}}$ from ${{\text{D}}^{0}}$ mesons reproduces reasonably well RHIC data in both centrality ranges, but we underestimate the ${\text{b}}\rightarrow{{\text{e}}^{\pm}}$ $R_{\text{AA}}$ data. Part of this discrepancy might originate from the difficult calibration of the bottom quark transport model coefficients using the heavy flavor data at RHIC (because of the sudden $R_{\text{AA}}$ increase and large error bars shown in Fig. 12). In Fig. 16 we show a comparison to the (low $p_{\text{T}}$) ${\text{b}}\rightarrow{{\text{e}}^{\pm}}$ LHC data at $\sqrt{s_{\text{NN}}}=$5.02\text{,}\mathrm{TeV}$$ where the heavy flavor electron data can be used with more confidence for the calibration. As for the ${{\text{D}}^{0}}$ data in Fig. 7, the Langevin models give a better agreement.

If one puts aside the effects of hadronization and initial/final hadronic stages on heavy meson production, a comparison between ${\text{c}}\rightarrow{{\text{e}}^{\pm}}$ and ${\text{b}}\rightarrow{{\text{e}}^{\pm}}$ $R_{\text{AA}}$ illustrates the mass dependence of the energy exchanges between the heavy quarks and the QGP medium. Within the usual energy loss framework, the experimental observation $R_{\text{AA}}^{\text{c}}<R_{\text{AA}}^{\text{b}}$ is in agreement with the theoretical prediction of the quark mass hierarchy for radiative energy loss, $\Delta E_{\text{c}}>\Delta E_{\text{b}}$, mainly due to the dead cone effect increasing with mass Dokshitzer and Kharzeev (2001). In DAB-MOD, the heavy quark mass plays a role in different parts of the simulation: 1) inherently in the formalism of the transport models: the evolution depends on mass even without any direct mass dependence in the Langevin and energy loss parametrizations, 2) in the values of the transport model free coupling factors for the charm and bottom quarks, although since they are in general of the same order (Tab. 1 and 2) this is only a small effect, 3) in the initial spectra from FONLL and 4) embedded in the hadronization process (both in the case of fragmentation and coalescence, see Sec. V).

III.2 Azimuthal anisotropy and multi-particle cumulants

In a previous work involving all charged high $p_{\text{T}}$ hadrons, multi-particle cumulants were predicted Betz et al. (2017) and later measured by CMS Sirunyan et al. (2018a). For a correlation of $m$ particles, because of the low statistics at high $p_{\text{T}}$, one hard particle is correlated with $m-1$ soft particles. In that case, when both particles come from the same $p_{\text{T}}$ bin, one can write the flow vectors $V_{n}$, accounting for both real and imaginary parts, as

[TABLE]

Here we use the same notation as in Gardim et al. (2012, 2015); Betz et al. (2017) for the sake of consistency. However, in the hard sector we are considering the particle of interest from separate $p_{\text{T}}$ bins so a $p_{\text{T}}$ dependence appears

[TABLE]

This leads to a $p_{\text{T}}$ dependent two particle correlation that uses the complex conjugate to produce a real-valued result as in

[TABLE]

We then consider the typical 2-particle correlations by taking one hard and one soft particle and the 4-particle correlation by taking 1 hard particle and 3 soft particles that are averaged over events within a fixed centrality window. They are defined by

[TABLE]

where

[TABLE]

Here, the outer bracket $\langle\ldots\rangle$ is an artifact of centrality rebinning where in experiments finer centrality bins are taken e.g. 0.5% centralities, which we indicate as $j$, that are then recombined using multiplicity weighing in a wider centrality bin of a width of 10%, for example. The inner brackets indicate averaging over the events within the $j$th fine centrality bin.

Initially, in Ref. Betz et al. (2017) it was expected that $v_{2}\{4\}/v_{2}\{2\}(p_{\text{T}})\sim\text{constant}$ at low $p_{\text{T}}$. This is because the ratio of $v_{2}\{4\}/v_{2}\{2\}(p_{\text{T}})$ encapsulates a non-trivial interplay between the event-by-event fluctuations of energy loss and the initial condition fluctuations, with the latter being typically associated only with the soft sector (though it was found in Betz et al. (2017) that these fluctuations also influence the multi-particle cumulants at high $p_{\text{T}}$ as well). In fact, in Eq. (A2) of Ref. Betz et al. (2017), it was shown that the exact interplay between event-by-event initial condition fluctuations and hard physics fluctuations is given by

[TABLE]

where the terms $\frac{\langle v_{n}^{4}\rangle}{\langle v_{n}^{2}\rangle^{2}}$ and $\frac{v_{n}\{2\}}{v_{n}\{4\}}$ are determined from initial condition fluctuations, which are translated into the final flow harmonics via linear+cubic response Noronha-Hostler et al. (2016d); Sievert and Noronha-Hostler (2019) within the soft sector, i.e., in the context of relativistic hydrodynamics. The only term that contains contributions from the hard sector is then $\frac{\langle v_{n}^{2}V_{n}V_{n}^{*}(p_{\text{T}})\rangle}{\langle v_{n}^{2}\rangle\langle V_{n}V_{n}^{*}(p_{\text{T}})\rangle}$, which can be interpreted as a measure of azimuthal anisotropy fluctuations associated with hard physics. When there are only soft fluctuations

[TABLE]

and the term in parentheses in Eq. (23) is nonzero. We should note that it is not known a priori if soft fluctuations are larger or smaller than hard fluctuations and, thus, the term $\left(\frac{v_{n}\{2\}}{v_{n}\{4\}}\right)^{4}\left(\frac{\langle v_{n}^{4}\rangle}{\langle v_{n}^{2}\rangle^{2}}-\frac{\langle v_{n}^{2}V_{n}V_{n}^{*}(p_{\text{T}})\rangle}{\langle v_{n}^{2}\rangle\langle V_{n}V_{n}^{*}(p_{\text{T}})\rangle}\right)$ may contribute positively or negatively to the total $\frac{v_{n}\{4\}(p_{\text{T}})}{v_{n}\{2\}(p_{\text{T}})}$. However, when a deviation from Eq. (25) occurs this implies that some other type of physics is occurring that does not stem from the soft sector. Therefore, this illustrates the importance of direct theory to experiment comparisons of $\frac{v_{n}\{4\}(p_{\text{T}})}{v_{n}\{2\}(p_{\text{T}})}$ in the hard sector.

Previous work had investigated 2-particle correlations coming from the same $p_{\text{T}}$ bin Heinz et al. (2013), for instance, 2-particle correlations up to $\sim$10\text{,}\mathrm{G}\mathrm{e}\mathrm{V} Aad *et al.* ([2012](#bib.bib104)); Zhou ([2014](#bib.bib105), [2019](#bib.bib106)). However, when it comes to energy loss models that are not fully integrated into hydrodynamics, they only typically become valid above $p_{\text{T}}>$10\text{\,}\mathrm{G}\mathrm{e}\mathrm{V} so those previous measurements require additional physics to interpret their results. More recently, multi-particle cumulant calculations have been extended to the heavy flavor sector Prado et al. (2017) by correlating heavy flavor particles with low $p_{\text{T}}$ charged particles. However, they have not yet been systematically studied in the heavy flavor sector nor have they been measured experimentally. In the following, we systematically compare the influence of different medium transport properties and initial conditions on heavy flavor multi-particle cumulants for the first time.555We note that in Xu et al. (2019) a systematic study was done with different energy loss models for $R_{\text{AA}}$ and $v_{2}$ but higher order flow harmonics and multi-particle cumulants were not yet considered.

III.2.1 Elliptic flow from two-particle cumulants

First we test the influence of different energy loss models. As a test bed we use PbPb collisions at $\sqrt{s_{\text{NN}}}=$5.02\text{,}\mathrm{TeV}$$ with MCKLN initial conditions, which are the same hydrodynamic backgrounds used in Betz et al. (2017). In Fig. 17 results for ${{\text{D}}^{0}}$ and ${{\text{B}}^{0}}$ mesons are shown for the same energy loss models previously shown in Fig. 3 for $R_{\text{AA}}$. We notice that the different models can lead to a variety of different curves for $v_{2}$. At the same time, models that are indistinguishable using only the $R_{\text{AA}}$ calculations are clearly separated when considering the $v_{2}$ results, as is the case for the two energy loss parametrizations $f=\alpha$ and $f=\xi T^{2}$. These two particular models differ from the other three in that they do not have an explicit dependence on the heavy quark momentum and lead to a bump in the low $p_{\text{T}}$ region as is expected from experimental data. This observation agrees with the previous choice of energy loss parametrizations made using $R_{\text{AA}}$ data.

The selected energy loss models are compared with the two different parametrizations within the Langevin formalism in Fig. 18 for a range of decoupling temperatures. On the left side we observe that the M&T parametrization seems to better capture the characteristic bump in the low $p_{\text{T}}$ regime present in the experimental data while other models seem to overlap in AuAu collisions at RHIC with $\sqrt{s_{\text{NN}}}=$200\text{,}\mathrm{GeV}. In the case of PbPb $\sqrt{s_{\text{NN}}}=$5.02\text{\,}\mathrm{TeV} collisions at LHC, the different models seem to be more distinguishable with the constant energy loss model leading to the largest values of $v_{2}$ at large $p_{\text{T}}$, while the Langevin model with M&T parametrization remains closest to experimental data at low $p_{\text{T}}$. Results for $v_{2}$ show better agreement with data for the lowest collision energy overall, while being slightly underestimated in the other case.

We take a special look at how the three parametrizations of energy loss fluctuations affect the results for the anisotropic flow coefficient. These results are shown in Fig. 19. The Gaussian fluctuation Ansatz leads to a more subtle effect by slightly lowering $v_{2}$ at low $p_{\text{T}}$. The other two models lead to stronger effects and, even though they have very different functional forms, the $v_{2}$ results do not show much difference among them. These results suggest that energy loss fluctuations are more relevant to the calculation of the low $p_{\text{T}}$ regime of the elliptic flow coefficient.

It is also possible to compare the simulation results with experimental data for central collisions, as it is shown in Fig. 20. The same behavior observed in Fig. 18 is also found in this case when comparing the different models. However, the experimental data comparison shows a better agreement with LHC collisions rather than RHIC’s, in which case the results from the simulation underestimate the data except for very low $p_{\text{T}}$ regime. Since these comparisons very much rely on the low $p_{\text{T}}$ region, we expect that the introduction of coalescence in the simulations to change these results.

In Fig. 21 we show predictions for $v_{2}$ of ${{\text{B}}^{0}}$ mesons obtained from our simulations, using only the constant energy loss parametrization and the M&T parametrization within the Langevin formalism, for the same collision setup that was used for the D meson results. We notice a similar behavior between both models with Langevin’s leading to a larger $v_{2}$ at low $p_{\text{T}}$ though it becomes lower at high $p_{\text{T}}$.

Having both results for ${{\text{B}}^{0}}$ and ${{\text{D}}^{0}}$ mesons, we can now obtain the elliptic flow coefficient for heavy flavor electrons from decays. The comparison between the different transport models for this case is shown in Fig. 22. We observe that the Langevin parametrizations tend to better reflect the features of the experimental data, especially for non-central collisions, shown in the right plot of the figure. In the case of central collisions, although the results for the ${{\text{D}}^{0}}$ meson underestimated the data, a good agreement with data is observed for the heavy flavor electrons for both Langevin parametrizations. Also, the difference between the models in the mid-$p_{\text{T}}$ range up to $\approx$30\text{,}\mathrm{G}\mathrm{e}\mathrm{V}$$ is not as pronounced for electrons as it is for the heavy mesons.

We further explore the role of each heavy quark, bottom or charm, in building up the elliptic flow of heavy flavor electrons in Fig. 23 where the results for AuAu collisions at $\sqrt{s_{\text{NN}}}=$200\text{,}\mathrm{GeV}$$ are shown compared to experimental data from the PHENIX collaboration. A good agreement is obtained considering the large uncertainties in the measured values. Unfortunately, the comparison is limited to the lower $p_{\text{T}}$ range where the physics is considerably more complex and other effects might be important.

III.2.2 Elliptic flow from multi-particle cumulants

Up to this point we only explored the correlation between the particles using 2-particle cumulants. We now test the influence of different energy loss and Langevin diffusion coefficients on the calculation of the 4-particle cumulants. For the 30–40% centrality class in the soft sector MCKLN initial conditions produce $v_{2}\left\{4\right\}/v_{2}\left\{2\right\}\approx 0.91$. In Fig. 24 the corresponding heavy flavor predictions are shown. Similar to the high $p_{T}$ particles result in Betz et al. (2017), the heavy flavor predictions are $v_{2}\left\{4\right\}/v_{2}\left\{2\right\}\approx 0.95$. There is a slight downward shift if a higher decoupling temperature is considered, which may provide some insight into the optimal values of the decoupling temperature parameter. Additionally, energy loss fluctuations have an opposite trend in the $p_{\text{T}}$ dependence compared to Langevin results. In the energy loss approach a peak in $v_{2}\{4\}/v_{2}\{2\}(p_{\text{T}})$ is seen at low $p_{\text{T}}$ whereas in the Langevin scenario $v_{2}\{4\}/v_{2}\{2\}(p_{\text{T}})$ increases at high $p_{\text{T}}$.

Furthermore, it is also instructive to investigate the centrality dependence of this observable. Thus, we also consider $v_{2}\{4\}/v_{2}\{2\}$ integrated in the range $8<p_{\text{T}}<40$ GeV. In Fig. 25 different types of energy loss fluctuations and Langevin diffusion transport coefficients are shown (left) for a fixed decoupling temperature. On the right, the decoupling temperatures are compared for our two best fit setups. One of the biggest takeaways from Fig. 25 is that the central collision region is the most sensitive regime in the description of heavy flavor dynamics. More peripheral collisions are predicted to have nearly identical results regardless of the underlying assumptions made in the heavy flavor modeling. However, peripheral collisions are more sensitive to the decoupling temperature (and less sensitive to the heavy flavor description) so by investigating 0–10% and 40–50% centrality classes one may be able to constrain both simultaneously. The caveat relies on the ability of experiments to obtain these results with reasonable error bars since the differences are small.

It is also interesting to consider possible mass effects on the $v_{2}$ fluctuations of heavy flavor mesons. Up to this point, our previous results only considered ${{\text{D}}^{0}}$ mesons but here we compare ${{\text{D}}^{0}}$ and ${{\text{B}}^{0}}$ mesons in Fig. 26. We find that ${{\text{B}}^{0}}$ mesons have a slightly larger $v_{2}\{4\}/v_{2}\{2\}$ in central collisions but the effect is very small. Finally, we show how the ${{\text{D}}^{0}}$ meson $v_{2}\{4\}/v_{2}\{2\}$ changes with collision energy and system size in Fig. 27. Again, the largest difference occurs in central collisions while, for peripheral collisions, results at LHC energies become similar though still distinct from the RHIC result.

III.2.3 Triangular flow with two-particle cumulants

As done previously for the elliptic flow coefficient, we first test how the different parametrizations of Langevin and the selected energy loss models affect the triangular flow coefficient $v_{3}$ using the correlation between two particles. From the plots of Fig. 28 one can observe much larger bands when studying the same $T_{\text{d}}$ range in comparison with the results for $v_{2}$. This wider band effect is consistent with a time hierarchy where $v_{2}$ is built up first in the evolution of the system and $v_{3}$ would be generated later. In this sense, a larger value of $T_{\text{d}}$ would hadronize the heavy quarks before they had enough time to build up a significant $v_{3}$ flow. Despite this difference, the general behavior is maintained among the different models with the constant energy loss model leading to higher values of $v_{3}$ in the high $p_{\text{T}}$ regime, while the low $p_{\text{T}}$ region is dominated by the M&T Langevin parametrization, leading to the best agreement with experimental data at $\sqrt{s_{\text{NN}}}=$200\text{,}\mathrm{GeV} despite the low number of data points. As with the $v_{2}$ results, the models tend to underestimate the data at $\sqrt{s_{\text{NN}}}=$5.02\text{\,}\mathrm{TeV}.

Energy loss fluctuations are also studied and the corresponding results are shown in Fig.  29 in which we observe a similar trend to that found for $v_{2}$. While the Gaussian fluctuations slightly decrease the flow coefficient values in the lower $p_{\text{T}}$ regime, stronger fluctuations such as the linear and the constant ones seem to have a much bigger effect. In addition, even though their functional form differ greatly, the overall effect of energy loss fluctuations on $v_{3}$ is not very different for the different cases considered. In the large $p_{\text{T}}$ regime the effect of the fluctuations seems to be negligible.

The corresponding results for $v_{3}$ in central collisions are shown in Fig. 30. We notice the scale between the different models remain unchanged, in agreement with all the calculations for $v_{2}$ and $v_{3}$ shown before. The simulation results still underestimate the experimental data for the largest beam energy collisions in agreement with the results observed in Fig. 28. One feature that is specific for these $v_{3}$ calculations is that there seems to be very little centrality dependence for this observable, differing from the $v_{2}$ observations. In that respect, all models seem to agree with the observed conclusion from data.

In Fig. 31 we show a prediction for the triangular flow of ${{\text{B}}^{0}}$ mesons using the constant energy loss model and the M&T Langevin parametrization. We observe a much lower $v_{3}$ in comparison with the results for the D meson, which also occurs for $v_{2}$, suggesting that the mass hierarchy may propagate throughout the higher flow harmonics.

III.2.4 4th order flow coefficient from two-particle cumulants

In Fig. 32 we compare the results for the 4th order flow coefficient from 2-particle correlations, $v_{4}\{2\}$, involving heavy flavor electrons and muons in the 0–10% (left) and $30-40\%$ (right) centrality classes of PbPb 2.76 TeV collisions, obtained using the constant energy loss and the M&T Langevin model. As already observed for the $R_{\text{AA}}$ in Fig. 14, one can first note that the electron and muon channels give equivalent results down to 1 GeV, so that comparing our electron results with muon data at mid-rapidities is fine. The $v_{4}$ results behave similarly to $v_{2}$ and $v_{3}$: they all decrease with increasing $p_{\text{T}}$ for $p_{\text{T}}\gtrsim 2$ GeV and, compared to the energy loss models, the Langevin framework leads to higher (smaller) values at low (high) $p_{\text{T}}$. Both models lie within the experimental data uncertainties although, as for the other flow coefficients, the $v_{4}$ results seem to slightly underestimate the data. We observe almost no variation of $v_{4}$ with centrality, as we did for $v_{3}$ in the previous section, which shows that $v_{4}$ and $v_{3}$ stem from geometrical fluctuations.

Despite the similarities with the other coefficients $v_{4}$ exhibits a new feature in $\sqrt{s_{\text{NN}}}=$2.76\text{,}\mathrm{TeV} PbPb collisions: for the decoupling temperature $T_{\text{d}}=160$ MeV the heavy flavor electron/muon $v_{4}$ is negative on the whole $p_{\text{T}}$ range for both transport models. The corresponding experimental data is either positive or negative but it is mostly compatible with zero within the error bands Aaboud *et al.* ([2018](#bib.bib17)). As shown in Fig.  [33](#S3.F33), this is also the case in 200 GeV AuAu collisions but not in $\sqrt{s_{\text{NN}}}=$5.02\text{\,}\mathrm{TeV} PbPb collisions. Additionally, as shown in Fig. 34 for 200 GeV AuAu collisions (note the different situation in Fig. 35 for 5.02 TeV PbPb collisions), both the ${{\text{D}}^{0}}$ and ${{\text{B}}^{0}}$ meson $v_{4}$ are negative when $T_{\text{d}}=160$ MeV with the ${{\text{D}}^{0}}$ meson $v_{4}$ being larger in absolute values.

The reason for these negative values is an anti-correlation between the heavy flavor and the bulk azimuthal anisotropy angles when $T_{\text{d}}=160$ MeV, $\psi_{4}^{\rm(heavy)}$ and $\psi_{4}^{\rm(soft)}$, respectively. Indeed, as shown in Fig. 36 the event-by-event $\psi_{4}^{\rm(heavy)}(\psi_{4}^{\rm(soft)})$ distributions are anti-correlated when $T_{\text{d}}=160$ MeV (centered on $\psi_{4}^{\rm(heavy)}=\pm\pi/2$ when $\psi_{4}^{\rm(soft)}=0$) whereas they are correlated when $T_{\text{d}}=120$ MeV (centered on $\psi_{4}^{\rm(heavy)}=0$ with a linear correlation to $\psi_{4}^{\rm(soft)}$). For $T_{\text{d}}=160$ MeV in 200 GeV AuAu and 2.76 TeV PbPb collisions the heavy quark interaction with the bulk medium is therefore not long enough for the heavy quark flow vector $V_{4}^{\rm(heavy)}$ to be in phase with the bulk $V_{4}^{\rm(soft)}$. This phenomenon was also observed in light quark models for $n>3$ Jia (2013). The $\psi_{4}$ observable will be further explored in the next section. Finally, by comparing the different flow harmonics we see a clear hierarchy $v_{2}>v_{3}>v_{4}$ for any collision energy as was also observed by the experiments Aaboud et al. (2018).

III.2.5 Correlation between the flow of heavy mesons and all charged particles

In addition to looking at the correlations between the heavy and the soft sectors from the point of view of the cumulants, it is also interesting to study the direct correlation between these sectors on an event-by-event basis using event-shape engineering techniques. Since the soft sector flow coefficients are directly related to the event eccentricities, this study can provide information on the role of the initial anisotropy on final observables for the heavy flavor sector. This work uses the same approach that has been previously used in Noronha-Hostler et al. (2016c) when investigating soft and hard sector correlations. We define the correlations by binning the distribution of integrated flow harmonics in the soft sector, $v_{n}^{\text{(soft)}}$, and then evaluating the corresponding flow coefficient of ${{\text{B}}^{0}}$ and ${{\text{D}}^{0}}$ mesons for each bin, $v_{n}^{\text{(heavy)}}$. The result can be understood as the probability that an event within a given soft sector flow class will correspond to a particular value of $v_{n}^{\text{(heavy)}}$. One can study the slope of these correlations with respect to different collision conditions leading to $v_{n}$ fluctuations. In that respect, if no $v_{n}$ fluctuations were to be observed, the plot would show a flat horizontal line.

Fig. 37 shows the correlation between the ${{\text{D}}^{0}}$ meson $v_{2}$ and the corresponding quantity for all charged particles computed using different transport models and different beam energies. We first notice that, regardless of the transport model used for the calculation, all results indicate that the correlation between the heavy and soft sectors is linear, although the slope of the lines might be different, reflecting the different $v_{n}$ already observed in previous sections for a given $p_{\text{T}}$ range. It is worth noticing that the hierarchy between the different transport models is not maintained for different energies, which may be due to the specific implementation details of each model.

The same approach can be used to investigate the same type of correlations involving higher order flow coefficients to check if the linear correlation observed for $v_{2}$ is still maintained at higher orders. In Fig. 38 it is possible to observe the same effect, although a deviation occurs when the soft $v_{3}$ becomes large. These deviations are related to the limited statistics of the event leading to a wider distribution of the heavy flavor sector $v_{3}^{\text{(heavy)}}$. Furthermore, we note that the hierarchy among the values obtained for the different transport models is similar to that found in the case of elliptic flow, except for the constant energy loss model at $200$ GeV.

One of the main questions involving heavy flavor quarks in the QGP concerns their coupling with the expanding medium. In our framework this is also encoded in the decoupling parameter $T_{\text{d}}$, which defines a temperature scale below which heavy quarks are considered to not be coupled with the medium anymore. This parameter affects the path length experienced by the heavy quark and can therefore greatly affect the results obtained for the flow coefficients. The lower the decoupling temperature, the longer the heavy quarks are under the influence of the medium in the transport models considered in this paper. In Fig. 39, a range of decoupling temperatures is studied for every transport model considering the flow of ${{\text{D}}^{0}}$ mesons. We notice that a variation of the decoupling temperature affects $v_{3}$ more than it affects $v_{2}$. Since $v_{2}$ is built up quickly, the difference in the slopes due to the decoupling temperature is very small, especially at very large collision energies which require more time to cool down below $T_{\text{d}}$. On the other hand, higher order harmonics and lower collision energies are more strongly affected by this parameter. In Fig. 40, the correlation between the $v_{4}$ of ${{\text{D}}^{0}}$ mesons and that of all charged particles is shown for the two collision energies considered. We obtain even wider bands compared with the $v_{3}$ correlation results, further indicating the hierarchy described before. These results also agree with previous observations regarding the hierarchy of flow harmonics determined using the different models.

All the analysis performed on ${{\text{D}}^{0}}$ mesons can also be performed on ${{\text{B}}^{0}}$ mesons and the results are very similar. The plots in Fig. 41 show results for ${{\text{B}}^{0}}$ mesons at different decoupling temperatures and for different models. In comparison with previous results for ${{\text{D}}^{0}}$ mesons we notice that the difference between transport models for ${{\text{B}}^{0}}$ mesons is less pronounced, though the decoupling temperature seems to play a bigger role in this case. Since this observable is integrated over $p_{\text{T}}$, both of these effects may be due to the mass difference between the mesons and a direct comparison in the same $p_{\text{T}}$ range may be misleading.

Not only the flow harmonics correlations are worth looking into. One can also study the correlations between their respective event plane angles $\psi_{n}$. This quantity gives information about the alignment between the event plane angles of the soft and the heavy sector. It is convenient to represent this correlation as a cosine term such as the ones that appear on the equations for the cumulants so they can be easily related to the $p_{\text{T}}$-differential results introduced earlier.

We first introduce a comparison for different collision energies of the event plane angles correlations in Fig. 42 in which both ${{\text{B}}^{0}}$ and ${{\text{D}}^{0}}$ mesons are studied using different values of the parameter $T_{\text{d}}$. The second order event plane angle is consistently aligned with \big{\langle}\cos\big{[}n\big{(}\psi_{2}^{\text{(soft)}}-\psi_{2}^{\text{(heavy)}}\big{)}\big{]}\big{\rangle}>0.95 across most of the centrality range. This is observed for all the different settings using the constant energy loss model. The triangular event plane correlations show a different picture, though. While the alignment is maintained for the low decoupling temperature $T_{\text{d}}=$120\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}, deviations for non-central collisions are found if this temperature is increased to $T_{\text{d}}=$160\text{,}\mathrm{M}\mathrm{e}\mathrm{V}. This effect is not observed for the highest collision energy, but it does not seem to depend on the collision energy otherwise. Furthermore, this result is consistent with the large suppression of $v_{3}$ observed previously for $T_{\text{d}}=$160\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}.

Turning to the M&T Langevin parametrization, the same analysis for the event plane angles is shown in Fig. 43. The elliptic flow event plane angles are shown to be heavily aligned with that of charged particles, in the same manner as the previous results for the constant energy loss model. In this case, however, the result for the triangular flow event plane angles is lower even for $T_{\text{d}}=$120\text{,}\mathrm{M}\mathrm{e}\mathrm{V}$$. For the larger decoupling temperature the correlation behaves similarly to the results for the constant energy loss in which the angles are less correlated for non-central collisions. The collision energy dependence is different in this case and the curves do not overlap in the same way we observe in Fig. 42. Looking back at Figs. 38, 39, and 41 we notice that this parametrization leads to the lowest values of $v_{3}$ and the low event plane angle correlation is therefore consistent with that observation.

It is interesting to check if the decorrelation keeps increasing when considering event plane angles of higher order flow harmonics. The plots in Fig. 44 summarize the results for the fourth order flow harmonics for both transport models studied. It is noticeable how these event plane angles are shown to be much less correlated with that of charged particles in comparison to the results found using lower order event plane angles shown before. In addition, the same general trend is maintained for ${{\text{D}}^{0}}$ mesons: the results are less correlated for the highest decoupling temperature and for the Langevin parametrization in contrast to the constant energy loss model. Results for ${{\text{B}}^{0}}$ mesons however do not discriminate the transport models in the same manner and there seems to be no collision energy dependence as well.

Our previous results show the presence of linear scaling between $v_{n}^{\text{(soft)}}$ and $v_{n}^{\text{(heavy)}}$. Considering that the former is also related to the eccentricities in the initial conditions, it is expected that this would also hold for the heavy sector. In Fig. 45 we show the results for the centrality dependence of the Pearson coefficients $Q_{2}$ and $Q_{3}$ Betz et al. (2017) associated with ${{\text{D}}^{0}}$ mesons computed using the constant energy loss model and the M&T parametrization of the Langevin model in PbPb collisions at $\sqrt{s_{\text{NN}}}=$5.02\text{,}\mathrm{TeV}$$ for $T_{\text{d}}=120$–$160\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$. We find that the centrality dependence of these quantities is generally weak. Concerning elliptic flow, we see that both models display large values of $Q_{2}$, with the constant energy loss model giving larger values than the Langevin description. In this case, the results do not vary appreciably when the decoupling temperature is varied in the range 120–$160\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$. However, while $Q_{3}$ displays large values $>0.8$ for the constant energy loss model, the same cannot be said about the Langevin result where $Q_{3}$ drops to $\sim 0.6$. This shows that in general $v_{n}^{\text{(heavy)}}$ is more strongly correlated to $v_{n}^{\text{(soft)}}$ in the constant energy loss model than in the Langevin description. This may be expected due to the noise term in the Langevin description. Moreover, we see that the results do not change significantly with the decoupling temperature, unless in the case of more peripheral collisions.

IV The effect of initial conditions: Trento vs. MCKLN

As discussed in Sec. II.1, while MCKLN and Trento can both reproduce $v_{2}\left\{2\right\}$ of the soft sector well, Trento generally does best when one considers multiple beam energies when it comes to the ratio $v_{2}\{4\}/v_{2}\{2\}$ (specifically at LHC energies, while for RHIC MCKLN does quite well). This implies that MCKLN initial conditions can capture the mean of the $v_{2}$ distribution well but they do not have a wide enough $v_{2}$ distribution compared to experimental data. On the other hand, Trento can capture both the mean and the width of the distribution well.

The question still remains how these differences translate into the heavy flavor sector. Previous studies have looked at the influence of the choice of initial conditions on $\langle v_{2}\rangle$ for smoothed, averaged initial conditions in Cao et al. (2013) and on $v_{2}\{2\}$ Noronha-Hostler et al. (2016c); Xu et al. (2019) defined using event-by-event initial conditions but the influence on the actual $v_{2}$ fluctuations has not yet been considered. However, as shown in Sec. II.1, this is precisely the sector where we expect the largest differences in the initial conditions.

In II.5 it was explained how that the parameters are tuned using the calculation of $R_{\text{AA}}$ in central collisions. Thus, there are slight differences in the parameters for Trento vs. MCKLN initial conditions. Finally, we note that in the figures below coalescence has not been included, which will be done in Sec.  V, so we do not anticipate a perfect match to experimental data. Rather, we focus here on the qualitative aspects of the results.

IV.1 Nuclear modification factor

Based on the previous study in Ref. Noronha-Hostler et al. (2016c), we did not expect large differences in $R_{\text{AA}}$ when comparing different initial conditions and, indeed, we find that to be true. We note, however, in  Noronha-Hostler et al. (2016c) only the energy loss scenario was considered, which we find to have essentially no dependence on the initial conditions. However, the Langevin model does see a slight enhancement in $R_{\text{AA}}$ at low $p_{\text{T}}$ when Trento initial conditions are considered. Since there is a tuning parameter to the $R_{\text{AA}}$, and Trento initial conditions do exhibit a slightly different behavior in the Langevin model, we anticipate that the Langevin model may see differences in the $v_{n}$’s as well.

IV.2 Two-particle $v_{n}$ cumulants

Fig. 47 shows $v_{2}\{2\}(p_{\text{T}})$ at a fixed decoupling temperature. In general, we find that MCKLN initial conditions produce a larger $v_{2}\{2\}$ than Trento initial conditions and that this effect is clearest in the low $p_{\text{T}}$ sector.

Given that MCKLN initial conditions have a slightly larger $\varepsilon_{2}$ than that found in Trento, one might conclude that this is why we see an enhancement in $v_{2}\{2\}(p_{\text{T}})$ for MCKLN initial conditions. We would like to point out, however, that the origin of this result is more complicated than that. In Fig. 1 we demonstrated that by varying the hydrodynamic parameters we were able to find similar $v_{n}\{2\}$ for MCKLN and Trento initial conditions in the soft sector. This affect is achieved by including a larger $\eta/s$ in the hydrodynamic backgrounds for the mckln initial conditions compared to those from TRENTO. Thus, initial conditions with similar results in the soft sector can lead to different results in the heavy flavor sector.

In Fig. 48 $v_{3}\{2\}$ results for Trento and MCKLN initial conditions are shown. While the energy loss scenario for heavy quarks shows no sensitivity to the initial conditions, we find that Trento enhances $v_{3}\{2\}$ for Langevin.

IV.3 Elliptic flow from multi-particle cumulants

While there are certain caveats to the 2-particle $v_{2}$ correlation when it comes to differences in the initial conditions, we do not expect those to strongly affect the ratio between the 4-particle and the 2-particle correlation. Thus, differences that arise in $v_{2}\{4\}/v_{2}\{2\}$ from the initial conditions may be very useful to constrain initial conditions (as was previously done in the soft sector in Giacalone et al. (2017)). In Fig. 49 a comparison between the integrated $v_{2}\{4\}/v_{2}\{2\}$ for MCKLN vs. Trento initial conditions is shown for our two best fitting dynamical models. As previously discussed in Sec. III.2.2, central collisions and peripheral collisions appear to be the best testing beds for different model parameters. This holds true for the initial conditions as well. In fact, central collisions are especially interesting because $v_{2}\{4\}/v_{2}\{2\}>1$ for MCKLN initial conditions (regardless of all dynamical parameters) and $v_{2}\{4\}/v_{2}\{2\}<1$ for Trento initial conditions. Thus, we find that one must incorporate the correct initial conditions first before being able to determine systematic differences from the dynamics of heavy quarks. Because of such a dramatic effect, we strongly encourage experimentalists to investigate these observables in upcoming runs.

V The effect of coalescence

In this section we explore the effect of heavy-light quark coalescence on D meson production. We use a hybrid coalescence plus fragmentation model Cao et al. (2015b) based on the widely used “instantaneous approach” of coalescence Dover et al. (1991); Lin and Molnar (2003); Fries et al. (2003b); Greco et al. (2004); Oh et al. (2009). Within this approach, the probability distribution ${\rm P}_{\rm coal}[q,Q\rightarrow M]$ that a heavy meson of momentum ${\bf p}_{M}$ is formed by coalescence of a heavy quark with momentum ${\bf p}_{Q}$ and a light quark of momentum ${\bf p}_{q}$ is given by

[TABLE]

where $N$ is a global normalization factor discussed below, $n_{q}$ is the momentum distribution of the light quarks at the time of hadronization, and $f_{M}$ is the probability density obtained from the projection of the two quark states onto the meson state. The different meson states are obtained from a simple harmonic oscillator model. The quarks being relativistic, the projection is performed in their center-of-mass frame following a Lorentz boost of their phase space coordinates from the global frame to the center-of-mass frame. Assuming that the light quarks have a uniform spatial distribution in the medium cell where the heavy quark is located, one can average the probability density over the spatial coordinates, which then gives

[TABLE]

where $g_{M}$ is the usual color-spin-isospin statistical factor for the two spin 1/2 quarks to form a color neutral meson $M$, and $h_{M}$ are new “thermal” factors discussed below. The width $\sigma=1/\sqrt{\mu\omega}$ is given by the angular frequency of the harmonic oscillator $\omega$ and the reduced mass of the two quark system $\mu=m_{q}m_{Q}/(m_{q}+m_{Q})$. The mass of the light quarks are assumed to be their constituent masses666If the effective or quasi-particle masses Levai and Heinz (1998); Plumari et al. (2011) for the light quarks is taken to be in the range $m_{u,d,s}\sim 300$-600 MeV, no significant differences in the final values of ${\rm P}_{\rm coal}$ are found. $m_{u,d}=300$ MeV and $m_{s}=460$ MeV Borka Jovanovic et al. (2010) to take into account the non-perturbative effects of QCD near the cross-over transition temperature $T_{c}$, while the heavy quark masses are taken to $m_{c}=1.27$ GeV and $m_{b}=4.19$ GeV. The value of the angular frequency is evaluated to be $\sim 0.3$ GeV for ${{\text{D}}^{0}}$, ${{\text{D}}^{+}}$, $\text{D}_{s}$ and ${\text{B}}^{-}$ mesons777Note that for ${{\text{B}}^{0}}$ and ${\text{B}}_{s}^{-}$ we obtain $\omega\sim 0.5$ GeV and $\sim 0.6$ GeV, respectively, which is quite different from the other meson ground states. from their vacuum charge radii obtained within the light cone model Hwang (2002). Nevertheless, it is known that at $T\sim T_{c}$ the charge radii must be (much) larger than in the vacuum. Motivated by a (limited) comparison to the case of quarkonia in which the radii of the $J/\psi$ and $\Upsilon(1S)$ are evaluated from lattice QCD potentials to be $\sim 3$ times larger at $T_{c}$ than in the vacuum, we set $\omega=0.1$ GeV ($\omega$ being inversely proportional to the radius) and assume this value to be valid for all the hadrons considered within this basic model. The relative momentum of the two quarks including relativistic corrections

[TABLE]

is defined via the momenta in the center-of-mass frame (denoted with primed coordinates). The momentum density distribution of the light quarks $n_{q}$ is assumed to be thermal in the local rest frame of the medium cell considered being given by the Fermi-Dirac distribution

[TABLE]

where ${\bf p}_{q}^{\rm cell}$ is the light quark momentum in the local rest frame of the medium cell, $g_{q}=6$ is the statistical factor that takes into account the spin and color degeneracy of the light quarks, $p_{q}.u$ is the 4-product between the light quark 4-momentum in the global frame and the fluid 4-velocity $u=(\gamma_{u},\gamma_{u}{\bf u})$ with $\gamma_{u}=\frac{1}{\sqrt{1-{\bf u}^{2}}}$. Thanks to this fluid velocity dependence of $n_{q}$, the derived coalescence probabilities will depend not only on the heavy quark momentum but also on the local flow and on the angle between them.

By extension, the probability distribution ${\rm P}_{\rm coal}[q_{1},q_{2},Q\rightarrow B]$ that a heavy baryon of momentum ${\bf p}_{B}$ is formed by coalescence of a heavy quark with momentum ${\bf p}_{Q}$ and two light quarks of momenta ${\bf p}_{q_{1}}$ and ${\bf p}_{q_{2}}$ is given by

[TABLE]

where $N$ is the global normalization factor, $f_{B}$ is the probability density obtained by combining the two light quarks first and then using their center of mass to recombine with the heavy quark,

[TABLE]

where $g_{B}$ is the usual color-spin-isospin degeneracy factor for the three spin 1/2 quarks to form a color neutral baryon $B$, and $h_{B}$ are the “thermal” factors discussed below. Note that the resulting ${\rm P}_{\rm coal}[q_{1},q_{2},Q\rightarrow B]$ is perfectly identical if we combine the two light quarks first or if we combine the heavy quark and one of the light quarks first. The widths $\sigma_{1}=1/\sqrt{\mu_{1}\omega}$ and $\sigma_{2}=1/\sqrt{\mu_{2}\omega}$ are given by the angular frequency of the harmonic oscillator $\omega$, the reduced mass of the two light quark system $\mu_{1}=m_{q_{1}}m_{q_{2}}/(m_{q_{1}}+m_{q_{2}})$, and the reduced mass of the system formed by the two light quarks center-of-mass and the heavy quark $\mu_{2}=(m_{q_{1}}+m_{q_{2}})m_{Q}/(m_{q_{1}}+m_{q_{2}}+m_{Q})$. The relative momenta with relativistic corrections are

[TABLE]

and

[TABLE]

defined888The relative momenta, as defined in Cao et al. (2015b), i.e. with $E_{q}^{\prime}$ replacing $m_{q}$ in the definitions, have the disadvantage of breaking the “combining order symmetry”, i.e. ${\rm P}_{\rm coal}[q_{1},q_{2},Q\rightarrow B]$ is different if we combine the two light quarks first or the heavy quark and a light quark. It is then difficult to justify why one type of ordering should be better than the other. via the momenta in the baryon center-of-mass frame (denoted with primed coordinates). Similarly to the mesons, the momentum density distributions of the light quarks $n_{q_{1}}$ and $n_{q_{2}}$ are assumed to be thermal in the local rest frame of the considered medium cell and given by the Fermi-Dirac distribution (30).

For the charmed mesons (baryons) we consider all the symmetric $J^{P}=0^{-}$ and $1^{-}$ states ($J^{P}=1/2^{+}$ and $3/2^{+}$ states), i.e. the ${{\text{D}}^{0}}$, ${{\text{D}}^{+}}$, ${\text{D}}_{s}^{+}$ ($g_{M}=1/36$ each), ${\text{D}}^{*0}$, ${\text{D}}^{*+}$, ${\text{D}}_{s}^{*+}$ ($g_{M}=1/12$ each) mesons and the $\Lambda_{c}^{+}$, $\Sigma_{c}^{0}$, $\Sigma_{c}^{+}$, $\Sigma_{c}^{++}$ ($g_{B}=1/108$ each), $\Sigma_{c}^{*0}$, $\Sigma_{c}^{*+}$, $\Sigma_{c}^{*++}$ ($g_{B}=1/18$ each), $\Xi_{c}^{0}$, $\Xi_{c}^{+}$, $\Xi_{c}^{{}^{\prime}0}$, $\Xi_{c}^{{}^{\prime}+}$ ($g_{B}=1/54$ each), $\Xi_{c}^{*0}$, $\Xi_{c}^{*+}$ ($g_{B}=1/27$ each), $\Omega_{c}^{0}$ ($g_{B}=1/108$), $\Omega_{c}^{*0}$ ($g_{B}=1/54$) baryons. The anti-symmetric states cannot be included as the considered Wigner distribution (28) is symmetric, being thus usually neglected in coalescence models.

The global normalization factor $N$ is usually chosen such that ${\rm P}_{\rm coal}[c\rightarrow{\rm any\,hadron}]({\bf p}_{Q}=0)=1$ assuming that a quark with zero momentum does not hadronize via fragmentation but via coalescence only999See a criticism of this assumption in Song et al. (2016).. In this work we extend this idea by assuming instead that ${\rm P}_{\rm coal}[c\rightarrow{\rm any\,hadron}]({\bf p}_{Q},{\bf u})=1$ when the heavy quark velocity vector is equal to the one from the local medium flow, i.e. when the heavy quark is moving together with the light quarks surrounding it. In this way we avoid the low ${\bf p}_{Q}$ “saturation” of ${\rm P}_{\rm coal}$ when $|{\bf u}|>0$ obtained in Cao et al. (2015b). This saturation seems unjustified to the present authors: a heavy quark with a velocity much lower than the typical velocities of the surrounding light quarks may not coalesce with unity probability.

Recently, the prompt heavy hadron ratios ${{\text{D}}^{+}}/{{\text{D}}^{0}}$, ${\text{D}}^{*+}/{{\text{D}}^{0}}$, ${\text{D}}_{s}^{+}/{{\text{D}}^{0}}$, ${\text{D}}_{s}^{+}/{{\text{D}}^{+}}$, and $\Lambda_{c}^{+}/{{\text{D}}^{0}}$ have been measured in heavy ion collisions and compared to proton-proton collisions Acharya et al. (2018d); Radhakrishnan (2018); Vermunt (2019); Acharya et al. (2018e). They provide important information about the hadronization process in heavy ion collisions and can serve as a benchmark to test and calibrate coalescence models. First, the low and intermediate $p_{\text{T}}$ enhancement of the $\Lambda_{c}^{+}/{{\text{D}}^{0}}$ ratio tends to confirm the presence of heavy-light quarks coalescence in heavy ion collisions when $p_{\text{T}}<7$ GeV. The $p_{\text{T}}$ and collision system independence of the ${{\text{D}}^{+}}/{{\text{D}}^{0}}$ and ${\text{D}}^{*+}/{{\text{D}}^{0}}$ ratios show a certain universality in the way the different hadronization processes distribute the quarks between the meson states. Finally, the enhancement of ${\text{D}}_{s}^{+}/{{\text{D}}^{+}}$ at low and intermediate $p_{\text{T}}$ confirms the thermal enhancement of strange quark production. Within the “basic” coalescence model (i.e. without the $h_{M}$ and $h_{B}$ factors) the prompt ${{\text{D}}^{+}}/{{\text{D}}^{0}}$ and ${\text{D}}^{*+}/{{\text{D}}^{0}}$ ratios are only based on the color-spin-isospin degeneracy factors, the rest of the probabilities canceling out, e.g.

[TABLE]

where “prompt” means including feed-downs from excited states and “dir.” means the direct production (without feed-downs) that we obtain from the coalescence model. The different ratio values obtained within the “basic” coalescence model are summarized in table 3 and compared to low $p_{\text{T}}$ experimental data in AA collisions. The ${{\text{D}}^{+}}/{{\text{D}}^{0}}$ and ${\text{D}}_{s}^{+}/{{\text{D}}^{+}}$ ratios particularly miss the experimental data. Additionally the $\Lambda_{c}^{+}/{{\text{D}}^{0}}$ ratio fits LHC preliminary data Vermunt (2019) but underestimate RHIC preliminaries Radhakrishnan (2018). To improve the fit to the ${{\text{D}}^{+}}/{{\text{D}}^{0}}$ ratio, one could naturally think of decreasing the angular frequency $\omega$ for the excited states — as their radii should be larger than the ground states — but the ${{\text{D}}^{+}}/{{\text{D}}^{0}}$ ratio is then observed to decrease, which is the opposite of what was expected.

To improve the fit to the ratios and, subsequently, the predictive power of coalescence in our study one can instead take a step aside from the basic model and include by hand a new element. One of the weak points of the basic model is the complete absence of hadron masses101010A consequence of this is the non conservation of the $4^{\text{th}}$ component of the system’s 4-momentum during the coalescence process. We do not intend to address this problem with these new elements. in the formalism whereas it is clear that the formation of excited states with larger masses requires more energy from the combining quark than to form the ground state. Inspired by the thermal model of hadronization Oh et al. (2009), we include the “thermal” factors $h_{H=M,B}=\exp[-(m_{H}-m_{{H}_{0}})/T_{\text{d}}]$, where $H$ is the hadron state for which we want to compute ${\rm P}_{\rm coal}$, $m_{H}$ is its mass, and $m_{{H}_{0}}$ is the mass of the corresponding ground state ${H}_{0}$ built with the same quark content as $H$, e.g. ${{\text{D}}^{0}}$ for ${\text{D}}^{*0}$ or $\Lambda_{c}^{+}$ for $\Sigma_{c}^{+}$ and $\Sigma_{c}^{+*}$. Therefore, for the same quark content, an excited state $H$ has $g_{H}\exp[-m_{H}/T_{\text{d}}]/(g_{{H}_{0}}\exp[-m_{{H}_{0}}/T_{\text{d}}])$ less probability to form than its corresponding ground state ${{H}_{0}}$. With these new factors the ratios are, thus, not anymore based on spin and color considerations alone but they also depend on the hadron masses, leading to a more relevant statistical hierarchy between the different energy states of a combination of quark flavors. Note that thanks to these new factors a natural justification arises for the non consideration of more excited states (such as the $J=2$ mesons and $J=5/2$ baryons) and antisymmetric states (such as the $J^{P}=0^{+}$ mesons and $J^{P}=1/2^{-}$ baryons), as they are now suppressed due to their larger mass. The resulting ratios are summarized in Tab. 3 and fit well the experimental data. Additionally, a larger $\Lambda_{c}^{+}/{{\text{D}}^{0}}$ ratio is obtained with the new factors, giving more confidence in the possibility of fitting both RHIC and LHC final data.

We emphasize again that the motivation for these new thermal factors is purely phenomenological, i.e., to improve the ratio fits and therefore the quality of the coalescence probabilities, and we do not assume that they are theoretically justified from the beginning. The probabilities of coalescence into any hadron and into ${{\text{D}}^{0}}$ mesons and $\Lambda_{c}^{+}$ baryons are shown in Fig. 50 for different medium cell velocities $|{\bf u}|$ and $\theta$ angles between ${\bf u}$ and ${\bf p}_{Q}$. Larger $|{\bf u}|$ and, thus, larger average $|{\bf p}_{q}|$ tend to shift the maximum of probability towards larger $|{\bf p}_{Q}|$. An increasing $\theta$ angle tends to lower the probability amplitudes and restrict the coalescence to heavy quarks with lower momenta. One can finally note that with increasing $|{\bf u}|$ the meson production from coalescence tends to decrease whereas the baryon production increases.

The choice between fragmentation and coalescence in heavy ion collisions is performed for each heavy quark through a simple Monte Carlo procedure using the probabilities of coalescence ${\rm P}_{\rm coal}$ described above and of fragmentation found in Cacciari and Nason (2003) for the rest of the probability $1-{\rm P}_{\rm coal}$. If the coalescence into a ${{\text{D}}^{0}}$ meson is drawn for a heavy quark of momentum ${\bf p}_{Q}$, a light quark of momentum ${\bf p}_{q}$ is generated according to the momentum space density distribution $n_{q}$ in the frame of the medium cell and then boosted back into the global frame where it can possibly combine with the given heavy quark according to the probability density $f_{M}({\bf p}_{q},{\bf p}_{Q})$ using a Monte Carlo procedure. If they do not recombine, another light quark is generated until the meson is formed. The ${{\text{D}}^{0}}$ meson momentum ${\bf p}_{{{\text{D}}^{0}}}$ is finally given by the momentum of the heavy-light quark system ${\bf p}_{Q}+{\bf p}_{q}$.

V.1 Nuclear modification factor

In Fig. 51 we show the effect of coalescence on the nuclear modification factor. The inclusion of coalescence in the M&T Langevin model remarkably improves the fit to the $p_{\text{T}}<5$ GeV data leading to a good description of the entire $p_{\text{T}}$ range. For the constant energy loss model the inclusion of coalescence only amounts to a very small increase in $R_{\text{AA}}$ for $p_{\text{T}}<10$ GeV, which does not fix the large discrepancy at low $p_{\text{T}}$ already observed using only fragmentation. Note however that these results correspond to direct $c\rightarrow{{\text{D}}^{0}}$ production whereas the experimental data in this section correspond to prompt ${{\text{D}}^{0}}$ production. The probabilities of $c\rightarrow{{\text{D}}^{*}}$ being different within the fragmentation and coalescence frameworks (see Tab. 3), the prompt results might be a bit different at low $p_{\text{T}}$. The effect of coalescence on $R_{\text{AA}}$ can be decomposed into a small $R_{\text{AA}}$ shift towards larger values - most probably due to a smaller direct $c\rightarrow{{\text{D}}^{0}}$ production with coalescence than with fragmentation (see Tab. 3) - and more importantly into a $p_{\text{T}}$ shift of the low $p_{\text{T}}$ $R_{\text{AA}}$ values of around $\sim 2$ GeV towards higher $p_{\text{T}}$. The latter occurs because fragmentation produces a hadron with less momentum than the heavy quark while in coalescence there is a momentum gain from the light quark “thermal” momentum and mass. With only fragmentation the low $p_{\text{T}}$ $R_{\text{AA}}$ for the energy loss model becomes almost $p_{\text{T}}$ independent and, thus, the $p_{\text{T}}$ shift cannot be observed. In contrast, with only fragmentation the low $p_{\text{T}}$ $R_{\text{AA}}$ of the Langevin model exhibits a large decrease in $R_{\text{AA}}$ such that the $p_{\text{T}}$ shift is visible and has a strong impact.

The difference between the $R_{\text{AA}}$ obtained with and without the thermal factors in the coalescence model is shown in Fig.  52. The extra thermal factors logically increase the low $p_{\text{T}}$ $R_{\text{AA}}$ due to a larger ${{\text{D}}^{0}}$ formation probability (see Tab. 3).

V.2 Elliptic and triangular flow coefficients with two-particle cumulants

Regarding $v_{2}$ and $v_{3}$, as shown in Fig. 53 we do not find a single model that can quantitatively capture all the experimental data though coalescence generally improves the description. The M&T Langevin model with coalescence has the best fit for $p_{\text{T}}\lesssim 4$ GeV but underpredicts the data at higher $p_{\text{T}}$, whereas the constant energy loss model with coalescence performs best from $p_{\text{T}}\gtrsim 5$ GeV but underpredicts the data at lower $p_{\text{T}}$. Similarly to what we observed for $R_{\text{AA}}$, coalescence tends to shift the low $p_{\text{T}}$ peaks in $v_{n}(p_{\text{T}})$ to higher $p_{\text{T}}$. Additionally, coalescence also tends to non-negligibly increase the low $p_{\text{T}}$ $v_{3}$, i.e., when combining the light quarks one communicates part of the medium’s triangular flow to the heavy quarks.

The difference between the $v_{2}$ obtained with and without the thermal factors in the coalescence model is shown in Fig.  54. The extra thermal factors are observed to change the $v_{2}$ except at $p_{\text{T}}\approx 1.5$ GeV where the coalescence probability is maximal (see Fig.50) and the fragmentation probability close to 0. The $v_{2}$ obtained through the coalescence and fragmentation processes being different, this effect can be explained by the different relative proportions of ${{\text{D}}^{0}}$ coming from coalescence and fragmentation in the two configurations of the coalescence model. More generally, these non-negligible variations of the observables with the coalescence model configuration underline the key role played by the $c\rightarrow{\text{D}}$ coalescence probability which value depends on various features of the coalescence models, e.g. the angular frequency $\omega$ and the considered set of hadron states111111Set of hadron states chosen differently for instance in Oh et al. (2009) and Cao et al. (2015b)..

V.3 Heavy mesons and all charged particles flow correlation

With the addition of coalescence one can explore the low $p_{\text{T}}$ soft-hard flow correlations with more confidence. As seen in Fig. 55, coalescence increases the linear correlation between the heavy meson anisotropies $v_{n}^{\rm(heavy)}$ and the all charged particles elliptic flow $v_{n}^{\rm(soft)}$ proportionally to the gain in $p_{\text{T}}$ integrated $v_{n}^{\rm(heavy)}$ due to coalescence. The increase is then more important for the M&T Langevin model in the $p_{\text{T}}$ range 0–8 GeV.

Fig. 56 shows how coalescence affects the Pearson coefficients associated with elliptic and triangular anisotropies, $Q_{2}$ and $Q_{3}$, computed using the M&T Langevin model. One can see that the correlation weakens with the addition of coalescence, in contrast to what was observed in Fig. 55. Coalescence tends to produce a somewhat more dispersed correlation on an event-by-event basis between the magnitudes $v_{n}^{\rm(heavy)}$ and $v_{n}^{\rm(soft)}$, and between the $\psi_{n}^{\rm(heavy)}$ and $\psi_{n}^{\rm(soft)}$ distributions leading to a slight event plane decorrelation.

V.4 Elliptic flow from multi-particle cumulants

We show in Fig. 57 the influence of the hadronization process on the low $p_{\text{T}}$ integrated $v_{2}\{4\}/v_{2}\{2\}$ cumulant ratio as a function of centrality. The addition of coalescence is observed to generally decrease the low $p_{\text{T}}$ integrated cumulant ratio. Remembering Eq. (23), it is indeed consistent for the coalescence mechanism to reduce the cumulant ratio as it enhances the heavy quark anisotropy fluctuations with respect to the soft fluctuations obtained via the heavy-light quark combination. As shown on the left panel, the extra thermal factors lead to a larger decrease of the cumulant ratio that can be explained by a larger proportion of ${{\text{D}}^{0}}$ coming from coalescence than in the “basic” configuration.

VI Conclusion

In summary, in this paper we have presented the details behind DAB-MOD, a new heavy flavor model code that allows for a modular description of the dynamics of heavy quarks in heavy-ion collisions. The code can be run on top of any event-by-event relativistic hydrodynamic background and it samples the initial heavy quark distribution according to pQCD FONLL calculations. Heavy quark evolution is done either by solving relativistic Langevin equations or employing parametrized energy loss models (which include energy loss fluctuations). Finally, the heavy quarks can either be fragmented or coalesce to produce heavy flavor mesons (heavy meson decay is implemented via Pythia 8 to produce semileptonic channels). To improve the fits to the recently observed heavy hadron ratios at low $p_{T}$ and therefore the predictiveness of our coalescence model, we proposed an “empirical” inclusion to the usual probabilities of some new thermal-like factors (based on hadron masses) between hadrons of identical quark content. They have the additional asset of giving a natural explanation for the non-consideration of anti-symmetric and higher spin states (because of their larger mass which then suppresses their probabilities).

In this paper we have thoroughly checked the influence of a wide range of assumptions/parameters involved in the modeling of the heavy flavor sector in heavy ion collisions such as the inclusion or not of energy loss fluctuations, energy loss models vs. Langevin descriptions, heavy quark decoupling temperatures (which define when the heavy quarks decouple from the medium), heavy flavor meson vs. electrons/muon observables, choice of initial conditions for the background hydrodynamic evolution, as well as coalescence vs. fragmentation hadronization mechanisms. Similarly to Betz and Gyulassy (2014), we find that energy loss fluctuations play a significant role in $v_{n}(p_{\text{T}})$ calculations at $p_{\text{T}}<10$ GeV. In general, either a no energy loss fluctuation scenario (which we admit to be unrealistic) or a Gaussian distribution fare best compared to experimental data. Our best fitting results at low $p_{\text{T}}$ stem from the Moore and Teaney inspired spatial diffusion coefficient within the Langevin formalism with the inclusion of coalescence. At high $p_{\text{T}}$ (roughly $p_{\text{T}}>5$ GeV), the energy loss model with a constant energy loss works the best when also coupled to coalescence. Thanks to the momentum brought by the light quark to the final meson, the addition of coalescence is observed to shift the low $p_{T}$ variations or lumps of the $R_{\text{AA}}$ and $v_{n}$’s of around $\sim 2$ GeV towards higher $p_{\text{T}}$.

Unlike in Andres et al. (2019) where the effects of the initialization time were explicitly studied, in this work we kept a fixed heavy quark initialization time, $\tau_{0}=0.6$ fm, but we note that that is compatible with the best fit in Andres et al. (2019). However, we remark that our two different setups for the initial conditions (MCKLN vs. Trento) used different equations of state (the results computed using MCKLN employed a now outdated equation of state from Huovinen and Petreczky (2010) while the runs using Trento employed a newer one which contains the most up-to-date Particle Data booklet resonances from  Alba et al. (2018)). In general, these correspond to different initialization temperatures $T_{0}$ and, thus, the MCKLN events had a lower $T_{0}$ (even though $\tau_{0}$ was fixed) than the Trento events. Due to the significant run time that these different backgrounds take, we did not check this difference explicitly in this paper. However, we do plan in a later paper to investigate possible equation of state effects, especially its underlying assumptions such as whether charm quarks are thermalized or not Borsanyi et al. (2016); Alba et al. (2018); Noronha-Hostler and Ratti (2018). Additionally, it would be interesting to check the effects of a pre-equilibrium stage as in Das et al. (2017). Keeping in mind these limitations, the type of initial conditions is observed to have a little impact on the heavy flavor $R_{\text{AA}}$ and $v_{n}\{2\}$’s within our model, with Trento leading to a slightly lower $v_{2}\{2\}$.

In this paper we investigated higher order flow coefficients involving heavy flavor. Similarly to the experimental data, the $v_{3}\{2\}$ is observed to be mostly centrality independent. An intriguing result was found concerning $v_{4}\{2\}$ of heavy flavor leptons in $\sqrt{s_{\text{NN}}}=$200\text{,}\mathrm{GeV} AuAu and $\sqrt{s_{\text{NN}}}=$5.02\text{\,}\mathrm{TeV} PbPb collisions: if the decoupling temperature parameter, $T_{\text{d}}$, is too large ($T_{\text{d}}=160$ MeV) then $v_{4}\{2\}$ becomes negative. This occurs when the heavy quark does not have enough time to interact with the medium, leading to an anti-correlation of the heavy flavor and bulk event planes. In Aaboud et al. (2018) $v_{4}\{2\}$ was measured by the ATLAS collaboration and, with the current error bars, it is not yet clear what the overall sign of $v_{4}\{2\}$ is. However, for some centrality classes it does appear that it could be positive and in others negative. With future upgrades to the experiments it may be eventually possible to use $v_{4}\{2\}$ as a way to constrain the decoupling temperature and obtain an estimate of how long heavy quarks remain coupled to the expanding medium. This would be very useful on the theoretical front since $T_{\text{d}}$ produces a large systematic uncertainty in our calculations.

We also make predictions for multi-particle cumulants that have not yet been measured experimentally in the heavy flavor sector. We find that the ratio $v_{2}\{4\}/v_{2}\{2\}$ is different from the soft sector, and especially interesting in the most central collisions (0–10%) and in peripheral collisions ($50-60\%$) because not only it is predominately sensitive to the initial conditions (which affects its overall magnitude) but it also has a different behavior across $p_{\text{T}}$ depending on the choice made to model heavy flavor evolution throughout the medium. Finally, $v_{2}\{4\}/v_{2}\{2\}$ has a rather small but non-trivial beam energy dependence so it would be especially interesting to see model comparisons to data from sPHENIX in a few years and also to data acquired at the top LHC energies.

Additionally, we perform event-shape engineering calculations in the heavy flavor sector, which is the theory analog of what was experimentally done in Acharya et al. (2018a). One caveat is that the results depend on the number of bins chosen within a centrality class (as well as the method of centrality binning). Thus, this must be set by an individual experimental collaboration depending on their available statistics before direct theory vs. experimental comparisons can be made. In fact, we suggest the experiments to always publish the corresponding integrated $v_{n}\{2\}$ of all charged particles for their event-shape engineered bins so more precise comparisons can be made. However, even with these caveats in mind, we clearly see linear correlations between soft and heavy $v_{n}$’s for both elliptical and triangular flow. We have observed as well the non-trivial correlations and decorrelations of the soft and heavy event planes $\psi_{n}$’s across centralities, colliding energies, heavy quark masses and models (e.g. the energy loss models leading to more correlations than Langevin dynamics).

The development of DAB-MOD allows for a systematic comparison of a variety of approaches to heavy quark dynamics using state-of-the-art hydrodynamic backgrounds that can describe the behavior of multi-particle cumulants in the soft sector of heavy-ion collisions. The predictions-data comparisons made throughout this paper tend to show that it is necessary to include more phenomena to the simulation, such as the initial shadowing or the final hadronic re-scattering, to improve its predictions (especially for the $v_{n}\{2\}$’s). With the upcoming era of sPHENIX the creation of DAB-MOD will allow for many future heavy flavor studies and constraints on our knowledge of heavy flavor dynamics. This will allow for many future comparisons involving different collisions systems, energies, and sizes (which will hopefully further constrain heavy flavor theoretical modeling) and become an invaluable tool to study heavy flavor in the upcoming RHIC and LHC runs.

Acknowledgments

We wish to thank Prof. Pol-Bernard Gossiaux and Shanshan Cao for their fruitful discussions and help. The authors thank Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) for support. R.K. is supported by the Region Pays de la Loire (France) under contract No. 2015-08473. C.A.G.P. is supported by the NSFC under grant No. 11521064, MOST of China under Project No. 2014CB845404. J.N.H. acknowledges the support of the Alfred P. Sloan Foundation, support from the US-DOE Nuclear Science Grant No. DE-SC0019175, and the Office of Advanced Research Computing (OARC) at Rutgers, The State University of New Jersey for providing access to the Amarel cluster and associated research computing resources that have contributed to the results reported here. J.N. is partially supported by CNPq grant 306795/2017-5 and FAPESP grant 2017/05685-2.