Statistical Production of $B_{c}$ Mesons in Heavy-Ion Collisions at the LHC Energy

The $B_{c}$ mesons as bound states of a bottom ($b$) quark with an anticharm ($\bar{c}$) quark, or vice versa, provide a unique avenue to the understanding of heavy quark dynamics. They are intermediate between charmonium ($c\bar{c}$) and bottomonium ($b\bar{b}$) states in terms of mass, size and binding energy [1, 2, 3, 4]. Being composed of two different heavy flavors, $B_{c}$ mesons cannot annihilate into gluons; consequently excited $B_{c}$ states lying below the open $BD$ threshold can only undergo radiative or hadronic transitions into the pseudoscalar (${}^{1}S_{0}$) ground state $B_{c}$ that decays weakly, leading to total widths less than a few hundred keV [1, 2], significantly smaller than their charmonia and bottomnia counterparts [5].

The production of $B_{c}$ mesons in hadronic (e.g. proton-proton, $pp$) collisions entails the simultaneous creation of a $c\bar{c}$ and a $b\bar{b}$ pair in a single collision, rendering it much rarer than that of other mesons containing a single $b$ quark. In $pp$ collisions at the LHC energy, the LHCb experiment has measured the production fraction of the ground state $B_{c}^{-}$ to be $\sim 0.26$% (with significant uncertainty due to the one from the branching fraction of the decay $B_{c}^{-}\rightarrow J/\psi\mu^{-}\bar{\nu}$ used in the measurement) relative to the total $b\bar{b}$ cross section [6], which is negligibly small compared to the fractions of $B$ mesons and $\Lambda_{b}$ baryons [7, 8].

However, in relativistic heavy-ion collisions where a deconfined medium (Quark-Gluon Plasma, QGP) is created, abundant heavy quarks are produced through primordial hard processes and then transported in the QGP [9, 10, 11, 12, 13], leading to a new production mechanism for $B_{c}$ mesons, namely the recombination of $b$ and $c$ quarks from different primordial nucleon-nucleon collisions that could significantly enhance the $B_{c}$ yield [14]. For example, in the most central 0-10% Pb-Pb collisions at the LHC energies, $dN_{c\bar{c}}/dy\sim 20$ [15] charm and $dN_{b\bar{b}}/dy\sim 0.9$ [8] bottom quarks per unit rapidity can be produced at mid-rapidity, respectively. In events in which there is a single $b\bar{b}$ pair along with a number of $c\bar{c}$ pairs, the $b$ quark can find any of the $c$ quarks near in phase space to recombine and produce a $B_{c}$ bound state. This recombination process is particularly strengthened given the high degree of thermalization of $c$ quarks in the QGP [13, 15, 16] as primarily evidenced by the large elliptic flow of $D$ mesons [17] measured in Pb-Pb collisions at the LHC energy.

This novel production mechanism makes $B_{c}$ an attractive part of the endeavor of using heavy quarkonia as a probe of the QGP properties [18, 19, 20, 21, 22, 23, 24]. In this context, the recombination production of $B_{c}$ mesons, since it was first proposed [14], has been studied in various transport [25, 26, 27] and instantaneous coalescence [28] models. The static in-medium properties of $B_{c}$ mesons such as binding energy and radius used as inputs of these models are sensitive to the employed heavy quark potential that is currently under hot debate [29, 30]. On the experimental side, pioneering measurement by the CMS experiment in $\sqrt{s_{\rm NN}}=5.02$ TeV Pb-Pb collisions, although restricted to relatively large transverse momenta ($p_{T}>6$ GeV), indeed gives a first hint that the $B_{c}$ production in the presence of QGP is enhanced relative to $pp$ collisions, as indicated by the nuclear modification factor well above unity in the lower $p_{T}$ bin accessed in the measurement [31]. One notes that, although the recombination of $B_{c}$ depends linearly on the $c$ quark densities, weaker than the quadratic dependence in the case of charmonia regeneration [32], the nuclear modification factor of $B_{c}$ is likely to to be significantly larger than that of $J/\psi$ mesons [33], owing to the too small production of the former in $pp$ collisions as the reference.

In this work, we treat $B_{c}$ as a member of the family of $b$-hadrons that contain a single $b$ quark and study its production in $\sqrt{s_{\rm NN}}=5.02$ TeV Pb-Pb collisions in the statistical hadronization model (SHM) following the spirit of the statistical production of charmonia and $c$-hadrons [34, 35, 36, 37]. For heavy quarkonia that feature deep binding, statistical coalescence (recombination) does not apply to their production in elementary collisions, in contrast to the case of heavy-light hadrons where SHM empirically works well [8, 38]. However, the high degree of thermalization for $c$ quarks facilitates the statistical description of not only charmonia [34, 39] but also $B_{c}$ mesons in Pb-Pb collisions at the LHC energy. For the latter, one can envisage that a single $b$ quark is surrounded by a number of highly thermalized diffusing $c$ quarks and stochastic (statistical) recombination satisfying phase space conditions occurs efficiently, which is reminiscent of the case of statistical $c$-light quark recombination to form open $c$-hadrons [36]. Given the small number of heavy quark production, it turns out important to adopt the canonical ensemble version of SHM to exactly conserve the $b$ and $c$ numbers for predicting the absolute yield of $B_{c}$ mesons. This statistical production description in particular allows us to make quantitative predictions for the production fraction of $B_{c}$ mesons relative to the total $b\bar{b}$ multiplicity and the enhancement of its relative production to $B$ mesons in the presence of QGP with respect to $pp$ collisions.

This article is organized as follows. In Sec. 2, we fix $B_{c}$’s $p_{T}$ differential cross section in $\sqrt{s}=5.02$ TeV $pp$ collisions by making use of the previously computed $B^{-}$’s differential cross section in combination with the $p_{T}$ dependent $B_{c}^{-}/B^{-}$ ratio measured by LHCb experiment. This will be used as the reference for the ensuing calculation of $B_{c}$’s nuclear modification factors and also enables us to make an estimate of its integrated cross section. In Sec. 3, we elaborate on the SHM calculation of $B_{c}$ production in $\sqrt{s_{\rm NN}}=5.02$ TeV Pb-Pb collisions, putting it in the context of the statistical production of the whole family of $b$-hadrons such that the $b$-number conservation through hadronization is taken into full account and serves as an constraint for the prediction of the production fraction of $B_{c}$ mesons. In Sec. 4, we compute the $B_{c}$’s $p_{T}$ distribution by normalizing its $p_{T}$ spectrum from resonance recombination of realistically transported $c$ and $b$ quarks to its statistical production yield, which is then supplemented by the component from fragmentation of $b$ quarks that dominates at high $p_{T}$. The total $p_{T}$ spectrum is converted into the $p_{T}$ dependent nuclear modification factors and compared to CMS data. We conclude in Sec. 5.

In this section, we construct the $p_{T}$ differential cross section $d\sigma/dp_{T}dy$ for the ground state $B_{c}^{-}$ in $\sqrt{s}=5.02$ TeV $pp$ collisions at mid-rapidity as the reference to measure its $p_{T}$ spectral modifications in Pb-Pb collisions. The production fraction of the ground state pseudoscalar $B_{c}^{-}$ mesons relative to the sum of those of $B^{-}$ and $\bar{B}^{0}$ ($f_{c}/(f_{u}+f_{d})$) has been measured by the LHCb experiment in $\sqrt{s}=7$ and 13 TeV $pp$ collisions in the kinematic regions of transverse momentum $4<p_{T}<25$ GeV and pseudo-rapidity $2.5<\eta<4.5$, showing almost the same results between these two energies [6]. The measured ratio shows no rapidity dependence and its rather weak $p_{T}$ dependence has been linearly fitted [6], which is converted into $B_{c}^{-}/B^{-}$ assuming isospin symmetry ($f_{u}=f_{d}$) and plotted in Fig. 2(b) (blue band), where the significant spread is due to the uncertainty from the branching fraction of $\mathcal{B}(B_{c}^{-}\rightarrow J/\psi\mu^{-}\bar{\nu})=1.95\%\pm 0.46\%$ [6] used in the measurement.

Assuming the $p_{T}$-dependent $B_{c}^{-}/B^{-}$ to remain unchanged, we construct $B_{c}^{-}$’s $d\sigma/dp_{T}dy$ at mid-rapidity in $\sqrt{s}=5.02$ TeV $pp$ collisions by multiplying the $B_{c}^{-}/B^{-}$ with $B^{-}$’s $d\sigma/dp_{T}dy$ that has been determined in [8]. The resulting $d\sigma/dp_{T}dy$ for $B_{c}^{-}$ and that for $B^{-}$ quoted from [8] (in comparison with CMS measurement [40]) are displayed in Fig. 1(a). As a byproduct, we’ve estimated the $B_{c}$’s integrated cross section to be $d\sigma^{B_{c}^{-}}/dy=0.104\pm 0.0245~{}\mu b$ at mid-rapidity and its relative ratio to $d\sigma^{B^{-}}/dy=13.16~{}\mu b$ for $B^{-}$ to be $B_{c}^{-}/B^{-}=0.00790\pm 0.00186$. Based on this ratio, the production fraction of $B_{c}^{-}$ was also estimated to be $f_{c}=0.00258\pm 0.00062$ [6] (only the dominant uncertainty from $\mathcal{B}(B_{c}^{-}\rightarrow J/\psi\mu^{-}\bar{\nu})$ is quoted).

The CMS experiment has measured the $d\sigma/dp_{T}dy$ for $B_{c}^{-}+B_{c}^{+}$ multiplied by the pertinent branching fractions through a two-step decay process $B_{c}\rightarrow(J/\psi\rightarrow\mu\mu)\mu\nu$ in two bins of the three-muons’ transverse momentum $6<p_{T}^{\mu\mu\mu}<11$ GeV at $1.3<|y^{\mu\mu\mu}|<2.3$ and $11<p_{T}^{\mu\mu\mu}<35$ GeV at $|y^{\mu\mu\mu}|<2.3$ [31]. To compare to the CMS data, we multiply the $B_{c}^{-}$’s $d\sigma/dp_{T}dy$ constructed above by $\mathcal{B}(B_{c}^{-}\rightarrow J/\psi\mu^{-}\bar{\nu})=1.95\%\pm 0.46\%$ [6] and $\mathcal{B}(J/\psi\rightarrow\mu\mu)=5.96\%$ (with negligible uncertainty) [5]. Taking into account the average $p_{T}$ shift of $p_{T}^{\mu\mu\mu}=0.85\cdot p_{T}^{B_{c}}$ [31], the result (the factor $2$ is to account for the sum of $B_{c}^{-}+B_{c}^{+}$) is plotted in Fig. 1(b) in comparison with the CMS data [31]. The central values of CMS data in two rather wide $p_{T}$ bins are better described by using the lower bound of $\mathcal{B}(B_{c}^{-}\rightarrow J/\psi\mu^{-}\bar{\nu})$.

In this section, we compute the absolute production yields of $B_{c}^{-}$ in different centralities of $\sqrt{s_{\rm NN}}=5.02$ TeV Pb-Pb collisions from the SHM. For large systems like heavy-ion collisions, quantum charges such as electric charge ($Q$), baryon number ($N$) and strangeness ($S$) can be treated grand-canonically in the statistical description. In contrast, the small number of $c$-hadrons ($C$) and $b$-hadrons ($B$) renders their exact conservation through the canonical ensemble description essential, such that the partition function reads [38, 41]

where $\lambda_{Q}$, $\lambda_{N}$ and $\lambda_{S}$ are fugacities associated with $Q$, $N$ and $S$, respectively; $C_{j}$ and $B_{j}$ are the charm and bottom number of the $j$-th particle, respectively; $\gamma_{s}$, $\gamma_{c}$ and $\gamma_{b}$ are fugacities that account for the deviation from chemical equilibrium for hadrons containing $N_{sj}$, $N_{cj}$ and $N_{bj}$ strange, $c$ and $b$ quarks plus antiquarks, respectively. In Eq. (3), $z_{j}$ denotes the one-particle partition function

which also specifies the chemical equilibrium multiplicity of the $j$-th hadron of mass $m_{j}$ and spin $J_{j}$ in a fireball of volume $V$ under the Boltzmann approximation at hadronization temperature $T_{H}$, with $K_{2}$ being the modified Bessel function of the second kind. The primary multiplicity of a heavy hadron produced from the SHM is then given by

where $Z(C-C_{j},B-B_{j})/Z(C,B)$ (known as chemical factor less than unity for $c$ or $b$ hadrons in a neutral system with $C=B=0$ [38]) characterizes the canonical suppression relative to the grand-canonical statistical production as a result of the exact conservation of $C$ and $B$.

To evaluate the partition function, hadrons to be summed over in the exponential of Eq. (3) are divided into five categories: light hadrons made of up, down and strange quarks with $C_{j}=B_{j}=0$, $N_{cj}=N_{bj}=0$ and $N_{sj}=0,1,2,3$; open $c$ hadrons with $C_{j}=\pm 1$, $N_{cj}=1$, $B_{j}=0$, $N_{bj}=0$ and $N_{sj}=0,1,2$; charmonia with $C_{j}=0$, $N_{cj}=2$, $B_{j}=0$, $N_{bj}=0$ and $N_{sj}=0$; open $b$ hadrons with $C_{j}=0$, $N_{cj}=0$, $B_{j}=\pm 1$, $N_{bj}=1$ and $N_{sj}=0,1,2$; $B_{c}$ mesons with $C_{j}=\pm 1$, $N_{cj}=1$, $B_{j}=\pm 1$, $N_{bj}=1$ and $N_{sj}=0$; bottomonia with $C_{j}=0$, $N_{cj}=0$, $B_{j}=0$, $N_{bj}=2$ and $N_{sj}=0$. We have safely neglected heavy baryons that contain two $c$ or $b$ quarks since their masses are much larger than other open $c$ or $b$ hadrons. Clearly the summations over light hadrons, charmonia and bottomonia factorize out of the exponential and for the purpose of computing the heavy hadrons’ multiplicity, the evaluation of the partition function reduces to

Now the summation in the exponential of Eq. (3) is restricted to be over positively charged $c$ and $b$ hadrons, such as $D^{0}$, $D_{s}^{+}$, $\Lambda_{c}^{+}$, $B^{0}$, $B_{s}^{0}$, $\bar{\Lambda}_{b}^{0}$ and $B_{c}^{+}$. In Eqs. (3) and (3), the strangeness fugacity is taken to be $\gamma_{s}=1$ (i.e., strangeness chemical equilibrium is reached down to the hadronization in Pb-Pb collisions). But the $c$ and $b$ fugacities $\gamma_{c}$ and $\gamma_{b}$ should be self-consistently determined from the balance equations that stipulate $c$ and $b$ conservation through hadronization:

where the summations run over the primary multiplicities $\langle N_{j}\rangle$ (cf. Eq. (3)) of positively charged open $c$/$b$ hadrons ($oc^{+}$ and $ob^{+}$, respectively), charmonia/bottomonia ($hc$ and $hb$, respectively) and $B_{c}^{+}$ mesons. Since the majority ($\sim 99$%) of $c$ content is carried by open $c$ hadrons, the $c$ fugacity $\gamma_{c}$ is essentially governed by the first term on the right hand side of Eq. (5); the same is true for $\gamma_{b}$.

To perform realistic calculations for $\sqrt{s_{\rm NN}}=5.02$ TeV Pb-Pb collisions, we focus on a fireball corresponding to a rapidity slice of one unit at mid-rapidity with vanishing net $c$ and $b$ number ($C=B=0$). We first determine the input $c$- and $b$-quark multiplicities $dN_{c\bar{c}}/dy$ and $dN_{b\bar{b}}/dy$ using the cross section $d\sigma_{c\bar{c}}/dy=1.165~{}{\rm mb}$ [42] and $d\sigma_{b\bar{b}}/dy=39.3~{}{\rm\mu b}$ [8, 43] in $\sqrt{s}=5.02$ TeV $pp$ collisions at mid-rapidity and the thickness functions ($\langle T_{\rm AA}\rangle$) for different centralities [44]. We have also applied 10-20% reduction to $c$ quark multiplicity because of nuclear shadowing effect [33, 45]. For the input open $c$- and $b$-hadron spectrum, we use the listings by the particle data group (PDG) [5], augmented by additional states that have been predicted by relativistic quark models (RQM) [46, 47] but not yet measured. The RQM predictions in particular feature many additional $c$ and $b$ baryons that are “missing” in the current PDG listings [47]. The inclusion of these “missing” states has proved pivotal for the successful reproduction of $c$ and $b$ baryon-to-meson ratios in $pp$ collisions via SHM calculations [8, 38, 48] that demonstrate significant enhancements relative to $e^{+}e^{-}$ collisions. For $B_{c}$ mesons, we have included all 14 states below the open $BD$ threshold [4]. The fireball volume is obtained by scaling the one $V_{\Delta y=1}=4997~{}{\rm fm^{3}}$ determined from SHM for light hadrons in the most central 0-10% centrality [36] to other centralities using the measured charged-particle multiplicities [49]. Finally, the hadronization temperature is taken to be $T_{H}=170$ TeV, which is higher than the pseudocritical chiral transition temperature ($T_{\chi}\sim 155$ MeV [50, 51]) but seems more appropriate for the hadronization (confinement) transition of heavy hadrons [8, 48, 52, 53]. We’ve also checked that lowering down the hadronization temperature by $\sim 10$ MeV does not cause significant change for the $b$-hadron production fractions [8].

The computed $c$ and $b$ fugacities for different centralities, alongside the corresponding fireball volumes and $c$- and $b$-quark multiplicities, are summarized in Table 1. The $\gamma_{c}$ exhibits a mild decrease from central to semicentral collisions, which then turns into an increase toward peripheral collisions, similar to the system size dependence of $\gamma_{c}$ found in [36]. In contrast, the $\gamma_{b}$ increases monotonously from central to peripheral centrality bins. Once the primary multiplicities from statistical production are computed, the total production yields of ground state open $b$ hadrons and $B_{c}^{-}$ mesons are obtained from the sum of the direct one and the feeddown contributions from excited states

where the branching fractions $\mathcal{B}$ for the strong decays of excited open $b$ hadrons have been estimated from a ${}^{3}P_{0}$ model [8] and those for the strong or radiative decays of excited $B_{c}$ mesons to the ground state $B_{c}^{-}$ are all taken to be 100% [1]. The resulting yields per unit rapidity of ground state $B^{-}$, $B_{c}^{-}$ and other open $b$ mesons and baryons, alongside the production yield ratios of $B_{c}^{-}$ to $B^{-}$ and to the total $b\bar{b}$ multiplicity (i.e. $B_{c}^{-}$’s production fraction $f_{c}$) in different centralities, are displayed in Table 2. While the production ratios between open $b$ hadrons (e.g., $\Lambda_{b}^{0}/B^{-}\sim 0.5$) previously computed from grand-canonical SHM [8] are well reproduced (all open $b$ hadrons containing a single $b$-quark suffer from the common canonical suppression not affecting their ratios), the $B_{c}^{-}/B^{-}$ reaches $\sim 0.043$ in central Pb-Pb collisions, amounting to a factor of $\sim 5$ enhancement with respect to the value in $pp$ collisions (cf. Sec. 2). Accordingly the $B_{c}^{-}$’s production fraction $f_{c}\sim 0.013$ in central collisions also represents a significant enhancement of the similar magnitude relative to $pp$ collisions, which is the prominent consequence of recombination production (implemented in a statistical fashion here) of $B_{c}$ mesons from abundant and highly thermalized $c$ and $b$ quarks in the QGP.

We also note that, while in central and semicentral collisions the exact conservation of $c$ number does not add to the canonical suppression for the production of $B_{c}$ mesons, strict conservation of $c$ number becomes important toward peripheral (especially 60-80%) collisions where the $dN_{c\bar{c}}/dy$ reduces to the order of unity (canonical limit). We’ve numerical checked that, should only the $b$ number conservation be implemented in the statistical partition function (cf. Eq. (3)), the $B_{c}^{-}$ yield in the 60-80% centrality would be $\sim 70\%$ greater than the one shown in Table 2. This mechanism of canonical suppression owing to strict conservation of the $c$ number largely explains the decrease of $B_{c}^{-}$’s production fraction toward peripheral collisions as indicated in Table 2.

We now convert the computed $B_{c}^{-}$’s statistical production yield into a $p_{T}$ distribution. Our strategy is to distribute the SHM yield according to the shape of the $p_{T}$ spectrum of $B_{c}^{-}$ calculated from the resonance recombination [15, 33, 54] of realistically transported $c$ and $b$ quarks. As an illustration, we focus on the 20-40% centrality which is used as a proxy for the minimum bias (0-90%) Pb-Pb collisions for which the CMS analysis was performed [31].

While recombination as a hadronization mechanism for $c$ or $b$ quarks in the QGP dominates at low $p_{T}$, it yields to vacuum-like fragmentation at high $p_{T}$ [8, 15]. It has been determined that more than 90% of $b$ quarks diffusing in the QGP in central and semicentral Pb-Pb collisions undergo recombination at hadronization [8]. Therefore, the absolute $p_{T}$ differential yield for the $B_{c}^{-}$’s recombination component is obtained by normalizing the $p_{T}$ spectrum of $B_{c}^{-}$ calculated from the resonance recombination of $c$ and $b$ quarks to its SHM yield corrected for the $b$-quark integrated recombination probability which is $\sim 92\%$ for the 20-40% centrality [8]. The resonance recombination model (RRM) [15, 33, 54] adopted here conserves 4-momentum and satisfies correct equilibrium limit [15]. The momentum distribution given by RRM for the recombination $B_{c}^{-}$ reads

where $f_{b}$ and $f_{\bar{c}}$ are the transported phase space distributions of $b$ and $\bar{c}$ quarks, and $v_{\rm rel}$ is their relative velocity. In the current $2\rightarrow 1$ formulation, recombination proceeds via a resonance cross section $\sigma_{B_{c}}(s)$ for $b+\bar{c}\rightarrow B_{c}^{-}$, taken to be of the Breit-Wigner form with the vacuum $B_{c}^{-}$ mass and a width $\Gamma_{B_{c}}\simeq 100$ MeV (variations by a factor of $\sim 2$ have practically no effect on the shape of the $B_{c}^{-}$’s $p_{T}$ spectrum), while the $C_{B_{c}}$ ensures normalization to the statistical production yield. The phase space distributions $f_{b}$ and $f_{\bar{c}}$ are constructed on the hadronization hypersurface from Langevin simulations of $b$- and $c$-quark diffusion in the QGP down to $T_{H}=170$ MeV and have been constrained by open $b$ and $c$ hadrons’ observables [8, 15]. Since the RRM has been carried out here with the full quark phase space distributions, space-momentum correlations between $b$ and $\bar{c}$ quarks built up from Langevin simulations have been incorporated, which help produce a significantly harder meson spectrum compared to recombination only in the momentum space [15, 33].

The high $p_{T}$ $B_{c}$ mesons are dominantly produced from fragmentation of $b$ quarks (the fragmentation probabilities for $\bar{c}\rightarrow B_{c}^{-}$ is two orders of magnitude smaller than those for $b\rightarrow B_{c}^{-}$ [55]). To calculate the absolute $p_{T}$ differential yield for the $B_{c}^{-}$’s fragmentation component, we take the fragmenting $b$-quark spectrum after Langevin diffusion (i.e., the $b$ quarks left over from recombination, whose integrated number accounts for $\sim 8\%$ of total $b$ quarks [8]), and simulate its fragmentation into $B_{c}^{-}$ using the fragmentation function [55]

where $z=p_{T}/p_{t}$ is the fraction of $B_{c}^{-}$’s $p_{T}$ relative to the parent $b$-quark’s transverse momentum $p_{t}$. In Eq. (4), the parameter $r$ controlling the slope of the resulting spectrum is tuned such that the latter resembles the $B_{c}^{-}$’s $p_{T}$ spectrum in $pp$ collisions as much as possible especially at high $p_{T}$, while the normalization constant $N$ is tuned to ensure that the resulting $B_{c}^{-}$’s integrated yield should account for the same fraction ($f_{c}=0.00258\pm 0.00062$, cf. Sec. 2) of the total fragmenting $b$ quarks as in $pp$ collisions.

The resulting $p_{T}$ differential yield $dN/dp_{T}dy$ for the recombination and fragmentation components, alongside their sum (i.e., the total one), are shown in Fig. 2(a). The widths of bands indicate uncertainties from the spread in $f_{c}$ (or equivalently $\mathcal{B}(B_{c}^{-}\rightarrow J/\psi\mu^{-}\bar{\nu})$, cf. Sec. 2) in $pp$ collisions (for the fragmentation component) or from the 10-20% reduction of $c$-quark multiplicity due to the shadowing effect (for the recombination component, cf. Sec. 3). An immediate observation is that the momentum reach of $B_{c}^{-}$’s recombination extends to rather high $p_{T}\sim 18$ GeV, beyond which the fragmentation takes over. In Fig. 2(a), we have also plotted the $B^{-}$’s $p_{T}$ differential yield for comparison which was previously computed from the $b$-quark transport approach [8] and is reproduced here; the integrated $B^{-}$ yield agrees with the value shown in Table 2. The total $p_{T}$ differential yields for $B_{c}^{-}$ and $B^{-}$ shown in Fig. 2(a) are converted into their $p_{T}$ dependent ratio $B_{c}^{-}/B^{-}$ which is shown in Fig. 2(b) and compared to its counterpart in $pp$ collisions [6] (cf. Sec. 2). The $B_{c}^{-}/B^{-}$ is enhanced by a factor of $\sim 5$ at low $p_{T}$ relative to $pp$ collisions, but tends to the value in the latter at very high $p_{T}$ within uncertainties. Accordingly, the $B_{c}^{-}$’s nuclear modification factor defined as

reaches the value of $\sim 5$-6 at $p_{T}<5$ GeV, in contrast to the $B^{-}$’s $R_{\rm AA}$ (quoted from [8]) around unity in the same $p_{T}$ range, as shown in Fig. 2(c). The enhancement of $B_{c}^{-}$’s production in Pb-Pb collisions as characterized by its $R_{\rm AA}$ above unity persists till $p_{T}\sim 10$ GeV, at which the $B^{-}$’s production is already suppressed by a factor of $\sim 2$ ($R_{\rm AA}\simeq 0.48$) relative to $pp$ collisions. The width of the band for $B_{c}^{-}$’s $R_{\rm AA}$ at low $p_{T}$ indicates uncertainties due to shadowing of $c$’s participating in the statistical recombination as well as the spread of the $pp$ reference spectrum, but is dominated by the latter. At high $p_{T}$ where the fragmentation component dominates, the uncertainties from the spread in $f_{c}$ (or $\mathcal{B}(B_{c}^{-}\rightarrow J/\psi\mu^{-}\bar{\nu})$) in the numerator spectrum and the denominator reference cancel, making the band become narrower.

Same as for the $p_{T}$ differential cross section of $B_{c}$ mesons in $pp$ collisions (cf. Fig. 1(b)), the CMS data [31] for their $p_{T}$ differential yield in minimum bias (0-90%) Pb-Pb collisions were presented for $B_{c}^{-}+B_{c}^{+}$, multiplied by the pertinent branching fractions of the two-step decay process $B_{c}\rightarrow(J/\psi\rightarrow\mu\mu)\mu\nu$ and normalized by the thickness function, in two bins of the three-muons’ transverse momentum $p_{T}^{\mu\mu\mu}$. To compare to the CMS data, the $dN/dp_{T}dy$ for $B_{c}^{-}$ shown in Fig. 2(a) is first doubled, then multiplied by $\mathcal{B}(B_{c}^{-}\rightarrow J/\psi\mu^{-}\bar{\nu})=1.95\%\pm 0.46\%$ [6] and $\mathcal{B}(J/\psi\rightarrow\mu\mu)=5.96\%$ [5], divided by the corresponding $\langle T_{\rm AA}\rangle$ [44], and re-plotted as a function of $p_{T}^{\mu\mu\mu}=0.85\cdot p_{T}^{B_{c}}$ [31] (cf. Sec. 2) in Fig. 3(a). The two data points measured by CMS are fairly described within uncertainties. This normalized $p_{T}$ differential yield is then divided by the $p_{T}$ differential cross section in $pp$ collisions shown in Fig. 1(b) to obtain the $B_{c}$’s $R_{\rm AA}$ as a function of $p_{T}^{\mu\mu\mu}$, which is shown in Fig. 3(b) and compared to the corresponding CMS data. While the data point in the lower $p_{T}$ bin is fairly reproduced within theoretical and experimental uncertainties, the suppression indicated in the higher $p_{T}$ bin is somewhat overestimated by our calculation.

Finally we calculate $B_{c}$’s integrated $R_{\rm AA}$ in the full $p_{T}$ range and in the $7<p_{T}<13$ GeV interval (corresponding to $6<p_{T}^{\mu\mu\mu}<11$ GeV in the CMS data [31]) as a function of participant numbers. As shown in Fig. 4, the integrated $R_{\rm AA}$ in the full $p_{T}$ range reaches $\sim 5$ in central and semicentral collisions but gradually drops off toward peripheral collisions, in line with the centrality dependence of $B_{c}^{-}$’s statistical production fraction shown in Table 2. The enhancement embodied in the calculated integrated $R_{\rm AA}$ in the higher $p_{T}$ interval is much milder, amounting to a factor of $\sim 50\%$ in central collisions and vanishing in peripheral collisions, fully consistent with the CMS data points in two centrality bins within uncertainties.

In this work, we have investigated the recombination production of $B_{c}$ mesons in Pb-Pb collisions at the LHC energy using a statistical hadronization approach. By treating $B_{c}^{-}$ as a member of the family of open $b$ hadrons that contain a single $b$-quark and implementing the strict conservation of $b$ and $c$ numbers, we were able to make quantitative predictions for $B_{c}^{-}$’s production fraction relative to the total $b{\bar{b}}$ multiplicity and its relative production to $B^{-}$ mesons in the presence of QGP, both demonstrating a factor of $\sim 5$ (3) enhancement in central (peripheral) Pb-Pb collisions with respect to $pp$ collisions. The statistical production yield of $B_{c}^{-}$ mesons is converted into a $p_{T}$ differential distribution according to the shape of $B_{c}^{-}$’s $p_{T}$ spectrum computed from resonance recombination of realistically transported $b$ and $\bar{c}$ quarks in the QGP. This was supplemented with the fragmentation of high $p_{T}$ $b$ quarks and enabled us to compute the $p_{T}$ dependent nuclear modification factors for $B_{c}^{-}$ mesons that reach $\sim$5-6 at low $p_{T}$ in central collisions, with the reference spectrum carefully constructed from previously determined differential cross section of $B^{-}$ mesons and the experimentally measured $B_{c}^{-}/B^{-}$ ratio. Pertinent data measured by the CMS experiment, including $B_{c}^{-}$’s normalized $p_{T}$ differential yield, nuclear modification factors in two $p_{T}$ intervals in the minimum bias collisions and the integrated ones in two centrality bins, can all be fairly described in our approach within theoretical and experimental uncertainties. Our study provides strong support for the formation of $B_{c}$ mesons through (statistical) recombination of abundant and highly thermalized $c$ and $b$ quarks in the deconfined QGP created in Pb-Pb collisions at the LHC energy.

Acknowledgments: This work was supported by the National Natural Science Foundation of China (NSFC) under Grant No.12075122.