Relativistic correction to gluon fragmentation function into pseudoscalar quarkonium

Inspired by the recent measurements of the $\eta_c$ meson production at LHC, we investigate the relativistic correction effect for the fragmentation function of the gluon into $\eta_c$, which constitutes the crucial nonperturbative elements to understand $\eta_c$ production at high $p_T$. Employing three distinct methods, we calculate the leading relativistic correction to the $g\to\eta_c$ fragmentation function in the NRQCD factorization framework, as well as verify the existing NLO result for the $c\to \eta_c$ fragmentation function. We also study the evolution behavior of these fragmentation functions with the aid of DGLAP equation.


Introduction
Heavy quarkonium production and polarization in various collider experiments has long been a topic of great interest in QCD, and has triggered intensive experimental and theoretical investigation in the past few decades (for a recent review, see [1]).
Thus far, the modern theoretical method to tackle heavy quarkonium (exemplified by J/ψ(ψ ) and Υ) production and decay is represented by the effective-fieldtheory approach dubbed nonrelativistic QCD (NRQCD) factorization [2]. For the production of a charmonium state H from colliding beams composed of particles of type A and B, NRQCD factorization allows one to, schematically, express the corresponding production rate as dσ[A+B → H+X] = n dσ n [A+B → cc(n)+X] O H n , (1) where n signifies the color/angular-momentum quantum number of the cc pair produced in the hard scattering. In (1), dσ n are the perturbatively-calculable short-distance coefficients, and O H n represent the nonperturbative, yet universal, vacuum expectation values of the NRQCD production operators that are sensitive to H and n. The power series in (1) is governed by the expansion in the characteristic velocity v of the c(c) quark inside a charmonium (the nonrelativistic nature of quarkonium implies that v 1), since each nonperturbative NRQCD production matrix element has definite power counting in v.
As an important step in theoretical progress in the past decade, various short-distance coefficients dσ n relevant to J/ψ production in virtually all the commissioning collider programs have been gradually made available to next-to-leading order (NLO) accuracy in the strong coupling constant, for both color-singlet and octet channels [1]. By comparing these NLO-accuracy NRQCD predictions with the various measurements conducted at B factorties, HERA, the Tevatron and the LHC, there is satisfactory agreement in some cases, but also alarming discrepancies in other cases, notoriously for J/ψ polarization in hadroproduction [1]. Unfortunately, current computational techniques limit our capability to further address the next-to-next-to-leading order (NNLO) perturbative corrections, so the viability of the NRQCD approach still awaits a sharper and more critical examination.
Very recently, for the first time, the LHCb collaboration measured the differential production rate of the pseudoscalar charmonium state η c in the range p T (η c ) > 6.5 GeV, tagged via the decay channel η c → pp [3]. This is an important supplement to our knowledge on charmonium production, since η c production is an even more ideal test-bed for NRQCD than J/ψ, owing to its simplicity as a spin-zero meson. Very recently, the NLO perturbative corrections to η c hadroproduction have been investigated in the NRQCD factorization approach [4][5][6].
Besides NRQCD factorization, there is another famous first-principle approach to tackle inclusive single hadron production, the so-called perturbative QCD (collinear) factorization, whose applicability is not confined merely to heavy quarkonium.
According to the collinear factorization theorem [7], at sufficiently high p T , the inclusive production rate of a specific hadron H is dominated by the following fragmentation mechanism: where i stands for a QCD parton (quark or gluon), and z is the light-cone momentum fraction carried by H with respect to the parent parton. dσ[A + B → i + X] is the perturbatively-calculable partonic hard cross section, and D i→H is the nonpertubative yet universal fragmentation function, characterizing the probability distribution for the parton i to hadronize into H carrying momentum fraction z. ⊗ indicates that the hard partonic cross section ought to be convoluted with the corresponding fragmentation function over z.
The fragmentation functions, such as a gluon fragmenting into π and p, are genuinely nonperturbative objects, which so far can only be extracted from experiments [8]. In contrast, the situation becomes greatly simplified if H is a heavy quarkonium. In this case, the fragmentation function D i→H (z, µ) contains several distinct energy scales: heavy quark mass m, typical threemomentum of quark mv, and even smaller scales such as mv 2 and Λ QCD . Owing to the fact m Λ QCD , it is conceivable that the hard scale m should be explicitly factored out from D i→H (z, µ). As a matter of fact, by demanding the equivalence of two factorization theorems (1) and (2), one concludes that the fragmentation function itself must be subject to the following NRQCD factorization theorem: Here d (n) (z) are the perturbatively calculable coefficient functions, and O H n are the same NRQCD matrix elements as appear in (1).
In passing, it is worth noting that, for heavy quarkonium production, the NLO power correction (the order-1/p 2 T contribution in (2)) has also recently been systematically developed. As a consequence, a new set of nonperturbative functions, dubbed double-parton fragmentation functions, must be introduced [9,10]. Analogous to (3), they are also subject to a similar NRQCD factorization procedure, with various LO short-distance coefficient functions having been recently calculated [11,12].
The physical picture underlying (3) was first elucidated and pursued in the NRQCD context by Braaten and collaborators in the early 1990s (note there was earlier work along this direction in the pre-NRQCD era [13]). In those studies, various quarkonium fragmentation functions were computed to lowest order in both α s and v, e.g. gluon/charm quark fragmentation into Swave quarkonium has been computed [14,15]. Recently, the relativistic corrections have been investigated for the g → J/ψ fragmentation function [16,17], as well as for the c → J/ψ, η c fragmentation functions [18]. Very recently, the order-α s correction has also been addressed to the g → η c fragmentation function [19].
The aim of the present work is to fill a missing gap, i.e., to compute the leading relativistic correction to the g → η c fragmentation function. This piece of knowledge, in supplement with the recently available radiative correction [19], might be helpful to interpret the recent LHC measurements on η c production at high p T .
From a theoretical perspective, there are many equivalent ways to calculate the quarkonium fragmentation functions. Originally, Braaten and collaborators invented a trick to directly extract the fragmentation functions in a process-independent fashion. This method is simple and efficient for a LO calculation, but may become cumbersome if one proceeds to higher order in α s and v. Soon after, Ma pointed out [15] that the NRQCD factorization of a quarkonium fragmentation function can be conveniently calculated starting from the operator definition of fragmentation function introduced by Collins and Soper [20]. This elegant approach has the advantage that it preserves manifest gauge invariance, and allows one to systematically address the higher-order corrections, as was illustrated in Refs. [16,19].
In this work, we will use three different approaches to compute the order-v 2 correction to the g → η c fragmentation function: the Collins-Soper definition, the Braaten-Yuan method, and extracting from a specific physical 023103-2 process involving η c production. After obtaining the desired expressions, we also attempt to study the evolution behavior of the fragmentation functions of g → η c and c → η c .
The rest of this paper is organized as follows. In Section 2, we present a short review of the factorization of the fragmentation function for g → η c in NRQCD framework, and briefly outline our matching strategy. In Section 3, we compute the short-distance coefficients for the gluon fragmenting into the spin-singlet state η c through relative order-v 2 by employing three different methods, and confirm that all of them yield the same answer. In Section 4, we revisit the c quark fragmenting into η c and verify the previous results through relative order-v 2 . We also employ the DGLAP equation to evolve both the fragmentation functions of g → η c and c → η c to higher energy scales. Finally we summarize in Section 5.

Fragmentation function in NRQCD
In line with the NRQCD factorization (3), through the relative order-v 2 , the fragmentation function of a gluon fragmenting into the pseudoscalar quarkonium η c reads and P ηc 1 are colorsinglet NRQCD production operators: where ψ, χ are Pauli spinor fields in NRQCD, and The involved nonperturbative matrix elements in (4) are the vacuum expectation values of those color-singlet NRQCD production operators as specified in (5b). Under the vacuum saturation approximation, which is accurate up to O(v 4 ), the LO matrix element can be approximated well by the Schrödinger wave functions at the origin for the η c in the potential model: where N c = 3 is the number of colors in QCD. Rather than cope with the order-v 2 matrix element itself, it is more convenient to introduce a dimensionless ratio of the following NRQCD matrix elements: where the second equality is again obtained by invoking the vacuum saturation approximation. It is often useful to estimate the v 2 ηc from the so-called Gremm-Kapustin relation [22]: If the charm quark mass is taken as the one-loop pole mass, m = 1.4 GeV, then v 2 ηc ≈ 0.13. Our central goal is to determine the short-distance coefficient functions d (0) (z) and d (2) (z). To this purpose, we will use the standard matching technique. Namely, since the short-distance coefficients are independent of the long-distance dynamics, we can freely replace the physical hadron η c by a free c(p)c(p) pair carrying the following momenta: where P 2 = 4E 2 , P · q = 0 and q 2 = m 2 − E 2 , so that p 2 =p 2 = m 2 . In the cc pair (quarkonium) rest frame, We can further enforce the cc pair to bear quantum number We then substitute this fictitious η c state into the factorization formula (4): Since both the left-hand side and right-hand side in (10) can now be computed in perturbation theory, we can readily solve for d (0) (z) and d (2) Firstly, one can trivially deduce the NRQCD matrix elements that appear in (10): Note that the c(c) state in the NRQCD matrix elements obeys the nonrelativistic normalization. We briefly describe our strategy of computing the left-hand side of (10). After writing down the QCD amplitude to produce the free c(p) andc(p),ū(p)Av(p), we have to project out the c(p)c(p) pair onto the desired 1 S (1) 0 state. We employ the standard covariant trace technique [23], with the aid of the following projector: Thereby we extract the singlet amplitude by the operation M =ū(p)Av(p) → Tr(AΠ (1) 1 ), where the cc pair is  in the spin/color-singlet state. We emphasize that, by using the above projection operator in (12), we have tacitly assumed that the quark and antiquark in the QCD side are normalized nonrelativistically. We need to further single out the S-wave orbital angular momentum contribution. We first truncate the amplitude M to quadratic order in q µ , then take the following procedure: where The M 0 and M 2 are then the desired QCD amplitudes to produce a cc pair in the 1 S (1) 0 state, accurate through order-v 2 .

Perturbative calculation of shortdistance coefficients
In this section, we will employ several different methods to ascertain the short-distance coefficients associated with the gluon fragmenting into η c through order v 2 . They all yield identical results, thus serving as a useful consistency check.

From Collins-Soper definition
The rigorous operator definition of the fragmentation function was introduced by Collins and Soper in 1981 [20]. It is convenient to adopt the light-cone coordinate, where x µ = (x + , x − , x ⊥ ) with x ± = (x 0 ± x 3 )/ √ 2, and x ⊥ = (0, x 1 , x 2 , 0). We assume the parent (virtual) gluon moves along the z axis, that is, k µ = (k + , k − , k ⊥ = 0), and the final-state hadron H carries the momen- where M H is the hadron mass. The corresponding g → H fragmentation function is defined as where d = 4 is the spacetime dimension, and G +µ a signifies the gluon field strength tensor. Specifically, the daughter hadron H carries the 4-momentum sents the gauge link: where A µ = A µ a T a is the matrix-valued gluon field, and T a is the generator of the SU (N c ) group in adjoint representation. P implies the path-ordering. The null vector n µ = (0, 1, 0 ⊥ ) defines the "minus" light-cone direction, so that n · A = A + . Following the above definition, together with the matching strategy outlined in Section 2, we can replace the physical η c by a fictitious free cc pair. For the gluon fragmenting into the cc( 1 S (1) 0 ) state, in the Feynman gauge one can draw four Feynman diagrams at the leading order in α s , one of which is depicted in Fig. 1. The relevant Feynman rules can be found in Ref. [20]. After some algebra, this perturbatively calculable fragmentation function can be expressed as where F c = Tr(T a T b )Tr(T a T b ) = N 2 c − 1 4 is the corresponding color factor, and l µ stands for the 4-momentum of the gluon recoiling against the cc pair. The rank-4 tensor G ρσαβ = g αβ g µν (g µρ k + −k µ n ρ )(g νσ k + −k ν n σ ) stems from the product of two vertices of the gluon interacting with the eikonal line, while L i (i = 1, 2) represents the remaining ordinary quark-gluon amplitudes: where the projector Π (1) 1 has been given in (12). The S-wave amplitudes are extracted following the recipe given in (13) and (14). After squaring the S-wave amplitudes, carrying out the phase-space integration to get rid of all the δ-functions in (17), and finally performing the trivial angular integration in P ⊥ , we end up with the following one-fold integral: where ρ = P 2 ⊥ is the modulus of the transverse momentum of the cc pair, and and The integral over ρ can be transformed into a more conventional phase-space integral through where s is the squared invariant mass of the final-state cc + g system. ρ 0 then implies that s 4E 2 /z. If we replace ρ by s in (19), and replace the lower boundary to s 4m 2 /z, we can fully recover the corresponding LO expression in Ref. [19]. After fulfilling the integration over ρ, and matching (10) onto (19), we can obtain the corresponding shortdistance coefficient functions accurate through order-v 2 : Equation (22a) recovers the well-known LO result [14,15]. Equation (22b) is the central result of this work. The relativistic correction tends to dilute the LO fragmentation contribution. A curious feature is that, although the integrand in (21) does not resemble (20) at all, d (2) (z) turns out to bear the exactly identical functional dependence on z as d (0) (z). It is interesting to observe, although likely to be only a coincidence, that for the short-distance coefficient function associated with g → cc( 3 S (8) 1 ), the ratio of the order-v 2 term and the order-v 0 term also turns out to be − 11 6 [16]. We have also redone the calculation in the light-cone gauge A + = 0, where the gauge link in Fig. 1 is absent. We again reproduce the results listed in (22). Hence, we have explicitly checked the gauge invariance of this fragmentation function according to the Collins-Soper definition.
Substituting these short-distance coefficients (22) into the NRQCD factorization formula (4), we then obtain the fragmentation function of g → η c . It is interesting to deduce the total fragmentation probability of g → η c (the 1st Mellin moment of the fragmentation function): Obviously, provided that v 2 ηc is positive, incorporating the relativistic correction decreases the fragmentation probability of g → η c .
For the latter use, we are also interested in the 2nd Mellin moment of the fragmentation function, which might be interpreted as the average momentum fraction of the η c meson in the g fragmentation process:

Extraction from Higgs boson decay
In this section, we attempt to extract the g → η c fragmentation function from a specific process, say, inclusive η c production from Higgs boson decay, h → g * g → η c +gg.
The Higgs coupling to two gluons plays a crucial role in discovering the Higgs boson at the LHC experiment. In the Standard Model, it can be represented by an effective operator − λ v hG a µν G a,µν , where v signifies the Higgs vacuum expectation value, and the effective coupling λ = α s 12π + O(α 2 s ) receives the major contribution from the top quark loop. From this effective operator, one can readily deduce the Higgs boson hadronic width We are interested in inferring the energy spectrum of the η c meson in h(K) → η c (P ) + g(k 1 )g(k 2 ). It is convenient to introduce the three dimensionless energy fraction variables: which is subject to the energy conservation condition . We then expect that the gluon fragmentation function can be read off from while holding M ηc fixed.
In accordance with the matching ansatz, our goal is again to compute the process h → cc(P, 1 S (1) 0 ) + gg. At the lowest order in α s , there are four Feynman diagrams, one of which is depicted in Fig. 2. After truncating the S-wave amplitude through order q 2 , we then obtain the squared amplitude according to . Accordingly, the decay rate of h → cc( 1 S (1) 0 ) + gg can be put in the form: where we have replaced r by 4E 2 /M 2 h . The integration boundaries of z are explicitly labeled. Obviously, its allowed range reduces to 0 z 1 as M h → ∞. The boundaries for the energy fraction of gluon 1, dubbed After carrying out the phase-space integration over x 1 and substituting (28), we can reshuffle the corresponding expression to the second order in q. In line with (26), dividing this expression by the hadronic decay width of h → gg, we only retain those terms that survive in the r → 0 limit. By solving the matching equation (10), we find the exactly identical expressions of d (0) (z) and d (2) (z) as given in (22), which were previous obtained via the Collins-Soper definition.
Alternatively, one can also apply the Braaten-Yuan trick [14] to extract the gluon fragmentation function, through order v 2 . In a similar spirit to the preceding Higgs boson decay example, this approach also aims to extract the fragmentation function from a physical process. Nevertheless, the virtue of this method is that it does not need specify any concrete process, and the only required knowledge is the off-shell amplitude g * → η c +g. After some trick in factorizing the phase space integration, one ends up with a one-dimensional integral which exactly resembles what is encountered in the Collins-Soper approach. Not surprisingly, we again reproduce the results for d 0 (z) and d 2 (z) as given in (22).

Evolution of fragmentation functions
Recently LHC has measured η c production with p T 16 GeV [3], which is already considerably greater than the charm quark mass. One may worry that the large collinear logarithms α s ln p 2 T m 2 c n could potentially ruin the fixed-order prediction. These large logarithms are most conveniently resummed by invoking the famous Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation: where P j←i is the splitting kernel for parton i splitting into parton j. Therefore, in order to understand the evolution of g → η c fragmentation function, inevitably we also need knowledge of the c → η c fragmentation function, due to the non-vanishing off-diagonal splitting kernel P j←i . For simplicity, we neglect the contribution from the light quark/antiquark fragmenting into η c , since they are suppressed by additional powers in α s .
Here the relativistic correction also tends to dilute the fragmentation probability. We also compute the 2nd Mellin moment (average momentum fraction of η c ) for the c → η c fragmentation function: For the nonperturbative input parameters, we take O ηc 1 ≈ 0.244 GeV 3 , which is obtained from |R ηc (0)| 2 = 0.512 GeV 3 [21] through (6). We also take v 2 ηc = 0.13, which is obtained through the Gremm-Kapustin relation [22] by choosing the one-loop charm quark pole mass to be 1.4 GeV. The QCD coupling constants at initial scales are set as α s (2m c ) = 0.266, α s (3m c ) = 0.233, associated with the gluon and the c quark fragmentation, respectively.
We utilize the elaborate FORTRAN/C++ package HOPPET [24] to numerically solve the DGLAP evolution equation. We consider two types of evolution equations: the evolution by only implementing the diagonal splitting kernels (c → c and g → g) as the non-singlet (NS) (without mixing), as well as the evolution incorporating the off-diagonal splitting kernel as the singlet (S) (with mixing effect). The starting scales for the non-singlet evolution are set as 3m c and 2m c , for c → η c and g → η c fragmentation functions, respectively, roughly corresponding to the invariant masses of the final states of c → η c +c and g → η c + g. We choose the evolved scales to be 15 GeV and 45 GeV, respectively.
The singlet evolution is a more delicate issue. In the beginning we carry out non-singlet evolution for the gluon fragmentation function from 2m c to 3m c . Subsequently, by incorporating the mixing effect, we perform the singlet evolution for both fragmentation functions from 3m c to the scales 15 GeV and 45 GeV.
In Fig. 3 and Fig. 4, we display the evolution of the fragmentation functions with different energy scales, also including relativistic correction effects. The effect of relativistic corrections at different scales is quantified by the following factor: where D (0) is the LO fragmentation function, and D represents the LO+NLO fragmentation function. The reason we choose the second Mellin moments in (34), rather than the first moments, is that upon evolution, the gluon fragmentation function diverges as z → 0. This renders the numerical extrapolation unreliable, so we deliberately multiply by a factor of z to suppress the contribution from the small z region.
The numerical values of ∆ are listed in Table 1. The magnitudes of the O (v 2 ) corrections range from 8% to 24% with respect to the LO results. It should be noticed that ∆ NS g→ηc (µ) is independent of the evolution scale, which can be attributed to the fact that the NLO correction to the g → η c fragmentation function is proportional to the LO result. In particular, one predicts ∆ NS g→ηc = − 11 6 v 2 ηc = −0.238. After including the mixing effect, we observe that at large evolution scale, a remarkable rise of the c → η c fragmentation function occurs in the small z region. This can be attributed to the singular behavior of the splitting kernel of c → g as z → 0. Incorporating the mixing effect leads to drastically different shapes of the fragmentation functions with respect to those obtained from the nonsinglet evolution, especially for the c → η c fragmentation function.

Summary
In this work, we have calculated the leading relativistic correction to the fragmentation function for g → η c . Curiously, the order-v 2 fragmentation function shares the same functional z dependence as the LO one, differing only by an overall factor − 11 6 v 2 . The relativistic correction appears to considerably suppress the g → η c fragmentation probability. The effect of relativistic correc-tion to this fragmentation function is opposite to the NLO QCD radiative correction, which notably enhances the fragmentation probability for g → η c [19].
We have also confirmed the previous results of the order-v 2 correction to the c → η c fragmentation function. We further study the DGLAP evolution of the fragmentation functions for g → η c and c → η c . We find that incorporating the mixing effect in evolution has a striking impact on the small-z behavior of the fragmentation function.
When more copious η c samples are collected at highp T in future LHC experiments, it will be interesting to conduct comprehensive phenomenological analysis to test our understanding of the g, c → η c fragmentation mechanism.