Quenching and flow of charm and bottom quarks via semi-leptonic decay of D and B mesons in Pb+Pb collisions at the LHC

A wealth of evidence indicates that a color deconfined quantum chromodynamic (QCD) matter, called quark-gluon plasma (QGP), is produced in high-energy nucleus-nucleus collisions performed at the BNL Relativistic Heavy-Ion Collider and CERN Large Hadron Collider (LHC) [1–6]. Heavy quarks, including charm and bottom quarks, have served as a successful probe of the QGP properties in these energetic collisions [7–9]. As their masses are much larger than $ \Lambda_\text{QCD} $ and the temperature scale of the QGP created in current RHIC and LHC experiments, they are mainly produced from the initial hard collisions between nuclei, whose cross sections can be calculated within the perturbative QCD (pQCD) framework. After being produced, heavy quarks travel through and interact with the QGP with their flavors conserved before forming heavy flavor hadrons upon exiting the QGP. At high transverse momentum ($ p_\text{T} $), heavy quarks are expected to lose energy inside the QGP via elastic and inelastic scatterings with the medium constituents, similar to the energy loss of light quarks and gluons but with mass effect introduced [10–12]. At low $ p_\text{T} $, heavy quarks are expected to diffuse inside the QGP and evolve towards thermal equilibrium with the QGP [13–16]. As a consequence, heavy quarks can inherit a significant amount of flow from the expanding medium, which is usually anisotropic in heavy-ion collisions. At intermediate $ p_\text{T} $, the hadronization process is crucial in forming the final heavy flavor hadron observables, including their spectra, collective flow coefficients, and chemical components [17–21]. Experimental measurements on heavy flavor mesons and their decay products, such as non-prompt D, $ J/\psi $, and leptons, have shown significant suppression of their yields [22–26] in nucleus-nucleus (A-A) collisions with respect to proton-proton (p-p) collisions. Strong anisotropic flow coefficients have also been observed [23, 27–29]. The combined study of heavy flavor particle yield modification and anisotropic flow coefficient can provide valuable constraints on the interaction strength (diffusion coefficient) and dynamics (contributions from elastic and inelastic processes) between heavy quarks and the QGP, as well as various thermodynamic properties of the QGP (e.g., viscosity and equation of state) [11, 14, 16, 30–57].

The yield suppression of heavy flavor particles can be quantified by the nuclear modification factor, $R_{AA}(p_\text{T}) = \dfrac{{\rm d}N_{AA}/{\rm d}p_\text{T}}{\langle N_{\rm coll}\rangle {\rm d}N_{pp}/{\rm d}p_\text{T}}$, where $ \langle N_{\rm coll} \rangle $ is the average number of binary nucleon-nucleon collisions per nucleus-nucleus collision for a given centrality class, which is usually calculated using the Glauber model [58]. If there are neither cold nor hot nuclear matter effects, $ R_{AA} $ should be unity, which has been verified by measurements of vector boson productions at the RHIC and LHC. Contrarily, clear suppression has been observed for the yield of high $ p_\text{T} $ hadrons and jets in relativistic heavy-ion collisions, which is stronger in more central collisions [59–63]. This can be well understood by the energy loss experienced by hard quarks and gluons inside the QGP, namely, jet quenching [64–70]. For the charm sector, the measurements of prompt D mesons show suppression similar to light flavor hadrons [22, 24, 25]; meanwhile, for the bottom sector, the measurements of B mesons, B-decayed D, and $ J/\psi $ also show substantial suppression, though less suppression than prompt D mesons at not very high $ p_\text{T} $ due to the mass effect in heavy quark energy loss [71–73].

The anisotropic flow of particles produced in high-energy nuclear collisions is quantified by performing Fourier expansion of the azimuthal angle distribution of particles as ${{\rm d}N}/{{\rm d}\phi} \propto 1 + 2 \sum_n v_n \cos[n(\phi - \Phi_n)]$, where $ v_n $ is n-th order anisotropic flow coefficient and $ \Phi_n $ is the corresponding event plane angle. The average almond shape of the collision zone in the transverse plane at non-zero impact parameter leads to the elliptic flow coefficient $ v_2 $. Meanwhile, event-by-event quantum fluctuations of nucleonic and sub-nucleonic degrees of freedom contribute to higher-order odd harmonic components. They are also the origin of anisotropic flow in ultra-central collisions [74, 75]. Heavy-ion experiments at the RHIC and LHC have shown strong anisotropic flows for charged and identified hadrons produced from the bulk matter [76–81], which have been successfully explained by relativistic hydrodynamics models [82–86]. The small specific viscosity extracted from these models suggests that the QCD matter produced in relativistic heavy-ion collisions is a strongly-coupled QGP that behaves like a perfect liquid. The elliptic flow of heavy flavor hadrons depends on the combination of several factors. At high $ p_\text{T} $, the heavy flavor flow is mainly driven by the anisotropic energy loss of heavy quarks through different path lengths in different propagation directions. At low $ p_\text{T} $, their flow depends on not only the interaction strength between heavy quarks and the QGP, which determines the degree of thermalization of heavy quarks with the medium background, but also the hadronization process, which further transfers the medium collectivity to the heavy flavor hadrons. Experimental data show that D mesons and light flavor hadrons share similar $ v_2 $ at low $ p_\text{T} $, indicating the thermal behavior of heavy quarks at low $ p_\text{T} $ [23, 27–29].

In this work, we focused on heavy flavor leptons, which are predominantly produced from the semi-leptonic decay of heavy flavor hadrons. In particular, we studied the separate contributions from D and B mesons to the suppression and flow of heavy flavor leptons, as recently measured by the ATLAS Collaboration [87, 88]. This allows us to study the mass dependence of quark interaction with a color deconfined medium. The remainder of this paper is organized as follows. In Sec. II, we give a brief introduction to our theoretical framework on heavy flavor production and evolution in heavy-ion collisions. In Sec. III, we present our numerical results on the nuclear modification and elliptic flow of heavy flavor electrons and muons produced in Pb+Pb collisions at $ \sqrt{s_{NN}} = 5.02 $ TeV. Finally, we summarize our work in Sec. IV.

In this work, the evolution of heavy quarks inside the QGP was simulated via our improved Langevin approach that solves the following equation [14]:

Here, both quasi-elastic scattering and medium-induced gluon bremsstrahlung processes are included. The first and second terms on the right hand side denote the drag force and thermal random force experienced by heavy quarks, respectively, when they scatter with the QGP, while the third term represents the recoil force $\vec{f}_g = {\rm d}\vec{p}_g/{\rm d}t$ exerted on heavy quarks when they emit medium-induced gluons with momentum $ \vec{p}_g $.

We assume the thermal random force $ \vec{\xi} $ to be white noise and satisfy the following correlation:

where κ is the momentum space diffusion coefficient of heavy quarks, characterizing the strength of the thermal force. This can be related to the drag coefficient $ \eta_\text{D} $ via the fluctuation dissipation theorem $ \eta_\text{D}(p)=\kappa/(2TE) $. The spatial diffusion coefficient is then given by $D_\text{s} = T/[M\eta_\text{D}(0)]=2T^{2}/\kappa$. By convention, we choose $ D_\text{s} $ as the model parameter in our Langevin approach.

For inelastic scattering, the medium-induced gluon radiation probability during a time interval $ \Delta t $ can be calculated as

where x and $ k_\perp $ are the energy fraction and transverse momentum of the bremsstrahlung gluon with respect to the heavy quark, respectively. Note that to interpret the average number of radiated gluons $ \langle N_{g}(t,\Delta t)\rangle $ as probability, we choose $ \Delta t $ to be sufficiently small such that $ \langle N_{g}(t,\Delta t)\rangle <1 $. In this study, the medium-induced gluon spectrum was taken from the higher-twist (HT) energy loss formalism [89–92]:

where the strong coupling strength $ \alpha_\text{s} = 4\pi /[9 \text{ln} (k_{\perp}^2/\Lambda^2)] $, with $ k_{\perp} $ being the transverse momentum of the radiated gluon in our calculation. $ \hat q $ is the (gluon) jet transport coefficient [93–95], which relates to κ via $ \hat q = 2\kappa C_{A}/C_{F} $ with $ C_A $ and $ C_F $ being color factors of gluon and quark. $ t_i $ is the initial time of forming the current gluon emission (or the time of the previous emission), and $ \tau _{f}=2Ex(1-x)/ (k_{\perp}^{2}+x^{2}M^{2}) $ is the average formation time of the gluon with E and M being the energy and mass of heavy quarks, respectively. We impose a cutoff for the radiated gluon energy ($ E_g = xE > \pi T $) to balance the gluon emission and absorption processes around the thermal scale. With this implementation, heavy quarks are able to approach thermal equilibrium at large time, though the extracted temperature parameter is approximately 5%−10% lower than the medium temperature [14]. Note that our improved Langevin approach has one free parameter, the dimensionless parameter $ D_\text{s} (2\pi T) $, which characterizes the strength of heavy quark interaction with the QGP. Other quantities, such as κ, $ \hat{q} $, and $ \eta_{D} $, can be calculated via their mutual relations. In this work, we follow our previous studies [54, 96] to take $ D_{s}(2\pi T)=4 $ for both charm and bottom quarks, which provides a reasonable description of the quenching and elliptic flow of D and B mesons produced in Pb+Pb collisions at the LHC. This value is 2 times larger than that from the recent Lattice QCD data [97]. Note that $ D_{s}(2\pi T) $ can depend on both the heavy quark energy and medium temperature and is not necessarily the same for charm and bottom quarks [98]. This will be left for our future improvement of this Langevin approach.

The initial momentum distribution of heavy quarks is initialized using the fixed-order-next-to-leading-log (FONLL) calculation [99–101], where the parton distribution functions are taken from the central sets of CT14NLO [102], and the nuclear shadowing effect is taken from the central sets of EPPS16 parametrization [103]. The initial spatial distribution of heavy quarks is taken from the binary collision vertices determined using the Monte-Carlo Glauber model [58]. The dynamical evolution of the QGP is simulated via the CLVisc (3+1)-dimensional viscous hydrodynamics code [104–106], where the initial condition is calculated using the TRENTo model [107]. In hydrodynamics simulation, the specific shear viscosity is set as $ \eta_\text{v}/s = 0.16 $, relaxation time is set as $ \tau_{\pi} = {3\eta_\text{v}}/{(sT)} $, and equation of state is taken from the $ \rm s95-pce-165 $ parametrization [108]. Heavy quarks start interacting with the QGP at the initial time of hydrodynamics evolution ($ \tau_0 = 0.6 $ fm/c), and the interaction stops when the local temperature of the medium drops below $ T_{c} = 160 $ MeV. These two parameters remain the same for different collision centralities. Upon exiting the QGP medium, both charm quarks and bottom quarks are converted to heavy flavor hadrons using a hybrid model that takes into account both fragmentation and coalescence processes [20]. In the end, the decay of D and B mesons into leptons is simulated using Pythia [109]. The main semi-leptonic decay modes of D meson we consider here are $ D^+ \rightarrow \overline{K}^0 l^+ \nu_l $, $ D^0 \rightarrow K^- l^+ \nu_l $, and $ D_s^+ \rightarrow \eta l^+ \nu_l $, and the main semi-leptonic decay modes of B meson are $ B^+ \rightarrow \overline{D}^0 l^+ \nu_l $, $ B^+ \rightarrow \overline{D}^* l^+ \nu_l $, $ B^0 \rightarrow \overline{D}^- l^+ \nu_l $, $ B^0 \rightarrow \overline{D}^* l^+ \nu_l $, $ B_s^0 \rightarrow \overline{D}_s^- l^+ \nu_l $, and $ B_s^0 \rightarrow \overline{D}_s^* l^+ \nu_l $. Because experimental measurements on b-decayed leptons include contributions from both $ b \rightarrow l $ and $ b \rightarrow c \rightarrow l $, leptons from $ B \rightarrow D \rightarrow l $ are categorized as b-decayed leptons in our work.

In this section, we present our numerical results on the nuclear modification factor $ R_{AA} $ and elliptic flow coefficient $ v_2 $ of heavy flavor electrons and muons produced from semi-leptonic decays of D and B mesons as functions of transverse momentum $ p_\text{T} $ in different centrality classes of Pb+Pb collisions at $ \sqrt{s_{NN}}=5.02 $ TeV.

First, Fig. 1 shows the $R_{AA}$ and $ v_2 $ for heavy flavor decayed electrons ($ c+b \to e $) as functions of $ p_\text{T} $ in various centrality classes of Pb+Pb collisions at $\sqrt{s_{NN}}=$ 5.02 TeV. For a better visual display, we present 0−10%, 20%−30%, and 30%−50% centralities in the left panels and 10%−20%, 30%−40%, and 60%−80% in the right panels. In general, our model calculation provides a good description of the existing $ R_{AA} $ data from the ALICE Collaboration. A strong centrality dependence of both $R_{AA}$ and $ v_2 $ of heavy flavor electrons can be observed here. The quenching of heavy flavor electrons is stronger (i.e., $R_{AA}$ is smaller) in central collisions than in peripheral collisions, which can be understood as the medium size effect on the heavy quark energy loss. Regarding $ v_2 $, it first increases from central to mid-central collisions but then decreases from mid-central to peripheral collisions. This results from the combined effect of medium anisotropy (eccentricity) and medium size: the former increases from central to peripheral collisions, while the latter decreases from central to peripheral collisions. Our results also show strong $ p_\text{T} $ dependence for both $R_{AA}$ and $ v_2 $ of heavy flavor electrons. As one moves from low to high $ p_\text{T} $, the value of $R_{AA}$ typically first decreases and then increases, with its minimum value appearing at intermediate $ p_\text{T} $ (around 5 GeV). In contrast, the value of $ v_2 $ first increases and then decreases, with its maximum value appearing at intermediate $ p_\text{T} $ (around 2 to 3 GeV). This may be understood as different mechanisms dominating the heavy flavor spectra at different $ p_\text{T} $ regions. At high $ p_\text{T} $, the quenching and elliptic flow of heavy flavor electrons mainly come from the quenching and elliptic flow of D and B mesons at high $ p_\text{T} $, which originate from the energy loss of charm and bottom quarks during their interaction with the (anisotropic) QGP. Meanwhile, at low and intermediate $ p_\text{T} $, the strong non-perturbative interaction and coalescence mechanism play more important roles in driving heavy quarks and hadrons towards thermal equilibrium with the medium background. Therefore, the quenching and elliptic flow of heavy flavor decayed electrons are more affected by the radial flow of the QGP medium.

Figure 1.  (color online) $R_{AA}$ (upper panel) and $ v_2 $ (lower panel) as functions of $ p_\text{T} $ for heavy flavor decayed electrons ($ c+b\to e $) in different centrality classes of Pb+Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV, compared to the ALICE data [110].

In Fig. 2, we present similar results as in Fig. 1, except for heavy flavor decayed muons ($ c+b \to \mu $). Similar to heavy flavor decayed electrons, strong centrality and $ p_\text{T} $ dependences are observed for both $ R_{AA} $ and $ v_2 $ of heavy flavor decayed muons here. For the same centrality and $ p_\text{T} $, the values of $ R_{AA} $ and $ v_2 $ are also similar between heavy flavor decayed electrons and muons, which is expected because the masses of both species of leptons are small compared to those of D and B mesons, and therefore, the decay functions to electrons and muons should be similar.

Figure 2.  (color online) $R_{AA}$ (upper panel) and $ v_2 $ (lower panel) as functions of $ p_\text{T} $ for heavy flavor decayed muons ($ c+b\to \mu $) in different centrality classes of Pb+Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV.

The $R_{AA}$ and $ v_2 $ shown in Figs. 1 and 2 are contributed by the mixture of electrons (muons) decayed from D and B mesons. To separate contributions from charm and bottom quarks and study the mass effect on heavy quark energy loss, the ATLAS Collaboration has recently measured the $ R_{AA} $ and $ v_2 $ of charm decayed and bottom decayed muons separately [87, 88]. Figure 3 shows our model calculation on the $ R_{AA} $ and $ v_2 $ of charm decayed muons as functions of $ p_\text{T} $ in Pb+Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV for various centrality classes, in comparison to the ATLAS data. Our model reasonably describes the ATLAS data in central, mid-central, and peripheral Pb+Pb collisions. Similar to Fig. 2, one can observe strong centrality and $ p_\text{T} $ dependences of the charm decayed muon $ R_{AA} $ and $ v_2 $, as a result of the combined effect of medium size, medium geometry, medium flow, and hadronization process. For the same centrality and $ p_\text{T} $, the values of $ R_{AA} $ here are smaller than those in Fig. 2, while the values of $ v_2 $ are larger. This is because the contribution from bottom decay is removed here, which yields larger $ R_{AA} $ and smaller $ v_2 $ than the contribution from charm decay within our model.

Figure 3.  (color online) $R_{AA}$ (upper panel) and $ v_2 $ (lower panel) as functions of $ p_\text{T} $ for charm decayed muons ($ c \to \mu $) in different centrality classes of Pb+Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV, compared to the ATLAS data [87, 88].

Finally, Fig. 4 presents the $R_{AA}$ and $ v_2 $ of bottom decayed muons as functions of $ p_\text{T} $ in Pb+Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV for different centrality classes, compared to the ATLAS data. One can see that while our model result reasonably agrees with the $ v_2 $ data, it overestimates $ R_{AA} $ or underestimates the energy loss effect for bottom quarks. Comparing Fig. 3 and Fig. 4, for the same centrality and $ p_\text{T} $, we see that our model calculation shows larger $ R_{AA} $ of bottom decayed muons than charm decayed muons due to weaker energy loss of heavier quarks within our model, although this hierarchy is not obvious in the current data. Regarding heavy flavor hadrons, this hierarchy is quite obvious for $ p_\text{T} \lesssim 10 $ GeV in both the experimental data and model calculations [12, 24, 72]. This discrepancy may result from the minimal assumption of the constant $ D_{s} $ parameter that we used in the present study. It is now generally acknowledged that $ D_{s} $ should depend on the medium temperature and heavy quark momentum, and its value is not necessarily the same between charm and bottom quarks. Including these temperature, momentum, and mass dependences of the diffusion coefficient should improve our model performance in the future.

Figure 4.  (color online) $ R_{AA} $ (upper panel) and $ v_2 $ (lower panel) as functions of $ p_\text{T} $ for bottom decayed muons ($ b \to \mu $) in different centrality classes of Pb+Pb collisions at $ \sqrt{s_{NN}}=5.02 $ TeV, compared to the ATLAS data [87, 88].

In this work, we studied charm and bottom quark energy loss and flow via heavy flavor leptons produced from semi-leptonic decays of D and B mesons. The CLVisc (3+1)-dimensional viscous hydrodynamics model was used to simulate the dynamical evolution of the QGP produced in Pb+Pb collisions at the LHC. The evolution of heavy quarks in the QGP was simulated using our improved Langevin model, including both collisional and radiative contributions to parton energy loss. The hadronization of heavy quarks after their escape from the QGP was simulated with our hybrid coalescence-fragmentation model. The semi-leptonic decay of D and B mesons was simulated using PYTHIA. Using our heavy quark model, we calculated the separate contributions from D and B mesons to the suppression and elliptic flow of heavy flavor decayed leptons. Both our model calculation and experimental data show strong centrality and $ p_\text{T} $ dependences of the $ R_{AA} $ and $ v_2 $ of heavy flavor leptons due to the combined effect of system size, medium anisotropy, medium flow, and the hadronization process. By using a constant $ D_{s}(2\pi T)=4 $ for both charm and bottom quarks, we found that bottom decayed leptons show larger $ R_{AA} $ and smaller $ v_2 $ than charm decayed leptons due to weaker energy loss of heavier quarks inside the QGP, though this hierarchy is not apparent in the current ATLAS data. While our model provides a reasonable description of the ALICE data on the mixed heavy flavor decayed electrons and ATLAS data on the charm decayed muons, it overestimates the $ R_{AA} $ of the bottom decayed muons observed at ATLAS. Future more detailed studies on heavy flavor hadrons and their decay products, especially with a more sophisticated temperature, momentum and mass dependences of heavy quark diffusion coefficient, can help to better constrain charm versus bottom quark energy loss, the collisional versus radiative contributions to parton energy loss, and transport properties of the QGP in relativistic heavy-ion collisions.