# Fragmentation functions for a quark into a spin-singlet quarkonium:
Different flavor case

Xu-Chang Zheng<EMAIL_ADDRESS>Ze-Yang Zhang<EMAIL_ADDRESS>Xing-
Gang Wu<EMAIL_ADDRESS>Department of Physics, Chongqing University,
Chongqing 401331, P.R. China.

###### Abstract

In the paper, we calculate the fragmentation functions for a quark to fragment
into a spin-singlet quarkonium, where the flavor of the initial quark is
different from that of the constituent quark in the quarkonium. The
ultraviolet divergences in the phase space integral are removed through the
operator renormalization under the modified minimal subtraction scheme. The
fragmentation function $D_{q\to\eta_{Q}}(z,\mu_{F})$ is expressed as a two-
dimensional integral. Numerical results for the fragmentation functions of a
light quark or a bottom quark to fragment into the $\eta_{c}$ are presented.
As an application of those fragmentation functions, we study the processes
$Z\to\eta_{c}+q\bar{q}g(q=u,d,s)$ and $Z\to\eta_{c}+b\bar{b}g$ under the
fragmentation and the direct nonrelativistic QCD approaches.

##  I Introduction

According to QCD factorization theorem, the cross section for the inclusive
production of a hadron $H$ with high transverse momentum ($p_{T}$) in a high-
energy collision is dominated by the single parton fragmentation
Collins:1989gx , i.e.,

$\displaystyle d\sigma_{A+B\to H+X}(p_{T})=$
$\displaystyle\sum_{i}d\hat{\sigma}_{A+B\to i+X}(p_{T}/z,\mu_{F})$ (1)
$\displaystyle\otimes D_{i\to H}(z,\mu_{F})+{\cal O}(m_{H}^{2}/p_{T}^{2}),$

where $\otimes$ denotes a convolution in the momentum fraction $z$, the sum
extends over all species of partons. $d\hat{\sigma}_{A+B\to i+X}$ indicates
the partonic cross section that can be calculated in perturbation theory,
while $D_{i\to H}$ indicates the fragmentation function for the parton $i$
into a hadron $H$. $\mu_{F}$ denotes the factorization scale which is
introduced to separate the energy scales of the two parts.

The factorization formula (1) was first derived by Collins and Soper for light
hadron production Collins:1981uw . This factorization formula can be equally
applied to the heavy quarkonium production. The proof of the factorization
formula (1) for the quarkonium production was presented by Nayak, Qiu, and
Sterman Nayak:2005rt . The factorization formula (1) is called leading power
(LP) factorization because it gives the LP contribution in the expansion in
powers of $m_{H}/p_{T}$. The factorization formula for the next-to-leading
power (NLP) correction was derived in Refs.Kang:2011zza ; Kang:2011mg ;
Fleming:2012wy ; Fleming:2013qu , and the NLP contribution comes from the
double-parton fragmentation.

Fragmentation functions play an important role in the calculation of the cross
sections under the LP factorization. Unlike the fragmentation functions for
the production of the light hadrons which are nonperturbative in nature, the
fragmentation functions for the heavy quarkonium production can be calculated
through the nonrelativistic QCD (NRQCD) factorization nrqcd . Under NRQCD
factorization, the fragmentation functions for a parton to fragment into a
quarkonium can be written as

$\displaystyle D_{i\to
H}(z,\mu_{F})=\sum_{n}d_{i\to(Q\bar{Q})[n]}(z,\mu_{F})\langle{\cal
O}^{H}(n)\rangle,$ (2)

where $d_{i\to(Q\bar{Q})[n]}$ are short-distance coefficients (SDCs) which can
be expanded as powers of $\alpha_{s}(m_{Q})$, and $\langle{\cal
O}^{H}(n)\rangle$ are long-distance matrix elements (LDMEs).

The fragmentation functions for quarkonia have been studied extensively. Most
of the fragmentation functions for the S-wave and P-wave quarkonia are known
up to $\alpha_{s}^{2}$ order Chang:1992bb ; Braaten:1993jn ; Braaten:1993mp ;
Braaten:1993rw ; Braaten:1994kd ; Braaten:1995cj ; Chen:1993ii ; Yuan:1994hn ;
Ma:1994zt ; Ma:1995ci ; Ma:1995vi ; Cho:1994gb ; Beneke:1995yb ;
Braaten:2000pc ; Hao:2009fa ; Sang:2009zz ; Jia:2012qx ; Bodwin:2014bia ;
Ma:2013yla ; Ma:2015yka ; Yang:2019gga , and a few fragmentation functions for
quarkonia were calculated up to $\alpha_{s}^{3}$ order Artoisenet:2014lpa ;
Sepahvand:2017gup ; Artoisenet:2018dbs ; Feng:2018ulg ; Zhang:2018mlo ;
Zheng:2019dfk ; Zheng:2019gnb ; Feng:2017cjk . Among these studies, the next-
to-leading order (NLO) corrections to the fragmentation functions for $g\to
Q\bar{Q}[^{1}S_{0}^{[1,8]}]$ have been calculated recently by three groups
Artoisenet:2014lpa ; Artoisenet:2018dbs ; Feng:2018ulg ; Zhang:2018mlo . The
NLO fragmentation functions for $g\to Q\bar{Q}[^{1}S_{0}^{[1,8]}]$ are
important in prediction of the production of the $\eta_{c,b}$ and the
$h_{c,b}$ at the LHC. However, the fragmentation functions for a quark into
the $\eta_{c,b}$ are only available up to $\alpha_{s}^{2}$ order. Those
fragmentation functions are also important to the precision prediction of the
$\eta_{c,b}$ production at the LHC. Moreover, for the production of the
$\eta_{c,b}$ in $e^{+}e^{-}$ collisions, the quark fragmentation contribution
is more important than the gluon fragmentation contribution due to the fact
that the cross section for a quark is $\alpha_{s}^{0}$ order but the cross
section for a gluon is $\alpha_{s}$ order. In this paper, we will calculate
the fragmentation functions for a quark $q$ into the $\eta_{Q}$, where $Q=c,b$
but $q\neq Q$.

In the early calculations of fragmentation functions for doubly heavy mesons
Chang:1992bb ; Chen:1993ii , the fragmentation functions are determined
through comparing the cross section calculated based on the NRQCD
factorization with that calculated based on the factorization formula (1) for
a process containing the doubly heavy meson being produced. In fact, the
fragmentation functions can be defined through the matrix elements of nonlocal
gauge-invariant operators Collins:1981uw . The operator definition for the
fragmentation functions was first applied to calculations of the fragmentation
functions for doubly heavy mesons by Ma Ma:1994zt . The calculations based on
the operator definition are particularly convenient to extend to higher
orders. Therefore, we will calculate the fragmentation functions based on the
operator definition suggested by Collins and Soper.

The paper is organized as follows. Following the Introduction, in Sec.II, we
present the definition and the analytical calculation for the fragmentation
functions. In Sec.III, we present the numerical results for the fragmentation
functions $D_{q\to\eta_{c}}(z,\mu_{F})(q=u,d,s)$ and
$D_{b\to\eta_{c}}(z,\mu_{F})$, and apply the fragmentation functions to
processes $Z\to\eta_{c}+q\bar{q}g(q=u,d,s)$ and $Z\to\eta_{c}+b\bar{b}g$.
Section IV is reserved for a summary.

##  II The analytical calculation for the fragmentation functions

###  II.1 The definition of fragmentation function

The fragmentation functions are usually defined in the light-cone coordinate
system. In this coordinate system, a d-dimensional vector $V$ is expressed as
$V^{\mu}=(V^{+},V^{-},{\bf V}_{\perp})$, with $V^{+}=(V^{0}+V^{d-1})/\sqrt{2}$
and $V^{-}=(V^{0}-V^{d-1})/\sqrt{2}$. Then the product of two vectors is
$V\cdot W=V^{+}W^{-}+V^{-}W^{+}-{\bf V}_{\perp}\cdot{\bf W}_{\perp}$. The
gauge-invariant fragmentation function for a quark $q$ to fragment into a
hadron $H$ is defined as Collins:1981uw

$\displaystyle D_{q\to H}(z)=$ $\displaystyle\frac{z^{d-3}}{2\pi}\sum_{X}\int
dx^{-}e^{-iP^{+}x^{-}/z}$ (3) $\displaystyle\times\frac{1}{N_{c}}{\rm Tr}_{\rm
color}\frac{1}{4}{\rm Tr}_{\rm Dirac}\left\\{\gamma^{+}\langle
0|\Psi(0)\bar{{\cal P}}{\rm
exp}\left[ig_{s}\int_{0}^{\infty}dy^{-}A_{a}^{+}(0^{+},y^{-},{\bf
0}_{\perp})t_{a}^{T}\right]|H(P^{+},{\bf 0}_{\perp})+X\rangle\right.$
$\displaystyle\left.\times\langle H(P^{+},{\bf 0}_{\perp})+X|{\cal P}{\rm
exp}\left[-ig_{s}\int_{x^{-}}^{\infty}dy^{-}A_{a}^{+}(0^{+},y^{-},{\bf
0}_{\perp})t_{a}^{T}\right]\bar{\Psi}(x)|0\rangle\right\\},$

where $\Psi$ is the field of initial quark, $A_{a}^{\mu}$ is the gluon field,
and $t_{a}(a=1\cdots 8)$ are $SU(3)$-color matrices. The longitudinal momentum
fraction is defined as $z\equiv P^{+}/K^{+}$, where $K$ is the momentum of the
initial quark. The fragmentation function is defined in a reference frame in
which the transverse momentum of the hadron $H$ vanishes. It is convenient to
introduce a lightlike momentum whose expression is $n^{\mu}=(0,1,{\bf
0}_{\perp})$ in the reference frame where the definition of the fragmentation
function carried out. Then, $z$ can be expressed as a Lorentz invariant, i.e.,
$z=P\cdot n/K\cdot n$. The Feynman rules can be derived from the definition
(3) directly, and we have presented the Feynman rules in a previous paper
Zheng:2019gnb .

###  II.2 The calculation of fragmentation function

The definition (3) is gauge invariant. However, for the practical calculation,
the gauge should be specified. We adopt the usual Feynman gauge throughout the
paper. There are ultraviolet (UV) divergences in the calculation. To deal with
the UV divergences, we adopt dimensional regularization with $d=4-2\epsilon$,
then the UV divergences appear as the pole terms in $\epsilon$.

In the calculation, we first calculate the fragmentation function for an on-
shell $Q\bar{Q}$ pair in ${}^{1}S_{0}^{[1]}$ state. Then the fragmentation
function $D_{q\to\eta_{Q}}$ can be obtained from
$D_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}$ through replacing the LDME $\langle{\cal
O}^{(Q\bar{Q})[^{1}S_{0}^{[1]}]}(^{1}S_{0}^{[1]})\rangle$ by $\langle{\cal
O}^{\eta_{Q}}(^{1}S_{0}^{[1]})\rangle$.

Figure 1: The cut diagrams for the fragmentation function
$D_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}$, where $q\neq Q$.

There are 16 cut diagrams for
$q(K)\to(Q\bar{Q})[^{1}S_{0}^{[1]}](p_{1})+g(p_{2})+q(p_{3})$ under the
Feynman gauge, which can be collectively represented by four diagrams in
Fig.1. The squared amplitudes, corresponding to four diagrams in Fig.1, can be
written as

$\displaystyle{\cal A}_{1}=$ $\displaystyle{\rm
tr}\left[\not{n}\frac{-i}{\not{p_{1}}+\not{p_{2}}+\not{p_{3}}-m_{q}-i\epsilon}(ig_{s}\gamma_{\nu}t^{b})(\not{p}_{3}+m_{q})\right.$
(4)
$\displaystyle\left.(-ig_{s}\gamma_{\mu}t^{a})\frac{i}{\not{p_{1}}+\not{p_{2}}+\not{p_{3}}-m_{q}+i\epsilon}\right]X_{ab}^{\mu\nu},$
$\displaystyle{\cal A}_{2}=$ $\displaystyle{\rm
tr}\left[\not{n}\frac{-i}{\not{p_{1}}+\not{p_{2}}+\not{p_{3}}-m_{q}-i\epsilon}(ig_{s}\gamma_{\nu}t^{b})(\not{p}_{3}+m_{q})\right.$
(5) $\displaystyle\left.\frac{i}{(p_{1}+p_{2})\cdot
n+i\epsilon}(ig_{s}n_{\mu}t^{a})\right]X_{ab}^{\mu\nu},$ $\displaystyle{\cal
A}_{3}=$ $\displaystyle{\rm
tr}\left[\not{n}(-ig_{s}n_{\nu}t^{b})\frac{-i}{(p_{1}+p_{2})\cdot
n-i\epsilon}(\not{p}_{3}+m_{q})\right.$ (6)
$\displaystyle\left.(-ig_{s}\gamma_{\mu}t^{a})\frac{i}{\not{p_{1}}+\not{p_{2}}+\not{p_{3}}-m_{q}+i\epsilon}\right]X_{ab}^{\mu\nu},$
$\displaystyle{\cal A}_{4}=$ $\displaystyle{\rm
tr}\left[\not{n}(-ig_{s}n_{\nu}t^{b})\frac{-i}{(p_{1}+p_{2})\cdot
n-i\epsilon}(\not{p}_{3}+m_{q})\right.$ (7)
$\displaystyle\left.\frac{i}{(p_{1}+p_{2})\cdot
n+i\epsilon}(ig_{s}n_{\mu}t^{a})\right]X_{ab}^{\mu\nu},$

where $\Pi_{1}$ is the spin-singlet projector

$\displaystyle\Pi_{1}=\frac{1}{(2\,m_{Q})^{3/2}}(\not{p}_{1}/2-m_{Q})\gamma_{5}(\not{p}_{1}/2+m_{Q}),$
(8)

$\Lambda_{1}$ is the color-singlet projector

$\displaystyle\Lambda_{1}=\frac{\textbf{1}}{\sqrt{3}},$ (9)

where 1 is the unit matrix of the $SU(3)_{c}$ group. There is a common factor
$X_{ab}^{\mu\nu}$ which arises from the annihilation of a virtual gluon into a
$(Q\bar{Q})[^{1}S_{0}^{[1]}]$ pair and a real gluon. The factor
$X_{ab}^{\mu\nu}$ can be expressed as

$\displaystyle
X_{ab}^{\mu\nu}=-g_{\rho\sigma}J^{\mu\rho}_{ac}J^{\nu\sigma*}_{bc},$ (10)

where

$\displaystyle J^{\mu\rho}_{ac}=$ $\displaystyle{\rm
tr}\left\\{\Pi_{1}\Lambda_{1}\left[(-ig_{s}\gamma_{\rho}t^{c})\frac{i}{\not{p}_{1}/2+\not{p}_{2}-m_{Q}}(-ig_{s}\gamma_{\mu}t^{a})\right.\right.$
(11)
$\displaystyle\left.\left.+(-ig_{s}\gamma_{\mu}t^{a})\frac{i}{-\not{p}_{1}/2-\not{p}_{2}-m_{Q}}(-ig_{s}\gamma_{\rho}t^{c})\right]\right\\}$
$\displaystyle\cdot\frac{-i}{(p_{1}+p_{2})^{2}+i\epsilon}.$

We employ the package FeynCalc Mertig:1990an ; Shtabovenko:2016sxi  to carry
out the Dirac and color traces, and then the total squared amplitude (${\cal
A}_{(1-4)}\equiv\sum_{i=1}^{4}{\cal A}_{i}$) can be written as

$\displaystyle{\cal A}_{(1-4)}=$
$\displaystyle\frac{c_{1}(s_{1},y,z)\,(p_{3}\cdot\tilde{p})^{2}}{s_{1}^{2}(s-m_{q}^{2})^{2}}+\frac{c_{2}(s_{1},y,z)\,p_{1}\cdot
p_{3}}{s_{1}(s-m_{q}^{2})^{2}}$ (12)
$\displaystyle+\frac{c_{3}(s_{1},y,z)\,p_{2}\cdot
p_{3}}{s_{1}(s-m_{q}^{2})^{2}}+\frac{c_{4}(s_{1},y,z)}{(s-m_{q}^{2})^{2}},$

where

$\displaystyle s_{1}=(p_{1}+p_{2})^{2},~{}~{}s=(p_{1}+p_{2}+p_{3})^{2},$
$\displaystyle y=\frac{(p_{1}+p_{2})\cdot n}{K\cdot
n},~{}\tilde{p}=p_{1}-\frac{z\,p_{2}}{(y-z)}.$ (13)

The coefficients $c_{i}(s_{1},y,z)$ can be easily extracted, and we do not
list their expressions here.

The differential phase space for the fragmentation function
$D_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}$ is111Here, we associate the scale factor
$\mu^{4-d}$ with each dimensionally regulated integration in $d$ space-time
dimensions. In our previous papers Zheng:2019dfk ; Zheng:2019gnb , this scale
factor was put in the squared amplitudes.

$\displaystyle d\phi_{3}(p_{1},p_{2},p_{3})=$ $\displaystyle
2\pi\delta\left(K^{+}-\sum_{i=1}^{3}p_{i}^{+}\right)\mu^{2(4-d)}$ (14)
$\displaystyle\times\prod_{i=2,3}\frac{\theta(p_{i}^{+})dp_{i}^{+}}{4\pi
p_{i}^{+}}\frac{d^{d-2}\textbf{p}_{i\perp}}{(2\pi)^{d-2}}.$

The contributions from the cut diagrams shown in Fig.1 can be calculated
through

$\displaystyle D^{(1-4)}_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}(z)=N_{CS}\int
d\phi_{3}(p_{1},p_{2},p_{3}){\cal A}_{(1-4)},$ (15)

where $N_{CS}\equiv z^{d-3}/(8\pi N_{c})$ is a factor from the definition of
the fragmentation function. The integral on the right-hand side of Eq.(15) is
UV divergent with $d=4$. This UV divergence is regularized by dimensional
regularization with $d=4-2\epsilon$, and the integral generates $1/\epsilon$
terms. To perform the integration in Eq.(15), it is important to choose proper
parametrization for the phase space. We present a parametrization for the
phase space in Appendix A.

The differential phase space given in Eq.(39) can be expressed as follows

$\displaystyle N_{CS}d\phi_{3}(p_{1},p_{2},p_{3})$
$\displaystyle=N_{g}(p_{1},p_{2})d\phi_{2}(p_{1},p_{2})d\phi^{(3)}(p_{1},p_{2},p_{3}),$
(16)

where $N_{g}(p_{1},p_{2})$ is defined as

$\displaystyle
N_{g}(p_{1},p_{2})=\frac{(z/y)^{1-2\epsilon}}{(N_{c}^{2}-1)(2-2\epsilon)(2\pi\,y\,K\cdot
n)},$ (17)

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

$\displaystyle d\phi_{2}(p_{1},p_{2})=$
$\displaystyle\frac{z^{-1+\epsilon}(y-z)^{-\epsilon}\mu^{2\epsilon}}{2(4\pi)^{1-\epsilon}\Gamma(1-\epsilon)K\cdot
n}$ (18)
$\displaystyle\times\left(s_{1}-\frac{y}{z}4m_{Q}^{2}\right)^{-\epsilon}ds_{1},$

where the range of $s_{1}$ is from $(4m_{Q}^{2}y/z)$ to $\infty$.
$d\phi_{2}(p_{1},p_{2})$ stands for the differential phase space for a gluon
with longitudinal momentum $yK\cdot n$ to fragment into a
$(Q\bar{Q})[^{1}S_{0}^{[1]}]$-pair with longitudinal momentum $zK\cdot n$ at
leading order (LO). According to Eqs.(16), (17), (18) and (39), the expression
of $d\phi^{(3)}(p_{1},p_{2},p_{3})$ can be derived

$\displaystyle d\phi^{(3)}(p_{1},p_{2},p_{3})$
$\displaystyle=\frac{(N_{c}^{2}-1)(2-2\epsilon)\mu^{2\epsilon}K\cdot
n}{16N_{c}(2\pi)^{3-2\epsilon}}y^{1-\epsilon}(1-y)^{-\epsilon}$
$\displaystyle~{}~{}~{}\times\left[s-s_{1}/y-m_{q}^{2}/(1-y)\right]^{-\epsilon}ds\,dy\,d\Omega_{3\perp}.$
(19)

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

The integrations over $\Omega_{3\perp}$ and $s$ of ${\cal A}_{(1-4)}$ can be
performed using the method introduced in Ref.Zheng:2019gnb . Then, we obtain

$\displaystyle D^{(1-4)}_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}(z)=$
$\displaystyle\frac{(4\pi\mu^{2})^{\epsilon}\Gamma(\epsilon)}{(4\pi)^{2}}\int_{z}^{1}dy$
(20) $\displaystyle\int N_{g}d\phi_{2}(p_{1},p_{2})f(s_{1},y,z),$

where $N_{g}d\phi_{2}(p_{1},p_{2})\equiv
N_{g}(p_{1},p_{2})d\phi_{2}(p_{1},p_{2})$. The expression of $f(s_{1},y,z)$ is
given in Appendix B.

The contribution $D^{(1-4)}_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}(z)$ contains a
UV pole, it should be removed through the operator renormalization
Mueller:1978xu . We carry out the renormalization using the $\overline{\rm
MS}$ procedure. Then the fragmentation function under the $\overline{\rm MS}$
scheme can be obtained through

$\displaystyle D_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}(z,\mu_{F})$
$\displaystyle=D^{(1-4)}_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}(z)-\frac{\alpha_{s}}{2\pi}\left[\frac{1}{\epsilon_{UV}}-\gamma_{E}+{\rm
ln}~{}(4\pi)+{\rm ln}\frac{\mu^{2}}{\mu_{F}^{2}}\right]$
$\displaystyle~{}~{}~{}\times\int_{z}^{1}\frac{dy}{y}P_{gq}(y)D_{g\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}^{\rm
LO}(z/y),$ (21)

where $\mu_{F}$ is the factorization scale, the expression of the splitting
function $P_{gq}(y)$ is

$\displaystyle P_{gq}(y)=C_{F}\frac{1+(1-y)^{2}}{y},$ (22)

and $D_{g\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}^{\rm LO}$ is the LO fragmentation
function in $d$-dimensional space-time. In the calculation, it is convenient
to use the unintegrated form of $D_{g\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}^{\rm
LO}$, i.e.,

$\displaystyle D_{g\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}^{\rm LO}(z/y)=\int
N_{g}d\phi_{2}(p_{1},p_{2}){\cal A}_{g\to(Q\bar{Q})[^{1}S_{0}^{[1]}]},$ (23)

where the expression of ${\cal A}_{g\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}$ is

$\displaystyle{\cal A}_{g\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}$
$\displaystyle=\frac{16g_{s}^{4}}{3}(1-2\epsilon)(y\,K\cdot
n)^{2}\left[\frac{(1-z/y)}{m_{Q}^{3}s_{1}}-\frac{(1-z/y)}{m_{Q}^{3}(s_{1}-4m_{Q}^{2})}\right.$
$\displaystyle~{}~{}~{}\left.+\frac{2(1-\epsilon)}{m_{Q}s_{1}^{2}}+\frac{4(1-z/y)^{2}}{m_{Q}(s_{1}-4m_{Q}^{2})^{2}}\right].$
(24)

The integral over $s_{1}$ in Eq.(23) can be carried out easily, and we have
checked that our result for $D_{g\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}^{\rm LO}$ is
consistent with that obtained in Refs.Braaten:1993rw ; Artoisenet:2014lpa ;
Feng:2018ulg ; Zhang:2018mlo .

Applying Eqs.(20) and (23) to Eq.(21), we can obtain the fragmentation
function under the $\overline{\rm MS}$ scheme. It is found that the UV pole of
$D^{(1-4)}_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}(z)$ is exactly canceled by the UV
pole of the counter term from the operator renormalization. The remaining
integrals no longer generate divergence, we can set $\epsilon=0$ before
carrying out the integrations. Multiplying the fragmentation function
$D_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}(z,\mu_{F})$ for the
$(Q\bar{Q})[^{1}S_{0}^{[1]}]$ pair by a factor $\langle{\cal
O}^{\eta_{Q}}(^{1}S_{0}^{[1]})\rangle/\langle{\cal
O}^{(Q\bar{Q})[^{1}S_{0}^{[1]}]}(^{1}S_{0}^{[1]})\rangle\approx|R_{S}^{(Q\bar{Q})}(0)|^{2}/(4\pi)$,
we obtain the fragmentation function $D_{q\to\eta_{Q}}(z,\mu_{F})$, i.e.,

$\displaystyle D_{q\to\eta_{Q}}(z,\mu_{F})$
$\displaystyle=\int_{z}^{1}dy\int_{4m_{Q}^{2}y/z}^{\infty}ds_{1}\,g(s_{1},\mu_{F},y,z),$
(25)

where

$\displaystyle g(s_{1},\mu_{F},y,z)=$
$\displaystyle\frac{\alpha_{s}^{3}\,|R_{S}^{(Q\bar{Q})}(0)|^{2}}{9\,\pi^{2}\,y^{4}m_{Q}s_{1}^{2}(s_{1}-4m_{Q}^{2})^{2}\left[(1-y)s_{1}+y^{2}m_{q}^{2}\right]}\bigg{\\{}(y-1)\Big{[}s_{1}^{3}\big{(}y^{4}-2y^{3}(z+1)+2y^{2}(z^{2}+6z+1)$
(26)
$\displaystyle-12y\,z(z+1)+12z^{2}\big{)}+s_{1}^{2}\big{(}2y^{4}m_{q}^{2}-4y^{3}(2m_{Q}^{2}(z+4)+z\,m_{q}^{2})+4y^{2}(4m_{Q}^{2}(3z+2)+z^{2}m_{q}^{2})$
$\displaystyle-48y\,z\,m_{Q}^{2}\big{)}+16s_{1}y^{2}m_{Q}^{2}\big{(}y^{2}m_{Q}^{2}-y(2m_{Q}^{2}+z\,m_{q}^{2})+2m_{Q}^{2}\big{)}+32y^{4}m_{Q}^{4}m_{q}^{2}\Big{]}-\left[(1-y)s_{1}+y^{2}m_{q}^{2}\right]$
$\displaystyle\left.\times\left[z^{2}(s_{1}-4m_{Q}^{2}y/z)^{2}+(y-z)^{2}s_{1}^{2}\right]\left[\big{(}(1-y)^{2}+1\big{)}{\rm
ln}\left(\frac{(1-y)s_{1}+y^{2}m_{q}^{2}}{\mu_{F}^{2}}\right)+y^{2}\right]\right\\}.$

Here, $R_{S}^{(Q\bar{Q})}(0)$ is the radial wave function at the origin for
the $(Q\bar{Q})$ bound state.

##  III Numerical results and discussion

In this section, we will present the numerical results for the fragmentation
functions and apply the fragmentation functions to the decay widths for the
$\eta_{c}$ production through Z boson decays.

The input parameters for the numerical calculation are taken as follows

$\displaystyle m_{c}=1.5\,{\rm GeV},\,m_{b}=4.9\,{\rm
GeV},\,m_{{}_{Z}}=91.1876\,{\rm GeV},$ $\displaystyle\alpha=1/128,{\rm
sin}^{2}\theta_{{}_{W}}=0.231,|R_{S}^{c\bar{c}}(0)|^{2}=0.810\,{\rm GeV}^{3}.$
(27)

The value of $|R_{S}^{c\bar{c}}(0)|^{2}$ is taken from the potential model
calculation pot . For the strong coupling constant, we adopt two-loop formula
as used in our previous paper Zheng:2019dfk , where
$\alpha_{s}(2m_{c})=0.259$.

###  III.1 The fragmentation functions

The fragmentation function for a light quark into $\eta_{c}$, where
$\mu_{F}=2m_{c}$, $4m_{c}$, and $6m_{c}$, is shown in Fig.2. In the numerical
calculation, the mass of the light quark is neglected, and the strong coupling
is taken as $\alpha_{s}(2m_{c})$.

Figure 2: The fragmentation function $D_{q\to\eta_{c}}(z,\mu_{F})$ as a
function of $z$ for $\mu_{F}=2m_{c}$, $4m_{c}$, and $6m_{c}$, where $q$
denotes a light quark.

From Fig.2, we can see that the fragmentation function is sensitive to the
factorization scale. When $\mu_{F}=2m_{c}$, the fragmentation function
increases first ($z<0.96$) and then decreases ($z>0.96$) with the increase of
$z$, and the fragmentation function is less than 0 in the most $z$ region;
when $\mu_{F}=4m_{c}$, the fragmentation function also increases first
($z<0.12$) and then decreases ($z>0.12$) with the increase of $z$, but the
fragmentation function is greater than 0 in the most $z$ region; When
$\mu_{F}=6m_{c}$, the fragmentation function decreases monotonically with the
increase of $z$, and the fragmentation function is greater than 0 for
$z\in(0,1)$. The fragmentation function has a singularity at $z=0$.

Figure 3: The coefficient of ${\rm ln}(\mu_{F}^{2}/m_{c}^{2})$ in the
fragmentation function $D_{q\to\eta_{c}}(z,\mu_{F})$, i.e.,
$F(z)=\frac{\alpha_{s}}{2\pi}\int_{z}^{1}\frac{dy}{y}P_{gq}(y)D_{g\to\eta_{c}}^{LO}(z/y)$.
The same curve is shown as a function of $z$ with a linear scale (upper one)
and with a logarithmic scale (lower one).  Figure 4: The fragmentation
function $D_{b\to\eta_{c}}(z,\mu_{F})$ as a function of $z$ for
$\mu_{F}=m_{b}$, $m_{b}+2m_{c}$ and $m_{b}+4m_{c}$.

In order to understand the dependence of the fragmentation function on the
factorization scale, we have calculated the coefficient of ${\rm
ln}(\mu_{F}^{2}/m_{c}^{2})$ in the fragmentation function
$D_{q\to\eta_{c}}(z,\mu_{F})$. The coefficient is shown as a function of $z$
in Fig.3 with a linear scale (upper one) and with a logarithmic scale (lower
one). We can see that, the coefficient decreases monotonically with the
increase of $z$ and is positive for $z\in(0,1)$. Like the fragmentation
function, the coefficient of ${\rm ln}(\mu_{F}^{2}/m_{c}^{2})$ also has a
singularity at $z=0$. When $\mu_{F}$ is very large, the fragmentation function
is dominated by the ${\rm ln}(\mu_{F}^{2}/m_{c}^{2})$ term. Therefore, the
behavior of the fragmentation function is similar to that of the coefficient
of ${\rm ln}(\mu_{F}^{2}/m_{c}^{2})$ when $\mu_{F}$ is large enough.

In Fig.4, the fragmentation function $D_{b\to\eta_{c}}(z,\mu_{F})$ for
$\mu_{F}=m_{b}$, $m_{b}+2m_{c}$ and $m_{b}+4m_{c}$ is presented. We can see
that the fragmentation function $D_{b\to\eta_{c}}(z,\mu_{F})$ is also
sensitive to the factorization scale. Actually, from Eqs.(25) and (26), we can
see that the coefficient of ${\rm ln}(\mu_{F}^{2}/m_{Q}^{2})$ in
$D_{b\to\eta_{c}}(z,\mu_{F})$ is the same as that in
$D_{q\to\eta_{c}}(z,\mu_{F})$.

###  III.2 Application to the decay widths

In this subsection, we will apply the obtained fragmentation functions to the
decay widths for the processes $Z\to\eta_{c}+q\bar{q}g$ and
$Z\to\eta_{c}+b\bar{b}g$.

Here, we only present the calculation formulas for $Z\to\eta_{c}+q\bar{q}g$,
the formulas for $Z\to\eta_{c}+b\bar{b}g$ are similar. Under the fragmentation
approach, the differential decay width $d\Gamma/dz$ for
$Z\to\eta_{c}+q\bar{q}g$ at LO can be written as

$\displaystyle\frac{d\Gamma^{\rm Frag,LO}_{Z\to\eta_{c}+q\bar{q}g}}{dz}=$
$\displaystyle 2\,\Gamma_{Z\to q+\bar{q}}D_{q\to\eta_{c}}(z,\mu_{F})$ (28)
$\displaystyle+\int^{1}_{z}\frac{dy}{y}\frac{d\hat{\Gamma}_{Z\to
g+q\bar{q}}(y,\mu_{F})}{dy}$ $\displaystyle\cdot D_{g\to\eta_{c}}(z/y),$

where the energy fraction is defined as $z\equiv E_{\eta_{c}}/E^{\rm
max}_{\eta_{c}}$, and $E_{\eta_{c}}$ and $E^{\rm max}_{\eta_{c}}$ are the
energy and the maximum energy of the $\eta_{c}$ in the rest frame of the Z
boson. $\Gamma_{Z\to q+\bar{q}}$ denotes the LO decay width for the process
$Z\to q+\bar{q}$, and the factor of 2 is due to the fact that the
contributions from the q fragmentation and $\bar{q}$ fragmentation are the
same. In Eq.(28), we have used the fact $d\hat{\Gamma}_{Z\to
q+\bar{q}}/dy=\Gamma_{Z\to q+\bar{q}}\delta(1-y)$ at LO. $d\hat{\Gamma}_{Z\to
g+q\bar{q}}/dy$ is the differential decay width for the inclusive production
of a gluon associated with a light-quark pair. Neglecting the light quark
mass, we obtain the differential decay width under the $\overline{\rm MS}$
factorization scheme

$\displaystyle\frac{d\hat{\Gamma}_{Z\to g+q\bar{q}}(y,\mu_{F})}{dy}$
$\displaystyle=\Gamma_{Z\to
q+\bar{q}}\frac{\alpha_{s}}{\pi}P_{gq}(y)\left[{\rm
ln}\frac{m_{{}_{Z}}^{2}}{\mu_{F}^{2}}+2{\rm ln}y+{\rm ln}(1-y)\right].$ (29)

$D_{g\to\eta_{c}}(z)$ and $D_{q\to\eta_{c}}(z,\mu_{F})$ are the LO
fragmentation functions. The expression of $D_{g\to\eta_{c}}(z)$ can be found
in Ref.Braaten:1993rw , and the expression of $D_{q\to\eta_{c}}(z,\mu_{F})$
has been given in Eq.(25). It is easy to check that the logarithm terms of
$\mu_{F}^{2}$ in Eq.(28) are canceled by each other which results in that
$d\Gamma^{\rm Frag,LO}_{Z\to\eta_{c}+q\bar{q}g}/dz$ is independent of
$\mu_{F}$.

The physical picture of Eq.(28) is as follows: The first term gives the
contribution from that the Z boson decays into a light quark and a light
antiquark with energies $m_{{}_{Z}}/2$ on a distance scale of order
$1/m_{{}_{Z}}$, and one of the light quark and the light antiquark decays into
an $\eta_{c}$ on a distance scale of order $1/m_{c}$. The second term gives
the contribution from that the Z boson decays into a light quark-antiquark
pair and a gluon on a distance scale of order $1/m_{{}_{Z}}$, and the gluon
decays into an $\eta_{c}$ on a distance scale of order $1/m_{c}$. The two
terms of Eq.(28) share the same Feynman diagrams which are shown in Fig.5, but
they come from different regions of the phase space. When the invariant mass
of the virtual light quark (antiquark) is very small compared to $\mu_{F}$,
the contribution is given by the first term of Eq.(28). When the invariant
mass of the virtual light quark (antiquark) is very large compared to
$\mu_{F}$, the contribution is given by the second term of Eq.(28)

Figure 5: The Feynman diagrams for $Z\to\eta_{c}+q\bar{q}g$ which are
responsible for the fragmentation mechanism.

There are logarithms of $m_{Z}/m_{c}$ in the decay widths of
$Z\to\eta_{c}+q\bar{q}g$, which may spoil the convergence of the perturbative
expansion. These large logarithms can be resummed through the evolution of the
fragmentation functions under the fragmentation approach. The decay width
after the resummation of the leading logarithms (LLs) can be written as

$\displaystyle\frac{d\Gamma^{\rm Frag,LO+LL}_{Z\to\eta_{c}+q\bar{q}g}}{dz}=$
$\displaystyle 2\,\Gamma_{Z\to q+\bar{q}}D^{\rm
LO+LL}_{q\to\eta_{c}}(z,\mu_{F})$ (30)
$\displaystyle+\int^{1}_{z}\frac{dy}{y}\frac{d\hat{\Gamma}_{Z\to
g+q\bar{q}}(y,\mu_{F})}{dy}$ $\displaystyle\cdot D^{\rm
LO+LL}_{g\to\eta_{c}}(z/y,\mu_{F}),$

where the factorization scale is set as $\mu_{F}=m_{{}_{Z}}$, and the
renormalization scale in the partonic decay widths is also set as
$\mu_{R}=m_{{}_{Z}}$, so as to avoid large logarithms appearing in the
partonic decay widths. The fragmentation functions $D^{\rm
LO+LL}_{q\to\eta_{c}}(z,\mu_{F}=m_{{}_{Z}})$ and $D^{\rm
LO+LL}_{g\to\eta_{c}}(z,\mu_{F}=m_{{}_{Z}})$ are obtained through solving the
Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations dglap1 ; dglap2 ;
dglap3  with LO splitting functions, where the initial fragmentation functions
$D_{q\to\eta_{c}}(z,\mu_{F0})$ and $D_{g\to\eta_{c}}(z,\mu_{F0})$ at
$\mu_{F0}=2m_{c}$ 222For $Z\to\eta_{c}+b\bar{b}g$ case, the initial
factorization scale is taken as $\mu_{F0}=m_{b}+2m_{c}$. are used as the
boundary condition. We solve the DGLAP equations by using the program FFEVOL
Hirai:2011si .

In addition to the fragmentation approach, we can also calculate the decay
width directly based on the NRQCD factorization, i.e.,

$\displaystyle d\Gamma^{\rm
Direct,LO}_{Z\to\eta_{c}+q\bar{q}g}=d\tilde{\Gamma}_{Z\to(c\bar{c})[^{1}S_{0}^{[1]}]+q\bar{q}g}\langle{\cal
O}^{\eta_{c}}(^{1}S_{0}^{[1]})\rangle,$ (31)

where we use “Direct” to denote the results from the direct calculation based
on the NRQCD factorization. The dominant contributions to the decay widths for
$Z\to\eta_{c}+q\bar{q}g$ and $Z\to\eta_{c}+b\bar{b}g$ come from four
fragmentation diagrams shown in Fig.5. The contributions from the
nonfragmentation diagrams are suppressed by powers of $m_{c}/m_{Z}$ compared
to the fragmentation contributions. For simplicity, under the direct NRQCD
calculation, we only consider the contributions from the four fragmentation
diagrams shown in Fig.5.

Figure 6: The differential decay width $d\Gamma/dz$ as a function of $z$ for
the process $Z\to\eta_{c}+q\bar{q}g$ under the fragmentation and the direct
NRQCD calculations. The contributions for $q=u,d,s$ are summed. Figure 7: The
differential decay width $d\Gamma/dz$ as a function of $z$ for
$Z\to\eta_{c}+b\bar{b}g$ under the fragmentation and the direct NRQCD
calculations. Figure 8: The ratios $(d\Gamma_{\rm Frag,LO}/dz)/(d\Gamma_{\rm
Direct,LO}/dz)$ as functions of $z$ for the processes $Z\to\eta_{c}+q\bar{q}g$
and $Z\to\eta_{c}+b\bar{b}g$.

The differential decay widths for $Z\to\eta_{c}+q\bar{q}g$ and
$Z\to\eta_{c}+b\bar{b}g$ under the fragmentation and the direct NRQCD
approaches are presented in Figs.6 and 7. In order to see the difference
between the “Frag,LO” results and the “Direct,LO” results more clearly, the
ratios $(d\Gamma_{\rm Frag,LO}/dz)/(d\Gamma_{\rm Direct,LO}/dz)$ are given in
Fig.8. For the “Frag,LO” and the “Direct,LO” calculations, the renormalization
scale is fixed as $\mu_{R}=2m_{c}$ for simplicity. For the “Frag,LO+LL”
calculation, the choice of the renormalization scale and the factorization
scale has been described below Eq.(30). From the figures, we can see that the
differential decay widths from the “Frag,LO” calculation are very close to
those from the “Direct,LO” calculation, especially for $0.2\leq z\leq 0.8$.

The total decay widths can be obtained through integrating the differential
decay widths $d\Gamma/dz$ over $z$. The total decay widths under the
fragmentation approach and the direct NRQCD approach are given in Table 1. We
can see that the total decay widths obtained from the “Frag,LO” calculation
and the “Direct,LO” calculation are also very close. Therefore, the fixed-
order fragmentation approach (i.e., the “Frag,LO” approach) provides a good
approximation to the direct NRQCD calculation.

The differential and total decay widths after the resummation of the large
logarithms under the fragmentation approach are also shown in Figs.6 and 7 and
Table 1. We can see that, after the resummation, the differential decay widths
are enhanced at smaller $z$ values but are reduced at larger $z$ values, and
the total decay widths are reduced after the resummation. In fact, the fixed-
order results have a big uncertainty caused by the choice of the
renormalization scale. If we set the renormalization scale as
$\mu_{R}=m_{{}_{Z}}$, the fixed-order results become
$[\alpha_{s}(m_{{}_{Z}})/\alpha_{s}(2m_{c})]^{3}=0.0958$ of those with
$\mu_{R}=2m_{c}$. On the contrary, the results after the resummation have a
smaller uncertainty caused by the choice of the renormalization scale. Because
the renormalization scale of the initial fragmentation functions should be
${\cal O}(m_{c})$, and the renormalization scale of the coefficient functions
should be ${\cal O}(m_{{}_{Z}})$. Moreover, the resummed results include the
leading logarithms up to all orders from the collinear radiation. Therefore,
the results after the resummation are more precise than the fixed-order
results.

|  $\Gamma_{Z\to\eta_{c}+q\bar{q}g}$ (keV) |  $\Gamma_{Z\to\eta_{c}+b\bar{b}g}$(keV)
---|---|---
Direct | 63.0 | 19.3
Frag,LO | 65.6 | 20.3
Frag,LO+LL | 40.3 | 15.3
Table 1: The decay widths of $Z\to\eta_{c}+q\bar{q}g$ and
$Z\to\eta_{c}+b\bar{b}g$ under the fragmentation and the direct NRQCD
approaches. For the $Z\to\eta_{c}+q\bar{q}g$ case, the contributions for
$q=u,d,s$ are summed.

##  IV Summary

In the present paper, we have calculated the fragmentation functions for a
(heavy or light) quark into a spin-singlet quarkonium, where the flavor of the
initial quark is different from that of the constituent quark in the
quarkonium. There are UV divergences in the phase-space integral, which are
removed through the operator renormalization of the fragmentation function. We
have carried out the renormalization under the $\overline{\rm MS}$ scheme. The
fragmentation function $D_{q\to\eta_{Q}}(z,\mu_{F})$ is given as a two-
dimensional integral, and this two-dimensional integral can be calculated
easily through numerical integration. Numerical results for a light quark or a
bottom quark into the $\eta_{c}$ with several factorization scales are
analyzed. The results show that these fragmentation functions are sensitive to
the factorization scale. Especially, when $\mu_{F}$ is small, the
fragmentation functions are negative at small $z$ values. There is a
singularity at $z=0$ for these fragmentation functions.

We have applied the obtained fragmentation functions to the decay widths for
the processes $Z\to\eta_{c}+q\bar{q}g(q=u,d,s)$ and $Z\to\eta_{c}+b\bar{b}g$.
The differential decay widths and total decay widths are calculated under the
fragmentation and the direct NRQCD approaches. It is found that the results
under the fixed-order fragmentation and the direct NRQCD approaches are close
to each other. Therefore, the fixed-order fragmentation approach provides a
good approximation to the direct NRQCD calculation. The more precise results
containing the resummation of the large logarithms under the fragmentation
approach are also presented. Moreover, the fragmentation functions obtained in
this paper can be used in the studies on the production of $\eta_{c}$ and
$\eta_{b}$ at high-energy colliders.

Acknowledgments: This work was supported in part by the Natural Science
Foundation of China under Grants No. 11625520, No. 12005028 and No. 12047564,
by the Fundamental Research Funds for the Central Universities under Grant No.
2020CQJQY-Z003, and by the Chongqing Graduate Research and Innovation
Foundation under Grant No. ydstd1912.

##  Appendix A The parametrization for the phase space

In order to extract the UV poles in the calculation analytically, the phase
space for the fragmentation function should be parametrized properly. In this
Appendix, we will present a parametrization for the phase space. The
differential phase space for the fragmentation function
$D_{q\to(Q\bar{Q})[^{1}S_{0}^{[1]}]}$ is

$\displaystyle d\phi_{3}(p_{1},p_{2},p_{3})=$ $\displaystyle
2\pi\delta\left(K^{+}-\sum_{i=1}^{3}p_{i}^{+}\right)\mu^{2(4-d)}$ (32)
$\displaystyle\times\prod_{i=2,3}\frac{\theta(p_{i}^{+})dp_{i}^{+}}{4\pi
p_{i}^{+}}\frac{d^{d-2}\textbf{p}_{i\perp}}{(2\pi)^{d-2}}.$

According to Ref.Zheng:2019gnb , the differential phase space for a single
parton with momentum $p_{i}$ and mass $m_{i}$ can be expressed as

$\displaystyle\frac{d^{d-1}\textbf{p}_{i}}{(2\pi)^{d-1}2p_{i}^{0}}=\frac{(\lambda_{i}p_{i}\cdot
n-m_{i}^{2})^{-\epsilon}}{4(2\pi)^{3-2\epsilon}}d\lambda_{i}d(p_{i}\cdot
n)\,d\Omega_{i\perp},$ (33)

where

$\displaystyle\lambda_{i}=2k_{i}\cdot p_{i}/k_{i}\cdot n,$ (34)

and $k_{i}$ is an arbitrary lightlike momentum which is not parallel to $n$.
$d\Omega_{i\perp}$ stands for the differential transverse solid angle, and the
total transverse solid angle
$\Omega_{i\perp}=2\pi^{1-\epsilon}/\Gamma(1-\epsilon)$.

Applying the parametrization (33) to the differential phase spaces for $p_{2}$
and $p_{3}$, we obtain

$\displaystyle d\phi_{3}(p_{1},p_{2},p_{3})=$
$\displaystyle\frac{2^{-2\epsilon}(K\cdot
n)^{1-2\epsilon}\mu^{4\epsilon}}{(4\pi)^{4-3\epsilon}\Gamma(1-\epsilon)}[(1-y)(y-z)]^{-\epsilon}$
(35) $\displaystyle\times\left[1-\frac{m_{q}^{2}}{\lambda_{3}(1-y)K\cdot
n}\right]^{-\epsilon}\lambda_{2}^{-\epsilon}\lambda_{3}^{-\epsilon}$
$\displaystyle\times dy\,d\lambda_{2}\,d\lambda_{3}\,d\Omega_{3\perp},$

where the integrations over $p_{2}\cdot n$ and $\Omega_{2\perp}$ have been
performed.

To obtain the phase-space parametrization used to extract the UV poles, we
choose the light-like momenta $k_{2}$ and $k_{3}$ as follows:

$\displaystyle k_{2}^{\mu}$ $\displaystyle=$ $\displaystyle
p_{1}^{\mu}-\frac{2m_{Q}^{2}}{p_{1}\cdot n}n^{\mu},$ $\displaystyle
k_{3}^{\mu}$ $\displaystyle=$
$\displaystyle(p_{1}+p_{2})^{\mu}-\frac{s_{1}}{2(p_{1}+p_{2})\cdot n}n^{\mu},$
(36)

then we obtain

$\displaystyle\lambda_{2}=\frac{1}{zK\cdot
n}\left(s_{1}-\frac{y}{z}4m_{Q}^{2}\right),$ (37)

and

$\displaystyle\lambda_{3}=\frac{1}{yK\cdot
n}\left(s-\frac{s_{1}}{y}-m_{q}^{2}\right),$ (38)

Changing variables in Eq.(35) from $\lambda_{2}$ and $\lambda_{3}$ to $s_{1}$
and $s$, we obtain

$\displaystyle d\phi_{3}(p_{1},p_{2},p_{3})=$
$\displaystyle\frac{2^{-2\epsilon}(z\,y)^{-1+\epsilon}\mu^{4\epsilon}}{(4\pi)^{4-3\epsilon}\Gamma(1-\epsilon)K\cdot
n}(1-y)^{-\epsilon}$ (39)
$\displaystyle\times(y-z)^{-\epsilon}[s-s_{1}/y-m_{q}^{2}/(1-y)]^{-\epsilon}$
$\displaystyle\times(s_{1}-4m_{Q}^{2}y/z)^{-\epsilon}dy\,ds\,ds_{1}\,d\Omega_{3\perp}.$

##  Appendix B The expression of $f(s_{1},y,z)$

The expression of $f(s_{1},y,z)$ which appeared in Eq.(20) can be written in
following form:

$\displaystyle
f(s_{1},y,z)=\frac{2^{8}g_{s}^{6}[(1-y)s_{1}+y^{2}m_{q}^{2}]^{-\epsilon}(K\cdot
n)^{2}}{9m_{Q}s_{1}^{2}(s_{1}-4m_{Q}^{2})^{2}}(f_{0}+\epsilon f_{1}),$ (40)

where

$\displaystyle
f_{0}=[(1-y)^{2}+1][s_{1}^{2}(y^{2}-2y\,z+2z^{2})-8m_{Q}^{2}s_{1}y\,z+16y^{2}m_{Q}^{4}],$
(41)

and

$\displaystyle f_{1}=$
$\displaystyle\frac{1}{[(1-y)s_{1}+y^{2}m_{q}^{2}]}\Big{\\{}s_{1}^{2}\Big{[}s_{1}(y-1)\big{(}5y^{4}-8y^{3}(z+1)+4y^{2}(z+2)(2z+1)-20y\,z(z+1)+20z^{2}\big{)}$
(42)
$\displaystyle-2y^{2}m_{q}^{2}(y^{2}-2y+2)(2y^{2}-3y\,z+3z^{2})\Big{]}-8ym_{Q}^{2}s_{1}\Big{[}s_{1}(y-1)(y^{3}+2y^{2}(2z+1)-2y(5z+1)+10z)$
$\displaystyle-y^{2}m_{q}^{2}(y^{2}-2y+2)(y+3z)\Big{]}+16m_{Q}^{4}y^{2}\Big{[}s_{1}(5y^{3}-13y^{2}+16y-8)-4m_{q}^{2}y^{2}(y^{2}-2y+2)\Big{]}\Big{\\}}.$

## References

  * (1) J. C. Collins, D. E. Soper and G. F. Sterman, Factorization of Hard Processes in QCD, Adv. Ser. Direct. High Energy Phys. 5, 1-91 (1989).
  * (2) J. C. Collins and D. E. Soper, Parton Distribution and Decay Functions, Nucl. Phys. B 194, 445-492 (1982).
  * (3) G. C. Nayak, J. W. Qiu and G. F. Sterman, Fragmentation, NRQCD and NNLO factorization analysis in heavy quarkonium production, Phys. Rev. D 72, 114012 (2005).
  * (4) Z. B. Kang, J. W. Qiu and G. Sterman, Factorization and quarkonium production, Nucl. Phys. B Proc. Suppl. 214, 39-43 (2011).
  * (5) Z. B. Kang, J. W. Qiu and G. Sterman, Heavy quarkonium production and polarization, Phys. Rev. Lett. 108, 102002 (2012).
  * (6) S. Fleming, A. K. Leibovich, T. Mehen and I. Z. Rothstein, The Systematics of Quarkonium Production at the LHC and Double Parton Fragmentation, Phys. Rev. D 86, 094012 (2012).
  * (7) S. Fleming, A. K. Leibovich, T. Mehen and I. Z. Rothstein, Anomalous dimensions of the double parton fragmentation functions, Phys. Rev. D 87, 074022 (2013).
  * (8) G. T. Bodwin, E. Braaten and G. P. Lepage, Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium, Phys. Rev. D 51, 1125 (1995) [Erratum-ibid. D 55, 5853 (1997)].
  * (9) C. H. Chang and Y. Q. Chen, The Production of B(c) or anti-B(c) meson associated with two heavy quark jets in Z0 boson decay, Phys. Rev. D 46, 3845 (1992), [erratum: Phys. Rev. D 50, 6013 (1994)].
  * (10) E. Braaten, K. m. Cheung and T. C. Yuan, Perturbative QCD fragmentation functions for $B_{c}$ and $B_{c}$* production, Phys. Rev. D 48, 5049 (1993).
  * (11) E. Braaten, K. Cheung and T. C. Yuan, Z0 decay into charmonium via charm quark fragmentation, Phys. Rev. D 48, 4230-4235 (1993).
  * (12) E. Braaten and T. C. Yuan, Gluon fragmentation into heavy quarkonium, Phys. Rev. Lett. 71, 1673-1676 (1993).
  * (13) E. Braaten and T. C. Yuan, Gluon fragmentation into P wave heavy quarkonium, Phys. Rev. D 50, 3176-3180 (1994).
  * (14) E. Braaten and T. C. Yuan, Gluon fragmentation into spin triplet S wave quarkonium, Phys. Rev. D 52, 6627-6629 (1995).
  * (15) Y. Q. Chen, Perturbative QCD predictions for the fragmentation functions of the P wave mesons with two heavy quarks, Phys. Rev. D 48, 5181-5189 (1993).
  * (16) T. C. Yuan, Perturbative QCD fragmentation functions for production of P wave mesons with charm and beauty, Phys. Rev. D 50, 5664-5675 (1994).
  * (17) J. P. Ma, Calculating fragmentation functions from definitions, Phys. Lett. B 332, 398-404 (1994).
  * (18) J. P. Ma, Gluon fragmentation into P wave triplet quarkonium, Nucl. Phys. B 447, 405-424 (1995).
  * (19) J. P. Ma, Quark fragmentation into p wave triplet quarkonium, Phys. Rev. D 53, 1185-1190 (1996).
  * (20) P. L. Cho, M. B. Wise and S. P. Trivedi, Gluon fragmentation into polarized charmonium, Phys. Rev. D 51, 2039-2043 (1995).
  * (21) M. Beneke and I. Z. Rothstein, Psi-prime polarization as a test of color octet quarkonium production, Phys. Lett. B 372, 157-164 (1996), [erratum: Phys. Lett. B 389, 769 (1996)].
  * (22) E. Braaten and J. Lee, Next-to-leading order calculation of the color octet 3S(1) gluon fragmentation function for heavy quarkonium, Nucl. Phys. B 586, 427-439 (2000).
  * (23) W. l. Sang, L. f. Yang and Y. q. Chen, Relativistic corrections to heavy quark fragmentation to S-wave heavy mesons, Phys. Rev. D 80, 014013 (2009).
  * (24) G. Hao, Y. Zuo and C. F. Qiao, The Fragmentation Function of Gluon Splitting into P-wave Spin-singlet Heavy Quarkonium, [arXiv:0911.5539 [hep-ph]].
  * (25) Y. Jia, W. L. Sang and J. Xu, Inclusive $h_{c}$ Production at $B$ Factories, Phys. Rev. D 86, 074023 (2012).
  * (26) G. T. Bodwin, H. S. Chung, U. R. Kim and J. Lee, Quark fragmentation into spin-triplet $S$-wave quarkonium, Phys. Rev. D 91, 074013 (2015).
  * (27) Y. Q. Ma, J. W. Qiu and H. Zhang, Heavy quarkonium fragmentation functions from a heavy quark pair. I. $S$ wave, Phys. Rev. D 89, 094029 (2014).
  * (28) Y. Q. Ma, J. W. Qiu and H. Zhang, Fragmentation functions of polarized heavy quarkonium, JHEP 06, 021 (2015).
  * (29) D. Yang and W. Zhang, Relativistic corrections of the fragmentation functions for a heavy quark to $B_{c}$ and $B_{c}^{*}$, Chin. Phys. C 43, 083101 (2019).
  * (30) P. Artoisenet and E. Braaten, Gluon fragmentation into quarkonium at next-to-leading order, JHEP 04, 121 (2015).
  * (31) R. Sepahvand and S. Dadfar, NLO corrections to $c$\- and $b$-quark fragmentation into $j/\psi$ and $\gamma$, Phys. Rev. D 95, 034012 (2017).
  * (32) P. Artoisenet and E. Braaten, Gluon fragmentation into quarkonium at next-to-leading order using FKS subtraction, JHEP 01, 227 (2019).
  * (33) F. Feng and Y. Jia, Next-to-leading-order QCD corrections to gluon fragmentation into ${}^{1}S_{0}^{(1,8)}$ quarkonia, [arXiv:1810.04138 [hep-ph]].
  * (34) P. Zhang, C. Y. Wang, X. Liu, Y. Q. Ma, C. Meng and K. T. Chao, Semi-analytical calculation of gluon fragmentation into1S${}_{0}^{[1,8]}$ quarkonia at next-to-leading order, JHEP 04, 116 (2019).
  * (35) X. C. Zheng, C. H. Chang and X. G. Wu, NLO fragmentation functions of heavy quarks into heavy quarkonia, Phys. Rev. D 100, 014005 (2019).
  * (36) X. C. Zheng, C. H. Chang, T. F. Feng and X. G. Wu, QCD NLO fragmentation functions for c or $\bar{b}$ quark to Bc or Bc* meson and their application, Phys. Rev. D 100, 034004 (2019).
  * (37) F. Feng, S. Ishaq, Y. Jia and J. Y. Zhang, Fragmentation function of gluon into spin-singlet $P$-wave quarkonium, Phys. Rev. D 102, 014038 (2020).
  * (38) R. Mertig, M. Bohm and A. Denner, FEYN CALC: Computer algebraic calculation of Feynman amplitudes, Comput. Phys. Commun. 64, 345-359 (1991).
  * (39) V. Shtabovenko, R. Mertig and F. Orellana, New Developments in FeynCalc 9.0, Comput. Phys. Commun. 207, 432-444 (2016).
  * (40) A. H. Mueller, Cut Vertices and their Renormalization: A Generalization of the Wilson Expansion, Phys. Rev. D 18, 3705 (1978).
  * (41) E.J. Eichten and C. Quigg, Mesons with beauty and charm: Spectroscopy, Phys.Rev. D49, 5845 (1994) and references therein.
  * (42) Y.L. Dokshitzer, Calculation of the Structure Functions for Deep Inelastic Scattering and $e^{+}e^{-}$ Annihilation by Perturbation Theory in Quantum Chromodynamics, Sov. Phys. JETP 46, 641 (1977);Zh.Eksp.Teor.Fiz. 73, 1216 (1977).
  * (43) V.N. Gribov and L.N. Lipatov, Deep inelastic ep scattering in perturbation theory, Sov. J. Nucl. Phys. 15, 438 (1972);Yad.Fiz. 15, 781 (1972).
  * (44) G. Altarelli and G. Parisi, Asymptotic Freedom in Parton Language, Nucl. Phys. B 126, 298 (1977).
  * (45) M. Hirai and S. Kumano, Numerical solution of $Q^{2}$ evolution equations for fragmentation functions, Comput. Phys. Commun. 183, 1002-1013 (2012).