Quarkonium TMD fragmentation functions in NRQCD

We study the transverse-momentum spectrum of quarkonium production from single light-parton fragmentation mechanism. In the case of semi-inclusive deep inelastic scattering, we observe that there are two possible initiating processes, namely photon-gluon fusion and light-quark fragmentation. For the second case we derive the factorization theorem, which involves a new hadronic quantity: the quarkonium transverse-momentum-dependent fragmentation functions in NRQCD. We calculate their matching onto the non-perturbative long distance matrix elements at the lowest order in the strong-coupling constant (${\mathcal O}(\alpha_s^2)$). Focusing on the case of the electron-ion collider, we make a comparative phenomenological study of the two production mechanisms and find the regions of the phase space where one is dominant over the other.

Despite the abundant interest in quarkonium TMDs, little has been done in the direction of TMD quarkonium fragmentation processes, by which we mean single parton fragmentation mechanism (see for example figure 1(c)). The fragmentation process becomes relevant when a hard scale, Λ, much larger than the quarkonium mass, M , exists. For example, in hadronic colliders the quarkonium TMD fragmentation can be accessed inside jets where quarkonia are found with relatively large transverse momenta. In this scenario, the hard scale of the problem is set by the transverse momentum of the jet, Λ = p jet T , and the TMD spectrum of quarkonium is then measured w.r.t. the jet axis. This type of studies can provide important insight into the heavy-quark hadronization mechanisms and, in terms of contamination from the underlying event (UE), the study of quarkonia is considered advantageous compared to light hadrons. 1 ⇤ < l a t e x i t s h a 1 _ b a s e 6 4 = " L Q G 5 N c S V T b q + k N N 1 O 6 D H i 6 / C q j Q = " > A A A B 7 3 i c d V D L S g M x F M 3 4 r P V V d e k m W B R x M W T a a q e 7 g h u X F e w D 2 r F k 0 k w b m m T G J C O U 0 p 9 w 4 0 I R t / 6 O O / / G 9 C G o 6 I E L h 3 P u 5 d 5 7 w o Q z b R D 6 c J a W V 1 b X 1 j M b 2 c 2 t 7 Z 3 d 3 N 5 + Q 8 e p I r R O Y h 6 r V o g 1 5 U z S u m G G 0 1 a i K B Y h p 8 1 w e D n 1 m / d U a R b L G z N K a C B w X 7 K I E W y s 1 O r 0 s R D 4 9 q y b y y O 3 U v F K f g E i 9 / y i 7 C P P E l Q s l v 0 S 9 F w 0 Q x 4 s U O v m 3 j u 9 m K S C S k M 4 1 r r t o c Q E Y 6 w M I 5 x O s p 1 U 0 w S T I e 7 T t q U S C 6 q D 8 e z e C T y 2 S g 9 G s b I l D Z y p 3 y f G W G g 9 E q H t F N g M 9 G 9 v K v 7 l t V M T + c G Y y S Q 1 V J L 5 o i j l 0 M R w + j z s M U W J 4

S N L M F H M 3 g r J A C t M j I 0 o a 0 P 4 + h T + T x o F 1 y u 6 h e t S v n q y i C M D D s E R O A U e K I M q u A I 1 U A c E c P A A n s C z c + c 8 O i / O 6 7 x 1 y V n M H I A f c N 4 + A R h i j + 4 = < / l a t e x i t >
g < l a t e x i t s h a 1 _ b a s e 6 4 = " d X 5 v e E E Z e V V B D o f X U 5 3 C L a J t h V c = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b F U 0 m q o M e C F 4 8 t 2 F p o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o r e N U M W y x W M S q E 1 C N g k t s G W 4 E d h K F N A o E P g T j 2 5 n / 8 I R K 8 1 j e m 0 m C f k S H k o e c U W O l 5 r B f r r h V d w 6 y S r y c V C B H o 1 / + 6 g 1 i l k Y o D R N U 6 6 7 n J s b P q D K c C Z y W e q n G h L I x H W L X U k k j 1 H 4 2 P 3 R K z q w y I G G s b E l D 5 u r v i Y x G W k + i w H Z G 1 I z 0 s j c T / / O 6 q Q l v / I z L J D U o 2 W J R m A p i Y j L 7 m g y 4 Q m b E x B L K F L e 3 E j a i i j J j s y n Z E L z l l 1 d J u 1 b 1 L q u 1 5 l W l f p 7 H U Y Q T O I U L 8 O A a 6 n A H D W g B A 4 R n e I U 3 5 9 F 5 c d 6 d j 0 V r w c l n j u E P n M 8 f x V e M 0 w = = < / l a t e x i t > q/q < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 B a n I D Q L i j J x T e c + B o G e I D t u z T A = " > A A A B 8 H i c b V B N T w I x E J 3 F L 8 Q v 1 K O X R q L x h L t g o k c S L x 4 x k Q 8 D G 9 I t X W h o u 0 v b N S E b f o U X D x r j 1 Z / j z X 9 j g T 0 o + J J J X t 6 b y c y 8 I O Z M G 9 f 9 d n J r 6 x u b W / n t w s 7 u 3 v 5 B 8 f C o q a N E E d o g E Y 9 U O 8 C a c i Z p w z D D a T t W F I u A 0 1 Y w u p 3 5 r S e q N I v k g 5 n E 1 B d 4 I F n I C D Z W e h x f d g O s 0 v G 0 V y y 5 Z X c O t E q 8 j J Q g Q 7 1 X / O r 2 I 5 I I K g 3 h W O u O 5 8 b G T 7 E y j H A 6 L X Q T T W N M R n h A O 5 Z K L K j 2 0 / n B U 3 R m l T 4 K I 2 V L G j R X f 0 + k W G g 9 E Y H t F N g M 9 b I 3 E / / z O o k J b / y U y T g x V J L F o j D h y E R o 9 j 3 q M 0 W J 4 R N L M F H M 3 o r I E C t M j M 2 o Y E P w l l 9 e J c 1 K 2 a u W K / d X p d p 5 F k c e T u A U L s C D a 6 j B H d S h A Q Q E P M M r v D n K e X H e n Y 9 F a 8 7 J Z o 7 h D 5 z P H 9 N 5 k F Y = < / l a t e x i t > ⇤ < l a t e x i t s h a 1 _ b a s e 6 4 = " L Q G 5 N c S V T b q + k N N 1 O 6 D H i 6 / C q j Q = " > A A A B 7 3 i c d V D L S g M x F M 3 4 r P V V d e k m W B R x M W T a a q e 7 g h u X F e w D 2 r F k 0 k w b m m T G J C O U 0 p 9 w 4 0 I R t / 6 O O / / G 9 C G o 6 I E L h 3 P u 5 d 5 7 w o Q z b R D 6 c J a W V 1 b X 1 j M b 2 c 2 t 7 Z 3 d 3 N 5 + Q 8 e p I r R O Y h 6 r V o g 1 5 U z S u m G G 0 1 a i K B Y h p 8 1 w e D n 1 m / d U a R b L G z N K a C B w X 7 K I E W y s 1 O r 0 s R D 4 9 q y b y y O 3 U v F K f g E i 9 / y i 7 C P P E l Q s l v 0 S 9 F w 0 Q x 4 s U O v m 3 j u 9 m K S C S k M 4 1 r r t o c Q E Y 6 w M I 5 x O s p 1 U 0 w S T I e 7 T t q U S C 6 q D 8 e z e C T y 2 S g 9 G s b I l D Z y p 3 y f G W G g 9 E q H t F N g M 9 G 9 v K v 7 l t V M T + c G Y y S Q 1 V J L 5 o i j l 0 M R w + j z s M U W J 4

S N L M F H M 3 g r J A C t M j I 0 o a 0 P 4 + h T + T x o F 1 y u 6 h e t S v n q y i C M D D s E R O A U e K I M q u A I 1 U A c E c P A A n s C z c + c 8 O i / O 6 7 x 1 y V n M H I A f c N 4 + A R h i j + 4 = < / l a t e x i t >
< l a t e x i t s h a 1 _ b a s e 6 4 = " 8 6 E s r N l B x y / z 2 G m D 9 M S v S / n T H H 8 = " < l a t e x i t s h a 1 _ b a s e 6 4 = " E x y Y 2 c R t M B l r 1 Q + v O y G K X 9 A h B d 4 = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b F U 0 m q o M e C F 4 8 V 7 Q e 0 o W y 2 k 3 T p Z h N 2 N 0 I p / Q l e P C j i 1 V / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k A q u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U 0 k m m G D Z Z I h L V C a h G w S U 2 D T c C O 6 l C G g c C 2 8 H o d u a 3 n 1 B p n s h H M 0 7 R j 2 k k e c g Z N V Z 6 i P p e v 1 x x q + 4 c Z J V 4 O a l A j k a / / N U b J C y L U R o m q N Z d z 0 2 N P 6 H K c C Z w W u p l G l P K R j T C r q W S x q j 9 y f z U K T m z y o C E i b I l D Z m r v y c m N N Z 6 H A e 2 M 6 Z m q J e 9 m f i f 1 8 < l a t e x i t s h a 1 _ b a s e 6 4 = " r 2 i B G M 4 d 9 9 d e j l 5 M U o p I c J E 3 w 1 E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b F U 0 m q o M e C F 4 8 V 7 Q e 0 o W y 2 k 3 T p Z h N 2 N 0 I p / Q l e P C j i 1 V / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k A q u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U 0 k m m G D Z Z I h L V C a h G w S U 2 D T c C O 6 l C G g c C 2 8 H o d u a 3 n 1 B p n s h H M 0 7 R j 2 k k e c g Z N V Z 6 i P q 1 f r n i V t 0 5 y C r x c l K B H I 1 + + a s 3 S F g W o z R M U K 2 7 n p s a f 0 K x a C 0 4 + c w x / I H z + Q P s I Y 1 4 < / l a t e x i t > ⇤ < l a t e x i t s h a 1 _ b a s e 6 4 = " L Q G 5 N c S V T b q + k N N 1 O 6 D H i 6 / C q j Q = " > A A A B 7 3 i c d V D L S g M x F M 3 4 r P V V d e k m W B R x M W T a a q e 7 g h u X F e w D 2 r F k 0 k w b m m T G J C O U 0 p 9 w 4 0 I R t / 6 O O / / G 9 C G o 6 I E L h 3 P u 5 d 5 7 w o Q z b R D 6 c J a W V 1 b X 1 j M b 2 c 2 t 7 Z 3 d 3 N 5 + Q 8 e p I r R O Y h 6 r V o g 1 5 U z S u m G G 0 1 a i K B Y h p 8 1 w e D n 1 m / d U a R b L G z N K a C B w X 7 K I E W y s 1 O r 0 s R D 4 9 q y b y y O 3 U v F K f g E i 9 / y i 7 C P P E l Q s l v 0 S 9 F w 0 Q x 4 s U O v m 3 j u 9 m K S C S k M 4 1 r r t o c Q E Y 6 w M I 5 x O s p 1 U 0 w S T I e 7 T t q U S C 6 q D 8 e z e C T y 2 S g 9 G s b I l D Z y p 3 y f G W G g 9 E q H t F N g M 9 G 9 v K v 7 l t V M T + c G Y y S Q 1 V J L 5 o i j l 0 M R w + j z s M U W J 4 S N L M F H M 3 g r J A C t M j I 0 o a 0 P 4 + h T + T < l a t e x i t s h a 1 _ b a s e 6 4 = " Z E D C m B q h C K o y w I r 1 R X h B p l 5 2 I k M = " > A A A C G H i c d V D L S s N A F J 3 U V 6 2 v q k s 3 g 0 W p I G n S i h b c F N 2 4 b N H a Q p K G y X T a D p 0 8 m J k I J e Q z 3 P g r b l w o 4 r Y 7 / 8 b p Q 9 C i B y 4 c z r m X e + / x I k a F N I On the other hand, in semi-inclusive deep inelastic scattering (SIDIS) the hard scale Λ is usually set by the invariant mass of the virtual photon, Λ = Q = −q 2 , with q the photon momentum. Thus, performing the TMD measurement with respect to the photon's direction in the Breit frame we can formulate a global TMD factorization theorem in a similar manner to light hadrons. Since in an electron-hadron collider the kinematical range of x and Q 2 is limited, in the present study we propose to explore the window in which single-quarkonium TMD fragmentation is relevant compared to the photon-gluon fusion quarkonium production.
In our analysis we employ the NRQCD factorization conjecture where quarkonium is produced at large distances through the hadronization of a heavy quark-antiquark pair, QQ(n). The pair can be found in any color and angular configuration n = 2S+1 L [col.] J but then the probability that the pair decays in the colorless quarkonium state scales with the relative velocity, v, of the quark-antiquark pair in the quarkonium rest frame. In this study we are primarily interested in the quarkonium state J/ψ, for which the four leading channels in v are: 3 S We are going to make our case for the forthcoming electron-ion collider (EIC). In order to start our discussion let us explore quarkonium production for this type of experiment. At leading order (LO) in the perturbative expansion, quarkonium production is possible through the photon-gluon fusion process as shown in the left panel of figure 1(a). In this case the heavy-quark pair is produced either in 1 S [8] 0 or 3 P [8] J=0,2 states, the transverse momentum of the heavy-quark pair vanishes and z ≡ p h · P/p h · q = 1 (where p h is the target momentum and P the heavy-quarkonium momentum). At higher orders a smearing around the Born kinematics will occur due the emissions of soft and collinear gluons from the incoming partons as well as the heavy quark-antiquark state. The cross section for these processes can be organized in terms of the TMDPDFs and the recently introduced TMD quarkonium shape-functions [14,15]. At O(a 2 s ) we also have a contribution from the 3 S (see figure 1(b)). Although this contribution is power suppressed in q 2 T /Q 2 , enhancements from the relative velocity scaling make it hard to argue that this channel can be neglected. At the same order in α s enters the quarkonium fragmentation process shown in figure 1(c), whose relative channel is To isolate the fragmentation process we take advantage of the fact that both 3 S 1 at Born level have continuous distributions in z and, thus, focussing on the region z 0.5 allows us to better isolate these two channels. In addition, to further suppress the non-fragmentation contributions we work in the region Q M ∼ P ⊥ (for details see discussion in section 4). We explore here how the competition of gluon-photon fusion and quark fragmentation depends on the kinematical conditions in which an experiment is run, in order to establish the study case of this process.
The factorization of the TMD cross section for the light-quark fragmentation mechanism is structurally very similar to the conventional SIDIS process [26][27][28][29], where the rapidity scale ζ and the renormalization scale µ are responsible of the evolution of the TMD matrix elements in initial and final states. This process has been extensively studied in the literature and we now have high-order QCD calculations for hard factors [30][31][32], evolution kernel [33][34][35] and TMDPDFs [36][37][38][39][40][41][42]. It is also possible then to use extraction of TMDs [43][44][45] to perform a phenomenological analysis. In this work we use the extractions done in ref. [44] with the so-called ζ-prescription [46] and the associated Artemide code [47].
The quarkonium TMD fragmentation function (TMDFF), D f →H , is the only piece which remains to be studied and included in our analysis. We decompose the TMDFF within NRQCD in terms of calculable short-distance matching coefficients and long-distance matrix elements (LDMEs). We proceed then to the LO calculation of this function, extracting the matching coefficient onto the corresponding LDMEs in the region q T ∼ M . Using this information we can evaluate the contribution to the cross section from light-quark fragmentation.
This paper is organized as follows: we introduce our notation and the form of the factorization theorem for light quark fragmentation in section 2. In section 3 we perform the calculation of the matching coefficient of the TMDFF onto the 3 S [8] 1 LDME at O(α 2 s ). We also include a short discussion on the RG evolution of LDMEs which becomes relevant for TMD calculations due to the additional soft scale q T . In section 4 we give an analytic comparison of the fixed-order results for channels for DIS in the small transverse momentum limit. In section 5 we demonstrate the phenomenological applicability of this formalism for the future EIC kinematics and provide for the first time numerical estimates of the relative importance of the fragmentation channel in various kinematic regions. Finally, in section 6 we conclude.

Notation
The momenta of quarkonium production in SIDIS are specified by where is a lepton, h and H are respectively the initial and the final hadrons, and X is the undetected final state. The masses of the hadrons are and we neglect lepton masses. For the moment we continue to consider hadron masses in order to control their effects. The differential cross-section in eq. (2.1) can be written as with q = l − l being the momentum of the intermediate photon and α em = e 2 /4π the QED coupling constant. The cross section in eq. (2.3) is proportional to the phase-space differentials for the detected lepton and heavy hadron, with E and E H being their energies. The leptonic and hadronic tensors (L µν and W µν ) are where e is the lepton charge, and J µ is the electro-magnetic current.

Kinematical variables for SIDIS
In SIDIS one makes use of different frames which specify the hadronic variables. In this section we follow closely [44], adapting the variables to the heavy quarkonium case. Although in our numerical evaluation we only consider some approximation of the kinematical variables, we prefer to write them down for future reference. The factorization theorem for the cross section is done in the Breit frame, where the momenta of hadrons are back to back and they are respectively with n 2 =n 2 = 0, (nn) = 1, and the vector decomposition The transverse component of a vector is defined by the projection We also use the convention that the bold font denotes vectors that have only transverse components and The kinematical scalar variables are defined as: Experimentally one needs to define a plane transverse to q and p, and the projector corresponding to the components transverse to this plane is given by the tensor g µν ⊥ defined as Using this notation we can distinguish the transverse components of v µ in the hadronic Breit frame as v µ T , see eq. (2.7), from the transverse components projected by g ⊥ , that is v µ ⊥ . The mass corrections are conveniently described using the following combinations The definition of ς 2 ⊥ in eq. (2.11) contains P 2 ⊥ = −P µ P ν g µν ⊥ = −P 2 ⊥ . Using these variables we can re-express g µν T in eq. (2.7) in terms of hadronic momenta With eq. (2.10) and eq. (2.12) one can derive the relation between q 2 T = q µ q ν g µν T and P 2 ⊥ , including mass corrections. Using these definition we can rewrite the elements of the SIDIS cross-section formula in terms of experimental variables, that is, the differential volumes of the phase space are where ψ is the azimuthal angle of scattered lepton, and ϕ is the azimuthal angle of the produced hadron.
Finally we introduce the variables x S and z S , that are the collinear fractions of parton momentum that include kinematic power corrections, which are invariant under boosts along the directions n andn, but are not invariant under a generic Lorentz transformation. The variables x S and z S in eq. (2.15) are (2.17) The kinematic corrections presented above are usually small when Q M, m. In this case the relation between observed and factorization variables simplifies: The next correction to this limit happens when Q M ∼ q T m, that is in the limit γ → 0 with ς ∼ O(1). In this limit Notice that in this case there is a shift effect in the variable z S which is independent of the transverse momentum. Nevertheless, the effect of hadron masses can be considered as a power correction to the process whose impact can be estimated numerically.

Factorization of the hadronic tensor
For Q M the proof of the factorization theorem for quarkonium fragmentation follows the same steps as for the light-hadron case, so we refer the reader to the literature [26][27][28][48][49][50][51][52] for a detailed treatment. The factorized hadronic tensor reads J 0 is the Bessel function, the index f runs through all quark flavors (including anti-quarks) and e f is the fractional charge of a quark measured in units of e. The function C V is the matching coefficient for the vector current. We also especify the factorization (µ) and rapidity (ζ) scales, and we have omitted the details of operator definitions, such as T (T )-ordering, color and spinor indices, and rapidity and ultraviolet renormalization factors.
The bare unsubtracted unpolarized TMDPDF and quarkonium TMDFF are defined as and similarly for antiquarks. Here, W v (x) are Wilson lines at x and pointing along vector v at infinity. In the case of SIDIS, the Wilson lines in the TMDPDF (TMDFF) are pointing respectively at past (future) infinity. Notice that the operator that appears in the fragmentation function is the same as in the case of light quarks and that the only change is the presence of a quarkonium in the final state.
The renormalized TMDs are defined in the usual way where Z 2 is the wave function renormalization constant for quarks and Z q the renormalization factor for the UV divergences of the TMD operator. The factor R q is the "rapidity renormalization factor". Its definition comes from the TMD factorization theorem and reads where S(b) is the soft function and Zb denotes the zero-bin contribution, i.e. the soft overlap of the collinear and soft sectors present in the factorization theorem [26][27][28]53]. The soft function in the present case is the same as in the light quark case; for SIDIS it reads (here b = |b T |) The Wilson lines are defined as usual The transverse gauge links T n(n) are to be used in singular gauges, like the light-cone gauge n · A = 0 (orn · A = 0) [54][55][56].
The definition of the zero-bin depends on the rapidity regularization used in loop calculations (see e.g. discussion in [28]). With the modified δ-regularization, the zero-bin subtraction is equal to the soft function: Zb = S(b). As a result in the modified δ-regularization the rapidity renormalization factor is Due to the process independence of the soft function [57], the factor R q is also process independent.

Leptonic tensor and cross section
The leptonic tensor for our case is the same as for the unpolarized SIDIS process: The contraction of the leptonic tensor with the hadronic tensor can be performed with the definition of the azimuthal angle of a produced hadron as [58]: As in the light-quark case, the kinematical rearrangements of the variables contribute to the cos φ and cos 2φ terms in the second line of eq. (2.30), which disappear when one integrates over angles [59].
Recollecting the equations for the cross-section in eq. (2.3), the differential phase-space volume in eq. (2.14), the hadronic tensor in eq. (2.21), the leptonic tensor in eq. (2.30), and integrating over the azimuthal angles we get where q 2 T , x S and z S are functions of P 2 ⊥ , x and z, and are defined respectively in eqs. (2.13), (2.16) and (2.17). In eq. (2.32) we have reported all calculable power corrections which come from the kinematics of the process. Other power corrections of the same order can be obtained from the g < l a t e x i t s h a 1 _ b a s e 6 4 = " d X 5 v e E E Z e V V B D o f X U 5 3 C L a J t h V c = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b F U 0 m q o M e C F 4 8 t 2 F p o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o r e N U M W y x W M S q E 1 C N g k t s G W 4 E d h K F N A o E P g T j 2 5 n / 8 I R K 8 1 j e m 0 m C f k S H k o e c U W O l 5 r B f r r h V d w 6 y S r y c V C B H o 1 / + 6 g 1 i l k Y o D R N U 6 6 7 n J s b P q D K c C Z y W e q n G h L I x H W L X U k k j 1 H 4 2 P 3 R K z q w y I G G s b E l D 5 u r v i Y x G W k + i w H Z G 1 I z 0 s j c T / / O 6 q Q l v / I z L J D U o 2 W J R m A p i Y j L 7 m g y 4 Q m b E x B L K F L e 3 E j a i i j J j s y n Z E L z l l 1 d J u 1 b 1 L q u 1 5 l W l f p 7 H U Y Q T O I U L 8 O A a 6 n A H D W g B A 4 R n e I U 3 5 9 F 5 c d 6 d j 0 V r w c l n j u E P n M 8 f x V e M 0 w = = < / l a t e x i t > q/q < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 B a n I D Q L i j J x T e c + B o G e I D t u z T A = " > A A A B 8 H i c b V B N T w I x E J 3 F L 8 Q v 1 K O X R q L x h L t g o k c S L x 4 x k Q 8 D G 9 I t X W h o u 0 v b N S E b f o U X D x r j 1 Z / j z X 9 j g T 0 o + J J J X t 6 b y c y 8 I O Z M G 9 f 9 d n J r 6 x u b W / n t w s 7 u 3 v 5 B 8 f C o q a N E E d o g E Y 9 U O 8 C a c i Z p w z D D a T t W F I u A 0 1 Y w u p 3 5 r S e q N I v k g 5 n E 1 B d 4 I F n I C D Z W e h x f d g O s 0 v G 0 V y y 5 Z X c O t E q 8 j J Q g Q 7 1 X / O r 2 I 5 I I K g 3 h W O u O 5 8 b G T 7 E y j H A 6 L X Q T T W N M R n h A O 5 Z K L K j 2 0 / n B U 3 R m l T 4 K I 2 V L G j R X f 0 + k W G g 9 E Y H t F N g M 9 b I 3 E / / z O o k J b / y U y T g x V J L F o j D h y E R o 9 j 3 q M 0 W J 4 R N L M F H M 3 o r I E C t M j M 2 o Y E P w l l 9 e J c 1 K 2 a u W K / d X p d p 5 F k c e T u A U L s C D a 6 j B H d S h A Q Q E P M M r v D n K e X H e n Y 9 F a 8 7 J Z o 7 h D 5 z P H 9 N 5 k F Y = < / l a t e x i t >  < l a t e x i t s h a 1 _ b a s e 6 4 = " E x y Y 2 c R t M B l r 1 Q + v O y G K X 9 A h B d 4 = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b F U 0 m q o M e C F 4 8 V 7 Q e 0 o W y 2 k 3 T p Z h N 2 N 0 I p / Q l e P C j i 1 V / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k A q u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U 0 k m m G D Z Z I h L V C a h G w S U 2 D T c C O 6 l C G g c C 2 8 H o d u a 3 n 1 B p n s h H M 0 7 R j 2 k k e c g Z N V Z 6 i P p e v 1 x x q + 4 c Z J V 4 O a l A j k a / / N U b J C y L U R o m q N Z d z 0 2 N P 6 H K c C Z w W u p l G l P K R j T C r q W S x q j 9 y f z U K T m z y o C E i b I l D Z m r v y c m N N Z 6 H A e 2 M 6 Z m q J e 9 m f i f 1 8 1 M e O N P u E w z g 5 I t F o W Z I C Y h s 7 / J g C t k R o w t o U x x e y t h Q 6 o o M z a d k g 3 B W 3 5 5 l b R q V e + y W r u / q t T P 8 z i K c A K n c A E e X E M d 7 q A B T W A Q w T O 8 w p s j n B f n 3 f l Y t B a c f O Y Y / s D 5 / A H q n Y 1 3 < / l a t e x i t > g 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " r 2 i B G M 4 d 9 9 d e j l 5 M U o p I c J E 3 w 1 E = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b F U 0 m q o M e C F 4 8 V 7 Q e 0 o W y 2 k 3 T p Z h N 2 N 0 I p / Q l e P C j i 1 V / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k A q u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U 0 k m m G D Z Z I h L V C a h G w S U 2 D T c C O 6 l C G g c C 2 8 H o d u a 3 n 1 B p n s h H M 0 7 R j 2 k k e c g Z N V Z 6 i P q 1 f r n i V t 0 5 y C r x c l K B H I 1 + + a s 3 S F g W o z R M U K 2 7 n p s a f 0 K V 4 U z g t N T L N K a U j W i E X U s l j V H 7 k / m p U 3 J m l Q E J E 2 V L G j J X f 0 9 M a K z 1 O A 5 s Z 0 z N U C 9 7 M / E / r 5 u Z 8 M a f c J l m B i V b L A o z Q U x C Z n + T A V f I j B h b Q p n i 9 l b C h l R R Z m w 6 J R u C t / z y K m n V q t 5 l t X Z / V a m f 5 3 E U 4 Q R O 4 Q I 8 u I Y 6 3 E E D m s A g g m d 4 h T d H O C / O u / O x a C 0 4 + c w x / I H z + Q P s I Y 1 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Z E D C m B q h C K o y w I r 1 R X h B p l 5 2 I k M = " > A A A C G H i c d V D L S s N A F J 3 U V 6 2 v q k s 3 g 0 W p I G n S i h b c F N 2 4 b N H a Q p K G y X T a D p 0 8 m J k I J e Q z 3 P g r b l w o 4 r Y 7 / 8 b p Q 9 C i B y 4 c z r m X e + / x I k a F N I x P L b O 0 v L K 6 l l 3 P b W x u b e / k d / f u R R h z T J o 4 Z C F v e 0 g Q R g P S l F Q y 0 o 4 4 Q b 7 H S M s b X k / 8 1 g P h g o b B n R x F x P F R P 6 A 9 i p F U k p s v N a D t I Z 4 0 U l j s m L e u 0 U m s q p P a V 7 S f l F L 7 s l O p u 4 l x W k 5 n + o m b L x i 6 M Q U 0 9 P M F Y s 6 t A p i j 7 u b H d j f E s U 8 C i R k S w j K N S D o J 4 p J i R t K c H Q s S I T x E f W I p G i C f C C e Z P p b C I 6 V 0 Y S / k q g I J p + r P i Q T 5 Q o x 8 T 3 X 6 S A 7 E o j c R / / K s W P a q T k K D K J Y k w L N F v Z h B G c J J S r B L O c G S j R R B m F N 1 K 8 Q D x B G W K s u c C u H 7 U / g / u S / r Z k U v N 8 4 K t e N 5 H F l w A A 5 B E Z j g A t T A D a i D J s D g E T y D V / C m P W k v 2 r v 2 M W v N a P O Z f f A L 2 v g L h 3 6 e I A = = < / l a t e x i t > QQ( 3 S Figure 2. Lowest order diagrams for quarkonium TMDFF from a light quark. expansion of the cross section. Since their study goes beyond the purpose of this work and given that they are still largely unexplored, we omit them in the present study.

Quarkonium TMDFF from a light quark
In the present section we report the calculation of the TMD fragmentation function as defined in eq. (2.23). In our approach and for q T ∼ M , we employ NRQCD factorization conjecture [18] to write the quarkonium TMDFF as a product of short distance coefficients, d q→QQ(n) , and the standard NRQCD long distance matrix elements (LDMEs): The short-distance coefficients can be calculated perturbatively through matching. Since they do not depend on the quarkonium state H, they can be calculated perturbatively by choosing an appropriate partonic state. On the other hand, the LDMEs O H (n) are non-perturbative numbers extracted from experimental studies, which scale with the relative velocity, v, of the heavy quark-antiquark pair in the quarkonium rest frame.
In this study we are primarily interested in the quarkonium state H = J/ψ for which the four leading channels in v are: 3 S 0 (∼ v 7 ) and 3 P where we use the shorthand notation ψ ≡ J/ψ. Contributions to light-quark fragmentation at LO in pQCD are shown in figure 2. Note that in our convention here the short distance matching coefficients are defined after renormalization and the inclusion of the soft functions and zero-bin subtraction. However, at the order we are working the soft function contribution is trivial and therefore we have only contributions from the direct matching of the collinear matrix element onto the NRQCD LDME. This is consistent with the fact that there are no rapidity divergent diagrams at this order. These will occur at higher orders in the perturbative expansion, where the soft function will enter into play to cancel them. In the calculation of the diagrams A, B1, B2, and C the only relevant master integral that we have to calculate is with B = b 2 /4 and K i is the i-th Bessel-K function. In the perturbative calculation we used the approximation m c M/2. The results for each diagram (d i ) reported in figure 2 are (we use the common notationz ≡ 1 − z and ∆ ≡ M 2z /z 2 ) where we have used .

(3.9)
Λ's are the Lorentz boosts from the heavy quarkonium center of mass frame to the boosted collinear frame and the expansion of the spinors is performed using the formulas found in the appendix of ref. [60]. With this we have completed the QCD part of the calculation. The pNRQCD part of the calculation is rather simple. We first rewrite the LDME in terms of the relativistically normalized matrix element We then perform the perturbative calculation by replacing

(3.12)
This then yields 2M 3 S From these we immediately have (see eq. (3.2)) (3.14) We also observe that the most divergent part of diagrams cancel in the sum. It is interesting to check also this result in the limit z → 1. In this case, neglecting the function F −1 , only the diagram A contributes so that we have , (3.15) which agrees with the collinear integrated fragmentation function in the same limit calculated in ref. [61]. For completeness we have performed the same calculation in momentum space. Here we only give the final result in d = 4 dimensions, , (3.16) which can also be computed from the Fourier transform of the impact parameter space result in eq. (3.14), We will use this result in the following section in order to make the comparison of the quarkfragmentation and the photon-gluon fusion processes in the large Q 2 regime, Q 2 M 2 ∼ q 2 T . ⇤ < l a t e x i t s h a 1 _ b a s e 6 4 = " L Q G 5 N c S V T b q + k N N 1 O 6 D H i 6 / C q j Q = " > A A A B 7 3 i c d V D L S g M x F M 3 4 r P V V d e k m W B R x M W T a a q e 7 g h u X F e w D 2 r F k 0 k w b m m T G J C O U 0 p 9 w 4 0 I R t / 6 O O / / G 9 C G o 6 I E L h 3 P u 5 d 5 7 w o Q z b R D 6 c J a W V 1 b X 1 j M b 2 c 2 t 7 Z 3 d 3 N 5 + Q 8 e p I r R O Y h 6 r V o g 1 5 U z S u m G G 0 1 a i K B Y h p 8 1 w e D n 1 m / d U a R b L G z N K a C B w X 7 K I E W y s 1 O r 0 s R D 4 9 q y b y y O 3 U v F K f g E i 9 / y i 7 C P P E l Q s l v 0 S 9 F w 0 Q x 4 s U O v m 3 j u 9 m K S C S k M 4 1 r r t o c Q E Y 6 w M I 5 x O s p 1 U 0 w S T I e 7 T t q U S C 6 q D 8 e z e C T y 2 S g 9 G s b I l D Z y p 3 y f G W G g 9 E q H t F N g M 9 G 9 v K v 7 l t V M T + c G Y y S Q 1 V

l a t e x i t s h a 1 _ b a s e 6 4 = " d X 5 v e E E Z e V V B D o f X U 5 3 C L a J t h V c = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b F U 0 m q o M e C F 4 8 t 2 F p o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o r e N U M W y x W M S q E 1 C N g k t s G W 4 E d h K F N A o E P g T j 2 5 n / 8 I R K 8 1 j e m 0 m C f k S H k o e c U W O l 5 r B f r r h V d w 6 y S r y c V C B H o 1 / + 6 g 1 i l k Y o D R N U 6 6 7 n J s b P q D K c C Z y W e q n G h L I x H W L X U k k j 1 H 4 2 P 3 R K z q w y I G G s b E l D 5 u r v i Y x G W k + i w H Z G 1 I z 0 s j c T / / O 6 q Q l v / I z L J D U o 2 W J R m A p i Y j L 7 m g y 4 Q m b E x B L K F L e 3 E j a i i j J j s y n Z E L z l l 1 d J u 1 b 1 L q u 1 5 l W l f p 7 H U Y Q T O I U L 8 O A a 6 n A H D W g B A 4 R n e I U 3 5 9 F 5 c d 6 d j 0 V r w c l n j u E P n M 8 f x V e M 0 w = = < / l a t e x i t > q/q < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 B a n I D Q L i j J x T e c + B o G e I D t u z T A = " > A A A B 8 H i c b V B N T w I x E J 3 F L 8 Q v 1 K O X R q L x h L t g o k c S L x 4 x k Q 8 D G 9 I t X W h o u 0 v b N S E b f o U X D x r j 1 Z / j z X 9 j g T 0 o + J J J X t 6 b y c y 8 I O Z M G 9 f 9 d n J r 6 x u b W / n t w s 7 u 3 v 5 B 8 f C o q a N E E d o g E Y 9 U O 8 C a c i Z p w z D D a T t W F I u A 0 1 Y w u p 3 5 r S e q N I v k g 5 n E 1 B d 4 I F n I C D Z W e h x f d g O s 0 v G 0 V y y 5 Z X c O t E q 8 j J Q g Q 7 1 X / O r 2 I 5 I I K g 3 h W O u O 5 8 b G T 7 E y j H A 6 L X Q T T W N M R n h A O 5 Z K L K j 2 0 / n B U 3 R m l T 4 K I 2 V L G j R X f 0 + k W G g 9 E Y H t F N g M 9 b I 3 E / / z O o k J b / y U y T g x V J L F o j D h y E R o 9 j 3 q M 0 W J 4 R N L M F H M 3 o r I E C t M j M 2 o Y E P w l l 9 e J c 1 K 2 a u W K / d X p d p 5 F k c e T u A U L s C D a 6 j B H d S h A Q Q E P M M r v D n K e X H e n Y 9 F a 8 7 J Z o 7 h D 5 z P H 9 N 5 k F Y = < / l a t e x i t >
⇤ < l a t e x i t s h a 1 _ b a s e 6 4 = " L Q G 5 N c S V T b q + k N N 1 O 6 D H i 6 / C q j Q = " > A A A B 7 3 i c d V D L S g M x F M 3 4 r P V V d e k m W B R x M W T a a q e 7 g h u X F e w D 2 r F k 0 k w b m m T G J C O U 0 p 9 w 4 0 I R t / 6 O O / / G 9 C G o 6 I E L h 3 P u 5 d 5 7 w o Q z b R D 6 c J a W V 1 b X 1 j M b 2 c 2 t 7 Z 3 d 3 N 5 + Q 8 e p I r R O Y h 6 r V o g 1 5 U z S u m G G 0 1 a i K B Y h p 8 1 w e D n 1 m / d U a R b L G z N K a C B w X 7 K I E W y s 1 O r 0 s R D 4 9 q y b y y O 3 U v F K f g E i 9 / y i 7 C P P E l Q s l v 0 S 9 F w 0 Q x 4 s U O v m 3 j u 9 m K S C S k M 4 1 r r t o c Q E Y 6 w M I 5 x O s p 1 U 0 w S T I e 7 T t q U S C 6 q D 8 e z e C T y 2 S g 9 G s b I l D Z y p 3 y f G W G g 9 E q H t F N g M 9 G 9 v K v 7 l t V M T + c G Y y S Q 1 V

L A o z Q U x C Z n + T A V f I j B h b Q p n i 9 l b C h l R R Z m w 6 J R u C t / z y K m n V q t 5 l t X Z / V a m f 5 3 E U 4 Q R O 4 Q I 8 u I Y 6 3 E E D m s A g g m d 4 h T d H O C / O u / O x a C 0 4 + c w x / I H z + Q P s I Y 1 4 < / l a t e x i t >
⇤ < l a t e x i t s h a 1 _ b a s e 6 4 = " L Q G 5 N c S V T b q + k N N 1 O 6 D H i 6 / C q j Q = " > A A A B 7 3 i c d V D L S g M x F M 3 4 r P V V d e k m W B R x M W T a a q e 7 g h u X F e w D 2 r F k 0 k w b m m T G J C O U 0 p 9 w 4 0 I R t / 6 O O / / G 9 C G o 6 I E L h 3 P u 5 d 5 7 w o Q z b R D 6 c J a W V 1 b X 1 j M b 2 c 2 t 7 Z 3 d 3 N 5 + Q 8 e p I r R O Y h 6 r V o g 1 5 U z S u m G G 0 1 a i K B Y h p 8 1 w e D n 1 m / d U a R b L G z N K a C B w X 7 K I E W y s 1 O r 0 s R D 4 9 q y b y y O 3 U v F K f g E i 9 / y i 7 C P P E l Q s l v 0 S 9 F w 0 Q x 4 s U O v m 3 j u 9 m K S C S k M 4 1 r r t o c Q E Y 6 w M I 5 x O s p 1 U 0 w S T I e 7 T t q U S C 6 q D 8 e z e C T y 2 S g 9 G s b I l D Z y p 3 y f G W G g 9 E q H t F N g M 9 G 9 v K v 7 l t V M T + c G Y y S Q 1 V   1 LDME associated with the mixing through renormalization group evolution with the P-wave matrix elements.

Renormalization of LDMEs and channel mixing
In this section we aim to rewrite eq. (3.1) considering resummation of logarithms of the ratio µ/µ f , where µ f is the scale at which the LDMEs are extracted. The most common choice for phenomenological extractions is µ f = M and µ ∼ 1/b. A complication to this attempt comes from the mixing of the 3 S [8] 1 and 3 P [8] J mechanisms under renormalization group equations (RGEs). At fixed order the mixing occurs when a real ultra-soft gluon is exchanged between heavy quarks (antiquarks) in the 3 S the order α s contributions to the 3 S [8] 1 LDME is scaleless and proportional to the P-wave LO matrix element. Using the 1/ UV − 1/ IR prescription we obtain for the sum of diagrams in figure 3: where B F = (N 2 c − 4)/(4N c ) and the ⊗ vertex in figure 3 corresponds to the bilinear operator: ψ † σT A χ. This UV divergence is removed through the renormalization of the LDMEs, which then leads to the following RGE: From eq. (3.1) this immediately implies that the matching coefficients d q→QQ(n) must also satisfy a non-diagonal RGE. A similar observation was made for the quarkonium TMD shape functions in ref. [15]. The major difference here is that the non-trivial evolution is performed within the TMD fragmentation function and therefore the form of the factorized cross section is not modified from the standard case, in contrast to what happens in [15]. Therefore the RGE satisfied by the short-distance matching coefficients is To confirm this equation we need to calculate the contributions to the quarkonium TMDFF at NLO, which is beyond the scope of this study. However, assuming that NRQCD factorization holds we may resum logarithms of the form ln(µ/µ f ) by solving eq. (3.19). Following the notation of [15] and up to next-to-leading-logarithmic (NLL) accuracy we have  .21) is the resummation of ln(µ/µ f ) and in the latter it is made explicit that TMD evolution is satisfied by the fragmentation function D q→H and not the short distance matching coefficients d q→QQ .

On quarkonium TMDFFs from a gluon
This study focuses primarily on quarkonium fragmentation in SIDIS and, for this process, the main interest is the light-quark fragmentation functions. However, a brief mention to the gluon fragmentation process is in order, since this formalism can also be applied in phenomenologically interesting processes beyond DIS. Particularly, gluon fragmentation will play an important role in quarkonium fragmentation within jets, which has attracted significant interest in recent years. This process can be accessed in hadron colliders and particularly at the LHC, where jets are produced in abundance. Measurements [63,64] of the momentum fraction z T ≡ p H T /p jet T for in-jet quarkonia have already revealed interesting aspects of quarkonium production and a striking deviation from Monte-Carlo implementations. An explanation for these deviations was given in ref. [65] and a better qualitative agreement with the data was provided through the fragmentation picture of quarkonium production at large transverse momenta. A multi-differential aspect of these measurements, such as including the transverse momentum w.r.t. the jet axis, will enable a three-dimensional picture of quarkonium fragmentation and an in-depth understanding of its production mechanisms. In the massless limit, i.e. q T M , the quarkonium TMD fragmentation within jets was studied in refs. [8,66]. The framework discussed here can be extended and applied to such studies in the small transverse momentum limit, q T M .
The LO contribution to gluon TMDFFs comes from the color-octet 3 S [8] 1 channel through a simple decay of an off-shell gluon to the heavy quark-antiquark pair. At NLO, contributions from both 1 S [8] 0 and 3 P [8] J channels are introduced. The color-singlet channel 3 S [1] 1 appears only at NNLO. Although these contributions are only relevant at higher orders in the perturbative expansion, recent phenomenological studies [65] of in-jet quarkonia suggest that numerical enhancements from the LDMEs require all four channels to describe the spectrum of quarkonia within jets.
Furthermore, for the 3 S 1 channel, light-quark fragmentation and gluon fragmentation are inevitably related, since at this order the heavy quark-antiquark pair can only be produced by a subsequent decay of an off-shell gluon (see figure 2). It is therefore implied that the channel mixing, through RG evolution as discussed in section 3.1, will also be relevant for gluon TMDFFs.

Quarkonium from photon-gluon fusion
The effect of photon-gluon fusion in the differential cross section that we study is highly complex in the TMD framework that we are using, because higher-twist gluon operators appear already at lowest order. To organize the discussion we consider two separate regions of z: region-1 where (1 − z) 1 and region-2 where (1 − z) ∼ 1. In region-1 one can study the small transverse momentum region q 2 T M 2 from channels 1 S [8] 0 and 3 P [8] 0,2 within the framework of shape-functions in SCET Q [14,15]. In addition, sub-leading operators in SCET Q can give contributions from the channels 3 S 1 . While these contributions are sub-leading in the q 2 T /Q 2 power-counting, the color-singlet channel can get significant enhancement of order 1/v 4 , where v is the relative velocity of the heavy quark-antiquark pair in the quarkonium rest frame. A proper description of region-1 is still in progress.
In this work we focus on region-2. In this case beyond the quark fragmentation that is discussed in section 3, we also have a contribution from the color-singlet channel from photon-gluon fusion process. In this section we give an estimate of this part of the cross section coming from perturbative QCD, that is, looking at the hard part of this process. In this way we get an order of magnitude understanding of both quark fragmentation and photon-gluon fusion (only color-singlet channel). The cross section that we are looking for is the one of the process depicted in figure 1 (b), already studied a long ago in [67].
Considering the result for the cross section in [67] integrated over all angles and in the limit Q 2 M 2 ∼ q 2 T , and expressing the result with the help of [68] one gets where and e H = 2/3 (e H = 1/3) for charmonium (bottomonium). We would like to compare this result with the quark fragmentation process in the same limit (Q 2 M 2 ∼ q 2 T ). The contribution to the cross section from this process at this limit can be expressed directly from eq. (2.32), where for the relative velocity v we have v 2 ∼ 0.3 for charmonium states and v 2 ∼ 0.1 for bottomonium. In table 1 we give the values of this estimate for different values of Q. We find that for charmonium sates and for Q > 20 GeV the photon-gluon fusion (only color-singlet channel) contributes less than 10% compared to the quark fragmentation process. The numerical relative contributions of the full fixed-order photon-gluon fusion in color-singlet channel and the light-quark fragmentation processes, which are sensitive to the values of the LDMEs and collinear PDFs, are discussed in the following section.

Numerics
The relative importance of photon-gluon fusion and quark fragmentation for different intervals of the kinematic variables is highly non-trivial. In order to plot the contribution to the cross section for the quark fragmentation we need some knowledge of the TMDPDFs as provided by other processes. We use here the ζ-prescription [46] and the phenomenological results of ref. [44]. The ζ-prescription allows to express the cross section in terms of scale-independent quantities: in this way we implement the TMD evolution kernel at N 3 LL, the TMDPDFs at NNLO, and the quarkonium TMDFF at LO, which we have calculated in this work. The non-perturbative parameters for the TMDPDF are chosen consistently with the PDF set NNPDF31 nnlo as 0118 [69]. The numerical values for the LDMEs we use are taken from ref. [70] and are given here in table 2.
We have considered two possible settings for the EIC, at √ s = 63 GeV and √ s = 140 GeV, which are expected to be typical for this collider [71]. In both settings we calculate the contribution to the cross section given respectively by photon-gluon fusion to color singlet and quark fragmentation processes. The results are shown in figures 4 and 5. The relative contribution of each input depends crucially on the values of Q, x and z. The general trend is that quark-fragmentation dominates for high Q, high x and low z. The enhancement at large Q is expected from the scaling relation in eq. (4.6), while the enhancement at large x is purely due to the values of the PDFs (gluon PDF in photon-gluon fusion is larger that quark PDFs in quark fragmentation at low x, and viceversa). We have plotted the cross section integrated over several intervals of these variables to show this effect more quantitatively.

Conclusion
Quarkonium production at the Electron-Ion Collider is generally considered as one of the most important processes for the study of gluon TMDs. We have shown that the physics for these heavy-quark events is highly non-trivial.
In this paper we have studied the quarkonium transverse-momentum-dependent fragmentation functions (TMDFFs) within the framework of non-relativistic QCD (NRQCD). We decompose the TMDFFs in a similar way as done in the corresponding collinear fragmentation functions [60]. We use the NRQCD factorization conjecture to write TMDFFs as a sum of products of short-distance matching coefficients and the standard NRQCD long-distance matrix elements (LDMEs).
Focusing on the case of J/ψ production in semi-inclusive deep-inelastic scattering (SIDIS) and for Q parametrically larger than the mass of the heavy quark, we use the standard TMD factorization to describe the cross section as a convolution of TMD parton distribution functions and the quarkonium TMDFFs. We provide the leading-order term in the perturbative expansion of the missing quarkonium TMDFF and use this result to make numerical comparisons against the competing photon-gluon fusion channel (through color-singlet intermediate QQ state). For this process, only  1 channel. We also included a short discussion on channel mixing through renormalization group evolution and phenomenological applications of the gluon TMDFFs at the LHC.
For our numerical implementation of the fragmentation cross section we use the public Artemide code [47], which we modified appropriately to incorporate the new quarkonium TMDFF. We discussed the kinematic regions of the future EIC and we provided two panels of plots for √ s = 63 GeV and √ s = 140 GeV. We showed that photon-gluon fusion is usually dominant at low values of Q and x, while quark fragmentation is the main initiating process in the opposite limit. The theoretical precision is at the moment limited only by the final quarkonium fragmentation part. In the case of photon-gluon fusion instead, a factorization theorem has not been yet formulated and its feasibility should be investigated in the future.  Figure 5. The quarkonium cross section coming from photon-gluon fusion in the color-singlet channel at leading order (γ * g) and TMD fragmentation (γ * q) for EIC kinematical settings √ s = 140 GeV.