Study of bottom quark dynamics via non-prompt $D^{0}$ and $J/\psi$ in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV

The quark-gluon plasma (QGP), which consists of deconfined quarks and gluons as predicted by Quantum Chromodynamics (QCD) on the lattice, has been created in relativistic heavy-ion collisions performed at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) Gyulassy:2004zy ; Muller:2012zq . Extensive studies have shown that the hot and dense QGP produced in these energetic nuclear collisions have two remarkable properties: small shear viscosity to entropy density ratio (or specific viscosity) Romatschke:2017ejr ; Rischke:1995ir ; Heinz:2013th ; Gale:2013da ; Huovinen:2013wma , and high opacity to the propagation of high-energy jet partons Wang:1991xy ; Gyulassy:2003mc ; Majumder:2010qh ; Qin:2015srf ; Blaizot:2015lma ; Cao:2020wlm ; Cao:2022odi .

Heavy quarks (charm and bottom quarks) are important hard probes of the QGP Dong:2019byy ; Andronic:2015wma ; He:2022ywp . They are mostly produced from the hard scatterings in the early stage of heavy-ion collisions, and then probe the entire history of the expanding QGP. During their propagation through the QGP, heavy quarks interact with the medium and lose energy via collisional and radiative processes Cao:2013ita ; Cao:2016gvr ; Xing:2019xae ; Liu:2021dpm . The energy loss of heavy quarks in the QGP can be quantified by the nuclear modification factor $R_{\mathrm{AA}}$, defined as the ratio of the particle yield in a given centrality class in nucleus-nucleus (AA) collisions ($dN_{\mathrm{AA}}/dp_{\mathrm{T}}$), scaled by the average number of binary nucleon-nucleon collisions ($N_{\rm coll}$), to the particle yield in proton-proton ($pp$) collisions ($dN_{pp}/dp_{\mathrm{T}}$). Experimental measurements have shown that the $R_{\mathrm{AA}}$ of heavy flavor particles produced in Au+Au collisions at RHIC and Pb+Pb collisions at the LHC is significantly below unity at high $p_{\mathrm{T}}$STAR:2014wif ; CMS:2017qjw ; STAR:2018zdy , indicating substantial energy loss of heavy quarks inside the QGP medium due to their interactions with the medium constituents Gossiaux:2006yu ; Qin:2009gw ; Uphoff:2011ad ; Young:2011ug ; Alberico:2011zy ; Fochler:2013epa ; Nahrgang:2013saa ; Djordjevic:2013xoa ; Cao:2015hia ; Das:2015ana ; Song:2015ykw ; Kang:2016ofv ; Prado:2016szr ; Xu:2017obm ; Liu:2017qah ; Rapp:2018qla ; Cao:2018ews ; Li:2018izm ; Ke:2018tsh ; Li:2019wri ; Katz:2019fkc ; Li:2020kax ; Chen:2021uar ; Yang:2023rgb ; Liu:2023rfi .

Another important observable to probe the geometrical and dynamical properties of the QGP is the azimuthal anisotropy of the momentum space distribution of final state (soft and hard) hadrons Ollitrault:1992bk ; STAR:2000ekf ; STAR:2001ksn ; ALICE:2010suc ; ALICE:2011ab , which can be quantified by the Fourier coefficients of the particle distribution $dN/d\phi$ in the azimuthal plane. These anisotropy coefficients can be calculated as $v_{n}=\langle\cos[n(\phi-\Phi_{n})]\rangle$, with $\Phi_{n}$ being the $n$-th order event plane angle. For example, $v_{2}$ is called elliptic flow, which mainly originates from the elliptic shape of the produced QGP matter. Quantum fluctuations of the initial state (or colliding nuclei) can also contribute to the anisotropic flow, especially for odd flow harmonics and in ultra-central collisions Alver:2010gr ; Qin:2010pf . As for heavy quarks, due to their interactions with the anisotropic QGP, the final momentum distributions of heavy quarks and their daughter hadrons are also anisotropic. At high $p_{\mathrm{T}}$, the heavy flavor hadron $v_{2}$ is sensitive to the energy loss difference of heavy quarks along different paths through the QGP. At low and intermediate $p_{\mathrm{T}}$, their $v_{2}$ is sensitive to the collective flow of the QGP medium, since low $p_{\mathrm{T}}$ heavy quarks can pick up the QGP flow either by direct interactions with the medium through diffusive process or by coalescence with thermal light partons inside the medium during hadron formation Moore:2004tg ; He:2011qa ; Xing:2021xwc ; Plumari:2017ntm ; He:2019vgs ; Cho:2019lxb ; Cao:2019iqs ; Zhao:2023nrz . Experimental data have shown that $D$ mesons have positive $v_{2}$ at RHIC and the LHC STAR:2017kkh ; CMS:2017vhp ; ALICE:2017pbx ; ALICE:2020iug , indicating charm quarks can build up significant collective flow through scatterings with an anisotropic QGP medium as well as the hadronization process. The magnitude of the $D$ meson $v_{2}$ is similar to the light flavor hadron $v_{2}$ at intermediate $p_{\mathrm{T}}$, indicating the relaxation time of low $p_{\mathrm{T}}$ charm quarks might be comparable to or even shorter than the lifetime of the QGP.

In this work, we focus on bottom quark evolution in relativistic heavy-ion collisions. Due to their even larger mass than charm quarks, bottom quarks provide a unique probe of the QGP properties. First, they suffer much smaller cold nuclear matter effect than charm quarks Kusina:2017gkz , thus offer clean observables to study the hot nuclear matter effect. In addition, bottom quarks provide a better tool to study the “dead cone” effect which strongly depends on the mass-to-energy ratio ($m_{Q}/E$) of heavy quarks Dokshitzer:2001zm ; Armesto:2003jh ; Zhang:2003wk ; Djordjevic:2003zk ; Zhang:2018nie . Therefore, at low to intermediate $p_{\mathrm{T}}$, bottom quarks are expected to experience less diffusion than charm quarks, and thus are harder to reach thermalization with QGP medium Moore:2004tg ; Xing:2021xwc ; Liu:2023rfi . Recently, ALICE, ATLAS and CMS Collaborations have measured the $R_{\mathrm{AA}}$ and $v_{2}$ of non-prompt $D^{0}$ and $J/\psi$ decayed from bottom hadrons CMS:2018bwt ; ALICE:2022tji ; ALICE:2023gjj ; ATLAS:2018xms ; ATLAS:2018hqe ; ATLAS:2020yxw ; ATLAS:2021xtw ; ALICE:2019nuy ; CMS:2023mtk ; ALICE:2023hou . Considering the large uncertainties of the current data directly on $B$ mesons, these measurements on non-prompt $D^{0}$ and $J/\psi$ provide an important supplementary opportunity for studying bottom quark interaction with the QGP in relativistic heavy-ion collisions. We will report our study on the $R_{\mathrm{AA}}$ and $v_{2}$ of bottom decayed $D^{0}$ and $J/\psi$, and compare them to results of prompt $D$ mesons directly produced from charm quark hadronization, in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV from low to high $p_{\mathrm{T}}$ regimes. The rest of this paper will be organized as follows. In Sec. II, we will present our theoretical framework to study heavy quark evolution in relativistic heavy-ion collisions. The numerical results for prompt $D$ mesons, $B$ mesons, and non-prompt $D^{0}$ and $J/\psi$ will be presented and compared to available data at the LHC in Sec. III. Section IV will contain our summary.

In this work, we use our linear Boltzmann transport model that combines perturbative and non-perturbative interactions (LBT-PNP) Xing:2021xwc to simulate heavy quark scatterings through a color-deconfined medium. In the LBT model He:2015pra ; Cao:2017hhk ; Cao:2016gvr , one solves the following Boltzmann equation for the evolution of the phase space distribution of heavy quarks (denoted by “$a$”) inside the QGP using the Monte-Carlo method:

$\displaystyle P_{a}\cdot\partial f_{a}(x,p)=E_{a}(\mathcal{C}_{\mathrm{el}}+\mathcal{C}_{\mathrm{inel}}).$

(1)

The right hand side includes the contributions from both elastic and inelastic scatterings between heavy quarks and constituent partons of the medium, as denoted by the collision integrals $\mathcal{C}_{\mathrm{el}}$ and $\mathcal{C}_{\mathrm{inel}}$, respectively.

To simulate elastic scatterings between heavy quarks and medium partons, one calculates the scattering rates $\Gamma_{ab\rightarrow cd}$ for a binary collision process $a+b\to c+d$ using the following formula:

	$\displaystyle\Gamma_{ab\rightarrow cd}$	$\displaystyle(\vec{p}_{a},T)=\frac{\gamma_{b}}{2E_{a}}\int\frac{d^{3}p_{b}}{(2\pi)^{3}2E_{b}}\frac{d^{3}p_{c}}{(2\pi)^{3}2E_{c}}\frac{d^{3}p_{d}}{(2\pi)^{3}2E_{d}}$
	$\displaystyle\times\,$	$\displaystyle f_{b}(\vec{p}_{b},T)[1\pm f_{c}(\vec{p}_{c},T)][1\pm f_{d}(\vec{p}_{d},T)]$
	$\displaystyle\times\,$	$\displaystyle\theta(s-(m_{a}+\mu_{d})^{2})$
	$\displaystyle\times\,$	$\displaystyle(2\pi)^{4}\delta^{(4)}(p_{a}+p_{b}-p_{c}-p_{d})|\mathcal{M}_{ab\rightarrow cd}|^{2}.$		(2)

In the above equation, $\gamma_{b}$ is the degeneracy factor of parton $b$, $f_{b}$ and $f_{d}$ are thermal distributions of the medium partons, $1\pm f$ is the Bose enhancement or Fermi suppression factor for the final states (neglected for heavy quark $c$ due to their dilute distribution in this work), and the $\theta$-function accounts for the thermal mass effect on medium partons, where $\mu_{d}$ represents the Debye mass. In this calculation, the masses of charm and bottom quarks are taken as $M_{c}=1.27$ GeV and $M_{b}=4.19$ GeV, and the medium partons are taken to be massless. The key information about the microscopic elastic scattering process is contained in the matrix element $\mathcal{M}_{ab\rightarrow cd}$.

In the LBT-PNP model Xing:2021xwc , the matrix element $\mathcal{M}_{ab\rightarrow cd}$ for elastic scatterings includes both perturbative and non-perturbative interactions between heavy quarks and medium constituents. More specifically, we use the following parameterized Cornell-type potential for the interaction between a heavy quark and a medium parton,

$\displaystyle V(r)=V_{\mathrm{Y}}(r)+V_{\mathrm{S}}(r)=-\frac{4}{3}\alpha_{\mathrm{s}}\frac{e^{-m_{d}r}}{r}-\frac{\sigma e^{-m_{s}r}}{m_{s}}.$

(3)

One can see that the above potential includes both short-range Yukawa interaction and long-range color confining interaction, the later of which is also called the string term. In these two terms, $\alpha_{\mathrm{s}}$ and $\sigma$ are the coupling strengths for the Yukawa and string interactions respectively, $m_{d}$ and $m_{s}$ are their corresponding screening masses, which are taken to be temperature dependent as $m_{d}=a+bT$ and $m_{s}=\sqrt{a_{s}+b_{s}T}$. In this work, the values of the parameters $\alpha_{\mathrm{s}},\sigma,a,b,a_{s},b_{s}$ are taken from Ref. Xing:2021xwc , which provides a reasonable description of the $D$ meson observables measured at RHIC and the LHC.

To calculate the matrix element $\mathcal{M}_{ab\rightarrow cd}$, one takes the Fourier transformation and obtains the above Cornell-type potential in the momentum space as

$\displaystyle V(\vec{q})=-\frac{4\pi\alpha_{\mathrm{s}}C_{F}}{m_{d}^{2}+|\vec{q}|^{2}}-\frac{8\pi\sigma}{(m_{s}^{2}+|\vec{q}|^{2})^{2}},$

(4)

with $\vec{q}$ being the momentum exchange between heavy quarks and medium constituents. To calculate the matrix element for two-body processes: $Qq\to Qq$ and $Qg\to Qg$, in which gluons with momentum $\vec{q}$ are exchanged, we treat the above in-medium Cornell-type potential as the effective gluon propagator (field). Assuming a scalar interaction vertex for the string term, the scattering amplitude can be written as:

	$\displaystyle i\mathcal{M}=\,$	$\displaystyle\mathcal{M_{\mathrm{Y}}}+\mathcal{M_{\mathrm{S}}}$
	$\displaystyle=\,$	$\displaystyle\overline{u}(p^{\prime})\gamma^{\mu}u(p)V_{\mathrm{Y}}(\vec{q})\overline{u}(k^{\prime})\gamma^{\nu}u(k)$
		$\displaystyle+\overline{u}(p^{\prime})u(p)V_{\mathrm{S}}(\vec{q})\overline{u}(k^{\prime})u(k).$		(5)

Here, $\mathcal{M_{\mathrm{Y}}}$ and $\mathcal{M_{\mathrm{S}}}$ represent the matrix elements for the Yukawa and string terms respectively. Note that we still use a vector interaction vertex for the Yukawa term, which can reproduce the leading-order perturbative QCD result (Combridge:1978kx, ). Since the in-medium potential represents the effective gluon propagator, the string term is only included for the $t$-channel scattering, i.e., by setting $|\vec{q}|^{2}=-t$ in the above potential. Note that the color information of interaction vertices has been included in the interaction potential. For the $Qq\rightarrow Qq$ scattering process, the final amplitude squared is given by:

	$\displaystyle|\mathcal{M}_{Qq}|^{2}$	$\displaystyle=\frac{64\pi^{2}\alpha_{\mathrm{s}}^{2}}{9}\frac{(s-m_{Q}^{2})^{2}+(m_{Q}^{2}-u)^{2}+2m_{Q}^{2}t}{(t-m_{d}^{2})^{2}}$
		$\displaystyle+\frac{(8\pi\sigma)^{2}}{N_{c}^{2}-1}\frac{t^{2}-4m_{Q}^{2}t}{(t-m_{s}^{2})^{4}};$		(6)

and for the $Qg\rightarrow Qg$ process, we have

	$\displaystyle|\mathcal{M}$	${}_{Qg}|^{2}=$
		$\displaystyle\frac{64\pi^{2}\alpha_{\mathrm{s}}^{2}}{9}\frac{(s-m_{Q}^{2})(m_{Q}^{2}-u)+2m_{Q}^{2}(s+m_{Q}^{2})}{(s-m_{Q}^{2})^{2}}$
	$\displaystyle+\,$	$\displaystyle\frac{64\pi^{2}\alpha_{\mathrm{s}}^{2}}{9}\frac{(s-m_{Q}^{2})(m_{Q}^{2}-u)+2m_{Q}^{2}(u+m_{Q}^{2})}{(u-m_{Q}^{2})^{2}}$
	$\displaystyle+\,$	$\displaystyle 8\pi^{2}\alpha_{\mathrm{s}}^{2}\frac{5m_{Q}^{4}+3m_{Q}^{2}t-10m_{Q}^{2}u+4t^{2}+5tu+5u^{2}}{(t-m_{d}^{2})^{2}}$
	$\displaystyle+\,$	$\displaystyle 8\pi^{2}\alpha_{\mathrm{s}}^{2}\frac{(m_{Q}^{2}-s)(m_{Q}^{2}-u)}{(t-m_{d}^{2})^{2}}$
	$\displaystyle+\,$	$\displaystyle 16\pi^{2}\alpha_{\mathrm{s}}^{2}\frac{3m_{Q}^{4}-3m_{Q}^{2}s-m_{Q}^{2}u+s^{2}}{(s-m_{Q}^{2})(t-m_{d}^{2})}$
	$\displaystyle+\,$	$\displaystyle\frac{16\pi^{2}\alpha_{\mathrm{s}}^{2}}{9}\frac{m_{Q}^{2}(4m_{Q}^{2}-t)}{(s-m_{Q}^{2})(m_{Q}^{2}-u)}$
	$\displaystyle+\,$	$\displaystyle 16\pi^{2}\alpha_{\mathrm{s}}^{2}\frac{3m_{Q}^{4}-m_{Q}^{2}s-3m_{Q}^{2}u+u^{2}}{(t-m_{d}^{2})(u-m_{Q}^{2})}$
	$\displaystyle+\,$	$\displaystyle\frac{C_{A}}{C_{F}}\frac{(8\pi\sigma)^{2}}{N_{c}^{2}-1}\frac{t^{2}-4m_{Q}^{2}t}{(t-m_{s}^{2})^{4}},$		(7)

with $s$, $t$, $u$ being the Mandelstam variables. To obtain the Yukawa parts $|\mathcal{M_{\mathrm{Y}}}|^{2}$ in the above two expressions, we calculate $|\mathcal{M}_{s}|^{2}$, $|\mathcal{M}_{u}|^{2}$, $|\mathcal{M}_{t}|^{2}$, and their interference terms $|\mathcal{M}_{s}\mathcal{M}_{t}^{*}|$, $|\mathcal{M}_{s}\mathcal{M}_{u}^{*}|$ and $|\mathcal{M}_{u}\mathcal{M}_{t}^{*}|$ separately, and then replace $t$ by $t-m_{d}^{2}$ in the denominators. Note that the form of the $|\mathcal{M_{\mathrm{Y}}}|^{2}$ part in Eq. (II) appears a little different from the form given by Ref. (Combridge:1978kx, ). This is because the amplitude squared for the $Qg\rightarrow Qg$ scattering process is simplified in Ref. (Combridge:1978kx, ) by extracting a term $-16\pi^{2}\alpha_{\mathrm{s}}^{2}(2m_{Q}^{2}+t)/t$ from both $|\mathcal{M}_{s}\mathcal{M}_{t}^{*}|$ and $|\mathcal{M}_{u}\mathcal{M}_{t}^{*}|$, and then combine them with $|\mathcal{M}_{t}|^{2}$ to obtain $32\pi^{2}\alpha_{\mathrm{s}}^{2}(s-m_{Q}^{2})(m_{Q}^{2}-u)/t^{2}$. If one introduces the Debye screening mass, i.e, replace $t$ by $t-m_{d}^{2}$ in the denominators, based on the form in Ref. (Combridge:1978kx, ), as shown in Refs. (Svetitsky:1987gq, ; GolamMustafa:1997id, ; Liu:2016ysz, ), the matrix element would appear different from our Yukawa part in Eq. (II). However, they are the same when $m_{d}$ is set as zero.

By summing over all possible scattering channels, one can obtain the total elastic scattering rate for a heavy quark propagating through the QGP according to Eq. (II). For a given time step $\Delta\tilde{t}$, the probability for a heavy quark to experience elastic scatterings with the QGP constituents can be calculated as: $P_{\mathrm{el}}^{a}=1-e^{-\Gamma_{\mathrm{el}}^{a}\Delta\tilde{t}}$.

For inelastic scatterings between heavy quarks and the QGP, the scattering rate at a given time $\tilde{t}$ can be obtained as follows:

$\displaystyle\Gamma_{\mathrm{inel}}^{a}(E_{a},T,\tilde{t})=\int dxdl_{\bot}^{2}\frac{dN_{g}^{a}}{dxdl_{\bot}^{2}d\tilde{t}}.$

(8)

In this study, we take the higher-twist energy loss formalism Guo:2000nz ; Majumder:2009ge ; Zhang:2003wk ; Zhang:2018nie for the gluon emission spectrum off a heavy quark inside a dense nuclear medium,

$\displaystyle\frac{dN_{g}^{a}}{dxdl_{\bot}^{2}d\tilde{t}}=\frac{2C_{A}\alpha_{\mathrm{s}}P_{a}(x)l_{\bot}^{4}\hat{q}_{a}}{\pi(l_{\bot}^{2}+x^{2}m_{a}^{2})^{4}}\sin^{2}\left(\frac{\tilde{t}-\tilde{t}_{i}}{2\tau_{f}}\right).$

(9)

In the above two equations, $E_{a}$ and $m_{a}$ are the energy and mass of heavy quarks, $x$ and $l_{\perp}$ represent the energy fraction and the transverse momentum of the radiated gluon with respect to the parent heavy quark, and $P_{a}(x)$ is the vacuum splitting function. The jet transport coefficient $\hat{q}_{a}$ represents the average transverse momentum squared exchanged between heavy quarks and the medium constituents per unit time during the elastic scattering process. Inside the sine function, $\tilde{t}-\tilde{t}_{i}$ is the time accumulated from the previous emission time ($\tilde{t}_{i}$), and $\tau_{f}=2E_{a}x(1-x)/(l_{\bot}^{2}+x^{2}m_{a}^{2})$ denotes the formation time of gluon emission. Based on the above formula, the probability for a heavy quark to experience inelastic scatterings with the QGP constituents during a time step $\Delta\tilde{t}$ can be calculated as: $P_{\mathrm{inel}}^{a}=1-e^{-\Gamma_{\mathrm{inel}}^{a}\Delta\tilde{t}}$.

To include both elastic and inelastic scatterings between heavy quarks and the QGP medium, the total scattering rate can be obtained as: $\Gamma_{\mathrm{tot}}=\Gamma_{\mathrm{el}}+\Gamma_{\mathrm{inel}}$. In terms of the total scattering probability, one may write out the following expression:

$\displaystyle P_{\mathrm{tot}}^{a}=1-e^{-\Gamma_{\mathrm{tot}}\Delta t}=P_{\mathrm{el}}^{a}+P_{\mathrm{inel}}^{a}-P_{\mathrm{el}}^{a}P_{\mathrm{inel}}^{a},$

(10)

where $P_{\mathrm{el}}^{a}(1-P_{\mathrm{inel}}^{a})$ can be understood as the pure elastic scattering probability, and $P_{\mathrm{inel}}^{a}$ is the inelastic scattering probability. Following our earlier work Cao:2017hhk ; Xing:2019xae ; Xing:2021xwc , we use different values of the coupling strength $\alpha_{\mathrm{s}}$ for different interaction vertices in our calculation of elastic and inelastic scatterings. For a vertex connecting to the propagating heavy quarks, we take $\alpha_{\mathrm{s}}=4\pi/[9\mathrm{ln}(2ET/\Lambda^{2})]$ with $\Lambda=0.2$ GeV; for a vertex connecting to the medium partons, we use the same value as in the interaction potential $V(r)$.

To evolve heavy quarks through a realistic medium, we utilize the (3+1)-dimensional viscous hydrodynamic model CLVisc Pang:2012he ; Pang:2018zzo ; Wu:2018cpc ; Wu:2021fjf to simulate the dynamical evolution of the QGP fireball produced in relativistic heavy-ion collisions at the LHC. In the present study, we apply smooth hydrodynamic profiles for investigating heavy quarks, whose specific shear viscosity is taken to be $\eta/s=0.16$ for the QGP produced in Pb+Pb collision at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV. The initial energy density distribution of the QGP is obtained from the Glauber model, which is also used to obtain the initial spatial distribution of heavy quarks. The momentum distribution of heavy quarks is initialized using the fixed-order-next-to-leading-log (FONLL) calculation Cacciari:2001td ; Cacciari:2012ny ; Cacciari:2015fta , where the parton distribution functions are taken from CT14NLO Dulat:2015mca , and the nuclear shadowing effect is taken from EPPS16 Eskola:2016oht at the next-to-leading-order. After their production, heavy quarks are assumed to stream freely before starting interaction with the QGP at the initial proper time of hydrodynamic evolution ($\tau_{0}=0.6$ fm/$c$). Interaction between heavy quarks and the QGP ceases when heavy quarks reach the QGP boundary, the hypersurface with $T_{\rm c}=165$ MeV in this work, on which they are converted into heavy flavor hadrons via a hybrid fragmentation-coalescence model Cao:2019iqs . The decay of $B$ mesons into non-prompt $D^{0}$ and $J/\psi$ is simulated by Pythia Sjostrand:2006za .

In this section, we present our numerical results on the nuclear modification factor ($R_{\mathrm{AA}}$) and elliptic flow coefficient ($v_{2}$) of non-prompt $D^{0}$ and $J/\psi$ in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV at the LHC.

Before showing results for non-prompt $D^{0}$ and $J/\psi$, we first show in Fig. 1 the $R_{\mathrm{AA}}$ and $v_{2}$ of prompt $D^{0}$ mesons as functions of transverse momentum $p_{\mathrm{T}}$ in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV for three different centrality classes: central (0-10%), mid-central (30-50%) and peripheral (60-80%) collisions. The ALICE and CMS data on 0-10% and 30-50% centralities are shown for comparison. Our model can provide a reasonable description of the LHC data on the $D^{0}$ meson $R_{\mathrm{AA}}$ and $v_{2}$. One can see that both $R_{\mathrm{AA}}$ and $v_{2}$ show strong centrality dependence. From central to mid-central to peripheral collisions, the quenching of the prompt $D$ mesons decreases ($R_{\mathrm{AA}}$ becomes larger at large $p_{\mathrm{T}}$) due to the decreasing system size, whereas the prompt $D$ meson $v_{2}$ first increases due to the increasing eccentricity of QGP medium, and then decreases due to the decreasing energy loss through a smaller system. Meanwhile, both $R_{\mathrm{AA}}$ and $v_{2}$ show strong $p_{\mathrm{T}}$ dependence. The prompt $D$ meson $R_{\mathrm{AA}}$ first exhibits a bump structure at low $p_{\mathrm{T}}$ and then increases with $p_{\mathrm{T}}$ at high $p_{\mathrm{T}}$, whereas the prompt $D$ meson $v_{2}$ first increases and then decreases. These structures originate from the combined effect of heavy quark spectrum, energy loss and hadronization, together with the QGP flow. At large $p_{\mathrm{T}}$, the quenching of prompt $D$ mesons mainly comes from charm quark energy loss inside the QGP, and the $v_{2}$ of prompt $D$ mesons mainly results from the anisotropic energy loss along different directions through a geometrically asymmetric medium. The decreasing fractional energy loss and the flatter $p_{\mathrm{T}}$ spectrum of charm quarks at high $p_{\mathrm{T}}$ lead to the increasing $R_{\mathrm{AA}}$ and decreasing $v_{2}$ of prompt $D$ mesons with $p_{\mathrm{T}}$. At low to intermediate $p_{\mathrm{T}}$, the strong non-perturbative interactions can quickly drive heavy quarks towards thermal equilibrium with the QGP, and therefore the motion of heavy quarks are strongly affected by the QGP flow. In addition, the coalescence between heavy quarks and the medium partons dominates the heavy flavor hadron formation at low to intermediate $p_{\mathrm{T}}$, which further enhances the QGP flow effect on the prompt $D$ meson spectrum. These lead to the flow bump of the prompt $D$ meson $R_{\mathrm{AA}}$, and the increasing $D$ meson $v_{2}$ at low $p_{\mathrm{T}}$, similar to the behaviors of light flavor hadron $R_{\mathrm{AA}}$ and $v_{2}$ observed at low $p_{\mathrm{T}}$.

Figure 2 shows our prediction for the $R_{\mathrm{AA}}$ and $v_{2}$ of $B$ mesons as functions of $p_{\mathrm{T}}$ in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV for 0-10%, 30-50% and 60-80% centrality classes. Here, results for $B$ mesons include contributions from $B^{+}$ and $B^{-}$ in our simulation. Similar to previous results for prompt $D$ mesons, the $B$ meson $R_{\mathrm{AA}}$ and $v_{2}$ also show strong centrality dependence. Moving from central to mid-central to peripheral collisions, the quenching of $B$ mesons becomes weaker due to the decreasing system size, while the $B$ meson $v_{2}$ first increases and then decreases due to the combined effect of medium eccentricity and size. For the same centrality class, the $B$ meson $R_{\mathrm{AA}}$ is larger than the prompt $D$ meson $R_{\mathrm{AA}}$ at low $p_{\mathrm{T}}$, while the former is comparable to the latter or even slightly smaller than the latter at high $p_{\mathrm{T}}$. This is due to the non-trivial mass hierarchy of quark energy loss within our LBT-PNP model, as discussed in an earlier study Dang:2023tmb : $b$-quarks lose less energy than $c$-quarks at low $p_{\mathrm{T}}$ due to the “dead cone effect” in the gluon emission process, while the former may lose more energy than the latter at high $p_{\mathrm{T}}$ due to the string interaction implemented in our model. Within the same centrality class, the $v_{2}$ of $B$ mesons is smaller than that of prompt $D$ mesons at low $p_{\mathrm{T}}$. At high $p_{\mathrm{T}}$, they are both small.

To further study the quenching and flow of $B$ mesons, Recently, ALICE, ATLAS and CMS Collaborations have measured the $R_{\mathrm{AA}}$ and $v_{2}$ of non-prompt $D^{0}$ and $J/\psi$, decayed from bottom hadrons CMS:2018bwt ; ALICE:2022tji ; ALICE:2023gjj ; ATLAS:2018xms ; ATLAS:2018hqe ; CMS:2023mtk ; ALICE:2023hou . In Fig. 3, we show the $R_{\mathrm{AA}}$ and $v_{2}$ for non-prompt $D^{0}$ as functions of $p_{\mathrm{T}}$ in central (0-10%), mid-central (30-50%) and peripheral (60-80%) Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, compared to the ALICE and CMS data for 0-10% and 30-50% centralities. One can see that our model provides a reasonable description of the non-prompt $D^{0}$ meson $R_{\mathrm{AA}}$ and $v_{2}$ at high $p_{\mathrm{T}}$. At low $p_{\mathrm{T}}$, our model underestimates both $R_{\mathrm{AA}}$ and $v_{2}$, which may result from deficiencies in evaluating non-perturbative processes at low $p_{\mathrm{T}}$, including the initial heavy quark spectra, string interactions between heavy quarks and the QGP, and hadronization of heavy quarks. Similar to prompt $D$ and $B$ mesons, the $R_{\mathrm{AA}}$ and $v_{2}$ of non-prompt $D^{0}$ mesons here also show strong dependence on centrality. In the same centrality class, the $R_{\mathrm{AA}}$ and $v_{2}$ of non-prompt $D^{0}$ show similar behaviors as those of $B$ mesons, except for some $p_{\mathrm{T}}$ shift during the decay from $B$ mesons to $D$ mesons.

In Fig. 4, we show the $R_{\mathrm{AA}}$ of non-prompt $J/\psi$ as a function of $p_{\mathrm{T}}$ in 0-80% Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV and the $v_{2}$ of non-prompt $J/\psi$ in 10-60% collisions. Our model provides a good description of the corresponding ATLAS data on $R_{\mathrm{AA}}$ and the CMS data on $v_{2}$ of non-prompt $J/\psi$. Both the $R_{\mathrm{AA}}$ and $v_{2}$ results here are obtained from averaging over a large interval of centrality. They do not show a strong dependence on $p_{\mathrm{T}}$.

In order to investigate the centrality dependence of quenching and flow of non-prompt $J/\psi$, we present in Fig. 5 the $R_{\mathrm{AA}}$ and $v_{2}$ of non-prompt $J/\psi$ as functions of $p_{\mathrm{T}}$ in central (0-10%), mid-central (30-50%) and peripheral (60-80%) Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, compared to the available data from the ALICE, ATLAS and CMS collaborations. One can see that our model provides a reasonable description of the corresponding data within comparable centrality bins. Similar to previous results of prompt $D$ mesons, $B$ mesons and non-prompt $D^{0}$, the $R_{\mathrm{AA}}$ and $v_{2}$ of non-prompt $J/\psi$ show strong centrality dependence.

In Fig. 6, we directly compare the quenching and flow between different species of heavy flavor particles presented above, left panel for the $R_{\mathrm{AA}}$ of prompt $D$, $B$, non-prompt $D^{0}$ and non-prompt $J/\psi$ in central (0-10%) Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, and right panel for their $v_{2}$ in mid-central (30-50%) collisions, in comparison to the available ALICE, ATLAS and CMS data. As discussed earlier, the $R_{\mathrm{AA}}$ of direct $D$ mesons is smaller than that of $B$ mesons at low $p_{\mathrm{T}}$, while the inverse order is seen at high $p_{\mathrm{T}}$. This is because of the opposite mass dependences of quark energy loss at low $p_{\mathrm{T}}$ and high $p_{\mathrm{T}}$ within our LBT-PNP model. In mid-central collisions, the $v_{2}$ of direct $D$ mesons appear larger than that of $B$ mesons. The $R_{\mathrm{AA}}$ of non-prompt $D^{0}$ and non-prompt $J/\psi$ are smaller than that of $B$ mesons at low $p_{\mathrm{T}}$, but larger at high $p_{\mathrm{T}}$, due to the $p_{\mathrm{T}}$ shift during the decay of $B$ mesons. The $v_{2}$ of non-prompt $D^{0}$ and non-prompt $J/\psi$ are smaller than that of $B$ mesons at low $p_{\mathrm{T}}$, but comparable at high $p_{\mathrm{T}}$. No apparent difference in $R_{\mathrm{AA}}$ and $v_{2}$ is observed between non-prompt $D^{0}$ and non-prompt $J/\psi$ due to their similar decay functions from $B$ mesons.

In Fig. 7, we compare the $p_{\mathrm{T}}$-integrated $R_{\mathrm{AA}}$ between prompt $D^{0}$, $B$, non-prompt $D^{0}$ and non-prompt $J/\psi$ in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, left panel for $1.5<p_{\mathrm{T}}<10$ GeV and right panel for $9<p_{\mathrm{T}}<40$ GeV. At higher $p_{\mathrm{T}}$ region (right panel), a clear participant number ($N_{\mathrm{part}}$), or centrality, dependence of $R_{\mathrm{AA}}$ can be seen for these heavy flavor particles: stronger energy loss of heavy quarks in more central collisions (larger $N_{\mathrm{part}}$) leads to smaller $R_{\mathrm{AA}}$. This trend is not apparent at low $p_{\mathrm{T}}$ (left panel), where in addition to parton energy loss, the hadronization process and the QGP flow also significantly affect the final state hadron spectra. Within the $p_{\mathrm{T}}$ range we explore here, the $R_{\mathrm{AA}}$ of direct $D$ mesons is smaller than that of $B$ mesons. No obvious difference is observed between $B$ mesons, non-prompt $D^{0}$ and non-prompt $J/\psi$ in these $p_{\mathrm{T}}$-integrated $R_{\mathrm{AA}}$. Since our current model underestimates the heavy flavor $R_{\mathrm{AA}}$ at low $p_{\mathrm{T}}$, as shown in previous $p_{\mathrm{T}}$-dependent $R_{\mathrm{AA}}$ results, our result on the $p_{\mathrm{T}}$-integrated $R_{\mathrm{AA}}$ of non-prompt $J/\psi$ here is also lower than the available data from the ALICE and ATLAS collaborations. The agreement becomes better as the $p_{\mathrm{T}}$ range becomes higher.

In the end, we compare in Fig. 8 the $p_{\mathrm{T}}$-integrated $v_{2}$ between prompt $D^{0}$, $B$, non-prompt $D^{0}$ and non-prompt $J/\psi$ in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, left panel for $1.5<p_{\mathrm{T}}<10$ GeV and right panel for $6.5<p_{\mathrm{T}}<50$ GeV. Moving from central to mid-central to peripheral collisions, or as $N_{\mathrm{part}}$ becomes smaller, the elliptic flow coefficients first increase and then decrease due to the combined effect of medium eccentricity and medium size. At low $p_{\mathrm{T}}$, the prompt $D^{0}$ mesons have much larger $v_{2}$ than $B$ mesons, non-prompt $D^{0}$ and non-prompt $J/\psi$. This difference becomes smaller at higher $p_{\mathrm{T}}$. Compared to the available data from the CMS collaboration, our model provides a reasonable description of the $v_{2}$ of non-prompt $J/\psi$.

Within the linear Boltzmann transport model that includes both string and Yukawa types of interactions between heavy quarks and the QGP, we study the dynamics of bottom quarks in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV at the LHC via the nuclear modification factors and elliptic flow coefficients of $B$ mesons, bottom decayed $D^{0}$ and $J/\psi$. Compared to direct $D$ mesons, $B$ mesons show larger $R_{\mathrm{AA}}$ at low $p_{\mathrm{T}}$ but slightly smaller $R_{\mathrm{AA}}$ at high $p_{\mathrm{T}}$, indicating weaker energy loss of heavier quarks at low $p_{\mathrm{T}}$ but a possible inverse order at high $p_{\mathrm{T}}$ within our LBT-PNP model. Both $R_{\mathrm{AA}}$ and $v_{2}$ of $D$ and $B$ mesons show strong $p_{\mathrm{T}}$ and centrality dependences. At high $p_{\mathrm{T}}$, the heavy meson $R_{\mathrm{AA}}$ increases with $p_{\mathrm{T}}$ due to the decreasing fractional energy loss and the flatter $p_{\mathrm{T}}$ spectra of heavy quarks at higher $p_{\mathrm{T}}$. A bump structure of $R_{\mathrm{AA}}$ can be observed at low $p_{\mathrm{T}}$, resulting from the QGP flow effect and the coalescence process of heavy quark hadronization. The heavy meson $v_{2}$ is driven by the anisotropic QGP flow at low $p_{\mathrm{T}}$, while driven by anisotropic energy loss through different directions at high $p_{\mathrm{T}}$, and therefore, first increases and then decreases as $p_{\mathrm{T}}$ becomes larger. From central to mid-central to peripheral collisions, the heavy meson $R_{\mathrm{AA}}$ becomes larger due to weaker energy loss through a smaller QGP system. On the other hand, their $v_{2}$ first increases due to larger medium eccentricity, and then decreases due to smaller medium size. Non-prompt $D^{0}$ and non-prompt $J/\psi$ show very similar $R_{\mathrm{AA}}$ and $v_{2}$ to $B$ mesons, except for a shift towards the lower $p_{\mathrm{T}}$ region. Compared to the available ALICE, ATLAS and CMS data, our model provides a reasonable description of the nuclear modification factors and elliptic flow coefficients of prompt $D$ mesons, non-prompt $D^{0}$ and non-prompt $J/\psi$, except for some deviation at low $p_{\mathrm{T}}$ due to possible inaccurate description of initial spectrum of heavy quarks, their interaction with the QGP and their hadronization process in this non-perturbative region. Therefore, studying non-prompt $D^{0}$ and $J/\psi$ provides a supplementary way for better understanding heavy quark dynamics in relativistic heavy-ion collisions.

This work is supported in part by the National Natural Science Foundation of China (NSFC) under Grant Nos. 12225503, 11890710, 11890711, 11935007, 12175122 and 2021-867. W.-J. X. is supported in part by China Postdoctoral Science Foundation under Grant No. 2023M742099. Some of the calculations were performed in the Nuclear Science Computing Center at Central China Normal University (NSC³), Wuhan, Hubei, China.