Physics case for quarkonium studies at the Electron Ion Collider
Daniël Boer Chris A. Flett Carlo Flore Daniel Kikoła Jean-Philippe Lansberg Maxim Nefedov Charlotte Van Hulse Shohini Bhattacharya Jelle Bor Mathias Butenschoen Federico Ceccopieri Longjie Chen Vincent Cheung Umberto D'Alesio Miguel Echevarria Yoshitaka Hatta Charles E. Hyde Raj Kishore Leszek Kosarzewski Cédric Lorcé Wenliang Li Xuan Li Luca Maxia Andreas Metz Asmita Mukherjee Carlos Muñoz Camacho Francesco Murgia Pawel Nadel-Turonski Cristian Pisano Jian-Wei Qiu Sangem Rajesh Matteo Rinaldi Jennifer Rittenhouse West Vladimir Saleev Nathaly Santiesteban Chalis Setyadi Pieter Taels Zhoudunmin Tu Ivan Vitev Ramona Vogt Kazuhiro Watanabe Xiaojun Yao Yelyzaveta Yedelkina Shinsuke Yoshida















1]Daniël BoerEditor
2]Chris A. FlettEditor
2,3,4]Carlo FloreEditor
5]Daniel KikołaEditor
2]Jean-Philippe LansbergEditor
2]Maxim NefedovEditor
2,6]Charlotte Van HulseEditor

[Editor]Editor

7]Shohini Bhattacharya
1,2]Jelle Bor 
7b]Mathias Butenschoen
2,8]Federico Ceccopieri
9,10]Longjie Chen
11]Vincent Cheung
4]Umberto D'Alesio
12]Miguel Echevarria
7,13]Yoshitaka Hatta
14]Charles E. Hyde
12,15]Raj Kishore
16]Leszek Kosarzewski
17]Cédric Lorcé
15,18]Wenliang Li 
19]Xuan Li
1,4]Luca Maxia
20]Andreas Metz
21]Asmita Mukherjee
2]Carlos Muñoz Camacho
4]Francesco Murgia
17,22]Pawel Nadel-Turonski 
4]Cristian Pisano
23]Jian-Wei Qiu
24]Sangem Rajesh
25]Matteo Rinaldi
26,27]Jennifer Rittenhouse West
28]Vladimir Saleev
29]Nathaly Santiesteban
1,30]Chalis Setyadi
31]Pieter Taels
7]Zhoudunmin Tu
19]Ivan Vitev
11,32]Ramona Vogt
33,34]Kazuhiro Watanabe
35,36]Xiaojun Yao
2,37]Yelyzaveta Yedelkina
9,10]Shinsuke Yoshida



[1]Van Swinderen Institute for Particle Physics and Gravity, University of Groningen, 9747 AG Groningen, The Netherlands
[2]Université Paris-Saclay, CNRS, IJCLab, 91405 Orsay, France
[3]Dipartimento di Fisica, Università di Torino, and INFN Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy
[4]Dipartimento di Fisica, Università di Cagliari, and INFN Sezione di Cagliari, Cittadella Univ., I-09042 Monserrato (CA), Italy
[5]Faculty of Physics, Warsaw University of Technology, plac Politechniki 1,00-661, Warszawa, Poland
[6]University of Alcalá, Alcalá de Henares (Madrid), Spain
[7]Physics Department, Brookhaven National Laboratory, Bldg. 510A, Upton, NY 11973, USA
[7b]II. Institut für Theoretische Physik, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany

[8]Université de Liège, B4000, Liège, Belgium
[9]Key Laboratory of Atomic and Subatomic Structure and Quantum Control (MOE), Guangdong Basic Research Center of Excellence for Structure and Fundamental Interactions of Matter, Institute of Quantum Matter, South China Normal University, Guangzhou 510006, China
[10]Guangdong-Hong Kong Joint Laboratory of Quantum Matter, Guangdong Provincial Key Laboratory of Nuclear Science, Southern Nuclear Science Computing Center, South China Normal University, Guangzhou 510006, China
[11]Nuclear and Chemical Sciences Division, Lawrence Livermore National Laboratory, Livermore, CA 94551 USA
[12]Department of Physics & EHU Quantum Center, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain
[13]RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA
[14]Department of Physics, Old Dominion University, Norfolk, VA 23529, USA
[15]Center for Frontiers in Nuclear Science, Stony Brook University, Stony Brook, NY 11794, USA
[16]Department of Physics, The Ohio State University, Columbus, Ohio 43210, USA
[17]CPHT, CNRS, Ecole polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France
[18]Stony Brook University, Stony Brook, NY 11794, USA
[19]Los Alamos National Laboratory, Los Alamos, NM 87545, USA
[20]Department of Physics,  Temple University,  Philadelphia,  PA 19122, USA
[21]Department of Physics, Indian Institute of Technology Bombay, Powai, Mumbai 400076, India
[22]University of South Carolina, Columbia, SC 29208, USA
[23]Theory Center, Jefferson Laboratory, Newport News, VA 23606, USA
[24]Vellore Institute of Technology, Vellore, Tamil Nadu 632014, India
[25]Dipartimento di Fisica. Università degli studi di Perugia, and INFN Sezione di Perugia. Via A. Pascoli, Perugia, 06123, Italy
[26]Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
[27]University of California at Berkeley, Berkeley, CA 94709, USA
[28]Joint Institute for Nuclear Research, Dubna, Russia
[29]University of New Hampshire, Durham, NH 03824, USA
[30]Department of Physics, Universitas Gadjah Mada, BLS 21 Yogyakarta, Indonesia
[31]Department of Physics, University of Antwerp, Groenenborgerlaan 171, 2020 Antwerpen, Belgium
[32]Department of Physics and Astronomy, University of California at Davis, Davis, CA 95616, USA
[33]SUBATECH UMR 6457 (IMT Atlantique, Université de Nantes, IN2P3/CNRS), 4 rue Alfred Kastler, 44307 Nantes, France
[34]Faculty of Science and Technology, Seikei University, Musashino, Tokyo 180-8633, Japan
[35]Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139 USA
[36]InQubator for Quantum Simulation, Department of Physics, University of Washington, Seattle, Washington 98195, USA
[37]School of Physics, University College Dublin, Dublin 4, Ireland

















































§ ABSTRACT 
The physics case for quarkonium-production studies
accessible at the  US Electron Ion Collider is described. 













§ INTRODUCTION


The Electron-Ion Collider (EIC) accelerator and detector systems are currently designed following the elaboration of an outstanding physics case aimed at further exploring the nucleon and nucleus partonic structure. The interested reader will find it useful to consult 
reviews of the EIC <cit.>. 








Bound states of heavy quark-antiquark pairs, QQ̅, i.e. quarkonia, allow for a detailed study of  basic properties of quantum chromodynamics (QCD), the theory of the strong interaction. 
Indeed, charmonia and bottomonia have played a crucial role in the establishment of QCD as the theory of the strong interaction, given the clean signature they provide in different observables.
On the theory side, the main origin of the simplifications is the hierarchy m_Q≫Λ_ QCD, with m_Q the mass of a heavy quark, meaning that for processes occuring at this scale (or higher), a perturbative expansion in α_s of QCD 
is allowed.

In parallel, the non-perturbative effects associated with the formation of the bound state can be factorised. 


Heavy quarkonia are multiscale systems. Besides m_Q and Λ_QCD, one needs to consider, in addition, the scale of the typical momentum transfer between heavy constituent quarks (m_Q v), v being the velocity of the heavy quarks in the rest frame of the bound state, and the scale of their binding energy (m_Q v^2), all of which become widely separated in the limit m_Q→∞.

At this point, the non-relativistic nature of the system comes into play, allowing for the development of different effective theories of QCD that attempt to more adequately describe the production of the bound state in the presence of different relevant scales as well as models such as the Colour-Singlet (CS) Model <cit.> or Colour-Evaporation Model (CEM) <cit.>. 





 

The multitude of existing theoretical approaches to describe quarkonium production reflects the fact that unfortunately, up to now, there is no universal physics picture of this process accepted by the community that would provide a satisfactory description of all available experimental data <cit.>. This complicates the 

use of quarkonia as tools for precision studies. Heavy quarkonia nevertheless remain useful to uncover new facets of the structure of nucleons and nuclei which we review in this document. 


In this context,  measurements of various quarkonium-production observables in electron-proton (ep) and electron-nucleus (eA) collisions at the EIC could provide crucial experimental clues to finally settle  the quarkonium-production-mechanism debate.


Important targets for the EIC experimental programme are vector-quarkonium-polarisation observables and cross-section measurements of C=+1 states, like the η_c and χ_c,b. These play a central role in the current debate about the heavy-quarkonium-production mechanism and yet corresponding precise data from ep collisions at HERA are simply lacking.


Such measurements would hopefully clear up the quarkonium-production debate and allow one to fully employ quarkonium data at the EIC as tools.



Before discussing how quarkonia can be used as tools to study nucleons, let us recall that the multi-dimensional structure of nucleons is parameterised by different hadronic functions, which encode the dynamics of partons at different levels of complexity. These span from the one-dimensional (1D) parton distribution functions (PDFs), to the five-dimensional (5D) Wigner distributions –or generalised transverse-momentum-dependent distributions (GTMDs)–, to mention a few. These also incorporate a variety of spin and momentum correlations between the parton (or partons) participating in the hard subprocess and its (their) parent hadron. 
Depending on the considered scattering process and the measured kinematics, different hadronic functions enter the relevant cross sections. 

Among them, let us cite the transverse-momentum-dependent PDFs (TMD PDFs or TMDs), arising from TMD factorisation <cit.> and which provide information on the distribution of partons inside the nucleon as a function of both their longitudinal 
and transverse momentum.

In the case of quarkonium production at transverse momenta, , small compared to their mass and for specific other kinematical end-point regions, new TMD functions related to the produced quarkonium and referred to as shape functions
,  are expected to enter the cross-section formula besides the TMD PDFs of the initial-state hadron(s). This reflects the interplay between radiation of soft gluons and effects of the formation of the QQ̅ bound state. 
Their impact on the phenomenology remains at present unknown.

Much progress has been made in the determination of the above-mentioned PDFs, achieving different levels of success. Currently, the gluon distributions in general remain 
much less explored than their quark analogues. In this context, quarkonia arise as a powerful handle to remedy this situation since, in the vast majority of 

situations,  the QQ̅ pair at the origin of the quarkonium comes from photon-gluon (gluon-gluon) fusion in ep (resp. pp) collisions, whereas deep inelastic scattering (DIS) or Drell-Yan-pair production are sensitive to gluon only through radiative QCD corrections.  
  

However, it has been shown <cit.> that factorisation of observables –cross sections, angular modulations, spin asymmetries,…– in terms of TMD PDFs is less universal than that in terms of 

standard (collinear)

PDFs and that consequently  such a factorisation could be violated in back-to-back-hadron production in proton-proton (pp) collisions. In ep and eA collisions, there is no anticipated violation of TMD factorisation, at least for 

inclusive single-hadron production, so quarkonium measurements  will likely be easier to interpret in terms of gluon TMDs at the EIC 

rather than at hadron colliders.    




Quarkonia are also key players in exclusive reactions. This is not surprising as exclusive meson production involving a hard scale 

is one of the main processes to access Generalised Parton Distributions (GPDs). 

Gluon GPDs are in particular accessible via exclusive 

heavy-quarkonium production <cit.>.
These GPDs provide information on the distribution of gluons inside the nucleon simultaneously as a function of their longitudinal momentum and their transverse position. They also provide information on the  angular momentum of the gluons inside the nucleon, about which very little is known to date. Furthermore, exclusive heavy-quarkonium production near the production threshold was suggested <cit.> as a tool for constraining the gluon condensate in the nucleon, itself linked to the nucleon mass, albeit with some unavoidable model dependence. 



At small momentum fraction, x, the differential exclusive electro- and photoproduction cross sections of quarkonia 


can be expressed in terms of particular products of integrals of GTMDs. 


In single-quarkonium production, when a collinear expansion is applied, the cross section reduces to expressions in terms of GPDs, see for example <cit.>. However, in general, especially beyond single-particle production, it 
provides additional information on GTMDs and offers an opportunity to learn more about the combined three-momentum and spatial distributions of gluons inside a nucleon. 

Moreover, while there is a direct relation between 

exclusive photoproduction case in ep collisions and 
in ultra-peripheral 
pp and pPb collisions (UPCs), 
studies at the EIC would allow one to probe in more details the transverse-momentum dependence of the GTMDs. 





To date, the detector simulations for the EIC physics case connected to quarkonium physics has been limited to   and Υ exclusive production 
as reported in the EIC Yellow report <cit.>. Whereas, as we discussed above, quarkonium production is still the object of intense debates within the community[We guide interested readers to the following reviews <cit.> which address HERA and Tevatron results, to more recent ones <cit.> addressing progress made thanks to the RHIC and LHC data and to to the HEPData database (https://www.hepdata.net/https://www.hepdata.net/), a dedicated repository of quarkonium measurements up to 2012 (http://hepdata.cedar.ac.uk/review/quarkonii/http://hepdata.cedar.ac.uk/review/quarkonii/ to <cit.> and to <cit.> for experimental quarkonium data.], there is no doubt that it can play a crucial role in the scientific success of the EIC. As 
was recently done for the High Luminosity LHC phase <cit.>, we gather in this review what we believe to be the most complete list of quarkonium studies that can be carried out at the EIC along with their motivation.


The document is organised as follows. In Section <ref>, the EIC accelerator system and the first EIC detector, ePIC, as currently envisioned are presented. After a description of the kinematics of lepton-hadron collisions, the importance and a theoretical treatment of QED radiation are discussed, to then end with a note on feed-down from b-quark production in the study of charmonium. In Section <ref>, 
the different theoretical descriptions and measurements related to the production mechanism of quarkonia are presented. First, the various existing theoretical formalisms in collinear factorisation are discussed. Then, the legacy of existing measurements 

and the potential of future measurements at the EIC in constraining these formalisms are presented. Subsequently, it is shown that TMD observables can also teach us about quarkonium formation and polarisation and predictions for the EIC are made. Finally, the effect of final-state interactions on the production of quarkonium in lepton-nucleus colissions is touched upon. Section <ref> focuses on the studies accessible in electron-proton collisions in order to advance our knowledge of the nucleon partonic structure 
and then moves on to studies with nuclear beams, which is a unique feature of the EIC, and which will allow us to make a giant leap forward into a new precision era of the partonic structure of nuclei. 
We underscore throughout this comprehensive review the diverse ways in which the EIC will utilise quarkonia to probe hadronic and nuclear physics and, conversely, will itself be a powerful tool for probing quarkonia. 






§ GENERALITIES ABOUT QUARKONIUM STUDIES AT THE EIC




 §.§ The proposed EIC accelerator system






The EIC

is an upcoming particle accelerator that will deliver intense beams of longitudinally polarised electrons and polarised light nuclei (p, d, ^3He) as well as unpolarised heavier ions, ranging up to uranium. It will produce electron-ion collisions at the highest energy and at the highest rate 
ever achieved. 

The EIC will be constructed in Brookhaven National Laboratory 
using a few key elements of the currently operating Relativistic Heavy Ion Collider (RHIC) <cit.>, such as the hadron ring
and the RHIC Electron-Beam-Ion-Source (EBIS) <cit.>.
The collider will be supplemented with a 
new electron ring, which will contain continuously injected,  polarised electrons with an energy from 5 GeV up to 18 GeV. 
The coverage in centre-of-mass energy will range from 
28 GeV to 141 GeV for 
 lepton-proton collisions, while for lepton-ion collisions an upper energy of 89 GeV/nucleon will be reached.



The 
expected instantaneous luminosity depends on the centre-of-mass energy and will range from 10^33 to 10^34 cm^-2s^-1 for electron-proton collisions,  with a maximum value expected for = 105 GeV. For electron-ion reactions, it will be on the order of 10^34 cm^-2s^-1. Such figures will correspond to integrated luminosities 
of the order of 10 to 100 fb^-1 per year. The designed average polarisation of electron, proton and ^3He beams is 
of the order of 70%.  
 








At present, the installation of a first EIC detector is foreseen 
at interaction point 6 (IP6).

A second interaction point (IP8) can, at 
any stage, host a second and complementary detector. 
The second interaction point could accommodate a design with a secondary focus, which in combination with forward spectrometry would
allow for an extension of the acceptance towards the detection of particles at very small polar angles.



The interaction points will re-use the existing large detector halls, currently occupied by the STAR and sPHENIX experiments. 








The first collisions at the EIC are expected in the early 2030s. 





 §.§ The proposed EIC detector









  §.§.§ Requirements for an EIC detector in the context of quarkonium studies





The specification of an EIC detector is determined by the kinematics of the electron-ion scattering (see fig:eic:generic:det) and the observables and processes of interest. It should address the full range of physics outlined in the 
EIC White Paper <cit.>, the NAS report <cit.> and the EIC Yellow Report <cit.>.
The basic requirements include <cit.> 4π hermeticity with large 
acceptance in pseudorapidity, η, of about -4 < η < 4, very good momentum resolution both in the central, forward and backward regions, 
very good energy resolution in electromagnetic calorimeters and particle identification capabilities up to 50 GeV 

in momentum. Such a setup allows 
one to study processes over a wider range of four-momentum transfer Q. In addition, measurements of heavy-flavour hadron production demand a microvertexing detector that provides good impact-parameter resolution. 


The detector technologies and configuration implementation will be known once the detector design is finalised.
However, existing high-energy experiments (for example ALICE at the LHC and STAR at RHIC) indicate that an EIC detector that fulfils the aforementioned requirements will have capabilities for  and Υ(nS) measurements via their e^+e^- decay channel  <cit.>. 
The precision of quarkonium reconstruction will strongly depend on the hardware configuration. For example, an internal silicon tracker could generate additional combinatorial background arising from conversions γ→ e^+e^-, limiting precision for low-mass quarkonia at low . Moreover, the energy loss of electrons due to Bremsstrahlung in the detector material deteriorates the mass resolution. It may complicate, if not make impossible, separation of the , 
and states. Measurements of other quarkonium states (for instance χ_c or χ_b) add constraints for the experimental apparatus. Studies of decays involving photon radiation (e.g., χ_c(1P) → + γ) would require an electromagnetic calorimeter able to isolate a soft photon and measure its energy with appropriate resolution. In addition, a muon detector would significantly extend capabilities for quarkonium studies. This is 
briefly discussed in Sec. <ref>.

Three different designs, ATHENA <cit.>, CORE <cit.> and ECCE <cit.>, were proposed.
The main difference between the ATHENA and ECCE design consists of the magnet, providing respectively a 3.0 T and 1.4 T magnetic field.
The distinguishing characteristic of the CORE detector is the compactness of the detector, obtained through exploitation of technological advances.
From the proposed designs, the ECCE proposal was selected as baseline for the first EIC detector, 
with improvements to the proposal 
at present under development.
This first EIC detector received the name electron-proton/ion collider (ePIC) detector.
A description of the ePIC detector in its current design state is given below.































































 
 
  
   
    
    
     
      
       
    
     
      
    
     
      
        




















































  §.§.§ The ePIC detector





The central barrel of the ePIC detector, as currently 
envisioned, is depicted in figure <ref>. Here, the hadron
beam comes in from the left and defines the forward-going direction.
The central barrel is around 
10 m long and 
5 m in diameter, providing a full coverage
in azimuthal angle and a coverage in polar angle between 0^∘ and 178^∘, corresponding to a pseudorapidity coverage
between -4 and 4.
In addition to the barrel detector, detectors in the far-forward and far-backward regions are foreseen.
The far-backward region will contain a luminosity monitor and two detectors to tag
low-Q^2 events. The far-forward region will contain a series of detectors aimed at detecting particles produced close to the beam line and as such will be instrumental 
to the reconstruction of an extensive set of diffractive processes and tagged measurements, such as proton reconstruction in exclusive processes, tagging of the two spectator protons when investigating 
the neutron structure through lepton-^3He interactions and tagging of respectively the neutron and Λ-baryon decay particles when probing the pion and kaon 
structure 
in lepton-proton interactions. 
The far-forward system will consist of a B0 spectrometer, 
containing an electromagnetic calorimeter and trackers for respectively the tagging of photons and reconstruction of charged particles,
Roman Pots 
and off-momentum detectors, 
performing charged-particle reconstruction, 
as well as 

Zero-Degree Calorimeters, 
capable of detecting photons and neutrons.

In the central barrel, track and vertex reconstruction will be performed by
silicon monolithic active pixel sensors placed close to the beam line and interaction point, while at a further distance
micro-pattern gaseous detectors (
micro-Resistive Well 
and Micro-Mesh Gaseous Structure) and AC-coupled low-gain
avalanche diodes will contribute to track reconstruction.
The tracking detectors will be embedded into a 1.7 T magnetic field.
Such a setup will provide the momentum resolution needed to 
fulfil the EIC physics programme.

Electromagnetic calorimeters cover the backward, central and forward regions of the central barrel, providing
electron and photon detection as well as hadron suppression.
In the backward region, a high-precision lead-tungstate calorimeter read out by silicon photo-multipliers is foreseen.
The detector will be 
critical to the reconstruction of (scattered) electrons,
improving the reconstruction precision over that obtained from tracking detectors only,
and in the identification of these electrons, by suppressing the background contribution strongly.
This contribution originates mostly from charged pions.
In the central region, 
a lead-scintillator imaging calorimeter is foreseen.

For the forward region, an electromagnetic calorimeter will be integrated with the forward hadronic calorimeter.
The system 
focuses on the containment of high-energetic particle showers while at the same time providing a
good energy resolution for lower-energetic particles.
Particle identification requires 
a good position resolution, in particular in the electromagnetic calorimeter. This will be provided by constructing the electromagnetic calorimeter
out of segments, of scintillating fibres embedded in tungsten powder, 
smaller than the Molière radius. This will also result in a good shower 
separation at high pseudorapidity.


In the central region, a hadronic calorimeter 
 will allow for the detection of neutral hadrons and as such will improve the resolution of jet reconstruction. 
Given the good momentum resolution of the central trackers, the central hadronic calorimeter system will not have
an impact on the reconstruction of charged particles.
The forward hadronic calorimeter, which forms an integrated system with the electromagnetic calorimeter, will
consist of layers of alternating tiles of scintillating material and steel, while towards the end of the detector
the steel is replaced by tungsten in order to serve as 
tail catcher of the shower and thus
maximise the interaction length within the available space.

Also in the backward region, an hadronic calorimeter will be installed, with the aim to serve as tail catcher of particle signals.

Detectors based on the detection of Cherenkov light will be used for the identification of
charged pions, kaons and protons, while also contributing to electron identification.
In addition, the 
aforementioned AC-coupled low-gain
avalanche diodes will provide particle identification in the low-momentum region, below ∼ 2 GeV, based on the detection of the
time of flight of a particle.
In the backward region, a proximity-focusing ring-imaging Cherenkov (RICH) detector with aerogel as radiator 
will be used. 
Because of the tight space constraints, a
DIRC – detection of internally reflected Cherenkov light – detector will be incorporated  in the central region.
The forward region will contain a dual RICH detector, with 
an aerogel radiator  for the low-momentum particles
and C_2F_6 for the high-energetic ones, covering the momentum range up to 50 GeV. 

No muon detectors are foreseen for the ePIC detector. While first studies, performed for the ATHENA and ECCE proposals, indicate that the reconstruction of 
 mesons from exclusive processes through their e^+e^- decay should be possible with the ePIC detector, there are neither studies for other quarkonium states nor for inclusive or semi-inclusive processes. 
Here, dedicated	muon detectors might be	needed.	This is	discussed in the following sub-section.




  §.§.§ The case for a muon detector for quarkonium studies at the EIC




Measurement of vector-quarkonium production using their di-muon decay provides significant benefits. The energy loss of muons due to interactions with detector material is much smaller than that of electrons. This leads typically to a better momentum resolution 
of the muons than of the electrons, and 
therefore the resolution of the quarkonium mass reconstructed in the μ^+μ^- channel is better compared to the e^+e^- one. The LHCb and CMS experiments provide a case in point as the performance of their muon detectors facilitated a rich and fruitful quarkonium physics program, which included that of 
, , and other quarkonium states such as the χ_c and χ_b via their radiative decays into vector quarkonia. Additional measurements via the μ^+μ^- decay channel would also essentially double the available statistics as the branching ratios into μ^+ μ^- and e^+e^- are 
nearly the same and enable analyses of rare decays (for example, χ_c →μμ). With a proper design,  

 studies via the di-muon channel benefit from a lower combinatorial background, thus improving the statistical precision of the measurement. In addition, they provide a cross check 
 of the e^+e^- results, which should in turn reduce  systematic uncertainties.  

In summary, a muon detector would significantly extend capabilities for quarkonium studies at the EIC. The present ePIC design does not consider muon-identification instrumentation, but

possibilities for an enhanced muon identification can be investigated for ePIC. Moreover, 
the incorporation of dedicated muon-identification detectors in
the design phase of the 2^ nd EIC detector can 


vastly extend quarkonium measurement capacities in the manner described above. 





 §.§ Kinematics and QED radiative corrections









  §.§.§ Kinematics of electron-hadron reactions





In this section, we collect 

basic kinematical definitions useful for the description of lepton-hadron reactions. The next section is devoted to how QED radiative corrections on the lepton side can affect the resolution on various kinematic variables and to possible ways to address this problem.


Let us consider the inclusive production of an identified hadron H, which in the context of this review is most likely to be a quarkonium, in electron-nucleon (eN) scattering:


    e(ℓ) + N(P_N) → e(ℓ^') +  H(P_ H) + X,

For electron-nucleus (eA) scattering, the momentum P_N usually denotes the average momentum of a single nucleon. 
 Depending on the experimental possibilities, one can  
tag the outgoing electron with the momentum ℓ^' or 

consider the reaction inclusive w.r.t. the final-state electron. If the momentum ℓ^' has been measured, 
one can define the momentum transfer q=ℓ-ℓ^' with q^2 = -Q^2 and the following Lorentz-invariant kinematic variables become experimentally accessible:

    = Q^2/2P_N· q ,      y = P_N· q/P_N·ℓ  ,      z = P_N· P_ H/P_N· q ,

where z is referred to as the elasticity and y as the inelasticity 
that should not be confused with the rapidity.



Among frame-dependent variables, one usually distinguishes the transverse momentum of the hadron P_T H in the laboratory frame (see 
fig:eP_frames (left)), in which the initial electron e(ℓ) and 
nucleon N(P_N) collide head on, defining the Z (collision) axis, from the transverse momentum P_T H^*  of the hadron H in the photon-hadron frame (see 
fig:eP_frames (right)), where three-momenta q and P_N are aligned with the Z axis of this frame[Different photon-hadron frames are related by 
a boost along the Z axis. In particular, one can adopt the photon-nucleon centre-of-momentum frame where q+ P_N=0. The transverse components of the momenta are the same in all photon-hadron frames.]. The word “photon” in the frame name 
specifically refers to the one-photon-exchange approximation between the electron and the hadronic part of the process.  In this review, we will often use the simplified notation  for the absolute value of the transverse momentum of the produced hadron: =| P_T H| or ^*=| P^*_T H|.

Even if the colliding particles are 
unpolarised, there could always be some dependence of the cross section on the azimuthal angle ϕ_T (or ϕ^*_T) formed 
by the vector 
P_T H (or P^*_T H) and the plane spanned by the initial () and final (')  lepton three-momenta (fig:eP_frames), due to the exchanged-photon polarisation. If the initial 
nucleon and/or electron have transverse polarisation, additional angular modulations of the cross section, related to the direction(s) of the transverse spin vector(s) of the colliding particles, can be generated. The transverse polarisation vector of the initial nucleon is denoted as S_T and the angle of this vector with respect to the lepton plane in the photon-hadron (resp. laboratory) frame is generally indicated as ϕ^*_S (resp. ϕ_S).   

If the recoil effects of the photons which can be emitted by the initial and final electrons during the scattering process (QED radiative corrections) 
are neglected, then the four-momentum of the exchanged photon is simply q=ℓ-ℓ^' as stated above. In such an approximation, the variables of eq:SIDIS:kin-var as well as the frame-dependent variables, such as P_T H^*, can be directly computed  from the measured energy and momentum of the scattered electron. However, such a QED Born approximation 
might be insufficient for precision studies. 
Section <ref> is devoted to this issue. 

The regime of the process of eq:process:e-pA, when the quasi-real-photon approximation can be applied to the exchanged photon,  
when Q is negligible compared to the hard scale (m_Q, P_T, P_T^*, ...), is commonly referred to as photoproduction, while the regime with Q being the hard scale, or among the potential hard scales, is called leptoproduction or (semi-inclusive) deep inelastic scattering (SIDIS). Experimentally, photoproduction is usually defined by a fixed cut on the photon virtuality, Q < 1 GeV. 

Beside the 
well-known regimes of photoproduction and leptoproduction (or SIDIS)
, which a priori require 
setting some constraints on Q^2, it appears very valuable for quarkonium studies to consider measurements where Q^2 is fully integrated over. Such yields then contain the contributions from both quasi-real and off-shell photons. This proposal is described in 
section <ref>. 




As it was mentioned in the introduction, polarisation observables play an important role in quarkonium physics. The polarisation parameters of a spin-1 heavy quarkonium λ_θ, μ_θϕ and ν_θϕ parametrise the angular distribution of decay leptons in the quarkonium rest frame:

    dσ dΩ∝ 1 + λ_θcos^2θ + μ_θϕsin 2θcosϕ + ν_θϕ 2sin^2 θcos 2ϕ.

These parameters depend on the orientation of the axes of the coordinate system chosen in the quarkonium rest frame with popular frame choices such as the Helicity, Collins-Soper, Gottfried-Jackson and target frames (see Section 2.3 of <cit.>). The same definition of polarisation parameters holds for the case of exclusive production of a vector quarkonium. 



  §.§.§ On the importance of QED corrections




The possibility to make a distinction between the photoproduction and 
electroproduction (or SIDIS) 
regimes, together with the rich phenomenology provided by measurements differential in the variables , y, z as well as 
^* and ϕ^*_T, has 
always been considered as an advantage of lepton-hadron reactions over hadron-hadron ones. 



However, the emission of photons by initial- and final-state leptons modifies the relation between the momentum ℓ^' of the final-state lepton measured in the detector and the four-momentum q of the photon exchanged between the leptonic and hadronic parts of the process in eq:process:e-pA (see the Fig. <ref>) which in turn modifies the Lorentz-invariant variables 
(eq:SIDIS:kin-var) as well as P_T H^* and ϕ^*_T. Beyond the Born approximation of QED, this relation is no longer simply q=ℓ-ℓ^' but includes the recoil from emitted photons. For strictly inclusive DIS measurements, as opposed to SIDIS, the application of QED radiative corrections boils down to an overall radiative correction factor to the cross section differential in x_B and the inelasticity y <cit.>. In the SIDIS case, fully differential Monte-Carlo computations have to be performed, using dedicated tools such as DJANGOH <cit.>.


Recently, it has been shown <cit.> that QED radiative corrections fundamentally limit the accuracy of 
SIDIS measurements, in particular for the kinematic regime where the TMD factorisation 

is needed. In this context, a new approach to their treatment has been proposed. The QCD factorisation for the SIDIS cross section was historically discussed in the photon-hadron frame. However, as it was mentioned above, the collision-induced photon radiations change both the direction and magnitude of the exchanged virtual photon, making the photon-hadron frame and the quantities related to it only approximately defined.  
The ambiguities in the definition of kinematic variables in photon-hadron frame can impact our ability to extract the TMDs and, in particular, to use the angular modulation in ϕ_T^* to separate contributions of different TMD PDFs and fragmentation functions (FFs). Since the QED radiations differently affect the determination of the angles ϕ^*_T and ϕ^*_S (fig:eP_frames), this can affect the determination of various azimuthal (spin) asymmetries <cit.>. In addition to the uncertainty of the “photon-hadron” frame, the collision-induced photon radiations also change 
the true values of x_B and Q^2. 


 Although the effects of the QED radiations could be calculated perturbatively, the main point of concern are those QED radiative correction effects 
 which are logarithmically enhanced due to the collinear and infrared sensitivity coming from the smallness of the electron mass m_e 
 compared to all the other scales of the process.  Omitting these effects may lead to significant uncertainties in some kinematic regimes 
 where 
 a wide phase space is available for collision-induced radiations, such as those relevant to the study of small-x physics or for 
 two-scale observables described by TMD factorisation. 

In Ref. <cit.>, it has been argued that a combined QCD+QED factorisation can be performed such that the exchanged photon momentum q is not fixed by the measured ℓ-ℓ', but rather has a range of values 
to be integrated over.  The range is determined by the observed momentum of the scattered lepton for inclusive DIS and the momenta of both the scattered lepton and the observed final-state hadron for 
SIDIS. 




The approach consists 
of using collinear factorisation to take into account the collision-induced-QED-radiation effects which are enhanced by large logarithms of either Q/m_e, | P_T H|/m_e or |_T'|/m_e, while either collinear or TMD factorisation can be used 
to account for QCD contributions depending on the hierarchy between the | P_T H-_T'| and the hard scale Q. For the SIDIS process 
of eq:process:e-pA on a proton target, the hybrid factorisation formula is given by <cit.>: 


    E_ℓ^' E_P_ Hdσ_ SIDIS/d^3^' d^3 P_ H   ≈   ∑_a,b ∫_ζ_ min^1dζ/ζ^2  D_e(ℓ^')/b(k')(ζ,μ_F^2)
    ∫_ξ_ min^1 dξ  f_a(k)/e(ℓ)(ξ,μ_F^2)  
    
          × [E_k^'E_p_ Hd σ^ap[a(k)+p(P) → b(k^')+ H(p_ H)+X]/d^3 k^'d^3 P_ H]_k=ξℓ,k^'=ℓ^'/ζ,

where a,b=e, e̅, γ
, and where the active lepton/photon momenta entering or leaving the hard collision 
are defined as k=ξℓ and k'=ℓ'/ζ with collinear momentum fractions ξ and ζ, and μ_F is the factorisation scale. The process-independent lepton distribution functions (LDFs) f_a/e(ξ) and lepton fragmentation functions (LFFs) D_e/b(ζ) 

in Eq. (<ref>) resum logarithmically-enhanced QED contributions in the limit when the hard scale, max (Q,| P_T H|,|_T'|), is much larger than m_e. The non-logarithmically-enhanced part of QED radiative corrections can be included into dσ^ap 
order by order in powers of α_ em. 

The differential cross section 
dσ^ap in the second line of eq:SIDIS-factor-QED can be further factorised by TMD or collinear factorisation in QCD depending on if the observed lepton and hadron are in the back-to-back regime or not.  As it has been demonstrated in Refs. <cit.>, the transverse-momentum broadening from the collision-induced QED radiations is much smaller than the TMD effects from QCD. Factorising out QED radiations using collinear LDFs and LFFs as done in the eq:SIDIS-factor-QED is therefore a good approximation. eq:SIDIS-factor-QED is valid up to Leading Power (LP), that is up to 
power corrections 
scaling as the inverse of the hard scale. Note that the same kind of equation holds in the case of e-A collisions. Note also that eq:SIDIS-factor-QED does not account for possible hadronic/resolved contributions from the photon.



Due to the smallness of α_ em, σ^ap in eq:SIDIS-factor-QED can be approximated  by its QED Born 
order, σ^ap,(0) 
with a=b=e.  This lowest order cross section is 
the same as the SIDIS cross section without QED radiation which can be parametrised in terms of the usual SIDIS structure functions <cit.> but with different kinematics: ℓ→ k=ξℓ and ℓ'→ k'=ℓ'/ζ. Consequently, the exchanged-virtual-photon momentum between the scattered lepton and the colliding hadron is modified as q=ℓ-ℓ' → k-k'=ξℓ - ℓ'/ζ. 
By neglecting higher order QED 
contributions to σ^ap in eq:SIDIS-factor-QED, the SIDIS cross section with the 
collision-induced QED radiation can thus be obtained from the same SIDIS cross section without QED radiation plus the knowledge of the universal LDFs and LFFs. 









 §.§ On the importance of b feed down



 An important and subtle  
 concept 
 needed to understand the quarkonium-production mechanism is the knowledge of feed downs. For instance, as shown in Ref. <cit.>, in the case of J/ψ photoproduction at HERA, not all the J/ψ are 
 produced by the hard scattering. Indeed, a non-negligible fraction of the J/ψ mesons produced at large P_T comes from the decay of a b quark. 
 fig:H1-data-ratio-b-FD shows the fraction of J/ψ coming from such a b feed down (also referred to as non-prompt yield) as a function of P_T^2 in the H1 kinematics. We guide the reader to Appendix A of Ref. <cit.> for more information about how it was estimated. 
One sees that the fraction of non-prompt J/ψ steadily grows to reach over 40% of the J/ψ yield at the highest reachable P_T ≲ 10 GeV.

Although the top energy of the EIC will be at most at √(s_ep) = 140 GeV, 
given the much higher luminosity of the EIC compared to HERA, the W_γ p reach[
W_γ p = √(s_γ p) designates the energy in the centre-of-mass of the photon-proton system.] might be such that, at high P_T, similarly large non-prompt fractions could be observed.


With this respect, further dedicated studies are necessary.




§ EIC TOOLS FOR QUARKONIUM STUDIES 





 §.§ Quarkonium-production mechanisms











As aforementioned, to justify the application of perturbative QCD to the studies of identified hadron production, the observables 
should 
involve some scale μ≫Λ_ QCD, such that α_s(μ)≪ 1. 
In such cases, the 
cross section can 
be factorised (up to power-suppressed corrections in μ) into a product or convolution of a short-distance part, which is meant to be computed perturbatively as a series in α_s(μ) and long-distance factors. The latter 
comprise (TMD-) PDFs of incoming hadrons and non-perturbative quantities which describe the hadronisation of partons produced at the short-distance/perturbative stage of the process into an observed final-state hadron. 

The treatment of hadronisation differs for hadrons containing light quarks in the “naive” quark-model picture of these states as opposed to quarkonia, the primary component of which is expected to be a 
QQ̅ Fock state with the same quantum numbers as 
quarkonium.  In the case of hadrons composed of light quarks or heavy-flavoured hadrons like D and B mesons, commonly denoted H_Q  in this review, in which relativistic (“light”) degrees of freedom play an important role, the hard-scale μ is ∼ p_T ≫ m_ H_Q and the “final-state” long-distance part of the cross section is usually 
encapsulated in a fragmentation function (FF).  Due to the importance of light degrees of freedom, the FFs of such hadrons can not be computed perturbatively and they are parametrised at some starting scale 
μ_0, on the order of 1 GeV, with parameters fitted to reproduce experimental data, see e.g. 
<cit.> and <cit.> for fits of respectively light and heavy-flavoured hadrons. 

For hadrons containing two tightly-bound heavy quarks, such as “standard” charmonia (η_c, J/ψ, χ_c, ψ(2S), …) and bottomonia (η_b, Υ(nS), χ_b,…), denoted hereafter by , a deeper understanding of hadronisation is believed to be possible. 

The overall success of non-relativistic potential models in the description of the mass spectrum of these states implies 
that the contributions of QCD Fock states containing gluons or light quarks is suppressed by the powers of the average velocity v of the heavy quarks in the bound state 
compared to that  of the simplest Fock state with only one heavy QQ̅ 
pair. 







The typical squared velocity v^2 is estimated in potential models to be ∼ 0.3 for charmonia and ∼ 0.1 for bottomonia, which turns it into an no additional n here useful small parameter with respect to which the observables can be expanded. 
The different existing models of quarkonium production <cit.> follow more or less closely the above observation which yields to somehow disparate predictions for some production observables. We review below the main features of three of the most popular ones which will follow us throughout this review.



  §.§.§ NRQCD & CSM




In the non-relativistic QCD (NRQCD) factorisation formalism <cit.>, the cross sections and decay rates are expanded in powers of α_s(μ) and v^2. At each order of the v^2 expansion, the short-distance part of the observable describes the production or annihilation of the QQ̅-pair in a colour-singlet or colour-octet state with a particular value of spin, orbital and total angular momentum.

 The hard scale, μ, for the short-distance part 
 can be the heavy-quark mass m_Q, 
 or any other larger scale not comparable to Λ_ QCD, justifying the perturbative calculation of this factor. The corresponding long-distance part 
 of the cross section is 
 a number called the  Long-Distance Matrix Element (LDME) 
 which, for the production case, can be written up to conventional colour and spin normalisation factors, which we omit for the sake of clarity, as:

    ⟨ O^ [i] ⟩∝∑_X_s 0( O_i^† Y_n^†)^ab(0)+X_s+X_s( Y_n  O_i)^ba(0)0,

where it is implied that any final state X_s containing light quarks and gluons can be produced together with the quarkonium . 
The factors Y_n in 
eq:LDME_def contain Wilson lines along the light-like direction n needed for the gauge invariance of the Colour-Octet (CO) LDMEs. The structure of the colour indices, ab, connecting the amplitude and complex-conjugate amplitude in eq:LDME_def reflects the process-dependent configuration of the Wilson lines in the factors Y_n. The local NRQCD operators O_i contain heavy-quark and antiquark fields[Denoted as χ and ψ in NRQCD.] and are labelled in the same way as the simplest Fock state | QQ̅[i] + X_s ⟩ which this operator can excite from the vacuum. The spectroscopic notation of the label i=^2S+1L_J^[1,8] is used to denote the total spin S, the orbital angular momentum L, the total angular momentum J and the singlet (CS, ^[1]) or octet (CO, ^[8]) colour quantum numbers of the heavy-quark pair. With these conventions, the complete traditional notation for the LDME becomes: ⟨ O^[ ^2S+1L_J^[1,8]] ⟩.  

NRQCD velocity-scaling rules <cit.> lead to the assignment of the O(v^m) suppression to LDMEs, thus allowing us to truncate the velocity 
expansion at some fixed order in v^2. Usually the contributions associated with the LDMEs up to Next-to-Next-to-Leading Order (NNLO) 

in v^2 (O(v^4) relative to the LDME of the ^3S_1^[1] state) are taken into account in phenomenological studies. 
This means that, besides the colour-singlet QQ̅ states, 
the colour-octet states ^1S_0^[8], ^3S_1^[8] and ^3P_J^[8] can contribute to J/ψ production, for example. 

For S-wave quarkonia, 
the expansion limited to the leading order of v^2 corresponds to the colour-singlet QQ̅-state with the same quantum numbers as those of . 
The colour-singlet model (CSM)
 <cit.>  for the production of these states is nothing but the truncation of the v^2 expansion at this order. The CS LDMEs can be estimated from potential-model wave functions <cit.>, while their accurate estimation from 
ℓ^+ℓ^- decay rates of  is rendered complicated by large NNLO QCD corrections <cit.> to the decay width. However, the CSM is not sufficient theoretically <cit.> 
for the description of the production of the P-wave quarkonia, such as χ_c,b, beyond LO in α_s and can not describe inclusive hadroproduction p_T 
spectra of charmonia and bottomonia at high 
 <cit.>. Nevertheless, the NNLO corrections in α_s to the short-distance part of the CSM cross section, only partially computed so far, may decrease 
the existing large discrepancy between the CSM and the data from Tevatron and the LHC  <cit.>. This point is still under debate <cit.>. 




 In contrast to the hadroproduction case described above, in 
 (prompt) inclusive photo- and electroproduction of heavy quarkonia, which are relevant for the EIC experimental program, the CSM has been expected <cit.> and proven to be able to account 
 for a large fraction of the observed cross section <cit.> even up to the highest reachable . Estimates varying from 50% <cit.> to almost 100% <cit.> can be found in literature. 






 
 


For the exclusive photo- and electroproduction of single J/ψ or Υ(nS), the CS contribution is 
also expected to be strongly dominating. In 
such exclusive reactions, no final-state radiation (X_s) is allowed and the NRQCD operators containing a CO QQ̅ pair can only couple 
to the higher Fock-state contributions in the expansion of the physical quarkonium eigenstate, which are velocity-suppressed, e.g. |J/ψ⟩= O(1)|cc̅[^3S^[1]_1]⟩ + O(v)|cc̅[^3P^[8]_J]+g⟩+… The matrix elements of the gauge-invariant CO operators which in principle can contribute to exclusive photoproduction, e.g. ψ^† (g_s E· D) χ where E is the chromoelectric field and D is the QCD covariant derivative, can be estimated[The scaling for D is O(v) and the scaling for g_s  E is O(v^3) so together with the O(v) suppression of the |cc̅[^3P^[8]_J]+g⟩ component of |J/ψ⟩, one obtains O(v^5).] to scale at least as O(v^5) at the level of the amplitude using the velocity scaling rules <cit.>. In Ref. <cit.>, the same conclusion 
has been made about the CO contributions to the matrix elements of the operator ψ^† D^2 χ = ψ^†∇^2 χ + ψ^† (g_s A·∇) χ + …, which are more suppressed 
than the CS relativistic corrections ⟨ J/ψ | ψ^†∇^2 χ |0⟩∼∇^2 Ψ(0)∼ O(v^2). Therefore, taking into account CS relativistic corrections to exclusive vector-quarkonium photoproduction is currently considered to be more important <cit.> than taking into account the CO corrections.    


Another success <cit.> of the CSM at NLO in α_s is the description of the prompt η_c hadroproduction, measured by LHCb <cit.>. 
However, such a success of the CSM to describe this data set, both at moderate ∼ m_η_c and for ≫ m_η_c is problematic for NRQCD. Indeed, from heavy-quark-spin-symmetry (HQSS) arguments, one expects the CO contributions to η_c cross section at ≫ m_η_c to be on the same order of magnitude as that previously found to describe J/ψ data at similar .




As 
aforementioned, at higher orders in the v^2 expansion, the CO LDMEs contribute, but at present they are treated as free parameters and are adjusted to describe experimental data. Besides order-of-magnitude constraints from O(v^n) scaling and HQSS constraints,  the progress on their theoretical calculation has been limited so far. Recently new expressions for LDMEs in terms of potential-model 
quarkonium wave functions and certain chromoelectric-field correlators have been proposed  in the potential-NRQCD (pNRQCD) formalism in the strongly coupled regime <cit.>. These relations can be used to reduce number of free parameters in the fit under the assumption m_Q v^2 ≪Λ_ QCD. Currently, the advantage of using pNRQCD compared to conventional NRQCD fits is still under debate as well as its applicability, since m_Q v^2 is naively not much smaller than Λ_ QCD . 

In 
Section <ref>, we describe existing phenomenological fits of LDMEs within collinear factorisation, commenting on their successes and shortcomings in more details. Unfortunately at present time there is no single set of LDMEs which can 
satisfactorily describe the charmonium e^+e^- annihilation, γγ fusion, hadro- and photoproduction data together with polarisation observables in the framework of NRQCD factorisation at NLO in α_s, which is a serious problem for the NRQCD factorisation approach. For the case of bottomonia, we lack 
photoproduction, e^+e^- annihilation and γγ fusion data, which prevents us from checking the process-independence of LDMEs for the bb̅ family. Another important task for the EIC, 
in connection with the clarification of the quarkonium-production mechanism, is to perform the first measurement of χ_c0,1,2 and η_c inclusive photoproduction cross sections. 
In this context, we discuss corresponding phenomenological predictions in Section <ref>. Such measurements will be complementary to 
those of χ_c and η_c hadroproduction to check the process-independence of the corresponding LDMEs. 

 
 Data at high ≫ m_, where CS and CO contributions behave differently, are potentially very discriminant for LDME fits. This calls for improvement of the perturbative accuracy of the short-distance part 
 since, at large ,
terms proportional to α_s^n+kln^n 
/m_ H appearing in the perturbative series for the short-distance part of the cross section both at LP in 
 and in  power-suppressed corrections at ≫ m_ need to be tackled.
These potentially large terms can be resummed using the formalism of FFs, perturbatively evolving with the scale μ∼. At LP, this formalism is analogous to the FFs for light hadrons mentioned in the beginning of this section, with 
a sole but important difference, namely that at the starting scale μ_0∼ m_ the FF 
is assumed to be factorised into a short distance part and a LDME. We refer to e.g. <cit.> as an example of the NLO study of this type as well as Refs. <cit.> at LO. At Next-to-Leading Power (NLP), new contributions with the QQ̅ pair as a whole participating in the fragmentation process 
appear <cit.>. These corrections 
seem to influence not only the cross section 
but also the evolution of leading-power FFs <cit.>. However, the effect of this corrections on cross sections and the polarisation is still under investigation in particular for the EIC phenomenology where the  reach, limited to roughly 15-20 GeV, might not be large enough for these to be relevant.




  §.§.§ CEM & ICEM



Given the above mentioned phenomenological problems along with others which we review later, NRQCD factorisation at fixed order in v^2 and α_s is not completely satisfactory. Due to its simplicity, the Colour Evaporation Model (CEM), introduced in Refs. <cit.> remains an attractive alternative mechanism to explain the formation of quarkonium. As the CEM is inspired from quark-hadron duality, one postulates that any QQ̅ pair produced at short distance with invariant mass M_QQ̅ less than the invariant mass of a pair of lightest mesons (H_Q) with open-heavy flavour Q (e.g. D^0 mesons in the case of charmonia) has to hadronise into one of the quarkonia below this heavy-flavour-production threshold with some universal probability. In the CEM, this probability, commonly denoted as F_ for the quarkonium state , is taken to be independent of spin, orbital momentum and colour quantum numbers of the pair, and is fit as a free parameter. 

In the improved CEM (ICEM) <cit.>, the kinematic effects  arising from the mass difference between the QQ̅-pair produced at short distance and the final-state quarkonium is taken into account, which roughly models the effects of soft-gluon emissions at hadronisation stage. This is done through the rescaling of the three-momentum of the pair by the mass ratio, so that the direct quarkonium-production cross section in collisions in the ICEM is given by <cit.>:

    σ   =    F_∑_i,j∫^2m( H_Q)_M_dM_QQ̅ dx_i dx_j  f_i(x_i,μ_F)f_j(x_j,μ_F)·σ̂_ij→ QQ̅(x_i,x_j, p_QQ̅,μ_R,μ_F) |_ p_QQ̅ = M_QQ̅/M_ P_ ,

where i and j are q, q̅ and g such that ij = qq̅, qg, q̅g or gg, 
x_i,j is the momentum fraction of the parton, f(x_i,j,μ_F) is the parton distribution function (PDF) in the proton as a function of x_i,j at the factorisation scale μ_F. Finally, σ̂_ij→ QQ̅ are the parton-level cross sections for the initial states ij to produce a QQ̅ pair of momentum p_QQ̅ at the renormalisation scale μ_R. 
In the ICEM, the invariant mass of the QQ̅ pair, M_QQ̅, is integrated from the physical mass of quarkonium M_ to two times the mass of the lightest open heavy Q-flavour meson m( H_Q). In the traditional CEM, see e.g. <cit.>, the value of 2m_Q is used as the lower limit of mass-integration instead of M_ and the momentum-shift due to the mass-difference between the QQ̅-pair and the quarkonium is neglected.

We emphasise  that the physical picture of the (I)CEM is opposite to NRQCD in the sense that the CS contributions play no special role at all. This assumption makes CEM incapable of describing observables where CS states are clearly dominating, e.g. the prompt hadroproduction of J/ψ pairs <cit.> and the e^+e^-→ J/ψ + cc̅ cross section <cit.>.
However, the (I)CEM still provides a reasonable description of single inclusive prompt quarkonium hadroproduction <cit.> although the model is not capable to describe ∼ m_ and ≫ m_ simultaneously even at NLO <cit.>. 
 
Recent ICEM calculations <cit.> have considered the polarisation in hadroproduction. Polarised production of quarkonium in these calculations restricts the final state quark-antiquark pair to be in the desired spin state, thus implicitly assuming that soft gluons are decoupled from heavy-quark spin. The polarisation parameters are then calculated in terms of the spin matrix elements σ_i_z,j_z. In these matrix elements, the quarkonium is assumed to have J_z = i_z when calculating the scattering matrix element, ℳ. The quarkonium is assumed to take J_z=j_z in calculating the conjugate, ℳ^*. The polar anisotropy (λ_θ), defined in the Eq. (<ref>), is given in this model by <cit.>

    λ_θ   =   σ_+1,+1-σ_0,0/σ_+1,+1+σ_0,0 .



As the ICEM is an alternative to NRQCD in hadroproduction, developments to extend it into other collision systems are still in progress. 
The authors of <cit.> anticipate that the value of λ_θ for production in ep collisions will also be very similar to the pp case, which they found to be compatible with the existing Tevatron and LHC data. 
In addition, they also find the free parameter F_𝒬 in photoproduction to be consistent with that in hadroproduction. The description of HERA H1 data <cit.> on photoproduction in the ICEM <cit.> is illustrated in Fig. <ref>.  However, the (I)CEM prediction introduces a parameter to keep the propagator at some minimum distance of M_ψ^2 from the pole. Thus, its prediction of the z-differential spectrum in photoproduction is likely to be complicated by large radiative corrections at 1-z≪ 1 if the parameter is removed, which was seen already in the LO analysis of Ref. <cit.> where the agreement with data at z→ 1 was reached only after introduction of an ad-hoc cut |t̂|>4m_c^2 on the partonic t̂ variable.

The observation that the CEM leads to unpolarised heavy-quarkonium hadroproduction at high- <cit.>, a result which is non-trivial to achieve with NRQCD fits, perhaps means that, in cases where CO LDMEs dominate, the dynamics of soft-gluon emissions should be taken into account more accurately than it is done in the fixed-order NRQCD factorisation. The recently proposed soft-gluon factorisation approach represents a progress in this direction <cit.>, whose phenomenological implications, however, remain to be investigated. 








 §.§ Legacy from HERA, the Tevatron and the LHC, and predictions for the EIC for cross-section and polarisation observables






 

  §.§.§ Status of  NRQCD LDME fits
 


 
 

  
A side note on the positivity of the LDMEs beyond LO. 
 Before discussing the NLO LDME fits, let us make a comment about the positivity of LDMEs. At LO in α_s, the LDMEs have a simple interpretation as “probabilities” of the transition of the QQ̅-pair in a certain colour, spin and angular-momentum state into an observed quarkonium. This physical interpretation follows from the operator definition of LDMEs (<ref>) in terms of “bare” fields <cit.> if QCD loop corrections are not taken into account and Wilson-line factors are ignored. Consequently, in LO calculations, LDMEs are typically assumed to be positive-definite. This is similar to the situation with LO (TMD)PDFs.

Already at NLO in α_s, both ultraviolet (UV)  and infrared (IR)

divergences appear in the operator definitions of LDMEs, see e.g. Appendix B of Ref. <cit.> as well as Section 6 of Ref. <cit.> and references therein. If NRQCD factorisation holds – which is yet to be proven beyond NNLO in α_s <cit.> – the IR divergences of the hard-scattering coefficients should cancel against the corresponding IR divergences of the LDMEs at all orders of the v^2 and α_s expansions, while the UV divergences appearing in LDMEs are removed by the operator renormalisation. The renormalised LDMEs then become non-perturbative fit parameters.  Therefore, these parameters do not necessarily have to be positive. Their definition involves the subtraction of the divergent part. In addition, the finite renormalised LDMEs are scheme- and scale-dependent, and mix with each 
other due to the NRQCD-scale evolution. The relation between the short-distance cross section and LDMEs, described above, is similar to the relation between NLO short-distance cross sections and QCD PDFs and/or fragmentation functions, which are also not necessarily positive-definite, at least if the calculation is truncated to a fixed order in α_s. This is the reason why there are usually no positivity constraints imposed in NLO LDME fits. One of the consequences of this is that the numerical values of LDMEs obtained in fits at NLO in α_s have limited physical significance outside the NLO context and should only cautiously be used in LO calculations, because this could create unjustifiable cancellation between some contributions.

In general though, it is not clear that negative NLO LDMEs would yield positive NLO cross sections for all possible measurable processes one could think of. Let us for instance mention the case of quarkonium-photon associated production for which it was shown <cit.> that some of the NLO LDME fits which we discuss below would yield negative NLO cross sections. Such a physical constraint on LDMEs at NLO has however not been systematically investigated as it requires the complete NLO computation of the hard scatterings for all the processes one wishes to consider.


 

  
Survey of existing NLO LDME fits. 

Several groups have performed fits of CO LDMEs for charmonia <cit.> and bottomonia <cit.> at NLO in α_s for the short-distance parts. We emphasise that the computation at NLO in α_s of short-distance cross sections for the production of NRQCD states (QQ̅[i]) is done in exactly the same framework of collinear factorisation by most of the groups with the exception of the fit of Bodwin et al. <cit.>. The latter computation includes, beside corrections at NLO in α_s,  the resummation of logarithms of P_T/m_ which become important at P_T≫ m_. Therefore the difference of the fits boils down mostly to the choice of different experimental data to fit and approximate (up to higher-orders in v^2) relations between different LDMEs which are assumed or not to hold exactly in the fitting procedure. For a detailed discussion, we refer to the recent review <cit.>. tab:LDME-fits_comp briefly compares phenomenological results of 
each fit for the case of charmonia using benchmark observables such as the 
cross sections and polarisation of inclusive prompt J/ψ produced in pp collisions as a function of  as well as  
photoproduction in ep collisions and the  total cross section of charmonium production in e^+e^- annihilation 
We also 
indicate in tab:LDME-fits_comp whether the corresponding set of LDMEs for J/ψ allows one to describe the prompt η_c hadroproduction 
-spectrum measured by LHCb <cit.> using 
heavy-quark-spin-symmetry relations between  η_c and J/ψ LDMEs which hold up to v^2 corrections.   
 



 
 The 
  spectra of the prompt inclusive 
 quarkonia produced in pp and pp̅ collisions at mid and large  at the Tevatron and the LHC are well described by all the fits mentioned in tab:LDME-fits_comp; this is 
 the major phenomenological success of NRQCD factorisation at NLO. Note, however,  that hadroproduction data with ≲ m_ (or integrated in  <cit.>) can not be simultaneously described by NLO NRQCD  fits of large  data. 
 In fact, most of the fits have been performed with even stronger  cuts, as indicated in tab:LDME-fits_comp. 
 
 The only existing global NLO LDME fit <cit.>, BK11, beyond hadroproduction, also provides a reasonable description of unpolarised charmonium production cross sections in e^+e^-, pp, pp̅ and ep collisions. The description of HERA H1 data <cit.> on J/ψ photoproduction by the BK11 fit is illustrated in the fig:LDME-fits_vs_H1-data_pT(a) and fig:LDME-fits_vs_H1-data_z(a).   However,  this fit is not able <cit.> to describe charmonium-polarisation observables, measured in hadroproduction at high-, see e.g. Ref. <cit.> for a global survey of heavy-quarkonium-polarisation data. 
 This situation is often referred to as the “heavy-quarkonium-polarisation puzzle” in the literature. Polarisation observables relevant for J/ψ production at the EIC will be discussed in Section <ref>. 



 
Two of the fits in tab:LDME-fits_comp, H14 and Z14, 
turned out to be able to simultaneously describe J/ψ and η_c hadroproduction data using heavy-quark-spin-symmetry 
 
relations between LDMEs. Remarkably, the J/ψ-polarisation observables in hadroproduction are also reasonably well reproduced by these fits but they significantly overestimate the HERA photoproduction cross section as can be seen in fig:LDME-fits_vs_H1-data_pT(b,c).  The same holds for all the other LDME fits (with the exception of BK11 discussed above), see fig:LDME-fits_vs_H1-data_pT(d-g). The discrepancies between the NRQCD NLO predictions with these fits range from 2 at ≃ 10 GeV up to 10 at ≃ 1-2 GeV in the case of pNRQCD and B14. This means that the yield predictions at the EIC using these LDMEs can be overestimated by up to one order of magnitude. Since the discrepancies remain at =10 GeV, which roughly corresponds to the maximum values which would be reached at the EIC, this should be kept in mind when considering predictions with CO contributions (except for the BK11 LDMEs) for the EIC case at any .
  




 As one can seen from fig:LDME-fits_vs_H1-data_z, all LDME fits except BK11 also strongly overestimate the z-differential cross section for z>0.6. The BK11 fit is consistent with the photoproduction data due to the cancellation between ^1S_0^[8] and ^3P_J^[8] channels. Other fits use this degree of freedom to accommodate the polarisation and/or η_c production data and therefore lose flexibility which is needed to achieve a global fit across different collision systems.

 

 In a recent study <cit.>, the NRQCD cross sections of J/ψ+Z and J/ψ+W hadroproduction have been completely calculated at NLO. Interestingly, the only set of LDMEs found to be marginally capable of reproducing the J/ψ+Z hadroproduction data from the LHC is the set of 
 Refs. <cit.>, referred to in the tab:LDME-fits_comp as “pNRQCD”. This fit uses potential-NRQCD relations between LDMEs to reduce the number of free parameters in the fit of the J/ψ  
 spectrum in hadroproduction and which also describes polarisation observables. However, this set of LDMEs is not able to describe 
 J/ψ photoproduction and e^+e^- annihilation data  and is  consistent with η_c hadroproduction data only within large uncertainties and with a  threshold for the  data large than for the η_c data. As just discussed, the pNRQCD fit, like all the hadroproduction, badly fails to account for the J/ψ-photoproduction data from the H1 collaboration at HERA as shown on fig:LDME-fits_vs_H1-data_pT(g) and fig:LDME-fits_vs_H1-data_z(g) which cast doubts on its relevant for EIC predictions.

Several fits of CO LDMEs for bottomonia have also been performed at  
NLO <cit.>. Only the most recent one <cit.> considered the Υ(1,2,3S) and χ_bJ(1,2P) LDMEs independently and systematically included the feed-down contributions from Υ(nS) and χ_bJ(nP) states with larger masses. These feed-down contributions constitute ∼ 40% of the Υ(1S, 2S) cross section, which is significant. In the case of Υ(3S), the feed down from χ_b(3P) states, which were discovered by ATLAS <cit.> and which lie just below the BB̅-threshold, also turns out to be significant (see <cit.> for a more detailed discussion of the feed-down impact). This was,  however, not taken into account in <cit.>. This may explain the difficulties of the corresponding fit to account for the Υ(3S) polarisation. The polarisation observables for Υ(1S,2S) states came out to be about consistent with data in this fit. We guide the reader to the recent Ref. <cit.> for a detailed discussion of the agreement with various polarisation observables. Note that there is no bottomonium data from inelastic photoproduction nor from e^+e^- annihilation. Hence, future measurements of Υ(nS) inclusive electro- and photoproduction at the EIC will serve as an excellent test of the LDME process-independence in the b-quark case, where it has more chances to hold due to smaller O(v^2) corrections.






  §.§.§ Recent developments regarding inclusive J/ψ photoproduction within the CSM




 
 

  
New P_T-enhanced contributions.
The recent study of Ref. <cit.>, performed within the CSM, is interesting regarding corrections which 
were not included in the NLO NRQCD analyses presented above, although they could become important at ≫ M_.  The study focused on the leading- leading-v next-to-leading-α_s corrections, within the NLO^⋆ approximation <cit.>. 

The latest HERA data from the H1 Collaboration <cit.> was first revisited, by including new contributions such as the pure QED one (γ + q →γ^⋆ + q→ + q  at 𝒪(α^3) where the off-shell photon γ^⋆ fluctuates into a ) 

and the associated + charm production (γ + g → + c +  c̅ and γ +{c,c̅}→+ {c,c̅}). The former involves quark PDFs in the initial state, while the latter is described within a LO Variable Flavour Number Scheme (LO-VFNS) <cit.>. It was shown that the CSM at 𝒪(αα_s^3) and 𝒪(α^3) is able to describe the latest HERA data at large . The NLO corrections to γ + g → + c +  c̅ were recently computed <cit.> and were found to increase the cross section a factor close to 2 in the HERA kinematics.

The corresponding predictions for the () spectrum in photoproduction at the EIC 
are shown in fig:EIC-CT14NLO-NLOstar 
with kinematical cuts on Q^2, the elasticity, z, and W_γ p≡√(s_γ p) inspired from the latest H1 measurements.  
The CT14NLO proton PDF set <cit.> was used. The factorisation and renormalisation scales were taken to be μ_F = μ_R = m_T = √(M_^2 + P_T^2), the transverse mass of the , later this is called m_T J/ψ and the corresponding uncertainties were evaluated by varying them in the interval μ_F, μ_R ∈ [1/2, 2]× m_T. The charm mass m_c was set to 1.5 GeV and the corresponding mass uncertainty was evaluated by varying it by ± 0.1 GeV. Moreover, the CS LDME ⟨𝒪^[^3S_1^[1]] ⟩ was taken to be 1.45 GeV^3. Finally, a 20% feed-down ψ^'→ was taken into account.



In fig:EIC-CT14NLO-NLOstar, predictions for two energy configurations are presented. At √(s_ep) = 45 GeV (fig:EIC-CT14NLO-NLOstar, left), as  increases, one enters the valence region. This makes the QED contribution become the dominant one at the largest measurable ≃ 11 GeV, with an integrated luminosity of ℒ = 100 fb^-1. Furthermore, γ + q fusion contributes more than 30% for > 8 GeV and the + unidentified charm contribution is comparable to the γ + g(q) fusion subprocesses. Hence, these so far overlooked contributions are going to be relevant at the EIC. At √(s_ep) = 140 GeV (right panel in fig:EIC-CT14NLO-NLOstar),  the yield is measurable up to ∼ 18 GeV. The QED contribution is the leading one at the largest reachable , while γ + g fusion is the dominant contribution up to ∼ 15 GeV. More generally, it turns out that the production of +2 hard partons (+{gg, qg, cc̅}) is dominant for ∼ 8 - 15 GeV. This could lead to the observation of + 2 jets with moderate , with the leading jet_1 recoiling on the +  jet_2 pair.
 



  
High-energy-enhanced contributions. 

In order to study the possible effects of higher-order QCD corrections enhanced by logarithms of the partonic centre-of-mass energy (ŝ), the Leading-Twist (LT) High-Energy Factorisation <cit.> (HEF) can be used. In many phenomenological studies, it is generalised to include, not only the resummation of ln (ŝ/M^2_)-enhanced effects in the leading-logarithmic approximation, but also the resummation of the “Sudakov” ln (M_/P_T) large logarithms at P_T≪M_ in the next-to-leading logarithmic approximation, assuming  CS state production, through the use of the Kimber-Martin-Ryskin-Watt (KMRW formula) <cit.>. However, the systematic study of the overlap between LT HE factorisation and the TMD factorisation 

usually employed to resum such transverse-momentum logarithms has been initiated only very recently <cit.>. The KMRW formula converts the set of usual collinear PDFs to the so-called unintegrated PDFs (uPDFs) of the LT HEF formalism. uPDFs depend not only on the longitudinal momentum fraction, x, but also on the transverse momentum of the parton. These objects can yield transverse momenta comparable to, or even larger than, M_ to the final state. This is indeed possible in the Regge limit ŝ≫M_T. For a more detailed review of the LT HEF and its connection to quarkonium physics, see 
Section 4.3 of 
Ref. <cit.>.

It has been shown earlier <cit.> that the phenomenological framework based on HEF with KMRW uPDF is capable of reproducing the J/ψ photoproduction data from HERA. This is already the case with the HEF coefficient function computed at LO in α_s and in the CS approximation of  NRQCD, as illustrated by the left panel of 
fig:HEF_H1-2010_EIC obtained with the version of KMRW uPDF introduced in the Ref. <cit.>. We note that the transverse-momentum integral of the uPDF exactly reproduces the input gluon PDF. The 
precise 
fulfilment of this normalisation condition both at x≪ 1 and x∼ 1 is important to avoid contradictions between LT HEF and NLO 
Collinear Factorisation (CF) 
predictions for the J/ψ prompt hadroproduction  spectrum in pp collisions at low energies, in particular at √(s_pp)=24 GeV for the planned Spin-Physics-Detector experiment at the NICA facility <cit.>.



From fig:HEF_H1-2010_EIC (left), one can see that there is still some room for additional contributions on top of the LO CS contribution from the fusion of a photon and a Reggeon γ(q)+R(x_1, q_T1)→ cc̅[^3S_1^[1]]+g
. These could be from cc̅[^3S_1^[1]]+c considered above (see fig:EIC-CT14NLO-NLOstar) and from CO contributions. The large scale uncertainty of the LO HEF prediction, shown in the fig:HEF_H1-2010_EIC, comes from the variation of μ_R and μ_F around their default value of M_TJ/ψ. Clearly, the uncertainty 
has to be reduced via the inclusion of the NLO corrections to make such predictions more precise.  

The comparison between LT HEF predictions for the EIC energy √(s_ep)=140 GeV 
and the full NLO CF 

CSM predictions (computed using FDC <cit.>) is shown in fig:HEF_H1-2010_EIC (right). The latter prediction is evaluated at 
a special value of the factorisation scale, μ_F=1.7 m_c, chosen <cit.> to minimise the NLO correction coming from the region of ŝ≫M_^2 (cfr. Section <ref>). There is 
a good agreement between these NLO predictions at the optimal scale and LO HEF predictions at the default scale μ_R=μ_F=M_T,J/ψ
: this indicates that the effects of the ln(ŝ/M_^2) 

resummation can be reproduced by the optimal factorisation scale choice at EIC energies and that the NLO CF prediction with the optimal scale is robust. At higher photon-nucleon collision energies, a matched calculation between LL HEF and NLO CF predictions, similar to that done in Ref. <cit.>, is necessary <cit.> to correctly capture the high-energy resummation effects at ŝ≫M_^2 while staying at NLO accuracy for ŝ∼M_^2.




  §.§.§ Testing NRQCD factorisation at the EIC

















  
Prompt J/ψ yields in photoproduction.

We plot in Figs. <ref> and <ref> the NLO NRQCD factorisation predictions for the 
- and z-differential photoproduction cross section of prompt J/ψ mesons in the EIC kinematic conditions. These predictions have been calculated using the short-distance cross sections of Refs. <cit.> and the LDME sets listed in tab:LDME-fits_comp. All LDME sets fitted only 
to 
the hadroproduction data predict a significantly (factor 3 to 6) higher J/ψ photoproduction cross section than the LDME set of 
Table 1 of 
Ref. <cit.>, denoted as “LDMEs Kniehl, Butenschön, fit # 1” 

in Figs. <ref> and <ref>, which includes the photoproduction data from HERA. We also plot the predictions performed with another set of LDMEs from the same paper, denoted as “LDMEs Kniehl, Butenschön, fit # 2”. The latter set of LDMEs had been fitted 
to 
the prompt J/ψ hadro- 
and photoproduction data corrected approximately for feed-down contributions from heavier charmonium states using constant feed-down fractions. For this fit, we calculate the feed-down contributions from χ_c0,1,2 and ψ(2S) decays to J/ψ using the χ_c LDMEs from Ref. <cit.> 
and the fit for ψ(2S) LDMEs performed in 
Ref. <cit.>. Calculating the feed-down contribution in this way is consistent with the treatment of feed-down in 
Ref. <cit.>.


 
 

 
 
 
 


 
 
 
 

 As expected, the predictions from both Kniehl-Butenschön LDMEs are reasonably close to each other. Yet, they differ from those obtained with the other LDME sets fit to hadroproduction data.


This is mostly because 
the latter sets predict a more pronounced z→ 1 growth of the cross section (see fig:Jpsi-photoprod-z_LDMEs) than the global fit LDME sets of Ref. <cit.> which, when integrated over z, translates into larger differential cross sections. This increase is due to both the ^1S_0^[8] and ^3P_J^[8] CO states.

It is important to note that such a rapid increase of the spectrum towards z→ 1  
is not a feature of the HERA data. Including these data in LDME fits calls for a 

compensation between contributions of ⟨ O^J/ψ[ ^1S_0^[8]] ⟩ and ⟨ O^J/ψ[ ^3P_J^[8]] ⟩ LDMEs 
resulting in different signs for these as in the LDME sets of Ref. <cit.>. Therefore the photoproduction data essentially fix the latter LDMEs and do not allow anymore to adjust them to describe the polarisation observables in hadroproduction, which leads to the polarisation puzzle discussed above. The EIC measurements will allow us to check the robustness of this feature of NRQCD predictions against variation of collision energy, since larger radiative corrections at z→ 1 could be expected at higher energies of the HERA collider.
 
The resolved-photon contribution manifests itself in the opposite region z≪ 1 (fig:Jpsi-photoprod-z_LDMEs) and EIC data are less sensitive to it than HERA data, again due to lower collision energies. Therefore the cleaner test of process-independence of LDMEs can be performed with EIC photoproduction data rather than with HERA data.  
 






  
Prompt  yields in Q^2-integrated lepton-nucleon interactions.














Another possibility to study the contributions of various LDMEs is to consider single-inclusive production of J/ψ in ep collisions, without detecting the final-state electron, as was pioneered recently in  Ref. <cit.>: e(ℓ)+h(p)→ J/ψ(P)+X. The rapidity (y) and transverse momentum () distributions of inclusive 
production at the EIC are promising observables for both studying the production mechanism of heavy quarkonia and extracting 
PDFs, in particular, the gluon PDF, complementary to other observables described in Section <ref>. 
 
When the transverse momentum of  defined relatively to the lepton-hadron collision axis 
≫ m_c 
the perturbative hard coefficient functions for producing the pair receive large higher-order QCD corrections that are enhanced by powers of ln(P_T^2/m_c^2).  Such logarithmically-enhanced higher-order corrections can be systematically resummed and factorised into 
FFs <cit.>.  
On the other hand, when ≳ m_c, the perturbative hard coefficients at a fixed order in α_s should be sufficient.

In addition, the occurrence of a hard partonic collision producing the with large transverse momentum ≫ m_e
necessarily induces multiple photon emissions from the incoming lepton, leading to large higher-order QED corrections enhanced by powers of ln(^2/m_e^2). As we discussed in Section <ref>, these QED corrections can also be systematically factorised and resummed into universal 
LDFs <cit.>.  
In order to predict the production rate of at the EIC, a new factorisation formalism, which takes into account both collision-induced QCD and QED radiation and provides a systematic transition from ≳ m_c to ≫ m_c, was introduced <cit.>. 
The 
factorisation formula for the inclusive 
production cross section is given by:

    E_J/ψdσ_eh→ J/ψ(P_)X/d^3 P_   =   ∑_a,b ∫ dx_a  f_a/e(x_a,μ_F^2)
    ∫ dx_b  f_b/h(x_b,μ_F^2)
    
          ×[ 
    E_J/ψdσ̃^ Resum_ab→ J/ψ(P_)X/d^3 P_
    +
    E_J/ψdσ̃^ NRQCD_ab→ J/ψ(P_)X/d^3 P_
    - E_J/ψdσ̃^ Asym_ab→ J/ψ(P_)X/d^3 P_]  ,

where indices a,b, in principle, run, respectively,  over all lepton and parton flavors, but in practice, as an approximation, a takes into account only (e,γ,e̅). The functions f_a/e(x_a,μ_F^2) and f_b/h(x_b,μ_F^2) are the LDFs of an electron and the usual parton PDFs 
respectively, depending on  partonic momentum fractions, x_a and x_b. The LDFs satisfy the DGLAP-like μ_F-evolution  equations mixing the QED and QCD splittings <cit.>.  In Eq. (<ref>), the  partonic cross sections σ̃_ab→ J/ψ(P_)X are computed with all the perturbative collinear singularities along the direction of  colliding lepton (a) and parton (b) removed. These singularities are absorbed into f_a/e and f_b/h, respectively.  


The cross section dσ̃^ Resum in eq:jpsi-lp-fac represents the partonic cross section with the ln(P_T^2/m_c^2) contributions being resummed to describe the production rate for ≫ m_c, as we have mentioned above. In 
σ̃^ NRQCD, the production of cc̅[^2S+1L^[1,8]_J]-state at the perturbative stage is computed at fixed order in α_s and the corresponding non-perturbative formation of a from a produced pair is taken care 
using the NRQCD velocity expansion and universal NRQCD LDMEs. This part of the cross section  should provide a good description of the production rate when ∼ m_c.  Finally, σ̃^ Asym is equal to a fixed-order expansion of σ̃^ Resum to the same order in α_s as in σ̃^ NRQCD. The latter part is needed to remove the double counting between σ̃^ Resum and σ̃^ NRQCD. 
By including all these three terms, this factorisation formalism 
can be applied to both lepton-hadron and hadron-hadron collisions, as well as e^+e^- collisions <cit.>, providing a smooth transition when observed ∼ m_c increases to ≫ m_c.  

The predictive power of eq:jpsi-lp-fac relies on the factorisation of each term and our ability to calculate them. Up to next-to-leading power corrections in m_c/, the σ̃^ Resum can be factorised as <cit.>, 

    E_J/ψdσ̃^ Resum_ab→ J/ψ(P_)X/d^3 P_   ≈   ∑_k∫dz/z^2 D_k→ J/ψ(z,μ_F^2) 
    	E_kdσ̂_ab→ k(p_k)X/d^3 p_k(z,p_k=P_/z,μ_F^2)
    	
           +
    	∑_κ∫dz/z^2D_[(κ)]→ J/ψ(z,μ_F^2) 
    	E_kdσ̂_ab→ [(κ)](p_k)X/d^3 p_c(z,p_k=P_/z,μ_F^2)
    	  ,

where k=q,g,q̅ and κ=v,a,t 
for   pairs respectively in a vector, axial-vector or tensor spin state <cit.>. The first and second terms are the factorised leading power (LP) and next-to-leading power (NLP) contributions to the cross section in its 1/ expansion. The corrections to eq:jpsi-lp-resum are suppressed by 1/^4 and cannot be further factorised <cit.>.  The universal single-parton and double-parton () FFs, D_c→ J/ψ(z,μ_F^2) and D_[(κ)]→ J/ψ(z,μ_F^2), respectively, satisfy a closed set of evolution equations with respect to 
changes of the factorisation scale μ_F <cit.>.  Solving these evolution equations one resums the logarithmic contributions scaling like ln(P_T^2/m_c^2) to these FFs. The universal FFs at an input scale μ_F = μ_0≃ 2m_c can be calculated assuming NRQCD factorisation <cit.> in terms of universal NRQCD LDMEs,

    D_c→ J/ψ(z,μ_0^2) 
    	   ≈   ∑_[^2S+1L_J]d̂_c→[^2S+1L_J](z,μ_0^2) ⟨ O_[^2S+1L_J]^(0)⟩  ,
    	
    D_[(κ)]→ J/ψ(z,μ_0^2) 
    	   ≈   ∫_-1^1 du ∫_-1^1 dv  D_[(κ)]→ J/ψ(z,u,v,μ_0^2) 
    	
       ≈   ∑_[^2S+1L_J]d̂_[(κ)]→[^2S+1L_J](z,μ_0^2) ⟨ O_[^2S+1L_J]^(0)⟩ .


eq:jpsi-ffs-cc involves further approximations, neglecting possible differences between the momentum fractions carried by the pair in the amplitude, u, and and its complex-conjugate, v, which can be taken into account through the more general FF D_[(κ)]→ J/ψ, defined in <cit.>. The approximation in the second line of eq:jpsi-ffs-cc reflects the fact that the 
integral of this function is dominated by the vicinity of u=v=1/2 <cit.>. 


 










The formalism described above has been already tested partially in the case of pp collisions, where instead of LDFs in eq:jpsi-lp-fac one substitutes  the proton PDFs. With perturbatively calculated short-distance matching coefficients for both single-parton and -pair FFs at the input scale <cit.> and solving the coupled evolution equations for these FFs, the factorised and resummed cross section in eq:jpsi-lp-resum 
describes 
the  distribution of production at the LHC and Tevatron <cit.> for  > 10 GeV, as we note in 
tab:LDME-fits_comp. 
At the LHC energies, the LP contributions, namely the first term in eq:jpsi-lp-resum, dominate when ≫ 20 GeV, while the NLP contributions, namely the second term in eq:jpsi-lp-resum, are comparable at ∼ 20 GeV and become dominant when  further decreases, which is critically important to describe the shape of the observed P_T distribution.  





Making predictions of the distribution of inclusive production at the EIC requires the knowledge of the universal LDFs.  In fig:LDFs of <ref>, the scale dependence of the LDFs with and without the mixing of QED and QCD evolution is shown.  

Like in any factorisation approach, the perturbatively calculated short-distance partonic cross section, such as σ̂_ab→ k(p_k)X in eq:jpsi-lp-resum, does not depend on the details of the hadronic state produced. It has been calculated for single hadron production at LO <cit.>, at NLO <cit.>, and at NNLO <cit.>. The fixed-order calculation for σ̃^ NRQCD has been carried out in NRQCD up to NLO <cit.>.  

In fig:jpsi-pt-eic, we present the predictions of the distribution of inclusive production in ep collisions at the EIC for √(s_ep)=140 GeV. For these predictions, only the σ̃^ Resum term in eq:jpsi-lp-fac is used and the same LDMEs that we used for describing the production at the LHC and Tevatron energies <cit.> are taken here. These LDMEs are close to those from the  Chao et al. <cit.> fit (H14) mentioned in tab:LDME-fits_comp. 
The CT18ANLO PDF central set <cit.> was used for the proton PDFs. Unlike production at the LHC and the Tevatron, the reach in the defined with respect to the lepton-hadron axis, is much smaller due to the smaller collision energy. The solid line in fig:jpsi-pt-eic refers to the total contribution, which is dominated by the subprocess  γ + g → [cc̅] + g (NLP Photon) and e + g → [cc̅] + e (NLP Lepton)  with the pair fragmenting into J/ψ.  The lepton or photon initiated LP contribution to the production cross section, namely the first term in eq:jpsi-lp-resum, is dominated by the lowest-order subprocesses, such as e+q→ e+q or γ+q→ g+q, respectively, with a produced parton fragmenting into the observed J/ψ, and is strongly suppressed by the single-parton FFs at the EIC energy. In summary, the LP contributions are essentially irrelevant in the EIC kinematics. Therefore, a matching to the fixed-order calculations (described above in this section), including the second and third terms in eq:jpsi-lp-fac, is awaited for.



  
Polarisation of J/ψ in photoproduction.



Since the prediction <cit.> in 1994 of a transversely-polarised  hadroproduction yield at high , much hope has been put in polarisation measurements to advance our understanding of quarkonium production, with a very limited success though <cit.>. NLO CSM computations of polarisation observables in photoproduction were performed in 2009 <cit.> and subsequently completed with the COM NRQCD contributions in 2011 <cit.> without clear conclusions owing to the large uncertainties in the H1 <cit.> and ZEUS <cit.> data and in the theory. 

In  Figs <ref> and  <ref>, we show the NLO NRQCD predictions for the  and z dependence of the polarisation parameter λ_θ of promptly photoproduced J/ψ mesons in the EIC kinematic conditions. These predictions include CS and CO contributions using the LDME sets discussed in  Section <ref> as well as direct and resolved-photon interaction contributions (see the Figs. <ref> and <ref> for the corresponding differential cross-section plots). As one can see from 
fig:Jpsi_photoprod_lambda-pT, the 
-dependent NRQCD predictions for all LDME sets are roughly consistent with unpolarised production (λ_θ=0 in all frames), unlike the predictions of the CSM, which leads to significant polarisation of photo 
produced J/ψ mesons. In the z-dependent case, the region of z→ 1 has the most 
discriminating power between different LDME sets. We however have reasons to doubt the relevance of these predictions given that all but the BK11 LDMEs are unable to describe the corresponding HERA data. From Figs. <ref> and  <ref>, one also observes that the detailed behaviour of λ_θ for different LDME sets is significantly different for different polarisation frames, which could be an important tool for additionally constraining the theory.








  
Polarisation of J/ψ in electroproduction.














The HERA collider experiments provided some results on the J/ψ polarisation, mostly for photoproduction <cit.>,  
but unfortunately these data do not allow to favour or disfavour different models and/or approaches.
The reasons behind this are twofold: data were not precise enough and they were collected in regions where theoretical predictions are very close to each other <cit.>. Furthermore in Ref. <cit.>, Yuan and Chao showed that the estimates for the λ_θ parameter in SIDIS, within both the CSM and NRQCD approaches, are overlapping for most of the values of the variable z.
In this respect EIC could play a crucial role: highly precise data are expected and other/extended kinematical regions could be explored. 



In the following, we present some predictions at LO, both in the CSM and NRQCD frameworks, adopting different NLO LDME sets. 
Some comments are therefore in order: (i) as previously discussed, the combined usage of NLO hard scattering with NLO LDMEs is subject to great caution. As of now, only the CSM part of the electroproduction cross section has been computed at NLO <cit.>. The only full NRQCD analysis has been performed at LO <cit.> and show mixed agreements between the different NRQCD predictions and HERA data; (ii) a number of quarkonium-production processes exhibit very large QCD corrections to polarisation observables <cit.>. The following LO results should therefore only be considered as a simple guidance for future measurements and certainly not as quantitative predictions to which future measurements should be confronted to. In this context, a NLO NRQCD analysis of electroproduction is eagerly awaited for.






fig: Jpsi lambda pol vs PT SIDIS shows some estimates for the λ_θ parameter at the centre-of-mass 
energy √(s_ep) = 45 GeV 
together with their uncertainty bands, visible mostly for the CSM and obtained by varying the factorisation scale in the range μ_0/2 < μ_F <2 μ_0, with μ_0 = √(M_^2 + Q^2). The integration regions are detailed in the legend box. 
No uncertainty bands from LDMEs are included, instead predictions for different sets are presented: C12 <cit.>, BK11 <cit.> and G13 <cit.>. This illustrates their impact on the results. From fig: Jpsi lambda pol vs PT SIDIS, it is clear that the λ_θ value can be significantly different if we consider different frames. In particular, the Gottfried-Jackson frame provides the better overall separation between CSM and NRQCD curves.





Another possibility offered by the EIC experiment is the collection of data at different energies. In fig: Jpsi lambda pol vs energy SIDIS, the impact coming from the energy variation on CSM and NRQCD predictions is shown. In this case, only the central values are presented (μ_F = μ_0); for the lower energy, √(s_ep) = 45 GeV, the integration region is the same as in fig: Jpsi lambda pol vs PT SIDIS, while for √(s_ep) = 140 GeV a wider W integration is considered (see legend).
Even focusing on one specific frame, like the helicity frame in fig: Jpsi lambda pol vs energy SIDIS, one clearly sees that the CSM is more affected by the energy shift. Note that moving to higher energies allows one to access contributions with higher virtuality, with an interesting effect: in the CSM these contributions are opposite to the lower virtuality ones (reducing the size of the estimates), while in NRQCD this phenomenon is less important. It however remains to be shown that such discriminant effects remain at NLO.






  
Prompt η_c and χ_c yields.


As mentioned in 
Section <ref>, the dominance of 
the CS mechanism in prompt-η_c hadroproduction at ≳ M_η_c was not expected by NRQCD factorisation. Therefore, from the point of view of studies of the heavy-quarkonium production mechanism, it is important to understand if this feature of η_c production persists also in ep collisions. If it is indeed the case, then η_c hadro-, photo- and leptoproduction can be used as 
a tool for hadron-structure studies with 
a reduced uncertainty stemming from the CO mechanism compared to production of other charmonium states.






In recent works <cit.>, η_c photo- and electroproduction cross sections were computed including all the CO and CS contributions  at LO in α_s. In the case of photoproduction <cit.>, both direct-photon and resolved-photon interactions were taken into account. The CS contribution had been assumed to be negligible in earlier studies <cit.>, because the corresponding direct-photon interaction subprocess appears at O(αα_s^3) due to the necessity of 
two-gluon radiation in the final state to produce a cc̅[^1S_0^[1]] pair and because resolved-photon contributions were assumed to be small. However, it was found <cit.> that the resolved-photon subprocesses make the CS contribution to the photoproduction cross section non-negligible. These predictions, updated with the use of CT14LO PDFs, are shown in Figs. <ref> and <ref>.  The CO contributions were computed by converting the J/ψ CO LDME sets listed in tab:LDME-fits_comp to the 
η_c LDMEs through HQSS relations valid up v^2 corrections. 
As one can see from these figures, the CO contributions are still important and the cross section at z>0.5 strongly depends on the LDME choice.




For 
electroproduction
 <cit.>, the CS contribution is also sizeable, but for  a different reason, namely an additional Q^2-dependent terms appearing in the short-distance cross section. Of course, the main problem of the predictions for η_c production in 
ep collisions is that they so far have been done only at 
LO in α_s. The NLO corrections could be particularly important for the CS ^1S_0^[1] state whose LO contribution is highly suppressed at ≳ M_η_c 
in photo- and leptoproduction in comparison to CO states, especially ^3S_1^[8]. As known from  J/ψ production, this suppression will be lifted by large NLO corrections <cit.>. NLO calculation, at least in the CS channel, should be done before drawing conclusions about the importance of the CS mechanism in η_c production at the EIC.






Besides J/ψ and η_c production, it is also essential to study χ_c0,1,2 states at the EIC.
Photo- or electroproduction of these mesons has not been observed experimentally yet.  NLO NRQCD predictions for the photoproduction cross sections of χ_c0,1,2 radiatively decaying to J/ψ are shown in fig:chic-photoprod. They are based on known calculations of short-distance cross sections for  photoproduction <cit.> and the χ_c0 LDME values obtained in hadroproduction fits by Ma et al. <cit.>, respectively Bodwin et al. <cit.>. We remark that the former χ_c0 LDME values are also those used by Gong et al. in Ref. <cit.> and in the LDME set denoted “Kniehl, Butenschoen, fit #2” in Figs. <ref> and <ref>. We remind the reader that for P-wave production at NLO in NRQCD, one cannot make a clear distinction between CO and CS contributions as they directly depend on the NRQCD factorisation scale, μ_Λ.

It is an expected feature that resolved photon contributions dominate photoproduction at low z. Interestingly, however, the predictions of fig:chic-photoprod are dominated by the resolved-photon contribution already for z below 0.5. Moreover, it is only due to the resolved photons that the χ_c cross sections are positive at low z after all. This feature of the theoretical predictions may indicate our poor understanding of χ_c 
photoproduction, but if confirmed, the photoproduction of these mesons could serve as a useful source of information about the poorly known gluon component of photon PDFs.










 §.§ Learning about quarkonia from TMD observables



  §.§.§ LDME constraints from TMD observables






One important reason to investigate quarkonium production at the EIC is the possibility to probe TMDs that have not been extracted from experiments yet. The semi-inclusive heavy vector quarkonium production process, e   p → e'  J/ψ  (Υ)   X at small transverse momentum, , is expected to offer a promising probe of gluon TMDs[Due to the presence of the large scale given by the quarkonium mass M_≈ 2m_Q, one can consider not only electroproduction, but in principle also the photoproduction case (Q^2≈ 0). A large photon virtuality is expected to suppress background from diffraction and higher-twist effects <cit.>. To our knowledge, at present there are no studies of the numerical impact of such background on the photoproduction process γ  p → J/ψ  (Υ)   X in the TMD regime. 
], as will be discussed extensively in section <ref>. Besides gluon TMD extractions, this process may also allow for improved determinations of certain LDMEs. In this way EIC can also improve our knowledge on NRQCD. 


At small P_T,

the differential cross section is expected to be described in terms of 
TMDs. 

As will be discussed in detail in the next subsection, for quarkonium production, this involves

TMD shape functions <cit.>, rather than TMD FFs  like for light hadron production. At the lowest order, α^2 α_s, the process e   p → e'  J/ψ  (Υ)   X at small transverse momentum is described by photon-gluon scattering producing a heavy quark-antiquark pair in the 
CO state. The transition from the heavy-quark pair into the bound state is then described by a shape function. If one assumes the shape function to be a delta function in transverse momentum, one can connect to the standard NRQCD expressions for this transition. 
To lowest order in the strong coupling, 

but with the inclusion of the NNLO in v^2 
^1S_0 and ^3P_J (J=0,1,2) CO intermediate states <cit.>, the resulting expression for the cross section involves 
two of the CO LDMEs which were discussed above,
⟨ O[^1S^[8]_0]⟩ and ⟨ O[^3P^[8]_0]⟩, for which constraints from new types of observables are clearly welcome. In this way, measurements of the transverse-momentum spectrum of e   p → e'  J/ψ (Υ)   X in the TMD regime can lead to improved determinations of these CO LDMEs.  However, inclusion of higher-order corrections, in particular from the leading v CS NRQCD contributions at α^2 α_s, and the proper shape functions will be required for a robust extraction of these LDMEs.







  §.§.§ TMD effects from quarkonia: shape functions


















The NRQCD factorisation approach can only be applied for transverse-momentum spectra when the quarkonium state is produced with a relatively large transverse momentum compared to its mass,  
≳ 2 m_Q. 

This is because the emissions of soft gluons from the heavy-quark pair cannot modify the large transverse momentum of the bound state. The large 
 is generated in the hard process through recoil off unobserved particles, while the infrared divergences are parametrised in terms of the well-known 
LDMEs, 
collinear PDFs and 
FFs, depending on the particular process under consideration.

On the contrary, when the quarkonium is produced with a small transverse momentum, all soft gluon effects can no longer be factorised in terms of standard TMD PDFs. In order to properly deal with soft-gluon radiation at small 
 in a transverse-momentum spectrum of quarkonium, it has recently been found that one needs to promote the LDMEs to so-called TMD shape functions (TMD ShFs) <cit.>. 
Earlier, similar shape functions had been introduced in quarkonium photo-/leptoproduction in the endpoint region <cit.>, which however are functions of z, but a more general form was discussed in <cit.>. 


The newly introduced non-perturbative TMD ShFs encode the two soft mechanisms present in the process at low 
: the formation of the bound state and the radiation of soft gluons. As a consequence, they parametrise the transverse-momentum smearing of the bound state, and carry a dependence on the factorisation and rapidity scale.

Schematically, for the production of a single quarkonium state  at the EIC, with mass M_,  we have:

    dσ   ∼
    F_g/P(b_T;μ,ζ) ∑_i∈{^1S_0^[1],…} 
    H^[i](M_, Q;μ) Δ^[i](b_T,μ,ζ)
     ,

where F_g/P stands for any of the eight leading-twist 
gluon TMDs <cit.>, H^[i]
are the process-dependent hard scattering coefficients and Δ^[i] are the 

quarkonium TMD ShFs <cit.>. 
The above formula is written down in coordinate space where 
b_T is Fourier-conjugate to the quarkonium transverse momentum P_T^* (to be specific, in the virtual photon-proton centre of mass frame). 

Moreover, μ and ζ are the factorisation/resummation and rapidity scales, respectively.
The summation is performed over the various colour and angular-momentum configurations (i) of the QQ̅ pair.
Similarly to LDMEs, the TMD ShFs are of a specific order in the relative velocity v of the heavy quark-antiquark pair in the quarkonium rest frame. 
Therefore, the factorisation formula is a simultaneous expansion in v and λ = ^*/M_.
The operator definition of a bare[It is understood that the TMD ShF in the factorised cross section in eq:TMD-shape is free from rapidity divergences, it has been divided by the relevant soft factor which has also been used to properly subtract rapidity divergences in the gluon TMD F_g/P.] 
TMD ShF with NRQCD quantum numbers i is:

    Δ^[i](b_T,μ,ζ) ∝ ∑_X_s 0( O_i^† Y_n^†)^ab(b_T)X_s   X_s( Y_n  O_i)^ba(0)0  ,
 
which is just the TMD generalisation of the LDME operator definition in eq:LDME_def. On the r.h.s., the usual LDME operators 𝒪 are evaluated at positions b_T and 0 and sandwiched between the vacuum |0⟩ and the state |  X_s⟩ of the produced quarkonium together with possible soft radiation carrying away color. 

Moreover, these operators are multiplied by 
Wilson lines 𝒴_n parametrising 
the resummation of gluons exchanged between the hard part and  
the 
state |  X_s⟩. 


The operator definition in eq:ShFoperator can be related to the NRQCD LDMEs by the first term in an operator product expansion (OPE) for b_T→ 0:


    Δ^[i](b_T,μ,ζ) 
    =
    ∑_n
    C^[i]_n(b;μ,ζ)
    ×⟨ O^ [n]  ⟩(μ)
    +𝒪(b_T)
     .



In order to extend this expression to larger b_T, one can introduce a prescription like b_T →  b_T^*≡ b_T/√(1+(b_T/b_T,max)^2)≤ b_T,max to ensure validity of this perturbative expression 
and include a nonperturbative overall factor Δ^[i]NP:

    Δ^[i](b_T,μ,ζ) 
     ≡Δ^[i]NP(b_T) ∑_n
    C^[i]_n(b_T^*;μ,ζ)
    ×⟨ O^ [n]  ⟩(μ)
     .





This expression involves the usual “collinear” LDMEs, multiplied by perturbatively calculable Wilson coefficients C_n^[i](b_T;μ,ζ) to match the expansion on pQCD, and a non-perturbative part Δ^[i]NP 

that needs to be modelled or extracted from experimental data. Note that, in principle, at higher orders in α_s, there might be operator mixing: the ^1S_0^[8] TMD ShF could become dependent on the ^3P_0^[8] LDME, hence the sum over NRQCD states n in eq:ShFOPE.

In Ref. <cit.>, the OPE  of eq:ShFOPE is implemented in a practical way by studying single-inclusive J/ψ electroproduction. In the regime 
^*2∼ Q^2∼ M_^2, 


with 
^* being the transverse momentum of the quarkonium in the virtual photon-proton centre-of-mass frame

and μ either given by Q or by the 
quarkonium mass M_, the cross section is computed as usual in collinear factorisation. On the other hand, when 
^*2≪μ^2, TMD factorisation eq:TMD-shape applies. By comparing both cross sections in the kinematical regime Λ^2_QCD≪^*2≪μ^2, one can 
match the relevant TMD ShF onto the collinear LDMEs, 

confirming the need for introducing shape functions. 
The analysis of Ref. <cit.> was revised in Ref. <cit.>, modifying the obtained expression for the shape function, but not its necessity. 


To summarise, the factorisation theorem in eq:TMD-shape contains a convolution of two non-perturbative hadronic quantities at low transverse momenta: the gluon TMD PDFs and the TMD ShFs. It is therefore possible to perform a phenomenological extraction of gluon TMDs from quarkonium production processes.
However, to do so, one also needs to model or extract the involved TMD ShFs. This is analogous to SIDIS where one observes a light hadron, where one needs information on the light-hadron TMD FFs in order to extract quark TMD PDFs.






























































































  §.§.§ Azimuthal cos 2 ϕ_T^* modulation in electroproduction









In (semi-inclusive) quarkonium electroproduction on an unpolarised proton target, 
an azimuthal cos 2 ϕ_T^* modulation (see Section <ref> for our kinematic definitions) of the differential cross section will arise from 
linearly 
polarised gluons inside 
the unpolarised proton. These are described by the TMD h_1^⊥  g <cit.>[Note that, in the photoproduction regime, one cannot determine the angle ϕ_T^* because the lepton plane is not defined, hence, also not the cos 2 ϕ_T^* modulation. In photoproduction, azimuthal modulations can only be seen for two-particle observables.]. 
In many studies the shape functions of the quarkonium are assumed  are not considered and then the differential cross section can be written as: 

    dσ=   1/2sd^3l'/(2π)^32E_l'd^3P_/(2π)^32E_P_∫ dx d^2𝐤_⊥(2π)^4δ(q+k-P_)
       ×1/Q^4ℒ^μμ'(l,q)Φ^νν'
      (x,𝐤_⊥) ℳ_μν(ℳ_μ'ν')^∗,

where ℳ_μν is the amplitude of  
production of the quarkonium  in  the subprocess 
γ^∗+g→

, ℒ^μμ' is the leptonic tensor, and the 
gluon correlator is given by <cit.>:

    Φ^νν'(x,𝐤_⊥)=-1/2x{g_⊥^νν'f_1^g(x,𝐤_⊥^2)-(k_⊥^νk_⊥^ν'/M_p^2+g_⊥^νν'𝐤_⊥^2/2M_p^2)h_1^⊥ g(x,𝐤_⊥^2)}.

Here, g_⊥^νν'=g^νν'-P^νn^ν'/P· n-P^ν'n^ν/P· n, x and 𝐤_⊥ are the light-cone momentum fraction and transverse momentum of the gluon. 

The asymmetry is defined as:

    ⟨cos(2ϕ_T^*)⟩   =   ∫ dϕ_T^* cos(2ϕ_T^*)dσ/∫ dϕ_T^* dσ,

where ϕ_T^* is the azimuthal angle of the production plane of J/ψ with respect to the lepton scattering plane. 



In a more complete picture of the P_T^*2≪ M_J/ψ^2 ∼ Q^2 
region, as explained in <ref>, the TMD factorisation applies and LDMEs are promoted to TMD ShFs. Hence, the differential cross section for this process can be recast in the following form: 


    dσ^UP/d y  dx_B  d^2 P^*_T =    N [∑_n A_UP^[n]  C[f_1^g  Δ^[n]] + ∑_n B_UP^[n]  C[w h_1^⊥  g Δ_h^[n]]  cos 2 ϕ_T^* ] ,





where the subscript UP on the amplitudes A_UP^[n] and B_UP^[n] denotes the polarisation state of the proton (U, since it is unpolarised) and of the quarkonium (P=U,L,T), respectively, and N denotes an overall normalisation factor. Here, the quarkonium polarisation is defined with respect to the direction of the quarkonium three-momentum in the virtual photon - proton 
centre-of-mass frame. 
Measurements of the transverse-momentum dependence of the above cross section at the EIC would 
allow one to gather information on the so-far unknown quarkonium shape functions. In particular, the cos2ϕ_T^*-weighted cross section would give access to a linear combination of the convolutions C[w h_1^⊥ g Δ^[n]_h], with n=^1S_0^[8] or n=^3P_0^[8]. 
Here the weight in the convolution expression 

    C[w h_1^⊥ g Δ^[n]_h]  (q_T) ≡∫ d^2p_T∫ d^2k_T δ^2(p_T+k_T-q_T)   w(p_T,q_T)   h_1^⊥ g(x,p_T^2)  Δ^[n]_h(k_T^2)  ,

is given by (in standard TMD notation, note however that q_T will correspond to P^*_T used here)


    w(p_T,q_T) = 1/M_p^2q_T^2[2(p_T·q_T)^2 - p_T^2q_T^2]   .



On the other hand, integrating over ϕ_T^* would single out a combination of the convolutions C[f_1^g Δ^[n]](P_T^*) for the same octet S- and P-waves, which could be in principle disentangled by looking at different values of the inelasticity y. Measurements of these observables should help to establish the relevance of smearing effects and, in case they turn out to be sizeable, to even perform a first extraction of the shape functions. In this way, it would be possible to compare Δ^[n] and Δ^[n]_h as well as determine some other properties, like their relations with the LDMEs and their dependence on n.

For unpolarised quarkonium production (P=U), 

applying the above expressions gives the following normalised asymmetry ratio:

    ⟨cos 2ϕ_T^* ⟩≡∫ d ϕ_T^* cos 2ϕ_T^* dσ^UU/d y  d x_B  d^2 P^*_T/∫ d ϕ_T^*  dσ^UU/d y  dx_B  d^2 P^*_T = 1/2∑_n B_UU^[n]  C[w h_1^⊥  g Δ_h^[n]]/∑_n A_UU^[n]  C[f_1^g  Δ^[n]].





As the matching analysis mentioned in Section <ref> suggests, it is expected that the shape functions are proportional to the LDMEs belonging to the [n] state, 
at least at LO: Δ^[n] ( k_^2;μ^2) ≃⟨ O^ [n]  ⟩ Δ( k_^2;μ^2) and Δ_h^[n] ( k_^2;μ^2) ≃⟨ O^ [n]  ⟩ Δ_h( k_^2;μ^2), for some 
Δ( k_^2;μ^2) and Δ_h( k_^2;μ^2). 

In this case the above asymmetry expression reduces to: 

    ⟨cos 2ϕ_T^* ⟩ = 1/2B_UU/A_UU  C[w h_1^⊥  g Δ_h ]/ C[f_1^g  Δ],

where A_UU = ∑_n A_UU^[n] ⟨ O^ [n]  ⟩ and 
B_UU= ∑_n B_UU^[n] ⟨ O^ [n] ⟩. 
At LO, the coefficients appearing in this expression are <cit.> :









    A_UU^[^1S_0^[8]] = 1+y̅^2,  
    A_UU^[^3P_0^[8]] = [2 y̅7+3Q̂^2/1+Q̂^2 + y^2 7+2Q̂^2+3Q̂^4/(1+Q̂^2)^2] 1/m_Q^2,  
    B_UU^[^1S_0^[8]] = - y̅,   
    B_UU^[^3P_0^[8]] =  3-Q̂^2/1+Q̂^2y̅/m_Q^2.

Here, we defined y̅=1-y, with y being the inelasticity variable (see eq:SIDIS:kin-var), and Q̂^2≡ Q^2/(4m_Q^2) and we approximated m_≃ 2 m_Q, where m_Q denotes the heavy-quark mass. 

At the EIC, one 
could try to determine the LDMEs together with the gluon TMDs. The Q^2 and y dependence of the P_T^*-independent pre-factor B_UU/A_UU can then be exploited
, as it makes the observable dependent on different linear combinations of the LDMEs. This can be paralleled to the slight rapidity dependence of the LDME linear combination appearing in the polarisation of the hadroproduction yield <cit.>. 
Another option is to consider ratios in which the gluon TMDs cancel out <cit.>, although that may only hold 
at LO in certain cases. An example of this will be discussed in Section <ref> where the quarkonium polarisation is used to cancel out the gluon TMDs.  

A further constraint on the LDMEs comes from the bound on the above asymmetry. At leading order, the bound q_T^2 | h_1^⊥  g(x,  q_T^2) | / (2M_p^2) ≤ f_1^g(x,  q_T^2) <cit.> and the fact that |⟨cos 2ϕ_T^* ⟩ | ≤ 1 leads to the condition |B_UU/A_UU| ≤ 1. The LDMEs that determine the ratio B_UU/A_UU will have to respect this bound. In this way one can find for instance that the CO LDMEs from Ref. <cit.> do not respect this bound at LO (and A_UU which should be positive becomes negative below the central value within the 1σ uncertainty range), but it has to be noted that these LDMEs were obtained at NLO from hadro- and photoproduction data. 


The ratio B_UU/A_UU at LO is shown in fig:ratioBUUoverAUUandratioAULoverAUU (left plot) for two different CO LDME sets: one obtained at LO <cit.> (SV), which is very similar 
to the NLO fit <cit.>, denoted C12 in tab:LDME-fits_comp, and another obtained at NLO with FF <cit.> (BCKL), denoted B14  in tab:LDME-fits_comp. The uncertainty bands are obtained assuming uncorrelated uncertainties on the ⟨ O^[^1S_0^[8]]⟩ and ⟨ O_8^[^3P^[8]_0]⟩  determinations. The ratio is shown as a function of Q̂^2. Here Q̂^2=0.01 is considered to be the minimum achievable value. 
Indeed, in order for ϕ_T^* to be determined, one needs to 
be in the electroproduction regime where Q^2 ≥ 1 GeV^2. In the bottomonium case, Q̂^2 should thus be larger than 1 GeV^2/(4 m_b^2) ≈ 0.01. 



The figure indicates that there is much uncertainty in the LO result. It also indicates the precision needed at the EIC in order to differentiate among the various fits and to improve on them. A determination of B_UU/A_UU at the 10% level would be an improvement of the current situation. Assuming h_1^⊥  g is 10% of its maximal value at EIC energies (for a more detailed analysis, see section <ref>), this translates into a percent level accuracy requirement on the measurement of ⟨cos 2ϕ_T^* ⟩. For other y values, similar conclusions hold. Needless to say, an NLO analysis of the asymmetry will be needed in order to arrive at more accurate predictions for the EIC and for a fully coherent NLO computation with NLO LDMEs. 


























For these measurements, a good P_T^*-resolution at small P_T^* is an important requirement. Small P_T^* applies to the range up to a few GeV for the EIC energies. Therefore, the transverse-momentum resolution in the small transverse-momentum region should be on the order of a few hundred MeV, such that sufficient bins can be selected to map out this region. For the determination of ⟨cos 2ϕ_T^*⟩, a sufficient angular resolution is needed. 




  §.§.§ Quarkonium polarisation in electroproduction within TMD factorisation









If the polarisation state P (L or T) of the produced quarkonium can be determined in the semi-inclusive quarkonium production process, e   p → e'  J/ψ  (Υ)   X at small transverse momentum, P_T^*, then that may offer a further possibility to improve our knowledge on LDMEs. As an illustration, here we consider the example of the ratio of the ϕ_T^*-integrated cross sections:

    ∫ d ϕ_T^*  dσ^UP/d y  dx_B  d^2 P^*_T/∫ d ϕ_T^*  dσ^UU/d y  dx_B  d^2 P^*_T = ∑_n A_UP^[n]  C[f_1^g Δ^[n]]/∑_n A_UU^[n]  C[f_1^g  Δ^[n]] = A_UP/A_UU.




Let us stress that eq:cos2phibis relies on the assumption that the shape functions are equal to the corresponding LDMEs times a universal shape function 
that is also polarisation independent. 
If so, the ratios A_UL/A_UU and A_UT/A_UU are independent of the value of 
P^*_T≡ |P^*_T| to all orders and hence not affected by TMD evolution. The ratio will only receive contributions from higher orders through modification of the amplitudes. Thus far only the 
LO expressions are known <cit.>: A_UU was already given in Section <ref>, and

    A_UL    = 1/3 [1+(1-y)^2]  ⟨ O^[^1S^[8]_0]⟩ +  [2(1-y) 1+10Q̂^2+Q̂^4/(1+Q̂^2)^2 +  y^2 1+2Q̂^2+Q̂^4/(1+Q̂^2)^2] ⟨ O^[^3P_0^[8]]⟩/m_Q^2 ,

where A_UT= A_UU - A_UL. Compared to A_UU, the ⟨ O^[^3P_0^[8]]⟩ term in A_UL has different inelasticity y and Q̂^2 dependences. This implies that there can be a significant deviation of A_UL from A_UU/3 (and of A_UT from 2A_UU/3), signalling the production of polarised quarkonia. Likewise, one could consider the ratios B_UL/B_UU or B_UT/B_UU which are similar, but different linear combinations of LDMEs. 


In fig:ratioBUUoverAUUandratioAULoverAUU (right plot) the ratio A_UL/A_UU at LO is shown for the LDME fits <cit.> at LO (SV) and <cit.> (here denoted BCKL, B14 in tab:LDME-fits_comp) at NLO, including uncertainty bands, assuming again uncorrelated uncertainties on the ⟨ O^[^1S_0^[8]]⟩ and ⟨ O^[^3P^[8]_0]⟩ determinations. 

In reality the uncertainties are correlated, which means that the bands are expected to be overestimations. The difference between the central values of the two different LDME sets could be viewed as another measure for the size of the involved uncertainties.
Although both fits are compatible with unpolarised production, 
both fits also allow, within their uncertainties, for values considerably different from 1/3. It is important to recall that quarkonium-polarisation observables are very sensitive to radiative corrections <cit.>. Computations of the NLO corrections to the hard parts entering the ratio A_UL/A_UU are therefore necessary to perform a reliable extraction of the LDMEs from these ratios.








In Ref. <cit.>, the fit C12 <cit.> is used to demonstrate the dominance of the ^1S_0^[8] cc̅ state in the inclusive process e h → J/ψ  X (which is dominated by Q^2 ≈ 0) described in collinear factorisation. As a result, it is concluded that the J/ψ will be approximately produced in an unpolarised state. However, the above results show that due to the large uncertainties in the CO LDMEs, one cannot draw the same conclusion for semi-inclusive J/ψ electroproduction in the TMD regime (where P_T^*2 is much smaller than the two hard scales M_J/ψ^2 and Q^2). 
Observation of 
a non-zero polarisation of the J/ψ yield would signal the relevance of the P-wave LDME ⟨ O^[^3P^[8]_0]⟩ or of higher-order contributions.

Again a 10% level precision of the determination of the ratio A_UL/A_UU at EIC would be sufficient to improve on the present situation. For this the polarisation state L or T of the quarkonium needs to be determined with sufficient precision of course.











 §.§ On the importance of final-state effects on quarkonium formation in electron-nucleus collisions 







Interest in quarkonium formation 
in reactions with nuclei goes back more than 30 years in the context of heavy-ion reactions.  The colour interaction between heavy quarks immersed in a high temperature quark-gluon plasma (QGP), such as produced in these reactions, was predicted to be screened, preventing quarkonium states from forming as well as dissociating them <cit.>.  

Excited, weakly-bound-state solutions to the Schrödinger equation, such as Υ(2S), ψ', χ_c, 
were expected  to melt  first in the QGP and provide a “thermometer” for determining the plasma temperature. Since their introduction, these ideas have evolved significantly. It was realised that dissociation and formation suppression of hadrons in QCD matter is not limited to quarkonia. Open heavy-
flavour mesons have short formation times and can also be destroyed by collisional interactions in the nuclear medium <cit.>, reducing the experimentally measured cross sections.  Importantly, the breakup of J/ψs and Υs is not exclusive to the QGP and can take place in different forms of strongly-interacting matter, for example a hadron gas or a large nucleus.  
Measurements of the modification of charmonium and bottomonium production in dAu, pAu, pAl and pPb collisions at RHIC <cit.> and at the LHC <cit.>, respectively, showed that production suppression increases with the multiplicity of hadrons recorded in a reaction <cit.>. 

Moreover, recent measurements of bottomonium yields in pPb collisions showed <cit.> that excited Υ states are more suppressed than the ground state, and the hierarchical pattern becomes more manifest in the negative rapidity region, which is the direction of the lead nucleus <cit.>. 
Studies indicate that final-state interactions can play a significant role in reducing the rates of quarkonium production at the EIC. This quenching effect has been demonstrated for light and heavy mesons (containing a single heavy quark) <cit.> and inclusive and heavy flavor-tagged  jets <cit.>. At forward rapidities and, especially at lower 
centre-of-mass energies, suppression in cold nuclear matter can be as large as a factor of two and serves as strong motivation to investigate these effects for quarkonium final states. 

These observations and predictions indicate that final-state effects (interaction of quarkonium with co-moving hadrons as well as the remnant of the nucleus) require careful treatment in order to extract information about nuclear PDF and transport properties of the nuclear matter. This 
can be addressed from both experimental and theoretical perspectives.

Experimentally, one can approach this concern by studying femtoscopic correlations (two-particle correlations at low relative momentum) between quarkonium and hadron in ep and eA reactions. Such observables are sensitive to interactions in the final state and strong interaction parameters can be measured directly (the scattering length and effective range) <cit.>.  Quarkonium-hadron elastic and inelastic scattering cross sections can be evaluated as a function of event multiplicity. Such information can be used to calculate the modification of the quarkonium yield in the hadronic environment.

 






From the theory point of view, in order to extract nuclear PDFs and constrain the transport properties of large nuclei using quarkonium production in eA collisions, we need to develop a theoretically well-controlled framework 
capable of describing final-state interactions. 
Below, we briefly present an example of such an attempt.



Since the remnant of the nucleus is a cold nuclear environment, we expect the energy transferred between the nucleus remnant and the heavy-quark pair traversing the nucleus to be small. With this assumption, one can use the open quantum system framework and the Boltzmann equation developed in Refs. <cit.> to study final-state interactions. In this approach the physical quantity that encodes the essential information of the nuclear remnant relevant for a final-state interaction is the chromoelectric field correlator, which is defined in a gauge invariant way:

    g_E^>(q)    =   ∑_i=1^3 ∫ d^4 (y-x)   e^iq · (y-x)  Tr_N( E_i(y) W E_i(x) ρ_N )  ,

where ρ_N is the density matrix of the remnant nucleus and W denotes a staple-shape Wilson line in the adjoint representation that connects the spacetime points y and x such that the correlator is defined gauge invariantly. For quarkonium dissociation and formation, the two time-like Wilson lines are connected at positive and negative infinite times separately, as shown in fig:eA_wilson. The Wilson lines involved here are similar to those involved in the definition of proton TMDs, with a 
difference in the orientation of the Wilson lines. 



In a nutshell, quarkonium production in eA collisions involves both initial-state and final-state effects. It will be important for the community to develop strategies for how to best separate these distinct contributions. The combination of both eA and pA experimental data will be useful to determine quarkonium-hadron interaction parameters, nuclear PDFs, and properties of the remnant nucleus such as the chromoelectric correlator strength.







§ QUARKONIA AS TOOLS TO STUDY THE PARTON CONTENT OF THE NUCLEONS







The goal of the present section is to show that quarkonium production in lepton-hadron collisions can be an excellent observable to probe the partonic content of the nucleon. 

First, we discuss 
how quarkonium production measurements at the EIC can contribute to our knowledge of collinear PDFs of the nucleon. Section <ref> is dedicated to accessing the gluon PDF from inclusive -quarkonium photoproduction processes. In Section <ref>, we emphasise how measurements of exclusive and electroproduction at the EIC, by extension of those from the HERA collider, can be used as an indirect probe of the gluon PDF at moderate values of the momentum fraction over a wide range of scales. 

Section <ref> is devoted to the sensitivity to light quark PDFs of inclusive photoproduction, while Section <ref> focuses on the charm PDF and the potential detection of intrinsic charm at the EIC. 
We then move to the multidimensional imaging of the partonic structure of nucleons through quarkonium-related measurements at the EIC.

Sections <ref> and
 <ref> are devoted to the possibility to extract information on TMD PDFs of unpolarised and polarised nucleons, respectively, from quarkonium electroproduction data at the EIC.
The systematic description of exclusive production processes is done in terms of GPDs (Section <ref>) and GTMDs (Section <ref>).  We stress that the relation between GPDs and PDFs used in Section <ref> is only an approximation, albeit a good one at the moderate and low values of the momentum fraction that we consider here. Furthermore, in Section <ref>, we touch on the 
possibility  to access the QCD trace anomaly through the measurement of exclusive electroproduction at the threshold. 

Finally, in Section <ref>, we concentrate on double-parton scattering (DPS), which is another interesting probe 
of nucleon structure. First estimates for -pair electroproduction at the EIC, which include DPS contributions, are presented.




 §.§ Unpolarised-nucleon PDFs





  §.§.§ Gluon PDF from inclusive quarkonium photoproduction









Inclusive photoproduction, when an almost real photon hits and breaks the proton producing a , is a useful tool to study the quarkonium-production mechanism and to 
probe the gluon PDF. 
This process has been the object of several studies at HERA <cit.>, and, in the future, it could be 
studied at the EIC. 

In Ref. <cit.> 
the inclusive photoproduction up to NLO in QCD for J/ψ and Υ(1S) at lepton-proton colliders was revisited, 
focusing on the P_T- and z-integrated yields.
Like for other charmonium-production processes <cit.>, 
the appearance  of negative hadronic cross sections was observed at increasing energies, due to large negative partonic cross sections. There can only be two sources of negative partonic cross sections:  the interference of the loop amplitude with the Born amplitude or the subtraction of the IR poles from the initial-emission collinear singularities to the real-emission amplitude. 
Here, the latter subtraction is the source of the negative cross sections.  
Conventionally, such divergences are removed by subtraction and included to the PDFs via Altarelli-Parisi counterterms (AP-CT). In principle, the negative term from the AP-CT should be compensated by the evolution of the PDFs according to the DGLAP equation. Yet, for the μ_F values on the order of the natural scale of these processes, the PDFs are not evolved much and can sometimes be so flat for some PDF parametrisations that the large ŝ region still significantly contributes. This results in negative values of the hadronic cross section.


To solve the negative cross-section issue, 
the μ̂_F prescription proposed in <cit.> was used, which, up to NLO, corresponds to a resummation of such collinear divergences in 
HEF <cit.>. According to this prescription, one needs to choose μ_F such that, for the partonic cross section σ̂_γ i (i = q, q̅, g), lim_ŝ→∞σ̂_γ i^ NLO=0
. 
It was found that, for z<0.9, the optimal factorisation scale is μ̂_F=0.86 M_ <cit.> which falls well within the usual ranges of used values. Like for η_c hadroproduction, such a factorisation-scale prescription indeed allows one to avoid negative NLO cross sections, but it of course in turn prevents one from studying the corresponding factorisation-scale uncertainties.


The NLO μ_R uncertainties become reduced compared to the LO ones but slightly increase around √(s_γ p) = 50 GeV
, because of rather large (negative[Let us stress that unless μ_R is taken very small with a large α_s(μ_R), these negative contributions are not problematic, unlike the oversubtraction by the AP-CT.]) interferences between the one-loop and Born amplitudes.
At NNLO, 
a further reduction of the μ_R uncertainties is expected. This is particularly relevant especially around √(s_γ p) = 50-100 GeV, which corresponds to the EIC region. This would likely allow us to better probe gluon PDFs using photoproduction data. Going further, differential measurements in the elasticity or the rapidity could provide a complementary leverage in x to fit the gluon PDF, even in the presence of the  v^4 

CO contributions. Indeed, these would likely exhibit a very similar dependence on x. 


The possibility to constrain PDFs using future J/ψ and Υ(1S) photoproduction data <cit.>
was investigated by comparing the PDF and μ_R uncertainties.

Unsurprisingly, the PDF uncertainties get larger than the (NLO) μ_R uncertainties with the growth of the γ p centre-of-mass energy, in practice from around 300 GeV,  for x below 0.01. Although this is above the reach of the EIC, with NNLO predictions at our disposal in the future, with yet smaller μ_R uncertainties, one could set novel constraints on PDFs with such EIC measurements.  Following the estimated counting rates for 100 fb^-1 of ep collisions given in <cit.>, 
a number of differential measurements (in , z and/or y) will be possible to reduce the impact of highly- or even partially-correlated theoretical uncertainties, including the contamination of higher-v^2 corrections, such as the 
CO contributions.



tab:numb_of_particles gathers estimates of the expected number of J/ψ and Υ(1S) possibly detected at the different ep 
centre-of-mass energies at the EIC. 

For Υ(1S), the yields should be sufficient to extract cross sections even below 
the nominal EIC luminosities. 

One can also estimate the expected number of detected ψ', Υ(2S)
and Υ(3S) using the following relations

    N_ψ'   ≃0.08 ×  N_J/ψ,
     
        N_Υ(2S)   ≃0.4 × N_Υ(1S),
     
        N_Υ(3S)   ≃0.35 × N_Υ(1S),

derived from the values of[These contributions were estimated using |R_ψ'(0)|^2 = 0.8 GeV^3, |R_Υ(2S)(0)|^2 = 5.0 GeV^3  and |R_Υ(3S)(0)|^2 = 3.4 GeV^3 and the corresponding measured branching fractions to J/ψ and Υ (1S) <cit.>.] |R_(0)|^2 (the quarkonium radial wave function at the origin, that is related to the ^3S_1^[1] LDME) and of the branching fractions to leptons. Using the values in tab:numb_of_particles and eq:yield_ratio <cit.>, one can see that the yield of ψ' should be measurable and the yields of Υ(2S) and Υ(3S) are close to about half of that of Υ(1S) and should be measurable as well at the EIC. 




  §.§.§ Gluon PDFs from exclusive quarkonium photo- and electroproduction






The exclusive production of heavy vector mesons has long been a fascinating observable to study, functioning as an enticing avenue to unravelling the small-x behaviour of the gluon PDF from low to moderate scales. Measured in the first instance in the fixed-target mode <cit.> and in DIS events at HERA, see e.g. <cit.>, and then more recently in ultra-peripheral collisions at the LHC <cit.>, they provide a means to explore the quarkonium production mechanism and act as sensitive probes at the frontier of small-x saturation physics.

The exclusive  electroproduction, γ^* p → p, has been measured via dilepton decays at HERA in a narrow range of photon virtualities, extending up to ⟨ Q^2 ⟩ = 22.4 GeV^2. The corresponding photoproduction has also been determined in ultraperipheral events at the LHC. 
There are not, as of yet however, any data from HERA and the LHC for exclusive  electroproduction,  γ^* p → p, away from the photoproduction limit. Going forward, the EIC will extend the kinematic reach in Q^2,  providing a lever arm up to larger virtualities and, moreover, allow for a measurement of the  electroproduction with off-shell photon kinematics for the first time, albeit with a projected lower Q^2 + M_^2 bin coverage and event count rates due to its heavier mass <cit.>. 

Recently, in Ref. <cit.>, the coefficient functions for exclusive heavy vector meson electroproduction were derived at NLO within the framework of collinear factorisation, with the transition from an open heavy quark-antiquark pair to a bound heavy vector meson made within LO NRQCD. 






Based on the above derivation of the coefficient functions, predictions for the exclusive electroproduction cross section have been made. They are shown in fig:chris_1 in bins of Q^2 at a fixed centre-of-mass energy, W = 90 GeV, of the γ^* p pair.

We use the Shuvaev transform <cit.> as a reliable means to obtain the GPD from input PDFs in the kinematic regions shown. We construct GPDs in such a way using MSHT20 <cit.>, NNPDF3.1 <cit.> and CT18 <cit.> input NLO PDFs and the predictions based on the former are shown in the figure. The choice of input PDF has the largest effect at the lowest Q^2, where the choice of the initial condition of the DGLAP evolution is felt, while for larger Q^2, this effect washes out and the predictions based on each PDF set agree at or below the percent level.
The central values of the prediction for low to moderate Q^2 are in good agreement with the experimental data from H1 and ZEUS, but for larger Q^2 there appears to be a downward shift of the prediction from the data. The prediction in the highest Q^2 bin exhibits a small factorisation-scale dependency and is essentially independent of the choice of the input PDF but, as shown, the deviation from the data is sizeable. 
Interestingly, in the large Q^2 limit, the gluon amplitude ∝ln(Q^2/m_Q^2)^2 while the quark amplitude ∝ln(Q^2/m_Q^2). This observation seems to necessitate a program of resummation for the exclusive electroproduction of heavy vector mesons for virtualities Q^2 ≫ m_Q^2,  those relevant for EIC kinematics, and may provide for the reconciliation of the theory prediction and experimental data at large scales.

The data statistics are currently limited for larger Q^2 and, in particular, there is a wide range 

where the EIC can provide a first coverage. This will help to ascertain on which front the difference between this prediction and the data at large Q^2 lies and if resummation effects are already needed. Other numerical effects in this framework such as the so-called `Q_0 subtraction' <cit.>, crucial for a fruitful description of the photoproduction data <cit.>, are not surmised to be important for electroproduction kinematics because the corresponding power correction 𝒪(Q_0^2/μ_F^2) is no longer of 𝒪(1). See also <cit.> for a recent baseline study of exclusive photoproduction in heavy-ion collisions in the collinear factorisation framework to NLO. 



Simulated event count projections were given for the exclusive electroproduction of the  and  in bins of Q^2 + M_^2 as a function of x in <cit.>.  In fig:chris_2 (left panel), we show predictions for the exclusive  electroproduction cross section as a function of W at a fixed scale ⟨ Q^2 ⟩ = 16 GeV^2 using Shuvaev-transformed MSHT20 input NLO PDFs, as well as the exclusive  electroproduction HERA data that lie in this bin for comparison purposes. The prediction agrees most favourably with the more up-to-date dataset, however the EIC will be able to provide more statistics and resolve the slight tension between (and discrepancies within) the datasets. In particular, the data point at W = 189 GeV is around a factor of two larger than other data lying in this bin.  We also show predictions for the exclusive electroproduction cross section as a function of W (right panel) for ⟨ Q^2 ⟩ = 0.001, 16, 22.4 GeV^2 and 47.3 GeV^2, which may ultimately be compared with data from the EIC.[Admittedly, the expected event count rate is a lot lower than that of the corresponding  production, even by three orders of magnitude in the photoproduction bin containing the most counts <cit.>. Any data will therefore likely be sparse and exhibit large uncertainties, but nonetheless complement those already existing from HERA and LHC, shown in the right panel of Fig. <ref>.] In each case the quark contribution to the total amplitude is negligible and so the forthcoming enhanced statistics and increased data coverage from the EIC will allow for refined and improved constraints on the gluon PDF at low to moderate scales. 



  §.§.§ Light quarks








At EIC energies, we also expect to be sensitive to quark-initiated partonic subprocesses. As shown in Ref. <cit.>, in inclusive quarkonium photoproduction, the quark-induced subprocesses γ + q → + q (+  g) will be a relevant contribution to the  cross section. Therefore, through quarkonium photoproduction, the EIC will also 
be 
partially sensitive to the light-quark PDFs.






To highlight the quark-induced contribution, we show in fig:EIC-CT14NLO-NLOstar-ratio the ratio to the 
CSM cross section for every partonic subprocess (up to 𝒪(αα_s^3)) depicted in fig:EIC-CT14NLO-NLOstar, at two different 
centre-of-mass energies, √(s_ep) = 45 GeV (left panel) and √(s_ep) = 140 GeV (right panel), as a function of transverse momentum. It is clear that 
the pure QED quark-initiated process at 𝒪(α^3) become dominant at high P_T, accounting for 
over half of the 
cross section at P_T ∼ 12 (16) GeV at √(s_ep) = 45 (140) GeV. The effect is larger at √(s_ep) = 45 GeV, where the valence region of the PDF is probed. The 𝒪(αα_s^3) contribution (γ + q → + q + g) is roughly 
5-15 % and 10-15 %  of the 
cross section at √(s_ep) = 45 GeV and √(s_ep) = 140 GeV, respectively. We then expect that, in photoproduction processes at the EIC, the produced at large P_T will be recoiling off of 
at least one quark jet.  The significant contribution of quark-induced subprocesses at high 
 of the J/ψ is also observed in the NLO NRQCD calculation, as shown in 
fig:quark_fractions. Moreover, this conclusion depends only mildly on the 
NRQCD LDMEs that were used. 








  §.§.§ Charm quark and intrinsic charm












The existence of a nonperturbative charm-quark content in the proton, referred to as intrinsic charm (IC), has long been postulated <cit.>. Intrinsic charm states are a fundamental property of hadronic bound-state wave functions <cit.>. They differ from extrinsic charm in perturbative QCD that arises from gluon splitting and contributes to the heavy-quark PDFs  (i.e., radiatively generated). The “intrinsic" label is due to the fact that a cc̅ pair formed by gluons from more than one quark line forces the cc̅ parameters to be dependent upon (i.e., reflective of) the hadron that creates it.  Therefore the c and c̅ distributions are “intrinsic” to the identity of the proton, or the meson, or whichever hadron contains the bound quarks that emit gluons.  “Extrinsic” means that the sea quark pairs come from a single quark line gluon and therefore do not reflect the bound state structure they exist in, at least not in the clear way that IC of the proton does, peaking at ∼ x_B=0.4 and imparting a difference in c and c̅ distributions, according to recent lattice calculations <cit.>.  

Since extrinsic charm contributions are due to a gluon emitted by a single quark line which then splits into a cc̅ pair, these charm distributions are soft, appear at low x and depend logarithmically on the mass of the heavy quark m_Q. On the other hand, IC contributions dominate at higher x and have a 1/m_Q^2 dependence. They come from five-quark (and higher) Fock-state configurations of the proton, |uud c c⟩, and are kinematically dominated by the regime where the state is minimally off-shell, leading to equal-rapidity constituent quarks. Thus, the charm quarks are manifested at large x. When the proton in this state interacts with its collision partner, whether a hadron or a lepton, the coherence of the Fock components is broken and the fluctuations can hadronise <cit.>. In hadroproduction, the state can be broken up by a soft gluon from the target interacting with the proton. In ep interactions, instead of a soft gluon, a low-energy photon can play the same role and bring the state on mass shell.

Several formulations of intrinsic charm in the proton wave function have been proposed. The first was proposed by Brodsky and collaborators in <cit.>:

    dP_ ic  5/dx_1 dx_2 dx_3 dx_c dx_c = P_ ic 5^0
    N_5 ∫ dk_x  1⋯ dk_x   5∫ dk_y  1⋯ dk_y   5δ(1-∑_i=1^5 x_i)δ(∑_i=1^5 k_x   i) δ(∑_i=1^5 k_y   i)/(m_p^2 - ∑_i=1^5 (m_i^2/x_i) )^2   ,

where i = 1, 2, 3 are the light quarks (u, u, d) and i = 4 and 5 are the c and c quarks.
Here, N_5 normalises the |uud c c⟩ probability to unity and P_ ic  5^0 scales the unit-normalised probability to the assumed intrinsic-charm content of the proton. The delta functions conserve longitudinal and transverse momentum. The denominator of eq:icdenom is minimised when the heaviest constituents carry the dominant fraction of the longitudinal momentum, ⟨ x_Q ⟩ > ⟨ x_q ⟩. 
In the
first papers, the c and c distributions were treated equally, but 
later studies showed 
 
an asymmetry in c and c distributions <cit.>.  The asymmetry is caused by QCD diagrams where, for example, two gluons from  two different valence quarks in the nucleon couple to a heavy-quark pair gg → QQ with charge conjugation value C=+1 <cit.>.
This amplitude interferes with QCD diagrams where an odd number of gluons attach to the heavy-quark pair, g→ QQ and ggg→ Q Q with C=-1.  The interference of amplitudes with the same final state but different charge conjugation symmetry for the Q Q produces the asymmetric distribution functions.  The analogous interference term is seen in the electron and positron distributions in e^+ e^- pair production <cit.>.

At leading order, the charm-quark structure function from this state can be written as

    F_2c^ ic(x_c) = 8/9 x_c c(x_c) = 8/9∫ dx_1 dx_2 dx_3
        dx_cdP_ ic  5/dx_1 dx_2 dx_3 dx_c dx_c .






  
Intrinsic-charm models Intrinsic-charm distributions in the proton have also been calculated using meson-cloud models where the proton fluctuates into a D(u c) Λ_c (udc) state <cit.>. A further development of this model examined all possible charm meson-baryon combinations  in the |uud c c⟩ state <cit.>, finding that charm mesons would predominantly be produced through D^* mesons. In these models the charm sea contribution would be asymmetric xc(x) ≠ x c(x).
In both the Brodsky et al. and the meson-cloud formulations, the intrinsic-charm contributions appear as an enhancement at large x.  On the other hand, a sea-like distribution <cit.> has also been considered. In this case, the intrinsic-charm distribution is represented simply as an overall enhancement to the light-quark-mass sea. These distributions are 
symmetric, x c(x) = x c (x).

Intrinsic-charm distributions from these models have been included in global analyses of the parton densities <cit.>. Earlier analyses <cit.> focused specifically on the European Muon Collaboration (EMC) high-x and high-Q^2 data <cit.>.  A range of values of P_ ic  5^0 were extracted, from 0.1% to 1%. For more details of these analyses, see <cit.>.  See also the recent review in <cit.> for more applications of intrinsic-heavy-quark states. New evidence for a finite charm-quark asymmetry in the nucleon wave function from lattice gauge theory, consistent with intrinsic charm, was published in <cit.>.  Further evidence for unequal c and c distributions in the proton has recently been presented along with proposed experimental tests with the EIC using flavour-tagged structure functions <cit.>.

Note that only the 5-particle intrinsic-charm state of the proton has been discussed.  However, one can also consider higher Fock components such as |uud cc q q⟩.  These will reduce the average momentum fraction of the charm quark and also have lower probability. See e.g. <cit.> for examples of charm hadron distributions from higher Fock states. Finally, the possibility for an enhanced IC component in the deuteron was studied in <cit.>.















  
Recent hints from the LHC 
A number of experimental measurements <cit.> over the last several decades have provided tantalising hints of intrinsic charm. 
Recently LHCb announced that their measurement of Z +  charm   jets relative to all Z +  jets is consistent with an intrinsic-charm component of the proton as large as 1% at large Z rapidity <cit.>. These results were recently confirmed by a phenomenological analysis made by the NNPDF Collaboration <cit.>. Measurements at lower scales than the Z-boson mass are therefore eagerly awaited for to advance our understanding of this higher-Fock-state phenomenon.














  
Intrinsic charm at the EIC The 
EIC will offer the possibility to probe the nonperturbative charm-quark content in the proton. Recent studies show that the EIC will be capable of precision studies of intrinsic-charm as well as gluon distribution functions in the nucleus and in the nucleon <cit.>. 

The associated production of a  and a charmed particle is 
an additional potential probe of intrinsic-charm related effects. A leading order VFNS study, first made in <cit.> for quarkonium hadroproduction, has been extended in <cit.> to the case of  photoproduction.  Such a scheme allows a proper merging of different partonic contributions, namely γ + g → + c + c̅ and γ + {c,c̅}→ +  {c,c̅}, respectively calculated with 3 and 4 flavours in the proton, using a counter term, dσ^ CT, that avoids double counting. When the charm-tagging efficiency ε_c is taken into account, the corresponding VFNS cross section is given by:


    dσ^ VFNS = dσ^ 3FS[1-(1-ε_c)^2] + (dσ^ 4FS - dσ^ CT)ε_c.


Based on such computations, the +charm yield has been calculated for two different EIC configurations: √(s_ep) = 45 (140) GeV, taking into account a 10% charm-tagging efficiency <cit.>. The calculation has been done with the CT14NNLO PDF set <cit.>, which includes different eigensets with some IC effects: a “sea-like” (in green in the following), 
a “valence-like” (in red) also called “BHPS,” and a central eigenset with no IC effects which we refer to as “no IC” (in blue). 



fig:EIC-VFNS-IC shows the result for the +charm yield at the 
EIC. First, we note that, at √(s_ep) = 45 GeV (left panel in fig:EIC-VFNS-IC), the yield is limited to low P_T values even with the largest estimated integrated luminosity. Nonetheless, it is clearly observable if ε_c = 0.1 with 𝒪 (5000, 500, 50) events for ℒ = (100,10,1) fb^-1. On the other hand, at √(s_ep) = 140 GeV (fig:EIC-VFNS-IC, right panel), the P_T range extends to ∼ 14 GeV and we expect 𝒪(10000) events at ℒ = 100 fb^-1. Such events could be observed by measuring a charmed jet. Finally, we note that, at √(s_ep) = 140 GeV, where the valence region 
is not probed, no clear IC effect is visible, while at √(s_ep) = 45 GeV we expect a measurable effect, where the BHPS valence-like peak is visible with a yield enhancement as large as 5-6 times the “no IC” yield. The EIC at √(s_ep) = 45 GeV will thus be the place to probe the nonperturbative charm content of the proton via associated +charm production.













































































 §.§ Unpolarised-nucleon TMDs










  §.§.§ Unpolarised gluons






Quark TMDs have now been extracted from data with reasonable precision <cit.>. On the contrary, phenomenological studies of gluon TMDs are still very much at the beginning stage. In Ref. <cit.>, a gluon TMD description of the Higgs-production transverse-momentum spectrum was compared to data,  

which, however, suffers from very large uncertainties. In Refs. <cit.>, a gluon TMD description of the LHCb -pair-production data <cit.>  was obtained. Like for Higgs-boson production, the experimental errors are 
large and require the subtraction of double parton scattering contributions (see <cit.> and Section <ref>), which adds an additional uncertainty. In Ref. <cit.>, it was discussed that back-to-back production of a heavy quarkonium, in particular of an Υ, and an isolated photon in proton-proton collisions at the LHC is a promising way to access the distribution of both the transverse momentum and the polarisation of gluons inside unpolarised protons. In a wide range of invariant masses of the quarkonium and photon system, gluon-gluon scattering into a photon plus a quarkonium in the CS state dominates. 

In the aforementioned processes, one 
however 
probes a convolution of  two gluon TMDs. At the EIC, one can probe gluon TMDs more directly through the P_T^* distribution, although upon the inclusion of 
ShFs (see section <ref>); this also deals with convolutions. 
At the LHC, with the consideration of such 
ShFs, one even folds three transverse-momentum-dependent distributions. Another possibility to study gluon TMDs at the EIC using quarkonia is to consider the transverse-momentum imbalance between the scattered lepton and the observed J/ψ in the electron-hadron centre-of-mass frame p_T = |ℓ'_T + P_T|. If large |ℓ'_T|≃ |P_T| ≫ | p_T| determines the hard scale of the process, then 
in quarkonium production at the EIC the leading subprocess is 
e+g→ (cc̅)^[8]+e with the octet pair 
hadronising into an observed J/ψ. Within the hybrid factorisation formalism for SIDIS discussed in Section <ref>, the p_T 
should be determined by the transverse momentum k_T of the colliding gluon (or its TMD distribution) and the quarkonium TMD 
ShF. 
Since gluon radiation from a heavy quark should be strongly suppressed compared to a light quark or a gluon, the observed momentum imbalance p_T is expected to be dominated by the k_T of the colliding gluon <cit.>. Therefore, the p_T
-distribution of production in SIDIS could be a more direct observable for the gluon TMD <cit.>.


It would be very interesting to compare the gluon TMD obtained at EIC to that from the J/ψ+J/ψ or Υ+γ process at LHC in the future. 

In principle, gluon TMDs are process dependent, even in the unpolarised case (see e.g. <cit.>).
However, provided that the CS final state dominates in 
J/ψ+J/ψ and Υ+γ production at the LHC, 
these processes involve the same gluon TMD. This then would provide a nice test of TMD factorisation in combination with NRQCD and of TMD evolution, if the processes are probed at different scales. Another comparison that seems worthwhile is the extraction of gluon TMDs from open heavy-quark pair production at the EIC <cit.> or from inclusive η_c or η_b production in proton-proton collisions <cit.>. Note that inclusive CS  or Υ production from two gluons is forbidden by the Landau-Yang theorem, while inclusive CO  or Υ production does not involve the same gluon TMD and may not even 
factorise to begin with.



















































  §.§.§ Linearly polarised gluons







As discussed in Section <ref>, linearly polarised gluons lead to a cos 2 ϕ_T^* asymmetry in  semi-inclusive electroproduction of J/ψ in unpolarised ep collisions <cit.>. In this section, we present some predictions for this asymmetry 
at low transverse momenta ^*. 








 
 
 
 













Within NRQCD, contributions to the asymmetry 
comes 
through the fusion of a virtual photon and a gluon <cit.> already at Born order, i.e. α_s α, but at NNLO in v^2 since via CO contributions. Such α_s α contributions however only sit at  z=1. 
As soon as z≠ 1, a recoiling particle against the quarkonium is needed and Born-order contributions are at α^2_s α both from CS and CO states. From a simple counting in v^2 the CS contributions <cit.> should be dominant at z≠ 1. However, the current LDME fits seem not to obey such a simple v^2 counting and, as a matter of fact, sometimes leads to an excess[It should be clear to the reader that such computations are as of now only carried at LO whereas the LDMEs are extracted at NLO. We refer to our introductory discussion at the beginning of Section <ref> regarding potential issues in doing so.] in describing the scarce data available from HERA <cit.>. 
In principle, the asymmetry thus receives contributions from both CS and CO states.


The first estimate we present here is based on 
a model expression for the cross section <cit.>:

    dσ=   1/2sd^3l'/(2π)^32E_l'd^3P_/(2π)^32E_P_∫d^3p_g/(2π)^32E_g∫ dx d^2𝐤_⊥(2π)^4δ(q+k-P_-p_g)
       ×1/Q^4ℒ^μμ'(l,q)Φ^νν'
      (x,𝐤_⊥) ℳ_μν(ℳ_μ'ν')^∗.

This expression is akin to the Generalised Parton Model employed to describe single-spin asymmetries in polarised proton collisions (to be discussed in section <ref>). It is not of TMD-factorisation form and differs from Eq. (<ref>) by considering the subprocess γ^∗+g→+g, where the additional hard gluon in the final state generates larger transverse momenta and elasticity z values below 1, 
while the dependence on the initial gluon transverse momentum is kept everywhere. In other words, no collinear expansion is performed and the obtained expression is thus not a CF expression either. 









In fig:lin, we 
show the cos 2 ϕ_T^* asymmetry as a function of P_T for =140 GeV, for fixed values of z and Q^2. Both  
CS and CO contributions are included. 

We show 
the results for two different models for the TMDs, 
the Gaussian <cit.> 
and the  McLerran-Venugopalan  model <cit.>, and for two different sets of LDMEs, CMSWZ <cit.> and BK <cit.>. The asymmetry is small and 
depends  on the chosen LDME set. 
Details of the calculation may be found in <cit.>.  

 






A second estimate – only relevant for z≃ 1 – is based on the 
TMD formalism involving shape functions. Although  the semi-inclusive quarkonium electroproduction is naturally described in TMD factorisation at small quarkonium transverse momentum  (^* ≪ M_∼ Q),  there is large uncertainty due to the non-perturbative part of the TMD description and due to the lack of knowledge on the TMD shape functions.  However, using the leading-order shape functions in terms of LDMEs and including leading-order TMD evolution, it is nevertheless possible to obtain rough predictions for the EIC (details on the shape function can be found in Ref. <cit.>). Using this approach, estimates for the cos 2ϕ_T^* asymmetry in J/ψ production as a function of ^* can be obtained. The results are shown in fig:EICpredictioninclTMDevolution for several LDME sets (for more predictions see Ref. <cit.>) and for kinematics similar to that of fig:lin (to be precise, for the same  and Q^2, and comparable x_B, but different values of z). Despite the large uncertainties in these TMD results (the uncertainty bands reflect the uncertainty in the non-perturbative Sudakov factor), it is clear that within these uncertainties it allows for significantly (by more than an order of magnitude) larger asymmetries than in fig:lin. 



Its measurement may thus be feasible at EIC such that further constraints on the LDMEs, and more generally on the TMD shape functions, can be obtained in this way. 



Observing a nonzero asymmetry would be a signal of linear polarisation of the gluons inside an unpolarised proton, which is expected theoretically but not established experimentally thus far. The range of predictions is currently too large to draw a definite conclusion about its observability at EIC, but that makes it all the more important to obtain first data on the cos 2 ϕ_T^* asymmetry. It would provide information on the distribution of linearly polarised gluons as well as on LDMEs. 















 §.§ Polarised-nucleon TMDs






Among the observables that can be measured at the EIC to access polarised nucleon TMDs (the Sivers function), 
the most common are probably the Single Transverse Spin Asymmetries (STSA), denoted A_N, or A_UT. 




Two theory approaches have been pushed forward to explain STSAs observed on polarised protons <cit.>. Both of them can in principle be extended to quarkonium production.

The first approach is referred to as collinear twist-3 (CT3) formalism <cit.> 

and, like CF, 
applies to single-scale processes.  
The STSA then arises from  quark-gluon-quark or triple-gluon correlators, which are the sub-leading (in the scale) twist-3 extensions of the usual collinear PDFs (putting aside for this discussion FF contributions). 
Some CT3 analyses for A_N in ep collisions have been performed in the past, see e.g. <cit.>, and only very recently this approach has been extended to STSAs in quarkonium production in polarised ep collisions <cit.>. 


The second approach is TMD factorisation, thus applicable when two very different momenta are measured, or when a small (yet perturbative) momentum is measured in a process involving a large mass ((Λ_ QCD) ≲ P_T^* ≪ Q in SIDIS, where P_T^* is the transverse momentum of the hadron in the final state and Q^2 is the photon virtuality). The STSA arises from the Sivers TMD PDF f_1T^⊥ <cit.>, i.e. the distribution of unpolarised partons inside the transversely-polarised hadron. In the case of quarkonium production in ep collisions, TMD factorisation has been assumed and used to compute the Sivers asymmetry in several cases <cit.>. 


In addition, a phenomenological approach, called the Generalised Parton Model (GPM) <cit.>, encapsulates the Sivers mechanism via the aforementioned TMD Sivers function, assumed to be universal, but also applied in single-scale processes. This is done by keeping track of the transverse-momentum exchanges in the partonic scattering. As such, it can be considered as a hybrid approach between strict CT3 and TMD factorisation. Its extension, called Colour Gauge Invariant GPM (CGI-GPM) <cit.>, allows one to recover the modified universality of the quark Sivers function between SIDIS and Drell-Yan <cit.>. Moreover, for the gluon Sivers effect, similarly to the CT3 approach case, two independent gluon Sivers functions (GSFs) appear <cit.>, dubbed as f- and d-type. 
This approach has proven to be quite successful in phenomenological analyses <cit.>. One should however be careful if one wishes to draw any conclusion about the properties of the used TMDs and the underlying phenomena. In any case, it is useful to get estimates of STSAs in single-scale processes where a CT3 analysis becomes challenging, like for quarkonium production, due to still unconstrained twist-3 functions appearing in its computation. It has been applied to the quarkonium cases in several studies <cit.>. 


Below STSAs in different quarkonium-production processes are discussed, in the context of the EIC, which could perform these measurements by polarising a target. In general, it is believed that quarkonium-related STSA would be key player to underpin the Sivers mechanism for gluons.










Experimentally, one defines the so-called transverse 
STSA as

    A_N = 1/ Pσ^↑ - σ^↓/σ^↑ + σ^↓ ,

where σ^↑ (↓) is the cross section of particles produced with the target nucleon spin orientation upwards (downwards), and P  is the average nucleon polarisation. In what follows, we present predictions and projections for 
STSA in inclusive photoproduction and for azimuthal weighted Sivers asymmetries in leptoproduction in SIDIS processes.







































  §.§.§ EIC reach for A^_N for inclusive photoproduction
















In this section, we study how to probe the GSF via the GPM approach by measuring the STSA in inclusive  photoproduction (γ + p^↑→ J/ψ+X) <cit.>. In such a process, only the f-type GSF contributes to the Sivers asymmetry.



In photoproduction, there are contributions from direct and resolved photons. Resolved photons mainly contribute in the  region of low elasticity z. 
At z close to unity, diffractive contributions become significant.  In inclusive photoproduction, the variable z can be measured using the Jacquet-Blondel method.   The differential cross section of inclusive J/ψ production in unpolarised ep collisions can be written as 

    E_dσ/d^3 P_=   1/2(2π)^2∫ dx_γ dx_g 
    d^2 k_⊥ g
    f_γ/e(x_γ)f_g/p(x_g, k_⊥ 
    g)δ(ŝ+t̂+û-M_^2)
       ×1/2ŝ|ℳ_γ+g→ +g|^2.


Here, x_γ and x_g are the light-cone momentum fractions of the photon and gluon, respectively; ŝ, t̂, û are the partonic Mandelstam variables; ℳ_γ+g→ +g is the matrix element for the partonic subprocess γ+g→ +g; f_g/p(x_g, k_⊥ 
g) is the unpolarised gluon TMD, while f_γ/e(x_γ) is the Weizsäcker-Williams distribution, giving  the density of photons inside the electron <cit.>.  






For theory predictions of measurements on a transversely polarised nucleon, the STSA, as introduced in eq:STSA,  is generally used. 





Some GPM predictions for STSA in inclusive  photoproduction at the EIC for =45 (140) GeV are shown in fig:Photo_Sivers, as a function of the transverse momentum, , as well as a function of the elasticity z. The amplitude for the J/ψ production is calculated in NRQCD. Details of the calculation can be found in Ref. <cit.>. 
The dominating channel of J/ψ production is γ g fusion. The contribution to the numerator of the 
STSA comes mainly from the GSF <cit.>, while the linearly polarised gluons do not contribute to the denominator for this specific process. Moreover, the numerator of the asymmetry only receives contributions from CO states <cit.>, whereas in the denominator, both 
CO and CS contributions are included. 

We have used the GSF parametrisations 
(SIDIS1, SIDIS2) from Ref. <cit.>. BV-a and BV-b are parametrisations of the GSF in terms of up and down quark Sivers functions <cit.>, where parameters from Ref. <cit.> are used. The effect of TMD evolution is not incorporated in the plot. The PDF set MSTW2008 <cit.> is used; the uncertainty bands have been obtained by varying the factorisation scale μ_F ∈[12 μ_0, 2μ_0], with μ_0 = m_T = √(M^2_ + P_T^2) being the  transverse mass. The value of α_s is calculated at the scale μ_0 and is taken from the MSTW set. The used cuts are the following: Q^2<1 GeV^2 and 0.3 < z < 0.9. 

Note that, in the photoproduction case, y coincides with x_γ. The corresponding cut is 0.01 < x_γ < 0.95. 
As shown in fig:Photo_Sivers, we expect A_N to be small and positive in the SIDIS1 and SIDIS2 cases, while it is larger (in size) but negative when the GSF is parametrised in terms of the up- and down-quark Sivers functions. 




Another estimate is shown in jpsi:AN:ep. Here, projections for statistical uncertainties for the  A_N measurement as a function of transverse momentum for ep collisions at =45 GeV and = 140 GeV for an integrated luminosity = 100 are presented.
We consider the reconstruction via its electron decay channel (→, B = 5.94 ± 0.06 %)
, and we assume the single-electron measurement efficiency to be 80% and constant 
with respect to its transverse momentum and in the pseudorapidity interval 
|η|<2. The measurement efficiency is calculated using decay kinematics simulated with PYTHIA8 <cit.> (see <ref> for details). Based on these results, we assume the measurement efficiency to be 64%. Furthermore, we assume the signal-to-background ratio S/B = 1, and use the same method as in Ref. <cit.> to estimate statistical uncertainties on A_N. 
For the expected cross section for prompt production in ep collisions at the EIC, 
we consider the CSM predictions from Ref. <cit.>, which were shown to approximately reproduce HERA data. 
For illustration, the projections are compared to results from  collisions reported by the PHENIX experiment <cit.>. At low , the statistical precision is at the per-cent level, exceeding the quality of the corresponding pp data. In this range, the final uncertainty will be dominated by systematic effects. The uncertainties increase fast with increasing  of because the  spectrum is predicted to be rather steep. Nonetheless, such a measurement would be valuable for constraining gluon TMDs at low transverse momentum. 




Finally, we suggest that the associated photoproduction of and a jet, having them back-to-back, 
can also probe the 
GSF <cit.>. In this case the produced can have large transverse momentum, and needs not be in the forward region. A wide kinematical region can be covered by varying the invariant mass of the J/ψ-jet pair.




  §.§.§ Azimuthal asymmetries for  production in SIDIS at the EIC







In this section we consider the Sivers effect in the SIDIS process
, e(l)+p^↑(P_N)→ e(l^')+ J/ψ (P_)+X, that represent a promising tool to probe the GSF. The weighted Sivers asymmetry for such a process is defined as



    A^sin(ϕ_T^*-ϕ_S^*)_N  ≡  2∫ϕ̣_S^*ϕ̣^*_T sin(ϕ^*_T -ϕ^*_S)
    (σ̣^↑ - σ̣^↓)/∫ϕ̣^*_S ϕ̣^*_T (σ̣^↑ + σ̣^↓)≡∫ϕ̣^*_S ϕ̣^*_T sin(ϕ^*_T -ϕ^*_S) Δ̣σ(ϕ^*_S, ϕ^*_T)/∫ϕ̣^*_S ϕ̣^*_T σ̣ ,

where σ̣^↑(↓)=σ̣^↑(↓)/Q̣^2  y ^̣2P_T^* ẓ is the differential cross section with the initial proton polarised along the transverse direction ↑(↓) with respect to the lepton plane in the γ^*p centre-of-mass  frame (at an angle ϕ^*_S). 



We start by presenting the predictions in the CT3 formalism. In Ref. <cit.>, the twist-3 contributions to the unpolarised and polarised cross sections (respectively denominator and numerator of Eq. (<ref>)) were computed in the CSM. Among the different contributions, one give access to gluon Sivers effect via the CT3 gluon Qiu-Sterman function, which at LO is realated via an integral relation to the GSF first k_⊥ moment. Predictions for the gluon Sivers asymmetry at the EIC at = 45 GeV are presented in Fig. <ref>. They are computed at Q^2 = 10 GeV^2, x_(B) = 0.005 and P_T = 2 GeV, and are presented as a function of z for two different models of the gluon Qiu-Sterman function. Both models are proportional to f_g/p(x), the unpolarised collinear gluon PDF, and read

    Model 1:    0.002 x f_g/p(x) ,
    Model 2:    0.0005 √(x) f_g/p(x) .




Notice that, as these CSM predictions are ratios of cross sections, they do not depend on the value of the CS LDME. Both models predict a sizeable Sivers asymmetry, with a steady increase as a function of z, reaching up to ∼ 13-14% at z = 0.8.

Another prediction for the Sivers asymmetry is performed within 
the GPM at  
αα_s^2. 

In order to study the effects of initial- and final-state interactions (ISIs and FSIs) on the Sivers asymmetry, the 
CGI-GPM approach <cit.> is employed.

In Ref. <cit.>, the same observable was studied at 𝒪(αα_s) within the GPM, which implies z=1. Here the analysis is extended to the region z<1.

Assuming TMD factorisation within the GPM framework, the unpolarised differential cross section, entering the denominator of eq:AN, can be written as

    σ̣/Q̣^2 ỵ ^̣2 P_T^* ẓ=1/(4π)^4 zs∑_a∫x̣_a/x_a ^̣2 k_⊥ a δ(ŝ+t̂+û-M_^2+Q^2) ∑_n1/Q^4 f_a/p(x_a,k_⊥ a) L^μν
    H^a,U_μν[n] ⟨ O^J/ψ[n] ⟩ ,


where a=g,q, q̅ and H^a,U_μν[n] is calculated at the perturbative order αα_s^2 using NRQCD. More precisely, it is the squared amplitude of the partonic process γ^∗+a→ cc̅[n] +a, averaged/summed over the spins and colours of the initial/final parton, with
n=^3S_1^[1,8], ^1S_0^[8], ^3P_J^[8], J=0,1,2.
L^μν is the standard leptonic tensor  and ⟨ O^J/ψ[n] ⟩ represents the LDME of the state indicated by n. 
The numerator in eq:AN is directly sensitive to the Sivers function and within the GPM reads


    Δ̣σ^GPM   =   1/(4π)^4 zs∑_a∫x̣_a/x_a ^̣2 k_⊥ a δ(ŝ+t̂+û-M_^2+Q^2) sin(ϕ^*_S-ϕ^*_a)
          ×∑_n1/Q^4 (-2  k_⊥ a/M_p) f_1T^⊥ a (x_a, k_⊥ a)  L^μν
    H^a,U_μν[n] ⟨ O^J/ψ[n] ⟩ ,

where f_1T^⊥ a(x_a, k_⊥ a ) is the Sivers function.


The numerator of the asymmetry in the CGI-GPM is 
given by 


    Δ̣σ^CGI   =   1/2s2/(4π)^4 z∫x̣_a/x_a ^̣2 k_⊥ a δ(ŝ+t̂+û-M_^2+Q^2)sin(ϕ^*_S-ϕ^*_a) (-2  k_⊥ a/M_p) 
          ×∑_n1/Q^4 L^μν{∑_q f_1T^⊥ q (x_a, k_⊥ a) H^q,Inc_μν[n]+
    
    
    
    f_1T^⊥ g(f) (x_a, k_⊥ a) 
    H^g,Inc(f)_μν[n]}⟨ O^J/ψ[n] ⟩ ,



where H^a,Inc_μν[n] is the perturbative square amplitude calculated by incorporating the FSIs within the CGI-GPM approach. 
Note that, in num:CGI-GPM, there is no contribution from the d-type GSF. In fact, in ep collisions, ISIs are absent due to the colourless nature of the virtual photon and only the f-type GSF is contributing to the Sivers asymmetry <cit.>. This means that quarkonium production in ep collisions is a powerful tool to directly access the process-dependent f-type GSF. Moreover, the modified colour factor associated with the ^3S_1^[1] state is zero in the CGI-GPM approach, which leads to a vanishing Sivers asymmetry in the 
CSM.

By adopting a Gaussian factorised form for the unpolarised TMD distribution, a Gaussian-like Sivers distribution 
and by maximising the latter we can give estimates for the upper bounds of the Sivers asymmetry (eq:AN) A^sin(ϕ^*_h-ϕ^*_S)_N at the EIC. 
Results are presented in fig:EIC45ptz_BK, and are computed using the following kinematical cuts: 2.5 GeV^2 < Q^2 < 100 GeV^2,  10 GeV< W_γ p < 40 GeV, 0.3 < z < 0.9 and < 5 GeV. The BK11 LDMEs set <cit.> is adopted.

 


The asymmetry is mostly dominated by the GSF, while the quark contribution is negligible. This indicates that such an observable is a powerful tool to probe the unknown GSF. The GPM predicts negative values around 20%. The asymmetry is drastically reduced in size in the CGI-GPM due to colour-factor relative cancellations and the absence of the ^3S_1^[1]-state contribution and is essentially driven by the f-type GSF.











 §.§ Generalised Parton Distributions





Information on the three-dimensional structure of the nucleon, correlating the transverse position of partons with their longitudinal momentum, 
is provided by 
GPDs. Processes to access GPDs include Deeply Virtual Compton Scattering (DVCS) and Deeply Virtual Meson Production (DVMP). 


A factorisation theorem has been proven for DVCS in the Bjorken limit <cit.>. It allows one to compute the DVCS amplitude as the product of some GPDs and corresponding coefficient functions that can be calculated perturbatively. GPDs are in very solid theoretical footing: at leading-twist level, all-order QCD-factorisation theorems directly relate the GPDs to particular hard exclusive scattering processes. GPDs are thus process-independent, universal quantities.
In the case of DVMP, factorisation applies in the case of longitudinally polarised photons. The hard-scattering process includes the exchange of hard quarks and gluons, involving the strong coupling constant α_s and a meson distribution amplitude, which is not completely understood to date.

The GPDs do not uphold a probabilistic interpretation like PDFs 
do, but are well-defined in quantum field theory as matrix elements of bilocal quark and gluon operators at a light-like separation. In the light-cone gauge at leading twist, the quark GPD is

    F^q(x,ξ,t)    = 1/2∫d z^-/2 π e^i x P^+ z^-⟨ p' | ψ̅^q (-z/2) γ^+ ψ^q (z/2) | p ⟩ |_z^+ = z_⊥ = 0
       =1/2 P^+[ H^q (x,ξ,t) u̅(p') γ^+ u(p) + E^q(x,ξ,t) u̅(p') iσ^+ μΔ_μ/2m_N u(p) ]

and the gluon GPD,

    F^g(x,ξ,t)    = 1/P^+∫d z^-/2 π e^i x P^+ z^-⟨ p' | F^+ μ(-z/2) F^+_μ(z/2) | p ⟩ |_z^+ = z_⊥ = 0
       =1/2 P^+[ H^g (x,ξ,t) u̅(p') γ^+ u(p) + E^g(x,ξ,t) u̅(p') iσ^+ μΔ_μ/2m_N u( p) ],

where z= (z^+, z_⊥, z^-) are the light-cone coordinates, P^+ is the light-cone plus-component of the average of the incoming- and outgoing-nucleon momenta, x is he fractional parton plus-component momentum of the nucleon, ξ the skewness variable and t the Mandelstam variable, which represents the four-momentum transfer squared to the nucleon.  The symbols γ and σ are the Dirac matrices,   
u and u̅ are nucleon spinors and m_N is the mass of the nucleon. Here, F^q and F^g are both expressed as a Fourier transform of a matrix element of a chiral-even operator formed from either quark fields ψ^q or the gluon-field strength tensor F^μ ν. The result is a decomposition into twist-2 parton-helicity conserving GPDs H and E.

GPDs cannot be directly extracted from experimental data. Indeed, in the expression of the cross section of exclusive electroproduction processes, GPDs appear in convolution integrals known as Compton Form Factors (CFFs). These CFFs are complex quantities, the real and imaginary parts of which provide complementary constraints on GPDs. The DVCS 
CFF ℋ, at leading-twist and leading-order (and at fixed momentum transfer t and skewness ξ), for example, is given by 

    ℋ = ∫_-1^1 dx F^q(x,ξ,t)/x-ξ + iϵ = 𝒫∫_-1^1 dx F^q(x,ξ,t)/x-ξ - iπ F^q(±ξ, ξ, t),

and with 

    σ(γ^* p →γ p) ∝ |ℋ|^2.



In addition, there are also spin-dependent GPDs

and are probed in measurements in which the spin or polarisation state is fully defined. If the spin states are averaged over, as in the description of an unpolarised measurement, then there is no way to have a direct dependence on, or be sensitive to, these objects. Moreover, there are also parton-helicity-flip GPDs (chiral odd), in which the initial- and final-state hadrons have different polarisations.

GPDs are also connected to the distribution of pressure and shear forces inside the nucleon <cit.> and, furthermore, the second moment of a particular combination of GPDs is related to the 

angular momenta of quarks and gluons via Ji's relation <cit.>. A comprehensive review on the phenomenology of GPDs in DVCS can be found in <cit.>.

















  §.§.§ Gluons





DVCS is sensitive to quarks and, at higher order and/or higher twist, also to gluons. On the other hand, the production of light mesons in DVMP probes quarks and gluons, depending on the energy scale at which the process is measured. However, production in exclusive photoproduction (or electroproduction) reactions is a golden channel for gluon GPDs. Indeed, in this case the quark exchange plays only a minor role and due to the large scale provided by the heavy-quark mass, perturbative calculations is expected to be applicable even for photoproduction <cit.>. 


At the EIC, precise measurements of exclusive cross sections will be possible in order to map out the dependence on the squared momentum transfer to the nucleon 
t=(P_N-P_N')^2  for , ϕ and K, among others. EIC will cover the region of 0 <|t| < 1.5 GeV^2, down to an impact parameter of ∼ 0.1 fm.



Figure <ref> shows the projected precision obtainable at the EIC in the exclusive electroproduction cross section as a function of the momentum transfer t to the proton, for different bins in x_=(Q^2+M_^2)/(2 p· q), 


the x-Bjorken equivalent scale variable for heavy mesons. The projections are produced using the LAGER <cit.> event generator and are based on the calculations presented in <cit.>. LAGER is described as a modular accept-reject generator, capable of simulating both fixed-target and collider kinematics, and has previously been used for vector-meson studies at EIC kinematics, with significant recent developmental effort in support of DVMP studies.
The transverse spatial distribution of partons can be obtained by a Fourier-transform
of the cross section as a function of t.

The key experimental feature of hard exclusive channels such as electroproduction is the detection of the recoil protons in the far-forward detectors, in particular in the  B0 spectrometer and the Roman Pots. This allows 
for accurate computation of the momentum transfer t, which is the Fourier conjugate variable to the impact parameter. A wide and continuous acceptance that extends to low-t is essential for a precision extraction of transverse-position distributions of partons.



On the other hand,  far-forward detectors can also help in detecting 
the process where the proton does not stay intact but breaks up. The dominance of this process over exclusive production increases with increasing t. In <cit.>, 
it has been shown that the cross-section measurement of dissociative diffractive photoproduction at large t as a function of the rapidity gap between the produced and the dissociated proton is possible at the EIC. The interest of this process lies in the presence of two comparable hard scales, the charm mass and the large t, and hence the possibility to probe the presence of Balitsky-Fadin-Kuraev-Lipatov (BFKL) dynamics.



  §.§.§ Light quarks



In <cit.>, it was shown that the rapidity differential cross section for exclusive photoproduction in heavy-ion ultra-peripheral collisions at NLO decomposes into a complicated interplay of contributions from both the quark and gluon sectors as well as their interference, over the whole region of rapidities accessible at the LHC. In particular, at mid-rapidities the quark contribution was shown to be the dominant player. While such a picture remains in place under a conservative factorisation and renormalisation scale variation, and is reflected in the original work of Ivanov et al. <cit.> in the context of the underlying hard scattering process, γ p → p, which drives the ultra-peripheral collisions, and indeed the eA collisions at the EIC, care must be taken to interpret such results. Indeed, it was shown that such a hierarchy arises from a coincidental cancellation of LO and NLO gluon contributions together with the positive-definite quark contribution at NLO. At NNLO, when there are also interference contributions wholly within the quark sector, one may anticipate a different final picture. Υ photoproduction on the other hand, sitting at a higher scale, does not exhibit such a complicated interplay of contributions at NLO, see <cit.>, with the gluon contribution dominating over all rapidities. The results are therefore indicative of the long-standing problem of the scale dependence and perturbative instability exhibited by low-scale processes. Indeed, after the so-called `Q_0 subtraction' <cit.> discussed in Sect. <ref>, the quark contribution to the amplitude becomes negligible. A new study <cit.> which includes the high-energy resummation effects in the coefficient function of exclusive photoproduction in the HEF formulism similar to one applied in the inclusive case <cit.> supports this conclusion.






 §.§ Generalised TMDs






The non-perturbative structure of the hadrons can be described in terms of parton correlation functions such as form factors, 1D PDFs and their 3D generalisations in terms of TMDs and GPDs. All these functions can be derived from more general objects called GTMDs <cit.>. Hence, GTMDs are also known as the “mother distributions". There are several compelling reasons to study GTMDs
. Firstly, GTMDs contain physics that outmatches the content encoded in the TMDs and GPDs.  Secondly, via Fourier transformation, GTMDs can be related to Wigner functions, a concept that spans across other branches of physics as well. Partonic Wigner functions may allow for a hadron tomography in 5D phase-space <cit.>. Thirdly, certain GTMDs can unravel unique correlations between parton orbital motion and spin inside hadrons <cit.>. In particular, the Wigner distribution can be used for a gauge-invariant definition of the canonical orbital angular momentum <cit.>, which makes this quantity also accessible for calculations in lattice QCD <cit.>.  Fourthly, there is a particular GTMD that is related to the Sivers TMD. By establishing a relation between GTMDs and the QCD odderon at small x, the authors in Ref. <cit.> have shown that one can access the gluon Sivers TMD through exclusive π^0 production in unpolarised ep scattering. This finding goes against our traditional belief that the Sivers function can only be measured with a transversely polarised target.





 
For a long time, it was questionable whether GTMDs could be measured at all. The authors in Ref. <cit.> were the first to propose addressing gluon GTMDs through exclusive diffractive dijet production in lepton-nucleon/nucleus collisions at small x (see left panel of fig:exclusive-dijet). The GTMDs depend on the average transverse parton momentum k⃗_⊥ and the transverse momentum transfer to the target Δ⃗_⊥, and it is possible to decompose the angular correlation between these two vectors into a Fourier series. The leading angular dependent term, known as the elliptic distribution, has a characteristic cos(2ϕ) angular modulation similar to the observed elliptic flow phenomenon in relativistic heavy-ion collisions <cit.>. It was shown that the cross section of this diffractive dijet process also exhibits such a cos(2ϕ) behavior where ϕ is now the angle between the dijet total and relative momenta. The pioneering work in Ref. <cit.> gave impetus to the field of GTMDs and subsequently many other interesting ideas were put forward; see, for instance, Refs. <cit.>.

An alternative idea <cit.> is to exclusively produce a single particle (instead of two jets) such as a J/ψ. The role of the second jet is now played by the scattered electron which must be detected. It has been shown that in this process the elliptic cos 2ϕ correlation of the gluon GTMD manifests itself in the angular correlation between the  scattered electron and the  J/ψ <cit.> (or the recoiling proton/nucleus <cit.>). For a proton target, a sizable v_2 of a few percent or larger has been predicted <cit.>. The same effect can also be seen in DVCS, but J/ψ production is more promising since there is no contamination from the Bethe-Heitler 
process. In the GPD-based approach to DVCS, the same angular correlation is known to be generated by the so-called gluon transversity GPD. The elliptic gluon GTMD is the mother distribution of the gluon transversity GPD <cit.>.

Quarkonium production processes are also useful to study other aspects of GTMDs. 
In Ref. <cit.>, it was shown that exclusive double production of pseudo-scalar quarkonia (η_c/b) in hadronic collisions could serve as a direct probe of GTMDs for gluons at moderate x (see right panel of fig:exclusive-dijet).  A similar idea came out in Ref. <cit.> where the authors proposed to access  the Weizsäcker-Williams gluon GTMD at small x via double  χ_cJ or η_c meson production in 
diffractive pp/pA collisions where (one of) the proton(s) stays intact.  


At the EIC, the primary process to look for gluon GTMDs is exclusive diffractive dijet production, as mentioned above. A challenge, however, is that due to the limited 
centre-of-mass energy, the transverse momenta of diffractively produced particles in the forward rapidity region are often not large enough to cleanly reconstruct jets. 
As a first step to test the underlying GTMD picture of exclusive diffractive production processes, like dijet or J/ψ electro- and photoproduction at small x, a GTMD model can be fitted to existing HERA data. Predictions can then be obtained for EIC in different kinematic regions. This has been considered for dijet production in <cit.>, where it was shown that a gluon GTMD model based on the impact-parameter-dependent McLerran-Venugopalan model can give a reasonably good description of diffractive dijet production data from 
H1 <cit.>. 


The same framework (slightly extended) can be applied to exclusive diffractive J/ψ production to describe the H1 and ZEUS data, as shown in fig:JPsi on the left (= 319 GeV). With the resulting GTMD parametrisation, predictions for exclusive  diffractive J/ψ production at EIC can be obtained. These are shown for = 45 and 140 GeV in fig:JPsi on the right. 
 Generally, at small x, and in particular for nuclear targets, a GTMD-based description becomes more appropriate for exclusive and diffractive processes.
 
 Exclusive quarkonium production at the EIC could be used to systematically study the transition between the collinear and k_⊥-dependent frameworks.




 




 §.§ Exclusive quarkonium production near threshold and the trace anomaly








It has been noticed long ago that the mass M of a hadronic system can be expressed in terms of the forward matrix element of the trace of the QCD energy-momentum tensor 
as <cit.>

    2M^2=⟨ p| β/2g F^2+(1+γ_m) ψmψ |p⟩,

where β and γ_m are anomalous dimensions and the operator β / (2g) F^2 + γ_m  ψmψ is the QCD trace anomaly <cit.>.  The decomposition of the r.h.s. of Eq. (<ref>) 
into quark and gluon contributions has been discussed in detail <cit.>. Other mass decompositions, based this time on the QCD Hamiltonian, have also been proposed in the literature <cit.>. The latter all require the knowledge of the same four quantities, combined in different ways for the physical interpretation <cit.>. Two of these quantities, namely the quark momentum fraction A_q(0)=⟨ x⟩_q and the gluon momentum fraction A_g(0)=⟨ x⟩_g, are already well known. The other two numbers C̅_q(0) and C̅_g(0) can be determined by measuring the quark and gluon contributions to Eq. (<ref>). While the quark condensate ⟨ p| ψmψ |p⟩ has already received a lot of attention over the last decades (see <cit.> and references therein), little is known so far about the gluon condensate ⟨ p | F^2 | p ⟩ from the experimental side. 

Four-momentum conservation implies that A_q(0)+A_g(0)=1 and C̅_q(0)+C̅_g(0)=0. From a phenomenological point of view, the knowledge of A_q(0) and the quark condensate  is therefore sufficient for specifying all the contributions to the various mass decompositions (see <cit.> for recent estimates). Measuring the gluon condensate is not expected to change much the current phenomenology of the nucleon mass, but it will provide a fundamental sanity check of the mass sum rules and the virial theorem <cit.>. Another motivation for measuring the gluon condensate is that it could shed light on the existence and nature of the recently discovered LHCb “pentaquark” states <cit.>. 

More than two decades ago, exclusive heavy-quarkonium production, near the production threshold, was suggested as a promising tool for constraining the gluon condensate in the nucleon <cit.>.  
This development together with the prospect to obtain through this process further information about the gravitational structure of the nucleon, which is contained in the form factors of the energy-momentum tensor 

(such as the mass radius and mechanical pressure distributions <cit.>), as well as the measurement of exclusive J/ψ photoproduction near threshold at Jefferson Lab <cit.> has stimulated a significant amount of activities in this area <cit.>.
Recently, it was argued that the extraction of the gravitational form factors through exclusive quarkonium photoproduction will necessarily retain model dependence <cit.>. Generally, access to the gravitational structure of the nucleon is expected to be cleaner for electroproduction <cit.>. 
At the EIC, one would have the unique opportunity to explore photo- and electroproduction of both J/ψ and Υ close to threshold  <cit.>.




 §.§ Probing double parton scattering at the EIC with quarkonium pairs


 



  §.§.§ A word of context

In this section, we study the possibility to observe double-J/ψ production at the EIC. In particular, we discuss both the single-parton-scattering (SPS)  and the double-parton-scattering (DPS) mechanisms, which could lead to the observation of a pair of . In fact, the cross section for the latter case would allow one to access new information on the so-called proton 
double-parton-distribution functions (dPDFs), which encode novel information on the partonic  structure of the proton.  



Let us recall the analysis of four-jet photoproduction  at HERA, which pointed out the relevance of  multi-parton interactions (MPIs) to account for the measured total cross section <cit.>.

In  Ref. <cit.>,  the DPS cross section for four-jet photoproduction was 
calculated. 
DPS are initiated by a quasi-real photon <cit.> splitting into a  q q̅ pair. The same strategy as for pp collisions <cit.> has been used to evaluate the photoproduction cross section. At this stage, the only missing quantity was σ_eff^γ p, the effective size of the photon-proton interaction, which is expected to be  process independent. It was estimated 
for the first time <cit.> and compared to that of the pp case from Refs. <cit.>.The four-jet DPS cross section has then been calculated for the HERA kinematics <cit.> to be σ_DPS^4j≥ 30 pb, while the total one was inferred from <cit.> to be  σ_tot^4j∼ 135 pb at x_γ < 0.75. This indicated that the DPS contribution  is  sizeable even in photon-induced reactions for the production of four jets 
and that it could also be so for other processes like quarkonium-pair production. Further analyses of the HERA data could lead to the extraction of σ_eff^γ p and, in turn, provide a first access to the mean transverse distance between two partons in the proton, 
an unknown property of the proton structure. To this aim, the needed luminosity was evaluated to be ℒ∼ 200 pb^-1  <cit.>. Double-J/ψ production from DPS at EIC will be 
presented below along the same lines.








  §.§.§ DPS at the EIC and -pair production 




Here we discuss J/ψ-pair photoproduction at the EIC.  In ep collisions, the radiated quasi-real photon can interact with the partons 
within the proton in two ways, namely as a “pointlike” particle and via its “resolved” hadronic content. In the first case, the photon “directly” interacts with the target while, in the latter case, the  photon splits into (colour charged) partons, which subsequently interact with partons in the proton. 



The treatment of the interaction between a proton and such a resolved photon is 
carried out by using a PDF describing the momentum distributions of these partons inside the photon. One of these is the  GRV <cit.> set, which is 
adopted here.


For what regards the quarkonium-production mechanisms, the CSM ( the leading v^2 contribution of NRQCD) is used. 

fig:DPS_SPS shows different Feynman graphs for SPS and DPS photoproduction. 
In the SPS case, the contributing channels at leading order,  αα_s^4, are shown in fig:DPS_SPS(a-c), namely, γ q→ J/ψ+J/ψ+q, g g→ J/ψ+J/ψ and q q̅→ J/ψ+J/ψ. However, the graph in fig:DPS_SPS(d) contributes at order αα_s^5, i.e. via the SPS γ g→ J/ψ+J/ψ+g+g. The gluon-initiated channel in DPS for di-J/ψ production at α_s^6 is shown in fig:DPS_SPS(e), while the quark-initiated channel does not contribute in the CSM at order α_s^6. The partonic channel 
gg→ J/ψ+g 
dominates for single-J/ψ production.
The SPS cross section,  the squared matrix elements convoluted with single-parton PDFs, can be calculated using 
HELAC-Onia
 <cit.>. In order to estimate the DPS cross section, we need to use the poorly known proton dPDFs, which

provide the number densities of a parton pair with a given transverse distance b_⊥ and carrying  the longitudinal momentum fractions (x_1,x_2) of the parent hadron <cit.>. 
Assuming that dPDFs can be factorised in terms of 
ordinary 1D PDFs and a transverse part, the DPS cross section can be expressed in terms of two SPS cross sections for the production of each 
of the observed particles among the pair:

    σ^(J/ψ ,J/ψ)_DPS=1/2σ^(J/ψ)_SPS σ^(J/ψ)_SPS/σ_eff^γ p ,

which is the so-called "DPS pocket formula", valid under the assumption of 

totally uncorrelated 

kinematics between both parton scatterings. 
The σ^(J/ψ)_SPS is the SPS contribution for single J/ψ production.  
In the present study, within the mentioned assumptions, one gets:

    σ_eff^γ p=  [∫ d^2k⃗_⊥(2 π)^2 F_2^γ(k⃗_⊥,Q^2) F_2^p(k⃗_⊥) ]^-1

where here F^p(γ)_2(k_⊥) parametrises the transverse structure of the proton (photon) <cit.>. For the photon, the only available calculation is that of Ref. <cit.> while, for the proton, there are several models based on the data for 
DPS in pp collisions.

Recently, several experimental analyses on DPS  have been carried out  for the production of J/ψ + W <cit.>, J/ψ+Z <cit.>, J/ψ+charm <cit.>  in pp and J/ψ+J/ψ <cit.> in pp̅ 
processes. A comprehensive comparison between theory and experiments for di-J/ψ production at the Tevatron and the LHC has been presented in <cit.>, and it 
was observed that DPS dominates the yield at large J/ψ-rapidity difference. 
DPS has been also studied for J/ψ-pair production for the LHC fixed-target (also referred to as AFTER@LHC) kinematics in <cit.>.

At the EIC, LO computations using HELAC-Onia show that measurements are possible at = 140 GeV with SPS contributions generally dominant over the DPS ones, but there are certain regions (low z and large Δ y) in the phase space where DPS cannot be disregarded. If σ_eff^γ p is not too small, DPS events could be measured. In these regions, there is thus a compelling opportunity to distinguish between the resolved and unresolved contributions in the cross section and thereby to gain valuable insight into the internal structure of photons and protons.









 



§ QUARKONIA AS TOOLS TO STUDY THE PARTON CONTENT OF NUCLEI










 §.§ Nuclear PDFs




Decades of experimental and theoretical studies showed that the distributions of partons in a nucleus are considerably modified compared to the nucleon ones. While significant progress has been made since the initial observation of the modification of PDFs in bound nucleons by the EMC Collaboration <cit.>, our understanding of nuclear PDFs (nPDFs) is still not satisfactory, most notably in the case of gluons. Measurements of quarkonium production in 
eA reactions can bridge this knowledge gap.

One of the main EIC goals is a high-precision survey of the 
partonic structure of the nucleus 
to significantly advance 
our quantitative understanding of 
nPDFs. The EIC 
will offer the possibility 
to study nPDFs over a broad range of momentum transfers


 <cit.>. An improved knowledge of nPDFs will enable more precise theoretical calculations for 
nuclear effects and increase the scientific benefit of already successful heavy-ion programmes at RHIC and LHC.     


A widely accepted approach to quantify nuclear effects in 
PDFs is to start with proton PDFs and use a function R(x,Q^2) that captures the modification of a given PDF in a nucleus. Experimentally, such a modification could be studied by a ratio of structure functions F_2 
or by the so-called nuclear modification factor as done by RHIC and LHC experiments. In the case of eA collisions, R(x,Q^2) is defined as  





    R_eA = 1/A(d)σ_eA/(d)σ_ep,

where (d)σ_eA and (d)σ_ep are the cross sections for 
the process under consideration, respectively, in eA and ep reactions, while mass number A serves as a normalisation factor. Note that these cross sections can be differential in different kinematical variables.  With the definition of eq:R_eA, = 1 in the absence of nuclear effects.  In the following, we review and quantify prospects for nuclear-PDF determination at EIC via   measurements.











  §.§.§ Gluons




In order to give an 
estimate of the potential impact of the EIC on nPDF determination, 
the nuclear modification factor 
R_eAu, 

which can be measured in inclusive J/ψ photoproduction in 
eAu reactions, 
is compared with projected statistical uncertainties. Such a prediction 
is shown in fig:R_eAu-45GeV and fig:R_eAu-90GeV at two different values of the centre-of-mass energy
, √(s_eN),  
45 GeV and 90 GeV, as a function of the J/ψ rapidity in the Nγ centre-of-mass frame
[Note that we adopt the same kinematical configuration as the EIC Yellow Report, with the proton(ion) moving along +ẑ and the electron along -ẑ (see also Fig. <ref>).] and as a function of W_γ N. Kinematical cuts are applied on the 
elasticity (0.2 < z < 0.9) and on the pseudorapidity of the electron pair coming from the J/ψ→ e^+e^- decay (|η_ee| < 3.5). Different cuts on W_γ N are applied for the rapidity spectra at the two different  energies. 



The nuclear-modification-factor predictions are calculated using HELAC-Onia <cit.>, adopting the CT14nlo set <cit.> as a proton PDF baseline and using two different nuclear PDF sets for the gold nucleus, namely EPPS16nlo <cit.> and nCTEQ15FullNuc <cit.>. Factorisation and renormalisation scales are taken to be the J/ψ transverse mass, μ_F = μ_R = m_T=√(M^2_J/ψ+^2). Note also that, since these predictions are calculated at LO in the CSM, where the only partonic subprocess is γ + g → J/ψ + g, they can be directly interpreted as R_g, the nuclear modification factor for the gluon nPDF. The statistical projections are calculated assuming 
R_eAu=1 (using the central value of CT14nlo) and assuming an integrated luminosity of 10 fb^-1/A. 
The branching ratio for the J/ψ→ e^+e^- decay was taken to be 5.94% and a  reconstruction efficiency of 64% was assumed (considering an average identification efficiency of the electrons from the  decay to be approximately 80%).    





Some comments are in order. First, as  can be seen in fig:R_eAu-45GeV and fig:R_eAu-90GeV,  
J/ψ is expected to be mostly produced in the backward region in the 
γ N centre-of-mass frame 
as the yield essentially vanishes at positive rapidities (see the increase of the statistical uncertainties of our projections). This happens for both energy configurations.  Second, 
the regions where shadowing (relative parton depletion at x smaller than 0.01), antishadowing  (relative parton excess at x around 0.11) and the EMC effect (relative parton depletion for 0.3 < x < 0.7) take place 
can be probed at the 
EIC via J/ψ photoproduction. The antishadowing peak is expected to be observed at moderate backward rapidity in the 
γ N centre-of-mass frame, while the shadowing region would be probed at larger negative rapidities. Such regions are also 
those where the projections point to a smaller statistical uncertainty compared to the PDF and scale uncertainties, 
 i.e. the gluon nPDFs would be the most constrained. 
 The W_γ N dependence of the nuclear modification factor would also be a very interesting tool to probe gluon nPDFs. A large shadowing tail is expected to be probed for larger values of W_γ N, while clear antishadowing peaks are expected in the region W_γ N∈ [10:20] GeV, in both energy configurations. The projected uncertainties are also 
 small, and seem to have an interesting constraining power for the gluon nPDFs. More detailed dedicated studies are surely required and would help in motivating new measurements to probe gluon nPDFs at the 
 EIC.





fig:jpsi-pt-eAu-eic presents predictions for 
the dependence of at =100 GeV 
by using the same factorisation formalism in eq:jpsi-lp-fac, with proton PDFs replaced by nuclear PDFs for the eA collision. The total, LP and NLP contributions are shown.  
The EPPS21nlo central set <cit.> is used as nPDF. 
Since the production rate is dominated by the γ+g→ [cc̅] + g subprocess, this ratio is directly sensitive to the nuclear dependence of the gluon PDF. At EIC energies, 
the distribution of production is sensitive to the gluon at a relatively large momentum fraction due to the soft-photon distribution in 
the incoming electron. The enhancement of the production rate in 
eAu over ep collisions in fig:jpsi-pt-eAu-eic is a direct consequence of the “antishadowing" behavior of the nuclear gluon distribution from the EPPS21nlo nuclear PDF set.  Since the quark-initiated subprocesses dominate the LP contribution, the ratio of the LP contribution (blue dashed and red dotted lines) shows the well-known EMC-type effect from nuclear quark PDFs. However, this feature of the LP contribution does not have a real impact on the observed nuclear dependence of the dependence of production at the EIC energies (the solid line), since the LP contribution is strongly suppressed; that is, the distribution of production at the EIC should also be an excellent observable for probing the nuclear gluon PDF.



























 §.§ Nuclear GPDs





In coherent diffractive production of vector mesons off a nucleus, the light (photon) generated by the electron interacts, similarly to optical experiments of diffraction, with the nucleus as a whole, resulting in the production of a vector meson in the final state. 
This process has been proposed as a tool to investigate gluon saturation dynamics <cit.>. Here, the production of lighter vector mesons, such as the ϕ meson, is expected to be sensitive to saturation effects. On the other hand, the production of quarkonia would because of the heavier quarkonium mass (and thus smaller size of the dipole formed by the quark–anti-quark pair that evolves into the vector meson) not be optimal to study gluon saturation and rather serve as a baseline free from saturation effects. Diffractive production also gives access to the spatial distribution of partons inside the nucleus. While coherent diffractive production provides information on the average spatial distribution of partons, incoherent production, where the nucleus does not stay intact, probes local fluctuations of this spatial distribution <cit.>.
For the study of the spatial distribution of gluons in heavy nuclei, in particular, the diffractive production of a quarkonium, such as a J/ψ, is most adequate. For the coherent process, the momentum transfer distribution 
√(|t|) from the photon to the target nucleus is expected to exhibit a diffractive pattern, where the details of the shape of this pattern encode information on the gluon GPD <cit.>. An example of such a diffractive pattern is shown in Fig. <ref>, as represented by the square symbols. The data points have been simulated using the Sartre Monte-Carlo event generator <cit.>. Results including (filled symbols)  and excluding (open symbols) saturation effects are shown. In addition to the diffractive coherent production, the expected incoherent contribution (circles) is shown.
As can be seen, apart from the very low |t| region, the incoherent contribution dominates
the coherent one.



















Elastic and inelastic diffractive quarkonium production off the proton has been studied at the HERA lepton-proton collider experiments H1 <cit.> and  ZEUS <cit.>, while a first measurement of exclusive J/ψ photoproduction at threshold has been performed in the fixed-target experiment GlueX at Jefferson Lab <cit.>.
At hadron-collider experiments, diffractive quarkonium production has been investigated in collisions <cit.> at the Tevatron, in  <cit.>,  <cit.> and  <cit.> collisions at the LHC and in 
dAu <cit.> and 
AuAu <cit.> collisions at RHIC.

The existing measurements off nuclei are at present restricted in statistical precision, while only offering a rough determination of the momentum transfer  √(|t|) and in general a limited separation of coherent and incoherent production. Hence, the knowledge on the gluonic structure of nuclei is at present poor, with many fundamental questions unanswered.













The EIC is expected to perform measurements of diffractive vector-meson production off light and nuclear ions with unprecedented precision. The two experimental challenges consist in determining t with high precision and in
distinguishing coherent from incoherent events <cit.>. Recently, the capability of proposed EIC detectors in reconstructing t and
their ability to suppress incoherent production have been examined <cit.>, <cit.>, <cit.>.
The variable t needs to be reconstructed from the scattered lepton and reconstructed vector meson, since in coherent production the trajectory of the ion after the interaction is nearly unmodified and thus the ion cannot be detected, while in the case of incoherent production not all fragments from the nuclear break up can be detected. The distribution in |t| for coherent diffractive J/ψ production off gold ions is shown in  Fig. <ref>, left.
Here, |t| is reconstructed  as the squared sum of the transverse momenta of the scattered lepton and of the lepton pair originating from the J/ψ decay.
It forms a good approximation for the true -t.

The data have been simulated again with Sartre and subsequently passed through a full simulation of the ePIC detector.
The histogram represented by the continuous line is the generated distribution, while the other curves represent the reconstructed distribution, with beam effects. 
The latter include an angular divergence originating from the focussing and defocussing quadrupoles in the interaction region and
a small angular kick from the crab cavities. The crossing angle from the beams in principle also influences the t distribution, but contrary
to the other effects it can be corrected for.
For the curve indicated by the open, blue circles only information from tracking detectors is used for the reconstruction of the scattered lepton, while for the curve indicated by the black, closed circles only information from the backward electromagnetic calorimeter is used for the reconstruction of the scattered lepton. The curve indicated by the red, open circles selects the best of the two methods. As can be seen, the quality of the reconstruction in t is strongly dependent on the quality of the reconstruction of the scattered beam lepton. In the diffractive process the beam lepton generally is scattered under a small angle and covers a region where the tracking performance is degraded. Using in addition the electromagnetic calorimeter in the backward region for the reconstruction of the scattered lepton improves the reconstruction in t vastly.

The spatial distribution of partons in impact-parameter space is
related to a Fourier transformation, with t going from 0 to infinity <cit.>. Experimentally, one is limited by a maximal momentum transfer, which preferably extends as far as possible. In practice,
studies have shown that it is necessary to resolve the minima up to the third one for the evaluation of the spatial distribution <cit.>. This dictates the needed level of suppression
of the incoherent contribution. The suppression of incoherent events includes the requirement of exactly three reconstructed lepton tracks with the correct charge in absence of any other signal in the main detector and various criteria corresponding to the absence of signal in a series of far-forward detectors, which can tag protons (Roman Pots for protons with energy close to the beam energy and the B0 spectrometer and off-momentum detectors for nuclear-breakup protons), neutrons (Zero-Degree Calorimeters) and photons (B0 and Zero-Degree Calorimeters). The capability to suppress incoherent production is illustrated in Fig. <ref>, right, which shows the -t distribution for coherent and incoherent production off gold nuclei.
The former is again simulated using Sartre, while for the latter the BeAGLE generator <cit.> is used. 




The generated coherent (incoherent) contribution is represented by the continuous (dotted) line. The generated data are passed through 
a full simulation of the ATHENA detector. The effect of data selection requirements  on the event activity in the main detector and on the absence of activity in the far-forward detectors,
based on the studies in Ref. <cit.>, is represented by the
blue, open circles. As can been seen, the obtained distribution lies close to the distribution from coherent events simulated by Sartre. The remaining contribution from incoherent events is given by the red, star symbols. The largest suppression of the incoherent process comes from the requirement on the absence of any neutron signal in the Zero-Degree Calorimeter, while the requirement on the absence of photon signals in this Zero-Degree Calorimeter 
also has an impact. 

Ways to further improve the reconstruction of t and the suppression of incoherent production are at present under investigation.













The study of light nuclei can offer additional insights into the 
internal structure of the nuclear medium.
In contrast to 
measurements with heavy nuclei, the total final state in incoherent diffractive production off light nuclei can be
unambiguously identified through tagging of the spectator nucleons.
Such measurements are of interest when studying the short-range correlation (SRC) of a nucleon pair, which is the temporal fluctuation of two nucleons into a strongly interacting pair in close proximity and large measured relative momentum <cit.>.
SRC pairs are suggested as a possible explanation for the nuclear modification of the momentum distribution of high-x partons, known as the EMC effect,
with a strong correlation between the two phenomena suggested by measurements by the CLAS experiment at Jefferson Lab <cit.> and a quark-level QCD basis for SRC has been proposed for the lightest nuclei <cit.> and A≥4 nuclei <cit.>.

The simplest nuclear system consists of deuteron and the first measurement of incoherent diffractive production with spectator tagging was performed in the measurement of incoherent diffractive J/ψ production in ultra-peripheral 
dAu collisions by the STAR experiment at RHIC <cit.>, with tagging of the spectator neutron in the Zero-Degree Calorimeter.
At the EIC, similar measurements can be performed with enhanced precision, and studies of incoherent diffractive J/ψ production off the deuteron at the EIC have been proposed to study the nuclear modification of the
gluon distribution and its possible link with the SRC <cit.>.
For the proposed measurement, the scattered lepton and J/ψ decay leptons are reconstructed in the main detector, while both the leading
and spectator nucleon (neutron and proton) can be detected in the far-forward detectors. The detection of both nucleons instead of only one 
offers certain advantages in the reconstruction of the event and some
kinematic variables <cit.>.

In Fig. <ref>, the three-momentum distribution of the tagged neutron (left) and tagged proton (right) in the deuterium rest frame is illustrated for incoherent diffractive production of J/ψ in the scattering of 18 GeV electrons off 110 GeV deuterons at the EIC, as
simulated with BeAGLE <cit.>.
The star symbols represent the generated distribution, the open circles represent the distribution including acceptance effects of the main and far-forward detectors, and the open squares also take 
the finite
detector resolution and beam effects into account. The momentum distribution of the tagged nucleon reflects the initial-state momentum
of the nucleons inside the deuteron.
The region above 300 MeV corresponds to the region of the SRC, and as visible in the figures, the EIC will be able to provide a good reconstruction of the tagged-nucleon momentum. A similar statement holds for the reconstruction of
other variables of interest <cit.>.

































 §.§ Study of transport properties of nuclear matter




The vital element of portraying nuclear matter is to get information on how the medium responds to a parton traversing the matter. It is characterised by transport coefficients, 
  a diffusion coefficient 
or 
q̂, which is the mean squared momentum transfer between the propagating particle and the medium per unit length. Transport coefficients are an essential ingredient in the modelling of nuclear reactions 
, and determining these parameters is one of the main goals of high-energy nuclear physics experimental and phenomenological efforts.

Measurements of hadron production in  collisions have shown a broadening of the transverse momentum distribution at 
intermediate 
hadron transverse momentum compared to  reactions. This phenomenon is visible over a wide range of hadronic collision energies, starting from collisions at ≈ 20 GeV <cit.> up to 200 GeV at RHIC <cit.>. The Cronin effect is also anticipated for quarkonium production in  collisions <cit.>. A similar effect was observed also in semi-inclusive deep-inelastic scattering off nuclei by the HERMES experiment <cit.>. One possible source of this effect is the multiple scattering of the struck parton while traversing the nucleus, which broadens the parton momentum . Under this assumption, the modification of can be related to the transport properties of matter, expressed by the transport coefficient q̂.  
Other effects, like nuclear absorption and parton energy loss, are also expected to contribute when studying particle production in nuclear matter. 

Additional measurements of the  spectrum in 
ep and 
eA at the EIC can help to discriminate between models and constrain their parameters, including the relative role of multiple scattering and nuclear absorption. Such a programme 
will greatly extend the studies pioneered by the HERMES collaboration.

We present here an example of the calculation of the expected modification of the quarkonium energy spectrum in 
eA collisions due to multiple scattering of the parton in the medium. The study is based on an earlier work <cit.>, where a microscopic approach was adopted for the calculation of the decay of J/ψ and Υ in the QGP. Here the QGP medium is replaced with cold nuclear matter, specifically with a large gold 
nucleus, and its properties are constrained taking into account various nuclear effects: nuclear shadowing <cit.>, coherent QCD multiple scattering <cit.>, initial- and final-state parton energy loss <cit.>, and initial and final scattering effects (including multiple scattering) <cit.>.  

To study the nuclear modification, the ratio of cross sections for quarkonium production in reactions that involve a nucleus 
and a proton baseline is used:

    R_AA = 1/⟨ N_ bin⟩dσ_AA/ dσ_pp  ,        
    R_eA =  1/Adσ_eA/ dσ_ep  .

Here, A and the average number of nucleon-nucleon collisions ⟨ N_ bin⟩ provide the relevant normalisation factors such that in the absence of  nuclear modification the ratios are unity. The R_AA presents suppression from QGP, including thermal dissociation in the QGP, while R_eA offers the cold nuclear-matter counterpart.  

A preliminary study demonstrates that most 
quarkonium states show a larger R_eA compared to R_AA, and thus a decreased suppression, with the exception of the state, which sees an increase in suppression by roughly 20%, and of the χ_b(1P) state, which sees a relatively low increase in suppression of roughly 10%. 
The χ_c state experiences a significant decrease of about 50% 
in the suppression factor. The 
Υ states follow 
an analogous trend, with  decreased suppression of around 25% for Υ(1S) and Υ(2S) and 90% for Υ(3S). Finally, χ_b(2P) and  χ_b(3P) show decreases in their suppression factors of roughly 55% and >95%, respectively. 

The overall trend seems to indicate that highly suppressed states see the largest decrease in suppression, while the least suppressed states show either a small decrease or a slight increase in their suppression factors. All states retain a similar amount of E dependence, which is not surprising given that it is assumed that the time for the onset of the interaction is τ_ form.=1 fm. We direct an interested reader to <ref> for more details. 

These preliminary results show that one can expect a significant modification due to cold nuclear-matter effects, which should allow for experimental investigation of these effects at the EIC. Thus, 
quarkonium studies in 
eA collisions at the EIC will help to understand the impact of different transport coefficients on quarkonium production in reactions that involve heavy nuclei and, in turn, help to calibrate quarkonium as a probe of the properties of matter created in high-energy 
pA and AA 
collisions.




























































































§ SUMMARY



Quarkonium is an extremely useful tool to probe the internal structure of matter, namely one of the main goals of the Electron Ion Collider. In this review, we argue that studies of quarkonium production and correlations in (polarised) electron-proton and electron-nucleus collisions can produce unprecedented insights into the 3D structure of the nucleon and into the partonic content of the nuclei as well as help to settle the long-lasting debate on how quarkonia form.

Section <ref> briefly introduced the EIC project, its key parameters, and requirements for an EIC detector. We also defined conventions and basic kinematical quantities useful for describing lepton-hadron reactions. Finally, we made a case for a muon detector for quarkonium studies at the EIC.


Studies of collinear PDFs, form factors, TMD PDFs, GPDs, GTMDs and even double-parton distribution functions can be done at EIC using quarkonium production on a nucleon. In Sections <ref>, <ref> and <ref>, we reviewed the physics case for quarkonium measurements at the EIC. Quarkonium production at large transverse momenta in proton-proton and electron-proton collisions has been studied extensively within the frameworks of NRQCD and collinear factorisation. As discussed in Sections 3.1 and 3.2, it remains a challenge to obtain a simultaneous description of all HERA, LHC and Tevatron data for J/ψ photo- and hadroproduction, η_c hadroproduction, J/ψ+Z hadroproduction, J/ψ polarisation as well as inclusive production in e^+e^- annihilation at B factories. 

Further data from the EIC can help but its p_T reach is limited to 10-15 GeV for charmonia and much less for bottomonia. The focus would then be on low-p_T data. The latter needs to be described within the framework of transverse momentum dependent parton distributions (TMDs) and requires the inclusion of so-called shape functions, which are the subjects of Section 3.3. In this way the EIC will provide new data to further unravel the quarkonium production mechanism, while at the same time offer new ways to employ quarkonium production as a tool to study TMDs and other parton distributions (the subjects of Section 4). This applies especially to gluon TMDs about which currently very little is known. 
Analogous studies can be performed in electron-nucleus collisions (including, among others, insights into transport properties of nuclear matter), which is the subject of section 5. J/ψ polarisation studies can be done, as well as various spin asymmetry measurements, where the electron, proton and light nuclei can be polarised.
All these observables can contribute to our understanding of hadron structure and hadron formation, in particular those involving heavy quarks.

Overall, the physics case for quarkonium physics at the EIC is very extensive and promising.









§ ACKNOWLEDGEMENTS


We thank M. Chithirasreemadam,  M.A. Ozcelik and H.F. Zhang for useful comments and inputs.
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the grant agreement No.824093 (STRONG-2020). 

This project has also received funding from the French ANR under the grant ANR-20-CE31-0015 (“PrecisOnium”).
This work was also partly supported by the French CNRS via the IN2P3 project GLUE@NLO, via the Franco-Chinese LIA FCPPL (Quarkonium4AFTER), via the IEA No.205210 (“GlueGraph") and “Excitonium”, by the Paris-Saclay U. via the P2I Department and by the GLUODYNAMICS project funded by the "P2IO LabEx (ANR-10-LABX-0038)" in the framework "Investissements d’Avenir" (ANR-11-IDEX-0003-01) managed by the Agence Nationale de la Recherche (ANR), France.




C.V.H. has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska–Curie grant agreement No 792684 and from the programme Atracción de Talento, Comunidad de Madrid (Spain), 
under the grant agreement No 2020-T1/TIC-20295.
M.N. has been supported by the Marie Skłodowska-Curie action “RadCor4HEF” under grant agreement No. 101065263.

The work of U.D. and C.P. is supported by Fondazione di Sardegna under the projects “Quarkonium at LHC energies”, No. F71I17000160002 (University of Cagliari) and  ”Proton tomography at the LHC”, No. F72F20000220007 (University of Cagliari).  The work of C.F. and C.P. is supported by the European Union ”Next Generation EU" program through the Italian PRIN 2022 grant n. 20225ZHA7W.

D.K. was supported by the National Science Centre, Poland, under the research grant no. 2018/30/E/ST2/00089.

P.T. is supported by a postdoctoral fellowship fundamental research of the Research Foundation Flanders (FWO) no. 1233422N.



The work of X.Y. was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics grant DE-SC0011090 and currently by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, InQubator for Quantum Simulation (IQuS) (https://iqus.uw.edu) under Award Number DOE (NP) Award DE-SC0020970 via the program on Quantum Horizons: QIS Research and Innovation for Nuclear Science.
The work of V.C. and R.V. was supported by the Office of Nuclear Physics in the U.S. Department of Energy at Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and the LLNL-LDRD Program under Project No. 21-LW-034 and No. 23-LW-036.
The work of I.V. was supported by the Laboratory Directed Research and Development program at Los Alamos National Laboratory. The work of C.S. was supported by the Indonesia Endowment Fund for Education (LPDP). 
S. B. and Y. H. are supported by the U.S. Department of Energy under Contract No. DE-SC0012704, and also by  Laboratory Directed Research and Development (LDRD) funds from Brookhaven Science Associates.
The work of A. M. was supported by the National Science Foundation under grant number PHY-2110472, and also by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the TMD Topical Collaboration. 
J.R.W. was supported by the EIC Center at Jefferson Lab, the LDRD programs of LBNL and by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract number DE-AC02-05CH11231.

L.C. and S.Y. were supported by the Guangdong Major Project of Basic and Applied Basic Research under the project No. 2020A1515010794.
The work of C.E.H. was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics grant DE-FG02-96ER40960.
J.W.Q. is supported in part by the U.S. Department of Energy (DOE) Contract No. DE-AC05-06OR23177, under which Jefferson Science Associates, LLC operates Jefferson Lab.
The work of M.B. was supported by the German Research Foundation DFG through Grant No. BU 3455/1-1.

















§ ESTIMATION OF MEASUREMENT EFFICIENCY




The measurement efficiency is calculated using the decay kinematic simulated with PYTHIA8 <cit.> and two cases for the minimum transverse momentum of the electron measurable in the experiment: ^ele > 0.2 GeV for a detector with a magnetic field B=1.5 T, and > 0.4 GeV for B = 3 T <cit.>. The single electron tracking efficiency is assumed to be 80%. jpsi:efficiency shows the efficiency as a function of  rapidity and : it approximately constant, and for the B = 3 T case, there is a mild decrease of efficiency with increasing due to decay kinematic. For high- , one of the electrons tends to carry the majority of the momentum; thus the of the other falls below the reconstruction threshold. Based on these results, we assume measurement efficiency to be 64%. 





§ NUMERICAL RESULTS FOR NUCLEAR MODIFICATION  R_AA AND R_EA FOR QUARKONIUM PRODUCTION WITHIN THE MICROSCOPIC MODEL PRESENTED IN SEC. <REF>

















§ THE LEPTON, PHOTON AND PARTON DISTRIBUTION IN AN UNPOLARISED ELECTRON










