Fragmentation function of $g\to Q\bar{Q}(^3S_1^{[8]})$ in soft gluon factorization and threshold resummation

We study the fragmentation function of the gluon to color-octet $^3S_1$ heavy quark-antiquark pair using the soft gluon factorization (SGF) approach, which expresses the fragmentation function in a form of perturbative short-distance hard part convoluted with one-dimensional color-octet $^3S_1$ soft gluon distribution (SGD). The short distance hard part is calculated to next-to-leading order in $\alpha_s$ and a renormalization group equation for the SGD is derived. By solving the renormalization group equation, threshold logarithms are resummed to all orders in perturbation theory. The comparison with gluon fragmentation function calculated in NRQCD factorization approach indicates that the SGF formula resums a series of velocity corrections in NRQCD which are important for phenomenological study.


Introduction
Heavy quarkonium physics has been on the focus of much experimental and theoretical attention since the discovery of the J/ψ in 1974. As the simplest bound state of strong interactions, heavy quarkonium is an ideal system to study both perturbative and nonperturbative aspects of QCD. Our current understanding of the decay and production of quarkonium is mainly based on the non-relativistic QCD (NRQCD) factorization [1], which factorizes physical quantities into perturbatively calculable short-distance coefficients (SDCs) multiplied by nonperturbative long-distance matrix elements (LDMEs).
However, recent studies shown that NRQCD factorization encounters some difficulties in describing inclusive quarkonium production data. In Ref. [2], it was argued that these difficulties may come from the bad convergence of velocity expansion in NRQCD. The velocity expansion suffers from large high order relativistic corrections due to ignoring the effects of soft hadrons emitted in the hadronization process. For this reason, the authors proposed a new factorization approach, called soft gluon factorization (SGF), to describe quarkonium production and decay [2]. It was argued in Ref. [3] that the SGF is equivalent to the NRQCD factorization, but with a series of important relativistic corrections originated from kinematic effects resummed. As a result, the SGF approach should has a much better convergence in the velocity expansion. Indeed, it was found in Ref. [2] that the lowest order result in velocity expansion in NRQCD factorization can deviate from the full SGF result by a factor of 4, but the lowest order velocity expansion in SGF can deviate from its full result by only a small value.
In the SGF approach, the quarkonium production cross section can be expressed in the following factorization formula  [2], and it can be simplified to the so called SGF-1d and SGF-0d formula by some further expansions. Using these SGF formulas, the authors studied the J/ψ hadroproduction via gluon fragmenting into polarization-summed 3 S [8] 1 intermediate state at tree level, and they found that the SGF-1d formula is a very good approximation of full SGF-4d formula. The SGF approach has also been applied to exclusive quarkonium processes [3,5], in which the SGDs are reduced to local matrix elements.
In this paper we study quarkonium inclusive production using gluon fragmentation function (FF) D g→H (z, µ), where z is the longitudinal momentum fraction of the quarkonium state and µ is the collinear factorization scale. To be concrete, we concentrate on the 3 S [8] 1 state production. In NRQCD framework, the SDC of gluon fragmenting to 3 S [8] 1 has been calculated up to next-to-leading order (NLO) in perturbative expansion [6][7][8] (but at the lowest order in velocity expansion). It was found that SDC in NRQCD suffers from large (threshold) logarithms of the form which spoils perturbative expansion in the region z → 1. Similar issues arise when analyzing the J/ψ photoproduction and electroproduction, where the large logarithms were resummed to all orders in α s [9][10][11][12] by combining the NRQCD with soft collinear effective theory (SCET) [13][14][15][16]. Although the same technique can also be used to resum the large logarithms appeared in FF, we choose to apply SGF instead of NRQCD+SCET to deal with this problem. The reason is that, as we will show, the SGF can not only resum the large logarithms but also resum a series of important relativistic corrections, which may result in better convergence for relativistic expansion. The rest of the paper is organized as following. In Sec. 2, we review the SGF formula for quarkonium production, including the definition of SGDs, the perturbative matching procedure and the velocity expansion of short distance hard parts. In Sec. 3, we present the perturbative calculation of the SGDs. In Sec. 4, we discuss the renormalization group equation (RGE) for the 3 S [8] 1 SGD. In Sec. 5, we study the gluon FF in SGF and calculate the short distance hard part in 3 S [8] 1 channel up to NLO. We conclude in Sec. 6. In appendix A we provide some details of solving the RGE of SGDs. In appendix B, we provide some integral formulas used in our calculation. Finally, we list the expressions of finite functions that enter the gluon FF in appendix C.

Soft gluon factorization
Before studying the gluon fragmentation function, we first briefly review the SGF formula for quarkonium production cross section. We denote M H and m Q as the mass of the heavy quarkonium H and the mass of the heavy quark Q, respectively. Lorentz-vector a, denoting as is sometimes also expressed in light-cone coordinates, a µ = (a + , a − , a 1 , a 2 ) = (a + , a − , a ⊥ ), The scalar product of two four-vector a and b then becomes We introduce a light-like vector l µ = (0, 1, 0 ⊥ ) so that a · l = a + . Following the discussion in [2], the factorization formula Eq. (1.1) for producing a quarkonium H with momentum P H can be simplified to the following SGF-1d formula which has similar good convergence in velocity expansion for our purpose, and at the same time is more suitable in practical use. In the above formula, dσ [nn ] (P H /z, m Q , µ f ) are the short-distance hard parts which, roughly speaking, produce a QQ pair with momentum P H /z and quantum number n in the amplitude and n in the complex conjugate of the amplitude, F [nn ]→H (z, M H , m Q , µ f ) are the one-dimensional SGDs which describe the hadronization from the QQ pair to quarkonium H, µ f is the factorization scale and z = P + H /P + is the longitudinal momentum fraction with P denoting the total momentum of the intermediate QQ pair.
SGDs are defined as where Ψ stands for Dirac field of heavy quark and the subscript "S" means the field operators in the above definition are the operators obtained in small momentum region. In additional, in this paper we define "S" to select only leading power terms in threshold (P − P H ) + = (1 − z)P + expansion, which is sufficient to factorize and then resum leading-power large logarithms as z → 1. In Eq. (2.2), the quarkonium state, created by a † H , has standard relativistic normalization and the projection operators K n can be decomposed to a spin operator, a color operator, and a gauge link [2]. For the case n = 3 S 1,λ which will be studied in this work, we have where D µ is the gauge covariant derivative with Ψ ← → D µ Ψ = Ψ(D µ Ψ) − (D µ Ψ)Ψ and λ are polarization vectors which satisfy the following relations The spin projection operator P αβ is defined as The color operator is defined as The gauge link Φ l (rb − )ā a is introduced to enable gauge invariance of SGD, which is defined along the l µ direction, where P denotes path ordering and A µ (x) is the gluon field in the adjoint representation: As the short distance hard parts dσ nn (P H /z, m Q , µ f ) do not depend on nonperturbative physics, they can be perturbatively calculated. To this end, we replace the quarkonium H in Eq. (2.1) by an on-shell QQ pair with certain quantum number m in the amplitude and m in the complex-conjugate amplitude, which results in We can determine the hard parts by computing both sides of Eq. (2.8) in perturbation theory. If the factorization Eq. (2.1) holds, the obtained hard parts will be free of infrared divergences. Momenta of the on-shell QQ pair are chosen as J,Jz ), we replace the spinors of QQ by the following projector 1,λ or 3 D The color operators and spin operator have similar definitions as those in Eq. (2.3), which are given by  . Finally, we note that perturbative calculation with analytical q dependence are not very easy for complicated processes. There are at least two independent hard scales, m Q and M H /z in the hard parts. As suggested in [2], we can further simplify the hard parts by expanding m 2 with v being the velocity of the heavy quark in the rest frame of the pair. Then the SGF formula Eq. (2.1) can be rewritten as which defines a velocity expansion series in SGF. The tree level calculation in Ref. [2] shows that the convergence of velocity expansion in SGF is much better than that in NRQCD. We will confirm this conclusion at one-loop level in this paper.

Definition
For the purpose of this paper, we are interested in the intermediate state n or n equals to state 3 S The notations in the second line denote polarization-summed intermediate states. In general, the SGD F [LL,λ]→H is λ dependent even for polarization-summed H. This is different to the polarized NRQCD LDMEs, which can be simplified to the usual unpolarized NRQCD LDMEs due to the rotation invariance [4,17]. But in the case of SGD, the z-axis direction needs to be specified because the longitudinal momentum of intermediate QQ pair is fixed, thus the rotation invariance is broken. Similar to the polarized NRQCD LDMEs defined in [4], it is convenient to construct the definitions of polarized SGDs as But for the polarization-summed S-wave quarkonium H, such as J/ψ, the polarized 3 S SGDs can be reduced to following unpolarized SGD: And we have We will explain this later.
After replacing H with QQ[LL, λ], the general form of F [L L ,λ ]→QQ[LL,λ] can be written as where dΦ is the final-state phase space, k i denotes momenta of the final-state gluons or light (anti-)quarks, and A L L (P H , P, m Q , k i , λ , λ) denotes the amplitude of transition from a QQ pair in 3 L 1,λ state which can be written as where A QQ L is the amplitude to produce an open QQ pair with spinors of the pair removed. As pointed out in the last section, in addition to force loop momenta and light-parton momenta k i in small-momentum region, the operator T S also selects only the leading power contributions in the threshold expansion.
Based on Eq. (3.5), we find that, for polarization-summed final state, SGDs can be decomposed as For general SGDs, F 2 and F 3 are nonzero and the existence of l µ breaks the rotation invariance. But for the S-wave case, F 2 and F 3 are zero thanks to keeping only LP term in threshold expansion, which results in Thus we argue that for the polarization-summed S-wave quarkonium H, we have the relations in Eq. (3.4).

LO calculation
The where Then we obtain which is consistent with the result in Eq. (2.14). Similarly we can derive To obtain the above results, we have used the relations |λ|=1 µ λ * ν The longitudinally polarized (λ = L) SGDs can be calculated similarly, but they are irrelevant for our purpose as we will show later, thus we do not consider their contributions in our calculation.

NLO calculation
At NLO, the corresponding Feynman diagrams are shown in Fig. 2 and Fig. 3. Let us first consider virtual corrections in Fig. 2. Following the calculation in ref. [3], we obtain with For real corrections in Fig.3, we have with Here k is the momentum of the emitted real gluon in the final state, which is restricted in the soft domain. Applying T S to I , we obtain Performing the Ω and k integration, we obtain the following result Combing the virtual and real correction contributions shown in Eqs. (3.15) and ( The ultraviolet divergences in the above result can be removed by using the MS renormalization procedure, and we get the renormalized SGD Here we distinguish the dimensional regularization scale µ c from the factorization scale µ f . Similarly we can calculate the other SGDs by projecting the initial and finial QQ pair into corresponding states, and we have Here T is a hypergeometric function which can be expended as

Renormalization group equation
Based the factorization formula in Eq. (2.1), RGEs for SGDs have the following general form The evolution equation with this kernel is difficult to be solved analytically. However, if we only want to resum large logarithms at leading order in v 2 expansion, we only need to keep the leading power term in the evolution kernel, which reads where Then the evolution equation for To solve the above evolution equation, we rewrite it as a function of the variable with ω = M H (1/x − 1). To derive the above equation, we have used the rescaling identity for plus functions [18] κ ln n (κω) κω The plus functions of the dimensionful variable ω are defined as where f (ω) is a smooth test function. To solve the RGE Eq. (4.7), it is convenient to perform a Laplace transformation [19][20][21], which, together with its inverse, is given bỹ where the constant c is chosen to be larger than the real part of the rightmost singularity off (s). Taking a Laplace transform in Eq. (4.7) we have a simple multiplicative RGE wheres = se γ E . Solving the RGE in Laplace space is straight forward, details of which are given in App. A. We eventually obtaiñ whereF fix (M H ) is given in Eq. (A13) which recovers fixed-order perturbative results at small s region, functions h 0,1 are given by Eq. (A9), and χ * is defined as where n f represents the number of light flavors. We have introduced a "frozen" scales * withs to prevent Landau singularity, where a (with a 1) is a parameter and its dependence should be compensated with a change of nonperturbative input. Therefore, we fix it as a = 1.3 in the rest of this paper. Finally all nonperturbative information is included in the modelF mod (s). Then SGD can be obtained by transforming Eq. (4.12) back to momentum space using where c is a constant that must be greater than the real part of all singularities ofF [SS]→H .

Sensitivity of nonperturbative model
We will demonstrate that SGD obtained above by solving RGE is only sensitive to two parameters of the nonperturbative model. To this end, we will study various models, with significant different shapes. We choose the model function used in [9][10][11]22] as our first class of models,  +1), respectively. Here N H determines the normalization,Λ characterizes the average radiated momentum, and b is related to the width of the model function. We varyΛ from 0.5GeV − 0.7GeV 1 and b from 1 − 3. Other models are chosen to cover various possible shapes, but with the same normalization (M H N H ) and average radiated momentum (0.6GeV). All models that we will study are listed in the following, Model-7:  These input models are shown in the first row in Fig. 4, and the corresponding SGDs obtained by solving RGE at LO are given in the second row, where we have taken M H = QCD = 0.25GeV and N H = 1/3. From the plots we find that, comparing with input models, SGDs are much broader and their peaks are shifted to smaller x. More importantly, SGDs at moderate and small x, i.e., x < 0.65, are almost independent of input models, as far as the models have the same normalization and average radiated momentum. This is reasonable because small-x region is dominated by perturbative effects. At larger x region, values of SGDs depend on input models, but the sensitivity after RGE running is much weaker than the original models. To see this more clearly, we introduce a function to describe the difference between two functions g 1 (x) and g 2 (x) (4.19) Differences between input models and that between corresponding SGDs are shown in Table. 1. It is clear that differences are significantly decreased after performing resummation using RGE. Especially for models with the same average radiated momentum, differences have been typically reduced by a factor of 1/5. The above phenomenon can be interpreted as following. The overall normalization of a SGD is fully determined by the normalization of the input model. The shape of a SGD is almost fully determined by perturbative calculation with all-order resummation, although its peak can be affected by the average radiated momentum of the input model. In this consideration, in practical use it is sufficient to choose a suitable parameter b in Eq. (4.17), and fit the parameters N H andΛ by comparing with experiment data.

SGF formula for gluon fragmentation function
The Collins-Soper definition for fragmentation function of a gluon fragmenting into a hadron (quarkonium) H is given by [23] where G µν is the gluon field-strength operator, P H and P c are the momenta of the hadron and initial virtual gluon, respectively, and z is the "+" momentum fraction defined as z = P + H /P + c . The projection operator P H(P H ) is given by where X sums over all unobserved particles. The gauge link Φ l (ξ − ) is the same as the one given in Eq. (2.7). The definition of Eq. (5.1) is gauge invariant and we use the Feynman gauge in our calculation. In Eq. (5.1), µ is the collinear factorization scale. The dependence of FFs on µ is controlled by the DGLAP evolution equations [24][25][26], where we have ignored contributions of quark fragmentation functions because they are usually less important. Evolving FFs from a low scale µ up to a high scale µ h by solving the DGLAP evolution, one can resum large logarithms of µ 2 /µ 2 h to all orders in perturbation theory. The gluon splitting function P gg has the perturbative expansion where, e.g., P LO gg (x) is given by In the parton (initial virtual gluon) frame, the gluon FF to quarkonium has the form Using the SGF-1d formula Eq. (2.18) and changing the variables from the parton to the hadron frame, we have By defining the short distance hard partsD [nn ] (ẑ, M H /x, m Q , µ, µ f ) as followinĝ where we have use the fact that polarized SGDs for S-wave states can be related polarizationsummed SGDs, as noted in Eq. (3.4). Similar to Eq. (2.8), we replace the quarkonium H in Eq. (5.10) by an on-shell QQ pair with specific quantum number and momenta P H . Then we have following matching relationŝ These relations enable us to calculateD [SS] perturbatively to NLO.

Perturbative calculation of gluon fragmentation functions
For the fragmentation functions we have where dΦ is the final state phase space, k i denotes momenta of the final-state gluons, and M is projected amplitude using projector introduced in Eq. (2.11), i.e.
where A ρ g→QQ is the amplitude from a virtual gluon to a QQ pair plus light partons with spinors of QQ removed, Ω is the solid angle of relative momentum q in the QQ rest frame, and N Ω is given by (5.14) At LO in α s , the phase space dΦ LO is given by Then the calculation of perturbative FFs is straightforward, which gives To compute the contributions of these diagrams, we rewrite the integral of solid angle as We then change all delta functions to propagator denominators [27] by using Then the obtained new Feynman integrals can be calculated analytically by using established multi-loop calculation techniques, like the IBP reduction [28][29][30][31][32] and differential equations [33][34][35][36][37][38], which results in =0. (5.20g) In the above results we have dropped imaginary parts that are irrelevant for our purpose.
where k is the momentum of the final-state gluon. In the calculation, we first integrate out the solid angle. To do this, we decompose the amplitude M ρ L into the following form where the functions M Li and M Li are q independent. Then the integral over Ω can be performed by using equations provided in Appendix. B. After that, we carry out the integration over the phase space dΦ (1) , which results in where R(z, M H , ∆) is infrared finite.
Combining the virtual and real corrections in Eq. (5.20) and Eq. (5.23), we obtain NLO correction for the FF D g→QQ [SS] , which contains ultraviolet and infrared divergences in the form of single poles in as shown below Ultraviolet divergences in the above equation are canceled by the renormalization of the coupling constant α s and the operator defining the FF. The renormalization of α s in the MS scheme can be carried out by making the following replacement in Eq. (5.16) 25) and the operator renormalization in the MS scheme can be carried out by further making the following replacement Then the renormalized D g→QQ [SS] reads Here we have distinguished between the dimensional regularization scale µ c and the collinear factorization scale µ. All remaining divergences in the above equation are infrared divergences, and they should be absorbed into the definition of SGDs if the SGF is valid at NLO. The infrared finite functions A, B, C, R, F in above equations are listed in App. C. They also can be read from the ancillary file. Based onD [SS] , we can easily obtain the short distance hard part at leading order in the velocity expansion, which readŝ

Numerical results
We now present our numerical results for the gluon FFs given in Eq. (5.29). We use the model function Eq. , and n f = 3 in β 0 which means contributions from virtual or initial heavy quarks are ignored. In the following, in not specified we refer to fully NLO FF D g→H (z, M H , m Q , M H ) calculated in SGF withΛ = 0.6GeV. We note that, as NLO evolution kernel for SGD is still not available, all results presented in this paper are obtained by using the LO evolution kernel (4.4). g→H is a good approximation to D g→H , with deviation smaller than 20%. This implies that the convergence of velocity expansion in SGF is good and that D g→H is insensitive to the input of heavy quark mass m Q . For comparison, we also provide plots of FFs obtained by setting x = 1 in hard part (5.28) everywhere except in the delta function δ(1 − z/x). This is clear a bad approximation because it overshoots original results by at least a factor of 2, which indicates big corrections at high order in 1 − x expansion. 2 This can be easily understood because the hard part are proportional to x 3 and SGD is peaked around x = 0.75. We note that, if one uses the shape function method [9][10][11][12] to calculate the gluon FF, one should take x → 1 in hard part, which will result in large relativistic corrections as discussed above.
In the right figure of Fig. 7 we show theΛ dependence by choosing it to 0.05GeV, 0.5GeV, 0.6GeV and 0.7GeV. We find that the gluon FF at small z is much less affected by the parameterΛ than those at large z, which indicates that gluon FF at small z is dominated by perturbative effects while that at large z is sensitive to nonperturbative dynamics. This sensitivity provides a possibility to extractΛ using experimental data. We note that the plot withΛ = 0.05GeV in Fig. 7 is to mimic the case withΛ → 0. Although nonperturbative input models withΛ = 0.6GeV and 0.05GeV have normalized to the same value, they result in significantly different FFs, mainly due to the large logarithm resummation. To quantify the difference, we calculate the n-th moment of FFs and define the following ratios, with numerical results given in Table. 2 for n = 2, 3, 4, 5, 6. By taking D g→H as "exact" result, the values of R indicate that the lowest order in velocity expansion in SGF (denoted as R (0) ) is a good approximation, while changingΛ from 0.6GeV to 0.05GeV (denoted as R 0.05 ) results in large deviation, about a factor of 3. As shown in Fig. 8, by perturbatively expanding DΛ =0.05GeV g→H to NLO in α s and replacing M H by 2m Q (denoted as "pert."), we can nicely recover NRQCD result [6][7][8]. Their deviations at end points of z are due to the nonzeroΛ, but they have almost the same moments as shown in Table. 2, where R NRQCD means that the numerator of the ratio is the NRQCD FFs. We thus find that the NRQCD results overshoot SGF results with our prefered treatment by a factor about 5 − 7, which is consistent with the conclusion in Ref. [2]. One of the main reasons for deviation is due to the large logarithm resummation, and the other one is due to the choice ofΛ = 0.6GeV instead of zero.

Summary
In this paper we studied the FF from the gluon to color-octet 3 S 1 heavy quark-antiquark pair in SGF approach, which is expressed as the convolution of perturbative short-distance hard part with the one-dimensional color-octet 3 S 1 SGD in Eq. (5.9). We calculated the color-octet 3 S 1 SGD to NLO in Eq. (3.24) and derived a RGE for it in Eq. (4.1). By solving the RGE, we resummed the threshold logarithms to all order in perturbation theory. Based on this, we perturbatively calculated the short distance hard part of the gluon FF up to NLO in Eq. (C6). With a natural scales choice, the hard part is free of threshold logarithms, which has been resummed by the RGE of SGD. This demonstrated the validity of SGF at NLO.
Our numerical results shown that velocity expansion in the SGF has a good convergence, and the lowest order approximation can already capture most physics. By a specific choice of parameters (but not prefered), the SGF can reproduce the NRQCD results. But then we get a big deviation, as large as a factor of 6, from our prefered results. This implies that velocity expansion in NRQCD has a very slow convergence. We also show that even the large threshold logarithms in NRQCD results are resumed by using shape function method, there should be still significant velocity corrections.
Our results may provide a new insight to understand the mechanisms of quarkonium production, especially for quarkonium production at high energy colliders and quarkonium production within a jet [39,40].
The functions h 0,1 are given by Note that the term ln(1 − 2χ) in above expressions gives rise to a cut singularity that start at the branch points This singularity, as a result of the divergence of the running coupling α s (µ) near the Landau pole at µ = Λ QCD , signals the onset of nonperturbative phenomena at very large s or, equivalently, when x is very close to the value x = 1. To deal with the Landau singularity, we introduce a cut-offs * above which the functions h 0,1 remained constant, i.e. "frozen" ats * . Additionally, a nonperturbative model was introduced to correct the formalism for nonperturbative effects at very large s. Following [41], we make the replacement s →s 1 + 1/(M Hs * ) 1 +s/s * ,s * =s L /a (A11) in Eq. (A9). Here a is a parameter not smaller than one, but of order one. Such a replacement prevents χ from entering the nonperturbative regime. Then we model the initial SGDF And the modelF mod (s) is introduced to describe the nonperturbative effects.

B Integration over solid angle
In this appendix, we give the expressions that integrate over the solid angle Ω of q in the QQ rest frame. For the momentum k which satisfies k 2 = 0, we have following equations, d d−2 Ω q µ q ρ q λ (P H − 2q) · k = d d−2 Ω −q µ q ρ q λ (P H + 2q) · k = F 3 (P µρ P λσ k σ + P µλ P ρσ k σ + P ρλ P µσ k σ ) Where