Gluon fragmentation into 3 P [1,8] J quark pair and test of NRQCD factorization at two-loop level
PPrepared for submission to JHEP
Gluon fragmentation into P [1 , J quark pair andtest of NRQCD factorization at two-loop level Peng Zhang, a,b
Ce Meng a Yan-Qing Ma, a,b,c and Kuang-Ta Chao a,b,c a School of Physics and State Key Laboratory of Nuclear Physics and Technology, Peking University,Beijing 100871, China b Center for High Energy Physics, Peking University,Beijing 100871, China c Collaborative Innovation Center of Quantum Matter,Beijing 100871, China
E-mail: [email protected] , [email protected] , [email protected] , [email protected] Abstract:
The next-to-leading order (NLO) ( O ( α s ) ) corrections for gluon fragmenta-tion functions to a heavy quark-antiquark pair in P [1 , J states are calculated within theNRQCD factorization. We use integration-by-parts reduction and differential equationsto semi-analytically calculate fragmentation functions in full-QCD, and find that infrareddivergences can be absorbed by the NRQCD long distance matrix elements. Thus, theNRQCD factorization conjecture is verified at two-loop level via a physical process, whichis free of artificial ultraviolet divergences. Through matching procedure, infrared-safe shortdistance coefficients and O ( α s ) perturbative NRQCD matrix elements (cid:104)O P [1 / J ( S [8]1 ) (cid:105) NLO are obtained simultaneously. The NLO short distance coefficients are found to have signif-icant corrections comparing with the LO ones.
Keywords:
Perturbative calculations, Factorization, Heavy quarkonia a r X i v : . [ h e p - ph ] N ov ontents Heavy quarkonium provides an ideal physical system to explore the strong interaction, asthe heavy quarkonium production contains both perturbative and non-perturbative effectsin QCD. Non-relativistic QCD (NRQCD) factorization [1] is the most widely used theoryto explain quarkonium production so far. When transverse momentum p T of the producedquarkonium is large, fragmentation mechanism dominates and factorization is easier tohold, because long-distance interactions between quarkonium and initial-state particles aresuppressed then.Inclusive production differential cross section of a specific hadron H at high p T can becalculated in collinear factorization [2], d σ A + B → H ( p T )+ X = (cid:88) i dˆ σ A + B → i ( p T /z )+ X (cid:48) ⊗ D i → H ( z, µ ) + O (1 /p T ) , (1.1)where i sums over all quarks and gluons, z is the light-cone momentum fraction carriedby H with respect to the parent parton i , and A and B are colliding particles whose ef-fect should be further factorized into partons if they are hadrons. d ˆ σ A + B → i ( p T /z )+ X areperturbatively calculable hard parts, while D i → H ( z, µ ) are non-perturbative but universalfragmentation functions (FFs) describing the probability of parton hadronizing to H withmomentum fraction z . For heavy quarkonium production, an important O (1 /p T ) contribu-tions are double parton FFs which can also be factorized [3–7]. In FFs, there is a collinearfactorization scale µ dependence, which cancels with the similar dependence in hard partsperturbatively order by order, leaving physical differential cross section independent of thescale. The evolution of single parton FFs with respect to µ are controlled by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation [8–10], and similar evolution– 1 –quations for double parton FFs are calculated in [4]. With these evolution equations, theonly unknown information for FFs are their values at any chosen factorization scale µ = µ f .When µ f is close to the quarkonium mass m H , it is natural to calculate FFs via NRQCDfactorization. For single parton FFs that will be considered in this paper, we have D i → H ( z, µ f ) = (cid:88) n d i → Q ¯ Q ( n ) ( z, µ f ) (cid:104) ¯ O Hn (cid:105) , (1.2)where d i → Q ¯ Q ( n ) represent the perturbative calculable short-distance coefficients (SDCs)to produce a heavy quark-antiquark pair Q ¯ Q with quantum number n , and (cid:104) ¯ O Hn (cid:105) arenormalized long-distance matrix elements (LDMEs) . The quantum number is usuallyexpressed in spectroscopic notation n = S +1 L [ c ] J , with c = 1 , respectively for color-singletstate and color-octet state. According to velocity scaling rule [1], (cid:104) ¯ O Hn (cid:105) is usually suppressedif L is too large. Therefore, the most important states for phenomenological purpose are S -wave and P -wave states. Because LDMEs are supposed to be process independent, they canbe determined by fitting experimental data, while SDCs can be calculated perturbativelythrough the matching procedure.All the SDCs for single parton FFs to both S -wave and P -wave states up to O ( α s ) areavailable [11–21] (see [22, 23] for a summary and comparison). At α s order, the leading order(LO) SDCs of g → Q ¯ Q ( S [1]1 ) + X and g → Q ¯ Q ( P [1]1 ) + X are calculated separately in Refs.[12, 24–27] and Ref. [28], and the next-to-leading order (NLO) SDCs of g → Q ¯ Q ( S [1 , )+ X have also been calculated in Refs. [29–31] recently. To perform phenomenological analysisat a consistent precision, such as that for J/ψ or χ cJ production, at this order, one alsoneeds the NLO SDCs of g → Q ¯ Q ( P [1 , J )+ X and the next-to-next-to-leading order (NNLO)SDCs of g → Q ¯ Q ( S [8]1 ) + X . In this paper, we will focus on the former.Another motivation to calculate the NLO FFs for processes g → Q ¯ Q ( P [1 , J ) + X is to test the NRQCD factorization conjecture at two-loop level with a physical process.The test at two-loop level exists in literature [32–35] but within the framework of eikonalapproximation for infrared (IR) divergences. The eikonal approximation introduces artificialultraviolet (UV) divergences which makes the extraction of IR divergences beyond one-looplevel highly nontrivial. Because we do not use eikonal approximation, the extraction of IRdivergences in our case is straightforward.In this paper, we aim to calculate NLO SDCs of g → Q ¯ Q ( P [1 , J ) + X to high precisionusing similar methods in our previous paper [29]. The rest of the paper is organized asfollowing. In Sec. 2, we calculate the process to LO and NLO in full-QCD. The calcu-lation method is very similar to that in Ref. [29]. A small improvement is the compu-tation of operator renormalization shown in Sec. 2.3 and Appendix. A. We argue thatthe NRQCD factorization at two-loop level without eikonal approximation is verified inSec. 2.4. The O ( α s ) perturbative NRQCD matrix elements (PMEs) (cid:104)O P [1] J ( S [8]1 ) (cid:105) NLO and (cid:104)O P [8] J ( S [8]1 ) (cid:105) NLO are matched, and the corresponding IR-safe SDCs are obtained simulta-neously in Sec. 3. Discussion will be presented in Sec. 4. (cid:104) ¯ O Hn (cid:105) can be related to the original definition of NRQCD LDME (cid:104)O Hn (cid:105) [1] by the following rules. Theyare the same if n is color-octet, and (cid:104) ¯ O Hn (cid:105) = (cid:104)O Hn (cid:105) / (2 N c ) if n is color-singlet. – 2 – Full-QCD Calculation
FF of a gluon to a hadron (quarkonium) is defined by Collins and Soper [36], D g → H ( z, µ ) = − g µν z D − πP + c ( N c − D − (cid:90) + ∞−∞ d x − e − iP + c x − × (cid:104) | G + µc (0) E † (0 , , ⊥ ) cb P H ( P ) E (0 , x − , ⊥ ) ba G + νa (0 , x − , ⊥ ) | (cid:105) , (2.1)where G µν is the gluon field-strength operator, P and P c are respectively the momentaof the produced hadron H and the initial-state fragmenting gluon g , and z = P + /P + c is the ratio of momenta along the “+” direction. It is convenient to choose the framein which the hadron has zero transverse momentum, P = ( zP + c , M / (2 zP + c ) , ⊥ ) , with P = 2 P + P − = M . The projection operator P H ( P ) is defined by P H ( P ) = (cid:88) X | H ( P ) + X (cid:105)(cid:104) H ( P ) + X | , (2.2)where X sums over all unobserved particles. The gauge link E ( x − ) is an eikonal operatorthat involves a path-ordered exponential of gluon field operators along a light-like path, E (0 , x − , ⊥ ) ba = P exp (cid:20) + ig s (cid:90) ∞ x − d z − A + (0 , z − , ⊥ ) (cid:21) ba , (2.3)where g s = √ πα s is the QCD coupling constant and A µ ( x ) is the matrix-valued gluonfield in the adjoint representation: [ A µ ( x )] ac = if abc A µb ( x ) .Since the SDCs d ’s in the NRQCD factorization formula (1.2) are independent of thefinal state H , we can replace H by on-shell Q ¯ Q states to extract d i → Q ¯ Q [ n ] through thematching procedure [See below (3.1) and (3.6)].Choosing the factorization scale µ ∼ m Q with m Q being the heavy quark mass,we can calculate FF for g → Q ¯ Q [ n ] in the full-QCD with the Feynman gauge. Most ofthe Feynman rules are the same as those of QCD, while the others related to the eikonalline can be found in Ref. [29]. Feynman amplitudes are denoted as M λ Q λ ¯ Q λ i ( q ) , where λ Q and λ ¯ Q are respectively spins of produced on-shell heavy quark and heavy antiquark, λ i ( i = 1 , , . . . ) are spins of the initial-state virtual gluon or final-state unobserved lightparticles, and q is the heavy quark momentum in the rest frame of heavy quark pair. Toproject the free Q ¯ Q pair to the state with spin-triplet and specific color, one can multiplyamplitudes by projection operators and obtain M αλ i ( q ) = Tr (cid:104) Γ c Γ αλ M λ Q λ ¯ Q λ i ( q ) (cid:105) , (2.4)where M λ Q λ ¯ Q λ i denotes amplitude with external heavy-(anti)quark spinors removed and– 3 –he projection operators are defined as Γ c =1 = 1 √ N c , Γ c =8 = √ T a (cid:112) N c − , Γ αλ = 1 √ E ( E + m Q ) ( /p − m Q ) 2 E − /P E γ α E + /P E ( /p − m Q ) , (2.5)where p = P/ q and p = P/ − q are momenta respectively for Q and ¯ Q , m Q is the massof heavy quark, and E is the energy of Q or ¯ Q in the rest frame of Q ¯ Q . With the on-shellconditions of Q and ¯ Q , we have P · q = 0 , P = 4 E = 4( m Q − q ) . (2.6)To project the amplitude into P-wave, we also need to take the derivation over q β , M αβλ i = ∂∂q β M αλ i ( q ) | q → . (2.7)And by summing over spin and color of initial-state and final-state particles, we get thesquared amplitude A J = N CS N J (cid:88) P αα (cid:48) ββ (cid:48) J Re (cid:104) M αβλ i M ∗ α (cid:48) β (cid:48) λ i (cid:105) , (2.8)where N CS = z D − ( N c − D − with D = 4 − (cid:15) is the space-time dimension. For different J , N J and P αα (cid:48) ββ (cid:48) J are respectively N J = , for J = 012 ( D − D − , for J = 112 ( D + 1)( D − , for J = 2( D − , for (cid:88) J , (2.9)and P αα (cid:48) ββ (cid:48) J = D − I αβ I α (cid:48) β (cid:48) , for J = 012 (cid:16) I αα (cid:48) I ββ (cid:48) − I αβ (cid:48) I α (cid:48) β (cid:17) , for J = 112 (cid:16) I αα (cid:48) I ββ (cid:48) − I αβ (cid:48) I α (cid:48) β (cid:17) − D − I αβ I α (cid:48) β (cid:48) , for J = 2 I αα (cid:48) I ββ (cid:48) , for (cid:88) J , (2.10)where I αβ = − g αβ + P α P β P . (2.11)– 4 –hen, fragmentation function in the full-QCD can be denoted as D [ g → Q ¯ Q ( P [1 / J ])] = (cid:90) dΦ A J . (2.12)Here, dΦ is the final-state phase space, dΦ = P + z S δ (cid:32) − zz P + − (cid:88) i k + i (cid:33) (cid:89) i d D k i (2 π ) D − δ + ( k i ) , (2.13)where S is the symmetry factor for final-state particles and k i ’s are momenta of final-statelight particles. In the rest of the paper, we will take the case of g → Q ¯ Q ( P [1]0 ) + X as anexplicit example to explain our calculations. The Feynman diagrams of gluon fragmenting into Q ¯ Q [ P [1]0 ] at the LO in α s are shown inFigure 1. These diagrams can be easily calculated, and the obtained FF can be expressed P + q P − qk Figure 1 . One of the two Feynman diagrams of gluon fragmenting into Q ¯ Q [ P [1]0 ] at LO in α s .Another diagram can be obtained by permuting the heavy quark and anti-quark. as D LO [ g → Q ¯ Q ( P [1]0 )] = N LO − (cid:15) (cid:2) (1 − z ) − (cid:15) f ( z ) + (1 − z ) − − (cid:15) f ( z ) (cid:3) , (2.14)where N LO = α s ( D − N c m Q (cid:32) πµ r m Q (cid:33) (cid:15) Γ(1 + (cid:15) ) (2.15)and f ( z ) = 3(5 − z ) (cid:15) + 32 (3 z − − z (cid:15) ,f ( z ) = 3(1 − z )(3 z − (cid:15) + 13 z − z + 15 z + 183 + 2 (cid:0) z − z − (cid:1) z(cid:15) + 2 z (cid:15) . (2.16)With the following decomposition (1 − z ) − − (cid:15) = − δ (1 − z )2 (cid:15) + 1[1 − z ] + − (cid:15) (cid:20) ln(1 − z )1 − z (cid:21) + + O ( (cid:15) ) , (2.17)the FF can be expressed in the limit of (cid:15) → as D LO [ g → Q ¯ Q ( P [1]0 )] = N LO (cid:16) − (cid:15) δ (1 − z ) + 13 δ (1 − z ) + 43 z [1 − z ] + + 16 (85 − z ) z + 3(5 − z ) ln(1 − z ) (cid:17) . (2.18)– 5 –t is clear that there is an IR divergence in the full-QCD result. This divergence can beabsorbed by the color-octet PME, which will be explained in Sec. 3. The NLO full-QCD calculations of this process resemble closely to the NLO calculations of g → Q ¯ Q ( S [1 , ) + X in Ref. [29]. The relevant Feynman diagrams are shown in Figure 2and Figure 3. In our calculation, we need to cut these diagrams at some specific positions.In Figure 4, we give an example to show how to obtain cut diagrams for the real and thevirtual corrections. Figure 2 . Typical Feynman diagrams without gluon eikonal vertex for g → Q ¯ Q ( P [1]0 ) + X . Theother diagrams can be obtained by permuting the heavy quark and anti-quark. The real and virtualdiagrams can be obtained by cutting at different positions. The calculation is almost the same as that in Ref. [29]. Firstly we generate theamplitude from the Feynman diagrams for real or virtual corrections. Then reduce theminto linear combinations of master integrals (MIs). There are 95 MIs for real correctionsand 66 MIs for virtual corrections. Fortunately, these MIs are the same as those in Ref.[29], which can be calculated by using differential equations (DEs) method [37–49]. We canestimate values of these MIs in regions ∼ / , / ∼ / and / ∼ respectively by theasymptotic expansions of them at z = 0 , / and . Finally, the real and virtual correctionscan be expressed by asymptotic expansions in different regions as (cid:88) s n s (cid:88) i =0 ( z − z ) s ln i ( z − z ) ∞ (cid:88) j =0 I s i j ( (cid:15) )( z − z ) j , (2.19)where s are linear functions of (cid:15) , n s are integers determined by s , and coefficients I s i j ( (cid:15) ) are functions of (cid:15) , which can be computed to sufficient high order with high precision.– 6 – igure 3 . Typical Feynman diagrams with gluon eikonal vertex for g → Q ¯ Q ( P [1]0 ) + X . The otherdiagrams can be obtained by permuting the heavy quark and anti-quark. The real and virtualdiagrams can be obtained by cutting at different positions. P + q P − q kl kl P + q P − q Figure 4 . Real and virtual cut diagrams by cutting at different positions for the last diagrams inFigure 3.
Renormalization has some differences from that in Ref. [29], because /(cid:15) pole alreadypresents at LO FF. Firstly, the on-shell renormalization constants of heavy quark field andmass ( δ and δ m ) should be expanded at least to O ( (cid:15) ) . To simplify the procedure, we useexpression with exact (cid:15) dependence, δ i = α s π Γ(1 + (cid:15) ) (cid:32) πµ r m Q (cid:33) (cid:15) ˆ δ i , (2.20)with ˆ δ = − C F − (cid:15) (cid:18) (cid:15) UV (cid:15) − (cid:15) + 2 (cid:15) IR (cid:19) , ˆ δ m = − C F − (cid:15) (cid:15) UV − (cid:15) − (cid:15) , (2.21)where C F = ( N c − / (2 N c ) .Apart from the calculations of the QCD counter terms, some tricks are needed in thecalculation of operator renormalization. In the MS scheme, the operator counter term isgiven by D OperatorNLO ( z ) = − α s π Γ(1 + (cid:15) ) (cid:15) (cid:18) πµ r µ (cid:19) (cid:15) (cid:90) z d yy P gg ( y ) D LO (cid:18) zy (cid:19) , (2.22)– 7 –here µ is the collinear factorization scale, and the Altarelli-Parisi splitting function P gg ( z ) is P gg ( z ) = b δ (1 − z ) + 2 N c (cid:18) z (1 − z ) + + 1 − zz + z (1 − z ) (cid:19) , (2.23)where b = (11 N c − n f ) / . If one takes the Eq. (2.18) as LO FF D LO ( z ) , there willbe two plus functions convolved in the integral in (2.22), which is difficult to define andcalculate. However, we can take Eq. (2.14) as the input form of the D LO ( z ) to overcome thisdifficulty, with corresponding integrated results given in Appendix A. In this way, operatorrenormalization results can be easily obtained in each region, which are convenient to beadded to the real and virtual corrections. In the full-QCD calculations, the FF can be expressed by the phase space integral of am-plitude squared as shown in Eq. (2.12). Before taking the derivative over q , the amplitudesquared at O ( q ) can be expressed as A ( q, q (cid:48) ) = ( q · q (cid:48) ) A + ( q · n )( q (cid:48) · n )( P · n ) A , (2.24)where n is the light-like vector that defines the direction of the eikonal line of fragmentationfunction. Because NRQCD LDMEs should be independent of the direction of the eikonalline of fragmentation function, and all remained IR divergences in full-QCD calculationmust be absorbed into the LDMEs if NRQCD factorization is valid, there should be no IRdivergences in A . On the other hand, if one can show that there is no IR divergence in A , then all IR divergences in full-QCD FF can be absorbed into the NRQCD LDMEs inprinciple, and the NRQCD factorization can be verified at two-loop level.Two-loop verification of NRQCD factorization was also done in Refs. [32–35], in whicheikonal approximation was used to simplify the calculation. The eikonal approximationhowever introduced artificial UV divergences that makes the extraction of IR divergenceshighly nontrivial.In this subsection, we verify the NRQCD factorization without applying eikonal ap-proximation. We take the process g → Q ¯ Q ( P [1] J ) + X with J summed as an example. Toextract the part A only, we use the following projection operator P αα (cid:48) ββ (cid:48) qn = I αα (cid:48) (cid:16) ( D − I ββ (cid:48) − ( D − d ββ (cid:48) (cid:17) , (2.25)where d ββ (cid:48) = − g ββ (cid:48) + P β n β (cid:48) + P β (cid:48) n β P · n − P n β n β (cid:48) ( P · n ) , (2.26)and we get P αα (cid:48) ββ (cid:48) qn ∂∂q β ∂q (cid:48) β (cid:48) A ( q, q (cid:48) ) = I αα (cid:48) ( D −
2) 1 P A . (2.27)Performing the calculations similar to those introduced in Sec. 2, we find that all diver-gences are canceled out after adding real corrections, virtual corrections as well as UVrenormalization contributions. As a result, IR divergences in A ( q, q (cid:48) ) are free of the gaugeline vector n . Thus the NRQCD factorization is verified to hold at two-loop level, which isconsistent with the observation in Refs. [32–35].– 8 – Matching procedure and the SDCs
As SDCs are independent of specifics of the hadronic states, we can match the O ( α s ) PMEsfrom the divergences in full-QCD FFs and obtain the IR-safe SDCs simultaneously, thanksto the fact that NRQCD factorization holds at this order.For g → Q ¯ Q ( P [1 / J ) process, the LO FF can be expressed by D LO [ g → Q ¯ Q ( P [1 / J )] = d LO [ g → Q ¯ Q ( P [1 / J )] (cid:104)O P [1 / J ( P [1 / J ) (cid:105) LO + d LO [ g → Q ¯ Q ( S [8]1 )] (cid:104)O P [1 / J ( S [8]1 ) (cid:105) LO , (3.1)where (cid:104)O P [1] J ( P [1] J ) (cid:105) LO = (cid:104)O P [8] J ( P [8] J ) (cid:105) LO = 1 is normalized. The LO SDC of gluon to S [8]1 in (3.1) has been calculated in Ref. [17, 18], which is given in D − dimension by d LO [ g → Q ¯ Q ( S [8]1 )] = πα s δ (1 − z )( D − N c − m Q , (3.2)and the O ( α s ) PMEs in the MS subtraction scheme can be expressed as (cid:104)O P [1] J ( S [8]1 ) (cid:105) LO = − C F α s πm Q Γ(1 + (cid:15) ) (cid:15) IR (cid:18) πµ r µ (cid:19) (cid:15) , (3.3) (cid:104)O P [8] J ( S [8]1 ) (cid:105) LO = − B F α s πm Q Γ(1 + (cid:15) ) (cid:15) IR (cid:18) πµ r µ (cid:19) (cid:15) , (3.4)where B F = ( N c − / (4 N c ) and µ Λ is the NRQCD factorization scale. With full-QCD FFin Eq. (2.18), the D − dimension SDC d LO [ g → Q ¯ Q ( P [1] J )] can be matched as d LO [ g → Q ¯ Q ( P [1]0 )] = α s Γ(1 + (cid:15) )( D − N c m Q (cid:32) πµ r m Q (cid:33) (cid:15) × (cid:32) (cid:32) −
23 ln (cid:32) µ m Q (cid:33)(cid:33) δ (1 − z )+ 43 z [1 − z ] + + 16 (85 − z ) z + 3(5 − z ) ln(1 − z ) (cid:33) , (3.5)which equals to the finite part of the full-QCD FF in Eq. (2.18) if factorization scale ischosen as µ Λ = 2 m Q . We find that the IR divergence in full-QCD FF D LO [ g → Q ¯ Q ( P [1]0 )] is absorbed into the O ( α s ) PME (cid:104)O P [1] J ( S [8]1 ) (cid:105) LO , and the SDCs are IR-safe up to LO. Thisresult is consistent with that in Ref. [22].Similarly, D NLO [ g → Q ¯ Q ( P [1 / J )] can be expressed in the NRQCD factorization as D NLO [ g → Q ¯ Q ( P [1 / J )] = d NLO [ g → Q ¯ Q ( P [1 / J )] (cid:104)O P [1 / J ( P [1 / J ) (cid:105) LO + d LO [ g → Q ¯ Q ( P [1 / J ))] (cid:104)O P [1 / J ( P [1 / J ) (cid:105) NLO + d NLO [ g → Q ¯ Q ( S [8]1 )] (cid:104)O P [1 / J ( S [8]1 ) (cid:105) LO + d LO [ g → Q ¯ Q ( S [8]1 )] (cid:104)O P [1 / J ( S [8]1 ) (cid:105) NLO , (3.6)where all LO terms are known. For NLO terms, D NLO [ g → Q ¯ Q ( P [1 / J )] are calculated inthis work, (cid:104)O P [1 / J ( P [1 / J ) (cid:105) NLO are well-known, which have only Coulomb divergences and– 9 –anish in dimensional regularization with MS subtraction scheme, and d NLO [ g → Q ¯ Q ( S [8]1 )] has been calculated in Refs. [18, 22], whose D − dimension form is given by d NLO [ g → Q ¯ Q ( S [8]1 )] = α s Γ(1 + (cid:15) )4( D − C F m Q (cid:32) πµ r m Q (cid:33) (cid:15) (cid:34) A ( µ r ) δ (1 − z ) + 1 N c P gg ( z ) × (cid:32) ln (cid:32) µ m Q (cid:33) − (cid:33) + 2(1 − z ) z − − z + z ) (cid:20) ln(1 − z )1 − z (cid:21) + (cid:35) , (3.7)with A ( µ r ) = b N c (cid:32) ln (cid:32) µ r m Q (cid:33) + 133 (cid:33) + 4 N c − π . (3.8)Therefore, SDCs d NLO [ g → Q ¯ Q ( P [1 / J )] and PMEs (cid:104)O P [1 / J ( S [8]1 ) (cid:105) NLO can be extractedfrom Eq. (3.6) thanks to the fact that the former one is finite and the later one is IR-divergent defined in MS subtraction scheme.The O ( α s ) PMEs (cid:104)O P [1 / J ( S [8]1 ) (cid:105) NLO are obtained numerically with high precision.Using PSLQ algorithm [50, 51], analytical expressions are obtained as (cid:104)O P [1] J ( S [8]1 ) (cid:105) NLO = − C F α s π m Q (cid:18) πµ r µ (cid:19) (cid:15) Γ(1 + (cid:15) ) (cid:18) b (cid:15) + 5 n f + (12 π − N c (cid:15) (cid:19) , (3.9) (cid:104)O P [8] J ( S [8]1 ) (cid:105) NLO = − B F α s π m Q (cid:18) πµ r µ (cid:19) (cid:15) Γ(1 + (cid:15) ) (cid:18) b (cid:15) + 5 n f + (3 π − N c (cid:15) (cid:19) . (3.10)Setting factorization scales µ = µ Λ = 2 m Q , the NLO SDC can be expressed as d NLO [ g → Q ¯ Q ( P [1]0 )] = α s πN c m Q (cid:32) p δ δ (1 − z ) + (cid:88) i =0 p i (cid:20) ln i (1 − z )1 − z (cid:21) + + p ( z )+ ln (cid:32) µ r m Q (cid:33) b d LO [ g → Q ¯ Q ( P [1]0 )] (cid:33) , (3.11)where the coefficients p δ , p i can be fitted with the PSLQ algorithm from high-precisionnumerical values as following p δ = 95 − π n f + 756 ζ (3) + 66 π − −
576 ln − N c − − N c ,p = − n f + − π + 119 + 64 ln 218 N c + 83 1 N c ,p = 49 n f − N c ,p = − N c , (3.12)– 10 –nd p ( z ) can be expressed as a piecewise function p ( z ) = − π z N c + (cid:88) i =0 ∞ (cid:88) j =0 ln i z (2 z ) j (cid:32) A fij n f + A ij N c + A Nij N c (cid:33) , for < z < ∞ (cid:88) j =0 (2 z − j (cid:32) B fj n f + B j N c + B Nj N c (cid:33) , for ≤ z ≤ (cid:88) i =0 ∞ (cid:88) j =0 ln i (1 − z ) (2 − z ) j (cid:32) C fij n f + C ij N c + C Nij N c (cid:33) . for < z < . (3.13)The coefficients A kij , B kj , C kij can be evaluated numerically with high precision. In the at-tached ancillary file, we give these coefficients calculated up to j = 500 to obtain the resultswith -digit precision at any value of z .The SDCs of g → Q ¯ Q ( P [1] J ) + X with J = 1 , and g → Q ¯ Q ( P [8] J ) + X with sum J are calculated similarly, with results given as attachments in the ancillary file. To see the effects of the NLO corrections, we choose parameters with m b = 4 .
75 GeV , N c = 3 , n f = 4 , and α s ( µ r = 2 m b ) = 0 . . In Figure 5, we plot the curves of LO SDCsand LO+NLO SDCs for the P -wave FFs with µ = µ Λ = 2 m b . The overall factor m b forcolor-singlet and m b / for color-octet comes from normalization factors at NLO.The sensitivities of LO and LO+NLO FFs with respective to the renormalization scale µ r are illustrated Fig. 5 with bands corresponding to vary µ r from m b to m b . We find thatthere are still large theoretical uncertainties after the NLO corrections.With µ r = µ = µ Λ = 2 m b , we also provide K-factors (the ratio of LO+NLO over LO)for some special values of z in Table 1 , where we find that the K-factors are negative formost values of z .As shown in Figure 5, LO FFs are divergent at z = 1 . These divergences come fromthe δ (1 − z ) and / [1 − z ] + terms in Eq. (3.5). Besides, the NLO FFs are negative anddivergent at both z = 0 and z = 1 . The leading behavior at z → is /z , which meansthat the total fragmenting probability obtained by integrating the NLO FF over z from to is divergent. But fortunately, physical cross sections are not sensitive to the behaviorat z → after convolving FFs with partonic hard parts, which behave as z n in the small z region with n usually larger than . As shown in Eq. (3.11), the divergences at z → behavior as δ (1 − z ) and plus functions with leading power [ln (1 − z ) / (1 − z )] + , whichcomes from the expansion of (1 − z ) − − (cid:15) (cid:15) − .Obviously, the distribution of FF will not give us direct information about physicalcross sections. To see the impact of the NLO results to the physical cross sections, weevaluate the n -th moments of the SDCs, (cid:90) z d z z n ( d LO ( z ) + d NLO ( z )) × c [1 / , (4.1)– 11 – → b ¯ b ( P [1]0 ) + X LOLO + NLO - - z d ( z ) × m P [ ] g → b ¯ b ( P [1]1 ) + X LOLO + NLO - - z d ( z ) × m P [ ] g → b ¯ b ( P [1]2 ) + X LOLO + NLO - - z d ( z ) × m P [ ] g → b ¯ b ( P [8] J ) + X LOLO + NLO - - z d ( z ) × m / P J [ ] Figure 5 . SDCs of the fragmentation functions of g → b ¯ b ( P [1 / J )) + X at LO and NLO. Thedashed line is for d LO ( z ) × c [1 / , and the solid line is for ( d LO ( z ) + d NLO ( z )) × c [1 / , with scalechoices µ r = µ = µ Λ = 2 m b , c [1] = 9 m b for color-singlet and c [8] = 144 m b / for color-octet. Thebands are obtained by varying the renormalization scale µ r by a factor of 2. where c [1] = 9 m b and c [8] = 144 m b / . We plot the numerical results of integration K-factorswith different z and n , which are shown in Figure 6. We find that for the same n , theK-factors are quite steady for different z , meaning that cross sections are dominated byFFs at large z . Besides, the K-factors decrease with the increase of n . Especially, the color-singlet K-factors become negative for large n . And there may be a cancellation betweenthe LO and the NLO FFs at about -th moment for color-singlet quarkonium productions.Further more, it is clear that though the K-factors and values of FFs at each position z = z are quite different for j = 0 , , as shown in Table 1 and in Figure 5, they are almost thesame in the sense of integration as shown in Figure 6.For the special case z = 0 , we can estimate the impact of NLO calculations on physicalcross sections. The numerical results are shown in Table 2. For large n , the influence ofthe behavior at z → becomes more and more prominent and the NLO correction becomesmore and more significant. In addition, the K-factors of color-octet process are much largerthan those of color-singlet process.The above discussion implies that resummation at large z region may be crucial. Thiscan be done within the framework of soft gluon factorization [52], which we leave for future– 12 – K [1]0 K [1]1 K [1]2 K [8] J . − . − . − . − . . − . − . − . − . . − . − . − . − . . − . − . − . − . . − . − . − . − . . − . − . − . − . . − . − . − . − . . − . . − . − . . − . . − . . . − . . − . . . − . . − . . . − . . . . . − . . . . . − . . − . . . − . . − . . . − . − . − . . . − . − . − . − . . − . − . − . − . Table 1 . K-factors at different values of z . K [ c ] j denotes the K-factor of g → Q ¯ Q ( P [ c ] j ) + X andthe symbol J means the sum of j = 0 , , . state SDCs ∗ c [1 / z z z z z P [1]0 LO − . − . − . − . − . LO+NLO − . − . − . . . K-factor . . . − . − . P [1]1 LO − . − . − . − . − . LO+NLO − . − . − . . . K-factor . . . − . − . P [1]2 LO − . − . − . − . − . LO+NLO − . − . − . . . K-factor . . . − . − . P [8] J LO − . − . − . − . − . LO+NLO − . − . − . − . − . K-factor . . . . . Table 2 . Moments and K-factors of SDCs. study. – 13 – → b ¯ b ( P [1]0 ) + X z - K P [ ] g → b ¯ b ( P [1]1 ) + X z - K P [ ] g → b ¯ b ( P [1]2 ) + X z - K P [ ] g → b ¯ b ( P [8] J ) + X z K P J [ ] Figure 6 . The K-factors calculated by the integration results with Eq. (4.1) for different z and n . A Operator renormalization integrals
The operator convolution integration is shown in Eq. (2.22) with the splitting functionsshown in Eq. (2.23). The splitting functions can be divided into three parts, while the LOFF have two regions. Thus there are types of integrals, (cid:90) z d yy δ (1 − y )(1 − zy ) − (cid:15) f (cid:18) zy (cid:19) = λ − (cid:15) f (1 − λ ) , (A.1) (cid:90) z d yy δ (1 − y )(1 − zy ) − − (cid:15) f (cid:18) zy (cid:19) = λ − − (cid:15) f (1 − λ ) , (A.2) (cid:90) z d yy g ( y )(1 − zy ) − (cid:15) f (cid:18) zy (cid:19) = λ − (cid:15) (cid:90) d yh ( y, (cid:15) ) λ − λy g (cid:18) − λ − λy (cid:19) f (1 − λy ) , (A.3)– 14 – z d yy g ( y )(1 − zy ) − − (cid:15) f (cid:18) zy (cid:19) = λ − − (cid:15) (cid:34) − (cid:15) λg (1 − λ ) f (1)+ (cid:90) d yy h ( y, (cid:15) ) (cid:18) λ − λy g (cid:18) − λ − λy (cid:19) f (1 − λy ) − λg (1 − λ ) f (1) (cid:19) (cid:35) , (A.4) (cid:90) z d yy − y ] + (1 − zy ) − (cid:15) f (cid:18) zy (cid:19) = λ − (cid:15) ln λ f (1 − λ ) + λ − (cid:15) (cid:90) d y − y (cid:18) h ( y, (cid:15) ) f (1 − λy ) − − λ − λy f (1 − λ ) (cid:19) , (A.5) (cid:90) z d yy − y ] + (1 − zy ) − − (cid:15) f (cid:18) zy (cid:19) = λ − − (cid:15) ln λ f (1 − λ )+ λ − − (cid:15) (cid:34) (cid:18) − (cid:15) − H ( − (cid:15) ) − ln(1 − λ ) (cid:19) f (1)+ (cid:90) d y − y (cid:18) h ( y, (cid:15) ) f (1 − λy ) − f (1) y − − λ − λy ( f (1 − λ ) − f (1)) (cid:19) (cid:35) , (A.6)where λ = 1 − z , g ( y ) = (1 − y ) /y + y (1 − y ) comes from the splitting function, H ( − (cid:15) ) is the Harmonic number which can be expanded in the limit of (cid:15) → and h ( y, (cid:15) ) = y − (cid:15) = (cid:88) i ( − (cid:15) ) i i ! ln i ( y ) (A.7)is the Taylor series at (cid:15) → . Based on these results, we find that two more behaviors (1 − z ) − (cid:15) ln(1 − z ) and (1 − z ) − − (cid:15) ln(1 − z ) are added to the region at z → . Acknowledgments
We thank Xiao Liu and Hao-Yu Liu for useful discussions. The work is supported inpart by the National Natural Science Foundation of China (Grants No. 11475005, No.11075002, and No. 11875071), the National Key Basic Research Program of China (No.2015CB856700), and High-performance Computing Platform of Peking University.
References [1] G. T. Bodwin, E. Braaten, and G. P. Lepage,
Rigorous QCD analysis of inclusiveannihilation and production of heavy quarkonium , Phys. Rev.
D51 (1995) 1125–1171[ hep-ph/9407339 ] [
InSPIRE ]. [Erratum: Phys. Rev.D55,5853(1997)].[2] J. C. Collins, D. E. Soper, and G. Sterman,
Factorization of Hard Processes in QCD , Adv.Ser.Direct.High Energy Phys. (1988) 1–91 [ hep-ph/0409313 ] [ InSPIRE ].[3] Z.-B. Kang, J.-W. Qiu, and G. Sterman,
Heavy quarkonium production and polarization , Phys.Rev.Lett. (2012) 102002 [ arXiv:1109.1520 ] [
InSPIRE ]. – 15 –
4] Z.-B. Kang, Y.-Q. Ma, J.-W. Qiu, and G. Sterman,
Heavy quarkonium production at colliderenergies: Factorization and Evolution , Phys.Rev.
D90 (2014) 034006 [ arXiv:1401.0923 ][ InSPIRE ].[5] Z.-B. Kang, Y.-Q. Ma, J.-W. Qiu, and G. Sterman,
Heavy Quarkonium Production atCollider Energies: Partonic Cross Section and Polarization , Phys.Rev.
D91 (2015) 014030[ arXiv:1411.2456 ] [
InSPIRE ].[6] S. Fleming, A. K. Leibovich, T. Mehen, and I. Z. Rothstein,
The Systematics of QuarkoniumProduction at the LHC and Double Parton Fragmentation , Phys.Rev.
D86 (2012) 094012[ arXiv:1207.2578 ] [
InSPIRE ].[7] S. Fleming, A. K. Leibovich, T. Mehen, and I. Z. Rothstein,
Anomalous dimensions of thedouble parton fragmentation functions , Phys.Rev.
D87 (2013) 074022 [ arXiv:1301.3822 ][ InSPIRE ].[8] V. Gribov and L. Lipatov,
Deep inelastic e p scattering in perturbation theory , Sov.J.Nucl.Phys. (1972) 438–450 [ InSPIRE ].[9] G. Altarelli and G. Parisi,
Asymptotic Freedom in Parton Language , Nucl.Phys.
B126 (1977)298 [
InSPIRE ].[10] Y. L. Dokshitzer,
Calculation of the Structure Functions for Deep Inelastic Scattering and e + e − Annihilation by Perturbation Theory in Quantum Chromodynamics. , Sov.Phys.JETP (1977) 641–653 [ InSPIRE ].[11] E. Braaten, K. Cheung, and T. C. Yuan, Z decay into charmonium via charm quarkfragmentation , Phys.Rev.
D48 (1993) 4230–4235 [ hep-ph/9302307 ] [
InSPIRE ].[12] E. Braaten and T. C. Yuan,
Gluon fragmentation into heavy quarkonium , Phys.Rev.Lett. (1993) 1673–1676 [ hep-ph/9303205 ] [ InSPIRE ].[13] P. L. Cho, M. B. Wise, and S. P. Trivedi,
Gluon fragmentation into polarized charmonium , Phys. Rev.
D51 (1995) R2039–R2043 [ hep-ph/9408352 ] [
InSPIRE ].[14] E. Braaten and T. C. Yuan,
Gluon fragmentation into P wave heavy quarkonium , Phys.Rev.
D50 (1994) 3176–3180 [ hep-ph/9403401 ] [
InSPIRE ].[15] M. Beneke and I. Z. Rothstein,
Psi-prime polarization as a test of color octet quarkoniumproduction , Phys. Lett.
B372 (1996) 157–164 [ hep-ph/9509375 ] [
InSPIRE ]. [Erratum: Phys.Lett.B389,769(1996)].[16] J. P. Ma,
Quark fragmentation into p wave triplet quarkonium , Phys. Rev.
D53 (1996)1185–1190 [ hep-ph/9504263 ] [
InSPIRE ].[17] E. Braaten and Y.-Q. Chen,
Dimensional regularization in quarkonium calculations , Phys.Rev.
D55 (1997) 2693–2707 [ hep-ph/9610401 ] [
InSPIRE ].[18] E. Braaten and J. Lee,
Next-to-leading order calculation of the color octet 3S(1) gluonfragmentation function for heavy quarkonium , Nucl. Phys.
B586 (2000) 427–439[ hep-ph/0004228 ] [
InSPIRE ].[19] G. Hao, Y. Zuo, and C.-F. Qiao,
The Fragmentation Function of Gluon Splitting into P-waveSpin-singlet Heavy Quarkonium , [ arXiv:0911.5539 ] [
InSPIRE ].[20] Y. Jia, W.-L. Sang, and J. Xu,
Inclusive h c Production at B Factories , Phys. Rev.
D86 (2012) 074023 [ arXiv:1206.5785 ] [
InSPIRE ]. – 16 –
21] G. T. Bodwin, H. S. Chung, U.-R. Kim, and J. Lee,
Quark fragmentation into spin-triplet S -wave quarkonium , Phys. Rev.
D91 (2015) 074013 [ arXiv:1412.7106 ] [
InSPIRE ].[22] Y.-Q. Ma, J.-W. Qiu, and H. Zhang,
Heavy quarkonium fragmentation functions from aheavy quark pair. I. S wave , Phys.Rev.
D89 (2014) 094029 [ arXiv:1311.7078 ] [
InSPIRE ].[23] Y.-Q. Ma, J.-W. Qiu, and H. Zhang,
Fragmentation functions of polarized heavy quarkonium , JHEP (2015) 021 [ arXiv:1501.04556 ] [ InSPIRE ].[24] P. Zhang, Y.-Q. Ma, Q. Chen, and K.-T. Chao,
Analytical calculation for the gluonfragmentation into spin-triplet S-wave quarkonium , Phys. Rev.
D96 (2017) 094016[ arXiv:1708.01129 ] [
InSPIRE ].[25] E. Braaten and T. C. Yuan,
Gluon fragmentation into spin triplet S wave quarkonium , Phys.Rev.
D52 (1995) 6627–6629 [ hep-ph/9507398 ] [
InSPIRE ].[26] G. T. Bodwin and J. Lee,
Relativistic corrections to gluon fragmentation into spin triplet Swave quarkonium , Phys. Rev.
D69 (2004) 054003 [ hep-ph/0308016 ] [
InSPIRE ].[27] G. T. Bodwin, U.-R. Kim, and J. Lee,
Higher-order relativistic corrections to gluonfragmentation into spin-triplet S-wave quarkonium , JHEP (2012) 020[ arXiv:1208.5301 ] [
InSPIRE ].[28] Q.-F. Sun, Y. Jia, X. Liu, and R. Zhu,
Inclusive h c production and energy spectrum from e + e − annihilation at a super B factory , Phys. Rev.
D98 (2018) 014039 [ arXiv:1801.10137 ][ InSPIRE ].[29] P. Zhang, C.-Y. Wang, X. Liu, Y.-Q. Ma, C. Meng, and K.-T. Chao,
Semi-analyticalcalculation of gluon fragmentation into S [1 , quarkonia at next-to-leading order , JHEP (2019) 116 [ arXiv:1810.07656 ] [ InSPIRE ].[30] P. Artoisenet and E. Braaten,
Gluon fragmentation into quarkonium at next-to-leading orderusing FKS subtraction , JHEP (2019) 227 [ arXiv:1810.02448 ] [ InSPIRE ].[31] F. Feng and Y. Jia,
Next-to-leading-order QCD corrections to gluon fragmentation into S (1 , quarkonia , [ arXiv:1810.04138 ] [ InSPIRE ].[32] G. C. Nayak, J.-W. Qiu, and G. F. Sterman,
Fragmentation, factorization and infrared polesin heavy quarkonium production , Phys. Lett. B (2005) 45–51 [ hep-ph/0501235 ].[33] G. C. Nayak, J.-W. Qiu, and G. Sterman,
Fragmentation, NRQCD and NNLO factorizationanalysis in heavy quarkonium production , Phys.Rev.
D72 (2005) 114012 [ hep-ph/0509021 ][ InSPIRE ].[34] G. C. Nayak, J.-W. Qiu, and G. F. Sterman,
NRQCD Factorization and Velocity-dependenceof NNLO Poles in Heavy Quarkonium Production , Phys. Rev. D (2006) 074007[ hep-ph/0608066 ].[35] G. T. Bodwin, H. S. Chung, J.-H. Ee, U.-R. Kim, and J. Lee, Covariant calculation of atwo-loop test of nonrelativistic QCD factorization , Phys. Rev. D (2020) 096011[ arXiv:1910.05497 ].[36] J. C. Collins and D. E. Soper,
Parton Distribution and Decay Functions , Nucl.Phys.
B194 (1982) 445 [
InSPIRE ].[37] A. V. Kotikov,
Differential equations method: New technique for massive Feynman diagramscalculation , Phys. Lett.
B254 (1991) 158–164 [
InSPIRE ]. – 17 –
38] Z. Bern, L. J. Dixon, and D. A. Kosower,
Dimensionally regulated one loop integrals , Phys.Lett.
B302 (1993) 299–308 [ hep-ph/9212308 ] [
InSPIRE ].[39] E. Remiddi,
Differential equations for Feynman graph amplitudes , Nuovo Cim.
A110 (1997)1435–1452 [ hep-th/9711188 ] [
InSPIRE ].[40] T. Gehrmann and E. Remiddi,
Differential equations for two loop four point functions , Nucl.Phys.
B580 (2000) 485–518 [ hep-ph/9912329 ] [
InSPIRE ].[41] J. M. Henn,
Multiloop integrals in dimensional regularization made simple , Phys. Rev. Lett. (2013) 251601 [ arXiv:1304.1806 ] [
InSPIRE ].[42] R. N. Lee,
Reducing differential equations for multiloop master integrals , JHEP (2015)108 [ arXiv:1411.0911 ] [ InSPIRE ].[43] L. Adams, E. Chaubey, and S. Weinzierl,
Simplifying Differential Equations for MultiscaleFeynman Integrals beyond Multiple Polylogarithms , Phys. Rev. Lett. (2017) 141602[ arXiv:1702.04279 ] [
InSPIRE ].[44] M. Caffo, H. Czyz, M. Gunia, and E. Remiddi,
BOKASUN: A Fast and precise numericalprogram to calculate the Master Integrals of the two-loop sunrise diagrams , Comput. Phys.Commun. (2009) 427–430 [ arXiv:0807.1959 ] [
InSPIRE ].[45] M. Czakon,
Tops from Light Quarks: Full Mass Dependence at Two-Loops in QCD , Phys.Lett.
B664 (2008) 307–314 [ arXiv:0803.1400 ] [
InSPIRE ].[46] R. Mueller and D. G. Öztürk,
On the computation of finite bottom-quark mass effects inHiggs boson production , JHEP (2016) 055 [ arXiv:1512.08570 ] [ InSPIRE ].[47] R. N. Lee, A. V. Smirnov, and V. A. Smirnov,
Solving differential equations for Feynmanintegrals by expansions near singular points , JHEP (2018) 008 [ arXiv:1709.07525 ].[48] X. Liu, Y.-Q. Ma, and C.-Y. Wang, A Systematic and Efficient Method to ComputeMulti-loop Master Integrals , Phys. Lett.
B779 (2018) 353–357 [ arXiv:1711.09572 ][ InSPIRE ].[49] X. Liu, Y.-Q. Ma, W. Tao, and P. Zhang,
Calculation of Feynman loop integration andphase-space integration via auxiliary mass flow , [ arXiv:2009.07987 ].[50] H. R. P. Ferguson, D. H. Bailey, and S. Arno,
Analysis of PSLQ, an integer relation findingalgorithm , Mathematics of Computation (Jan., 1999) 351–369.[51] D. H. Bailey and D. J. Broadhurst, Parallel integer relation detection: Techniques andapplications , Math. Comput. (2001) 1719–1736 [ math/9905048 ].[52] Y.-Q. Ma and K.-T. Chao, New factorization theory for heavy quarkonium production anddecay , Phys. Rev. D (2019) 094007 [ arXiv:1703.08402 ].].