The relativistic corrections to the fragmentation functions for a heavy quark to $B_c$ and $B_c^{*}$

TL;DR

This paper calculates relativistic corrections to the fragmentation functions of a heavy quark into B_c and B_c* mesons using NRQCD, including their evolution via DGLAP equations.

Contribution

It provides the first detailed computation of relativistic effects on heavy quark fragmentation functions into B_c mesons within the NRQCD framework.

Findings

Relativistic corrections significantly affect fragmentation functions.

Presented DGLAP evolution for both singlet and non-singlet cases.

Results improve the precision of heavy quark fragmentation modeling.

Abstract

In this paper, we compute the relativistic corrections to the fragmentation functions (FFs) of a heavy quark to $B_c$ and $B_c^*$ within the framework of non-relativistic QCD (NRQCD) factorization. The non-singlet and singlet DGLAP evolution are also presented.

Equations

$$E_p\frac{d{\sigma }_{A+B\to H+X}}{d^{3}p}\approx \sum _f\int \frac{dz}{z^2}{D}_{f\to H}(z;\mu )E_c\frac{d{\sigma }_{A+B\to f(p_c)+X}}{d^{3}p}\left(p_c=\frac{1}{z}p\right),$$

$${\mu }^2\frac{d}{d{\mu }^2}{D}_{f\to H}(z;\mu )$$

$$D_{b\to B_c^{(\ast )}}(z;\mu )$$

$$W(x^{-})\equiv \mathrm{P}\exp \left(+i{g}_s{\int }_{x^{-}}^{\infty }ds{A}^{+}(sn)\right),$$

$$D_{b\to B_c{(}^{2s+1}{S}_J)}(z,\mu )$$

$$⟨{\mathcal{O}}_1^{B_c}{(}^1{S}_0)⟩$$

$$⟨{\mathcal{O}}_8^{B_c}{(}^1{S}_0)⟩$$

$$⟨{\mathcal{P}}_1^{B_c}{(}^1{S}_0)⟩$$

$$⟨{\mathcal{O}}_1^{B_c^{\ast }}{(}^3{S}_1)⟩$$

$$⟨{\mathcal{O}}_8^{B_c^{\ast }}{(}^3{S}_1)⟩$$

$$⟨{\mathcal{P}}_1^{B_c^{\ast }}{(}^3{S}_1)⟩$$

$$⟨{\mathcal{O}}_1^{B_c}{(}^1{S}_0)⟩\simeq \frac{N_c}{2\pi }|R_S(0){|}^2,⟨{\mathcal{O}}_1^{B_c}{(}^3{S}_1)⟩\simeq (d-1)\frac{N_c}{2\pi }|R_S(0){|}^2.$$

$$|[b\bar{c}{(}^{2s+1}{S}_J^{[1]})](p)⟩=\sum _{a,b=1}^{N_c}\frac{{\delta }_{ab}}{\sqrt{N_c}}|b^a(p_b),{\bar{c}}^b(p_c)⟩,$$

$$|[b\bar{c}{(}^{2s+1}{S}_J^{[8]})](p)⟩=\sum _{a,b=1}^{N_c}\frac{T_{ab}^A}{\sqrt{T_F}}|b^a(p_b),{\bar{c}}^b(p_c)⟩,$$

$$p^{\mu }$$

$$p_b$$

$$D_{b\to [b\bar{c}]}(z;\mu )=\frac{z}{4N_c}\sum _{\text{colors\&spins}}\int \frac{d{q}^{+}{d}^2{q}_{\perp }}{(2\pi {)}^{3}2q^{+}}\theta (q^{+})\delta (p^{+}/z-p^{+}-q^{+}){\left|\mathcal{M}\right|}^2,$$

$${\left|\mathcal{M}\right|}^2=\mathrm{Tr}\left[n/⟨0|W(0)b(0)|[b\bar{c}](p),c(q)⟩⟨[b\bar{c},c(q)](p),c(q)|\bar{b}(0)W^{†}(0)|0⟩\right],$$

$$⟨b^a(p_b),{\bar{c}}^b(p_c),c^c(q)|{\left[\bar{b}(0)W^{†}(0)\right]}_d|0⟩$$

$$k^{\mu }=\left(\frac{{\mathbf{k}}^2}{2m_b},\mathbf{k}\right),{\tilde{k}}^{\mu }=\left(\frac{{\mathbf{k}}^2}{2m_c},-\mathbf{k}\right).$$

$$v(p_c)\bar{u}(p_b)$$

$$⟨{\mathcal{O}}_1^{b\bar{c}}{(}^1{S}_0)⟩$$

$$⟨{\mathcal{O}}_8^{b\bar{c}}{(}^1{S}_0)⟩$$

$$⟨{\mathcal{O}}_1^{b\bar{c}}{(}^3{S}_1)⟩$$

$$⟨{\mathcal{O}}_8^{b\bar{c}}{(}^3{S}_1)⟩$$

$$⟨{\mathcal{P}}_1^{b\bar{c}}(n)⟩$$

$$\sum _{\lambda =0,\pm 1}{\epsilon }^{\ast \mu }(p,\lambda ){\epsilon }^{\nu }(p,\lambda )\to -g^{\mu \nu }+\frac{p^{\mu }{p}^{\nu }}{p^2}=-g^{\mu \nu }+v^{\mu }{v}^{\nu },$$

$$\sum _{\lambda =\pm 1}{\epsilon }^{\ast \mu }(p,\lambda ){\epsilon }^{\nu }(p,\lambda )\to -g^{\mu \nu }+\frac{p^{\mu }{n}^{\nu }+p^{\nu }{n}^{\mu }}{p^{+}}-\frac{p^2{n}^{\mu }{n}^{\nu }}{(p^{+}{)}^2}=-g^{\mu \nu }+\frac{v^{\mu }{n}^{\nu }+v^{\nu }{n}^{\mu }}{v^{+}}-\frac{n^{\mu }{n}^{\nu }}{(v^{+}{)}^2},$$

$$D_{b\to b\bar{c}{(}^1{S}_0^{[1]})}$$

$$D_{b\to b\bar{c}{(}^3{S}_1^{[1]})}$$

$$D_{b\to b\bar{c}{(}^1{S}_0^{[8]})}$$

$$D_{b\to b\bar{c}{(}^3{S}_1^{[8]})}$$

$$d_{{}^1{S}_0^{[1]}}^{(0)}(r,z)$$

$$d_{{}^3{S}_1^{[1]}}^{(0)}(r,z)$$

$$d_{{}^1{S}_0^{[8]}}^{(0)}(r,z)$$

$$d_{{}^3{S}_1^{[8]}}^{(0)}(r,z)$$

$$d_{{}^1{S}_0^{[1]}}^{(2)}(r,z)$$

$$d_{{}^3{S}_1^{[1]}}^{(2)}(r,z)$$

Full text

††preprint: To be published in Chinese Physics C

I Introduction

One of main fields for precision examination of the perturbative Quantum Chromodynamics (QCD) is the study of the hadron production with the momentum $p$ ($p^{2}=m_{H}^{2}$) at the large transverse momentum $p_{T}$ where $p_{T}^{2}>>m_{H}^{2}$. At the leading power of $1/p_{T}$, the fragmentation mechanism may be dominant so that the differential cross-section can be factorized as (for reviews, see Altarelli:1981ax ; Collins:1989gx ; Ellis:1991qj ; Collins:2011zzd )

[TABLE]

where parton-level differential cross-section $d\sigma_{A+B\to f(p_{c})+X}$ contains the short-distance dynamics, and $D_{f\to H}(z;\mu)$ is the so-called single parton fragmentation function (FF) of the hadron $H$ from parton $f$ parametrizing the universal hadronization effects Collins:1981uk ; Collins:1981uw . Generally, the fragmentation functions are non-perturbative, but their renormalization scale-dependences are perturbatively calculable, and can be described by the famous Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations which also govern the renormalization group (RG) running of the parton distribution functions (PDFs)

[TABLE]

where $P_{f\leftarrow f^{\prime}}(z,\alpha_{s}(\mu))$ are the celebrated Altarelli-Parisi splitting functions DGLAP ; LOAPkernel .

For the FFs of quarkonia, which are different from the FFs for the light mesons which rely completely on the dynamics in the non-perturbative regime of QCD, one believes that they can be further factorized into the products of the perturbatively calculable parts and non-perturbative behavior of the wave-functions of quarkonia at origin, due to the nature of quarkonium as a non-relativistic bound state of heavy quark and anti-quark. The standard theoretical tool to deal with the heavy quark bound state system is the NRQCD factorization Bodwin:1994jh , in which all information of hadronization of quarkonium are encoded in the NRQCD matrix elements. In the literature Braaten:1996pv ; Braaten:1993rw ; Braaten:1994kd ; Braaten:1996rp ; Ma:1994zt ; Sang:2009zz ; Hao:2009fa , various single-parton FFs quarkonia have been calculated. Recently, some of such fragmentations for quarkonia have been calculated to the higher order of $\alpha_{s}$ and $v$ analytically or numerically Ma:2013yla ; Artoisenet:2014lpa ; Zhang:2017xoj ; Feng:2017cjk ; Artoisenet:2018dbs ; Feng:2018ulg ; Zhang:2018mlo ; Bodwin:2003wh ; Bodwin:2002hg ; Ma:2015lka by implementing the state-of-art multi-loop calculation techniques.

The $B_{c}$ mesons are unique mesons which consist of two different flavor heavy quarks. Their production at high energy collision is a hot topic of perturbative QCD since 1990s Chang:1991bp ; Chang:1992bb ; Chang:1992pt ; Braaten:1993jn ; Cheung:1993pk ; Chen:1993ii ; Yuan:1994hn ; Chang:2001pm ; Chang:2014jca ; Zheng:2015ixa ; Zheng:2017xgj ; Zheng:2019gnb ; Qiao:2011yk ; Qiao:2011zc ; Jiang:2015jma ; Jiang:2015pah . Among these works, a lot of efforts have been paid on the FFs of $B_{c}$ mesons, and most of calculations were done within the framework of NRQCD factorization similar to corresponding calculations for quarkonia. It is worth noting that the first NLO calculation of the FF for a heavy-quark to $B_{c}^{(*)}$ meson, which is quite non-trivial, has been carried out in Zheng:2019gnb .

In this paper, we compute the relativistic corrections to the FFs of a heavy-quark to $B_{c}$ and $B_{c}^{*}$ at the leading order of $\alpha_{s}$. These corrections are as equally important as the NLO QCD radiative corrections due to the NRQCD power counting rule. Similar results had appeared in the appendix of Sang:2009zz as a by-product of calculation of FFs for quarkonia. The new things in this paper, that distinguish from results of Sang:2009zz , are: 1) we consider a different scheme to calculate the relativistic corrections to S-wave heavy-quark bound system, in which the mass of $B_{c}^{(*)}$ meson appears as an overall factor different from the scheme adopted in Sang:2009zz that the binding energy is treated as a part of non-relativistic corrections; 2) we consider the contribution from the color-octet NRQCD operators; 3) we calculate the FF for a heavy quark to the transversely polarized $B_{c}^{*}$; 4) we further investigate the DGLAP evolution effects.

The paper is organized in the following way. We first present the definition of the FFs for $b$ quark to $B_{c}^{(*)}$ following Collins and Soper, and the desired NRQCD factorization formula in Sect.II. In Sect.III, we illustrate the matching procedure, present the explicit expressions for the short-distance distributions, and compare them with those in literature. In Sect.IV, we present the non-singlet and singlet DGLAP evolution of the obtained FFs. Finally, we summarize in Sect.V.

II Fragmentation Function and its NRQCD facotrization

In this section, we briefly review the definition of the fragmentation function which was given by Collins and SoperCollins:1981uk ; Collins:1981uw , and present the NRQCD factorization formula for FFs of $b$-quark to $B_{c}^{(*)}$.

We adopt the following notations for the decompositions of momenta: the 4-velocity of $B_{c}$ is $v^{\mu}$ with $v^{2}=1$. We also use the same notation $v$ for the non-relativistic expansion parameter, which is typical size of the relative velocity of quark and anti-quark inside a $B_{c}$ meson. In the context, one should not confuse these two. We also introduce two light-like vectors $n^{\mu}$ and $\bar{n}^{\mu}$ such that $n^{2}=\bar{n}^{2}=0$ and $\bar{n}\cdot n=1$. Any 4-vector $a^{\mu}$ can be decomposed as $a^{\mu}=a^{+}\bar{n}^{\mu}+a^{-}{n}^{\mu}+a_{\perp}^{\mu}$($a^{+}\equiv n\cdot a$ and $a^{-}\equiv\bar{n}\cdot a$) with $n\cdot a_{\perp}=\bar{n}\cdot a_{\perp}=0$. Thus, $a\cdot b=a^{+}b^{-}+a^{-}b^{+}+a_{\perp}\cdot b_{\perp}$. For convenience, we set $v^{\mu}=v^{+}\bar{n}^{\mu}+v^{-}n^{\mu}$. Apparently $v^{+}v^{-}=1/2$.

The $d$-dimensional fragmentation function for $b$ quark to $B_{c}^{(*)}$ meson is defined in the way given by Collins and Soper in Collins:1981uk ; Collins:1981uw ,

[TABLE]

where $b(x)$ is the $b$ quark field in QCD, and the Wilson line along a light-like path is defined as

[TABLE]

in which $g_{s}$ is the SU(3) gauge coupling and $A_{\mu}(x)\equiv A^{A}_{\mu}(x)T^{A}$ ($T^{A}$s are the generators of SU(3) group in the fundamental representation). Since in this paper we never encounter the situation that we need the extra dimension to regulate the ultraviolet or infrared divergences, we always set $d=4$.

Since $B_{c}^{(*)}$ can be regarded as a non-relativistic bound-state of $b$ and anti-$c$ quark, the contribution around the heavy-quark mass may be further factorized. Therefore, up to ${\cal O}(v^{2})$, we may have the factorization formula within the NRQCD factorization framework,

[TABLE]

where the mass scale $M\equiv m_{b}+m_{c}$ is introduced to balance the mass dimensions so that the perturbatively calculable short-distance distributions $d(z,\mu)$ are dimensionless, and matrix-elements of the singlet and octet NRQCD operators are defined as

[TABLE]

where $\chi_{c}$ and $\psi_{b}$ are the two-component effective fields in NRQCD for $\bar{c}$ and $b$ quark, respectively, $\sigma^{i}$ is the $i$-th Pauli matrix, and $\mathbf{D}=\mathbf{\nabla}-ig_{s}\mathbf{A}$ is the covariant derivative. Note that all the $B_{c}^{(*)}$ states defined in (6-11) are non-relativistically normalized, meanwhile such states are relativistically normalized in the definition of FF in (II). The factor of $2m_{B_{c}}$ in the factorization formula (II) is just from the relativistic normalization condition for single particle states.

By using the vacuum saturation approximation, at the leading order of $\alpha_{s}$, the matrix elements $\langle\mathcal{O}_{1}\rangle$ can be related to the radial wave function of S-wave $B_{c}$ meson at the origin $R_{S}(0)$ in the color-singlet model through

[TABLE]

However, there is no similar simple relation for the matrix-elements $\langle\mathcal{O}_{8}\rangle$.

III Computation of the Fragmentation Functions

III.1 Matching procedure

In the practical computation of the short-distance distributions $d(z,\mu)$, we replace the $B_{c}^{(*)}$ with a color-singlet quark pair $[b\bar{c}(^{2s+1}S_{J}^{[1]})]$,

[TABLE]

or a color-octet quark pair $[b\bar{c}(^{2s+1}S_{J}^{[8]})]$

[TABLE]

where $a,b$ are color indices, $T_{F}=1/2$ is the Dynkin index of SU(3) group in the fundamental representation. And we set the momenta as

[TABLE]

where $k^{\mu}$ and $\tilde{k}^{\mu}$ are the residual momenta of $b$-quark and $\bar{c}$-quark, respectively.

At leading order of $\alpha_{s}$, we have

[TABLE]

with

[TABLE]

where $q^{2}=m_{c}^{2}$. Schematically, it has the similar NRQCD factorization formula as in (5) by just replacing the $B_{c}^{(*)}$ state by corresponding $b\bar{c}$ pair state.

In the light-cone gauge,

[TABLE]

where we have used various equations of motion to simplify the spin structures. We also checked the expression of this matrix-element in the covariant gauge, which is the same as the above expression in the light-cone gauge as it should be.

Since the FF is defined in a Lorenz invariant way, we can do the calculation in any reference frame. For convenience, we choose the rest frame of the $b\bar{c}$ pair, in which we have $v^{\mu}=(1,\mathbf{0})$ and

[TABLE]

Generally, one applies the replacements

[TABLE]

to project out the spin-singlet and spin-triplet parts, respectively. Here $\varepsilon^{\mu}$ denotes the polarization vector for the spin-triplet quark pair, and the energies $E_{b}=\sqrt{m_{b}^{2}+\mathbf{k}^{2}}$ and $E_{c}=\sqrt{m_{c}^{2}+\mathbf{k}^{2}}$. After expanding the resulting matrix-element (18) to the second order in the momentum ${\bf k}$, one can extract the S-wave contributions by neglecting the first order terms in ${\bf k}$ (which contribute to $P$-wave state only) and making the replacement $k^{i}k^{j}\to{\bf k}^{2}\delta^{ij}/3$. This leads to the corresponding matrix-elements of NRQCD operators at leading order of $\alpha_{s}$,

[TABLE]

It is worth noting that we calculate the FF for $B_{c}^{*}$ by using the replacement rules for the summation of polarization vectors

[TABLE]

for the unpolarized case, and

[TABLE]

for the transversely polarized case. We will invoke a superscript T for FF and corresponding short-distance distributions of transversely polarized $B_{c}^{*}$.

Through the matching equations,

[TABLE]

we can extract the short-distance distributions straightforwardly.

III.2 Results of short distance distributions

At the lowest order of $v$,

[TABLE]

in which we denote $r\equiv m_{c}/(m_{b}+m_{c})$ and $\bar{r}\equiv 1-r$ for shortening the expressions. The color factor $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ with $N_{c}=3$. Here we fix $\alpha_{s}$ at $\mu=m_{b}+m_{c}$. $d^{(0)}_{{}^{1}{}S_{0}^{[1]}}(r,z)$ and $d^{(0)}_{{}^{3}{}S_{1}^{[1]}}(r,z)$ agree with the results obtained by Braaten et al in Braaten:1993jn by expressing $\langle O_{1}\rangle$ in term of the Schrödinger wave function at origin.

At the order of $v^{2}$, we get

[TABLE]

and

[TABLE]

For transversely polarized fragmentation function for $B_{c}^{*}$, we have the short-distance distributions as

[TABLE]

and

[TABLE]

In Braaten:1993jn , the authors attribute $\int_{0}^{z}D_{b\to B_{c}^{(*)}}(z)$ to be a kind of fragmentation probability. It is quite interesting to investigate that the ratios among such fragmentation probabilities of $b$ to $B_{c}$, unpolarized $B_{c}^{*}$ and transversely polarized $B_{c}^{*}$, such as

[TABLE]

One can see that such ratios in the limit $r\to 0$ agree with the rough estimation based on arguments of counting the degrees of freedom on polarizations. Even when ${\cal O}(v^{2})$ corrections are taken into account, the above limit of ratios when $r$ goes to zero still hold.

The fragmentation functions $D_{\bar{c}\to B_{c}^{(*)}}$ can be easily obtained by replacing $r$ to $1-r$ in the above results. Namely, one has

[TABLE]

which is just a consequence of a kind of charge charge-conjugation Braaten:1993jn .

III.3 Comparisons and discussions

As mentioned in Sect.I, in Sang:2009zz , the relativistic corrections to the FF have been calculated. Neglecting the color-octet contribution, and taking into account the Gremm-Kapustin relation Gremm:1997dq , that is for $B_{c}^{(*)}$,

[TABLE]

where

[TABLE]

we can derive the same results as in Sang:2009zz .

We note that in many phenomenological applications of NRQCD factorization, the matrix-element $\langle{\cal P}^{H}\rangle$ is treated as a free parameter which is usually not bounded by the Gremm-Kapustin relation. In this sense, it is inconsistent with attributing the binding energy as part of relativistic corrections, which is equivalent to apply the Gremm-Kapustin relation, in the matching procedure. Our scheme to calculate the relativistic corrections by taking out the meson mass as an overall factor is more relevant and consistent in phenomenological applications.

IV Evolution of Fragmentation Function

In this section, we solve the DGLAP evolution equations numerically by using the 4-th order classical Runge-Kutta method. We consider both the non-singlet and singlet evolution. We neglect the contributions from the color-octet operators.

The singlet DGLAP equations read as

[TABLE]

where the Altarali-Parisi kernels $P_{f^{\prime}\leftarrow f}(z/z^{\prime},\alpha_{s}(\mu))$ are

[TABLE]

where $C_{A}=N_{c}$ and $\beta_{0}=11C_{A}/3-4T_{F}n_{f}/3$. The non-singlet DGLAP equation can be obtained from the singlet case by neglecting the off-diagonal Altarali-Parisi kernels.

Here the sliding scale $\mu$ is rather as factorization scale than renormalization scale since eventually the such scale dependence of fragmentation functions will be canceled by the corresponding factorization scale dependence in the parton-level cross-sections or partial decay widths after convolution. Of course, the fragmentation function of heavy-quark to $B_{c}$ will be also renormalization scale dependent. However, such dependence should be canceled by the $\beta$-function dependent terms (mainly from virtual vertex corrections and self-energy corrections) in the higher order radiative corrections to these fragmentation functions, which is not the main concern of this work.

We take the strong coupling $\alpha_{s}(m_{Z})=0.118$ with the leading logarithm renormalization running with $n_{f}=5$, the heavy quark pole masses $m_{b}=4.78$ GeV and $m_{c}=1.67$ GeV, $m_{B_{c}}=6.275$ GeV Tanabashi:2018oca . And we also take the Schrödinger wave function at origin for $B_{c}^{(*)}$ from the potential model in Eichten:1994gt , which is $|R_{S}(0)|^{2}=1.642\,\text{GeV}^{3}$. This implies that

[TABLE]

We use the Gremm-Kapustin relation and heavy-quark spin symmetry to get the values for

[TABLE]

We set the FFs at initial renormalization scale $\mu=m_{b}+m_{c}=6.00$ GeV as what we present at the end of Sect.III. The corresponding plots are depicted in Fig.1. One can see that the relativistic corrections enhance the distribution around $15\%-20\%$.

Since the behavior of the FFs for $b$ to $B_{c}$, unpolarized $B_{c}^{*}$ and transversely polarized $B_{c}^{*}$ are similar, for illustration, we just show the effects of the DGLAP evolution for FFs of $B_{c}$. We show the non-singlet DGLAP evolution in Fig.2. One can see that the evolution effects are sizable, but does not change the behavior of $D_{b\to B_{c}}$ very much. The jump around $z=0$ region is mainly due to the instability of the numerical method. With improving the accuracy of the numerical method, such deficiency would be demmished.

The singlet DGLAP evolution of $D_{b\to B_{c}}$, $D_{\bar{c}\to B_{c}}$ and $D_{g\to B_{c}}$ are shown in Fig. 3, Fig. 4 and Fig. 5, respectively. $D_{g\to B_{c}}(z)$ starts from ${\cal O}(\alpha_{s}^{3})$ which means that $D_{g\to B_{c}}(z)=0$ at ${\cal O}(\alpha_{s}^{2})$, therefore, for the consistency of the initial condition setting for the RG running, here we take $D_{g\to B_{c}}(z)=0$ at $\mu=6$ GeV as the initial value. From these three diagrams, one can observe that

• •

all three FFs at $z=0$ tend to diverge due to the singlet DGLAP evolution;

• •

the value of $D_{\bar{c}\to B_{c}}$ is almost one order of magnitude smaller than $D_{b\to B_{c}}$, which is understandable since it is harder to get a heavier quark from parton fragmentation;

• •

at $\mu=50.0$ GeV and $125.0$ GeV, the values of $D_{g\to B_{c}}$ are comparable with $D_{\bar{c}\to B_{c}}$, which implies that the contribution from gluon fragmentation is needed for more accurate theoretical prediction on $B_{c}$ production at large transverse momentum region.

V Summary and outlook

In this paper, we calculate relativistic corrections to the FFs for a heavy quark to $B_{c}$ and unpolarized and transversely polarized $B_{c}^{*}$ within the framework of NRQCD factorization at the leading order of $\alpha_{s}$. We show that by adopting the Gremm-Kapustin relation, we can reproduce the results in the literature, however we argue that our scheme to count the relativistic corrections is more consistent in the phenomenological applications.

We show the ${\cal O}(v^{2})$ corrections is around $15-20\%$ as we expect, which is as important as the NLO radiative corrections in Zheng:2019gnb . We also show the DGLAP evolution is equally important for phenomenological application. The singlet DGLAP evolution can drastically change the behavior of the FFs of a heavy quark to $B_{c}^{(*)}$ at the small $z$ region. Furthermore, the singlet DGLAP evolution implies that the FFs of gluon to $B_{c}^{(*)}$ could be as important as both relativistic corrections and NLO radiative corrections to the FFs of heavy-quark to $B_{c}^{(*)}$, especially in the production at hadron collider at the large transverse momentum region. Therefore, this FF deserves the further investigation beyond the level of DGLAP evolution.

When the transverse momentum goes down, the so-called double parton fragmentation mechanism becomes more and more important to increase the theoretical predictive power Ma:2013yla . The color-singlet double parton fragmentation function is closely related to the light-cone distribution amplitude (LCDA) which appears mostly in the exclusive production. The LCDAs of $B_{c}$ meson have been studied up to NLO in both $\alpha_{s}$ and $v$-expansion in Xu:2016dgp ; Wang:2017bgv . However, the color-octet double fragmentation functions for $B_{c}^{(*)}$ are still missing, and valuable to be carried out.

Acknowledgement

We thank Yu Jia and Wen-Long Sang for valuable discussions. The work of Deshan Yang is supported in part by the National Natural Science Foundation of China under Grants No. 11275263 and 11635009.