The NLO fragmentation functions of heavy quarks into heavy quarkonia

TL;DR

This paper derives the next-to-leading order fragmentation functions for heavy quarks into heavy quarkonia, providing a tool for precise predictions of quarkonium production at high-energy colliders.

Contribution

It presents the first NLO calculation of heavy quark to quarkonium fragmentation functions, including operator renormalization and initial scale determination.

Findings

Fragmentation functions are obtained at NLO accuracy.

The functions can be evolved to any scale using DGLAP equations.

Application to $J/ ext{}/ ext{Psi}$ and $ ext{Upsilon}$ production at $Z$-pole colliders.

Abstract

In the paper, we derive the next-to-leading order (NLO) fragmentation function for a heavy quark, either charm or bottom, into a heavy quarkonium $J/\Psi$ or $\Upsilon$. The ultra-violet divergences in the real corrections are removed through the operator renormalization, which is performed under the modified minimal subtraction scheme. We then obtain the NLO fragmentation function at an initial factorization scale, e.g. $\mu_{F}=3 m_c$ for $c\to J/\Psi$ and $\mu_{F}=3m_b$ for $b\to \Upsilon$, which can be evolved to any scale via the use of Dokshitzer-Gribov-Lipatov-Altarelli-Parisi equation. As an initial application of those fragmentation functions, we study the $J/\Psi$ ($\Upsilon$) production at a high luminosity $e^+e^-$ collider which runs at the energy around the $Z$ pole and could be a suitable platform for testing the fragmentation function.

Tables

Table 1. Table 1: The fragmentation probability and the average value of the energy fraction for c → J / Ψ → 𝑐 𝐽 Ψ c\to J/\Psi . μ F = 3 ​ m c subscript 𝜇 𝐹 3 subscript 𝑚 𝑐 \mu_{F}=3m_{c} .

	$P^{\mathrm{LO}}\times {10}^4$	$P^{\mathrm{NLO}}\times {10}^4$	${⟨z⟩}^{\mathrm{LO}}$	${⟨z⟩}^{\mathrm{NLO}}$
${\mu }_R=2m_c$	1.88	1.86	0.62	0.63
${\mu }_R=3m_c$	1.40	1.72	0.62	0.63

Table 2. Table 2: The fragmentation probability and the averaged energy fraction for b → Υ → 𝑏 Υ b\to\Upsilon . μ F = 3 ​ m b subscript 𝜇 𝐹 3 subscript 𝑚 𝑏 \mu_{F}=3m_{b} .

	$P^{\mathrm{LO}}\times {10}^5$	$P^{\mathrm{NLO}}\times {10}^5$	${⟨z⟩}^{\mathrm{LO}}$	${⟨z⟩}^{\mathrm{NLO}}$
${\mu }_R=2m_b$	2.09	2.05	0.62	0.63
${\mu }_R=3m_b$	1.73	1.99	0.62	0.63

Table 3. Table 3: Total cross sections (in unit: p ​ b 𝑝 𝑏 pb ) of production channels e + ​ e − → J / Ψ + X → superscript 𝑒 superscript 𝑒 𝐽 Ψ 𝑋 e^{+}e^{-}\to J/\Psi+X and e + ​ e − → Υ + X → superscript 𝑒 superscript 𝑒 Υ 𝑋 e^{+}e^{-}\to\Upsilon+X at the Z 𝑍 Z pole.

	Frag, LO	Frag, NLO	Frag, NLO+NLL
$J/\Psi$	2.58	2.77	2.65
${\rm Y}$	0.368	0.382	0.377

Equations

$$d\sigma (A+B\to H+X)$$

$$=\sum _{n}d\tilde{\sigma }(A+B\to (Q\bar{Q})[n]+X)⟨{\mathcal{O}}_n^H⟩,$$

$$d\sigma (e^{+}{e}^{-}\to H(p)+X)=$$

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

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

$$=\frac{{\alpha }_s({\mu }_F)}{2\pi }\sum _j{\int }_z^1\frac{dy}{y}{P}_{ji}(y,{\alpha }_s({\mu }_F))D_{j\to H}(z/y,{\mu }_F),$$

$$P_{ji}(y,{\alpha }_s({\mu }_F))=$$

$$P_{QQ}^{(0)}(y)=C_F\left[\frac{1+y^2}{(1-y{)}_{+}}+\frac{3}{2}\delta (1-y)\right],$$

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

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

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

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

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

$$\Pi =-\frac{1}{2\sqrt{2m_Q}}/ϵ(p_1)(/p_1+2m_Q),$$

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

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

$$a_2=$$

$$a_3=$$

$$a_4=$$

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

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

$$D_{Q\to Q\bar{Q}{[}^3{S}_1^{[1]}]}^{\mathrm{LO}}(z)=N_{CS}\int d{\varphi }_{\mathrm{Born}}{\mathcal{A}}_{\mathrm{Born}},$$

$$D_{Q\to Q\bar{Q}{[}^3{S}_1^{[1]}]}^{\mathrm{LO}}(z)$$

$$=\frac{C_F^2{\alpha }_s^{2}z(1-z)(4\pi {)}^{ϵ}\Gamma (1+ϵ)}{2N_c(2-z{)}^{4+2ϵ}{m}_Q^{3+2ϵ}}[a_2+a_3\frac{(1+ϵ)z(1-z)}{2(2-z{)}^2}$$

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

$$D_{Q\to Q\bar{Q}{[}^3{S}_1^{[1]}]}^{\mathrm{LO}}(z)$$

$$=\frac{32{\alpha }_s^{2}z(1-z{)}^2}{27(2-z{)}^6{m}_Q^3}(5z^4-32z^3+72z^2-32z+16)$$

$$\times \frac{⟨{\mathcal{O}}^{Q\bar{Q}{[}^3{S}_1^{[1]}]}{(}^3{S}_1^{[1]})⟩}{6N_c}.$$

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

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

$$D_{Q\to H}^{\mathrm{LO}}(z)$$

$$=\frac{8{\alpha }_s^{2}z(1-z{)}^2|R_S(0){|}^2}{27\pi (2-z{)}^6{m}_Q^3}(5z^4-32z^3+72z^2-32z+16),$$

$$D_{Q\to Q\bar{Q}{[}^3{S}_1^{[1]}]}^{\mathrm{virtual}}(z)=N_{CS}\int d{\varphi }_{\mathrm{Born}}{\mathcal{A}}_{\mathrm{virtual}}.$$

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

$$D_{Q\to Q\bar{Q}{[}^3{S}_1^{[1]}]}^{\mathrm{real}}(z)=N_{CS}\int d{\varphi }_{\mathrm{real}}{\mathcal{A}}_{\mathrm{real}},$$

Full text

I Introduction

Since the observation of $J/\Psi$ meson, the heavy quarkonium has attracted great interests from theorists and experimentalists. Due to the constituent quark ($Q$) and antiquark ($\bar{Q}$) of a heavy quarkonium are heavy, i.e, $m_{Q}\gg\Lambda_{QCD}$, its production rate and decay width involve both perturbative and non-perturbative aspects of quantum chromodynamics (QCD). It then provides a good platform for testing various QCD factorization theories. Within the framework of the nonrelativistic QCD (NRQCD) factorization theory nrqcd , the production cross section of a heavy quarkonium via the collision of two incident particles $A$ and $B$ can be written as

[TABLE]

where $H$ denotes the produced heavy quarkonium, $d\tilde{\sigma}$ is the production cross section for the perturbative state $(Q\bar{Q})[n]$ with quantum numbers $n$, which is calculable and can be expanded in powers of the strong coupling constant $\alpha_{s}$. $\langle{\cal O}_{n}^{H}\rangle$ denotes the non-perturbative but universal NRQCD matrix element, which is proportional to the transition probability of the perturbative state $(Q\bar{Q})[n]$ into the hadron state $H$.

The NRQCD factorization formulism has been used to deal with the quarkonium production at $e^{+}e^{-}$ and hadronic colliders, and most of the calculations have been performed up to next-to-leading order (NLO) accuracy ybook1 ; ybook2 . In some cases, there are large logarithms in the short-distance part $d\tilde{\sigma}$. For instance, there are large logarithms in powers of $\ln(\sqrt{s}/m_{Q})$ for the heavy quarkonium production at a $e^{+}e^{-}$ collider with $\sqrt{s}$ being the central of mass collision energy, or there are large logarithms in powers of $\ln(p_{T}/m_{Q})$ in high $p_{T}$ region ($p_{T}$ being the transverse momentum of the quarkonium) at a hadronic collider. Those large logarithms may spoil the convergence of the perturbative expansion, leading to the unreliable pQCD predictions. As a solution, it has been pointed out that those logarithms mainly come from the emission of collinear gluons, which can be systematically resummed through the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation dglap .

In the present paper, we shall take the heavy quarkonium production at an $e^{+}e^{-}$ collider as an explicit example to explain this idea, especially, we shall give the results for the heavy quark to heavy quarkonium fragmentation function up to NLO accuracy. The fragmentation function that gives the probability for the splitting of a parton into the desired hadron plus other partons. For a hadron composed of only heavy quarks, its fragmentation function can be calculated using pQCD theory. The Z boson decays $Z\to J/\Psi+X$ and $Z\to\Upsilon+X$ have been studied up to NLO level using the complete pQCD approach JXWang , and in the present paper we will study these processes but using the fragmentation approach. Compared with the complete pQCD approach, even though the fragmentation approach neglects some power-suppressed terms, it resums the large collinear logarithms more conveniently via the use of the DGLAP evolution equation and achieves more reliable predictions in specific kinematic regions.

We shall also study these processes but using the fragmentation function approach. Compared with the complete pQCD approach, the fragmentation function approach can resum the large collinear logarithms, however, it neglects some power suppressed terms. One may combine these two approaches to give a good prediction.

Under the pQCD factorization theory, the production cross section of heavy quarkonium at an $e^{+}e^{-}$ collider can be factorized as

[TABLE]

where the sum extends over all the parton types, $d\hat{\sigma}$ denotes the partonic cross section (coefficient function), $D_{i\to H}(z,\mu_{F})$ denotes the fragmentation function (decay function) for the parton $i$ into heavy quarkonium $H$ with longitudinal momentum fraction $z$. $\mu_{F}$ is the factorization scale which separates the energy scales of two parts. In order to avoid large logarithms appearing in $d\hat{\sigma}$, $\mu_{F}$ is usually set as $\mu_{F}={\cal O}(\sqrt{s})$.

The pQCD factorization formula (2) was firstly suggested by Collins and Soper for the inclusive production of a light hadron Collins , and the proof of the pQCD factorization formula to the case of quarkonium production was given by Nayak, Qiu, and Sterman fragnrqcd . The recent progress of the pQCD factorization for quarkonium production is the derivation of the next-leading-power (NLP) contribution, which comes from the double-parton fragmentation nlp1 ; nlp2 ; nlp3 ; nlp4 . The fragmentation function $D_{i\to H}(z,\mu_{F})$ contains nonperturbative information, which is calculable through the NRQCD factorization (or the Mandelstam formulation Mandelstam under the instantaneous approximation) fragbc1 ; fragbc2 , e.g. $D_{i\to H}(z,\mu_{F})$ can be factorized as

[TABLE]

where $d_{i\to(Q\bar{Q})[n]}(z,\mu_{F})$ is the short-distance coefficient, which contains the logarithms of $\mu_{F}/m_{Q}$. To avoid such kind of large logarithms, one can first calculate the fragmentation function at some initial factorization scale which is of order ${\cal O}(m_{Q})$, and then evolve it to a higher factorization scale by using the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation dglap1 ; dglap2 ; dglap3 :

[TABLE]

where $P_{ji}$ are splitting functions, which can be expanded in perturbative series as:

[TABLE]

For the quark to quark case, the LO coefficient

[TABLE]

where $C_{F}=4/3$ for $SU(3)_{c}$ group. The NLO coefficient $P^{(1)}_{QQ}(y)$ is too lengthy to be presented here, which can be found in Refs.nlospfun1 ; nlospfun2 ; nlospfun3 .

The LO fragmentation function for a heavy quark to $J/\Psi$ or $\Upsilon$ was firstly calculated by Braaten and Cheung in 1993 Braaten1 , where the LO fragmentation function is derived from a LO calculation of the process $Z\to J/\Psi+c+\bar{c}$. Subsequently, Ma calculated the LO fragmentation functions Ma by using the gauge invariant definition suggested by Collins and Soper. In Ref.Sepahvand , the authors calculated the NLO corrections to the transverse-momentum dependent fragmentation functions for a heavy quark to $J/\Psi$ and $\Upsilon$, which is however not convenient for practical applications. Recently, the NLO fragmentation function for a gluon into heavy quarkonium has been finished by Refs.Braaten2 ; Braaten3 ; YJia ; YQMa . In the present paper, we shall give the fragmentation functions for a heavy quark into $J/\Psi$ and $\Upsilon$ up to NLO level.

The remaining parts of the paper are organized as follows. In Sec.II, we present the LO fragmentation function for a heavy quark into the heavy quarkonium $J/\Psi$ or $\Upsilon$. In Sec.III, we present the NLO fragmentation functions $D_{c\to J/\Psi}$ and $D_{b\to\Upsilon}$, in which the renormalization is carried out by using the conventional $\overline{\rm MS}$-scheme. In Sec.IV, we apply those NLO fragmentation functions to $J/\Psi$ and $\Upsilon$ production at a super $Z$ factory. Sec.V is reserved for a summary.

II The LO fragmentation function

Before carrying out the calculation for the fragmentation function, we first give a brief introduction of the gauge invariant fragmentation function suggested by Collins and Soper Collins . We adopt the dimensional regularization to regularize the infrared (IR) and ultraviolet (UV) divergences, and shall work in $d=4-2\epsilon$ dimensional space-time.

The light-cone coordinates are conventionally adopted to define the fragmentation function, where a vector $V^{\mu}$ in $d$ dimensions is expressed as $V^{\mu}=(V^{+},V^{-},\textbf{V}_{\perp})=((V^{0}+V^{d-1})/\sqrt{2},(V^{0}-V^{d-1})/\sqrt{2},\textbf{V}_{\perp})$. The scalar product of two vectors $V$ and $W$ then becomes $V\cdot W=V^{+}W^{-}+V^{-}W^{+}-\textbf{V}_{\perp}\cdot\textbf{W}_{\perp}$. The gauge-invariant fragmentation function for the heavy quark $Q$ into a spin-triplet and color-singlet quarkonium $H$ in $d=4-2\epsilon$ dimension is defined as

[TABLE]

where $\Psi$ and $A_{a}^{\mu}$ are quark field and gluon field, respectively. $t^{a}$ is the color matrix, ${\cal P}$ implies the path ordering, $z\equiv P^{+}/K^{+}$ is the longitudinal momentum fraction carried from the incident heavy quark $Q$. This definition is carried out in a reference frame where the quarkonium $H$ carries no transverse momentum, i.e, $P^{\mu}=(P^{+},P^{-}=m_{H}^{2}/2P^{+},0_{\perp})$. It is convenient to introduce a light-like vector $n^{\mu}$, which has the value of $n^{\mu}=(0,1,0_{\perp})$ in the reference frame where the definition of the fragmentation function is carried out. In this frame, the plus component of a momentum $p$ can be expressed as $p^{+}=p\cdot n$, and $z=P\cdot n/K\cdot n$. The Feynman rules for the fragmentation function can be directly derived from the above definition Bcfragnlo .

The fragmentation function (7) is gauge independent. For practical treatment, we work under the usual Feynman gauge. To derive the fragmentation function $D_{Q\to H}$, we first calculate the fragmentation function for the production of a free on-shell $Q\bar{Q}$ state with the quantum numbers ${}^{3}S^{[1]}_{1}$ (i.e, $D_{Q\to Q\bar{Q}[^{3}S_{1}^{[1]}]}$), where the superscript “[1]” denotes the $Q\bar{Q}$ state is in color singlet. Then, $D_{Q\to H}$ is obtained by replacing the NRQCD matrix element of a free $Q\bar{Q}$ state ($\langle{\cal O}^{Q\bar{Q}[^{3}S_{1}^{[1]}]}(^{3}S_{1}^{[1]})\rangle$) to the matrix element of the quarkonium ($\langle{\cal O}^{H}(^{3}S_{1}^{[1]})\rangle$).

At the LO level, we need to deal with process $Q(K)\to Q\bar{Q}[^{3}S_{1}^{[1]}](p_{1})+Q(p_{2})$. As shown by Fig.1, there are four LO cut diagrams which contribute to $D_{Q\to Q\bar{Q}[^{3}S_{1}^{[1]}]}$. The four squared amplitudes from four cut diagrams are

[TABLE]

[TABLE]

[TABLE]

Here $\Pi$ denotes the spin projector, and for ${}^{3}S_{1}$ state

[TABLE]

and $\bar{\Pi}\equiv\gamma^{0}\Pi^{\dagger}\gamma^{0}$. $\Lambda_{1}$ is color-singlet prejector

[TABLE]

where 1 is the ${\rm SU}_{c}(3)$ unit matrix.

Then we obtain the total squared amplitude at the LO level,

[TABLE]

where $s_{1}=(p_{1}+p_{2})^{2}$ is the invariant mass of final $Q\bar{Q}[^{3}S^{[1]}_{1}]+Q$, and the coefficients $a_{j=2,3,4}$ are

[TABLE]

The differential phase space for the LO fragmentation function is

[TABLE]

where the $\delta$-function comes from the final cut of the eikonal line. The integration over $p_{2}^{+}$ can be carried out with the $\delta$-function. The integration over the angles of $\textbf{p}_{2\perp}$ is trivial and can be carried out easily. Then we have

[TABLE]

The range of $s_{1}$ is from $[4m_{Q}^{2}/z+m_{Q}^{2}/(1-z)]$ to $+\infty$.

The LO fragmentation function for $Q\to Q\bar{Q}[^{3}S^{[1]}_{1}]$ can be obtained through

[TABLE]

where $N_{CS}=z^{1-2\epsilon}/8\pi N_{c}$ is an overall factor. Performing the integration over $s_{1}$, we obtain

[TABLE]

Setting $d=4$, we obtain

[TABLE]

Here, the LO fragmentation function for the free $Q\bar{Q}[^{3}S_{1}^{[1]}]$ state has been written as the factorization form by the use of the fact that

[TABLE]

at the order $\alpha_{s}^{0}$ with the normalization of the NRQCD matrix element in Ref.nrqcd . Then the LO fragmentation function for the quarkonium can be obtained through replacing $\langle{\cal O}^{Q\bar{Q}[^{3}S_{1}^{[1]}]}(^{3}S_{1}^{[1]})\rangle$ by $\langle{\cal O}^{H}(^{3}S_{1}^{[1]})\rangle$. Under the leading non-relativistic approximation, the matrix element $\langle{\cal O}^{H}(^{3}S_{1}^{[1]})\rangle$ can be expressed by the radial wave function at the origin for the quarkonium $H$, i.e,

[TABLE]

Finally, the LO fragmentation function for the ${}^{3}S_{1}$ quarkonium state takes the form

[TABLE]

which is exactly the same as that of Ref.Braaten1 .

III The NLO correction to the fragmentation function

At the NLO level, we need to deal with the virtual and real corrections to the LO terms. It is hard to give the analytic expressions for those NLO-terms. In the following subsections, we shall give some explanations on how to deal with the virtual and real corrections.

In doing the calculations, the FeynCalc package feyncalc1 ; feyncalc2 is adopted to carry out the color and Dirac traces, the $Apart package apart and the FIRE package fire are used to do partial fraction and integration-by-parts (IBP) reduction. The master integrals are calculated by using the LoopTools package looptools . As a subtle point, there are some master integrals, which contain an eikonal propagator and can not be calculated by using the LoopTools, and we adopt the method introduced in Ref.Braaten2 to deal with those master integrals.

III.1 The virtual corrections

At the NLO level, the virtual corrections come from the cut diagrams containing a loop on either side of the cut. Four typical virtual corrections are shown in Fig.2. The differential phase space of the virtual corrections is the same as that of the LO one, e.g. Eq.(15). The virtual corrections can be obtained through

[TABLE]

where ${\cal A}_{\rm virtual}$ is the squared amplitude for the virtual corrections.

We adopt the method of regions region to calculate the fragmentation function. Within this method, we only need to calculate the contributions from the hard region, the Coulomb divergences which come from the potential region do not appear in the calculations. The hard region contributions are then obtained by expanding the relative velocity between the produced $Q$ and $\bar{Q}$ before the loop integration. Thus in the calculation, we only need to set the relative momentum to zero before the loop integration region .

III.2 The real corrections

At the NLO level, the real corrections come from the fragmentation process emitting an extra gluon, i.e, we need to deal with the process, $Q(K)\to Q\bar{Q}[^{3}S_{1}^{[1]}](p_{1})+Q(p_{2})+g(p_{3})$. The cut diagrams for the real corrections can be obtained from the LO cut diagrams by adding a gluon line crossing the cut. Four typical real correction cut diagrams are shown in Fig.3.

The differential phase space for the real correction is

[TABLE]

The real corrections can be obtained through

[TABLE]

where ${\cal A}_{\rm real}$ denotes the squared amplitude for the real corrections. As for the real corrections, the IR divergences come from the limits $p_{3}\to 0$ and $p_{3}\cdot n\to 0$, and the UV divergences come from the limit $|\textbf{p}_{3\perp}|\to\infty$. We adopt the method of Ref.Braaten2 to extract those divergences. Following this method, the real corrections can be represented as

[TABLE]

where ${\cal A}_{S}$ denotes the subtraction term, which has the same singularities as the squared amplitude of the real corrections. The first term on the right hand side of Eq.(26) is finite and can be directly calculated in $4$-dimensions. The integral of the subtraction term is divergent and should be calculated in $d$-dimensions.

The subtraction term can be constructed according to the singularity behavior of the squared amplitude of the real corrections. More explicitly, the squared amplitude for the real corrections can be written as

[TABLE]

where the Lorentz invariants are defined as follows

[TABLE]

The ${\cal A}^{\rm finite}_{\rm real}$ stands for the terms which are finite after the phase-space integration. There are neither $1/t_{1}^{2}$ term nor $1/t_{1}t_{2}$ term in ${\cal A}_{\rm real}$ because those terms are canceled after summing all the terms of the real corrections. The subtraction term ${\cal A}_{\rm S}$ can then be constructed as follows

[TABLE]

where $\tilde{s}$ is defined as

[TABLE]

with

[TABLE]

Numerically, we have found that the integration of $({\cal A}_{\rm real}-{\cal A}_{\rm S})$ over the phase space is finite in $4$-dimensions, which confirms our present choice of the subtraction term ${\cal A}_{\rm S}$. Due to its much simpler structure, following the method of Ref.Bcfragnlo , the phase-space integration over the subtraction term ${\cal A}_{\rm S}$ can be done analytically. The interesting reader can turn to Ref.Bcfragnlo for details.

III.3 Renormalization

The UV divergences in the virtual and real corrections should be canceled through renormalization. The counter-term approach is adopted to carry out the renormalization, where the fragmentation function is calculated with the renormalized quark mass $m_{Q}$, the renormalized field $\Psi_{r}$, the renormalized gluon field $A^{\mu}_{r}$ and the renormalized coupling constant $g_{s}$. The renormalized quantities are related to their corresponding bare quantities as

[TABLE]

where the renormalization constants $Z_{i}=1+\delta Z_{i}$, with $i=m,2,3,g$, respectively. The quantities $\delta Z_{i}$ are fixed by the renormalized conditions which define a renormalization scheme. The quark field, quark mass and gluon field are renormalized in the on-mass-shell scheme (OS), whereas the strong coupling constant $g_{s}$ is renormalized in the $\overline{\rm MS}$ scheme. The expressions of $\delta Z_{i}$ can be derived:

[TABLE]

where $\mu_{R}$ is the renormalization scale, $\beta_{0}=11C_{A}/3-4T_{F}n_{f}/3$ is the one-loop coefficient of the $\beta$-function in QCD, and $n_{f}$ is the number of active quark flavors. $\beta^{\prime}_{0}=11C_{A}/3-4T_{F}n_{lf}/3$, and $n_{lf}=3$ is the number of the light-quark flavors. For $SU_{c}(3)$ group, we have $C_{A}=3$, $C_{F}=4/3$ and $T_{F}=1/2$.

The operator products in the definition of the fragmentation functions also require renormalization Collins ; Mueller , whose counter-terms in $\overline{\rm MS}$ scheme can be written as Braaten2

[TABLE]

where $D_{Q\to Q\bar{Q}[^{3}S_{1}^{[1]}]}^{\rm LO}(z)$ denotes the LO fragmentation function in $d$-dimensional space time.

IV Numerical results

In doing the numerical calculations, the input parameters are taken as follows:

[TABLE]

where the values of $|R_{J/\Psi}(0)|^{2}$ and $|R_{\Upsilon}(0)|^{2}$ are taken from the potential-model calculations pot . For the 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 the $\beta$-function. According to $\alpha_{s}(m_{{}_{Z}})=0.1185$ pdg , we obtain $\Lambda_{\rm QCD}^{n_{f}=5}=0.233{\rm GeV}$ and $\Lambda_{\rm QCD}^{n_{f}=4}=0.337{\rm GeV}$. Then we have $\alpha_{s}(2m_{c})=0.259$, $\alpha_{s}(3m_{c})=0.223$, $\alpha_{s}(2m_{b})=0.180$, and $\alpha_{s}(3m_{b})=0.164$.

IV.1 The fragmentation functions

The fragmentation functions for $c\to J/\Psi$ and $b\to\Upsilon$ are shown in Figs.4 and 5, respectively. The initial factorization scale has been set as the minimal invariant mass of the initial $c$ quark or $b$ quark. Two typical values for the renormalization scale $\mu_{R}$ are adopted, one is the threshold energy to create a $c\bar{c}$ or $b\bar{b}$ pair, and the other is the minimal invariant mass of the initial $c$ quark or $b$ quark. Figs. 4 and 5 indicate that the renormalization scale dependence for the LO fragmentation functions are large, and the NLO corrections are very important to reduce the renormalization scale dependence.

For future applications, we use polynomials to fit the NLO fragmentation functions. We write the NLO fragmentation functions as the form

[TABLE]

where $D^{\rm LO}_{Q\to H}(z)$ is the LO fragmentation function given in Eq.(22). For $c\to J/\Psi$, we have

[TABLE]

For $b\to\Upsilon$,

[TABLE]

Using the fragmentation function, we can obtain two useful quantities, i.e. the fragmentation probability ($P$) and the averaged energy fraction ($\langle z\rangle$), which are defined as

[TABLE]

The fragmentation probabilities and the averaged energy fractions for $c\to J/\Psi$ and $b\to\Upsilon$ are shown in Tables 1 and 2, respectively.

The fragmentation functions shown in Figs. 4 and 5 are for $\mu_{F}=3m_{Q}$. The fragmentation functions at any other factorization scales can be obtained through DGLAP equation. To apply the NLO fragmentation functions to the production of heavy quarkonia at a $Z$ factory, we present the fragmentation functions for $\mu_{F}=m_{{}_{Z}}$ in Figs. 6 and 7, which are obtained by using the fragmentation functions at the initial value $\mu_{F}=3m_{Q}$. For definiteness, we set the renormalization scale $\mu_{R}=2m_{Q}$. In doing the numerical calculation, the DGLAP equation is solved through the Mellin transformation evolution1 ; evolution2 , and the NLO expression for $P_{QQ}$ is used as the evolution kernel. Figs. 6 and 7 show that the peaks of the fragmentation functions become lower and shift to a smaller $z$ value for a larger factorization scale, leading to a smaller value for the averaged energy fraction $\langle z\rangle$. For examples, we obtain $\langle z\rangle^{\rm NLO+NLL}|_{\mu_{F}=m_{Z}}=0.47$ for $D_{c\to J/\Psi}$ and $\langle z\rangle^{\rm NLO+NLL}|_{\mu_{F}=m_{Z}}=0.54$ for $D_{b\to\Upsilon}$.

IV.2 The $J/\Psi$ and $\Upsilon$ production at a $Z$ factory

As an application, we apply the NLO fragmentation functions of $c\to J/\Psi$ and $b\to\Upsilon$ to the production of $J/\Psi$ and $\Upsilon$ at a super $Z$ factory.

A Chinese group has raised the proposal of constructing a high-luminosity $e^{+}e^{-}$ collider in China, the so-called super $Z$ factory zfactory , which is similar to the Gigaz program suggested by the Internal Linear Collider Collaboration gigaz1 ; gigaz2 but with a even higher luminosity. As for this purposed super $Z$ factory, an $e^{+}e^{-}$ collider shall run at the energies around the mass of $Z^{0}$-boson resonance and with a high luminosity up to ${\cal L}=10^{34}-10^{36}{\rm cm}^{-2}{\rm s}^{-1}$. Due to the $Z^{0}$-boson resonance effect, large numbers of $J/\Psi$ and $\Upsilon$ events can be produced, thus providing a good platform studying the $J/\Psi$ and $\Upsilon$ properties.

In this case, the factorization formula for the production of $J/\Psi$ or $\Upsilon$ can be written as

[TABLE]

where the energy fraction $z$ is defined as $z\equiv 2p\cdot k/k^{2}$, $p$ is the momentum of the produced quarkonium, $k$ is the sum of the momenta of the initial electron and positron. Up to NLO level, the parton $i$ may be a heavy quark or a heavy antiquark. For the quarkonium production, the fragmentation function $D_{\bar{Q}\to H}$ is the same as $D_{Q\to H}$.

Due to that the coefficient functions $d\hat{\sigma}/dy$ are independent to the species of the produced hadron, they can be extracted by applying the pQCD factorization formula to the production of an on-shell heavy quark ($Q$) or heavy antiquark ($\bar{Q}$) pqcdfac1 ; pqcdfac2 . The expression for the coefficient function $d\hat{\sigma}/dy$ in $\overline{\rm MS}$-scheme up to NLO level has been given in Refs.coefun1 ; coefun2 , e.g.

[TABLE]

where $\sigma_{e^{+}e^{-}\to Q\bar{Q}}^{\rm LO}$ is the LO cross section for the $Q\bar{Q}$ production. And in the massless limit $m_{Q}\to 0$ 111Because the coefficient function $d\hat{\sigma}/dy$ is infrared safe, so we can take the limit $m_{Q}\to 0$ to do our calculation, which shall introduce a small error of ${\cal O}(m_{Q}^{2}/s)$., we have

[TABLE]

where $e_{f}$ is the electric charge of fermion $f$,

[TABLE]

are the vector and axial-vector couplings of fermion $f$ to the $Z$ boson, $T_{3f}$ is the third component of weak isospin of the fermion $f$, $\theta_{w}$ is the weak mixing angle, and the propagator functions are:

[TABLE]

The function $C(y)$ in the massless limit $m_{Q}\to 0$ takes the form

[TABLE]

For comparison, we adopt three strategies to calculate the differential cross sections $d\sigma/dz$ under the fragmentation approach. We denote them as “Frag, LO”, “Frag, NLO” and “Frag, NLO+NLL”, respectively. For the case of “Frag, LO”, the differential cross sections are calculated through

[TABLE]

where the factor $2$ comes from that the contribution of the $\bar{Q}$ fragmentation which is the same as that of $Q$ fragmentation. $D^{LO}_{Q\to H}(z)$ denotes the LO fragmentation function which has been given in Eq.(22). In the second equation, we have used the fact that $d\hat{\sigma}^{\rm LO}_{e^{+}e^{-}\to Q+\bar{Q}}/dy=\sigma^{\rm LO}_{e^{+}e^{-}\to Q+\bar{Q}}\delta(1-y)$. For the case of “Frag, NLO”, the differential cross sections are calculated through

[TABLE]

In the above calculation, the factorization and renormalization scales are set as $\mu_{F}=3m_{Q}$ and $\mu_{R}=2m_{Q}$. For the case of “Frag, NLO+NLL”, the differential cross sections are calculated through

[TABLE]

where the factorization scale and renormalization scale are set as $\mu_{F}=\mu_{R}=m_{{}_{Z}}$, and the fragmentation functions $D_{Q\to H}(z,\mu_{F}=m_{{}_{Z}})$ are obtained through solving the DGLAP evolution equation, i.e., Eq.(4) and the initial fragmentation functions $D^{\rm NLO}_{Q\to H}(z,\mu_{F}=3m_{Q})$ with $\mu_{R}=2m_{Q}$ are used as the boundary condition. This way, the large log-terms such as $\ln(m_{Q}^{2}/m_{Z}^{2})$ are resummed up to next-to-leading-log (NLL) accuracy.

We present the differential cross sections $d\sigma/dz$ for the production of $J/\Psi$ and $\Upsilon$ via the fragmentation approach in Figs. 8 and 9. In drawing the figures, we have omitted the $\gamma-\gamma$ and $\gamma-Z$ contributions, which are quite small compared with the dominant $Z-Z$ contribution around the $Z$ pole. Figs. 8 and 9 show how the NLO-terms and the leading and next-to-leading logarithms affect the predictions. For the $J/\Psi$ production, the NLO contribution is significant, and after including the NLO-terms, the distribution becomes softer and the value of $z$ corresponding to the peak value of the distribution becomes smaller. For the $\Upsilon$ production, the NLO contribution is relatively small compared with the $J/\Psi$ case.

Integrating the differential cross sections $d\sigma/dz$ over $z$, we can obtain the total cross sections for $J/\Psi$ and $\Upsilon$ production at the $Z$ factory, which can be simplified as

[TABLE]

where $P_{Q\to H}$ is the fragmentation probability for $Q$ into the quarkonium $H$. The results are presented in Table 3, which shows that the NLO corrections enhance the total cross section $\sim 7\%$ for $J/\Psi$ and $\sim 4\%$ for $\Upsilon$.

IV.3 A comparison of $J/\Psi$ and $\Upsilon$ production via the $Z$ decays up to NLO level

As a final remark, we compare the NLO results under the fragmentation function approach and the complete pQCD approach. The complete NLO pQCD calculations of $Z\to J/\Psi+X$ and $Z\to\Upsilon+X$ have been done by Ref.JXWang . We present the differential decay widths $d\Gamma/dz$ for $Z\to J/\Psi+X$ and $Z\to\Upsilon+X$ in Figs. 10 and 11, in which all the input parameters are taken as those of Ref.JXWang and the results denoted by “Li and Wang” are results under the complete NLO calculation. The results denoted by “Combined, NLO+NLL” are results combining the results from the complete pQCD approach which are taken from Ref.JXWang and the results from the fragmentation approach, i.e, $d\Gamma/dz|_{\rm Combined,NLO+NLL}=d\Gamma/dz|_{\rm Li\;and\;Wang}+(d\Gamma/dz|_{\rm Frag,NLO+NLL}-d\Gamma/dz|_{\rm Frag,NLO})$. For the case of $J/\Psi$ production at the $Z$ factory, Fig.10 shows the fragmentation contributions dominant the decay width, since the “Frag, NLO” shape is very close to the complete NLO one. Moreover, the total decay widths are $\Gamma(Z\to J/\Psi+X)|_{\rm Frag,NLO}=136\,{\rm keV}$, $\Gamma(Z\to J/\Psi+X)|_{\rm Frag,NLO+NLL}=130\,{\rm keV}$, $\Gamma(Z\to J/\Psi+X)|_{\rm Li\;and\;Wang}=136\,{\rm keV}$ and $\Gamma(Z\to J/\Psi+X)|_{\rm Combined,NLO+NLL}=130\,{\rm keV}$. Fig.11 shows that the “Frag, NLO” result is larger than the result under the complete NLO calculation at small $z$-region, thus the fragmentation approximation for the $\Upsilon$ production is not as good as the $J/\Psi$ case. This is because the $b$-quark mass is larger than $c$-quark mass, and the power correction in $m_{b}^{2}/s$ for the case of $\Upsilon$ is larger than the case of $J/\Psi$. However, the combined result counts both the large power correction and the large logarithms, then give a good prediction. Moreover, the total decay widths are $\Gamma(Z\to\Upsilon+X)|_{\rm Frag,NLO}=20.9\,{\rm keV}$, $\Gamma(Z\to J/\Upsilon+X)|_{\rm Frag,NLO+NLL}=20.6\,{\rm keV}$, $\Gamma(Z\to J/\Upsilon+X)|_{\rm Li\;and\;Wang}=17.38\,{\rm keV}$, and $\Gamma(Z\to J/\Upsilon+X)|_{\rm Combined,NLO+NLL}=17.08\,{\rm keV}$.

V Summary

In the present paper, we have calculated the fragmentation function for a heavy quark into heavy quarkonium, e.g. $c\to J/\Psi$ or $b\to\Upsilon$, up to NLO level. Our present results are complementary to the previous works on the fragmentation function of a gluon into heavy quarkonia done in the literature, which is pQCD calculable due to the fact that the gluon should be hard enough to form a heavy quark-and-antiquark pair.

Our results show that the NLO correction is important to suppress the renormalization scale uncertainty and to achieve a reliable fragmentation prediction. Our calculations are based on the gauge-invariant definition of the fragmentation function suggested by Collins and Soper. To avoid large logarithms appearing in the perturbative series of the fragmentation function, we first derive the fragmentation function at an initial (reasonable) factorization scale $\mu_{F}=3m_{Q}$, and then run to any factorization scale with the help of the DGLAP evolution equation. This treatment, in effect, resums the large logarithms and forms a reliable prediction. Thus for the cases when the fragmentation dominants the quarkonium productions or decays, our present calculated fragmentation functions shall be of great help for a more precise pQCD prediction. As an application, we have applied the obtained fragmentation functions to the production of $J/\Psi$ and $\Upsilon$ at the super $Z$ factory. The shape of the $J/\Psi$ distribution changes significantly by introducing the NLO corrections, and the total cross section increases by $\sim 7\%$. The shape of the $\Upsilon$ distribution changes slightly by introducing the NLO corrections, and the total cross section increases only $\sim 4\%$.

Acknowledgments: This work was supported in part by the Natural Science Foundation of China under Grant No.11625520, No.11847222, No.11847301, No.11675239, No.11535002, and by the Fundamental Research Funds for the Central Universities under Grant No.2019CDJDWL0005.