$J/\psi$ meson production in SIDIS: matching high and low transverse momentum

We consider the transverse momentum spectrum and the $\cos 2\phi$ azimuthal distribution of $J/\psi$ mesons produced in semi-inclusive, deep-inelastic electron-proton scattering, where the electron and the proton are unpolarized. At low transverse momentum, we propose factorized expressions in terms of transverse momentum dependent gluon distributions and shape functions. We show that our formulae correctly match with the collinear factorization results at high transverse momentum, obtained in the framework of collinear, nonrelativistic QCD (NRQCD). We conclude that the shape functions at high transverse momentum are independent of the quantum numbers of the produced quarkonium state, except for the overall magnitude given by the appropriate NRQCD long distance matrix element.


I. INTRODUCTION
The production of a light hadron h with a specific transverse momentum in semi-inclusive, deep-inelastic electronproton scattering (SIDIS), e p → e h X, is in general characterized by three different scales: the hard scale of the process Q, given by the virtuality of the gauge boson exchanged in the reaction, the nonperturbative QCD scale Λ QCD , and the magnitude of the hadron transverse momentum q T in a suitable reference frame. Depending on the value of q T , two different factorization frameworks can be adopted for the description of this process. Both of them enable to separate the short-distance from the long-distance contributions to the cross section. While the former can be perturbatively calculated through a systematic expansion in the strong coupling constant, the latter has to be parametrized in terms of parton distributions (PDFs) and fragmentation functions (FFs), which need to be extracted from data.
More explicitly, collinear factorization is applicable in the so-called high-q T region, namely for q T Λ QCD , where the transverse momentum in the final state is generated by perturbative radiation and the cross section is expressed in terms of collinear (i.e., integrated over transverse momentum) PDFs and FFs. The other framework, based on transverse momentum dependent (TMD) factorization [1][2][3] is valid at low q T , q T Q, and involves TMD PDFs and FFs (or TMDs for short). The high-and low-q T regions overlap for Λ QCD q T Q, where both descriptions can therefore be applied. If the two results describe the same dynamics, characterized by the same power behavior, they have to match in this intermediate region. Conversely, if the two results describe competing mechanisms, they should be considered independently and added together [4].
The SIDIS cross section differential in q T and integrated over the azimuthal angle of the final hadron can be expressed in terms of unpolarized, twist-two TMDs in the small-q T region, and its matching with the collinear description has been demonstrated in Ref. [4]. The matching for the analogous observable in Drell-Yan (DY) dilepton production, p p → X, has been proven as well in Refs. [5,6]. Azimuthal asymmetries in both SIDIS [7] and DY [8][9][10] processes have also been widely investigated within the TMD framework. In particular, the matching of the cos φ modulations, which are suppressed by a factor q T /Q with respect to the φ-integrated cross sections and involve twist-three TMDs, has been shown very recently in Ref. [11].
In this paper, we analyze J/ψ production in SIDIS, e p → e J/ψ X, along the same lines of Refs. [4,11]. There are several reasons that motivate our study. First of all, since its discovery in 1974, the J/ψ meson, a charm-anticharm quarkonium bound state with odd charge parity, has always attracted a lot of attention as a probe of the perturbative and nonperturbative aspects of quantum chromodynamics (QCD) and their interplay. Moreover, J/ψ production in both ep [12][13][14][15] and pp collisions [12,[16][17][18][19] has been proposed lately as a tool to access gluon TMDs. Similar studies have also been carried out within the so-called generalized parton model approach [20][21][22][23][24][25][26][27]. From the experimental point of view, these reactions should have a very clean signature due to the large branching ratio of the J/ψ leptonic decay modes. Hence our findings could be in principle verified at the future Electron-Ion Collider (EIC) planned in the United States [28,29].
As compared to e p → e h X, the study of e p → e J/ψ X presents an additional complication, namely a second hard scale given by the J/ψ mass, M ψ . Since we want to avoid contributions from photoproduction processes, we focus on the kinematic region where Q M ψ . In principle, each of the two scales, or any combination of them, can be chosen as the factorization scale µ in the calculation of the cross section. From our analysis it will turn out that the choice µ = Q 2 + M 2 ψ allows for a smooth transition of the cross section from the high-to low-q T region. When the J/ψ meson is produced with a large transverse momentum, q T M ψ , nonrelativistic QCD (NRQCD) can be applied. This rigorous theoretical framework implies a separation of short-distance coefficients, which can be calculated perturbatively as expansions in the strong-coupling constant α s , from long-distance matrix elements (LDMEs), which must be extracted from experiment [30]. The relative importance of the latter can be estimated by means of velocity scaling rules, i.e. the LDMEs are predicted to scale with a definite power of the heavy-quark velocity v in the limit v 1. In this way, the theoretical predictions are organized as double expansions in α s and v. The main feature of this formalism is that the charm-anticharm quark pair forming the bound state can be produced both in a color-singlet (CS) configuration, with the same quantum numbers as the J/ψ meson, and as an intermediate color-octet (CO) state with different quantum numbers. In the latter case, the pair subsequently evolves into a colorless state through the emission of soft gluons.
In the small-q T region, q T M ψ , TMD factorization has not yet been proven in a rigorous way for the process e p → e J/ψ X. There are however strong arguments in favor of its validity, if we consider the analogy with e p → e h X, for which TMD factorization holds at all orders [1]. The only difference from the color point of view is that the dominant partonic subprocess is now γ * g → cc instead of γ * q → q . Hence final state interactions will be resummed in the gauge link of the gluon correlator, which will be in the adjoint representation, rather than in a quark correlator in the fundamental representation. Since the J/ψ mass does not affect the gauge link structure, we do not expect any TMD factorization breaking effects due to color entanglement [13].
Motivated by these arguments, in the present analysis we put forward factorization formulae, valid at the twisttwo level, for the transverse momentum spectrum and the cos 2φ azimuthal distribution of J/ψ mesons produced in SIDIS. In addition to the usual TMD PDFs, we consider the recently proposed shape functions [31,32], which are the generalization of the collinear LDMEs in NRQCD. Alternatively, they can be seen as the analog of the TMD FFs for light hadron production in SIDIS. By requiring a proper matching with the collinear results, in complete analogy with the TMD cross sections for e p → e h X and p p → X, we are able to assess the role of the shape functions in the TMD formalism for quarkonium production. This will have important implications for a recent suggestion to extract poorly known CO LDMEs from a comparison between quarkonium production and open heavy quark pair production in SIDIS at the EIC [13].
The paper is organized as follows. In Section II we define the variables and the reference frames that are adopted in our calculation. Parametrizing the cross section in terms of different structure functions, we compute the cross section in the collinear framework after which we take the small-q T limit. Section III is devoted to the computation of the cross section in the TMD regime, under the approximation that the J/ψ meson is collinear with the outgoing heavy-quark pair. The large-q T limit of the result is then taken and compared with the small-q T limit of the collinear calculation. In Section IV, both results are shown to match after including the smearing of the transverse momentum of the quarkonium in its hadronization, which is encoded in the appropriate shape functions. Conclusions are given in Section V.

II. FROM HIGH TO INTERMEDIATE TRANSVERSE MOMENTUM
In this section the framework of collinear NRQCD is adopted for the description of the process e( ) + p(P ) → e( ) + J/ψ(P ψ ) + X , where all the particles are unpolarized and their four-momenta are given within brackets. This reaction is described by the conventional variables with q ≡ − and q 2 = −Q 2 . We denote by M p and M ψ the masses of the proton and the J/ψ meson, respectively. In the one-photon exchange approximation, at leading order (LO) in NRQCD, the transverse momentum of the J/ψ is due to parton emission in the hard scattering process where parton a can be either a gluon, a quark or an antiquark. The charm-anticharm quark pair is produced in a Fock state n = 2S+1 L [c] J with spin S, orbital angular momentum L, total angular momentum J and color configuration c = 1, 8.
We also introduce the scaling variablesx If we neglect the proton mass and any smearing effects both in the initial and in the final state, we can take and thereforex which impliesx A. Reference frames A convenient reference frame for the calculation of the structure functions for this process is defined by adopting light-cone coordinates with respect to the directions of the relevant hadron four-momenta, P and P ψ . We introduce the light-like vectors n + and n − such that n + · n − = 1. Neglecting the proton mass, Hence the four-momentum of the virtual photon can be written as In the calculation of the cross section, we can perform the following replacement Moreover, from the kinematical constraintŝ we obtain and, from momentum conservation,t By comparing the above relation with the second of Eqs. (12), we get Alternatively, this process can be studied in a frame where the three-momenta P and q are collinear and lie on the z-axis. In this frame, the virtual photon has obviously no transverse momentum. The four-momenta of the particles can be decomposed using two new vectors κ µ + and κ µ − , such that κ 2 + = κ 2 − = 0 and κ − · κ + = 1, where the following relations between the light-like vectors of the two frames hold The partonic Mandelstam variables in this frame read By comparing the above expression fort with the one in Eq. (12), we obtain the relation between the transverse momentum of the photon q T w.r.t. the hadrons, and the transverse momentum of the hadron |P ψ⊥ | w.r.t. the photon and the proton, FIG. 1: Representative leading order diagrams for the partonic process γ * (q) + a(pa) → J/ψ(P ψ ) + a(p a ), with a = g, q,q. The only diagrams contributing to the CS production mechanism are of the type (a), and there are six of them. There are two diagrams for each type (b), (c), (d). The dominant diagrams in the small-qT limit are those of type (c) and (d).
Moreover, using the above expression together with the Sudakov decomposition of P µ ψ in the two frames, namely and using Eq. (19), we obtain

B. Calculation of the cross section
Within the collinear NRQCD approach and in a frame where the longitudinal directions are fixed by the proton and the photon, the cross section for the process under study can be written as follows, where φ ψ is the azimuthal angle of the final J/ψ meson, defined with respect to the lepton plane according to the conventions of Ref. [33], and where H a[n] µ is the amplitude for the hard scattering subprocess γ * a → cc[n] a, with a = g, q,q. The corresponding Feynman diagrams, at the perturbative order αα 2 s , are depicted in Fig. 1. The dominant Fock states included in the calculation are n = 3 S J , with J = 0, 1, 2 [34]. Furthermore, in Eq. (24), f a 1 is the unpolarized PDF, which depends on the light-cone momentum fraction ξ = x B /x of parton a and on the a priori arbitrary hard factorization scale µ, while 0|O(n)|0 is the NRQCD LDME.
The leptonic tensor L µν can be written as [7] L µν = e 2 − g µν Q 2 + 2( µ ν + ν µ ) where the second equality can be obtained from the first one by replacing the expression for µ in Eq. (18), and where the transverse projector g µν ⊥ is given by Furthermore, we have introduced the longitudinal polarization vector of the exchanged virtual photon, which fulfills the relations 2 L (q) = 1 and µ L (q) q µ = 0. From Eq. (21) it follows that the cross section differential in q 2 T can be obtained by simply multiplying Eq. (24) by a factor z 2 , Along the same lines of Ref. [7], the final result can be expressed in terms of four independent structure functions: where the first and second subscripts of the structure functions F denote the polarization of the initial electron and proton, respectively, while the third one, when present, specifies the polarization of the exchanged virtual photon. The expressions of the structure functions in the limit q 2 T Q 2 can be calculated from Eq. (28), replacing the Dirac delta with its expansion in the small-q T limit derived in Appendix A, namely By using the relations we obtain the leading power behavior of the structure functions where a sum over i = q,q is understood. These results are valid up to corrections of the order of O(Λ QCD /q T ) and O(q T /Q). The structure function F cos φ ψ U U is suppressed by a factor q T /Q with respect to the other ones and will not be considered in the following. In Eq. (34), we have defined and where n f refers to the number of active flavors, T R = 1/2, C A = N c , with N c being the number of colors. The symbol ⊗ denotes a convolution in the longitudinal momentum fractions: The well-known LO unpolarized splitting functions read while the splitting functions of an unpolarized parton into a linearly polarized gluon are [35,36] δP with C F = (N 2 c − 1)/2N c . The plus-prescription on the singular parts of the splitting functions is defined, as usual, such that the integral of a sufficiently smooth distribution G is given by and Assuming the validity of the common heavy-quark spin symmetry relations [30] 0|O( 3 P the cross sections for the partonic processes γ * g → cc[n] in Eq. (34) read where e c is the electric charge of the charm quark in units of the proton charge. We note that the relevant partonic subprocesses, contributing to the above cross sections in the small-q T region, are only the n = 1 S [8] 0 , 3 P [8] J ones that correspond tot-channel Feynman diagrams of the type (c) and (d) in Fig. 1.
The structure functions F U U,T and F U U,L in Eq. (34) exhibit logarithmic (collinear) singularities as q T → 0. Their behavior is similar to the analogous structure functions for light hadron production in SIDIS discussed in Ref. [4], where the dominant underlying partonic process is γ * q → q. There are some differences though. In the latter case, the logarithmic term L is given by C F (2 ln Q 2 /q 2 T − 3) instead of Eq. (36). The color factor C F clearly corresponds to a quark initiated process, while C A corresponds to a gluon initiated one. The two different finite terms originate from the virtual corrections to the splitting functions P qq and P gg , respectively. Moreover, in light hadron production extra terms appear, containing convolutions of FFs with the P qq and P gq splitting functions, which cannot be present in our calculation for quarkonium production within the NRQCD framework. We also point out that the structure function F cos 2φ ψ U U does not contain any large logarithm in the region q T Q, whereas the corresponding observable for light hadron production diverges logarithmically and is suppressed by an overall factor q 2 T /Q 2 . Finally, the appearance of a logarithm ln(Q 2 + M 2 ψ )/q 2 T , instead of ln Q 2 /q 2 T , suggests Q 2 + M 2 ψ as the natural choice for the hard scale in our process.

III. FROM SMALL TO INTERMEDIATE TRANSVERSE MOMENTUM
The process γ * g → cc[n] has been calculated in Ref. [13] within the TMD framework, taking into account the intrinsic transverse momentum effects of the gluons inside the proton. If we neglect smearing effects in the final state, i.e. if we assume that the final J/ψ meson is collinear to the cc pair originally produced in the hard scattering process, the cross section can be cast in the following form where f g 1 and h ⊥ g 1 are, respectively, the unpolarized and linearly polarized gluon TMDs inside an unpolarized proton [37][38][39][40][41].
We note that, beyond the parton model approximation, for those processes where TMD factorization is valid, soft gluon radiation to all orders is included into an exponential Sudakov factor, which can be split and its parts absorbed into the TMD PDFs and FFs involved in the reaction, whereas the remaining perturbative corrections are collected into a hard factor H. As a consequence of the regularization of their ultraviolet and rapidity divergences, TMDs depend on two different scales, not explicitly shown in the above equations. In the following we will take these two scales to be equal to each other and denote them by µ, to be identified with a typical hard scale of the process.
TMDs can be calculated perturbatively in the limit p T ≡ |p T | Λ QCD . This can be better achieved in the impact parameter space. To this aim, we focus first on f g 1 (x, p 2 T ) and its Fourier transform, which is defined as with J n being the Bessel function of the first kind of order n. The perturbative part of the gluon TMD, valid in the limit b T 1/Λ QCD , reads where f a 1 (x, µ) are the collinear parton distribution functions for a specific (anti)quark flavor or a gluon a, and The coefficient functions C g/a and the perturbative Sudakov exponent S A , which resums large logarithms of the type ln(b T µ), are calculable in perturbative QCD. The coefficient functions can be expanded in powers of α s as follows while the perturbative Sudakov factor at LO reads Therefore, in the small-b T limit, namely b T 1/Λ QCD , and at LO in α s we find Using the DGLAP equations we can evolve f g 1 from a scale µ down to another scale µ b < µ and obtain where a sum over i = q,q is understood. By substituting the above expression into Eq. (50) we get Transforming back to momentum space, we find the transverse momentum distribution in the region p T Λ QCD , where we have used the following integrals By substituting the expression for f g 1 (x, p 2 T ) given in Eq. (53), evolved to the scale µ 2 = Q 2 + M 2 ψ , into Eq. (45), we find that the TMD structure functions F U U,T and F U U,L do not exactly match the corresponding collinear ones in the small-q T limit given in Eq. (34). In the intermediate region Λ QCD q T Q, we get This suggests that one needs to include smearing effects in the final state as well, through the inclusion of a suitable shape function [31,32], to be convoluted with the TMD in momentum space. Imposing the validity of the matching will give us the LO expression of the shape function, as we shall see in the next section. We now turn to the polarized gluon distribution h ⊥ g 1 (x, p 2 T ) and the structure function F cos 2φ ψ U U . From Ref. [35], we know the perturbative tail of h ⊥ g 1 (x, p 2 T ) at LO in terms of the unpolarized collinear PDFs f a 1 (x, µ 2 ), The above expression, together with Eq. (45), leads to . This is achieved without the need of any shape function because of the absence of a logarithmic term at the perturbative order we are considering.

IV. TMD FACTORIZATION AND MATCHING WITH THE COLLINEAR FRAMEWORK
On the basis of the above considerations, TMD factorized expressions for the structure functions F U U,T and F U U,L have to take into account smearing effects [13], encoded in the shape function ∆ [n] [31,32], which can be thought as a generalization of the long distance matrix elements of NRQCD in collinear factorization. We start by assuming the validity of the following formulae, where the H represent the hard parts, which can be calculated pertubatively. Moreover, we have introduced the transverse momentum convolution As expected, in absence of smearing, ∆ [n] (k 2 T ; µ 2 ) = 0|O(n) |0 δ 2 (k T ), and the convolution in Eq. (59) reduces to the product of the LDME 0|O(n) |0 with the gluon TMD f g 1 (x, q 2 T ). Furthermore, this convolution can be expressed as follows where we have introduced the Fourier transform of the shape function, We are now able to show that the following expression of ∆ [n] , valid at LO in α s , which leads to in the limit |k T | Λ QCD , will solve the matching issue of the TMD and collinear results in the region Λ QCD q T Q. In fact, by plugging it together with Eq. (52) into Eq. (60), with the choice µ 2 = Q 2 + M 2 ψ , we get Substituting the last formula in the LO expressions for F U U,T and F U U,L in Eqs. (58), we recover the correct results for F U U,T and F U U,L in Eqs. (34). Along the same lines of Ref. [11], in which a TMD factorized formula has been proposed at the twist-three level for the cos φ asymmetry for light-hadron production in SIDIS, we can conjecture that the formulae in Eq. (58) are valid to all orders in α s , provided one includes also the nonperturbative contributions of the TMD gluon distribution and shape function.
At this point, it may be good to stress the difference between the shape function and the wave function of the J/ψ meson. In the lowest-order picture a quarkonium state consists of a heavy quark and an antiquark. In the center-of-mass frame of the quarkonium the momenta k 1 and k 2 of the heavy quark and antiquark add up to zero. Their difference defines the relative velocity v = |v|: k ≡ k 1 − k 2 = m Q v, and NRQCD involves an expansion in v 1. The wave function Ψ(k) of the quarkonium in momentum space is expected to be positronium-like, with a tail that depends on L: the orbital angular momentum of the quark-antiquark pair. What enters in the expressions for L = 0 states is dk Ψ(k) = R(0), the wave function at the origin, and for L = 1 states its derivative R (0). On the other hand, the shape function ∆ [n] of the QQ system is a function of k 1 + k 2 and, in LO NRQCD, is equal to a delta function in k 1 + k 2 . Upon radiating additional gluons and quarks, this becomes smeared out. That is the reason why in Ref. [13] the shape function was referred to as smearing function. For lack of better input, the model adopted there for this function was based on the expectations for the wave function Ψ(k). However, what the results presented in this paper show is that the shape function (or at least its perturbative tail) is actually independent of L. An L-independent shape or smearing function would imply that it would not be an obstacle to the extraction of the CO matrix elements from a comparison between quarkonium production and open heavy quark pair production in SIDIS as proposed in Ref. [13], implying a much more robust result. Note that the transverse momentum dependence of the shape function in Eq. (63) not only implies independence from the L quantum number, but actually from any quantum number of the produced quarkonium. Only the overall magnitude of the shape function is a function of these quantum numbers determined by the relevant LDME.
For the angular dependent structure function we expect the following result to hold, even if due to the absence of a collinear divergence a shape function is not strictly needed: where with w(p T , k T ) being a transverse momentum dependent weight function. The shape function ∆ [n] h could be in general different from ∆ [n] : the determination of its perturbative tail would require a similar study at higher order in α s . We note however that a full calculation of the cross section for J/ψ production in SIDIS at the order α 2 α 3 s , within the collinear NRQCD framework, is still missing [42].
Based on the fact that the p T dependence for h ⊥ g 1 in the gluon correlator has a rank-two tensor structure in the noncontracted transverse momentum, and unpolarized vector-meson production generally has a rank-zero structure, we consider a shape function ∆ [n] h of rank zero (ignoring a possible contribution from a linearly polarized quark-pair state) and a weight function expression: Furthermore, the convolution in Eq. (66) can be rewritten as where we have introduced the second derivative w.r.t. b 2 T of the Fourier transform of the linearly polarized gluon TMD distribution, Concerning the L-independence of ∆

[n]
h at high transverse momentum, we cannot draw any conclusion, due to the absence of a logarithmic dependence. Moreover, we cannot conclude any L-independence for small transverse momentum for either ∆ [n] or ∆ [n] h , but the suggestions of Ref. [13] and the proposed cross-checks, do allow an experimental investigation of this. Shape functions need to be experimentally extracted, just like the LDMEs and FFs have to. The EIC can play an important role in this regard.

V. CONCLUSIONS
Let us recapitulate the main points of this work. Our starting point is the assumption that transverse momentum dependent factorization is valid for J/ψ production in SIDIS at small q T . This Ansatz is a very reasonable one, since SIDIS for light hadrons is one of the few processes for which TMD factorization is proven at all orders in α s . Building further on this premise, we calculate the cross section for e p → e J/ψ X in two different regimes. At low q T , the cross section is factorized in terms of TMD PDFs and generic shape functions (first introduced in Refs. [31,32] as the generalization of the NRQCD LDMEs to TMD factorization), while at high q T the factorization involves collinear PDFs and NRQCD matrix elements. Our Ansatz of TMD factorization then requires the consistency condition that both descriptions match in the intermediate region Λ 2 QCD q 2 T Q 2 . The perturbative calculations for the leadingpower SIDIS structure functions F U U,T and F U U,L allow, at least at LO accuracy in the strong coupling constant, to deduce from the matching the specific form of the color-octet shape function at large transverse momentum. Moreover, the angular structure function F cos 2φ ψ U U matches without any dependence on a shape function whatsoever, due to the absence of a logarithmic divergence.
We therefore conclude that the assumption of TMD factorization for J/ψ production in SIDIS, and the necessary (but not sufficient) condition of matching at intermediate q T , imposes certain properties to the perturbative structure of the shape functions. In particular, we find that their perturbative tails are independent of the quantum numbers of the produced quarkonium state, except for their overall magnitude given by the collinear NRQCD LDMEs. One consequence of this is that the feasibility of extracting the CO LDMEs by comparing quarkonium and open heavyquark production at an EIC, proposed in Ref. [13], is not hampered.
We note that our perturbative result agrees with the findings in recent work based on the soft-collinear effective theory (SCET) approach [32], although we are unable to draw conclusions on the nonperturbative structure of the shape functions with our method. Corroborating the results in Ref. [32] from a different and arguably simpler approach, we thus believe that our study could contribute towards a full proof of TMD factorization for quarkonium production in SIDIS.

VI. ACKNOWLEDGEMENTS
We thank Werner Vogelsang for useful discussions on the derivation of Eq. (A25). The work of U.D., F.M. and C.P. is supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No.