Transverse momentum dependent shape function for J/ψ production in SIDIS

,


I. INTRODUCTION
In recent years, heavy quarkonium production in various inclusive processes has attracted great interest as a way to probe the transverse momentum dependent (TMD) gluon distributions [1][2][3][4][5][6][7][8][9][10][11].In this paper, we focus on J/ψ production in semi-inclusive deep inelastic scattering (SIDIS), e (ℓ) + p (P ) → e ′ (ℓ ′ ) + γ * (q) + p (P ) → e ′ (ℓ ′ ) + J/ψ (P ψ ) + X , where the particle momenta are given between brackets and the virtual photon momentum is given by q = ℓ − ℓ ′ .The J/ψ mass M 2 ψ = P 2 ψ and the photon virtuality Q 2 = −q 2 > 0 are considered hard scales in the process, i.e. they are considered much larger than the nonperturbative QCD scale Λ QCD , although most results will also be valid for photoproduction (Q 2 = 0).The electron and proton masses will be neglected w.r.t.M ψ and Q whenever possible.The virtual photon transverse momentum is denoted by q T and can be directly related to the J/ψ transverse momentum P ψ⊥ .The distinct subscripts used for the transverse momentum components, specifically "T " and "⊥", serve to emphasize the different frames in which they are measured.In particular, we consider q T when both the target proton and the J/ψ have no transverse components and P ψ⊥ when the photon and the proton have only longitudinal components.
Depending on the value of |q T |, we can identify two different transverse momentum regions, see Fig. 1.The high transverse momentum (HTM) region is given by the condition |q T | ≫ Λ QCD , while the low transverse momentum (LTM) region corresponds to |q T | ≪ µ H .Here µ H = f (Q, M ψ ) with f (Q, M ψ ) ≳ M ψ generically denotes the hard scale of the process.The cross section can be evaluated within the two transverse momentum regions by adopting the proper factorization that enables to separate the short-distance from the long-distance contributions.The collinear factorization is applicable at HTM, while the TMD factorization is expected to be valid at LTM [12], for which the cross section is sensitive to TMD quantities.In addition, we can identify an intermediate transverse momentum (ITM) region, namely Λ QCD ≪ |q T | ≪ µ H , where both factorizations are valid.Since our attention will be mostly directed towards this overlapping region, where |q T | (or equivalently |P ψ⊥ |) becomes small compared to the hard scale, we will neglect any transverse momentum dependence in f (Q, M ψ ).To describe J/ψ hadronization we employ nonrelativistic QCD (NRQCD) [13], in which the heavy-quark pair forms a Fock state, specified by n = 2S+1 L [c] J : S denotes the spin, L the orbital angular momentum, J the total angular momentum and c the color state of the pair.Note that the pair can couple either as a color-singlet (CS), with c = 1, or as a color-octet (CO) state, with c = 8.The (low-energy) transition from this general state to the J/ψ is encoded in the nonperturbative Long-Distance Matrix Elements (LDMEs) that are distinct for each quarkonium Fock state.States with different quantum numbers n do not interfere as the cross section is proportional to a direct sum of LDMEs, up to a required precision in the expansion w.r.t.v, which corresponds to the (non-relativistic) relative velocity of the heavy quark-antiquark pair in the quarkonium rest frame.In the following we will truncate the expansion up to the relative order v 4 , including the 3 S [1] 1 CS state and the 1 S J CO states.Note that, in the following we will not consider the interference among P -wave states since it is not necessary in the evaluation of the unpolarized differential cross section.However, we have taken them into account in our brief digression on the production of polarized J/ψ mesons in SIDIS (see Sec. IV).
In Refs.[14,15] it was found that the TMD factorized expressions have to take into account final state smearing effects that are encoded in the TMD shape function (TMDShF).This nonperturbative hadronic quantity describes the transition from the heavy quark pair to a bound quarkonium state, which not only contains the formation of the bound state in terms of an LDME, but also the transverse momentum effects that arise from the soft-gluon radiation.
In Refs.[16,17] the matching procedure in SIDIS has been investigated, according to which the TMD and collinear expressions are compared in the ITM region.It was found that the introduction of TMDShFs solves the mismatch between the collinear and TMD expressions, by resumming |q T | divergences in the Sudakov factor.However, this term is in contradiction with other studies, as it has been demonstrated that no double logarithms in the nonperturbative Sudakov factor associated to heavy quark production are present for pp → (J/ψ or Υ) + X [18] and for open heavyquark pair production, both in ep [19] and pp [20] collisions.The absence of the double logarithms in J/ψ production can also be seen in Ref. [21].Due to this discrepancy, universality was assumed in [22] using the explicit result of [18]; however, as we will see, this only holds for photoproduction, not electroproduction.
We found that the discrepancy in the ep matching study arises from the presence of discontinuities in the structure functions that appear in the small-q T limit of collinear factorized expressions.These structure functions contain a Dirac delta function for which a small-q T approximation is applied.The approximation employed in Refs.[16,17], which is an extension of a well-known expression [23] to the heavy quarkonium case, would be valid when multiplied by a continuous function, but that turns out to be invalid in the present case of discontinuous hard scattering factors.In this article we show how to properly treat these expressions to resolve this discrepancy.
In addition, we extend our analysis to single quarkonium production in pp collisions, where the hard scale is a function of the quarkonium mass only: µ H = f (M ψ ) with f (M ψ ) ∼ M ψ .This allows us to test the connection of TMDShFs obtained in different cases, i.e. to study their universal properties.Even if we expect that the LDMEs are process independent in the collinear description, the same is not necessarily true for the TMDShFs.Indeed in the latter process dependences may arise due to the transverse momentum exchange with other colored objects.
The paper is organized as follows.In Sec.II we revise the matching procedure.In particular, in Sec.II A we discuss the pole structure of the collinear cross section in the small transverse momentum limit in detail, while the new TMDShF results for SIDIS are presented in Sec.II B. In Sec.III we address the aforementioned process dependence, comparing the TMDShFs in SIDIS and in pp collisions.Conclusions are given in Sec.IV, together with a summary of our findings.In addition, there are two appendices at the end of this paper.In Appendix A we present a more complete derivation of our method to include the pole structure contributions in our results.In Appendix B we derive the soft gluon emission from the Born amplitude obtained through the eikonal approximation.

II. THE MATCHING PROCEDURE
The SIDIS reaction in Eq. ( 1) is described by the conventional kinematical SIDIS variables We consider a frame where the virtual photon has no transverse momentum component, and we identify two light-cone directions n + and n − , for which n + • n − = 1.With these, the Sudakov decomposition of the relevant momenta can be written as where P 2 ψ⊥ = −P 2 ψ⊥ is the squared J/ψ transverse momentum (w.r.t. the photon and proton), while M ψ⊥ = M 2 ψ + P 2 ψ⊥ is the J/ψ transverse mass.In particular, we will consider the fully unpolarized differential cross section dσ/(dx B dy dz dq 2 T dϕ ψ ), where ϕ ψ is the J/ψ azimuthal angle measured w.r.t. the lepton plane.Moreover, we replaced the transverse momentum of the J/ψ with that of the photon q T (evaluated w.r.t. the hadrons); this replacement is achieved via The differential cross section can be parameterised in the HTM region as follows [16] dσ dx B dy dz dq 2 T dϕ ψ where the first two subscripts of the structure functions F refer to the polarization of the initial (unpolarized) proton and electron.The last subscript in F U U,P with P =⊥, refers to the virtual photon polarization (transverse or longitudinal), while for F Φ U U with Φ = cos ϕ ψ , cos 2ϕ ψ the superscript refers to the angular term that accompanies it.Henceforth, we will refer to the aforementioned hard scattering structure functions via the general notation F Φ U U,P .On the other hand, the same differential cross section evaluated in the LTM region is given by where the structure function , being subleading power/twist, has not been included.Note the difference in the structure functions: F Φ U U,P are evaluated in collinear factorization, while the calligraphic F Φ U U,P are calculated within transverse momentum factorization.

A. From high to intermediate transverse momentum
In this section we provide a systematic method to investigate the small-q T limit of SIDIS observables at HTM (|q T | ≫ Λ QCD ).Adopting the parton model, the production of a J/ψ possessing a high transverse momentum component is possible at the lowest order in α s via where a can be either a quark, antiquark, or gluon.In this kinematical regime we adopt collinear factorization, for which The perturbative amplitude squared |M| 2 for the hadronic process in Eq. ( 1) is obtained by contracting the leptonic tensor L µν with the amplitude H (a) [n] µ , describing the partonic process in Eq. (7), and its conjugate H In particular, the lepton tensor can be written as follows: where we introduced the tensors Moreover, g µν ⊥ is the transverse projector while ϵ µ L (q) is the longitudinal polarization vector and lµ ⊥ is the unit vector along the transverse component of ℓ, w.r.t. the photon-proton axis.Henceforth, we refer to one of the tensors in Eq. ( 10) via the general notation ϵ µν P;Φ , where P =⊥, and Φ = cos ϕ ψ , cos 2ϕ ψ .Employing this, the structure functions introduced in Eq. ( 5) can be evaluated via where the sum n runs over the dominant LDMEs ⟨O[n]⟩ and a runs over the parton types.Furthermore, we introduced the partonic scaling variables together with In Ref. [16] the Dirac delta present in Eq. ( 13) was expanded at small-q T as follows (see its Appendix B for the derivation) where Note that on the right-hand side of Eq. ( 16) the coefficient in front of the double delta logarithmically diverges with q T .However, as we have found this is not sufficient to obtain the correct behavior of structure functions in the ITM region, which is restored by adding a constant term to the double-delta coefficient. 1 The need to include this subdominant term can also be understood from the following argument.The Dirac-delta expansion in Eq. ( 16) was obtained in Ref. [16] by applying the full Dirac delta to two continuous test functions.However, the structure functions defined in Eq. ( 13) contain discontinuities that come from the soft gluon radiation associated with the CO final state in the NRQCD calculations (see Appendix B).These contributions are made explicit via the decomposition into poles through a Laurent expansion, namely2 1 (4π) where H  18) is in agreement with Refs.[24,25].Note that to get the pole structure on the right-hand side of Eq. ( 18) we are explicitly writing the amplitude squared in terms of x′ and ẑ (Eqs.( 14) and ( 17)).We have found that for J/ψ production in SIDIS the poles are present only for the gluon-initiated process γ * g with the expansion running up to k = 2. Instead, the quark-initiated processes γ * q are fully described by the k = 0 finite term.Moreover, up to the precision considered in this work, the poles contribute only to the structure functions F U U,P introduced in Eq. ( 5).
These poles are under control when the amplitude squared is evaluated at high-q T values, as the transverse momentum forces the phase space to deviate from ẑ = 1 and x′ = 1.Solely when we consider the small-q T limit they have a significant impact.The Dirac-delta expansion in Eq. ( 16) is applicable only to the first term (H ), while all the others require a different approach.In particular, we can split the differential cross section in Eq. ( 5) into three parts in the HTM region, namely with where the function G(x ′ , ẑ) is given by The difference in the lower integration limit for both x′ and ẑ between Eqs. ( 13) and ( 20) has been introduced for the convenience of the calculation.This modification is possible since the added integration range does not contribute to the final result (see Appendix A of this work and Appendix B of [16]).As mentioned, one can apply directly the delta expansion in Eq. ( 16) to evaluate the small-q T behavior of dσ A .Instead, the expansions of dσ B and dσ C are obtained by considering the integral w.r.t.dx and dẑ of those terms that are truly indeterminate in the limit q T → 0, with the indeterminacy solved by the presence of the full Dirac delta (see Eq. ( 21)).Therefore, it is legitimate to approximate H [n]; (k) P (ẑ) → H [n]; (k) P (1) which gives Hence from Eq. ( 20) we obtain and where x ≡ x B /x max .More details on the previous results can be found in Appendix A. Since the small-q T limit of these quantities is proportional to H [n]; (0) (1, 1), we can effectively add these terms to the double delta coefficient, obtaining that with Considering the contributions from the various terms in Eq. ( 26), we found that the small-q T limit is dominated by the first two terms for this γ * g channel (and similarly the first two terms of Eq. ( 16) for the quark and antiquark channels).Instead, the contribution coming from the "+"-distribution of ẑ is subdominant and can be neglected in the following.In principle, this last term will lead to a fragmentation-like contribution to the process considered here.We will further comment on the connection to a fragmentation description below.Hence, the leading power behavior of the structure functions in the ITM region is given by while is suppressed by a factor of |q T |/µ H w.r.t. the other structure functions.This is in accordance with the TMD formula in Eq. ( 6) which does not show any cos ϕ ψ contribution too.The logarithmic function L(q 2 T ) reads and the quantities σ U U,P and σ are related to the partonic process γ * g → cc[n] and they correspond to (see Ref. [16]) Moreover, P ab in Eq. ( 27) denotes the leading order, fully unpolarized splitting functions, that can be found in Ref. [26], while δP ab are the splitting function of an unpolarized parton into a linearly polarized gluon, which can be found in Refs.[27,28].The convolution (denoted by the "⊗" symbol) between these splitting functions and the parton distribution functions is defined as where P ab denotes either P ab or δP ab .
The logarithmic function defined in Eq. ( 28) is our most important difference compared to Ref. [16], where the logarithmic function contains twice the logarithm log[(M 2 ψ + Q 2 )/q 2 T ] compared to Eq. ( 28).This is due to the presence of the poles, not considered in Ref. [16].Indeed, it is through the inclusion of Eqs. ( 23) and ( 24) that in Eq. ( 28) one of the logarithms has been removed.The price to pay corresponds to the novel q T -independent terms found, namely Clearly, Eq. ( 28) has an impact on the TMDShF derivation too, as will be discussed in Sec.II B. Besides, Eq. ( 28) implies the presence of divergences related to soft gluon emission from the leading order γ * g → cc[n] process.It is then possible to check the validity of this expression by investigating the soft-limit of Eq. ( 7) via the eikonal method, as done in Appendix B.
Although our work is based on the J/ψ production in SIDIS, we expect that the presence of the poles as in Eq. ( 18) is an intrinsic feature of any inclusive quarkonium production, and they apply to different processes and observables too.Hence, the suppression of the q T -logarithm in Eq. ( 28) is not an exclusive outcome of the specific process under consideration, but rather a general statement.Thus, these discontinuities may be connected to other regularization procedures associated with CO contributions to heavy quarkonium productions.While it is worthwhile to further pursue these connections in the NRQCD factorization, we consider such a study to be beyond the scope of the current paper but we hope to address it in the future.However, to emphasize the importance of further investigation, we will briefly comment on the similarities of our findings with those obtained by adopting the fragmentation function description.
The same cross section in the HTM region can be expressed in terms of fragmentation functions, as shown in Refs.[29][30][31].Hence, the TMDShF may also be seen as a fragmentation-like function of a cc into a J/ψ evaluated at ITM.The evolution of the latter has been studied in Ref. [30], which includes real contributions having a component proportional to the (1 − ẑ) + distribution and another one to δ(1 − ẑ).Hence, our analysis is related to the latter term.However, an important difference concerns the integration range of the outgoing gluon.The integration of the softgluon momentum in our case has a lower limit set by the J/ψ transverse momentum (see Eq. (B11) of Appendix B), whereas no lower limit is present in Ref. [30] causing infrared divergences.Therefore, the connection between our work and the fragmentation-function description cannot be carried out further without the inclusion of next order (real and virtual) contributions and, as previously stated, we leave this discussion to further studies.

B. From low to intermediate transverse momentum
In this section we evaluate the evolution of the structure functions valid in the LTM region (|q T | ≪ µ H ) up to the ITM region.Even if not formally proven, there are strong arguments in favor of TMD factorization [12].Therefore, the differential cross section for the semi-inclusive production of a J/ψ with a small transverse momentum component is given by Eq. ( 6).In this case, the structure functions F can be calculated from the partonic process where, contrarily to the HTM case, the initial gluon has a non-negligible transverse momentum component w.r.t. the parent proton, namely with p 2 T = −p 2 T .Hence, Eq. ( 31) leads to Following Refs.[14][15][16][17]22], we consider from the beginning the TMD factorized formula that includes the presence of a TMDShF, ∆ [n] or ∆

[n]
h , which is related to the production of a J/ψ with a small transverse momentum component w.r.t. the photon and proton.As commented below Eq. ( 26), in this paper we focus on the δ(1 − ẑ) contribution from the TMDShF.Subdominant terms and higher order corrections are also expected to contribute away from z = 1.In that case the description of the heavy quark pair that hadronizes into the heavy quarkonium state will be even more similar to a single-parton TMD fragmentation functions description applied in light hadron production.See also Ref. [32] for a description involving both single parton and parton-pair fragmentation processes.
In the above equations convolutions between the TMD distributions and the TMDShFs appear, namely where in the last line we introduced the weight function w defined as (with M p the proton mass) Beyond the parton model approximation, soft gluon radiation to all orders is included into an exponential Sudakov factor.One can relate its logarithmic divergences to the TMD objects (both PDFs and shape function) involved in the reactions, whereas the remaining perturbative q T -independent corrections are collected into the hard term.As a consequence of the regularization of their ultraviolet and rapidity divergences, TMD-PDFs depend on two different scales, respectively µ and √ ζ.We take these two scales to be equal and denote them by µ.In contrast, there are no rapidity divergences associated to the TMDShF.Thus, we can impose for its rapidity parameter ζ ∆ = 1, in line with Ref. [33].
Implementing TMD evolution is more easily done in impact parameter space, where convolutions in the cross section become simple products.Besides, in Ref. [16] it was found that, up to the precision considered, from the matching procedure it is possible to deduce only the naive order-α 0 s part of Note that in reality smearing effects will be involved, but that small-k T behavior cannot be obtained from a perturbative matching calculation, at least up to the perturbative order considered here.As a consequence, in the following we focus on the convolution C[f g 1 ∆ [n] ].We define the Fourier transform of f a 1 (x, p 2 T ) as and the Fourier transformed TMDShF as from which where we fixed the factorization scale so that the convolutions are evaluated at the hard scale.The perturbative tail of the fully unpolarized gluon TMD f g 1 , valid in the limit |b T | ≪ 1/Λ QCD , is given by [26] where 123.Note that the coefficient function C a/b in Eq. ( 39) can be expanded in powers of α s and can be explicitly found in Ref. [26,34].Nevertheless, the coefficient g/a in the right-hand side of Eq. ( 40) will not enter in the following (leading order) discussion since they are independent of the parameter b T .Consequently, their explicit expression at all orders is not required.Furthermore, the (leading order) Sudakov factor S A present in Eq. (39) reads where in the last line the running of the coupling has been neglected.By inserting Eqs. ( 40) and (41) in Eq. ( 39) and using the DGLAP equations to evolve the PDF from a scale µ H down to the scale µ b < µ H , we find that up to order α s the perturbative tail of the gluon TMD-PDF reads [34] where once again P ab denotes the leading order splitting functions [26].Employing this, and by requiring that the TMD expressions evolved to the scale µ 2 with the expansion of the collinear ones obtained in Eq. ( 27), we deduce the TMDShF perturbative tail which in momentum space becomes valid in the |k T | ≫ Λ QCD limit.Inserting Eq. (43) in Eq. ( 38), we find that the convolution in momentum space is given by where L(q 2 T ) is the logarithmic function defined in Eq. ( 28).Hence, with the choice µ = Q, the first two lines of Eq. ( 33) and the first line of Eq. ( 27) match.Note how the modification of Eq. ( 28) compared to Ref. [16] has a significant impact on the TMDShF expression.Indeed, the TMDShF perturbative tail in Eq. ( 44) does not contain any kind of logarithmic divergence in k T , being tamed by the presence of the heavy mass.We emphasized that the absence of k T -divergent terms associated to the quarkonium is in accordance with other works in the literature, e.g.Refs.[14,15,[18][19][20][21].
For completeness, we remark that the matching of , which involves the second convolution in Eq. (34), is fulfilled by taking the perturbative tail of h ⊥g 1 [27] up to α s order and the leading order naive shape function from which Since the h ⊥g 1 expansion starts at order α s , we notice that to get the non-trivial perturbative tail of ∆ h it is required that the SIDIS cross section within NRQCD is evaluated at order αα3 s .However, this calculation is currently unavailable.

III. UNIVERSALITY
In the previous section, we found that an extra factor ∆ is needed to absorb all the q T -divergent terms coming from the collinear limit, and we identified it as the dominant TMDShF perturbative tail.However, it has been obtained at the particular scale Q, whereas for more general application it needs to be considered at a general scale µ H .This can be obtained by tracing back the µ H dependence in Eq. (28), that is related to the full Sudakov factor for J/ψ production in SIDIS in terms of this general scale and up to order α s , namely where We checked that Eq. (49) (and subsequently Eq. ( 50)) agrees in the kinematic limit corresponding to a bound pair with the Sudakov factor obtained in the open heavy-quark pair production in electron-proton collisions, which can be found in Ref. [19].
It is not natural to fully include Eq. ( 50) into something that we identify as the TMDShF.Indeed, being a quarkonium-related object, its complete dependence is given by ∆ , while it may depend on the process-related hard-quantity Q only via the µ H choice.Thus, the Q 2 dependence deriving from Eq. (50) must stem from a process dependent part, which can be incorporated into an extra process-dependent factor ep into these two terms: 3 The ∆

[n]
ShF is what we truly identify as the TMDShF and is universal because it solely depends on M ψ .Instead, the S ep is an extra soft factor which incorporates the specific process dependence and it can be removed by a proper choice of the factorization scale µ = µ H .This implies that at that scale the full ∆ ep is equivalent to the TMDShF.At this level, the simplest way to perform the splitting in b T -space is to take With this splitting convention and by taking µ H ≡ Q, the full ∆ [n] ep reduces to the TMDShF, implying that the latter is given by Eq. (44).
To test the proposed factorization, one may consider another process and check if it is possible to identify the same TMDShF in Eq. ( 52).We take into account J/ψ production in hadron collisions, namely pp → J/ψ + X. 4 For this process the small-q T behavior of the cross section evaluated in the HTM region has been calculated in Ref. [18].The corresponding Sudakov factor can be written as where in which the first term of Eq. ( 55) is directly related to the δ 8c term in Ref. [18].Also in this case we checked that previous equations agree in the kinematic limit corresponding to a bound pair with the open heavy-quark pair production Sudakov factor, which can be found in the literature (for instance Ref. [20]).Moreover, even if it is possible to produce quarkonia in a CS state (e.g.η c ), for pp our perturbative tail only applies to CO states.Despite this, we cannot exclude that a non trivial TMDShF perturbative tail applies to the CS channel too, if one goes to next orders in perturbation theory.
Although the full ∆ [n] ep , we can still identify the same ∆ [n] ShF in Eq. ( 52), which is now combined with a different (extra) soft factor S pp , namely with Interestingly, for S pp the coefficient in front of the log is "3", whereas the same coefficient for S ep is "2", which corresponds to the number of TMD quantities (PDFs and shape functions) involved.Hence, even if process dependent, these terms are the same apart from the number of TMDs involved.This may allow to guess the required term for other processes, such as for di-quarkonium production in pp collisions (if that factorizes at all for CO-CO production).
The factor S pp reduces to ShF (M 2 ψ ).For this scale choice, ∆ [n] ShF (M 2 ψ ) is compatible with the corresponding one presented in Ref. [15] for χ c decay into light-quarks, where the NLO TMDShF up to corrections of O(|k T | −1 ) is given by a constant too.
According to our findings, in principle one may obtain the value of . Hence, we propose a strategy for the extraction of the TMDShF from different processes, relying on their factorizability.For processes where we have a dominant hard scale it is reasonable to expect that by setting µ H equivalent to it we reduce our uncertainties in the extraction of the TMDShF. 5Then, this term can be re-used for every process involving J/ψ by evolving ∆ For completeness, we mention that the soft factor derived for the open heavy-quark pair production also involves an additional process-dependent factor [35][36][37] (which is sometimes denoted by ∆, but should not to be confused with ours).This additional factor stems from soft radiation in the Q Q production and can in principle even depend on the angle of q T .Hence, it is natural to expect an additional process-dependent soft term in the quarkonium case too.In that sense we expect that our extra soft term S will acquire azimuthal and rapidity dependences if one goes beyond the order and approximation we have considered, as they are present in the ∆ quantity of Refs.[35,37].

IV. CONCLUSIONS
In this work, we revised the procedure to derive the leading order TMDShF perturbative tail for heavy quarkonium production.We focused on the SIDIS unpolarized cross section, which is parameterized in terms of structure functions.In particular, we considered the cross section evaluated at low q T and order αα s , which involves the convolution between the gluon TMD-PDF and a general TMDShF, taking the reasonable assumption that factorization holds.This description should match the collinear one at high q T and order αα 2 s when both are evaluated at intermediate q T , namely Λ QCD ≪ |q T | ≪ µ H .We emphasize that, although the exact choice of µ H is important from a phenomenological point of view where it is advantageous to extend the intermediate-q T region, our findings hold for any choice of µ H .
We show that in the high transverse momentum region, these structure functions present poles when the small-q T limit is taken.We expect that these poles will be contained in other hard amplitudes concerning inclusive quarkonia production.Therefore, we presented a systematic way to deal with them, showing how they provide non-negligible terms in the expansion at small q T .These terms, neglected in [16,17], significantly alter our findings of the TMDShF perturbative tail.At variance with previous works, it does not present a logarithmic dependence on the transverse momentum (double-logarithm in b T -space), which makes them different from usual TMD fragmentation functions for light hadron production.However, this non-logarithmic dependence is in agreement with other works [14,15], and with the Sudakov factors obtained for open heavy-quark pair production in electron-proton and proton-proton collisions.
We remark that our results on the transverse momentum dependence of the TMDShFs hold for every CO quarkonium state with the same quantum numbers as the J/ψ we considered, e.g.Υ(nS) and ψ(2S).The magnitude of TMDShFs can be different though and is determined by the LDMEs.This conclusion holds up to the precision considered, corresponding to the αα 2 s and v 4 orders in the NRQCD double expansion.Moreover, the same considerations apply if we take into account the polarization of the J/ψ, since the kinematics is the same.Namely, we have the same TMDShF perturbative tail for both the longitudinal and transverse J/ψ polarization states.Besides, to check that the same form of the TMDShF applies for observables involving h ⊥g 1 we would require the computation of the cross section within NRQCD at higher order in α s , both for polarized and unpolarized J/ψ productions.However, this calculation is still unavailable.
Furthermore, we showed that if we consider the evolution w.r.t. the factorization scale µ, the TMDShFs would have to depend on the hard scale Q too.As it is not reasonable to include this dependence into a quantity that is related to the quarkonium formation solely, we considered a split into two terms: a process-independent quantity that we identify as the universal TMDShF, and an extra process-dependent soft factor.This then allows to make a connection between ep and pp processes, without losing predictability completely.It is also in line with results for open heavy quark production, where extra process dependent soft factors are also required, at least in pp collisions [35,37].
Despite the process dependence, we showed that it is possible to extract the universal TMDShFs by appropriate choices of scales, which allows to relate different processes.Hence, we expect that with the upcoming Electron-Ion Collider and more data provided by pp facilities (e.g.LHC in fixed target mode) extractions of the TMDShFs will become available in the future and new features of heavy quarkonium production will be uncovered.
where we recall that x = x B /x max .Note how the last lines in Eqs.(A3) and (A4) are respectively equivalent to what is presented in Eqs. ( 23) and (24).
if the bound state is produced in a CO configuration (the relation is still independent of the other quantum numbers).By combining Eqs.(B4) and (B5) we obtain the full amplitude that includes the soft gluon radiation from both the incoming gluon and the outgoing (CO) cc pair, namely Averaging over colors and using Eq.(B1), we then find that where Considering a frame where q and p a are along the z axis, we can choose two light-cone vectors κ µ + and κ µ − such that where from the momentum conservation we have that P ψ⊥ = −p g⊥ , while by considering the softness of the gluon in the final state (q + p a ) 2 ≈ M 2 ψ .We can introduce the variable x g defined by which is also constrained by the momentum conservation (B11) Then, the phase space of the emitted (on-shell) soft gluon is given by

ΛFigure 1 .
Figure 1.Schematical overview of matching in ep → e ′ J/ψ + X to obtain the leading order shape function.
are finite.Despite the different notation, the amplitude squared on the left-hand side of Eq. ( scale µ ′ H and combining it with the proper process-dependent extra soft factor S(µ ′ 2 H ).