The QCD NLO fragmentation functions for $c$ or $\bar{b}$ quark to $B_c$ or $B_c^*$ meson and their application

TL;DR

This paper derives and numerically computes the QCD NLO fragmentation functions for charm and bottom quarks to Bc mesons, applying them to predict Bc production at Z-factories with improved accuracy.

Contribution

It provides the first derivation of NLO fragmentation functions for c and b quarks to Bc mesons and applies them to collider predictions.

Findings

Fragmentation functions derived up to NLO accuracy.

Numerical results presented with high precision.

Application to Bc production at Z-factories.

Abstract

The fragmentation functions for $c$ or $\bar{b}$ quark to $B_c$ or $B_c^*$ meson are derived up to QCD next-to-leading order (NLO). They are further computed numerically and presented in figures precisely. In order to reach a higher accuracy, we have also tried properly to apply them to estimating the $B_c$ and $B_c^*$ production at a $Z$-factory (an $e^+e^-$ collider running at the energy of the $Z$-boson pole).

Tables

Table 1. Table 1: The fragmentation probability and average value of z 𝑧 z for D b ¯ → B c ( z , μ F ​ 0 = m b + 2 m c , μ R ) D_{\bar{b}\to B_{c}}(z,\mu_{F0}=m_{b}+2m_{c},\mu_{R}) with two typical renormalization scales.

${\mu }_R$	$P\times {10}^4$(LO)	$P\times {10}^4$(NLO)	$⟨z⟩$(LO)	$⟨z⟩$(NLO)
$2m_c$	3.82	3.14	0.68	0.70
$m_b+2m_c$	2.05	2.73	0.68	0.69

Table 2. Table 2: The fragmentation probability and average value of z 𝑧 z for D b ¯ → B c ∗ ( z , μ F ​ 0 = m b + 2 m c , μ R ) D_{\bar{b}\to B^{*}_{c}}(z,\mu_{F0}=m_{b}+2m_{c},\mu_{R}) with two typical renormalization scales.

${\mu }_R$	$P\times {10}^4$(LO)	$P\times {10}^4$(NLO)	$⟨z⟩$(LO)	$⟨z⟩$(NLO)
$2m_c$	5.36	2.91	0.73	0.77
$m_b+2m_c$	2.89	3.25	0.73	0.74

Table 3. Table 3: The fragmentation probability and average value of z 𝑧 z for D c → B c ( z , μ F ​ 0 = 2 m b + m c , μ R ) D_{c\to B_{c}}(z,\mu_{F0}=2m_{b}+m_{c},\mu_{R}) with two typical renormalization scales.

${\mu }_R$	$P\times {10}^6$(LO)	$P\times {10}^6$(NLO)	$⟨z⟩$(LO)	$⟨z⟩$(NLO)
$2m_b$	4.95	8.07	0.51	0.51
$2m_b+m_c$	4.63	7.72	0.51	0.51

Table 4. Table 4: The fragmentation probability and average value of z 𝑧 z for D c → B c ∗ ( z , μ F ​ 0 = 2 m b + m c , μ R ) D_{c\to B^{*}_{c}}(z,\mu_{F0}=2m_{b}+m_{c},\mu_{R}) with two typical renormalization scales.

${\mu }_R$	$P\times {10}^6$(LO)	$P\times {10}^6$(NLO)	$⟨z⟩$(LO)	$⟨z⟩$(NLO)
$2m_b$	4.28	5.75	0.55	0.54
$2m_b+m_c$	4.00	5.57	0.55	0.54

Table 5. Table 5: The total cross sections (in pb pb {\rm pb} ) for the production of B c subscript 𝐵 𝑐 B_{c} and B c ∗ superscript subscript 𝐵 𝑐 B_{c}^{*} mesons at the Z pole. Here, “LO” and “NLO” denote the results from the complete LO and NLO calculations, while “Frag” denotes the results from the fragmentation approach.

States	LO	NLO	Frag
$B_c$	1.76	2.53	2.51
$B_c^{\ast }$	2.46	3.07	2.98

Equations

$$d{\sigma }_{e^{+}{e}^{-}\to B_c+X}=\sum _{n}d{\tilde{\sigma }}_{e^{+}{e}^{-}\to (c\bar{b})[n]+X}⟨{\mathcal{O}}^{B_c}(n)⟩$$

$$d{\sigma }_{e^{+}{e}^{-}\to B_c^{\ast }+X}=\sum _{n^{\prime }}d{\tilde{\sigma }}_{e^{+}{e}^{-}\to (c\bar{b})[n^{\prime }]+X}⟨{\mathcal{O}}^{B_c^{\ast }}(n^{\prime })⟩$$

$$\frac{d{\sigma }_{e^{+}{e}^{-}\to B_c+X}}{dz}=$$

$$D_{i\to B_c}(z,{\mu }_F)=\sum _n{d}_{i\to c\bar{b}[n]}(z,{\mu }_F)⟨{\mathcal{O}}_n^{B_c}⟩,$$

$$\frac{d}{d\ln {\mu }_F^2}{D}_{i\to B_c}(z,{\mu }_F)$$

$$P_{qq}(y)$$

$$P_{gq}(y)$$

$$P_{qg}(y)$$

$$P_{gg}(y)$$

$$D_{Q\to H}(z)=$$

$$D_{\overline{Q}\to H}(z)=$$

$${\mathcal{A}}_1=$$

$${\mathcal{A}}_2=$$

$${\mathcal{A}}_3=$$

$${\mathcal{A}}_4=$$

$$p_{11}=\frac{m_c}{M}{p}_1+q,p_{12}=\frac{m_b}{M}{p}_1-q,$$

$$\Pi =\frac{-\sqrt{M}}{4m_b{m}_c}(/p_{12}-m_b){\gamma }_5(/p_{11}+m_c)$$

$$\Pi =\frac{-\sqrt{M}}{4m_b{m}_c}(/p_{12}-m_b)/ϵ(p_1)(/p_{11}+m_c)$$

$${\Lambda }_1=\frac{1}{\sqrt{3}}\text{1},$$

$${\mathcal{A}}_{\mathrm{LO}}$$

$$d{\varphi }_{\mathrm{LO}}=\frac{\theta (p_2^{+})d{p}_2^{+}}{4\pi p_2^{+}}\frac{d^{d-2}{\text{p}}_{2\perp }}{(2\pi {)}^{d-2}}2\pi \delta (K^{+}-p_1^{+}-p_2^{+}),$$

$$d{\varphi }_{\mathrm{LO}}=$$

$$D_{\bar{b}\to c\bar{b}[n]}^{\mathrm{LO}}(z)=N_{CS}\int d{\varphi }_{\mathrm{LO}}{\mathcal{A}}_{\mathrm{LO}}.$$

$$D_{\bar{b}\to c\bar{b}[n]+\cdots }^{\mathrm{LO}}(z)$$

$$=\frac{C_F^2{\alpha }_s^2(1-z)(4\pi {)}^{ϵ}\Gamma (1+ϵ)}{4N_c{r}_c^{2}z(1-r_{b}z{)}^{4+2ϵ}{M}^{3+2ϵ}}[a_2+a_3\frac{(1+ϵ)z(1-z)}{2(1-r_{b}z{)}^2}$$

$$+a_4\frac{(2+ϵ)(1+ϵ)z^2(1-z{)}^2}{6(1-r_{b}z{)}^4}].$$

$$D_{\bar{b}\to c\bar{b}{[}^1{S}_0^{[1]}]+\cdots }^{\mathrm{LO}}(z)$$

$$=\frac{8{\alpha }_s^{2}z(1-z{)}^2}{81r_c^2(1-r_{b}z{)}^6{M}^3}[6-18(1-2r_c)z$$

$$+(21-74r_c+68r_c^2)z^2-2r_b(6-19r_c+18r_c^2)z^3$$

$$+3r_b^2(1-2r_c+2r_c^2)z^4]\frac{⟨{\mathcal{O}}^{c\bar{b}{[}^1{S}_0^{[1]}]}{(}^1{S}_0^{[1]})⟩}{2N_c},$$

$$D_{\bar{b}\to c\bar{b}{[}^3{S}_1^{[1]}]+\cdots }^{\mathrm{LO}}(z)$$

$$=\frac{8{\alpha }_s^{2}z(1-z{)}^2}{27r_c^2(1-r_{b}z{)}^6{M}^3}[2-2(3-2r_c)z$$

$$+3(3-2r_c+4r_c^2)z^2-2r_b(4-r_c+2r_c^2)z^3$$

$$+r_b^2(3-2r_c+2r_c^2)z^4]\frac{⟨{\mathcal{O}}^{c\bar{b}{[}^3{S}_1^{[1]}]}{(}^3{S}_1^{[1]})⟩}{6N_c},$$

$$⟨{\mathcal{O}}^{c\bar{b}{[}^1{S}_0^{[1]}]}{(}^1{S}_0^{[1]})⟩=2N_c,$$

$$⟨{\mathcal{O}}^{c\bar{b}{[}^3{S}_1^{[1]}]}{(}^3{S}_1^{[1]})⟩=2(d-1)N_c$$

$$⟨{\mathcal{O}}^{B_c}{(}^1{S}_0^{[1]})⟩\approx N_c|R_S(0){|}^2/(2\pi ),$$

$$⟨{\mathcal{O}}^{B_c^{\ast }}{(}^3{S}_1^{[1]})⟩\approx (d-1)N_c|R_S(0){|}^2/(2\pi ),$$

$$D_{\bar{b}\to B_c}^{\mathrm{LO}}(z)$$

Full text

Keywords:

I Introduction

$B_{c}$ and $B_{c}^{*}$ are the ground states of the $c$,$\bar{b}$ binding system with spin [math] and $1$ respectively. Carrying two different heavy flavors, they are unique doubly heavy mesons in the Standard Model. Thus they attract a lot of attentions, particularly after the $B_{c}$ meson was first observedCDF . Their components, being of heavy flavor quarks, move nonrelativistically inside the mesons, so the effective theory —nonrelativistic quantum chromodynamics (NRQCD) nrqcd —is applicable, and the Mandelstam formulation of the Bethe-Salpeter equationmand under the instantaneous approximation also works well.

The production of $B_{c}$ or $B_{c}^{*}$ in $e^{+}e^{-}$ collisions at the Z-boson resonance [i.e., $e^{+}e^{-}\to Z/\gamma\to B_{c}(B_{c}^{*})+X$] under the framework of NRQCD or the Mandelstam formulation under the instantaneous approximation can be factorized as follows doublyhadron ; ybook1 ; ybook2 :

[TABLE]

where $d\tilde{\sigma}$ denotes the cross section for the perturbative production of the two-quark state $(c\bar{b})[n]$ [or $(c\bar{b})[n^{\prime}]$] with proper quantum numbers $n$ (or $n^{\prime}$), which can be calculated using perturbative QCD (pQCD), and the nonperturbative matrix element $\langle{\cal O}^{B_{c}}(n)\rangle$ [or $\langle{\cal O}^{B^{*}_{c}}(n^{\prime})\rangle$] representing the transition probability from the perturbative two-quark state $(c\bar{b})[n]$ [or $(c\bar{b})[n^{\prime}]$] into the hadronic state (a $B_{c}$ or $B^{*}_{c}$ meson) can be related to the wave function at origin of the $(c\bar{b})$ binding system squared in the potential model framework, and can also be calculated using lattice QCD.

Since the $B^{*}_{c}$ meson is similar to the $B_{c}$ meson [the difference is that the spin of the diquark $(c\bar{b})$ inside $B_{c}$ is $S=0$ but the spin of the diquark $(c\bar{b})$ inside $B^{*}_{c}$ is $S=1$], throughout the paper we often use $B_{c}$ to represent both $B_{c}$ and $B_{c}^{*}$ for simplicity.

However, when the center-of-mass energy of a collision is larger than the heavy-quark mass and the terms in ${\cal O}(m^{2}_{Q}/s)$ can be neglected, according to the factorization formulation of pQCD the production can also be calculated in terms of the fragmentation approach:

[TABLE]

where $z\equiv 2p\cdot q/q^{2}$ is the energy fraction (e.g. here $p$ is the momentum of $B_{c}$, and $q$ is the momentum of $e^{+}$ and $e^{-}$ collision), $d\hat{\sigma}_{e^{+}e^{-}\to i+X}$ is the cross section (coefficient function) for the inclusive production of a parton $i$($i=c,\bar{b}$, etc.) and can be calculated using pQCD, $\mu_{F}$ denotes the factorization scale for the production, and $D_{i\to B_{c}}$ is the fragmentation function (FF) from a parton $i$ to a $B_{c}$ meson, which is universal and can be extracted experimentally. The authors of Refs.fraglo1 ; fraglo2 realized that the production is calculable in terms of QCD factorization as shown in Eq.(I) and the leading-order (LO) FFs can be extracted by comparing Eqs.(I) and (2), i.e., the FFs are theoretical calculable, and they were first extracted in Refs.fraglo1 ; fraglo2 . The authors of Ref.comparlo applied the obtained FFs to the production to the QCD leading logarithm approach and made comparisons between their results and those obtained using the complete LO QCD approach, which gives us a understanding of the two approaches.

In order to obtain a better theoretical estimation on $B_{c}$ production, etc., at a Z factorychangch (an $e^{+}e^{-}$ collider running at the energy of the Z-boson pole), we would like to adopt the factorization approach (2) but with the FFs from the $c$ or $\bar{b}$ quark to a $B_{c}$ meson which are of next-to-leading order (NLO), because the NLO QCD calculations are generally more accurate. The NLO FFs cannot be extracted from the complete NLO calculation of the relevant $B_{c}$ production as easily as those for the LO ones, although $B_{c}$ production at a Z factory has been studied using the “complete computation approach”bcnlo . Therefore we must start with the definition given in Ref.Collins to derive them up to the NLO of QCD. In addition, the QCD NLO FFs have many applications, so we would like to derive them precisely here, although the derivation is complicated.

Note that in Refs.braaten ; braaten1 ; YJia ; YQMa the QCD NLO FFs for a gluon to heavy quarkonium were derived, but here the FFs from a quark $i$ to a $B_{c}$ meson involve two heavy quarks of different flavors, so they are quite different from the ones for a gluon to heavy quarkonium.

According to NRQCD, the FFs $D_{i\to B_{c}}$ ($i=c,\bar{b}$), which depict the hadronization and contain nonperturbative effects, can be factorized as follows:

[TABLE]

where the first factor $d_{i\to c\bar{b}[n]}$ denotes a parton $i$ generating a $c\bar{b}$ quark pair with matched quantum number $n$, and being perturbative it can be calculated using pQCD; the factor $\langle{\cal O}_{n}^{B_{c}}\rangle$ denotes the “long-distance matrix elements”, and being nonperturbative they may be related to the wave functions at the origin in the potential model framework or computed using lattice QCD. The nonperturbative factors are reduced to a few long-distance matrix elements $\langle{\cal O}_{n}^{B_{c}}\rangle$ under the required accuracy111The relevant discussions about the accuracy of applying NRQCD to the FFs of heavy quarkonia can be found in Refs.fragnrqcd1 ; fragnrqcd2 , and the conclusions also apply to the FFs of the $B_{c}$ meson.. With the normalization $\int_{0}^{1}dzD_{i\to B_{c}(B^{*}_{c})}(z)=1$, the LO FFs $D_{i\to B_{c}}$(where $i=c,\bar{b}$) were first obtained in Ref.fraglo1 . The LO FFs were extracted from the LO calculations of the processes $Z\to B_{c}+b+\bar{c}$ and $Z\to B_{c}^{*}+b+\bar{c}$ with the approximation $m_{B_{c}}\ll m_{{}_{Z}}$. Subsequent calculationsfraglo2 ; fraglo3 confirmed the results. The LO FFs for the production of the P-wave and D-wave excited states of the $B_{c}$ were calculated in Refs.fraglo4 ; fraglo5 ; fraglo6 . So far there is no NLO calculation for the FFs $D_{i\to B_{c}}$ ($i=c,\bar{b}$). Thus, in the present paper, we devote ourselves to calculating the QCD NLO corrections to $D_{i\to B_{c}}$ (and $D_{i\to B^{*}_{c}}$).

Since the FFs $D_{i\to B_{c}}(z,\mu_{F})$(where $\mu_{F}$ is the factorization energy) generally contain terms like ${\rm ln}(\mu_{F}/m_{Q})$, in order to properly take into account the possible large-logarithm terms the FFs $D_{i\to B_{c}}(z,\mu_{F})$[$\mu_{F}={\cal O}(\sqrt{s})$] will be obtained by solving the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations dglap1 ; dglap2 ; dglap3 with the NLO QCD FFs $D_{i\to B_{c}}(z,\mu_{F0})$[$\mu_{F0}={\cal O}(m_{Q})$] being the “initial FFs”,

[TABLE]

where $P_{ji}(y,\alpha_{s}(\mu_{F}))$ are splitting functions for parton $i$ into parton $j$222In fact, here they are of LO.:

[TABLE]

where $C_{F}=4/3,T_{F}=1/2,C_{A}=3$ for QCD and $P_{\bar{q}\bar{q}}$ is equal to $P_{qq}$. Note that in order to focus on the consequences of NLO QCD corrections for FFs, we restrict ourselves to evaluating the evolution of the FFs from $\mu_{F0}$ to $\mu_{F}$ only to leading-logarithm (LL) accuracy so that here the “splitting functions” in Eq.(5) are of leading order.

The paper is organized as follows. Following the Introduction, in Sec.II we present the definition of the FFs which was given by Collins and SoperCollins , and with this definition we calculate the LO FFs for $i\to B_{c}(B^{*}_{c})+\cdots(i=\bar{b},c)$. In Sec.III we describe the adopted method for calculating the virtual and real corrections to the FFs, and how to carry out the renormalization, so as to obtain the “initial FF” $D_{i\to B_{c}(B^{*}_{c})}(z,\mu_{F0})$. Then, we present the numerical results for the FFs $D_{i\to B_{c}}$ and $D_{i\to B^{*}_{c}}$ up to QCD NLO. In Sec.IV we apply the obtained QCD NLO FFs to the production of $e^{+}e^{-}\to B_{c}(B^{*}_{c})+\cdots$ at a Z factory and compare the results with those obtained from the complete QCD NLO calculations. Section V is devoted to discussions and a conclusion.

II The fragmentation functions

II.1 The definition of fragmentation functions

The FFs may be defined as the hadron matrix elements of certain quark-field operators, and the light-cone coordinate is conventionally adopted. In the light-cone coordinate a vector in d-dimensional space-time333In this work, we adopt dimensional regularization with $d=4-2\epsilon$ to regularize UV and IR divergences, and adopt the reading point prescriptiongamma5 to handle $\gamma_{5}$ in $d$ dimensions. is represented as $V^{\mu}=(V^{+},V^{-},V_{T})=((V^{0}+V^{d-1})/\sqrt{2},(V^{0}-V^{d-1})/\sqrt{2},V_{T})$. The gauge-invariant definition of the FFs for a quark $Q$ fragmenting into a hadron $H$ in $d=4-2\epsilon$-dimensional space-time isCollins

[TABLE]

where $\Psi$ is the quark field and $A_{a}^{\mu}$ is the gluon field. ${\cal P}$ denotes path ordering, $t^{a}$ is the color matrix, $z$ is the longitudinal momentum fraction $z=P^{+}/K^{+}$, and $K$ is the momentum of the initial quark $Q$. The FFs are defined in the reference frame where the hadron $H$ carries the momentum $P^{\mu}=(P^{+},P^{-}=m_{H}^{2}/2P^{+},0_{T})$. It is convenient to introduce a light-like vector $n^{\mu}=(0,1,0_{T})$ in the reference frame where the FFs are defined. Then, the plus component of a momentum p can be written as $p^{+}=p\cdot n$, and $z=P\cdot n/K\cdot n$.

The definitionCollins of FFs for an antiquark $\bar{Q}$ into a hadron $H$ is

[TABLE]

Given the Feynman rules and the definition of the FFs (6)-(7), the relevant Feynman diagrams can be drawn. The part to the left of the cut line in the Feynman diagram corresponds to the right part of the definition, and the part to the right of the cut corresponds to the left part of the definition (which is just the complex conjugate of the right part of the definition). Note that for the FFs of an antiquark into a hadron we have the following:

• •

The vertex for a gluon line attached to an eikonal line contributes a factor $ig_{s}n^{\mu}t^{a}_{ij}$, where $\mu$ and $a$ are the Lorentz index and color index of the gluon, respectively.

• •

The eikonal propagator, which carries momentum $q$ flowing from the operator to the cut side, is $i\delta_{ij}/(q\cdot n+i\epsilon)$.

• •

The cut of final-state eikonal line carrying momentum $q$ contributes $2\pi\delta(q\cdot n)$.

An overall factor of $N_{CS}=z^{1-2\epsilon}/8\pi N_{c}$ from the definition should also be taken into account. The Feynman rules of the FFs of a quark into a hadron are the same as those in the antiquark cases except that the color matrix for the eikonal line–gluon vertex should be $t^{a}_{ij}$ instead of $-t^{a}_{ji}$. Thus, given the Feynman diagrams the LO and NLO FFs $D_{i\to B_{c}}$ and $D_{i\to B^{*}_{c}}$ can be derived.

II.2 LO fragmentation functions

To understand the definition (6)-(7) and to present the conventions used in this paper, here we derive the LO FFs, $D_{i\to B_{c}}$ and $D_{i\to B^{*}_{c}}$, where $i=c,\bar{b}$, although they have been obtained in the past using other approachesfraglo1 ; fraglo2 ; fraglo3 .

In this section and the next one we will show the derivations of the FFs $D_{\bar{b}\to B_{c}}$ and $D_{\bar{b}\to B^{*}_{c}}$ from the definition (6)-(7). The FFs $D_{c\to B_{c}}$ and $D_{c\to B^{*}_{c}}$ can be derived in the same way and the results are the same as those for $D_{\bar{b}\to B_{c}}$ and $D_{\bar{b}\to B^{*}_{c}}$ with the replacement $m_{b}\leftrightarrow m_{c}$, so we will not repeat the derivation for them.

According to the factorization (3), as the first step we derive the “FFs” $D_{\bar{b}\to c\bar{b}[n]}$ with the diquark $c\bar{b}$ states with quantum numbers ${}^{1}S^{[1]}_{0}$ and ${}^{3}S^{[1]}_{1}$, where the superscript $[1]$ denotes the color singlet. Then the second step is to derive the FFs for a heavy quark ($\bar{b}$ or $c$) into a $B_{c}$ or $B^{*}_{c}$ meson, where the “free diquark” $c\bar{b}$ state is replaced by the NRQCD matrix element (the wave function at the origin), which depicts QCD nonperturbative effects in the formation of a $B_{c}$ or $B^{*}_{c}$ meson from the relevant diquark state $c\bar{b}[n]$. (In this paper we assume that the QCD NLO matrix element is the same as the QCD LO one.444The matrix element appears as an overall factor, so its correction(s) (if any) can be considered easily.)

Based on the definition (6)-(7), there are four cut diagrams (Fig.1) for the LO FF $D_{\bar{b}\to c\bar{b}[n]}$. The squared Feynman amplitudes, corresponding to the four diagrams with a “cut”, can be written as follows:

[TABLE]

where $p_{11}$ and $p_{12}$ are the momenta of the $c$ quark and $\bar{b}$ quark inside the $c\bar{b}$ pair and

[TABLE]

where $M\approx m_{b}+m_{c}$ is the mass of the $c\bar{b}$ pair. $\Pi$ is the spin projector: for the spin singlet it is

[TABLE]

and for the spin triplet it is

[TABLE]

$\bar{\Pi}$ is defined as $\bar{\Pi}=\gamma^{0}\Pi^{\dagger}\gamma^{0}$. The color-singlet projector is

[TABLE]

where 1 is the unit matrix of the color $SU_{c}(3)$ group. Note that throughout the paper we work in the Feynman gauge.

Having taken traces, the squared amplitudes corresponding to the LO FFs can be written as follows:

[TABLE]

where $r_{c}=m_{c}/M$ and $r_{b}=m_{b}/M$. $s_{1}=(p_{1}+p_{2})^{2}$ is the invariant mass of the lowest (LO) final states $(c\bar{b}+\bar{c})$. The coefficients $a_{i}$ can be found in the Appendix A.

The differential phase space for the LO FFs can be written as

[TABLE]

where the $\delta$ function comes from the cut through the eikonal line. The integration over $p_{2}^{+}$ can be carried out precisely due to the $\delta$ function. The integrand does not depend on the angles of $\textbf{p}_{2\perp}$, so the integration over the angles of $\textbf{p}_{2\perp}$ is trivial, and can be carried out too. Thus, now the differential phase space is reduced to

[TABLE]

The range of $s_{1}$ is from $(M^{2}/z+m_{c}^{2}/(1-z))$ to $\infty$. The LO FFs can be represented as

[TABLE]

The integration over $s_{1}$ can be carried out with Eq.(16). Integrating over $s_{1}$, we obtain

[TABLE]

Setting $d=4$, we obtain

[TABLE]

and

[TABLE]

where the LO FFs for the $(c\bar{b})$ states have been written in the factorization form, and at order $\alpha_{s}^{0}$,

[TABLE]

with the normalization for the NRQCD matrix elements as that in Ref.nrqcd . $\langle{\cal O}^{c\bar{b}[^{1}S_{0}^{[1]}]}(^{1}S_{0}^{[1]})\rangle$ and $\langle{\cal O}^{c\bar{b}[^{3}S_{1}^{[1]}]}(^{3}S_{1}^{[1]})\rangle$ denote the NRQCD matrix elements for the states $(c\bar{b})$.

Thus, the LO FFs for the $B_{c}$ and $B_{c}^{*}$ mesons are obtained by replacing $\langle{\cal O}^{c\bar{b}[^{1}S_{0}^{[1]}]}(^{1}S_{0}^{[1]})\rangle$ and $\langle{\cal O}^{c\bar{b}[^{3}S_{1}^{[1]}]}(^{3}S_{1}^{[1]})\rangle$ with $\langle{\cal O}^{B_{c}}(^{1}S_{0}^{[1]})\rangle$ and $\langle{\cal O}^{B^{*}_{c}}(^{3}S_{1}^{[1]})\rangle$, respectively. The NRQCD matrix elements $\langle{\cal O}^{B_{c}}(^{1}S_{0}^{[1]})\rangle$ and $\langle{\cal O}^{B^{*}_{c}}(^{3}S_{1}^{[1]})\rangle$ can be estimated as follows:

[TABLE]

where $R_{S}(0)$ is the radial wave function at the origin for the $B_{c}(B_{c}^{*})$ meson. Replacing the NRQCD matrix elements in Eqs.(21) and (22) with the NRQCD matrix elements in Eq.(24), we obtain

[TABLE]

and

[TABLE]

The LO FFs $D^{\rm LO}_{\bar{b}\to B_{c}}(z)$ and $D^{\rm LO}_{\bar{b}\to B^{*}_{c}}(z)$ obtained here are exactly the same as those obtained in Refs.fraglo1 ; fraglo2 , although the authors of Refs.fraglo1 ; fraglo2 derived them in a different way.

III QCD NLO corrections to the FFs for a $\bar{b}$ quark to a $B_{c}$ or $B_{c}^{*}$ meson

In this section we will derive the NLO FFs as defined by Eqs.(6) and (7), and divide the derivation of the NLO corrections into virtual corrections, real corrections and renormalization for convenience. Finally we will compute them numerically and present them in figures.

III.1 The virtual NLO corrections

The virtual NLO corrections to the FFs $D_{\bar{b}\to c\bar{b}[n]}$ come from the “cut diagrams” with one loop on either side of the cut. Four typical cut diagrams for the virtual corrections are shown in Fig.2.

There are Coulomb divergences in the conventional matching procedure. These Coulomb divergences may be regularized by a small relative velocity $v$ between the $\bar{b}$ quark and $c$ quark inside the produced $c\bar{b}$ pair. The Coulomb divergences also appear in the virtual corrections to the NRQCD matrix elements $\langle{\cal O}^{c\bar{b}[^{1}S_{0}^{[1]}]}(^{1}S_{0}^{[1]})\rangle$ and $\langle{\cal O}^{c\bar{b}[^{3}S_{1}^{[1]}]}(^{3}S_{1}^{[1]})\rangle$, while the NRQCD short-distance coefficients are free from Coulomb divergences at all. However, in dimensional regularization, we can avoid the divergence and extract the NRQCD short-distance coefficients by using the so-called region methodregion . In the method, one may calculate the contributions from the hard region directly by expanding the relative momentum $q$ of the $c\bar{b}$ pair before performing the loop integration, and under the lowest nonrelativistic approximation one just needs to take $q=0$ before the loop integration. Thus the Coulomb divergences, which come from the potential region, do not appear in the calculations of the FFs for the free $c\bar{b}$ states and the NRQCD matrix elements. With this method, the NRQCD matrix elements $\langle{\cal O}^{c\bar{b}[^{1}S_{0}^{[1]}]}(^{1}S_{0}^{[1]})\rangle$ and $\langle{\cal O}^{c\bar{b}[^{3}S_{1}^{[1]}]}(^{3}S_{1}^{[1]})\rangle$ at NLO are the same as the LO ones.

The squared amplitudes of the virtual corrections can be read off from the virtual-correction cut diagrams with the Feynman rules in Section II. The Dirac traces are carried out using the Mathmatica packages FeynCalcfeyncalc1 ; feyncalc2 and FeynCalcFormlinkformlink . Then Apart[apart ](#bib.bib31) and FIRE[fire ](#bib.bib32) are adopted to do the partial fraction and integration-by-parts (IBP) reduction. After the IBP reduction, all one-loop integrals in amplitudes are reduced to master integrals. The master integrals include the common scalar one-loop integrals(A_{0}$,$B_{0}$and$C_{0}$functions) and the scalar one-loop integrals with one eikonal propagator. The$A_{0}$,$B_{0}$and$C_{0}$ functions are calculated numerically using LoopToolslooptools . The scalar one-loop integrals with one eikonal propagator can be calculated using the method introduced in the Appendix of Ref.braaten .

The differential phase space for the virtual corrections is the same as that for the LO FFs. The virtual corrections to the FFs can be expressed as

[TABLE]

where ${\cal A}_{\rm virtual}$ denotes the squared amplitudes for the virtual corrections.

III.2 The real NLO corrections

The real corrections to the FFs $D_{\bar{b}\to c\bar{b}[n]}$ come from the fragmentation processes in which an additional gluon is emitted in comparison with the corresponding LO ones. We denote the momenta of the initial and final particles as $\bar{b}(K)\to c\bar{b}[n](p_{1})+\bar{c}(p_{2})+g(p_{3})$. The cut diagrams can be obtained from the LO cut diagrams in Fig.1 by adding a gluon line crossing the cut and connecting two of the lines on each side of the cut. Four typical cut diagrams for the real corrections are shown in Fig.3.

The differential phase space for the real corrections to the FFs can be written as

[TABLE]

The real corrections to the FFs can be written as

[TABLE]

where ${\cal A}_{\rm real}$ denotes the squared amplitudes for the real corrections.

There are UV and IR divergences in the real corrections. These divergences come from the phase-space integration over the momentum of the final gluon $p_{3}$, and yield UV and IR poles in $\epsilon$ in dimensional regularization. However, it is impractical to do the phase-space integration for ${\cal A}_{\rm real}$ analytically. We follow the strategy used in Ref.braaten to calculate the real corrections to the FF $D_{g\to\eta_{Q}}$, in order to extract the UV and IR poles. Namely, we construct the subtraction terms ${\cal A}_{S}$ which have the same singularities as ${\cal A}_{\rm real}$ in the phase space, but the subtraction terms are simpler than ${\cal A}_{\rm real}$, and can be analytically integrated out over the phase space. Then the real corrections can be expressed as

[TABLE]

Therefore, the first term on the right-hand side of Eq.(30) is finite and can be calculated directly in four-dimensional space-time.

The UV divergences in the real corrections arise from the integrations over the phase-space region $p_{3\perp}\to\infty$. The IR divergences arise from the regions $p^{+}_{3}\to 0$ and $p_{3}\to 0$. The squared amplitudes for the real corrections can be expressed as

[TABLE]

where the Lorentz-invariant parameters are defined as follows:

[TABLE]

Since we only consider the production of color-singlet S-wave $c\bar{b}$ states, the $1/t_{1}^{2}$ and $1/t_{1}t_{2}$ terms cancel in ${\cal A}_{\rm real}$bcnlo . The coefficients $b_{i},c_{i},d_{i},g$ and $h$ can be obtained from the squared real-correction amplitudes and the results are very lengthy, so we do not present them here. The integrals of the $b_{i}$ ($i=1,2,3$) terms are UV and IR divergent, and yield double poles $\frac{1}{\epsilon_{{}_{\rm UV}}}\frac{1}{\epsilon_{{}_{\rm IR}}}$ or $\frac{1}{\epsilon^{2}_{{}_{\rm IR}}}$. The integrals of the $c_{i}$ terms are UV divergent, and yield a UV pole $\frac{1}{\epsilon_{{}_{\rm UV}}}$. The integrals of the $d_{i}$ terms are IR divergent, and yield a double pole $\frac{1}{\epsilon^{2}_{{}_{\rm IR}}}$. The integrals of the $g$ and $h$ terms are IR divergent, and yield an IR pole $\frac{1}{\epsilon_{{}_{\rm IR}}}$. The term ${\cal A}^{\rm finite}_{\rm real}$ represents the remaining terms in ${\cal A}_{\rm real}$ which do not contribute divergences.

Now the subtraction terms can be constructed as follows:

[TABLE]

where $\tilde{s}$ is defined as

[TABLE]

where

[TABLE]

One can check that the integration of $({\cal A}_{\rm real}-{\cal A}_{\rm S})$ over the phase space is finite in four space-time dimensions.

In order to analytically extract the UV and IR poles in $\epsilon$ in the real corrections, it is better to choose proper phase-space parametrizations for the terms in ${\cal A}_{\rm S}$. Various phase-space parametrizations can be found in Appendix B.

To integrate the subtraction terms that contain $s$, we use the parametrization in Eq.(LABEL:eqb10) for the differential phase space. The expression of the differential phase space in Eq.(LABEL:eqb10) can be decomposed as

[TABLE]

where the prefactor $N_{\rm LO}(p_{1},p_{2})$ is defined as

[TABLE]

and $d\phi_{\rm LO}(p_{1},p_{2})$ is defined as

[TABLE]

$d\phi_{\rm LO}(p_{1},p_{2})$ represents the differential phase space for a $\bar{b}$ quark with longitudinal momentum $yK\cdot n$ to fragment into a $B_{c}$($B^{*}_{c}$) meson with longitudinal momentum $zK\cdot n$ at LO. $N_{\rm LO}(p_{1},p_{2})$ and $d\phi_{\rm LO}(p_{1},p_{2})$ reduce to $N_{\rm LO}$ and $d\phi_{\rm LO}$ respectively, if $y=1$. Then $d\phi^{(3)}(p_{1},p_{2},p_{3})$ can be expressed as

[TABLE]

The range of $y$ is from $z$ to 1, the range of $s_{1}$ is from $[M^{2}/(z/y)+ym_{c}^{2}/(y-z)]$ to $\infty$, and the range of $s$ is from $s_{1}/y$ to $\infty$.

Thus, with Eqs.(36)-(39), we can obtain

[TABLE]

where $Nd\phi_{\rm LO}(p_{1},p_{2})\equiv N_{\rm LO}(p_{1},p_{2})d\phi_{\rm LO}(p_{1},p_{2})$. The remaining integral in this equation does not generate poles in $\epsilon$.

We can also obtain

[TABLE]

The integration over $y$ will diverge if $\epsilon=0$ in the limit $y\to 1$, and contribute an IR pole. In order to extract this IR pole, we use the plus prescription, where

[TABLE]

Inserting Eq.(LABEL:eqplus) into Eq.(41), we obtain

[TABLE]

where $d\phi_{\rm LO}$ is the LO differential phase space given by Eq.(18).

The integration of the terms with the coefficient $c_{2}(s_{1},z,y)$ over $\Omega_{3\perp}$ is not trivial. We first calculate the integration of the vector $p_{3}^{\mu}$ over $\Omega_{3\perp}$. According to Lorentz invariance, this integral can be expressed as

[TABLE]

We can determine the coefficients $A$ and $B$ by contracting both sides of Eq.(44) with $n_{\mu}$ [and contracting with $(p_{1}+p_{2})_{\mu}$]. Then we obtain

[TABLE]

where $\Omega_{\perp}$ is the total transverse solid angle and $\Omega_{\perp}=2\pi^{1-\epsilon}/\Gamma(1-\epsilon)$.

Inserting Eqs.(45)-(46) into Eq.(44) and contracting both sides of Eq.(44) with $p_{1\mu}$, we obtain

[TABLE]

Carrying out the integration over $s$, we obtain

[TABLE]

The method used to extract the poles from the subtraction terms involving $s$ integration can also be used to extract the poles from the integrations over the $s_{2}$ and $s_{3}$ of the subtraction terms.

For the $s_{2}$ integration, we adopt the parametrization in Eq.(LABEL:eqb13). The expression in Eq.(LABEL:eqb13) can also be decomposed to the form of Eq.(36), but the expression for $d\phi^{(3)}(p_{1},p_{2},p_{3})$ becomes

[TABLE]

The ranges of $y$ and $s_{1}$ are the same as above. The range of $s_{2}$ is from $(1-y+r_{b}z)m_{b}^{2}/r_{b}z$ to $\infty$. Then we can readily obtain

[TABLE]

and

[TABLE]

For the subtraction terms involving $s_{3}$, the parametrization in Eq.(LABEL:eqb16) is adopted and the expression for $d\phi^{(3)}(p_{1},p_{2},p_{3})$ in the form of Eq.(36) is

[TABLE]

The ranges of $y$ and $s_{1}$ are the same as above. The range of $s_{3}$ is from $r_{c}(1-r_{b}z)(s_{1}-m_{b}^{2})/(y-r_{b}z)$ to $\infty$. Then we obtain

[TABLE]

and

[TABLE]

To integrate the subtraction terms that contain $t_{1}$, we adopt the parametrization in Eq.(LABEL:eqb23) for the differential phase space. The expression of the differential phase space in Eq.(LABEL:eqb23) can be written as

[TABLE]

where the prefactor $N_{\rm LO}(p_{1},\tilde{p})$ is defined as

[TABLE]

and $d\phi_{\rm LO}(p_{1},\tilde{p})$ is defined as

[TABLE]

Then the expression of $d\tilde{\phi}^{(3)}(p_{1},p_{2},p_{3})$ can be written as

[TABLE]

The range of $\tilde{s}$ is from $[M^{2}/z+m_{c}^{2}/(1-z)]$ to $\infty$, the range of $u$ is from 0 to 1, and the range of $t_{1}$ is from $(1/z-1)M^{2}u$ to $\infty$.

After integrating over $\Omega_{3\perp}$, $t_{1}$ and $u$, we obtain

[TABLE]

[TABLE]

and

[TABLE]

For the subtraction terms that contain $t_{2}$, the parametrization in Eq.(LABEL:eqb26) is adopted and the expression for $d\tilde{\phi}^{(3)}(p_{1},p_{2},p_{3})$ in the form of Eq.(57) is

[TABLE]

The ranges of $\tilde{s}$ and $u$ are the same as above. The range of $t_{2}$ is from $u\,m_{c}^{2}/(1-u)$ to $\infty$. Then we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Since the remaining integrals in these expressions do not generate poles in $\epsilon$, we can expand these expressions in powers of $\epsilon$ before performing the integration.

III.3 The renormalization

There are UV divergences remaining after summing the contributions from virtual and real corrections, while they are removed by renormalization. We adopt the counterterm approach to carry out the renormalization, where the FFs are calculated with the renormalized coupling constant $g_{s}$, the renormalized quark mass $m$, field $\Psi_{r}$555Here the mass $m$ and field $\Psi_{r}$ may be the mass and field of a $b$ quark or $c$ quark., and the renormalized gluon field $A^{\mu}_{r}$. The renormalized quantities are related to their corresponding bare quantities as

[TABLE]

where $Z_{i}=1+\delta Z_{i}$ with $i=g,m,2,3$ are renormalization constants. The quantities $\delta Z_{i}$ are fixed by the precise definitions of the renormalized quantities. The renormalized quark field, quark mass, and gluon field are defined in the on-mass-shell scheme (OS), whereas the renormalized strong coupling constant $g_{s}$ is defined in the modified-minimal-subtraction scheme ($\overline{\rm MS}$). The expressions of the corresponding renormalization constants in this scheme are obtained as follows:

[TABLE]

where $\mu_{R}$ is the renormalization scale, $\beta_{0}=\frac{11}{3}C_{A}-\frac{4}{3}T_{F}n_{f}$ is the one-loop coefficient of the $\beta$ function in QCD, $n_{f}$ is the number of active quark flavors, $\beta^{\prime}_{0}=\frac{11}{3}C_{A}-\frac{4}{3}T_{F}n_{lf}$ and $n_{lf}=3$ is the number of the light-quark flavors.

Then the contribution from these counterterms can be expressed as

[TABLE]

where ${\cal A}_{\rm counter}$ denotes the squared amplitudes for the counterterms from the renormalization of the quark field, the gluon field, the quark mass and the strong coupling.

Obviously the NLO FFs defined as in Ref.Collins by operator products require renormalizationMueller . We carry out the operator renormalization in the $\overline{\rm MS}$ scheme. The expression for the counterterms for the operator products in this scheme is

[TABLE]

where $\mu_{F}$ is the factorization scale for the FFs and $D_{\bar{b}\to c\bar{b}[n]}^{\rm LO}(z)$ denotes the LO FFs in $d$-dimensional space-time.

III.4 The numerical results

Canceling the pole terms in $\epsilon$, the NLO FFs can be obtained by summing the finite parts from virtual and real corrections and counterterms:

[TABLE]

where the terms on the right-hand side of the equation are defined in Eqs.(19),(27),(29),(30),(72) and (III.3), and the renormalization and factorization scales are written explicitly here. The FFs $D_{\bar{b}\to B_{c}(B_{c}^{*})}^{NLO}(z,\mu_{F},\mu_{R})$ can be obtained by multiplying the matrix element $\langle{\cal O}^{B_{c}(B_{c}^{*})}(n)\rangle/\langle{\cal O}^{c\bar{b}[n]}(n)\rangle\approx|R_{S}(0)|^{2}/4\pi$ to $D_{\bar{b}\to c\bar{b}[n]}^{NLO}(z,\mu_{F},\mu_{R})$, where $n=\,^{1}S_{0}$ or ${}^{3}S_{1}$ accordingly. In the numerical calculations, the integrations over phase space are performed numerically with the help of the program Vegasvegas .

The necessary input masses in the numerical calculations are taken as follows:

[TABLE]

The value of $|R_{S}(0)|^{2}$ may be extracted from the experimental widths of the $B_{c}$ pure leptonic decays, potential-model calculations and lattice QCD calculations, etc., whereas now there is no very accurate value of $|R_{S}(0)|^{2}$. In fact, due to the fact that for the FFs $|R_{S}(0)|^{2}$ is an overall factor, the numerical results obtained in this paper with a given value of $|R_{S}(0)|^{2}$ can be easily updated with a more accurate value. Thus, as an approximation, in numerical calculations we just take the value from the potential-model calculationspot :

[TABLE]

For strong coupling constant, we adopt the two-loop formula

[TABLE]

where $\beta_{1}=\frac{34}{3}C_{A}^{2}-4C_{F}T_{F}n_{f}-\frac{20}{3}C_{A}T_{F}n_{f}$ is the two-loop coefficient of $\beta$ function in QCD. According to $\alpha_{s}(m_{{}_{Z}})=0.1185$pdg , we obtain $\Lambda_{QCD}^{n_{f}=5}=0.233{\rm GeV}$ and $\Lambda_{QCD}^{n_{f}=4}=0.337{\rm GeV}$. Then we have $\alpha_{s}(2m_{c})=0.259$, $\alpha_{s}(m_{b}+2m_{c})=0.190$, $\alpha_{s}(2m_{b})=0.180$, and $\alpha_{s}(2m_{b}+m_{c})=0.174$.

The LO FFs $D^{LO}_{\bar{b}\to B_{c}}(z,\mu_{F0},\mu_{R})$, $D^{LO}_{\bar{b}\to B^{*}_{c}}(z,\mu_{F0},\mu_{R})$ and the NLO FFs $D^{NLO}_{\bar{b}\to B_{c}}(z,\mu_{F0},\mu_{R})$, $D^{NLO}_{\bar{b}\to B^{*}_{c}}(z,\mu_{F0},\mu_{R})$(the latter is that in Eq.(III.4) are presented in Figs.4 and 5, respectively. In order to keep the logarithm terms ${\rm ln}(\mu_{R}/m_{Q})$ and ${\rm ln}(\mu_{F0}/m_{Q})$ ($m_{Q},Q=c,b$) in higher-order corrections “becoming” large and to have better accuracy, here we set the renormalization scale $\mu_{R}$ and factorization scales $\mu_{F0}$ to ${\cal O}(m_{Q})$, i.e., we set $\mu_{R}$ and $\mu_{F0}$ to be $2m_{c}$ and $m_{b}+2m_{c}$(the minimum invariant mass of the initial off-shell $\bar{b}$ quark), respectively. In Figs.4 and 5 the results for $\mu_{R}=\mu_{F0}=m_{b}+2m_{c}$ are also presented.

From Figs.4 and 5, one can see that the QCD NLO corrections to the FFs of $\bar{b}$ quark are quite large with a normalization scale $\mu_{R}=2m_{c}$ or $\mu_{R}=m_{b}+2m_{c}$. The maximum points of the FFs are shifted to smaller values of $z$ when the NLO corrections are involved. Moreover, the QCD NLO FFs are scheme and scale dependent, and the FFs in this paper are defined in the $\overline{MS}$ scheme.

There are two useful quantities which can be easily computed from the numerical results for the FFs: the fragmentation probability $P$ and the average value of $z$, $\langle z\rangle$. They are defined as follows:

[TABLE]

where $D(z)$ denotes an FF at a given energy scale. The numerical results for the obtained FFs are presented in Tables 1 and 2. From the two tables, one can see that the NLO corrections to the fragmentation probabilities are sizable with the two choices of the renormalization scale. However, due to the QCD NLO corrections the average values $\langle z\rangle$ change by only a small amount.

The FFs of a c quark to the meson $B_{c}$ or $B^{*}_{c}$ can be derived out by applying the method presented in Secs. II and III precisely. Whereas, the contributions to the FFs from the cut diagrams without a heavy-quark loop on either side of the cut can be obtained by the alternation of $m_{b}$ and $m_{c}$. The NLO QCD FFs of a $c$ quark into $B_{c}$ and $B^{*}_{c}$ mesons are presented in Figs.2 and 3 with two possible renormalization scales, $\mu_{R}=2m_{b}$ and $\mu_{R}=2m_{b}+m_{c}$, and the factorization scale is set $\mu_{F0}=2m_{b}+m_{c}$, which is the minimum invariant mass of the initial off-shell $c$ quark.

From Fig.6 and Fig.7, one can see that the NLO QCD corrections to the FFs of $B_{c}$ and $B^{*}_{c}$ are also sizable with the renormalization scales $\mu_{R}=2m_{b}$ or $\mu_{R}=2m_{b}+m_{c}$, and the difference between the FFs at the two renormalization scales is quite small. This is because these two scales are quite close to each other.

The fragmentation probabilities and average values of $z$ for $c$-quark fragmentation are presented in Tables 3 and 4. One can see that the fragmentation probability of $c\to B_{c}$($c\to B^{*}_{c}$) is smaller than that of $\bar{b}\to B_{c}$($\bar{b}\to B_{c}$) by about 2 orders of magnitude.

The FFs at a large factorization scale such as at $\mu_{F}\gg m_{Q}$ can be obtained by solving the DGLAP evolution equations from the FFs at a smaller $\mu_{F0}$ ($\sim m_{Q}$). Note that for convenience in the paper we call the FFs at a smaller factorization scale as “initial FFs”.

Here to solve the DGLAP evolution equations the approximation method introduced in Ref.APsolve is adopted, and as stated in the Introduction, the evolution of FFs are from a low energy scale to a high energy scale is restricted to the LL QCD level, namely, only the LO splitting function $P_{ij}$ ($i,j=g,q$, where $g$ is gluon and the quarks $q=b,c$) in Eq.(5) is considered.

Although for solving the DGLAP equations, the QCD NLO FFs at comparatively low energy scale of $\bar{b}$ and $c$ quarks to the mesons $B_{c}$ or $B_{c}^{*}$ provide main parts of the necessary “initial FFs”, due to the mixing of the gluon’s and flavor-singlet quarks’ FFs, the FF of a gluon at the low energy scale also is a necessary part of the initial condition, so we need to calculate out the FF of a gluon at the comparatively low energy scale $\mu_{F}=2m_{b}+2m_{c}$, where is the threshold of the $B_{c}$ or $B^{*}_{c}$ production by a gluon. Now the “initial FFs” of b̄ and c quarks as well as a gluon to the meson $B_{c}$ or $B^{*}_{c}$ all at the low energy scale $\mu_{F0}=2m_{b}+2m_{c}$, which as “initial condition” are needed for solving the DGLAP equations, are shown in Figs.8 and 9, where the FFs of $\bar{b}$ and $c$ quarks to the meson $B_{c}$ or $B^{*}_{c}$ at this energy scale are obtained by solving the DGLAP equations from $\mu_{R}=\mu_{F0}=m_{b}+2m_{c}$(for $\bar{b}\to B_{c}(B_{c}^{*})$) or $\mu_{R}=\mu_{F0}=2m_{b}+m_{c}$(for $c\to B_{c}(B_{c}^{*})$). In order to show the FF curves in one figure, in Figs.8 and 9, the gluon and $c$-quark FFs are artificially multiplied by a factor of 30. From Figs.8 and 9, one can see that the FFs for $g\to B_{c}(B^{*}_{c})$ and $c\to B_{c}(B_{c}^{*})$ are about 2 orders of magnitude smaller than the FF for $\bar{b}\to B_{c}(B_{c}^{*})$.

For the following application in the next section, with the initial FFs at $\mu_{F}=2m_{b}+2m_{c}$ which are obtained by means of this work on the QCD NLO FFs, we calculate out the FFs at the energy scale $\mu_{F}=m_{{}_{Z}}$ by solving DGLAP evolution equations Eq.(4); and the results are shown in Figs.10 and 11. For the same reason as in Figs.10 and 11, we artificially multiply the gluon and $c$ quark FFs by a factor of 30. From Figs.10 and 11, one can see that the FFs are changed due to the evolution. The average values of $z$ for the $\bar{b}$-quark and $c$-quark fragmentation are shifted to smaller values. For the $\bar{b}$-quark fragmentation,

[TABLE]

For the $c$-quark fragmentation,

[TABLE]

The fragmentation probabilities for the gluon fragmentation are increased compared to the gluon fragmentation at $\mu_{F}=2m_{b}+2m_{c}$. However, the fragmentation probabilities for $g\to B_{c}(B^{*}_{c})$ are small compared to the fragmentation probabilities for $\bar{b}\to B_{c}(B^{*}_{c})$.

IV Application to $B_{c}(B_{c}^{*})$ production at a Z factory

The production of $B_{c}$ and $B_{c}^{*}$ mesons at a Z factory is the simplest case where only the fragmentation from $\bar{b}$ and $c$ quarks should be considered, and the fragmentation from light quarks and gluon (being high-order processes) can be ignored. Moreover, this production is a typical process for doubly heavy flavored hadron production at a Z factory, which can be a good reference for doubly heavy hadron production at a Z factory. Thus, to try to have a higher accuracy for the fragmentation approach in computing the production of the $B_{c}$ and $B_{c}^{*}$ at a Z factory, we would like to apply the FFs, which are accurate up to the QCD NLO at a low factorization energy scale $\mu_{F0}$ and evolved with the DGLAP equations to the proper and higher energy scale (here it is $\mu_{F}=m_{{}_{Z}}$ ), to computing the production of the $B_{c}$ and $B_{c}^{*}$ at a Z factory, and to compare the results with those obtained by the approaches of complete LO and NLO QCD.

With the pQCD factorization (2), the differential cross sections of $B_{c}$ and $B_{c}^{*}$ production at a Z factory can be calculated straightforwardly. The expressions for the coefficient functions $d\hat{\sigma}_{e^{+}e^{-}\to i+X}/dy$ ($i=\bar{b},c$) in the limit $m_{Q}\to 0$ can be found in Refs.coefun1 ; coefun2 .

For the numerical calculations, the additional and relevant input parameters are taken as follows:

[TABLE]

where $\alpha=\alpha(m_{{}_{Z}})$ is the electromagnetic coupling constant renormalized at $m_{{}_{Z}}$.

The differential cross sections $d\sigma/dz$ for the production of $B_{c}$ and $B^{*}_{c}$ mesons at the Z pole are presented in Figs.12 and 13, where the contributions from $\gamma$ exchange and $\gamma-Z$ interference are neglected due to the fact that they are much smaller in comparison to the contributions from Z exchange at the Z polebcnlo . In Figs.12 and 13, “LO” and “NLO” denote the results of the complete LO and NLO calculations, respectively, “Frag” denotes the results of the fragmentation approach and the leading logarithms (LLs) being resummed through DGLAP evolution equations. For the results of the complete LO and NLO approaches, the renormalization is set at $\mu_{R}=m_{b}+2m_{c}$666In order to maintain compatibility with the results of the fragmentation calculations, for the complete NLO calculations we adopt the renormalization scale $\mu_{R}=m_{b}+2m_{c}$, although the complete NLO results in our previous paperbcnlo are those with the renormalization at $\mu_{R}=2m_{b}$..

It is interesting to compare the total cross sections for the production of $B_{c}$ and $B_{c}^{*}$ mesons obtained using the fragmentation approach with those from the complete LO and NLO calculations. The obtained total cross sections are presented in Table 5. We believe that the fragmentation approach provides better results for the values of the total cross sections of $B_{c}$ and $B_{c}^{*}$ production at the Z pole.

V Discussions and Conclusion

In this paper, by means of the general operator definition of the FFs, we have derived the FFs for a $\bar{b}$ or $c$ quark fragmenting to $B_{c}$ and $B^{*}_{c}$ mesons in LO and NLO QCD, and the numerical results with reasonable input parameters are presented in figures.

In the derivation of the NLO “real corrections” to the FFs, the difficulty in extracting the singularities is overcome by the fact that certain proper subtraction terms are constructed, which contain the exact same singularities ($1/\epsilon$) as those in the real corrections under dimensional regularization, but they can be computed almost analytically [see Eqs.(30), (31), and (33)]. Then, with the constructed auxiliary terms for subtractions, the singular and finite contributions from the real corrections can be computed separately and the finite contributions can be calculated numerically. Note that here the integrations of the subtraction terms over the phase space are carried out under suitable parametrizations, which are very similar to those introduced in Ref.braaten , and the expressions for the subtraction terms and the phase-space parametrizations may be useful in calculating the real QCD NLO corrections for other FFs.

It is known that the choices for the factorization scale $\mu_{F}$ and the renormalization scale $\mu_{R}$ are very important in QCD calculations. For NLO corrections of FFs, one may set them equal to each other or different from each other according to convenience. As a typical case, here we set the “(initial) factorization scale” to $\mu_{F0}=m_{b}+2m_{c}$ for the QCD NLO “initial FFs”, and the results show that the NLO corrections are significant with two possible choices of the renormalization scale. Moreover, for an important application specifically discussed in this paper, to gain a higher accuracy FFs with the factorization $\mu_{F}=m_{{}_{Z}}$ were used. We obtained FFs by solving the DGLAP evolution equation, starting with the “initial QCD NLO FFs” at a low energy scale $\mu_{F}=2m_{b}+2m_{c}$. Since the solution of the DGLAP evolution equation shows certain shifts of the average value of the energy fraction $z$ in a small region, we hope that future experiments can test this effect(s).

Finally, the production of $B_{c}$ and $B_{c}^{*}$ mesons at a Z factory is the simplest case where only the fragmentation from a $\bar{b}$ or $c$ quark should be considered, and this production is a typical process for doubly heavy hadron production at a Z factory, which can be used as a reference to estimate doubly heavy hadron production at a Z factory. Thus we applied the FFs at the energy scale $\mu_{F}=m_{{}_{Z}}$, which were obtained by evolving the QCD NLO ones at a low energy scale $\mu_{F0}$ and shown in Figs.10 and 11, to computing the production of $B_{c}$ and $B_{c}^{*}$ mesons at a Z factory, and we suspect that the results presented in Figs.12 are 13 are comparatively accurate. For comparison, the results from the fully complete LO and NLO calculations were also presented in these figures.

In summary, we derived the QCD NLO FFs of $\bar{b}$ and $c$ quarks to $B_{c}$ and $B_{c}^{*}$ mesons, and the physical picture for the production of $B_{c}$ and $B_{c}^{*}$ mesons to QCD leading logarithm (LL) order at a Z factory was described as follows: the $\bar{b}$ and $c$ quarks are produced at high energy ($\sqrt{s}=m_{{}_{Z}}$), then the produced $\bar{b}$ and $c$ quarks are evolved to the lower invariant mass (${\cal O}(m_{Q})$) by emitting real and virtual collinear gluons and quarks (that are summed by the LO DGLP equations), at last they fragment into the meson $B_{c}$ or $B_{c}^{*}$, that is described by the QCD NLO FFs. Therefore one may reasonably understand why the physics picture summarized here has more solid QCD foundation and works better in estimating the $B_{c}$ and $B_{c}^{*}$ production at a Z factory.

Acknowledgments: This work was supported in part by the National Natural Science Foundation of China (NSFC) under Grants No. 11745006, No. 11535002, No. 11675239, No. 11821505, No. 11275036, No. 11625520, No. 11705045, and No. 11847222. It was also supported in part by the Key Research Program of Frontier Sciences, CAS, Grant No. QYZDY-SSW-SYS006.

Appendix A The expressions for the coefficients $a_{i}$ in the LO squared amplitudes

For the production of the $c\bar{b}[^{1}S_{0}^{[1]}]$ state, the expressions for the coefficients $a_{i}$ in Eq.(16) are

[TABLE]

For the production of the $c\bar{b}[^{3}S_{1}^{[1]}]$ state,

[TABLE]

Appendix B Phase space for the real corrections

The differential phase space for the real corrections to the FFs is

[TABLE]

The different parameterizations are required in order to extract the poles in $\epsilon$ in the real corrections. We adopt similar parametrizations as those used in Ref.braaten . In Ref.braaten , the authors derived the phase space for two massless partons in final state. In our case, there is one massive parton and one massless parton in the final state, so we derive the formulas for this case.

The differential phase space for a single parton with momentum $p$ and mass $m$ can be expressed as

[TABLE]

where $\theta$ denotes the polar angle and $d\Omega_{\perp}$ denotes the differential transverse solid angle. The total transverse solid angle $\Omega_{\perp}=\int d\Omega_{\perp}=2\pi^{1-\epsilon}/\Gamma(1-\epsilon)$.

It is useful to introduce a light-like momentum $k$, and define the variable

[TABLE]

Then, the differential phase space for a single parton can be expressed as

[TABLE]

Here, $d\Omega_{\perp}$ is Lorentz invariant due to the fact that the differential phase space, $\lambda$ and $p\cdot n$ are Lorentz invariant. This expression can be easily derived from Eq.(85) in a Lorentz frame where the spatial parts of the light-like vectors $n$ and $k$ are back to back. The differential phase space for a massless parton can be easily obtained from Eq.(87) by setting $m=0$.

We can apply the parametrization Eq.(87) to the differential phase spaces for $p_{2}$ and $p_{3}$ in Eq.(84). Two light-like vectors $k_{2}$ and $k_{3}$ corresponding to the parametrizations of $p_{2}$ and $p_{3}$ are introduced. The integral over $p_{2}\cdot n$ can be carried out through the $\delta$ function, and the integral over $\Omega_{2\perp}$ is trivial. Then we obtain the expression

[TABLE]

where

[TABLE]

We have converted the integral variable $p_{3}\cdot n$ to $u$ by using the definition of $u$ in Eq.(32). If we set $m_{c}=0$ in Eq.(LABEL:eqb5), we obtain an expression that is the same as Eq.(A.6) in Ref.braaten .

We need to choose proper light-like vectors $k_{2}$ and $k_{3}$ in order to extract the poles in $\epsilon$. For the subtraction terms that contain $s$, we choose

[TABLE]

then

[TABLE]

and

[TABLE]

Changing variables in Eq.(LABEL:eqb5) from $u$, $\lambda_{2}$ and $\lambda_{3}$ to $y$, $s_{1}$ and $s$, we obtain

[TABLE]

For the subtraction terms that contain $s_{2}$, we choose

[TABLE]

Then

[TABLE]

After changing variables in Eq.(LABEL:eqb5) from $u$, $\lambda_{2}$ and $\lambda_{3}$ to $y$, $s_{1}$ and $s_{2}$, we obtain

[TABLE]

For the subtraction terms that contain $s_{3}$, we choose

[TABLE]

Then,

[TABLE]

After changing variables in Eq.(LABEL:eqb5) from $u$, $\lambda_{2}$ and $\lambda_{3}$ to $y$, $s_{1}$ and $s_{3}$, we obtain

[TABLE]

To derive the differential phase space for the subtraction terms that contain $t_{1}$ or $t_{2}$, we multiply Eq.(84) by

[TABLE]

which is equal to 1 and does not change the phase space. After integrating over $p_{2}$, the differential phase space can be expressed as

[TABLE]

Using the parametrization Eq.(87) on the differential phase spaces for $\tilde{p}$ and $p_{3}$ in Eq.(101), we obtain the expression

[TABLE]

where

[TABLE]

For the subtraction terms that contain $t_{1}$, we choose the light-like vectors $\tilde{k}$ and $k_{3}$ as follows

[TABLE]

then

[TABLE]

After changing variables in Eq.(LABEL:eqb19) from $\tilde{\lambda}$ and $\lambda_{3}$ to $\tilde{s}$ and $t_{1}$, we obtain

[TABLE]

For the subtraction terms that contain $t_{2}$, we choose the light-like vectors $\tilde{k}$ and $k_{3}$ as follows:

[TABLE]

then

[TABLE]

[TABLE]

After changing variables from $\tilde{\lambda}$ and $\lambda_{3}$ to $\tilde{s}$ and $t$, we obtain

[TABLE]