$J/\psi$ polarization in semi-inclusive DIS at low and high transverse momentum

We study the polar and azimuthal decay angular distributions of $J/\psi$ mesons produced in semi-inclusive, deep-inelastic electron-proton scattering. For the description of the quarkonium formation mechanism, we adopt the framework of nonrelativistic QCD, with the inclusion of the intermediate color-octet channels that are suppressed at most by a factor $v^4$ in the velocity parameter $v$ relative to the leading color-singlet channel. We put forward factorized expressions for the helicity structure functions in terms of transverse momentum dependent gluon distributions and shape functions, which are valid when the $J/\psi$ transverse momentum is small with respect to the hard scale of the process. By requiring that such expressions correctly match with the collinear factorization results at high transverse momentum, we determine the perturbative tails of the shape functions and find them to be independent of the $J/\psi$ polarization. In particular, we focus on the $\cos 2\phi$ azimuthal decay asymmetry, which originates from the distribution of linearly polarized gluons inside an unpolarized proton. We therefore suggest a novel experiment for the extraction of this so-far unknown parton density that could be performed, in principle, at the future Electron-Ion Collider.

phenomenological analyses within the NRQCD framework are not able to consistently account for all cross sections and polarization measurements for J/ψ mesons produced both in proton-proton and in electron-proton collisions.
In this paper, we show how the study of J/ψ polarization in semi-inclusive, deep-inelastic electron-proton scattering (SIDIS), i.e. e p → e J/ψ X, can shed light on the still open puzzles of quarkonium production mechanism and polarization. Such an analysis is complementary to the one performed in Ref. [14], where the transverse momentum spectrum and the azimuthal distribution for the production of unpolarized J/ψ mesons were considered. A satisfactory theoretical description will be achieved only if all these observables can be simultaneously described within a unique framework. Their experimental determination could be reached at the future Electron-Ion Collider (EIC) [15][16][17], which will be built in the United States.
Along the lines of Ref. [14], in order to avoid contributions from photoproduction processes, we only consider the kinematic region where the virtuality of the photon exchanged in the reaction, Q, is equal or greater than the J/ψ mass, M ψ , namely Q ≥ M ψ . Moreover, we denote by q T the transverse momentum of the photon with respect to the J/ψ and the proton four-momenta. Hence, in the calculation of the polarization parameters, we have to deal with three relevant scales: the above-defined transverse momentum q T , a hard scale µ (to be identified with Q, or M ψ , or any combination of them) and a soft scale (the nonperturbative QCD scale Λ QCD or, alternatively, the proton mass). Depending on the value of q T we can adopt two different factorization frameworks, both of them allowing for a separation of the short-distance from the long-distance contributions to the observables under study.
In the high-q T region, namely q T Λ QCD , in a frame where the J/ψ meson is at rest, the photon transverse momentum is generated by perturbative radiation. NRQCD and collinear factorization can be applied and the resulting polarization parameters will depend on collinear (i.e. integrated over transverse momentum) parton distribution functions (PDFs) and LDMEs. In the small-q T region, q T µ, the photon transverse momentum is nonperturbative instead, and transverse momentum dependent (TMD) factorization [18][19][20] is expected to be appropriate. Observables should then depend on TMD PDFs (or TMDs for short) and shape functions [21,22], which are the generalization of the NRQCD LDMEs. Alternatively, the latter can be seen as the analog of the TMD fragmentation functions for light hadrons. We note that in the overlapping region, Λ QCD q T µ, both frameworks can be applied and the results obtained within the two formalisms have to match, provided they describe the same underlying mechanism [23]. This property has been proven for several observables for which TMD factorization at the twist-two level has been demonstrated, such as the unpolarized cross sections, differential in q T and integrated over the azimuthal angles of the final particles, for the SIDIS process e p → e h X, where h is a light hadron, and for Drell-Yan (DY) dilepton production, p p → X [24,25]. Although there is not yet a rigorous proof of TMD factorization for the J/ψ polarization parameters in SIDIS, strong arguments exist in favor of its validity, if we consider the analogy with e p → e h X, for which TMD factorization holds at all orders [18]. As already pointed out in Ref. [14], these two processes are essentially equivalent from the point of view of the color flow, which determines the gauge-link structure of the TMD parton correlators [26]. Since neither the J/ψ mass nor its spin can affect such structure, we do not expect any factorization breaking effects due to color entanglement. We therefore propose reasonable factorized expressions for the J/ψ polarization parameters in terms of twist-two TMDs and shape functions, which properly match with the collinear results in the intermediate region Λ QCD q T µ. Finally, we point out that the theoretical TMD framework we have devised will also have important implications for the extraction of gluon TMDs at the EIC. Indeed, although several proposals for the extraction of these distributions have been put forward, both in ep [26][27][28][29] and in pp collisions [27,[30][31][32][33], as well as within the more phenomenological generalized parton model approach [34][35][36][37][38][39][40][41], they are still basically unkown. In particular, we find that the so-called ν parameter, related to the cos 2φ azimuthal asymmetry of the leptons from the J/ψ decay, can give direct access to the distribution of linearly polarized gluons inside unpolarized protons.
The paper is organized as follows. In Section II we review in detail the main properties of the cross section using only kinematic considerations and the symmetries of the strong and electromagnetic interactions, without referring to any specific model concerning quarkonium formation. In particular, the well-known result of the angular structure of the cross section, expressed in terms of four independent helicity structure functions, is obtained along the lines of the derivation presented in Ref. [42] for the DY process. In Section III we compute the polarized structure functions within the framework of collinear factorization and NRQCD at the order α 2 s . Moreover, we investigate their small-q T limit and show that the dominant terms do not depend on the choice of the coordinate axes. Section IV is devoted to the calculation of the structure functions in the TMD regime at the order α s and to the study of their matching with the collinear results in the common region of validity. Summary and conclusions are gathered in Section V. Details on the reference frames and on the transformations which connect them can be found in Appendix A. Finally, the explicit expressions for the partonic helicity structure functions in the Gottfried-Jackson frame are collected in Appendix B.

II. MODEL INDEPENDENT PROPERTIES OF THE CROSS SECTION
We study the process with X being an undetected hadronic system, and the subsequent leptonic decay where the four-momenta of the particles are given within brackets. Throughout this work the mass of the leptons, both in the initial and in the final state, are neglected, while we denote by M p and M ψ the proton and the J/ψ masses, respectively. Moreover, we will sum over the lepton and J/ψ helicities. The virtual photon exchanged in the reaction carries four-momentum with Q 2 ≡ −q 2 > 0. In the deep-inelastic limit, where Q 2 , P · q and P ψ · q are large while the variables are finite (with values between zero and one), the square of the invariant amplitude for this process can be split into a leptonic part and a purely hadronic part, where the leptonic tensor L µν reads and H µν is the product of the hadronic current matrix elements, The corresponding cross section reads where S = (k + P ) 2 ≈ 2 k · P if also the proton mass is neglected, while B is the branching ratio for the decay process J/ψ → + − . The integration of H µν over P X leads to the usual hadronic tensor Moreover, it is convenient to express the phase space of the final leptons as follows where φ k is the azimuthal angle of the scattered electron, while Ω = (θ, ϕ) refers to the solid angle of lepton + in a reference frame where the system formed by + and − is at rest. Hence, after integration over φ k , Eq. (8) can be cast in the form The hadronic tensor W µν is a function of the four-vectors q µ , P µ , P µ ψ and contains the information on the proton structure as well as on the J/ψ formation process and polarization. It has to fulfill certain constraints due to electromagnetic gauge invariance, parity, and hermiticity. More specifically, the gauge-invariance condition, is implied by conservation of the hadronic current J µ (0) in Eq. (7) and limits W µν to be a tensor in the threedimensional space orthogonal to q µ . Moreover, as a consequence of the hermiticity condition, the symmetric part of the hadronic tensor in the indices (µ, ν) is real, while its antisymmetric part is imaginary. Since the leptonic tensor in Eq. (6) is symmetric, only the symmetric part of W µν contributes to the cross section. Finally, parity conservation implies . By imposing Eqs. (12) and (14), the hadronic tensor can be decomposed into a set of four basis tensors multiplied by scalar functions, the so-called structure functions, which have to be real because of Eq. (13). To this end, it is useful to introduce the projector operator to the space orthogonal to q µ , satisfying the relation g µν q µ = g µν q ν = 0 and such that, when contracted with any four-vector, yields a vector orthogonal to q µ . In particular, with P · q = P ψ · q = 0. Hence, the hadronic tensor can be expressed as in terms of the common Lorentz invariant structure functions W i = W i (Q 2 , q · P, q · P ψ ), with i = 1, 2, 3, 4. An alternative parametrization of the hadronic tensor is provided by four, frame dependent, helicity structure functions, which, as already pointed out in Ref. [42] for the DY process, have the advantage of allowing for an explicit factorization of the leptonic and hadronic variables in the cross section. In the following, we will show that this is the case also for the process under study.
As a first step of our derivation, it might be useful to write separately the contributions coming from transversely and longitudinally polarized virtual photons. This can be achieved by expressing the leptonic tensor in Eq. (6) in the form which is valid upon integration over the azimuthal angle of the electron scattering plane, taken with respect to the plane formed by the proton and the virtual photon, and neglecting the proton mass. In Eq. (18), we have introduced the projector to the space orthogonal to both q and P g µν and the longitudinal polarization vector of the exchanged photon, which fulfills the relations 2 γ (q) = 1 and µ γ (q) q µ = 0. Hence we obtain with the definitions which in turn can be expressed in terms of either the invariant structure functions or the helicity ones. It still remains to clarify the dependence of the cross section on the angles of the decaying lepton + . To this end, we write the sum over the final J/ψ helicities λ, λ = −1, 0, 1 in an explicit form, where the superscript P =⊥, refers to the photon polarization and we have introduced the frame-dependent helicity structure functions with µ λ (P ψ ) being the polarization vectors of the spin-1 J/ψ meson defined with respect to a covariant set of coordinate axes (T µ , X µ , Y µ , Z µ ) with normalizations T 2 = 1 and X 2 = Y 2 = Z 2 = −1. These four-vectors can be defined as linear combinations of the physical momenta q µ , P µ , P µ ψ , in such a way that X µ , Y µ , Z µ become three-vectors in the quarkonium rest frame [42,45]. An illustration of the commonly used frames is given in Fig. 1. In all of them, with ε 0123 = +1, while they differ for the choice of X µ and Z µ , see also Appendix A for details. Once the frame is fixed, the polarization vectors are given by µ 0 (P ψ ) = Z µ = (0, 0, 0, 1) , and fulfill the orthogonality and completeness relations, * The constraints of parity conservation and hermiticity imposed by QCD on the hadronic tensor imply the following relations for the helicity structure functions We can therefore decompose W αβ and W ⊥ αβ as [42,44] in terms of the eight independent helicity structure functions where the subscripts refer to the J/ψ polarization: W P T and W P L are respectively the structure functions for transversely and longitudinally polarized J/ψ mesons, W P ∆ are the single-helicity flip structure functions, while W P ∆∆ are the double-helicity flip ones.
A further simplification occurs because the leptonic decay conserves helicity in the massless limit we are considering, i.e. the lepton spins do not flip in the coupling with the J/ψ meson. In the quarkonium rest frame, where + and − are produced back-to-back along the direction identified by the three-vector l = (sin θ cos ϕ, sin θ sin ϕ, cos θ) , the component along this direction of the total angular momentum of the dilepton system, and therefore of the J/ψ, can be either +1 or −1, but not zero. Hence, the sum in Eq. (23) can be performed by choosing a specific set of polarization vectors α σ (P ψ ), where σ is the J/ψ helicity along the direction l, and the transverse projector with respect to the directions of the four-vectors P ψ and l is given by By contracting the above expression for g αβ l⊥ with the one for W P αβ in Eq. (29), according to Eq. (32) we obtain Hence, by using this result and the one in Eq. (21), the angular structure of the cross section in Eq. (11) can be made explicit. We find with the definitions We point out that Eq. (35) represents the main result of this section and shows, as anticipated above, that the helicity structure functions defined in Eqs. (30) and (36) allow for an explicit factorization of the leptonic and hadronic variables in the cross section. It is also convenient to introduce the following ratio of differential cross sections [44], which is therefore given by Alternatively, this ratio can be expressed in terms of the largely used polarization parameters λ, µ, ν, where, by comparison with Eq. (38), we have Notice that Eqs. (38) and (39) describe the typical angular distribution of the leptons that originate from the decay of a spin-one particle in its rest frame, for instance a virtual photon produced in the Drell-Yan process [17] or a J/ψ meson produced in photon-proton collisions [45].

III. HELICITY STRUCTURE FUNCTIONS WITHIN NRQCD AND COLLINEAR FACTORIZATION
The helicity structure functions W P λλ can be calculated within NRQCD and collinear factorization in the kinematic region q T Λ QCD . One first needs to evaluate the partonic structure functions w P(a) λλ for each of the underlying hard scattering subprocesses, which at the order αα 2 s are where a is either a gluon, a quark or an antiquark with helicity λ a , and the charm-anticharm quark pair is produced in the intermediate Fock state J , with S, L, J being the spin, orbital and total angular momenta of the cc pair, respectively, while c = 1, 8 specifies its color configuration. In addition to the leading color-singlet production channel with n = 3 S 1 , we also include the subleading color-octet channels that are relatively suppressed by at most a factor of v 4 , namely 1 S J , with J = 0, 1, 2. The corresponding Feynman diagrams are shown in Fig. 2. The polarized hadronic structure functions are then obtained as convolutions of the partonic ones with suitable parton distribution functions f a 1 . If we denote by M µα (γ * a → J/ψ X) the amplitude for the processes γ * a → J/ψ X, where we assume that the J/ψ meson is produced in a definite helicity state λ, the partonic structure functions are given by with P =⊥, ; the tensors ε µν P being defined as ε µν ⊥ = −g µν ⊥ and ε µν = µ γ ν γ , see also Eqs. (19)- (20). The J/ψ polarization vectors α λ (P ψ ) are presented in Appendix A in terms of the four-momenta q, P , P ψ in the four different frames considered.
According to NRQCD, the J/ψ polarization state λ can be reached from the cc pairs produced in various orbital and spin angular momentum states in the subprocesses of Eq. (41) [45]. Because of charge and parity conjugation, there are no interference effects between intermediate states with different orbital angular momentum L and spin S. However, NRQCD does not forbid interferences of different 3 P J states, for which S = L = 1 [46][47][48]. The following decomposition of w which expresses the coherent sum of the partonic structure functions for the four intermediate cc states. Since the 1 S [8] 0 state is rotational invariant in the quarkonium rest frame, it leads to a random orientation of the J/ψ spin. In other words, each J/ψ helicity state λ amounts to one third of the unpolarized cross section. Furthermore, for the states, NRQCD spin symmetry implies that the perturbatively calculable third component of the spin of the cc pair, namely the quantum number S z , is not modified in the subsequent nonperturbative evolution into the observed J/ψ meson. Hence, S z is equal to the helicity of the J/ψ, i.e. S z = λ. The only nontrivial term of the sum in Eq. (43) Representative diagrams of the partonic subprocesses contributing to J/ψ production in SIDIS at the order αα 2 s , γ * (q) + a(pa) → J/ψ(P ψ ) + a(p a ) with a = g, q,q. The six diagrams of type (a) are the only ones corresponding to the CS production mechanism. Moreover, there are two diagrams for each type (b), (c), (d).
is the last one. This can be calculated explicitly by projecting the hard scattering amplitude onto states of definite S z = λ and L z , squaring the amplitude and summing over L z [45][46][47]. In formulae, where we have denoted the quantum numbers of the cc pair by (L, L z ; S, S z ). At the order in v we are considering and by further exploiting the symmetries of NRQCD, it can be shown that the partonic structure functions can be written in a factorized form, in which the nonperturbative contributions, encoded in the same four LDMEs 0|O(n)|0 that describe unpolarized J/ψ production [14,[50][51][52][53], are multiplied by perturbative short-distant coefficients K P(a) αβ [n]. These can be calculated from the hard scattering partonic subprocesses γ * a → cc[n] a. Hence we can write where the sums run over the Fock states n = 3 S 0 . In the actual calculation we neglect the proton mass and any smearing effects both in the initial and in the final states, hence partons can be taken to be collinear to the parent proton, p µ a = ξP µ . The w P(a) λλ can therefore be expressed as functions of the scalar products of the four-vectors q, p a , P ψ , or, equivalently, of the usual partonic Mandelstam variablesŝ It is then straightforward to obtain the structure functions w in analogy to Eq. (30). Their expressions in the GJ frame are collected in Appendix B, while the corresponding formulae in the other three reference frames considered (TF, CS, HX) can be obtained by applying the transformations listed in Appendix A. We note that in the photoproduction limit (Q → 0) our results agree with the ones in Ref. [45]. Furthermore, we are able to reproduce the SIDIS unpolarized cross section, related to the combination 2W T + W L , presented in Ref. [50]. The structure functions W T and W L for J/ψ production in SIDIS have been calculated also in Ref. [49], however we are in disagreement with the relative normalizations of most of the contributions to the structure functions for longitudinally polarized virtual photons, w (g) Λ [n], presented in Ref. [49]. Since w (g) Λ [n] → 0 as Q → 0, this might explain why the authors of Ref. [49] recover the photoproduction results in Ref. [45], but not the unpolarized SIDIS cross section in Ref. [50]. Conversely, to the best of our knowledge, the SIDIS helicity flip structure functions W ∆ and W ∆∆ have been explicitly computed for the first time in the present work.
In order to write down the expressions for the hadronic structure functions, it is convenient to introduce the scaling variablesx By comparison with the hadronic variables in Eq. (4), we obtain the relationsẑ = z andx = x B /ξ. The latter impliesxp a = x B P , withx ≥ x B . These variables allow us to perform, in the J/ψ rest frame, the following Sudakov decomposition of the four-momenta of the particles involved in the reaction, where n + and n − are two light-like vectors such that n + · n − = 1 and q 2 T = −q 2 T . Hence the Mandelstam variables in Eq. (46) can be written aŝ and the hadronic structure functions are given by the convolutions where the final result has been obtained by substituting the expression for w P(a) λλ given in Eq. (45), whilê Moreover, in Eq. (50) µ 2 is the hard factorization scale, on which also the partonic structure functions w P(a) λλ [n] depend, even if not explicitly shown.
The helicity structure functions discussed so far are expected to describe the azimuthal decay distributions of J/ψ mesons for large values of the transverse momentum, namely for q 2 T Λ 2 QCD . Their behavior in the small-q T region, Λ 2 QCD q 2 T Q 2 , can be obtained along the lines of Ref. [14], by replacing the Dirac delta in Eq. (50) with its expansion in the small-q T limit, wherex =x/x max . Up to corrections of the order of O(Λ QCD /|q T |) and O(|q T |/Q), we find the following leading power behavior of the structure functions for transversely polarized photons where a sum over i = q,q is understood. The single-helicity flip structure function W ⊥ ∆ is suppressed by a factor |q T |/Q relatively to the ones in Eq. (53) and therefore we do not consider it in the following. Similarly, for longitudinally polarized photons, while all the other structure functions are suppressed. We note that the above results are independent of the chosen reference frame. The quantities w P (g) Λ in Eqs. (53) and (54) are the gluon helicity structure functions at LO in α s , i.e. they refer to the 2 → 1 subprocess γ * g → J/ψ in Fig. 3 and read where, as before, the superscripts (subscripts) refer to the photon (J/ψ) polarization states and e c is the electric charge of the charm quark in units of the proton charge. Furthermore, in Eqs. (53) and (54) we have used the definitions and where T R = 1/2, n f is the number of active flavors and C A = N c , with N c being the number of colors. The symbol ⊗ stands for a convolution in the longitudinal momentum fractions of the splitting functions and the parton distributions, The leading-order (LO) unpolarized splitting functions are explicitly given by with C F = (N 2 c − 1)/2N c . As well known, the plus-prescription on the singular parts of the splitting function P gg is defined so that the integral of a sufficiently smooth distribution G is given by Moreover, the splitting functions of an unpolarized parton into a linearly polarized gluon read [54,55] δP Finally, we point out that the partonic subprocesses which contribute to the structure functions at small q T , see Eqs.  Fig. 2 (a) and the quark-induced subprocesses in Fig. 2 (c), are suppressed and vanish in the limit q 2 T → 0. This means that they are not relevant in our analysis of the matching of the collinear and TMD results. Hence they will not be considered in the following discussion.

IV. TMD FACTORIZATION AND MATCHING
Within the TMD framework at the order αα s in the strong coupling constant and v 4 in the NRQCD velocity parameter, the underlying partonic subprocesses are given by where the cc pairs are produced perturbatively in one of the color-octet bound states n = 1 S [8] 0 , 3 P

[8]
J with J = 0, 1, 2, see also Fig. 3. The leading-twist helicity structure functions for incoming transversely polarized photons read while for incoming longitudinally polarized photons we find where we have introduced the transverse momentum convolutions In the above equations, f g 1 and h ⊥ g 1 are, respectively, the unpolarized and linearly polarized gluon TMDs inside an unpolarized proton [56][57][58][59][60][61], while ∆ In absence of smearing effects, the final J/ψ meson would be collinear to the cc pair originally produced in the hard scattering process and the shape functions would be given by ∆ Λ (k 2 T ; µ 2 ) = 0|O(n) |0 δ 2 (k T ). Hence the convolutions in Eq. (65) reduce to the products of a LDME 0|O(n) |0 with a gluon TMD, namely TMDs can be calculated perturbatively in the limit |q T | Λ QCD , where their soft parts can be safely neglected. In particular, the unpolarized gluon TMD distribution at the scale µ 2 and order α s can be expressed in terms of collinear parton distributions as follows [14], where i = q,q, the LO splitting functions P gg and P gi are given in Eq. (59) and the convolution ⊗ is defined in Eq. (58). Similarly, the perturbative tail of the linearly polarized gluon distribution reads [54] with the LO polarized splitting functions δP gg and δP gi given in Eq. (61). By substituting Eqs. (67) together with the expressions for f g 1 (x, q 2 T ; µ 2 ) and h ⊥g 1 (x, q 2 T ; µ 2 ) in Eqs. (68)-(69), evaluated at the scale µ 2 = Q 2 + M 2 ψ , into Eqs. (63)-(64), we find that, if smearing is neglected in the final state, only the double-helicity flip structure function W ⊥ ∆∆ exactly matches the corresponding collinear one in the small-q T limit given in Eqs. (53), (55). This is possible without the need of any shape function because of the absence of a logarithmic term in W T ∆∆ at the perturbative order we are considering. Conversely, in the intermediate region This suggests that smearing effects need to be taken into account in those helicity structure functions which depend on the unpolarized gluon TMD distribution, through the inclusion of a suitable shape function [21,22]. Along the lines of the calculation for unpolarized J/ψ production presented in Ref. [14], to which we refer for details, we are able to find the perturbative tails of the shape functions ∆ L , i.e. they are independent of the J/ψ polarization and equal to the shape function ∆ [n] calculated for unpolarized J/ψ production in Ref. [14] by applying the same matching procedure, This can be checked directly by plugging into the first line of Eq. (65) the full perturbative expansions of the gluon TMD and the shape function up to the order α s , which are obtained by adding the order-α 0 s terms to the expressions in Eqs. (68) and (71), respectively. We get Hence, in the limit q 2 If we substitute this convolution, with the choice µ 2 = Q 2 + M 2 ψ , into Eqs. (63)-(64), we recover the correct formulae for W ⊥ T , W ⊥ L and W L in Eqs. (53)- (54). In this way we are able to solve the matching issue of the TMD and collinear results for all the structure functions in the region Λ 2 QCD q 2 T Q 2 . We note that, at the perturbative order we are considering, nothing can be said about the tail of the shape function ∆ [n] ∆∆ , which could be in principle different from ∆ [n] . Its determination would require a study of the matching of the helicity structure function W ⊥ ∆∆ at higher order in α s , because the perturbative expansion of h ⊥ g 1 starts at the order α s and not α 0 s like for f g 1 . However, a full calculation of the cross section for J/ψ production in SIDIS at the order α 2 α 3 s within the NRQCD framework, is still missing.

V. SUMMARY AND CONCLUSIONS
We conclude by summarizing our main results. We have calculated the structure functions W T , W L , W ∆ , W ∆∆ , and the corresponding polarization parameters λ, µ, ν, for J/ψ production in SIDIS. In particular, the helicity-flip structure functions W ∆ and W ∆∆ , which generate, respectively, cos φ and cos 2φ azimuthal asymmetries in the angular distribution of the decaying leptons, have been derived here for the first time. In general, the explicit expressions depend on the choice of the coordinate axes. We have presented our results in the Gottfried-Jackson frame and we have shown how to transform them to the Target, Collins-Soper and Helicity frames.
In the kinematic region where the transverse momentum q T of the virtual photon exchanged in the reaction is large, q T Λ QCD , the theoretical frameworks of collinear factorization and NRQCD have been adopted. By performing an analysis up to the order α 2 s accuracy, with the inclusion of CO contributions up to the order v 4 with respect to the CS one, we obtain the small-q T behavior of the structure functions W T , W L and W ∆∆ in terms of the LO DGLAP splitting functions. Furthermore, W T and W L receive also large logarithmic corrections, which are absent in W ∆∆ and need to be resummed. In contrast, the structure function W ∆ turns out to be suppressed by a factor q T /Q with respect to the other ones. Quite remarkably, these features at small q T do not depend on the choice of the reference frame.
The large logarithms of the structure functions at small q T can be resummed within the TMD framework. Although a rigorous proof of TMD factorization only exists for light-hadron production in SIDIS, it is reasonable to assume its validity also for the production of heavy quarkonium states, just changing from the fundamental to the adjoint representation in the gauge-link structure underlying the process [26]. Hence we propose a factorized expression of the polarized structure functions in terms of transverse momentum dependent parton distributions and shape functions, the latter being a generalization of the NRQCD LDMEs to TMD approach first introduced in Refs. [21,22]. Our assumption has to fulfill the consistency condition that both descriptions match in the intermediate region Λ QCD q T µ, where µ is a hard scale typical of the process. At the order α s , we find that a smooth transition from low q T to high q T is possibile if we choose µ = Q 2 + M 2 ψ . Moreover, while W ∆∆ matches without depending on any shape function because of the absence of logarithmic divergences, from the matching of W L and W T we are able to deduce the specific form of the color-octet shape functions at large transverse momentum. Very interestingly, the perturbative tails of the shape functions turn out to be independent of the J/ψ polarization, as well as of the quantum numbers of the intermediate color-octet Fock states, except for their overall magnitude given by the NRQCD LDMEs. These conclusions are in full agreement with the ones obtained in the analysis of unpolarized J/ψ production in SIDIS [14]. They of course hold for any other quarkonium state with the same quantum numbers of the J/ψ meson, such as the ψ(2S) and Υ(nS) states.
Finally, our findings show that a combined analysis of the J/ψ polarization parameters and production yields in the large-q T region will be very helpful in testing NRQCD and will likely improve our knowledge of the CO LDMEs. A similar study in the small-q T region will provide important information on the shape functions, especially on their nonperturbative part which we cannot predict. Moreover, we suggest a novel experiment for the extraction of the distribution of linearly polarized gluons inside an unpolarized proton, by looking at the double helicity-flip parameter ν, related to the cos 2φ azimuthal decay asymmetry of the J/ψ meson. The high-luminosity Electron-Ion Collider to be build in the U.S. would be the ideal facility for carrying out this proposal.
In this appendix, we present the explicit expressions of the covariant four-vectors (T µ , X µ , Y µ , Z µ ), which form a set of coordinate axes in the four dimensional Minkowski space. These vectors are given as linear combinations of the physical momenta q µ , P µ , P µ ψ in four typical reference frames, such that X µ , Y µ , Z µ become three-vectors in the quarkonium rest frame [42,45]. As already mentioned in Section III, the normalization conditions are given by T 2 = 1 and X 2 = Y 2 = Z 2 = −1, with T µ and Y µ fixed by Eq. (25) in all frames, namely T µ = P µ ψ /M ψ , Y µ = ε µναβ X ν Z α T β . On the other hand, the definition of X µ and Z µ is not unique. For their determination, we proceed along the lines of Ref. [45], introducing the following auxiliary four-vectors Sinceq · T =P · T = 0, X µ and Z µ can be easily written as a linear combination ofq µ andP µ . One starts by writing with one of the two coefficients fixed by the condition Z 2 = −1. Similarly, where A x and B x can be determined by imposing X 2 = −1 and X · Z = 0. The remaining sign ambiguity in A x and B x is fixed by requiring that the three-vector Y points in the direction of q × (−P ) in the J/ψ rest frame. The four commonly used reference frames described in the following are specified by the choice of Z µ . The covariant expressions for X µ and Z µ are given in terms of the invariants s = 2q·p a , t = M 2 ψ −Q 2 −2q·P ψ , u = M 2 ψ +Q 2 −2p a ·P ψ , with p a = ξP in the collinear approach. These invariants are related to the partonic Mandelstam variables in Eq. (46) as follows and fulfill the equation 1. Target frame. The target frame (TF) is defined in such a way that the vector Z points along the opposite direction of the proton momentum in the J/ψ rest frame, i.e. Z = −P /|P |. Hence the four-vectors Z µ and X µ are given by 2. Collins-Soper frame. In the Collins-Soper frame (CS), the direction of Z is chosen as the bisector of the angle between the proton and the virtual photon three-momenta in the J/ψ rest frame, Z = q/|q| − P /|P |. Therefore, with 3. Gottfried-Jackson frame. In the Gottfried-Jackson frame (GJ), also known as u-channel helicity frame, Z is chosen to be along the direction of the three-momentum of the virtual photon in the J/ψ rest frame. Namely, Z = q/|q| and with 4. Helicity frame. In the helicity (HX) or recoil frame, the Z-axis is chosen to be along the direction of the J/ψ three-momentum in the hadronic center of mass frame, that is Z = −(P + q)/|P + q| in the J/ψ rest frame. Moreover, with The transition from one frame to the other is given by a rotation around the Y -axis, where, from the GJ to the TF frame: from the GJ to the CS frame: form the GJ to the NX frame: As a consequence, the hadronic structure functions W Λ (as well as the partonic ones) in two different frames F and F are connected through the linear transformation Similarly, for the polarization parameters we have Appendix B: Partonic helicity structure functions We collect here the analytic expressions for the partonic helicity structure functions in the GJ frame. For each subprocess initiated by a parton a, and for each Fock state n, they are given in the form where, as already stated, P =⊥, refers to the virtual photon polarization, while Λ = T, L, ∆, ∆∆ to the J/ψ one. Moreover, F [n] is a common prefactor independent of the polarizations. In order to render the formulae more compact, they are given in terms of the Lorentz invariants s, t, u in Eq. (A4). Moreover, for simplicity we will drop the superscript (a) on w P (a) Λ , as well as the dependences of F and w P Λ on n when not needed.