Higgs boson decays to $B_c$ meson in the fragmentation-function approach

In the paper, we present a calculation of the decay widths for the Higgs boson decays to the $B_c$, $B_c^*$, $B_c(2^1S_0)$ and $B_c^*(2^3S_1)$ mesons using the fragmentation-function approach. In the calculation, the fragmentation functions up to order $\alpha_s^3$ based on the nonrelativistic QCD factorization theory are used, and the decay widths for $H\to Q+X$ and $H \to g+X$ at the partonic level are calculated up to order $\alpha_s$. The large logarithms of $m_H^2/m_{Bc}^2$ are resummed up to next-to-leading logarithmic accuracy by solving the evolution equations for the running quark masses and the fragmentation functions. Compared to the leading-order decay widths based on the nonrelativistic QCD approach, the decay widths based on the fragmentation-function approach that include the higher-order QCD corrections are reduced significantly. Our numerical results show that there are about $1.2\times 10^5$ $B_c$ events via the Higgs decays to be produced at the HL-LHC with $3ab^{-1}$, and about $1.6\times 10^6$ $B_c$ events via the Higgs decays to be produced at the HE-LHC with $15ab^{-1}$.


I. INTRODUCTION
The discovery of the Higgs boson at the LHC [1,2] in 2012 is an important breakthrough in our understanding of fundamental interactions. After that, an important task is to accurately study the properties of the Higgs boson, including the Higgs couplings to the fundamental fermions and the gauge bosons, as well as the Higgs self-coupling, and test whether these couplings are completely consistent with those predicted by the Standard Model (SM). Any measurement that deviates from the SM prediction may be a signal of new physics.
The LHC has achieved great success in discovering the Higgs boson, and has studied some couplings of the Higgs boson, e.g., the couplings to heavy vector bosons [3][4][5] and the charged fermions of the third generation [6][7][8][9][10][11][12][13][14]. However, the precision of these measurements is restricted due to the limited Higgs events and complicated hadronic background. After a period of shut down, the LHC has just upgraded to Run 3. During Run 3, more data will be collected than the first two runs combined. Futhermore, the LHC is planned to upgrade to the High-luminosity LHC (HL-LHC) and the Highenergy LHC (HE-LHC) after Run 3. At the HL-LHC ( √ s = 14 TeV), with an integrated luminosity of 3 ab −1 , about 1.6 × 10 8 Higgs boson events will be produced; At the HE-LHC ( √ s = 27 TeV), with an integrated luminosity of 15 ab −1 , about 2.2 × 10 9 Higgs boson events will be produced [15]. In addition, several lepton colliders are under consideration, e.g., the Circular Electron-Positron Collider (CEPC) [16], the International Linear Collider (ILC) [17], and the e + e − Future Circular Collider (FCCee) [18], and the Muon Collider [19,20]. One of the advantages of lepton colliders is that the background is * zhengxc@cqu.edu.cn † wuxg@cqu.edu.cn ‡ zhanxj@cqu.edu.cn § wanggy@cqu.edu.cn ¶ liht@cqu.edu.cn clean, thus they are suitable for the precision measurements of the properties of the Higgs boson. With these collider platforms, some rare decays of the Higgs boson, such as the Higgs decays to quarkonium, may be measured [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]. These rare decays can be used to determine the magnitude of the Yukawa couplings of the Higgs boson to the heavy quarks, and they have distinguished signals to be detected at the high-luminosity or high-energy colliders. The B c meson carries two different heavy flavors and provides a unique bound-state system for testing the SM. As a combination, in Ref. [39], the authors studied the Higgs decays to the B c meson at the leading order (LO) accuracy, and they found that about 1.4 × 10 5 B c events can be produced through the Higgs decays at the HL-LHC. In addition to studying the Higgs properties, this decay process can also be used to study the production mechanism of the B c meson. Thus, it is attractive to present a more precise study on this decay process. In the present paper, we devote ourselves to reanalyzing this decay process with higher accuracy in the fragmentation-function approach.
There are large logarithms of the form ln(m 2 H /m 2 Bc ) in the perturbative series of the decay width of the Higgs boson into a B c meson, which come from two sources: the renormalization of the Yukawa couplings and the emission of the collinear gluons. These large logarithms may spoil the convergence of the perturbative expansion, thus it is important to sum them to all orders (in α s ) in the calculation. It is noted that under the fragmentationfunction approach, the large logarithms from these two sources can be resummed simultaneously. More explicitly, the large logarithms from the renormalization of the Yukawa couplings can be resummed by using the running quark masses for the heavy (b and c) quarks [40,41]; while the large logarithms from the collinear gluon emission can be resummed through solving the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations for the fragmentation functions of the B c production [42,43]. In this paper, we will resum these large logarithms up to next-to-leading logarithmic (NLL) accuracy under the fragmentation-function approach.
Because of carrying two different heavy flavors, the excited states of the B c meson below the BD threshold 1 will decay to the ground state B c meson with almost 100% probability through electromagnetic or strong interaction. Thus, the excited states are important sources of the ground state B c meson. Furthermore, the production of these excited states via the Higgs boson decays is also interesting by itself. Therefore, besides the decay width for the Higgs boson into the ground state B c meson, we will also calculate the decay widths for the Higgs boson decays into the S-wave excited states, e.g. B * c , B c (2 1 S 0 ), and B * c (2 3 S 1 ). The paper is organized as follows. In Sec.II, we present useful formulas for the decay width of the Higgs boson into the B c meson under the fragmentation-function approach. In Sec.III, numerical results and discussions are presented. Section IV is reserved as a summary.

II. CALCULATION FORMALISM
In this section, we present useful formulas for the considered decay width under the fragmentation-function approach. For simplicity, we only give the formulas for the ground state B c meson, the formulas for the excited states [i.e., B * c , B c (2 1 S 0 ), and B * c (2 3 S 1 )] are similar to the ground state B c case.
Under the fragmentation-function approach, the differential decay width for the decay channel H → B c + X can be written as (1) where dΓ H→i+X (y, µ F )/dy stands for the differential decay width of H → i + X at the partonic level, D i→Bc (z/y, µ F ) stands for the fragmentation function for a parton i into the B c meson, z = 2p Bc · P H /m 2 H denotes the energy fraction carried by the B c meson from the Higgs boson, µ F denotes the factorization scale, and the sum extends over the patron species.
The decay widths for the Higgs boson into a parton can be calculated through perturbation theory. Up to now, the decay width for the Higgs boson into bottom quarks has been calculated up to order α 4 s [40,41,[45][46][47][48][49][50][51][52]. However, the expressions for the differential decay widths dΓ/dz of the Higgs boson into a quark or gluon are not given in those references. We calculate the differential decay widths for H → Q + X and H → g + X up to order α s in this work. In the calculation, we neglect the quark mass in the amplitudes and phase space integrals except the quark mass in the Yukawa coupling. This approximation will only lead to an error of O(m 2 Q /m 2 H ). Then, we have where Q can be a quark or an antiquark, N c = 3 is the number of quark colors, C F = (N 2 c − 1)/(2N c ) is the quadratic Casimir operator, G F is the Fermi constant, and m Q (µ R ) is the running quark mass defined in the modified-minimal-subtraction scheme (MS). The expressions of the LO splitting functions are 1 The Bc excited states above the BD threshold will decay mainly into a pair of B and D mesons [44].
The expressions of C i (y) functions in Eqs. (2) and (3) are C g (y) = C F 1 + (1 − y) 2 y [2ln y + ln(1 − y)] + y . (7) In the calculation, the ultraviolet (UV) divergences are removed by renormalization, and the renormalization of the quark mass is carried out in the usual MS scheme. Besides the UV divergences, there are infrared (IR) soft and collinear divergences appearing in the virtual and real corrections. The IR soft divergences are canceled between the virtual and the real corrections, while the IR collinear divergences (which should be absorbed into the bare fragmentation functions) are subtracted according to the MS scheme. To avoid large logarithms appearing in dΓ H→i+X /dz, we will set the renormalization and factorization scales as µ R = µ F = m H in the following calculation.
The fragmentation functions for a parton into the B c meson can be calculated based on the nonrelativistic QCD (NRQCD) factorization theory [53], i.e., where d i→cb[n] (z, µ F ) is the short-distance coefficient (SDC) for the (cb)[n] pair production, which can be calculated through perturbative QCD. O Bc (n) is the longdistance matrix element (LDME) for the transition of a (cb)[n] pair into the B c meson, which can be estimated through phenomenological models, e.g. the potential models. The sum extends over intermediate Fock states. In the lowest nonrelativistic approximation, only the Fock state n = 1 S  [54,55]. They extracted the fragmentation functions from the processes Z → B c (B * c ) + b +c by taking the approximation of m Bc /m Z → 0. Their results were confirmed by the subsequent calculations in Refs. [56,57] using different methods. For a long time, the NLO fragmentation functions for the B c production are absent. Recently, with the development of loop-diagram calculation techniques, the NLO fragmentation functions forb → B c (B * c ) and c → B c (B * c ) has been given in Ref. [58]. Furthermore, the fragmentation functions for g → B c (B * c ), which start at order α 3 s , have been obtained in Refs. [59,60]. In this work, we will adopt the fragmentation functions up to order α 3 s obtained in Refs. [58,59]. In order to avoid large logarithms appearing in dΓ H→i+X /dz, we have set the factorization scale as µ F = m H in Eq.(1). However, the large logarithms of m 2 H /m 2 Bc will appear in the fragmentation functions. To resum the large logarithms in the fragmentation functions, we first calculate the fragmentation functions up to order α 3 s with initial scales µ R0 = µ F 0 = m b + m c using the codes developed in our previous works [58,59]. Then the fragmentation functions with µ F = m H can be obtained through solving the DGLAP equations, i.e., where P ji (y, α s (µ F )) are the splitting functions, which can be expanded in powers of α s : The LO splitting functions for Q → Q and Q → g have been given in Eqs. (4) and (5), and the LO splitting functions for g → Q and g → g are where C A = N c and T F = 1/2. The NLO corrections to these splitting functions have been obtained in Refs. [61][62][63][64][65], they are too lengthy to be replicated here.
It is nontrivial to solve these integro-differential equations. In this work, we adopt the program FFEVOL [66] to solve the DGLAP equations numerically. In solving the DGLAP equations, the NLO fragmentation functions at µ F 0 = m b + m c are used as the boundary conditions and the NLO spiliting functions are used as the evolution kernel. After the evolution of the fragmentation functions from the initial factorization scale

III. NUMERICAL RESULTS AND DISCUSSION
To do the numerical calculation, the input parameters are adopted as follows: where the values for the Fermi constant and the masses are taken from the Particle Data Group (PDG) [67]. R 1S (0) and R 2S (0) are the radial wave functions at the origin for the (cb) bound states, which are taken from the calculation based on the Buchmüller-Tye potential model [68]. The running masses at µ R = m H can be obtained through solving the renormalization-group equation, i.e., where the first two coefficients [69,70] are and n f is the number of active flavors. We solve this renormalization-group equation by using the Mathematica package RunDec [71], and only the first two coefficients of the right hand of Eq. (14) are preserved (i.e., the obtained running masses at µ R = m H reach the NLL accuracy). Then we have In the calculation of the fragmentation functions in Ref. [58], the heavy quark masses are renormalized in the on-shell (OS) scheme. The OS (pole) masses for the heavy quarks can be obtained from the MS masses [72][73][74][75][76], and we have For the strong coupling constant, we adopt the two-loop formula, i.e., In the fragmentation-function approach, some terms which are suppressed by powers of m 2 Bc /m 2 H are neglected. In order to see the magnitude of those neglected higher-power (in m 2 Bc /m 2 H ) contributions, we compare the decay widths calculated by the "direct" NRQCD and the fragmentation-function approaches. Here, we only present the comparison at the LO level.
The differential decay width for the decay H → B c +X under the (direct) NRQCD approach [53] can be written as dΓ H→Bc +X = n dΓ H→(cb)[n]+X O Bc (n) , (19) where dΓ H→(cb)[n]+X is the SDC for the (cb)[n] pair production, which can be calculated through perturbation theory. At the LO level, there are four Feynman diagrams responsible for the decay H → B c +X. The details about the LO calculation based on the direct NRQCD approach can be found in Ref. [39]. In the calculation, we adopt the package FeynArts [77] to generate the Feynman diagrams and the amplitudes, and use the FeynCalc [78,79]  The partial decay widths for H → B c + X and H → B * c + X under the direct NRQCD approach and the fragmentation-function approach 2 are presented in Tables I and II. Here, for consistency, the quark masses are taken as the corresponding pole masses and the strong coupling is taken as α s (m b + m c ) = 0.204 under the two approaches. In the tables, the contributions from thē b-fragmentation and the c-fragmentation as well as the total contribution are presented explicitly. For the direct NRQCD approach, the first two Feynman diagrams in Fig.1 are responsible for theb-fragmentation contribution, while the last two Feynman diagrams are responsible for the c-fragmentation contribution 3 . The interference contribution comes from the interference of the first two Feynman diagrams and the last two Feynman diagrams, and the fragmentation-function approach can not give the interference contribution.
From Tables I and II, we can see that the differences between the decay widths under two approaches are very small for both theb-fragmentation and the cfragmentation. 4 Moreover, the interference contributions are very small in both the B c and B * c cases. There are two reasons for the small interference contributions. First, the interference contributions come from the interference of theb-fragmentation diagrams and the c-fragmentation diagrams. The c-fragmentation diagrams are suppressed compared to theb-fragmentation diagrams. This leads to the interference contributions are suppressed compared to theb-fragmentation contributions. Second, the dominant contribution of theb-fragmentation diagrams comes from the phase-space region where the relative velocity of the B c (B * c ) meson and thec-quark is small and back-to-back with the b-quark, while the dominant contribution of the c-fragmentation diagrams comes from the phase-space region where the relative velocity of the B c (B * c ) meson and the b-quark is small and back-to-back with thec-quark, i.e., the dominant contributions of thē b-fragmentation diagrams and the c-fragmentation diagrams come from different phase-space regions, which further suppresses the contributions of the interference. These two reasons lead to the very small interference contributions. Therefore, the higher-power terms neglected in the fragmentation-function approach are very small in the two decay processes, the fragmentation-function approach can give a good approximation to the direct NRQCD approach.
The results show that the c-fragmentation contributions are about 3 orders of magnitude smaller than the correspondingb-fragmentation contributions. There are two reasons for the very small c-fragmentation contributions: One is that the magnitude of the Yukawa coupling of Hcc is smaller than that of Hbb, and the other is that the fragmentation probability of c → B c is smaller than that ofb → B c .
The differential decay widths dΓ/dz of H → B c + X and H → B * c +X under the direct NRQCD approach and the fragmentation-function approach are shown in Figs. 2 and 3. From the figures, we can see that the curves from the two approaches are very close. The differences constitute two gauge-invariant subgroups, respectively. Hence, in an arbitrary covariant gauge (e.g. the Feynman gauge), the first (last) two Feynman diagrams in Fig.1 should be simultaneity taken into account for theb-fragmentation (c-fragmentation) contribution. 4 The results show that the difference between the two approaches for B * c is larger than that for Bc. This indicates that the accuracy of the fragmentation-function approximation depends not only on the large energy scale involved in the process (e.g. m H in H → Bc(B * c ) + X), but also on the quantum number of the produced state. This point is also shown in the process gg → Bc(B * c ) + b +c [80]. between the two approaches are relatively small at large z values, and relatively large at small z values. The reason is that the fragmentation contribution mainly comes from the phase-space regions with large z value. More explicitly, the contribution ofb-fragmentation comes from the phase-space region where the B c (B * c )-c system has small invariant mass and large momentum. The contribution of c-fragmentation comes from the phase-space region where the B c (B * c )-b system has small invariant mass and large momentum. In these phase-space regions, the fragmentation-function approach can well describe the process with the B c (B * c ) production. In the phase-space regions where the momentum of the B c (B * c ) meson is very small, the non-fragmentation contribution becomes relatively important. approach. The upper one shows the contributions ofb-fragmentation and cfragmentation respectively, the lower one shows total contribution. In order to put the results from theb-fragmentation and the c-fragmentation into one figure, the curves for the c-fragmentation are multiplied by a factor of 50.

FIG. 4. Feynman diagrams induced by the triangle top-quark loop.
From the comparison of the LO results under the two approaches given in the last subsection, we found that the fragmentation-function approach (without the resummation of large logarithms) can give a good approximation to the direct NRQCD approach, i.e., the power corrections in the fragmentation-function approach are negligi-ble at the LO level. At the NLO, there are nonfragmentation Feynman diagrams induced by a triangle top-quark loop as shown in Fig.4. Compared with the fragmentation contribution, the contribution from these Feynman diagrams is suppressed by powers of m 2 Bc /m 2 H but enhanced by the Htt coupling. Therefore, before giving the results from the fragmentation-function approach up to the NLL accuracy, it is important to see how much these triangle top-loop diagrams contribute. In fact, the authors of Ref. [39] have calculated the contribution from the triangle top-quark loop diagrams. They obtained very strange results, i.e., the triangle top-loop contribution in the B * c case is one order of magnitude smaller than that in the B c case. To further illustrate the reason of the smallness of the triangle top-loop contribution in the B * c case, we recalculate the contribution from the triangle top-loop diagrams here. Furthermore, the authors of Ref. [39] found that the interference contribution between Fig.4 and Fig.1 is very small for both the B c and B * c cases. However, we can not conclude from the small interference contribution that the contribution from the square of Fig.4 is also small. Because the topologies of Fig.4 and Fig.1 are significantly different, they may be dominated by different phase-space regions. Hence, in addition to the contribution from the interference between the Fig.4 and Fig.1, we also calculate the contribution from the square of Fig.4.

Contributions
Bc B * c Interference of Fig.1 Fig.1 and Fig.4, and "Square" denotes the contribution coming from the square of Fig.4.
In Table III and Fig.5, the contributions from the tri-angle top-quark loop are presented 5 . In the calculation, the top-quark mass is taken as m t = 172.8 GeV [67], and the other parameters are taken the same values as those in the last subsection. From Table III, we can see that the interference contribution in the B * c case is one order of magnitude smaller than that in the B c case. This can be understood by the distribution shown in Fig.5. The distribution of the interference contribution in the B * c case is negative at large z values, this indicates that there is a large cancellation between the contributions from different phase-space regions. Furthermore, we can also see that the contributions (the interference contribution as well as the squared contribution) from the triangle top-quark loop are very small compared with the LO contributions.   S0) and B * c (2 3 S1), where the contributions from three fragmentation channels and the triangle top-loop contribution are given explicitly.
From the analysis presented in the above two subsections, we believe that the neglected power suppressed terms in the fragmentation-function approach are small in the decays H → B c + X and H → B * c + X even up to the NLO level. In this subsection, we present the decay widths calculated based on the the fragmentationfunction approach up to the NLL accuracy.
In Table IV, the partial decay widths for the Higgs decays to B c , B * c , B c (2 1 S 0 ) and B * c (2 3 S 1 ) are presented, where the contributions from the different fragmentation channels and the contribution from the triangle top-quark loop are given explicitly. The calculation method for these fragmentation contributions have been described in Sec.II, i.e., the factorization scale is taken as µ F = m H in Eq.(1) and the fragmentation functions at µ F = m H are obtained through the DGLAP evolution from the initial factorization scale µ F 0 = m b + m c . Hence, the large logarithms of m 2 H /(m b + m c ) 2 which arise from the collinear gluon emission and the renormalization of the Yukawa couplings have been resummed in the numerical results for these fragmentation contributions. The input parameters have been given at the beginning of this section.
From Table IV, we can see that the decay widths up to NLL accuracy are significantly smaller than the corresponding LO decay widths shown in Tables I and II. 5 Adopting the same input parameters as in Ref. [39], we are able to reproduce the numerical results for the contribution of the interference between the triangle top-quark loop diagrams and the LO diagrams given in Table IV of Ref. [39]. FIG. 6. The differential decay width dΓ/dz for H → Bc + X, where the contributions from the three fragmentation channels and the triangle top-quark loop are shown explicitly. In order to put the results into one figure, the curve for the cfragmentation is multiplied by a factor of 50, and the curve for the g-fragmentation is multiplied by a factor of - 50. This indicates that the higher-order corrections, especially the large logarithmic terms, are important in these decay processes. Therefore, the resummation of large logarithms should be taken into consideration for giving high-precision predictions. The contributions from the c-fragmentation and the g-fragmentation are very small compared to theb-fragmentation contribution. Moreover, the g-fragmentation contribution is negative for these processes. In Figs.6, 7, 8 and 9, the differential decay widths dΓ/dz for the B c , B * c , B c (2 1 S 0 ) and B * c (2 3 S 1 ) states are shown. In the figures, the contributions of the different fragmentation channels and the triangle top-quark loop to dΓ/dz are shown explicitly. From these figures, we can see that the differential decay widths are dominated by theb-fragmentation contribution for all the z values. The c-fragmentation and the g-fragmentation contributions mainly come from the small z values. However, even for small z values, these contributions are also small compared to theb-fragmentation contribution. FIG. 9. The differential decay width dΓ/dz for H → B * c (2 3 S1) + X, where the contributions from the three fragmentation channels and the triangle top-quark loop are shown explicitly. In order to put the results into one figure, the curve for the c-fragmentation is multiplied by a factor of 50, and the curve for the g-fragmentation is multiplied by a factor of -50.

D. Uncertainty analysis
In this subsection, we will estimate the theoretical uncertainties for these partial decay widths. The main uncertainty sources for these decay widths include the factorization and renormalization scales, the heavy quark masses, the Higgs boson mass, and the cb radial wave functions at the origin. The dependence of the decay widths on the Higgs boson mass mainly comes from the partonic decay widths dΓ H→i+X (y, µ F )/dy(i = Q, g), which contain a global factor m H . The uncertainty for the world average value of the Higgs boson mass given by the PDG is about 0.2 GeV [67], which is only about 0.2% of the Higgs boson mass. Therefore, the uncertainties for the decay widths caused by the Higgs boson mass are only about 0.2% of their central values. Since the uncertainties caused by the Higgs mass are very small, we will neglect this uncertainty source in the following uncertainty estimaton.
There are several factorization and renormalization scales involved in the calculation based on the fragmentation-function approach: the initial (lower) factorization and renormalization scales (µ F 0 and µ R0 ) for the initial fragmentation functions; the final (upper) factorization and renormalization scales (µ F and µ R ). In the calculation presented in the last subsection, these scales are set as µ F 0 = µ R0 = m b + m c and µ F = µ R = m H . In the uncertainty estimation, we vary them by a fac- where the first uncertainty is caused by the lower factorization and renormalization scales, and the second uncertainty is caused by the upper factorization and renormalization scales.
where the first uncertainty is caused by the bottom quark mass, while the second uncertainty is caused by the charm quark mass. For the wave functions at the origin, we have adopted the values based on the Buchmüller-Tye potential model given in Ref. [68]. However, the authors of Ref. [68] did not give an error estimate to the wave functions. In order to give an estimate to the uncertainties from the wave functions, we take the values based on the Buchmüller-Tye potential model as the central values, while take the values based on the logarithmic potential and the Cornell potential as the boundary values for the wave functions. The values for the radial wave functions based on the three potential models are shown in Table V. We obtain the uncertainties caused by the wave functions as Adding the uncertainties from different sources in quadrature, we obtain the total theoretical uncertainties for these partial decay widths as follows:

IV. SUMMARY
In the present paper, we have calculated the partial decay widths for the Higgs boson decays to the B c , B * c , B c (2 1 S 0 ), and B * c (2 3 S 1 ) mesons based on the fragmentation-function approach. The decay widths and the differential distributions are obtained, and the theoretical uncertainties for the decay widths are estimated. In the calculation, the fragmentation functions up to order α 3 s for the B c (B * c ) production calculated in the previous works are used as the initial fragmentation functions. The large logarithms that arise from the renormalization of the Yukawa couplings and the collinear gluon emissions are resummed up to NLL accuracy through solving the evolution equations of the running heavy-quark masses and the fragmentation functions. After including these higher-order contributions, the decay widths are reduced significantly compared to the LO predictions.
In order to have a feel on the size of the higher-power terms (in m 2 Bc /m 2 H ) that neglected in the fragmentationfunction approach, we compared the decay widths under the direct NRQCD approach with those under the fragmentation-function approach at the LO level. The results show that the decay widths under the two approaches are very close to each other, i.e., those higherpower terms are very small at the LO level. We also studied the contributions induced by the triangle topquark loop, which are enhanced by the Htt coupling. The results show that these contributions are very small compared to the fragmentation contributions. Moreover, we found that the interference between the triangle topquark loop diagrams and the LO Feynman diagrams for H → B * c + X has a strong cancellation between different phase-space regions. This leads to a much smaller contribution from the triangle top-quark loop in the B * c case than that in the B c case.
Since the B c excited states below the BD threshold will decay to the ground state B c with almost 100% probability, the total decay width for the Higgs boson decay to the B c meson is approximately equal to the sum of the decay widths for the Higgs boson decays to the cb meson states below the BD threshold. Adding the decay widths for the S-wave states (B c , B * c , B c (2 1 S 0 ), and B * c (2 3 S 1 )) shown in Eq. (23), and using the Γ H ≈ 3.2 MeV [68], we obtain that the total branching fraction for the Higgs boson decays into the B c meson is about 7.47×10 −4 . According to the total branching fraction, we estimate that there are about 1.2 × 10 5 B c events will be produced via the Higgs boson decays at the HL-LHC with 3 ab −1 , and about 1.6 × 10 6 B c events will be produced via the Higgs boson decays at the HE-LHC with 15 ab −1 . Therefore, these decay processes may be studied at the future HL-LHC and HE-LHC, and provide a complementary method for measuring the bottom-quark Yukawa coupling.