Extracting Color Octet NRQCD matrix elements from $J/\psi$ production at the EIC

Recently unpolarized and polarized $J/\psi \,(\Upsilon)$ production at the Electron-Ion Collider (EIC) has been proposed as a new way to extract two poorly known color-octet NRQCD long-distance matrix elements: $\langle0\vert{\cal O}_{8}^{J/\psi}(^{1}S_{0})\vert0\rangle$ and $\langle0\vert{\cal O}_{8}^{J/\psi}(^{3}P_{0})\vert0\rangle$. The proposed method is based on a comparison to open heavy-quark pair production ideally performed at the same kinematics. In this paper we analyze this proposal in more detail and provide predictions for the EIC based on the available determinations of the color-octet matrix elements. We also propose two additional methods that do not require comparison to open heavy-quark pair production.

( 1 S 0 )|0 and 0|O J/ψ 8 ( 3 P 0 )|0 . In Ref. [4], taking ratios with spin-dependent asymmetries was suggested to cancel out these quantities. It was also pointed out in Ref. [4] that ratios with the analogous expressions for unpolarized open heavy-quark pair production could be used to cancel out the gluon TMDs so as to determine the CO matrix elements experimentally. This can hopefully help reducing their uncertainty and, as a consequence, forgo the need for comparison to spin asymmetries. In this paper we look into the determination of CO matrix elements in more detail, providing the general expressions and their current estimates, which can be of help in the study of these quantities at the EIC.
Following the calculations of Ref. [4], we study the semi-inclusive deep-inelastic (SIDIS) process of quarkonium production e( ) + p(P ) → e( ) + Q (P Q ) + X , where Q is either a J/ψ or a Υ meson. Besides unpolarized quarkonium production we also examine the cases in which the quarkonium is longitudinally or transversely polarized w.r.t. its direction of motion in the γ * p center-of-mass frame. The reference frame is such that both the exchanged virtual photon γ * and the incoming proton move along theẑ-axis. Azimuthal angles are measured w.r.t. to the lepton scattering plane. In order for TMD factorization to apply, the component of the quarkonium momentum transverse to the lepton plane, i.e. q T ≡ P QT , should be small compared to the virtuality of the photon Q and to the mass of the quarkonium M Q 1 The differential cross section can then be written as dσ [e p → e Q P X] = 1 2s where Γ g is the gluon correlator and L( , q) is the lepton tensor. For details of the intermediate calculation we refer to Ref. [4]. Sticking to the notation of Ref. [4], the differential cross section will be denoted by which depends on the variable z = P · P ψ /P · q, on the inelasticity y = P · q/P · , on Bjorken-x x B = Q 2 /2P · q, and on the transverse momentum q T = P QT of the quarkonium state Q (P 2 QT = −P 2 QT ) with an azimuthal angle φ T w.r.t. the lepton plane (φ = φ = 0). The variables y and x B are related to the total invariant mass squared s = ( + P ) 2 ≈ 2 · P and to the photon virtuality Q 2 = −q 2 ≡ −( − ) 2 through the relation s = 2 P · q/y = Q 2 /x B y. The superscript U P = U U, U L or U T denotes an unpolarized proton and a polarization state P for the quarkonium which can be either unpolarized (U ), longitudinally polarized (L) or transversely polarized (T ) with respect to the direction of its three-momentum in the photon-proton center-of-mass frame.
For the calculation of the cross section, we adopt the TMD factorization framework in combination with nonrelativistic QCD (NRQCD) [9][10][11]. In NRQCD, the nonperturbative hadronization process of the heavy-quark pair into a quarkonium bound state is encoded in long-distance matrix elements (LDMEs) [12]. Observables are evaluated by means of a double expansion in the strong coupling constant α s and in the average velocity v of the heavy quark in the quarkonium rest frame [13], where v 2 0.3 for charmonium and v 2 0.1 for bottomonium. In SIDIS to leading order (LO) in α s , the heavy quark-antiquark pair can only be produced in a color-octet (CO) state. The CO LDMEs are typically determined from fits to data on J/ψ and Υ yields [14][15][16][17][18]. The extracted values are not compatible with each other, even within the large uncertainties. Therefore, any new method to determine the CO matrix elements with better precision is worth exploring. Moreover, the methods to be discussed partly rely on experimental determination of the polarization states of the produced quarkonia, which is not only beneficial for the extraction of the CO matrix elements but also to the understanding of the production mechanism, which is still incomplete. Although NRQCD successfully explains many experimental observations, describing all cross sections and polarization measurements for charmonia in a consistent way still poses challenges [19][20][21].
In the process under study, at LO in the strong coupling constant α s , the QQ pair forms a bound state with spin S, orbital angular momentum L, and total angular momentum J, for which we employ spectroscopic notation: 2S+1 L J . Since to LO accuracy there is no gluon emission in the final state, z = 1. The relevant CO matrix elements for J/ψ production are 0|O 1 state together with a gluon, hence 0 ≤ z ≤ 1. The CS contribution is, therefore, suppressed relatively to the CO one by a perturbative coefficient of the order α s /π [22,23] and vanishes at z = 1. On the other hand, 0|O by v 3 and v 4 , respectively. As a result the CO contribution should be enhanced by about a factor v 3 π/α s ≈ 2 w.r.t. the CS contribution. In the analysis of Ref. [22] this factor turned out to be larger (≈ 4) for values of Q 2 > 4 GeV 2 . Based on these findings, we restrict our study to the CO contributions at leading order in α s .
The cross section for e p → e Q P X in Eq. (3) can be cast in the following form: where the unpolarized gluon TMD f g 1 and the linearly polarized gluon TMD h ⊥ g 1 depend, besides on the transverse momentum, also on the momentum fraction x given by The normalization factor N reads with e Q denoting the fractional electric charge of the quark Q. The expressions for A U P and B U P are where the subscripts U + L, L, T refer to the specific polarization of the photon [24,25]. Using heavy-quark spin symmetry relations [12] 0|O J/ψ 8 one obtains to leading order in v: Similarly, the expressions for longitudinally polarized quarkonium production are Since dσ U U = dσ U L + dσ U T , this determines the corresponding expressions for transverse polarization of the quarkonium: One can isolate the angular independent term of the cross section and the azimuthal modulation by (weighted) integration: The hard scales of the process e p → e Q X are the photon virtuality Q and the quarkonium mass M Q ≈ 2M Q . To avoid ratios of the two it is convenient to choose Q = M Q , but there are also advantages to varying Q. In the case of P = U and writing Q 2 = cM 2 Q we find that (as a function of the heavy-quark mass M Q rather than the quarkonium mass M Q ) where These two long-distance matrix elements have been determined with quite large uncertainty from fits to data (see below), but in all cases Using the Q 2 -and, to a lesser extent, y-dependence offers one method to extract the two matrix elements or at least their ratio if the gluon TMDs are not yet known. In A U U , the prefactor multiplying O P 8 is largest when y is smallest, and has a stronger y-and c-dependence than the one multiplying O S 8 . Considering larger Q 2 will reduce the contribution from O P 8 . More specifically, for y = 0.1 and c = 1 the prefactor multiplying O P 8 is approximately 5 times the one multiplying O S 8 , whereas for large c it reduces to a factor 3 independently of y. In B U U , the contribution from O P 8 is completely independent of y and flips sign when c = 3. Therefore, especially the Q 2 -dependence could be exploited to change the relative contributions from O S 8 and O P 8 to the terms in Eqs. (19) and (20).
If one is able to determine the polarization of the quarkonium state, then there is another method to extract the two LDMEs. For longitudinally polarized quarkonium P = L one obtains Compared to A U U , the O P 8 term in A U L has a different y and c dependence, which has quite important implications as we will discuss now. For y = 0.1 and c = 1 the prefactor multiplying O P 8 in A U L is around 9 times larger than the one multiplying O S 8 , whereas for large c it reduces to a factor 3 independently of y. This means that the prefactor by as much as an order of magnitude), but it can also lead to a significant deviation of A U L from A U U /3 (and of A U T from 2A U U /3), signaling the production of polarized J/ψ mesons. This would be in contradiction to the recent conclusion of Ref. [23], i.e. that the dominance of the contribution of the 1 S [8] 0 cc state to the cross section for the electron-hadron scattering process e h → J/ψ X implies that the produced J/ψ meson will likely be unpolarized. This is presented as a robust test of NRQCD factorization and as a way to shed light on the J/ψ production mechanism. As opposed to this statement from a collinear factorization analysis, our TMD factorization study indicates that the production of polarized J/ψ mesons is not necessarily in contradiction with NRQCD, as it could simply signal a kinematic enhancement factor in front of the matrix element O P 8 . For B U L there is no c-dependence, only the overall y-dependence of 1 − y. Analogous expressions can be obtained when the J/ψ meson is transversely polarized, namely For c = 1, we obtain for the production of unpolarized quarkonium states, from Eqs. (19), (20), (21), and (22): while for longitudinally polarized quarkonium states, from Eqs. (19), (20), (23), and (24): To our present LO accuracy in α s , the ratios D Q L /D Q U and N Q L /N Q U are independent of the TMDs: These ratios offer another method to extract the two CO matrix elements if the polarization of the quarkonium state can be determined. Analogous ratios can be defined for transversely polarized quarkonium states, see Eqs. (19), (20), (25), and (26). Notice that N Q T depends only on one of the two matrix elements, namely O S 8 . Again, the ratios D Q T /D Q U and N Q T /N Q U are independent of the TMDs, as are the ratios D Q T /D Q L and The third and last method to be discussed is the one first suggested in Ref. [4]. For this, we will compare the above expressions to the analogous ones for the process e p → e Q Q X, where Q is either a charm or bottom quark. The differential cross section for the process in which the quark-antiquark pair is almost back to back in the plane orthogonal to the direction of the proton and the exchanged virtual photon, is written as [24,26,27] In the γ * p center-of-mass frame, the difference of the transverse momenta of the outgoing quark and antiquark, K ⊥ ≡ (K Q⊥ − K Q⊥ )/2, should be much larger than the vector sum q T ≡ K Q⊥ + K Q⊥ . The angles φ T and φ ⊥ are the azimuthal angles of q T and K ⊥ w.r.t. the lepton plane, respectively. Furthermore, z = K Q · P/q · P .
The process e p → e Q Q X depends on three large scales: where In order to avoid the necessity for TMD evolution in the comparison to the quarkonium production case, we fix c = 1, i.e. we take K ⊥ = Q = 2M Q . Assuming furthermore that the energy is approximately equally distributed to the heavy quark and antiquark, we choose z = 1/2. Therefore, from Eqs. (37) and (38), we get [26] The comparison to quarkonium production then yields the following two independent observables Since they appear in two distinct combinations, the LDMEs O S 8 and O P 8 can be extracted from experimental measurements of the above ratios. It should be stressed that the presented expressions are LO in α s and order v 3 in the NRQCD velocity parameter v. Next-to-leading order (NLO) corrections in α s will reintroduce sensitivity to the gluon as well as quark TMDs, which in turn will introduce a q T -dependence in R and R cos 2φ T . Since the LO CO contribution, which we consider here, is expected to dominate over the NLO CS and CO contributions, parametrically by a factor v 3 π/α s ≈ 2 and in practice by a larger factor for high Q 2 [22], we expect this dependence to be small. Another source of q T -dependence in R and R cos 2φ T , however, may come from final state smearing effects. Indeed, in NRQCD the outgoing quarkonium state has the same transverse momentum as the QQ pair, which is of course an idealization. While adding some additional uncertainty, the effects from smearing are expected to be mild [4], entering merely as an overall q T -dependent prefactor which is independent of the orbital angular momentum (L) of the CO state, at least in the perturbative region [28]. Thus, whenever a q T -dependence is observed in R and/or R cos 2φ T , final state smearing effects in terms of shape functions [28][29][30] should be considered next to the higher-order contributions from TMDs. Differences in the q T -dependence of the O S 8 and O P 8 contributions would signal higher order terms or an L-dependence in the shape functions.
In case the quarkonium polarization state can be determined, additional expressions can be obtained which can be used for consistency checks as they involve different combinations of the same CO matrix elements. For longitudinally polarized quarkonium production, one obtains: and for transversely polarized quarkonium production: These expressions satisfy: Note that the measurement of R cos 2φ T T directly probes the matrix element O S 8 = 0|O Q 8 ( 1 S 0 )|0 . Selecting photon virtualities such that Q M Q ensures that this relation is not spoiled by large logarithmic corrections. The polarized quarkonium expressions also allow to assess the impact of higher order contributions, because these are different for the unpolarized and polarized cases and do not cancel in the ratios.
In order to assess the potential of this third method we perform a numerical investigation of the various ratios in the next section.

II. NUMERICAL RESULTS
In Tables I and II we list various available determinations of the CO matrix elements. These have been used in the calculation of the cross section σ CO for the production of unpolarized J/ψ and Υ mesons presented in Tables III and IV, respectively. The results are obtained within the collinear factorization framework, in the kinematic region defined by y = 0.1, Q = 2M Q and √ s = 100 GeV, which imply x B = 0.01 for J/ψ production (M c = 1.5 GeV) and x B = 0.09 for Υ production (M b = 4.88 GeV). Moreover, we have adopted the LO MSTW parametrization for the gluon distribution function [31], with the hard scale taken to be equal to Q = 2M Q . The ratio with the CS background, σ CO /σ CS , is shown as well. Note that all the LDME extractions were performed at NLO accuracy, except the Sharma & Vitev set which is LO. The negative cross section corresponding to the NLO LDMEs of Butenschön & Kniehl is due to the dominant contribution of the negative matrix element O P 8 for the process and the kinematic region under investigation. Although σ CO /σ CS strongly depends on the choice of the J/ψ LDME set, we find that it is possible to suppress the CS contribution by requiring z > 0.9. Such a cut will not affect the CO contribution, for which z = 1 at LO and neglecting smearing effects. Butenschön & Kniehl [14] 4.50 ± 0.72 −0.54 ± 0.16 ×10 −2 GeV 3 1.5 GeV 2 Chao et al. [15] 8.9 ± 0.98 0.56 ± 0.21 ×10 −2 GeV 3 not specified 3 Sharma & Vitev [16] 1.8 ± 0.87 1.8 ± 0.87 ×10 −2 GeV 3 1.4 GeV 4 Bodwin et al. [17] 9.9 ± 2.2 0.49 ± 0.44 ×10 −2 GeV 3 1.5 GeV   The values for σ CO and σ CO /σ CS for J/ψ production, calculated at y = 0.1, Q = 2Mc = 3 GeV and xB = 0.01, with different cuts on the variable z. The results are obtained using the fits for the CO LDMEs given in Table I. For the CS matrix element, we adopted the values 1.32 GeV 3 (fit 1), 1.16 GeV 3 (fit 2), and 1.2 GeV 3 (fits 3, 4). The values for σ CO and σ CO /σ CS for Υ production, calculated at y = 0.1, Q = 2M b = 9.76 GeV [16], and xB = 0.09, with different cuts on the variable z. The results are obtained using the fits for the CO LDMEs given Table II. For the CS matrix element, we adopted the value 10.9 GeV 3 .
The results for R and R cos 2φ T are displayed in Fig. 1, where for fit 2 we assumed M c = 1.5 GeV. Possible correlations between the errors of the two fitted LDMEs have not been taken into account, i.e. the errors are simply added in quadrature. This may imply an overestimate of the uncertainty. From these results it is clear that different fits give quite different results. For instance, fit 1 gives a result compatible with zero for R but not for R cos 2φ T , whereas for fit 3 it is the other way around. These results give the general impression that R cos 2φ T may be significantly larger than R, perhaps even by as much as a factor of 10.
Note that the ratios are not normalized to [0, 1] for R or to [−1, 1] for R cos 2φ T . However, one can relate their ratio  to a ratio of asymmetries that are normalized: where cos 2φ T X = N X /D X for X = Q or X = QQ. From this one can see that even though cos 2φ T Q and cos 2φ T QQ are not known, a rough average of the fits shows that the cos 2φ T asymmetry could be substantially larger in J/ψ production than in open charm production. This would suggest that a study of the linearly polarized gluon TMD may be even more promising in J/ψ production than in open charm production. Indeed, this seems to be supported by studies in a small-x model, where cos 2φ T QQ was found to be at the 5 − 10% level [26] and cos 2φ T Q around 10 − 20% [4], in both cases increasing with increasing values of q T . Since we consider gluon induced processes, smaller x values may be beneficial because the gluon TMDs would be enhanced. On the other hand, nonlinear QCD evolution was seen to have a suppressing effect on the cos 2φ T Q asymmetry as x becomes smaller [4]. At the EIC, the smaller the x value, the smaller the Q values covered, and the smaller the q T -range for which TMD factorization is expected to hold. Therefore, one has to keep a balance between the x-and Q-ranges. Similarly, the results for the P = L and P = T cases are given in Figs. 2 and 3. The figures show only a moderate y-dependence, most pronounced for large y-values. Note that although R cos 2φ T P has no y-dependence, both numerator and denominator have a prefactor (1 − y), hence vanish at y = 1. It is also worth remarking that, for the Υ meson, the error bars for R cos 2φ T T are significantly smaller than those for R cos 2φ T L and R cos 2φ T . This is due to the fact that in the P = T case, this ratio only depends on the LDME O S 8 .

III. SUMMARY AND CONCLUSIONS
In this paper we have presented several methods to use J/ψ and Υ production in unpolarized SIDIS at the EIC to obtain improved determinations of two CO LDMEs, 0|O is independent of the gluon TMDs. The third method is based on the comparison of the process e p → e J/ψ (Υ) X with e p → e D D (B B) X. Estimates based on the available fits of the CO matrix elements were presented for this third method. Since the available fit values are not compatible with each other, any new method to determine the CO matrix elements with better precision will be worth exploring. Despite the large uncertainty the estimates show sizable ratios. These results suggest that the cos 2φ T asymmetry that arises from linearly polarized gluons inside the unpolarized proton could be substantially larger in J/ψ production than in open charm production. For typical EIC kinematics, this is corroborated by the small-x model studies in Refs. [26] and [4], where asymmetries around 5-10% and 10-20% were found for open charm and J/ψ production, respectively.
The experimental determination of the polarization states of the produced quarkonia will not only be very helpful in the determination of the CO matrix elements, it could also help improve the understanding of the quarkonium production mechanism in NRQCD, which still poses challenges [19][20][21]. In contrast to the recent conclusion of Ref. [23] that J/ψ mesons produced in electron-hadron collisions will likely be unpolarized and thus can provide a rigorous test of NRQCD, we find, from our LO TMD factorization analysis, that polarized production of J/ψ mesons is not necessarily a stringent NRQCD prediction.
The robustness of our results is ensured by the following considerations. First, TMD factorization is applicable to SIDIS in the kinematic region where q 2 T Q 2 . Considering an outgoing quarkonium state instead of a light meson is not expected to spoil the factorization, provided one allows for a q 2 T -dependence of the LDMEs (the shape functions, Refs [28][29][30]). However, in Ref. [28] it was shown that for large q 2 T , the final state smearing is independent of the orbital angular momentum of the bound state, hence the same for the S-and P -wave CO states, and thus drops out in the ratios at LO. We have assumed this to hold for all q 2 T Q 2 , but otherwise the restriction M 2 p q 2 T Q 2 could be considered. Second, based on the velocity scaling rules, NLO CS contributions are suppressed by a factor v 3 π/α s ≈ 2 w.r.t. the LO CO contributions which we consider here. In an actual numerical study (in the collinear framework), we showed that this suppression can be even stronger, especially for the LDME sets of Ref. [16], and the CS contributions can be made negligible if one imposes a lower cut on z. In turn, the possible diffractive background that becomes important at high z, as well as potential corrections from higher-twist effects, can be suppressed by looking at sufficiently high virtualities Q 2 [22]. Third, the asymmetries that we propose are always the ratios of two cross sections. Hence, besides being independent of the normalizations of the cross sections, they are also expected to be less sensitive to higher-order corrections or to other sources of uncertainties such as the exact value of the charm and bottom mass or the final state smearing effects. The importance of higher order corrections and final state smearing effects will make themselves apparent in the possible q T -dependence of the various ratios that were presented, and can be included as needed.
We conclude that our findings show that heavy-quark final states at EIC are very promising tools to provide improved determinations of CO LDMEs.