#  $J/\psi$ polarization in large-$P_{T}$ semi-inclusive deep-inelastic
scattering at the EIC

Umberto D’Alesio<EMAIL_ADDRESS>Dipartimento di Fisica,
Università di Cagliari, Cittadella Universitaria, I-09042 Monserrato (CA),
Italy  INFN, Sezione di Cagliari, Cittadella Universitaria, I-09042 Monserrato
(CA), Italy  Luca Maxia<EMAIL_ADDRESS>Dipartimento di Fisica,
Università di Cagliari, Cittadella Universitaria, I-09042 Monserrato (CA),
Italy  INFN, Sezione di Cagliari, Cittadella Universitaria, I-09042 Monserrato
(CA), Italy  Francesco Murgia<EMAIL_ADDRESS>INFN, Sezione di
Cagliari, Cittadella Universitaria, I-09042 Monserrato (CA), Italy  Cristian
Pisano<EMAIL_ADDRESS>Dipartimento di Fisica, Università di
Cagliari, Cittadella Universitaria, I-09042 Monserrato (CA), Italy  INFN,
Sezione di Cagliari, Cittadella Universitaria, I-09042 Monserrato (CA), Italy
Sangem Rajesh<EMAIL_ADDRESS>Department of Physics, School of
Advanced Sciences, Vellore Institute of Technology, Vellore, Tamil Nadu
632014, India  INFN, Sezione di Perugia, via A. Pascoli snc, 06123, Perugia,
Italy

(August 30, 2024)

###### Abstract

We present a detailed phenomenological study of $J/\psi$ polarization in semi-
inclusive deep inelastic scattering processes, focusing on the kinematics
accessible at the future Electron-Ion Collider. We show theoretical estimates
for the standard polarization parameters for different frames usually adopted
in the literature, in the large $P_{T}$ region, namely
$P_{T}\gg\Lambda_{\text{QCD}}$, where collinear factorization is expected to
hold. We adopt both the Color Singlet Model and the Nonrelativistic QCD
approach, paying special attention to the role of different sets of Long
Distance Matrix Elements. Finally we present a preliminary analysis of some
frame independent polarization invariants.

##  I Introduction

Our understanding of the $J/\psi$ production mechanism at high energies has
improved significantly since its discovery almost 50 years ago [1, 2], thanks
to the combined efforts from both the theoretical and experimental
communities. However, there are still major problems in the theoretical
analyses of the available data, such as the long-standing $J/\psi$
polarization puzzle. Namely, $J/\psi$ polarization measurements cannot yet be
explained in a way entirely consistent with the world experimental results for
the unpolarized $J/\psi$ yields.

The present theoretical frameworks all agree in providing a perturbative
description of the creation of the charm quark-antiquark ($c\bar{c}$) pair.
The charm mass $m_{c}$ plays the role of the hard scale, since it is much
larger than the asymptotic scale parameter of QCD, $\Lambda_{\text{QCD}}$.
These approaches nonetheless differ in the treatment of the subsequent
nonperturbative transition to the hadronic bound state. For instance, in the
traditional Color-Singlet Model (CSM) [3] the $c\bar{c}$ pair is produced at
short distances directly with the quantum numbers of the $J/\psi$ meson, i.e.
in a color-singlet (CS) state with spin one and no orbital angular momentum.
This is possible by the emission of an additional hard gluon, which implies
the suppression of the cross section by one power of the strong coupling
constant $\alpha_{s}$. However, the CSM cannot be considered as a complete
theory, since at the next-to-leading order (NLO) $P$-wave quarkonia are
affected by uncanceled infrared singularities.

These singularities are properly removed in the effective field theory
approach of nonrelativistic QCD (NRQCD), based on a rigorous factorization
theorem, which was assumed in the original paper by Bodwin, Braaten, and
Lepage [4], and later explicitly proven to next-to-next-to-leading order
(NNLO) [5]. NRQCD therefore implies a separation of process-dependent short-
distance coefficients, to be calculated perturbatively as expansions in
$\alpha_{s}$, from long-distance matrix elements (LDMEs), which are expected
to be universal and have to be extracted from experiments. Scaling rules [6]
predict each of the LDMEs to scale with a definite power of the relative
velocity $v$ of the heavy quark-antiquark pair in the quarkonium rest frame in
the limit $v\ll 1$. Observables are hence evaluated by means of a double
expansion in $\alpha_{s}$ and in $v$, with $\alpha_{s}\simeq 0.2$ and
$v^{2}\simeq 0.3$ for charmonium states. An essential feature of this approach
is that the $c\bar{c}$ pair at short distance can be produced in any Fock
state $n=\,^{2S+1}L_{J}^{[c]}$ with definite orbital angular momentum $L$,
spin $S$, total angular momentum $J$ and color configuration $c=1,8$. NRQCD
hence predicts the existence of intermediate color-octet (CO) states, which
subsequently evolve into physical, CS quarkonia by the emission of soft
gluons. For $S$-wave quarkonia, the CSM is recovered in the limit $v\to 0$. In
the specific case of $J/\psi$ production, the CSM prediction is based only on
the ${}^{3}S_{1}^{[1]}$ CS state, while NRQCD includes the leading
relativistic corrections as well, which at the relative order ${\cal
O}(v^{4})$ are given by the CO states ${}^{1}S_{0}^{[8]},$
${}^{3}S_{1}^{[8]},$ and ${}^{3}P_{J}^{[8]}$ with $J=0,1,2$.

The values of the CO LDMEs extracted from different fits to data on $J/\psi$
and $\Upsilon$ yields [7, 8, 9, 10, 11] are not compatible with each other,
even within the large uncertainties [12, 13, 14]. Therefore, any new method to
determine them with better precision is worth exploring [15, 16, 17]. In this
paper we propose to look at the $J/\psi$ polarization parameters in semi-
inclusive deep-inelastic scattering (SIDIS), $e\,p\to e^{\prime}\,J/\psi\,X$,
in a kinematic region where the transverse momentum of the $J/\psi$ meson
$P_{T}$ is large, namely $P_{T}\gg\Lambda_{\text{QCD}}$, and collinear
factorization is expected to hold. Analysing SIDIS at finite values of the
exchanged photon virtuality $Q^{2}$ has certain experimental and theoretical
advantages as compared to photoproduction. Namely, as $Q^{2}$ increases
theoretical uncertainties in the different frameworks decrease and resolved
photon contributions are expected to be negligible. Moreover, background from
diffractive $J/\psi$ production is expected to decrease with $Q^{2}$ faster
than the SIDIS cross section. The distinct signature of the scattered lepton
makes the process particularly easy to detect. Clearly, cross sections are
smaller than those expected in the photoproduction case, however, considering
the achievable high luminosities, this study should be feasible at the future
Electron-Ion Collider (EIC) planned in the United States [18, 19, 20].

So far, only a single experimental study of $J/\psi$ polarization in SIDIS has
been performed, by the H1 Collaboration at HERA [21]. Such a measurement is
limited to the polarization parameter $\lambda$ in the helicity frame. This
result turns out to be compatible with the predictions provided in Refs. [22,
23], but it can hardly discriminate among the different models. In analogy
with Refs. [22, 23], our phenomenological analysis has been carried out at the
perturbative order $\alpha_{s}^{2}$, which has to be considered as the state
of the art for these observables. Higher-order effects have been calculated
very recently only for the unpolarized cross section within the CSM [24].
Anyway, we expect these effects (at least in the large $Q^{2}$ region) to be
small for the observables we are investigating, because they are ratios of
cross sections. We point out that our estimates include also the polarization
parameters $\mu$ and $\nu$, not addressed in Refs. [22, 23], which are studied
in different reference frames. Furthermore, we perform a preliminary study of
rotational invariant combinations of these parameters.

The remainder of the paper is organized as follows. In section II we recall
the standard SIDIS variables and collect the expressions of the differential
cross section for quarkonium production and its leptonic decay in terms of the
helicity structure functions and the polarization parameters. In section III
we discuss the three polarization parameters $\lambda$, $\mu$, $\nu$, showing
their estimates in two reference frames and paying special attention to their
energy, $z$ and $P_{T}$ dependences as well as to the impact of the LDME set
adopted. To overcome the intrinsic frame dependence of the polarization
parameters, in section IV we present two classes of the so-called rotational
invariant quantities, and show, as a case of study, some results for one of
them. Finally in section V we gather our conclusions.

##  II Kinematics and Formalism

In this section we provide the main analytic expressions needed to carry out
the phenomenological analysis. For more details and the complete formalism we
refer the reader to Ref. [25]. We consider the SIDIS process

$e(k)+p(P)\to e^{\prime}(k^{\prime})+J/\psi(P_{\psi})+X(P_{X})\,,$ (1)

with the subsequent $J/\psi$ decay into a lepton pair

$J/\psi(P_{\psi})\to l^{+}(l)+l^{-}(l^{\prime})\,,$ (2)

where, in brackets, we have shown the four-momenta of each particle. The
$J/\psi$ meson is produced via the partonic subprocess

$\gamma^{*}(q)+a(p_{a})\to c\bar{c}[n](P_{\psi})+a(p_{a}^{\prime})\,,$ (3)

with $q^{2}=-Q^{2}$ and $P_{\psi}^{2}=M^{2}_{\psi}=(2m_{c})^{2}$. The initial
parton momentum, $p_{a}$, is related to the parent proton one, $P$, as

$p_{a}=\xi P\,.$ (4)

We adopt the following three standard invariant quantities, defined in terms
of the photon and hadron momenta

$x_{\scriptscriptstyle B}={Q^{2}\over 2P\cdot q},\quad y={P\cdot q\over P\cdot
k},\quad z={P\cdot P_{\psi}\over P\cdot q},$ (5)

where $x_{\scriptscriptstyle B}$ is the Bjorken variable, $y$ is the
inelasticity and $z$ is the energy fraction carried out by the $J/\psi$ (in
the proton rest frame). All these variables are constrained in the region
$0\leq x_{\scriptscriptstyle B},y,z\leq 1$ and they are connected to other
kinematical quantities of the system, like the total center-of-mass (cm)
energy $\sqrt{s}$ and the virtual photon-proton cm energy, $W$.

The cross section that describes the $J/\psi$ formation and its decay into a
lepton pair can be written as

${1\over B_{ll}}{{\rm d}\sigma\over{\rm d}x_{B}\,{\rm d}y\,{\rm d}z\,{\rm
d}^{2}\bm{P}_{T}\,{\rm
d}\Omega}=\frac{\alpha}{8\,y\,z\,Q^{2}}\,\frac{3}{8\pi}\,\left[{\cal
W}_{T}(1+\cos^{2}\theta)+{\cal W}_{L}(1-\cos^{2}\theta)+{\cal W}_{\Delta}\sin
2\theta\cos\phi+{\cal W}_{\Delta\Delta}\sin^{2}\theta\cos 2\phi\right]\,,$ (6)

where $\bm{P}_{T}$ is the $J/\psi$ transverse momentum in the cm frame of the
virtual photon and the proton, $B_{ll}$ is the branching ratio for the decay
process $J/\psi\to\ell^{+}\ell^{-}$ and $\Omega(\theta,\phi)$ refers to the
solid angle spanned by the lepton $\ell^{+}$ in a reference frame where the
system formed by $\ell^{+}$ and $\ell^{-}$ is at rest. Moreover, we have
introduced the following helicity structure functions

$\displaystyle{\cal W}_{T}$ $\displaystyle\equiv{\cal W}_{11}={\cal
W}_{-1,-1}\,,$ $\displaystyle{\cal W}_{L}$ $\displaystyle\equiv{\cal
W}_{00}\,,$ $\displaystyle{\cal W}_{\Delta}$
$\displaystyle\equiv\frac{1}{\sqrt{2}}\,({\cal W}_{10}+{\cal
W}_{01})=\sqrt{2}\,\text{Re}\,[{\cal W}_{10}]\,,$ $\displaystyle{\cal
W}_{\Delta\Delta}$ $\displaystyle\equiv{\cal W}_{1,-1}={\cal W}_{-1,1}\,,$ (7)

where the subscripts refer to the $J/\psi$ polarization states. More
specifically, ${\cal W}_{T}$ and ${\cal W}_{L}$ are respectively the structure
functions for transversely and longitudinally polarized $J/\psi$ mesons,
${\cal W}_{\Delta}$ is the single-helicity flip structure function, and ${\cal
W}_{\Delta\Delta}$ is the double-helicity flip one. Notice that in Eq. (6) we
have introduced a proper overall constant factor w.r.t. Eq. (2.35) of Ref.
[25] to ensure the normalization when integrated over the solid angle, see Eq.
(8) below. This does not affect any conclusion of Ref. [25], where all
relevant quantities are defined as ratios of helicity structure functions.

As shown in Ref. [25], the structure functions in Eq. (7) can be further
decomposed in terms of the contributions coming from the longitudinal (  ) and
transverse ($\perp$) polarizations of the virtual photon. Moreover, within a
collinear factorization scheme, they are given as convolutions of collinear
parton distribution functions (PDFs) with partonic helicity structure
functions (weighted by proper LDMEs). These, in turn, can be expressed as
functions of the partonic Mandelstam invariants.

The unpolarized cross section is obtained by integrating Eq. (6) over the
solid angle $\Omega$,

${1\over B_{ll}}{{\rm d}\sigma\over{\rm d}x_{\scriptscriptstyle B}\,{\rm
d}y\,{\rm d}z\,{\rm d}^{2}\bm{P}_{T}}={\alpha\over{8\,y\,z\,Q^{2}}}\,(2{\cal
W}_{T}+{\cal W}_{L})\,.$ (8)

It is then useful to introduce the ratio of polarized and unpolarized cross
sections

$\displaystyle\frac{{\rm d}N}{{\rm d}\Omega}\equiv\frac{{\rm d}\sigma}{{\rm
d}x_{\scriptscriptstyle B}\,{\rm d}y\,{\rm d}z\,{\rm d}^{2}\bm{P}_{T}\,{\rm
d}\Omega}\left(\frac{{\rm d}\sigma}{{\rm d}x_{\scriptscriptstyle B}\,{\rm
d}y\,{\rm d}z\,{\rm d}^{2}\bm{P}_{T}}\right)^{-1}\,,$ (9)

which can be expressed as follows

$\displaystyle\frac{{\rm d}N}{{\rm
d}\Omega}=\frac{3}{4\pi}\,\frac{1}{\lambda+3}\,\left[1+\lambda\cos^{2}\theta+\mu\sin
2\theta\cos\varphi+\frac{1}{2}\,\nu\,\sin^{2}\theta\cos 2\varphi\right]\,,$
(10)

where we have defined the polarization parameters

$\lambda={{\cal W}_{11}-{\cal W}_{00}\over{\cal W}_{11}+{\cal
W}_{00}},\quad\mu={\sqrt{2}\,{\rm Re}\left[{\cal W}_{10}\right]\over{\cal
W}_{11}+{\cal W}_{00}},\quad\nu={{\cal W}_{1,\,-1}\over{\cal W}_{11}+{\cal
W}_{00}}\,,$ (11)

or alternatively adopting Eq. (7),

$\lambda={{\cal W}_{T}-{\cal W}_{L}\over{\cal W}_{T}+{\cal
W}_{L}},\quad\mu={{\cal W}_{\Delta}\over{\cal W}_{T}+{\cal
W}_{L}},\quad\nu={2\,{\cal W}_{\Delta\Delta}\over{\cal W}_{T}+{\cal
W}_{L}}\,.$ (12)

The parameterizations shown in Eqs. (6) and (10) are standard for the study of
the angular distribution of a spin-one particle decay into a lepton pair and,
indeed, they are commonly adopted in Drell-Yan processes [26] and in $J/\psi$
photoproduction [27].

Among the polarization coefficients, $\lambda$, $\mu$ and $\nu$, the most
investigated experimentally is $\lambda$. Moreover, from the phenomenological
point of view it has a very intuitive interpretation, with $\lambda=+1(-1)$
describing a transverse(longitudinal) polarization state for the $J/\psi$
(i.e. a $J/\psi$ helicity equal to $\pm 1$ or 0), while $\lambda=0$ for an
unpolarized one.

The main goal of our study is to present estimates for these polarization
quantities, within both the CSM and the NRQCD frameworks, focusing on the
kinematic region accessible at the future EIC. As we will show in the
following, such a detailed phenomenological study could help in disentangling
among the production mechanisms.

LDME Set |  $\braket{{\cal O}_{1}[\,^{3}S_{1}]}$ $\big{[}{\rm GeV}^{3}\big{]}$ |  $\braket{{\cal O}_{8}[\,^{1}S_{0}]}$ $\big{[}{\rm GeV}^{3}\big{]}$ |  $\braket{{\cal O}_{8}[\,^{3}S_{1}]}$ $\big{[}{\rm GeV}^{3}\big{]}$ |  $\braket{{\cal O}_{8}[\,^{3}P_{0}]}$ $\big{[}{\rm GeV}^{5}\big{]}$
---|---|---|---|---
C12 | $1.16$ | $0.089$ | $0.003$ | $0.0126$
G13 | $1.16$ | $0.097$ | $-0.0046$ | $-0.0214$
BK11 | $1.32$ | $0.0304$ | $0.00168$ | $-0.00908$
Table 1: LDME set (central) values for the $J/\psi$ state: C12 [8], G13 [28]
and BK11 [29]. For the other ${}^{3}P_{J}$ states we use the standard spin-
symmetry relation $\braket{{\cal O}_{8}[\,^{3}P_{J}]}=(2J+1)\,\braket{{\cal
O}_{8}[\,^{3}P_{0}]}$.

##  III Angular distributions

In this section we analyze the polarization parameters defined in Eq. (11)
showing both their $z$ and $P_{T}$ distributions. The explicit analytic
expressions of the underlying partonic structure functions, calculated at the
perturbative order $\alpha_{s}^{2}$, are presented in Ref. [25] for the so-
called Gottfried-Jackson frame, together with all prescriptions needed to
transform them in the other relevant frames. For the predictions based on the
NRQCD approach, in addition to the CS contribution, given by a pure gluon
fusion channel, we consider the CO channels up to the order $v^{4}$, which
involve both gluon and quark final states. The CTEQ6L1 set [30] is used for
the unpolarized parton distribution functions. Moreover, in order to assess
the stability of our results against higher order corrections, we produce
uncertainty bands by varying the factorization scale $\mu_{F}$ in the range
$\mu_{0}/2<\mu_{F}<2\mu_{0}$, around the central value
$\mu_{0}=\sqrt{Q^{2}+M_{\psi}^{2}}$.

Concerning the CO LDME values, three different sets are adopted, see Table 1.
Here we only recall their main features: the C12 set [8] has been extracted
simultaneously from both polarized and unpolarized $J/\psi$ production data in
$pp$ collision at $P_{T}>7$ GeV, measured by the CDF (Run II) Collaboration;
the G13 set [28] is obtained including only $P_{T}>7$ GeV unpolarized data
from the CDF and LHCb Collaborations and then used to predict the $J/\psi$
polarization in $pp$ collisions; it is in agreement with the C12 set if feed-
down contribution is negligible; the BK11 set [29] is based on a fit without
any polarization data, but starting from a lower $P_{T}$ value, around $3$
GeV, and including both photoproduction and hadroproduction data.

The high cm energy kinematical set-ups expected at the EIC are an ideal
environment to study $J/\psi$ polarization in electroproduction. Moreover,
they will allow to better explore high photon virtualities ($Q$), avoiding the
competing contributions from photoproduction. Furthermore, since we are
interested in the region where collinear factorization holds, our results will
be shown only for $P_{T}$ values above $P_{T\rm{min}}=1$ GeV. Notice that
around this value we actually enter the region where the transverse momentum
dependent (TMD) factorization could be applied and therefore our estimates are
pushed down to the overlapping region of validity of the two factorization
schemes.

###  III.1 The $\lambda$ parameter

In Fig. 1 we present our predictions for $\lambda$ at $\sqrt{s}=140$ GeV, as a
function of both the $J/\psi$ energy fraction $z$ (left panels) and its
transverse momentum $P_{T}$ (right panels). Two quarkonium rest frames are
explicitly considered: the Gottfried-Jackson (upper panels) and the Helicity
(lower panels) ones. In this and in the following figures, the kinematical
ranges explored are indicated in the legend boxes. For completeness we report
here the corresponding regions explored in $x_{B}$ and $y$ at $\sqrt{s}=140$
GeV, $10^{-3}\lesssim x_{B}\lesssim 0.2$ and $y\lesssim 0.5$ respectively,
even if the effectively probed maximum value in $x_{B}$ is around 0.07.

Concerning other typical frames, like the Target and Collins-Soper ones, we
only notice that the first one give estimates very close to those in the
Helicity frame, while predictions obtained in the second one, at least for the
kinematics considered, are in general much smaller than those in the
Gottfried-Jackson frame or even close to zero.

Notice that for such observable, defined as a ratio of cross sections, the
dependence on the scale $\mu_{F}$ in the range $[\mu_{0}/2,2\mu_{0}]$ is
barely appreciable and therefore is not shown.

Figure 1: Estimates for $\lambda$ at $\sqrt{s}=140$ GeV as a function of $z$
(left panels) and $P_{T}$ (right panels) for different models and LDME sets
and two reference frames: Gottfried-Jackson (upper panels) and Helicity (lower
panels) frames. Integration ranges are given in the light-blue legend box.
Figure 2:  Estimates for $\lambda$ at cm energy $\sqrt{s}=45$ GeV. The
integration region, different with respect to the higher-energy case, is given
in the red legend box, while curves and panels have the same meaning as in
Fig. 1. The scale error bands are sizable and explicitly shown only for the
CSM prediction as a function of $P_{T}$.

The study of the $\lambda$ parameter as a function of $z$ presents very
interesting features from the phenomenological point of view. The reasons are
manifold: first of all its expected relative large size as compared to the
$\mu$ and $\nu$ parameters. Moreover, it is experimentally under more active
investigation. On the other hand, theoretical estimates for $\lambda$ as a
function of $z$ (for small and moderate values) do not vary significantly
adopting different frameworks (Fig. 1, left panels), which implies that, in
order to get information on the quarkonium formation mechanism, one would need
highly precise measurements. The same problem was found in different analyses
performed by the HERA Collaborations, Refs. [21, 23].

The situation changes considerably at $z>0.6$, which represents a very
interesting region from the phenomenological point of view. As is well known,
NRQCD estimates for the unpolarized cross section manifest a divergent
behavior as $z\to 1$, due to the corresponding $\hat{t}\to 0$ singularities.
This can potentially spoil the validity of NRQCD factorization. As shown in
Ref. [31], in order to extend the region of applicability of NRQCD up to
$1-z\sim v^{2}$, one can introduce a new set of functions, the so-called shape
functions [32], that allow to improve noticeably the convergence for
photoproduction. We expect such quantities to be relevant also for the SIDIS
process, together with their TMD extensions, which have been adopted in the
study of $pp$ collisions in Refs. [33, 34] and whose perturbative tails have
been derived in Ref. [35] for unpolarized and in Ref. [25] for polarized
$J/\psi$ SIDIS. On the other hand, the impact of the shape functions on
$\lambda$ is expected to be strongly reduced since $\lambda$ is a ratio of
cross sections. This can be tested with future available data.

A much more powerful tool to assess the relevance of the CO contributions is
the study of the $P_{T}$ distribution (Fig. 1, right panels). In the
Gottfried-Jackson frame (upper panel) we see a clear separation as well as a
different behavior between the CSM and NRQCD curves, in particular in the
region $4<P_{T}<7$ GeV; similarly in the Helicity frame there is a wide
separation between the CSM and the NRQCD curves, while different LDME sets
give predictions much closer to each other and closer to $\lambda=0$. It is
worth noticing that, even if the unpolarized cross section decreases as
$P_{T}$ increases, a good separation can be found already around $P_{T}\simeq
5$ GeV, which is also far away from the TMD region.

Before concluding the analysis of $\lambda$ at large cm energies, a comment on
the contributions from different partonic channels and/or different NRQCD
waves can be useful. Concerning the $z$ distribution, we find that the main
contribution to the numerator of $\lambda$ comes from the (gluon) CS wave,
while the differences among NRQCD predictions, especially around $z\to 0.9$,
are due to the gluon P-wave, modulated by the corresponding LDME parameter.
For the $P_{T}$ distribution we find, similarly, that the CS term is on the
whole the most relevant contribution, followed again by the gluon P-wave one.
In particular at $P_{T}\to 1$ GeV the size of the gluon P-wave contribution
becomes comparable to (or even larger than) the CS one; moreover, since the
low-$P_{T}$ region dominates the integration over $P_{T}$, one can also
understand why the gluon P-wave is so relevant in our estimates vs. $z$, with
the most visible effects for $z\to 0.9$.

At medium $P_{T}$ values the quark P-wave starts becoming important and at
even higher $P_{T}$ values it is similar in size to the gluon one; this means
that in this region, the full P-wave contribution (gluon+quark) dominates over
the CS one.

Another interesting possibility given by the future EIC facility is the
corresponding analysis at smaller energies: in the following we will adopt
$\sqrt{s}=45$ GeV. In this case, different integration ranges have been
considered for $W$ and $Q^{2}$, as reported in the legend box of Fig. 2.
These, in turn, correspond to $10^{-3}\lesssim x_{B}\lesssim 0.5$ (with an
effective upper limit around $x_{B}\simeq 0.2$) and $y\lesssim 0.8$, a more
valence-like region w.r.t. the previous case. Moreover, since at lower
energies it is more difficult to reach high photon virtualities, we get
contributions mostly from moderately low $Q^{2}$. Consistently we adopt a
lower limit, $Q_{\rm min}\simeq 1.6$ GeV, in the integration. Notice that in
this kinematic region, at least for the high $P_{T}$ dependence of $\lambda$
within the CSM, the scale error bands are once again sizeable enough.

From Fig. 2 (left panels) we can see that the $z$ distribution does not depend
significantly on the energy for $z\leq 0.6$, while at higher $z$ values the
estimates are closer to zero, at variance with those at higher cm energy. As
said, a polarization study pushed up to this regime can suffer from
factorization breaking effects in NRQCD even if data in this region could be
relevant from the phenomenological point of view. We also observe a rapid
variation of all curves in the Helicity frame at $z\sim 0.1$. This is due to
geometrical factors which are energy dependent (see also Eq. (A16) of Ref.
[25]). The same variation is also present at higher cm energy, but for $z<0.1$
(outside the range shown in the lower-left panel of Fig. 1).

Concerning the $P_{T}$ dependence, Fig. 2 (right panels), we notice that the
CSM results are very different with respect to the corresponding ones in Fig.
1, while the same is not true for the NRQCD cases. This is related to the
different virtualities explored, on which the CSM estimates depend heavily.
This difference can be considered as an extra tool in the quest of discerning
among different frameworks.

Finally, we briefly comment on how the parton and/or wave contributions vary
with the energy. While the $z$ distribution manifests almost no energy
dependence, the $P_{T}$ spectrum presents interesting features in the two
frames considered. For the Gottfried-Jackson one the relative contribution
from the quark P-wave is widely increased at this lower energy, making it the
leading term in the numerator at medium/high $P_{T}$. Regarding the Helicity
frame the situation is, potentially, even more interesting, since the CSM and
P-wave (both gluon and quark) contributions are highly suppressed at this
energy, especially at large $P_{T}$. The main role is then played by the
${}^{3}\\!S_{1}^{(8)}$ quark wave, which is responsible for the difference
among the predictions based on the LDME sets considered. Even if in this
region it is quite hard to expect precise enough data to discriminate between
models, it is nevertheless worth stressing that it could be very useful in
constraining the nonperturbative physics.

###  III.2 The $\mu$ parameter

Estimates for the $\mu$ parameter are again provided both in the Gottfried-
Jackson and in the Helicity frames, as a function of $z$ and $P_{T}$ at
$\sqrt{s}=140$ GeV, Fig. 3, and $\sqrt{s}=45$ GeV, Fig. 4.

Figure 3: Estimates for the parameter $\mu$ at $\sqrt{s}=140$ GeV. Paneling
order is the same as in Fig. 1. Integration ranges are given in the blue
legend box. Figure 4: Estimates for the parameter $\mu$ at $\sqrt{s}=45$ GeV.
Paneling order is the same as in Fig. 1. Integration ranges are given in the
red legend box.

From these figures we see that the Gottfried-Jackson frame is the best choice
to discern among the CSM and NRQCD approach. A similar conclusion holds for
the parameter $\nu$ as well, see the discussion in Sec. III.3. Indeed, in Fig.
3 the separation between the CSM estimates and the corresponding NRQCD ones
are remarkably sizeable for $z\gtrsim 0.5$ and $P_{T}\gtrsim 5$ GeV. On the
contrary, estimates in the Helicity frame both with respect to $z$ and $P_{T}$
are so close to each other that one cannot draw any conclusion.

The wave/parton decomposition of the ${\cal W}_{\Delta}$ helicity function,
that is directly related to the $\mu$ numerator, allows us to get some further
insights. The main CO contribution comes from the $P$-wave term. In
particular, differences in NRQCD predictions as a function of $z$ (left panels
of Fig. 3) are driven by the gluon $P$-wave LDMEs. Moreover, the gluon
$P$-wave dominates the numerator behavior with respect to $P_{T}$ too (right
panels of Fig. 3). In addition, we find that the NRQCD predictions in the
Gottfried-Jackson frame receive a significant contribution from the gluon
$P$-wave also at low-$P_{T}$, namely $P_{T}\lesssim 3$ GeV. At variance with
the behavior in $z$, here the quark $P$-wave channel is relevant at high
$P_{T}$, especially when considering the Helicity frame.

Moving to the lower cm energy, we see that the CSM $\mu$ estimates in the
Gottfried-Jackson frame, Fig. 4 (upper panels), vary significantly for
$z\gtrsim 0.5$ and $P_{T}\gtrsim 5$ GeV, as compared with what happens at
$\sqrt{s}=140$ GeV. We remark that this variation can also appear via a proper
$Q$-binning in the higher cm energy case ($\sqrt{s}=140$ GeV). In contrast,
estimates within the Helicity frame at lower energies (lower panels of Fig. 4)
do not present the same energy/$Q$-binning dependence. The only remarkable
exception resides in the $P_{T}$ distribution, where CSM predictions increase
up to $\sim 40\%$, to be compared with the $\sqrt{s}=140$ GeV case where the
CSM result is at most $\sim 25\%$. Despite this, $\mu$ estimates in the
Helicity frame do not differ enough to discern among different models.

Looking at the wave/parton decomposition, we confirm that also for the $\mu$
numerator the role of quarks is enhanced at lower energies. This is
particularly true for the $P_{T}$ dependence. Here we find that NRQCD
predictions at the higher $P_{T}$ values, namely $P_{T}\gtrsim 6$ GeV, are
mostly driven by the quark $P$-wave; moreover, in the same $P_{T}$ region we
observe that the ${}^{3}S_{1}^{[8]}$ quark wave is non-negligible.

###  III.3 The $\nu$ parameter

We now discuss the parameter $\nu$, which is particularly important in the TMD
framework, since it is directly related to the TMD distribution of linearly
polarized gluons inside an unpolarized proton, $h_{1}^{\perp g}$. This could
play a role in the region of moderately low $P_{T}$, where the two
factorization schemes overlap.

Again, we focus initially on the higher cm energy ($\sqrt{s}=140$ GeV), Fig.
5, and then we describe the main differences with respect to the smaller cm
energy ($\sqrt{s}=45$ GeV), Fig. 6.

Figure 5: Estimates for the parameter $\nu$ at $\sqrt{s}=140$ GeV. Paneling
order is the same as in Fig. 1. Integration ranges are given in the blue
legend box.

Starting from the $z$-dependent distribution in Fig. 5 (left panels), we see
once again that even if the estimated $\nu$ values are potentially sizeable,
at least in the Helicity frame, the separation among the different approaches
is in general very poor. Nevertheless, it is worth remarking that at high $z$
we find more sensitivity to the LDME sets in the NRQCD framework. The
situation is slightly different for the $P_{T}$ case (right panels): if the
Helicity frame does not show a promising scenario, in the Gottfried-Jackson
case the differences in the medium/high-$P_{T}$ region between the two
approaches are sizeable.

As said, results at high $z$ and/or small $P_{T}$ are in general promising for
future analyses regarding the $h_{1}^{\perp g}$ gluon distribution in the TMD
region. Nevertheless, it is important to remark that for the $\nu$ parameter
the shape functions and their TMD extensions enter, potentially, in a
different way in the numerator and the denominator, and their role could be
important. This requires further investigation, together with a full higher-
order description in $\alpha_{s}$, which is not available at present.

It is once again interesting to look into the parton and wave decomposition.
The $z$-dependent ${\cal W}_{\Delta\Delta}$ is dominated, for almost all $z$
values, by the CS wave; only for $z\to 0.9$ the CS contribution becomes
negligible, and the results are driven by the CO P-wave, in particular by the
gluon term. Moving to the $P_{T}$ dependence, we find again some similarities
with the $\lambda$ case: the CS term is the relevant contribution to the
numerator over the whole $P_{T}$ spectrum, together with the gluon P-wave. At
variance with the $\lambda$ parameter case, the quark contribution to the
P-wave term starts becoming important already at small-$P_{T}$ values.

Figure 6: Estimates for the parameter $\nu$ at $\sqrt{s}=45$ GeV. Paneling
order is the same as in Fig. 1. Integration ranges are given in the red legend
box.

Moving to the lower cm energy, from Fig. 6 we see that the $z$ distribution is
sensitive to the energy change in the whole spectrum, at variance with the
$\lambda$ case. The differences, particularly noticeable in the Gottfried-
Jackson frame, are mostly in size and not in the general behavior, implying
that even in this case it would be difficult to extract any information.
Again, we remark that the rapid variation of $\nu$ estimates at low-$z$ values
is due to a geometrical factor (Eq. (A16) of Ref. [25]). The $P_{T}$-dependent
distributions, instead, have a quite different behavior for the two frames
displayed. The Gottfried-Jackson estimates vary significantly in size,
especially if one considers the CSM; moreover all the LDME sets give similar
predictions, compatible with zero, for $P_{T}>5$ GeV, while predictions, in
both approaches, are sizeable (up to $\sim 20\%$) at low-$P_{T}$ values. This
could be very promising for further extensions to the TMD region. The curves
in the Helicity frame, instead, do not show the same dependence on the energy.
In general, we conclude that the study of the $\nu$ parameter, at least in
this frame, is not very effective. Nevertheless it becomes more interesting
when its information is combined with other parameters, as done in the study
of the invariant quantities in the next section, Sec. IV.

Concerning the wave decomposition, we find that both quark and gluon P-wave
contributions to the $P_{T}$ and $z$ distributions are enhanced at lower
energies, even if for the latter this is true only at large $z$ values. Notice
that the different (larger) size of the $\nu$ parameter at $z\to 0.9$ could
also affect the TMD region, increasing the possibility of extracting
information on the linearly polarized gluon distribution. The main source of
this enhancement at $\sqrt{s}=45$ GeV is related once again to the lower
photon virtualities explored. In this sense, very similar predictions might be
expected at higher cm energy via a binned analysis with
$1.6\leavevmode\nobreak\ {\rm GeV}<Q<M_{\psi}$.

##  IV Rotational invariants

The polarization parameters $\lambda$, $\mu$ and $\nu$, as widely discussed in
the previous sections, are frame dependent by definition, since they are
expressed with respect to the solid angle $\Omega$ spanned by the $l^{+}$
particle in the $J/\psi$ decay and in its rest frame. As already pointed out,
the frame choice is not unique and the results appear different from frame to
frame. On the other hand, the relations among the most used reference frames
are computable, since they differ only in the $Z$-axis direction.

A complementary and powerful tool to study $J/\psi$ polarization, both from
the experimental and the phenomenological points of view, is the use of
rotational invariant parameters, that are rest-frame independent by
construction. These can be defined taking into account what follows.

For all the most common choices, the $Z$\- and $X$-axes, lying in the $J/\psi$
production plane, are defined in terms of physical momenta in the quarkonium
rest frame (see Appendix A of Ref. [25]), with the $Y$-axis always
perpendicular with respect to this plane and always pointing in the same
direction. This implies that two frames ($F,F^{\prime}$) can be connected by a
simple rotation of an angle $\psi$ around the $Y$-axis, and the corresponding
polarization parameters can be directly related as111Here $\mu_{F}$ stands for
the $\mu$ parameter in a specific frame $F$, not to be confused with the
factorization scale $\mu_{F}$ defined in the previous sections.

$\begin{pmatrix}\lambda\\\ \mu\\\ \nu\end{pmatrix}_{F^{\prime}}={1\over
1+\rho}\begin{pmatrix}1-{3\over 2}\sin^{2}\psi&{3\over 2}\sin 2\psi&{3\over
4}\sin^{2}\psi\\\ -{1\over 2}\sin 2\psi&\cos 2\psi&{1\over 4}\sin 2\psi\\\
\sin^{2}\psi&-\sin 2\psi&1-{1\over
2}\sin^{2}\psi\end{pmatrix}\begin{pmatrix}\lambda\\\ \mu\\\
\nu\end{pmatrix}_{F}\,,$ (13)

with

$\rho={\sin^{2}\psi\over 2}\left(\lambda_{F}-{\nu_{F}\over 2}\right)-\sin
2\psi\,{\mu_{F}\over 2},$ (14)

as given in Eqs. (A.18) and (A.19) of Ref. [25], where we have changed the
rotation angle from $\theta$ to $\psi$ to avoid any confusion with the polar
angle of the final lepton $l^{+}$. Notice that the quantity $\rho$ depends on
the kinematics, since the rotation angle itself depends on the partonic
Mandelstam variables (see Eqs. (A.14)-(A.16) of Ref. [25] for details).

From Eq. (13), one can construct several quantities which do not change upon
rotation around the $Y$ direction. The following relations are extremely
useful in this respect:

$3+\lambda_{F^{\prime}}={1\over 1+\rho}\left(3+\lambda_{F}\right),\qquad
1-{\nu_{F^{\prime}}\over 2}={1\over 1+\rho}\left(1-{\nu_{F}\over 2}\right)\,.$
(15)

A group of rotational invariants, as initially proposed in Ref. [36], can be
defined in terms of two polarization parameters, namely $\lambda$ and $\nu$,

${\cal F}_{(c_{i})}={c_{0}(3+\lambda)+c_{1}(1-\nu/2)\over
c_{2}(3+\lambda)+c_{3}(1-\nu/2)}\,,$ (16)

where $c_{i}$ are suitable free constants.

Among all possible combinations, two of them play an important role and have
received special attention [37, 38, 39, 40, 41]

${\cal F}\equiv{\cal F}_{(1,-2,1,0)}={1+\lambda+\nu\over 3+\lambda}$ (17)

and

$\tilde{\lambda}\equiv{\cal F}_{(1,-3,0,1)}={2\,\lambda+3\,\nu\over 2-\nu}\,.$
(18)

These invariants have been widely studied in $pp$ and heavy-ion processes [42,
43].

It is worth noticing that both invariants can be similarly defined for Drell-
Yan processes, where they acquire a constant value if the Lam-Tung relation
($1-\lambda=2\nu$) holds [26]: ${\cal F}_{\rm DY}=1/2$ and
$\tilde{\lambda}_{\rm DY}=+1$, as pointed out in Refs. [38, 41]. Another
interesting feature is that $\tilde{\lambda}=+1(-1)$ is related to a natural
transverse (longitudinal) polarization [36]. It is important to stress that
the constant behavior is purely dynamical, and in particular for the Drell-Yan
case is a consequence of rotational invariance and helicity conservation [44].
Since $J/\psi$ couples differently in SIDIS processes, the Lam-Tung relation
is expected to be broken in this case.

Figure 7: Estimates for the invariant ${\cal F}$, Eq. (17), as a function of
$z$ (left panels) and $P_{T}$ (right panels) at two cm energies,
$\sqrt{s}=140$ GeV (upper panels) and $\sqrt{s}=45$ GeV (lower panels), for
different approaches and LDME sets. Kinematic ranges are given in the legend
boxes.

Not all the invariants belong to the previous family. Indeed, one can exploit
another relation that involves all polarization parameters in two frames and
that, upon rotation around the $Y$-axis, reads

$(\lambda_{F^{\prime}}-\nu_{F^{\prime}}/2)^{2}+4\mu_{F^{\prime}}^{2}={(\lambda_{F}-\nu_{F}/2)^{2}+4\mu_{F}^{2}\over(1+\rho)^{2}}\,.$
(19)

From this, one can construct an invariant quantity involving the polarization
parameters squared, as first pointed out in Ref. [45]. As an example, we
recall

$\tilde{\lambda}^{\prime}={(\lambda-\nu/2)^{2}+4\mu^{2}\over(3+\lambda)^{2}}\,,$
(20)

as introduced in Ref. [41].

The study of rotational invariants has not only a theoretical interest, but it
is relevant also from the experimental point of view, since their expected
equality among different frames is an important check of experimental
acceptances and systematics as shown, for instance, by the ATLAS Collaboration
[46].

For these reasons, we consider, as a case of study, one of these quantities at
the kinematics explored by the EIC. In Fig. 7 we show the theoretical
estimates in the collinear framework, for the invariant ${\cal F}$, Eq. (17),
as a function of $z$ (left panels) and $P_{T}$ (right panels). Once again we
compute this quantity at two energies, $\sqrt{s}=140$ GeV (upper panels) and
$\sqrt{s}=45$ GeV (lower panels) for different approaches and LDME sets.

From Fig. 7 we clearly see that ${\cal F}$ is not equal to $1/2$, as expected
from the Lam-Tung relation. Moreover, it is neither a constant, since its
value depends on both $z$ and $P_{T}$ variables. In principle, for some LDME
sets a constant behavior could accidentally appear, but this would be limited
to a specific kinematic region.

Another interesting remark is that, while the denominator of ${\cal F}$ is
proportional to the unpolarized cross section, its numerator is controlled by
the relative size of the $\lambda$ and $\nu$ parameters. This can vary
significantly, depending on the frames and approaches adopted, as discussed in
the previous Section.

From this preliminary study we can conclude that, even if not easily
accessible from the experimental point of view, these invariant quantities
could represent an invaluable tool to learn on the $J/\psi$ polarization
mechanism.

##  V Conclusions

The study of quarkonium polarization, interesting by itself, is also a
powerful tool to explore the still challenging issue of its formation
mechanism within QCD. In this spirit, we have presented a phenomenological
analysis of $J/\psi$ polarization in SIDIS at large $P_{T}$. More
specifically, we have looked at the dilepton angular distribution in the
$J/\psi\to\ell^{+}\ell^{-}$ decay in terms of the associated polarization
parameters, that could be accessed at the future EIC. By exploiting the
theoretical results of Ref. [25], we have computed the parameters, $\lambda$,
$\mu$ and $\nu$, in different frames, trying to emphasize whether one can use
these observables to discriminate among two well consolidated frameworks,
still under investigation: the Color Singlet Model and the NRQCD approach.
Moreover, for the latter we have employed three different LDME sets, based on
different extractions and assumptions, highlighting their impact on quarkonium
polarization estimates.

We have shown results both as a function of $z$ and $P_{T}$, adopting two
quite different cm energies, for standard kinematics at the EIC, together with
a detailed analysis in terms of parton and NRQCD wave contributions.

The main findings of our study can be summarized as follows: $i)$ concerning
the $\lambda$ parameter, the large-$z$ region, both in the Gottfried-Jackson
and the Helicity frame, turns out to be very promising, with the only caveat
of possible contributions from (TMD) shape functions (even if expected to be
reduced being $\lambda$ a ratio of helicity structure functions); similarly
its $P_{T}$ distribution, at medium-large values, could be an ideal ground to
disentangle the formation mechanisms, both at high and low energies. $ii)$ The
$\mu$ parameter displays some interesting features when studied in the
Gottfried-Jackson frame, namely: a clear separation among the estimates in
different frameworks at medium-large $z$ or as a function of $P_{T}$ in the
high-energy set-up; a different behavior with respect to the corresponding
lower-energy estimates at medium-large $z$ or at moderate $P_{T}$. Moreover,
in the Helicity frame at low energies one could extract important information
by looking in the large $P_{T}$ region. $iii)$ Similarly, for the $\nu$
parameter, relevant also in the context of the TMD framework, medium-large
$P_{T}$ values in the Gottfried-Jackson frame are certainly worth to be
explored.

Finally, we have discussed a selection of frame-independent (rotational
invariant) polarization parameters, relevant not only from the theory point of
view, but extremely useful as an important check of experimental acceptances
and systematics. In particular, we have focused on the invariant ${\cal{F}}$,
controlled by the relative weight of the $\lambda$ and $\nu$ parameters, that
strongly depend on the frames and frameworks adopted. As shown, this
observable could clearly help in getting information on the $J/\psi$ formation
mechanism, both at large $z$ (high- and low-energy set-ups) and as a function
of $P_{T}$ (at large energy).

We can certainly conclude that a study of the dilepton angular distribution in
$J/\psi$ decay in SIDIS at the EIC could be an invaluable tool to shed light
on the $J/\psi$ polarization as well as on its formation mechanism.

###### Acknowledgements.

We thank P. Faccioli, T. Stebel and R. Venugopalan for clarifying some aspects
concerning the rotational invariants. This project has received funding from
the European Union’s Horizon 2020 research and innovation programme under
grant agreement STRONG 2020—No 824093. U.D. and C.P. also acknowledge
financial support by Fondazione di Sardegna under the project “Proton
tomography at the LHC”, project number F72F20000220007 (University of
Cagliari).

## References

  * [1] E598 Collaboration, J. J. Aubert et al., “Experimental Observation of a Heavy Particle $J$,” Phys. Rev. Lett. 33 (1974) 1404–1406.
  * [2] SLAC-SP-017 Collaboration, J. E. Augustin et al., “Discovery of a Narrow Resonance in $e^{+}e^{-}$ Annihilation,” Phys. Rev. Lett. 33 (1974) 1406–1408.
  * [3] R. Baier and R. Ruckl, “Hadronic collisions: A quarkonium factory,” Z. Phys. C 19 (1983) 251.
  * [4] G. T. Bodwin, E. Braaten, and G. P. Lepage, “Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium,” Phys. Rev. D 51 (1995) 1125–1171, arXiv:hep-ph/9407339. [Erratum: Phys.Rev.D 55, 5853 (1997)].
  * [5] G. C. Nayak, J.-W. Qiu, and G. F. Sterman, “NRQCD factorization and the velocity dependence of NNLO poles in heavy quarkonium production,” Phys. Rev. D 74 (2006) 074007, arXiv:hep-ph/0608066.
  * [6] G. P. Lepage, L. Magnea, C. Nakhleh, U. Magnea, and K. Hornbostel, “Improved nonrelativistic QCD for heavy-quark physics,” Phys. Rev. D 46 (1992) 4052–4067, arXiv:hep-lat/9205007.
  * [7] M. Butenschoen and B. A. Kniehl, “Reconciling $J/\psi$ production at HERA, RHIC, Tevatron, and LHC with NRQCD factorization at next-to-leading order,” Phys. Rev. Lett. 106 (2011) 022003, arXiv:1009.5662 [hep-ph].
  * [8] K.-T. Chao, Y.-Q. Ma, H.-S. Shao, K. Wang, and Y.-J. Zhang, “$J/\psi$ Polarization at Hadron Colliders in Nonrelativistic QCD,” Phys. Rev. Lett. 108 (2012) 242004, arXiv:1201.2675 [hep-ph].
  * [9] R. Sharma and I. Vitev, “High transverse momentum quarkonium production and dissociation in heavy ion collisions,” Phys. Rev. C 87 no. 4, (2013) 044905, arXiv:1203.0329 [hep-ph].
  * [10] G. T. Bodwin, H. S. Chung, U.-R. Kim, and J. Lee, “Fragmentation contributions to $J/\psi$ production at the Tevatron and the LHC,” Phys. Rev. Lett. 113 no. 2, (2014) 022001, arXiv:1403.3612 [hep-ph].
  * [11] H.-F. Zhang, Z. Sun, W.-L. Sang, and R. Li, “Impact of $\eta_{c}$ hadroproduction data on charmonium production and polarization within the NRQCD framework,” Phys. Rev. Lett. 114 no. 9, (2015) 092006, arXiv:1412.0508 [hep-ph].
  * [12] N. Brambilla et al., “Heavy Quarkonium: Progress, Puzzles, and Opportunities,” Eur. Phys. J. C 71 (2011) 1534, arXiv:1010.5827 [hep-ph].
  * [13] A. Andronic et al., “Heavy-flavour and quarkonium production in the LHC era: from proton–proton to heavy-ion collisions,” Eur. Phys. J. C 76 no. 3, (2016) 107, arXiv:1506.03981 [nucl-ex].
  * [14] J.-P. Lansberg, “New observables in inclusive production of quarkonia,” Phys. Rept. 889 (2020) 1–106, arXiv:1903.09185 [hep-ph].
  * [15] H.-F. Zhang, W.-L. Sang, and Y.-P. Yan, “Statistical analysis of the azimuthal asymmetry in the $J/\psi$ leptoproduction in unpolarized $ep$ collisions,” JHEP 10 (2019) 234, arXiv:1908.02521 [hep-ph].
  * [16] J.-W. Qiu, X.-P. Wang, and H. Xing, “Exploring $J/\psi$ Production Mechanism at the Future Electron-Ion Collider,” Chin. Phys. Lett. 38 no. 4, (2021) 041201, arXiv:2005.10832 [hep-ph].
  * [17] D. Boer, C. Pisano, and P. Taels, “Extracting color octet NRQCD matrix elements from $J/\psi$ production at the EIC,” Phys. Rev. D 103 no. 7, (2021) 074012, arXiv:2102.00003 [hep-ph].
  * [18] R. Abdul Khalek et al., “Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report,” Nucl. Phys. A 1026 (2022) 122447, arXiv:2103.05419 [physics.ins-det].
  * [19] A. Accardi et al., “Electron Ion Collider: The next QCD frontier - Understanding the glue that binds us all,” Eur. Phys. J. A 52 no. 9, (2016) 268, arXiv:1212.1701 [nucl-ex].
  * [20] D. Boer et al., “Gluons and the quark sea at high energies: Distributions, polarization, tomography,” arXiv:1108.1713 [nucl-th].
  * [21] H1 Collaboration, C. Adloff et al., “Inelastic leptoproduction of $J/\psi$ mesons at HERA,” Eur. Phys. J. C 25 (2002) 41–53, arXiv:hep-ex/0205065.
  * [22] S. Fleming and T. Mehen, “Leptoproduction of $J/\psi$,” Phys. Rev. D 57 (1998) 1846–1857, arXiv:hep-ph/9707365.
  * [23] F. Yuan and K.-T. Chao, “Polarized $J/\psi$ production in deep inelastic scattering at DESY HERA,” Phys. Rev. D 63 (2001) 034017, arXiv:hep-ph/0008301. [Erratum: Phys.Rev.D 66, 079902 (2002)].
  * [24] Z. Sun and H.-F. Zhang, “QCD corrections to the color-singlet $J/\psi$ production in deeply inelastic scattering at HERA,” Phys. Rev. D 96 no. 9, (2017) 091502, arXiv:1705.05337 [hep-ph].
  * [25] U. D’Alesio, L. Maxia, F. Murgia, C. Pisano, and S. Rajesh, “$J/\psi$ polarization in semi-inclusive DIS at low and high transverse momentum,” JHEP 03 (2022) 037, arXiv:2110.07529 [hep-ph].
  * [26] C. S. Lam and W.-K. Tung, “Systematic approach to inclusive lepton pair production in hadronic collisions,” Phys. Rev. D 18 (1978) 2447.
  * [27] M. Beneke, M. Kramer, and M. Vanttinen, “Inelastic photoproduction of polarized $J/\psi$,” Phys. Rev. D 57 (1998) 4258–4274, arXiv:hep-ph/9709376.
  * [28] B. Gong, L.-P. Wan, J.-X. Wang, and H.-F. Zhang, “Polarization for Prompt $J/\psi$ and $\psi(2s)$ Production at the Tevatron and LHC,” Phys. Rev. Lett. 110 no. 4, (2013) 042002, arXiv:1205.6682 [hep-ph].
  * [29] M. Butenschoen and B. A. Kniehl, “World data of $J/\psi$ production consolidate NRQCD factorization at next-to-leading order,” Phys. Rev. D 84 (2011) 051501, arXiv:1105.0820 [hep-ph].
  * [30] J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky, and W. K. Tung, “New generation of parton distributions with uncertainties from global QCD analysis,” JHEP 07 (2002) 012, arXiv:hep-ph/0201195.
  * [31] M. Beneke, G. A. Schuler, and S. Wolf, “Quarkonium momentum distributions in photoproduction and $B$ decay,” Phys. Rev. D 62 (2000) 034004, arXiv:hep-ph/0001062.
  * [32] M. Beneke, I. Z. Rothstein, and M. B. Wise, “Kinematic enhancement of non-perturbative corrections to quarkonium production,” Phys. Lett. B 408 (1997) 373–380, arXiv:hep-ph/9705286.
  * [33] M. G. Echevarria, “Proper TMD factorization for quarkonia production: $pp\to\eta_{c,b}$ as a study case,” JHEP 10 (2019) 144, arXiv:1907.06494 [hep-ph].
  * [34] S. Fleming, Y. Makris, and T. Mehen, “An effective field theory approach to quarkonium at small transverse momentum,” JHEP 04 (2020) 122, arXiv:1910.03586 [hep-ph].
  * [35] D. Boer, U. D’Alesio, F. Murgia, C. Pisano, and P. Taels, “$J/\psi$ meson production in SIDIS: matching high and low transverse momentum,” JHEP 09 (2020) 040, arXiv:2004.06740 [hep-ph].
  * [36] P. Faccioli, C. Lourenco, and J. Seixas, “New approach to quarkonium polarization studies,” Phys. Rev. D 81 (2010) 111502, arXiv:1005.2855 [hep-ph].
  * [37] P. Faccioli, C. Lourenco, and J. Seixas, “Rotation-invariant relations in vector meson decays into fermion pairs,” Phys. Rev. Lett. 105 (2010) 061601, arXiv:1005.2601 [hep-ph].
  * [38] P. Faccioli, C. Lourenco, J. Seixas, and H. K. Wohri, “Towards the experimental clarification of quarkonium polarization,” Eur. Phys. J. C 69 (2010) 657–673, arXiv:1006.2738 [hep-ph].
  * [39] P. Faccioli, C. Lourenco, J. Seixas, and H. K. Wohri, “Quarkonium polarization in $pp$ and $p$-nucleus collisions,” Nucl. Phys. A 855 (2011) 116–124.
  * [40] P. Faccioli, C. Lourenco, J. Seixas, and H. K. Wohri, “Quarkonium polarization measurements,” Nucl. Phys. B Proc. Suppl. 214 (2011) 97–102.
  * [41] J.-C. Peng, D. Boer, W.-C. Chang, R. E. McClellan, and O. Teryaev, “On the rotational invariance and non-invariance of lepton angular distributions in Drell–Yan and quarkonium production,” Phys. Lett. B 789 (2019) 356–359, arXiv:1808.04398 [hep-ph].
  * [42] Y.-Q. Ma, T. Stebel, and R. Venugopalan, “$J/\psi$ polarization in the CGC+NRQCD approach,” JHEP 12 (2018) 057, arXiv:1809.03573 [hep-ph].
  * [43] ALICE Collaboration, S. Acharya et al., “Measurement of the inclusive $J/\psi$ polarization at forward rapidity in $pp$ collisions at $\sqrt{s}=8$ TeV,” Eur. Phys. J. C 78 no. 7, (2018) 562, arXiv:1805.04374 [hep-ex].
  * [44] P. Faccioli, C. Lourenco, J. Seixas, and H. K. Wohri, “Model-independent constraints on the shape parameters of dilepton angular distributions,” Phys. Rev. D 83 (2011) 056008, arXiv:1102.3946 [hep-ph].
  * [45] S. Palestini, “Angular distribution and rotations of frame in vector meson decays into lepton pairs,” Phys. Rev. D 83 (2011) 031503, arXiv:1012.2485 [hep-ph].
  * [46] ATLAS Collaboration, G. Aad et al., “Measurement of the differential cross-sections of inclusive, prompt and non-prompt $J/\psi$ production in proton-proton collisions at $\sqrt{s}=7$ TeV,” Nucl. Phys. B 850 (2011) 387–444, arXiv:1104.3038 [hep-ex].