Gluon TMD fragmentation function into quarkonium

We compute the gluon transverse-momentum-dependent fragmentation function (TMDFF) at next-to-leading order (NLO) into heavy quarkonium in the color-octet $^3S_1^{[8]}$ channel, based on the NRQCD factorization approach. The spurious rapidity divergences are explicitly shown to cancel in a well-defined TMDFF, which incorporates the needed soft factor. We also compute the integrated gluon FF at NLO in the same $^3S_1^{[8]}$ channel, and show that the matching coefficient of the TMDFF onto the FF at large transverse momentum is the expected one. These results are relevant to perform precise and sensible phenomenological studies of transverse-momentum spectra of quarkonium production, for which the production mechanism through fragmentation plays a relevant role, like in the future Electron-Ion Collider.

In more recent years, with the formulation of soft collinear effective theory (SCET) [32][33][34][35] new details have been included in this description, providing a robust TMD factorization theorem for the production of quarkonium at small transverse momentum right from the hard reaction, which is given in terms of the so-called and newly introduced TMD shape functions [16,17].
Despite the abundant interest in quarkonium TMDs, little has been done so far in the direction of TMD quarkonium fragmentation processes, by which we mean single parton fragmentation mechanism.In processes like ep → J/ψ + jet, e + e − → J/ψ + jet or e + e − → (J/ψ)(J/ψ) and even more energetic ones where the J/ψ is traded for an Υ, the quarkonium states can also be produced via the fragmentation of a single parton.This mechanism becomes relevant when a hard scale, much larger than the quarkonium mass, exists.Such processes can certainly be observed at an Electron Ion Collider (EIC) as the ones expected to be built in several locations [36,37].
While the light-quark TMD fragmentation function to quarkonia has been considered in [38] (see also the recent works [39,40] for polarized TMDFFs), in this work we concentrate on the gluon TMD fragmentation function to quarkonia.
In our analysis we employ the NRQCD factorization conjecture, where the quarkonium state is produced at large distances through the hadronization of a heavy quark-antiquark pair, Q Q(n).The pair can be found in any color and angular configuration n = 2S+1 L [col.] J , but then the probability that the pair decays in the colorless quarkonium state scales with the relative velocity, v, of the quark-antiquark pair in the quarkonium rest frame.We decompose the TMDFF within NRQCD in terms of calculable short-distance matching coefficients and long-distance matrix elements (LDMEs).We proceed then to the NLO calculation of this function, extracting the matching coefficient onto the corresponding LDMEs in the region q T ∼ M .Using this information we can evaluate the contribution to the cross section from gluon fragmentation.We also compute the matching of the TMDFF onto the corresponding integrated FF at q T > M .
The so-called short-distance coefficient in NRQCD for a gluon fragmenting into a heavy-quark pair in the 3 S [8] 1 channel has been calculated up to next-to-leading order (NLO) in the strong coupling in [41][42][43], with some discrepancies in the finite terms.We have checked this result as a by-product of our analysis, and we do agree with the result obtained in [43].
We regularize the rapidity divergences using a δ-regulator whenever necessary (see [44,45] for more details) and we check explicitly the cancellation of the rapidity divergences with the needed soft factor, which is the same as the one appearing e.g. in Higgs boson production at small transverse momentum [4,46].This paper is organized as follows.In Sec. 2 we start by setting the notation and the relevant definitions.Then in Sec. 3 we present the main results, for which more technical details can be found in the Appendix.Finally in Sec. 4 we conclude.

Notation
In this Section we establish the notation that we are going to use in the calculations, as well as the general definition of the gluon TMD fragmentation function and its matching onto the LDMEs.

Kinematics
We use the coordinates of the light cone set by the following scalar products in the spacetime dimensions d = 4 − 2ε: with n 2 = n2 = 0, (nn) = 1.We denote the momentum of the fragmenting gluon by q and the momentum of the heavy-quark Q Q pair by P .Also, we work in the frame in which the transverse momentum of P vanishes.

Gluon TMDFF definition
The operator for the gluon fragmentation function follows from [44].The unsubtracted TMDFF is the hadronix matrix element of that operator , with iD µ n⊥ = ∂ µ n⊥ + g s A µ n⊥ where A µ n⊥ is the SCET n-collinear field.In (2.5), W n is the collinear Wilson lines defined as The renormalized gluon TMD is defined as follows where Z g is the usual renormalization factor for UV divergences and R g is the rapidity renormalization factor, µ is the scale of UV subtraction and ζ is the scale of rapidity subtraction.Here, R g is the following expression which describes the ratio between the soft function denoted as S(b ⊥ ) and the soft overlap of the collinear and soft sectors through the term Z b [47], denoting the zero-bin contribution.The soft function is defined as a expectation value of soft Wilson lines: where the soft Wilson lines are defined as These soft Wilson lines, with calligraphic typography, are in the adjoint representation, where the color generators are given by (t a ) bc = −if abc .

Gluon TMDFF factorization
We employ, for q T ∼ M , the NRQCD formalism [28] to write the gluon TMDFF as a product of short distance coefficients and the long distance matrix elements (LDMEs): Here, n = L 2S+1 [col.]J describes the color and angular momentum configuration of the heavy-quark pair.All relativistic effects are absorbed in d g→Q Q(n) , which can be calculated as a perturbative series in the strong coupling constant α s through matching, and the LDMEs are defined as follows where a J/ψ and a † J/ψ are the operators of annihilation and creation of the state describing the J/ψ, K n and K ′ n are products of a color matrix, a spin matrix and other fields, and χ and ψ are the field operators for the heavy quarks in NRQCD.

Results
In the present section we show the results of the calculation of the gluon TMDFF for quarkonium at NLO, defined in (2.7).The details of the calculation can be found in the Appendix.We have used the δ-regularization [46] in order to regularize the rapidity divergences, defined at operator level as follows: and similarly for the rest of the collinear and soft Wilson lines.

Leading Order
On the one hand, we calculate the left hand side of the equation (2.11), i.e. the matrix element (2.4), around the threshold (q = 0): This result is already the TMDFF, since the soft function at LO is just 1.The spinorial structure which we obtain describes the configuration n = S 3 [8] 1 .
On the other hand, we know from pNRQCD that the LDME describing this configuration is the following, Therefore, by matching both sides of the equation (2.11) the SDC is (3.4)

Next to Leading Order
The virtual contributions to the gluon TMDFF at NLO are shown in figure 1.
) where δ ≡ δ + /P + and β 0 = 11C A /3 − 2n f /3.In order to obtain a well-defined hadronic quantity it is necessary to renormalize the divergences according to the equations (2.7) and (2.8).We have used the δ-regularization where the subtractions related with Z b are equal to the soft function [44]: The virtual contribution to the SF at one loop is as follows [4,46] S vir.(δ At the end, after the renormalization of the rapidity divergences, the virtual contribution of the SDC at NLO is Since real diagrams will not contain UV divergences, because the transverse momentum (or distance) is finite, we can extract the evolution of the TMDFF solely from its virtual part.The renormalization of the TMDFF at NLO requires the renormalization of both α s , which appears already at LO, and the operator itself.In the MS-scheme it is well-known that the coupling is renormalized as follows: And combining this result with the TMDFF at LO in (3.4) and the NLO result in (3.8) we obtain the complete renormalization factor for the UV divergences: The anomalous dimension of the TMDFF [48], which gives the evolution in the renormalization scale µ to O(α s ), is the following This result obviously agrees with the one that can be found e.g. in [44], since the UV behavior of the operator does not depend on the nature of the hadronic state.
We now turn to the real contribution to the gluon TMDFF at NLO, coming from the diagrams shown in figure 2. We separate the diagrams into two groups.Those with rapidity divergences, i.e. c, d1 and d2, and those without, i.e. a, b, e1, e2, f 1 and f 2.
The contribution of diagrams c, d1 and d2 is The contribution of diagrams a, b, e and f is where B(z, b T ) is defined in the Appendix.
Adding the results (3.11) and (3.12) we obtain the real contribution before renormalizing the rapidity divergences: where The real contribution of the Soft Function, which can be found e.g. in [4,46], is Therefore, the real contribution to the SDC at NLO after the renormalization of the rapidity divergences is the following: Finally, by adding the virtual part in (3.8) and the real part in (3.16), we get the short-distance matching coefficient for the gluon TMDFF at NLO into heavy quarkonium in the color-octet S 3 [8] 1 channel: (3.17) We have obtained a well-defined quantity, since it is free from rapidity divergences and is well-behaved in the limit z → 1 (in terms of distributions), remaining only the UV divergences which are removed by standard renormalization.
We end this section by focusing on the ζ-evolution of the gluon TMDFF.The TMD evolution equation in ζ [49] is the following: where D g (b T ; µ) is called the rapidity anomalous dimension (RAD), also called Collins-Soper (CS) kernel.Therefore, from our results, the RAD which gives the evolution in the rapidity scale with L T (b T ; µ) = ln(µ 2 b 2 T e 2γ E /4).This result, es expected, is the same as for the TMDPDF [4], since being the TMD operator the same, the structure of the rapidity divergences is the same.
Finally, the evolution in the plane (µ, ζ) of the gluon TMDFF is ) with the UV anomalous dimension γ F (µ, ζ) and the RAD D g (b T ; µ) obtained at O(α s ) in (3.11) and (3.19), respectively.In the integral, P denotes any path connecting the points (µ f , ζ f ) and (µ i , ζ i ).

Matching onto integrated FF
Following the notation of [44], the small-b T matching between the gluon TMDFF and its corresponding integrated function is described by the OPE of the gluon TMDFF onto the standard FF: where ⊗ is the Mellin convolution in variable z, and both hadronic matrix elements are understood to be renormalized.All the dependence on the transverse coordinate b T and rapidity scale is in the OPE Wilson coefficient.The integrated FF is defined as Notice the different prefactor of z as compared to the TMDFF (apart from the obvious difference in the separation of the fields of the operator, which in this case is just in the collinear direction).In Section 5.3 we have calculated the collinear unpolarized FF in order to obtain the OPE Wilson coefficient of the perturbative expansion of the unpolarized TMDFF at large transverse momentum: With this result we have extracted the matching coefficient in the g/g channel at NLO: Here, notice that we have expanded in b T M the functions in the TMDFF in (3.17), in order to obtain the Wilson matching coefficient in (3.24) via equation (3.21).The obtained coefficient is the same as in [44], as expected, since the matching is done at the operator level, regardless of the hadronic state.

Conclusions
The calculation of the 3 S [8] 1 gluon TMD fragmentation function at NLO that we have performed in this paper shows how the TMDFF behaves at the threshold of quarkonium production.For this we have combined both TMD formalism and NRQCD factorization at leading power.As expected, the rapidity divergences typical of TMD functions are not affected by the threshold, because the light-cone and heavy-mass regimes do not interfere with each other.This means that one can use all our knowledge about the UV and rapidity evolution kernel also in the case of quarkonia production.We recall that the rapidity evolution kernel has been recently extracted from experiment with a N 4 LL analysis [50].The results of this paper provide necessary ingredients to perform phenomenological studies of transverse-momentum spectra of quarkonia production, for which the fragmentation production mechanism is important.Real diagrams contributing to NLO.Only diagrams c, d1 and d2 have rapidity divergences.Hermitian conjugates of diagrams c, d1, d2, e1, e2, f1 y f2 are not shown.

Appendix
In this section we present the details for the calculation of the short-distance matching coefficient of the unpolarized gluon TMDFF at NLO into heavy quarkonium in the color-octet 3 S 1 channel, which is shown in (3.17).We have used dimensional regularization for ultraviolet (UV) divergences in the MS-scheme, i.e. µ 2 → µ 2 e γ E /4π, and δ-regularization for infrared (IR) and rapidity divergences.Our goals are to obtain the short-distance coefficient d NLO g→J/ψ around the threshold by matching both sides of (2.11), to show the cancellation of the spurious rapidity divergences and to extract the Wilson matching coefficients of the TMDFF onto its integrated function, as shown in (3.21).
With δ-regulator the relevant Feynman rules for the collinear gluon field in (2.5), needed in the calculation at NLO, become where the superscripts (0) and (1) denote the power of the strong coupling constant g.The next order is zero.We assume that the relative momentum of the heavy-quark pair is small compared with the mass m c , so The Dirac spinors for the c an c in the frame in which the pair has total momentum P are where ξ and η are the Pauli spinors and they are normalized such that ξ † ξ = 1 and η † η = 1.
First, we calculate the left hand side of (2.11), that is the matrix element defined in (2.4), around the threshold (q = 0).To do this, we need the following nonrelativistic expansions in powers of q, which we will obtain when we make the calculation of equation (2.11) and which will allow us to compare with the LDMEs when matching: where σ denotes the Pauli matrices and L µ i is a boost matrix defined as (5.4) For example, equation (2.4) at Leading Order (LO), which is shown in fig. 1, is (5.5) By using the expressions of (5.3) we get We can average that factor over rotations if we also average the projection operator on the right side of (2.11).The average of the spin factor is so, by using the second relation in (5.4) we get where we have defined (5.9) Therefore, the left hand side of (2.11) at LO is After calculating that side of (2.11) we note that we obtain the spin factor ξ † T a σ k η × η † T a σ k ξ that defines the color and angular momentum configuration n = S 3 [8] 1 .We now focus on the NRQCD side of the matching equation.The spinor structure ξ † T a σ k η × η † T a σ k ξ can be defined as the expansion to LO in α s of the following NRQCD matrix element: and, in turn, (5.12) Finally, by matching between the two sides of equation (2.11) we can obtain the shortdistance matching coefficient.We can conclude that by following the equations (5.11) and (5.12), the short-distance matching coefficient will be the result of the left-hand side without the spin factor ξ † T a σ k η × η † T a σ k ξ and with a 1/m c factor: (5.13)

Details for the calculation of the virtual contribution
Diagrams contributing to the gluon TMDFF at NLO are shown in figures 1 and 2. In this section we calculate the diagrams shown in fig. 1 which are the so-called virtual diagrams.The diagram 1a is the propagator correction for a gluon with invariant mass M 2 = 4 m 2 c , diagrams 1b1 and b2 describe the vertex correction factor and diagrams 1c1 and c2 are the wave renormalization factor (WFR) for a heavy quark and anti-quark: ∆ 1a+h.c.g→J/ψ = 2 Re n f ∆ 1a,q g→J/ψ + ∆ 1a,g g→J/ψ + ∆ 1a,G In order to calculate the contribution of the renormalization of the heavy-quark wave function we need to define the following: (5.21) The renormalization factor of the heavy-quark wave function is then (5.22) The UV pole comes from B(m 2 c ) and the IR pole from the other term.
Diagrams 1(d1+d2) and its hermitian conjugate give (5.24) Here, we have a different spinorial structures from the previous one: where we have used that the rotational average of Diagrams 1e and its hermitian conjugate give where V αβγ (p, q, k) = g αγ (k − p) β + g βγ (q − k) α + g αβ (p − q) and we have used eq.(5.8).
The integrals that appear in the calculation of the virtual contribution are as follows We note that we obtain the same spin factor ξ † T a σ k η × η † T a σ k ξ from all the virtual diagrams, which in turn, is the same as the one we obtain in eq.( 5.10) at LO.As we said in the calculation of the LO contribution, this spin factor defines the configuration n = S 3 [8] 1 .
In order to obtain the SDC of each diagram, matching translates into removing the spin factor and adding the factor 1/m c to the previous results: (5.28) Here 1b denotes 1b1 + 1b2 and so on.

Details for the calculation of the real contribution
In this section we calculate the contributon of the diagrams shown in fig. 2 which are the so-called real diagrams.Diagram 2a gives (5.29)In the second equality, we have used that where L µν is defined in eq.(5.9).Diagram 2b gives (5.31) Diagrams 2c and its Hermitian conjugate give (5.32) Diagrams 2d and its Hermitian conjugate give (5.33) Here, we have a different spinorial structures from the previous ones: ) where we have used that the rotational average of ξ † σ i η × η † {[σ j , σ k ], σ l }ξ vanishes.Diagrams 2e1 + 2e2 and its Hermitian conjugate give Diagrams 2f1 + 2f2 and its Hermitian conjugate give (5.36)Where we have used that (5.37) Finally, we note that we obtain the same structure ξ † T a σ k η ×η † T a σ k ξ for all diagrams.As we mentioned in the calculation at LO, this spin factor defines the configuration n = 3 S [8] 1 .In order to obtain the SDCs from each diagram, matching translates into removing the spin factor and adding the factor 1/m c to the previous results.

Fourier transform
We show the calculation of the Fourier transform of the results of the previous section.We separate the calculation into the diagrams that have rapidity divergences and those that do not.
The real contribution of diagrams c and d as a function of k ⊥ is (5.38) We need the following integrals in order to obtain that contribution as a function of b ⊥ , where B = b 2 ⊥ /4 and a = 2 BM 2 /z 2 .Therefore, the Fourier transform of the equation (5.38) is the following (5.41) The real contribution of diagrams a, b, e and f as a function of k ⊥ is (5.42)We need the following integrals in order to obtain that contribution as a function of b ⊥ , (5.45)

Behaviour at z → 1
In this section we make explicit the infrared divergences associated with the limit z → 1 for the real diragams.We want to extract those divergences and reexpress them as poles of ε and a combination of plus distributions of z.
We know that the limit form when x → 0 of the modified Bessel function of second kind of order n, K n (x), is so K n (x) diverges as a logarithm when n = 0 and as 1/x n when n > 0. We can extract this behaviour from K n (x) as follows where f (x) and g n (x) are regular at x → 0. On the other hand, (5.48) Since the order O(ε) of K ε (x) is zero we can define (5.49) with f (x) and g 1+ε (x) regular at x → 0. If we proceed in the same way in eq. ( 5.41) we obtain the following that are regular at z → 1.And for the eq.(5.45): ( (5.56) which is regular at z → 1.
Note that after substituting the expansions of this section in equations (5.41) and (5.45) we obtain the final result of equations (3.11) and (3.12).

Integrated FF
In this section we show the details for the calculation of the equation (3.23).The virtual contribution for the integrated function is the same as for the TMDFF in eq.(3.5): (5.57)In the case of the real contribution, we start from equations (5.38) and (5.42), and instead of computing the Fourier transform to obtain the TMDFF, we integrate over k ⊥ .
For the contribution of the diagrams 2(c + d1 + d2) + h.c.we need the following integrals  (5.60) In that equation there is an infrared divergence associated with the limit z → 1 which can be made explicit by using the expansion For the contribution of the diagrams 2(a + b + e1 + e2 + f1 +f2) + h.c.we need the following integrals  (5.65) In that equation there is an infrared divergence associated with the limit z → 1 which can be made explicit by using the expansion (5.68) and b T is the conjugate variable to the transverse momentum.The color normalization factor and the normalization of the number of physical gluon polarizations in d = 4 − 2ε are encoded in the prefactor.In this equation, B µ n⊥ is the gluon field strength defined as follows ) where L T = ln µ 2 b 2 T e 2γ E /4 and A 1 (z, b T ) and A 2 (z, b T ) are defined in the Appendix.

Figure 1 .
Figure 1.Diagram contributing to LO and virtual diagrams contributing to NLO.Only diagrams d1, d2 and e have rapidity divergences.Hermitian conjugates of diagrams a, b, c, d1, d2 and e are not shown.

) where B = b 2 ⊥ / 4
and a = 2 BM 2 /z 2 .Therefore the Fourier transform of the equation (5.42) is the following d a,b,e,f g→J/ψ (z, b ⊥