NLO fragmentation functions for a quark into a spin-singlet quarkonium: Same flavor case

In the paper, we calculate the fragmentation functions for $c \to \eta_c$ and $b \to \eta_b$ up to next-to-leading-order (NLO) QCD accuracy. The ultraviolet divergences in the real corrections are removed through operator renormalization under the modified minimal subtraction scheme. We then obtain the fragmentation functions $D_{c \to \eta_c}(z,\mu_F)$ and $D_{b \to \eta_b}(z,\mu_F)$ up to NLO QCD accuracy, which are presented as figures and fitting functions. The numerical results show that the NLO corrections are significant. The sensitives of the fragmentation functions to the renormalization scale and the factorization scale are analyzed explicitly.


Introduction
Since the discovery of the J/ψ in 1974, the heavy quarkonium production has been a focus of theoretical and experimental interest. It is because the production of a heavy quarkonium involves both perturbative and nonperturbative aspects of QCD, and it provides a good platform to study QCD. The most successful effective theory to describe the quarkonium production is the nonrelativistic QCD (NRQCD) effective theory [1]. Many important quarkonium production processes have been studied up to next-to-leading order (NLO) accuracy under the NRQCD factorization at various colliders [2,3]. However, there are still challenges in understanding the quarkonium production within the NRQCD factorization, such as the J/ψ polarization puzzle [4][5][6] and the large differences among various sets of long-distance matrix elements (LDMEs) extracted by several groups, c.f. refs. [5][6][7][8][9]. Therefore, it is important to further study the quarkonium production mechanism.
The production of a quarkonium at high transverse momentum p T region is simpler than other cases, because the long-distance interactions between the produced quarkonium and initial particles are suppressed. Therefore, in order to explore the quarkonium production mechanism, it is important to study the quarkonium production at high transverse momentum (p T ) region.
According to QCD factorization theorem, the production cross section of a hadron (H) at high p T region is dominated by the single parton fragmentation [10], i,e., where the symbol ⊗ represents a convolution integral over the momentum fraction z, dσ A+B→i+X (p T /z, µ F ) are partonic cross sections, D i→H (z, µ F ) are fragmentation functions for a parton into a hadron H. The sum extends over all species of partons. µ F is the factorization scale which separates the energy scales of two parts. The factorization formula (1.1) is also called the leading-power (LP) factorization, since it gives the LP contribution in the expansion in powers of m H /p T . For the quarkonium production, the next-to-leading power (NLP) contribution can be factorized to double-parton fragmentation [11][12][13][14]. Unlike the fragmentation functions for light hadrons, the fragmentation functions for quarkonia contain perturbatively calculable information. In fact, the fragmentation functions for the production of quarkonia can be refactorization through the NRQCD factorization, i.e.,  , and a few fragmentation functions for quarkonia have been calculated up to order α 3 s [36][37][38][39][40][41][42][43][44][45]. Among these studies, the fragmentation functions for q → η Q (q = Q), which are of order α 3 s , have been calculated in our recent paper [45]. However, the fragmentation functions for c → η c and b → η b are only available up to order α 2 s . For the precision prediction of the production rate of the η c,b at high-energy colliders such as LHC and etc., it is also important to know the fragmentation functions of c → η c and b → η b up to order α 3 s . In this paper, we will calculate the fragmentation functions D c→ηc (z, µ F ) and D b→η b (z, µ F ) up to NLO QCD accuracy.
The paper is organized as follows. In Sec.II, we present the definition of fragmentation function and the calculation for the LO fragmentation functions. In Sec.III, we sketch the method used in the calculation of the NLO corrections to the fragmentation functions. In Sec.IV, we present the numerical results for the fragmentation functions D c→ηc (z, µ F ) and D b→η b (z, µ F ). Section V is reserved for a summary.

LO fragmentation function
Before calculating the fragmentation function D Q→η Q (z, µ F ), we first present the definition of the fragmentation function for a quark into a hadron. To give the definition for the fragmentation function, it is convenient to work in the light-cone coordinate system. In this coordinate system, a d-dimensional vector is expressed as We adopt a gauge-invariant definition of fragmentation function which was introduced by Collins and Soper in ref. [46].
For a quark into a hadron H, the fragmentation function is defined as where Ψ is the initial quark field, A µ a is the gluon field, t a is the color matrix, and P indicates the path ordering. The longitudinal momentum fraction is defined as z ≡ P + /K + , and K is the momentum of the initial quark. This definition of the fragmentation function is carried out in a reference frame where the transverse momentum of the produced hadron H vanishes. It is convenient to introduce a lightlike vector, and in the reference frame, we have n µ = (0, 1, 0 ⊥ ). Then, the "+" component of a vector can be expressed as V + = V · n, and z can be expressed as a Lorentz invariant z = P · n/K · n. According to the definition, the Feynman rules for the fragmentation function can be derived, and the Feynman rules can be found in our previous paper [42].
To obtain the fragmentation function for Q → η Q , we first calculate the fragmentation function for an on-shell QQ pair in 1 S [1] 0 state. Then the fragmentation function D Q→η Q can be obtained from D Q→(QQ)[ 1 S [1] 0 ] through replacing the LDME (4) (3) 0 ](p 1 ) + Q(p 2 ), as shown in Fig.1. The squared amplitudes, corresponding to four diagrams in Fig.1, can be written according to the Feynman rules, i.e,

4)
where Π 1 is the spin projector for the 1 S 0 state where 1 is the unit matrix for the SU c (3) group. Carrying out the color and the Dirac traces, we obtain the expression of the total squared amplitude ( where C F = 4/3, s 1 = (p 1 + p 2 ) 2 , and The differential phase space for the fragmentation function at LO is where the δ function is due to the cut through the eikonal line. The integration over p + 2 can be performed directly using the δ function. The squared amplitude does not depend on the angles of p 2⊥ , thus the integration over the angles of p 2⊥ is trivial, and can be performed easily. Then the differential phase space reduces to where the relation (2.14) The fragmentation function D Q→η Q can be obtained from which is the same as that obtained in ref. [17].

NLO corrections to the fragmentation function
At NLO, there are virtual and real corrections contributing to the fragmentation function.
In this section, we will briefly introduce the methods used to calculate the virtual and the real corrections.
In the calculation, the package FeynCalc [47,48] is used to carry out the Dirac and color traces, the packages $Apart [49] and Fire [50] are used to do partial fraction and integration-by-part (IBP) reduction. After the IBP reduction, all the one-loop integrals are reduced to master integrals which include the common A 0 , B 0 and C 0 functions, and the scalar one-loop integrals with one eikonal propagator. The common A 0 , B 0 and C 0 functions are calculated using LoopTools [51] numerically, while the scalar integrals with one eikonal propagator are calculated using the method which was introduced in ref. [36].

Virtual corrections
The virtual corrections come from the cut diagrams which contain a loop on either side of the cut line. Six sample cut diagrams for the virtual corrections are presented in Fig.2. The virtual corrections can be calculated through where A virtual denote the squared amplitude for the virtual corrections. (1) (3)  There are ultraviolet (UV) and infrared (IR) divergences in the NLO calculations, dimensional regularization is adopted to regularize these divergences. In dimensional regularization, γ 5 should be noted. We adopt a practical prescription, which was introduced in ref. [52], to handle γ 5 in dimensional regularization. There are Coulomb divergences in the virtual corrections. Conventionally, these Coulomb divergences are regularized by the relative velocity of the produced (QQ) pair, and they should be absorbed into the LDMEs. In this paper, we adopt the threshold expansion method [53] to extract the NRQCD SDCs. In this method, we expand the relative momentum q of the produced (QQ) pair before carrying out the loop integration. In the leading nonrelativistic approximation, we just set q = 0 before the loop integration. Then, the Coulomb divergences, which are power divergences, vanish in the calculation. (1) (2) The real corrections come from the fragmentation process Q(K) → (QQ)[ 1 S [1] 0 ](p 1 ) + Q(p 2 ) + g(p 3 ). Six sample cut diagrams are given in Fig.3. The real corrections can be obtained through

Real corrections
where A real denotes the squared amplitude for the real corrections, and the differential phase space for the real corrections is There are UV and IR divergences in the real corrections. These divergences come from the phase space integrals and should be regularized by dimensional regularization as that in the virtual corrections. In order to isolate the divergent and the finite terms, we adopt the subtraction method. Under this method, the real corrections are expressed as where A S denotes the subtraction term, which has the same singularities as A real . The first term on the right-hand side of Eq.(3.4) is finite, thus it can be calculated in 4 dimensions directly. The second term on the right-hand side of Eq. (3.4) is divergent, and it should be calculated in d dimensions analytically. The methods of constructing and integrating the subtraction terms can be found in our previous paper [42]. It should be noted that there are new types of cut-diagrams (e.g, the fifth and the sixth diagrams in Fig.3) in the η Q case compared to the J/ψ case. There are additional subtraction terms arise from these new-type cut diagrams, and the integration of these additional subtraction terms can be found in our another paper [45]. With those formulas presented in refs. [42,45], the real corrections can be calculated directly.

Renormalization
The calculation is carried out with the renormalized fields Ψ R and A µ R , the renormalized quark mass m Q , and the renormalized strong coupling constant g s . The relations between the renormalized and the bare quantities are where Z i ≡ 1 + δ i are renormalization constants. The renormalization constants are fixed by the renormalization scheme. In this paper, the renormalization of the strong coupling constant is carried out in the modified minimal-subtraction scheme (MS), while the renormalization of the quark field, the gluon field and the quark mass m Q are carried out in the on-mass-shell scheme (OS). The expressions of the quantities δ i are where C F = 4/3, T F = 1/2, β 0 = 11 − 2 n f /3, β 0 = 11 − 2 n lf /3, n f is the number of the active-quark flavors and n lf is the number of the light-quark flavors. The contributions from these counterterms can be obtained through where A counter denotes the squared amplitude for the counterterms from the strong coupling constant, the quark field, the gluon field and the quark mass.
The operator used to define the fragmentation function is also need to be renormalized [54]. We carry out the operator renormalization in the MS scheme, and and P gQ (y) = C F 1 + (1 − y) 2 y . (3.10)

Numerical results and discussion
Summing the contributions from the virtual and real corrections and the counter terms, the UV and IR divergences are canceled, and the finite fragmentation function up to NLO QCD accuracy is obtained. The fragmentation function for the η Q can be obtained from the fragmentation function for the (QQ)[ 1 S [1] 0 ] pair by multiplying a factor O η Q ( 1 S In doing the numerical calculation, the input parameters are taken as follows: S (0) are the radial wave functions at the origin for the (cc) and (bb) systems, and the input values for them are taken from the potential-model calculation [55]. For the strong coupling constant, we adopt the two-loop running formula, i.e., , where β 1 = 102 − 38 n f /3 is the two-loop coefficient of the QCD β-function. According to α s (m Z ) = 0.1185 [56], we obtain Λ in the LO fragmentation process; the factorization scale is set as µ F = 3m c , i.e., the minimal invariant mass of the initial c quark. From the figure, we can see that the effect of the NLO corrections is significant. The fragmentation function is decreased at small z values and increased at large z values after including the NLO corrections. The NLO fragmentation function behaves like 1/z as z → 0, while the LO fragmentation function tends to 0 as z → 0. This is because that there are new type cut diagrams (e.g., the sixth diagram in Fig.3) contributing to the real correction, and these new type cut diagrams are the same as those of a light quark into the η c . As shown in our previous paper [45], the fragmentation function D q→ηc (z, µ F ) behaves like 1/z as z → 0, then the 1/z behavior of the NLO fragmentation function for c → η c is expected. The sensitivity of the LO and NLO fragmentation functions to the renormalization scale µ R is also shown in Fig.4. The bands in Fig.4 are obtained by varying the renormalization scale µ R by a factor of 2 around the center value 2m c . From the figure, we can see that the sensitivity of the NLO fragmentation function to µ R is decreased at the moderate and the large z values, but increased at the small z values, compared to the LO fragmentation function. The reason for the large sensitivity of the NLO fragmentation function to µ R at small z values is that the NLO correction is large compare with the LO contribution at the small z values.
The sensitivity of the LO and NLO fragmentation functions to the factorization scale µ F is shown in Fig.5. The band in Fig.5 is obtained by varying the factorization scale µ F by a factor of 2 around the center value 3m c . From the figure, we can see that the LO fragmentation function is independent of µ F , the µ F dependence of the fragmentation function starts at the NLO. The NLO fragmentation function is very sensitive to µ F when z is very small, which is similar to the fragmentation function of q → η c [45]. It is noted that the dependence of the fragmentation function on µ F is not a theoretical uncertainty, because the fragmentation function is always defined at a specified scale µ F . The fragmentation function D c→ηc (z, µ F ) contains logarithms of µ F /m Q , and these logarithms may spoil the the convergence of the perturbative expansion of the fragmentation function when µ F m c . These large logarithms can be resummed through solving the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [57][58][59], where the fragmentation function with µ F ∼ m c is used as the boundary condition.
The LO and NLO fragmentation functions for b → η b are shown in Fig.6, the sensitivities of the LO and NLO fragmentation functions to µ R and µ F are shown in Fig.6 and Fig.7 respectively. The sensitivities of the fragmentation functions for b → η b to µ R and µ F are similar to the c → η c case.
Since the NLO fragmentation functions for c → η c and b → η b behave like 1/z as z → 0, the fragmentation probabilities for c → η c and b → η b at the NLO level are infinite. However, the physical cross sections shall be finite, because the phase space limitations present a lower bound for z.
To see more about the effect of the NLO corrections, we calculate the moments of the fragmentation functions. The moments of the fragmentation functions D Q→η Q (z, µ F ) can be defined as those moments are moderate although the K-factors of the fragmentation probabilities are divergent.
For future applications, we give fitting functions to the NLO fragmentation functions. The NLO fragmentation functions can be written as following form  The difference between the f (z) functions for c → η c and b → η b arises from the heavy quark loop in the gluon vacuum polarization. We only consider one heavy flavor (i.e., c) contributing to the gluon vacuum polarization for c → η c , while we consider two heavy flavors (i.e., c and b) contributing to the gluon vacuum polarization for b → η b .

Summary
In the present paper, we have calculated the fragmentation functions for c → η c and b → η b up to NLO QCD accuracy. The results obtained in this paper are complementary to the previous works on the fragmentation functions for q → η Q and g → η Q at order α 3 s . The most difficult part in this work is the calculation of the real corrections. We adopt the subtraction method to calculate the real corrections. The construction of the subtraction terms and the parametrization of the phase space have been developed in our previous works.
The fragmentation functions D c→ηc (z, µ F ) and D b→η b (z, µ F ) with µ F = 3m Q and µ R = 2m Q under the MS factorization scheme are presented in the forms of figure and fitting function. The results show that the effect of the NLO corrections is significant. The fragmentation functions are decreased at small z values and increased at large z values after including the NLO corrections. Moreover, the NLO fragmentation functions have a singularity at z = 0, while the LO fragmentation functions are zero at z = 0. The sensitives of these fragmentation functions to the renormalization scale µ R and the factorization scale µ F are analyzed explicitly. The NLO fragmentation functions obtained in this paper can be applied to the precision predictions of the η c and η b production at high-energy colliders.