Large-P_T inclusive photoproduction of J/psi in electron-proton collisions at HERA and the EIC

We study the inclusive J/psi production at large transverse momenta at lepton-hadron colliders in the limit when the exchange photon is quasi real, also referred to as photoproduction. Our computation includes the leading-P_T leading-v next-to-leading alpha_s corrections. In particular, we consider the contribution from J/psi plus another charm quark, by employing for the first time in quarkonium photoproduction the variable-flavour-number scheme. We also include a QED-induced contribution via an off-shell photon which remained ignored in the literature and which we show to be the leading contribution at high P_T within the reach of the EIC. In turn, we use our computation of J/psi+charm to demonstrate its observability at the future EIC and the EIC sensitivity to probe the non-perturbative charm content of the proton at high x.


Introduction
Quarkonium inclusive reactions have been widely studied at colliders since the discoveries of the first quarkonia in the 70's but they remain unsatisfactorily understood. We guide the reader to the following reviews covering the HERA and Tevatron legacy [1][2][3][4] and the current evolution of the field at RHIC and the LHC [5,6].
Inclusive J/ψ photoproduction, when a near on-shell photon hits and breaks a proton to produce the J/ψ, has been the object of several studies at HERA [7][8][9][10][11][12][13] in order to understand the quarkonium-production mechanisms and with the hope to learn more about the gluon content of the proton (see e.g. [14]). Indeed, the replacement of one hadron by a photon, in principle, simplifies the computations of such direct photon-proton interactions with respect to hadron-hadron ones. This also allows one to safely use quarkonia to extract Transverse Momentum Dependent gluon distributions [15][16][17][18].
Otherwise, this process is very similar to the widely studied hadroproduction where two hadrons collide to produce the quarkonium. Yet, the photoproduction cross sections are smaller and require significant luminosities to produce objects like quarkonia. The statistical samples collected at HERA are indeed more limited with quasi no data on ψ and none on the bottomonia.
Besides, at high energies, the hadronic content of the photon can be "resolved" by these reactions. The compu-tations relevant for these interactions between a resolvedphoton and a proton are then similar to those for hadroproduction. Moreover, the inclusive contribution is not always much larger than the exclusive or diffractive ones. Indeed, for vector mesons like the J/ψ, the corresponding cross sections can be extremely large. Both resolved-photon and diffractive contributions can be avoided with specific kinematical cuts as we shall see later. A last possible complication comes from the possible decays -also referred to as feed down (FD)-of excited quarkonia or b quarks, which should be analysed case by case like for hadroproduction.
In the present analysis, we focus on the large-P T region (P T M J/ψ ) which was seldom studied at HERA and which could be experimentally studied in more details at the future US Electron-Ion Collider (EIC) [19]. We will use the Colour-Singlet Model (CSM) [20], which is the leadingv contribution of Non-Relativistic QCD (NRQCD) [21]. Since the seminal work by Krämer [22] when he found that real-emission αα 3 s corrections associated to t-channelgluon-exchange topologies were large, we know that such NLO CSM corrections are in general large when P T increases. This can be traced back to the more favourable P T scaling of these specific real-emission contributions. The same has been observed in several hadroproduction processes of spin-triplet vector quarkonia [23][24][25][26][27][28]. The conclusion made by Krämer that the HERA data could be accounted for by the CSM after the inclusion of the NLO corrections were then questioned by 3 studies [29][30][31] in particular at large P T . We will revisit this for HERA in order to provide predictions for the EIC.
is that fragmentation cannot be initiated by g + g → g + g partonic reactions. Quarks are necessarily involved, which reduces their overall impact at HERA energies. In 1996, Godbole et al. [32] further showed that the charm fragmentation contribution was taking over the gluon fragmentation for P T > 10 GeV (see also [33]) which leads us to reconsider the contribution of J/ψ + (an unobserved) c to the inclusive production. Another J/ψ + c channel, where a cquark comes from the proton, had already been addressed in the early 80's by Berger and Jones [34] but was essentially forgotten. We also revisit it to show that it can be a relevant contribution to the inclusive yield but also that it can be studied for itself by further detecting a charmed hadron in the final state.
We also consider another overlooked contribution, namely from the pure QED process where the J/ψ is produced by an off-shell photon, which happens to be relevant for the large-P T domain. Owing to the presence of a quark line interacting with the initial photon, such a contribution should be as large as the typical fragmentation Colour-Octet (CO) contributions [35]. As a matter of fact, we see that this QED contribution matters when one reaches the valence region in the proton.
Thanks to different possible running energies, the EIC will indeed allow one to probe rather large x. This also motivate us to see to which extent extracting the J/ψ + c yield in this region could help measuring a valence-like non-perturbative charm content in the proton, also referred to as Intrinsic Charm (IC) [36].
The structure of our Letter is as follows. Section 2 outlines our methodology to compute the contributions of the aforementioned reactions, whereas Section 3 focusses on the FD from excited quarkonium states or b decays. Section 4 gathers our results for the inclusive cross section for the HERA and EIC kinematics, while Section 5 is devoted to the discussion of J/ψ+charm associated production. Section 6 gathers our conclusions.
2. J/ψ photoproduction at finite P T As announced, we will limit our study to the CSM. As such, at Born order, or LO in α s , J/ψ photoproduction proceeds via the partonic process γ + g → J/ψ + g thus at αα 2 s (Fig. 1a). Strictly speaking, γ+{c,c} → J/ψ+{c,c} (Fig. 1b) appears at the same order in α s as far as the matrix element is concerned [34].
Even though the contribution of the leading-P T fragmentation contribution have been discussed in a couple of studies [32,33], the QED fragmentation contribution (Fig. 1g), Figure 1: Representative diagrams contributing to direct J/ψ (unresolved) photoproduction via CS channels at orders αα 2 s (a,b), αα 3 s (c,d,e,f) and α 3 (g). The quark and antiquark attached to the ellipsis are taken as on-shell and their relative velocity v is set to zero.
whereby the photon strikes a quark off the proton which then radiates an off-shell photon turning into a J/ψ, has never been studied. It is indeed suppressed by the third power of α but should give a similar contribution than a similar graph with an off-shell radiated gluon into a 3 S [8] 1 CO. As it is expected to scale like P −4 T , it can be relevant at large P T . In addition, in the valence region, it may be enhanced by the larger quark densities compared to that of the gluons.
Qiu et al. recently reported [37] on the first NLO NRQCD study of the Q 2 -integrated cross section for J/ψ production at the EIC. Since the yield is expected to be dominated by photoproduction (very low Q 2 ) as opposed to leptoproduction (finite Q 2 ), we find it useful to stress the differences between both their and our studies. In fact, they started with Born-order (LO) contributions from 2 → 1 partonic processes, at αα 2 s , from CO states. Such LO contributions only lie at P T = 0 1 . As such, the corresponding virtual corrections only contribute at P T = 0 and the real-emission corrections generate the P T spectrum. In the photoproduction region, on which we focus here, P T P T . Out of their NLO contributions, the virtual corrections are thus cut away by the low-P T cut usually applied to the data (see section 2.2), whereas their real-emission CS corrections are in fact our αα 2 s γ+g → J/ψ+g contribution. Going further, our αα 3 s γ +{g, q,q} → J/ψ+{g, q,q}+g contributions would be part of a NNLO analysis of their Q 2 -integrated observable.

The CSM partonic amplitude
In the CSM [20,38,39], one computes the matrix element to create a 3 S 1 quarkonium Q with a momentum P Q and a polarisation λ accompanied by other partons, commonly noted { j}, by building the product of the amplitude to photoproduce the corresponding heavy-quark pair, M(γ + a → QQ + { j})), a spin projector N(λ|s 1 , s 2 ) and R(0), the radial wave function at the origin in the configura-tion space. More precisely, we have where P Q = p Q + pQ, p = (p Q − pQ)/2, s 1 and s 2 are the heavy-quark spins, and δ ii / √ N c is the projector onto a CS state. N(λ|s i , s i ) is a spin projector. In the non-relativistic limit, N(λ|s i , When one then sums over the quark spins, one obtains usual traces which can be evaluated without any specific troubles. In fact, such a computation can be automated at tree level as done by HELAC-Onia [40,41]

Photoproduction cross section in ep collisions
Let us start with some elements of kinematics. We first define s ep = (P e + P p ) 2 = 4E e E p (E e(p) is the electron (proton) beam energy) and s γp = W 2 γp = (P γ + P p ) 2 . Introducing x γ as P γ = x γ P e , one has s γp = x γ s ep .
Diffractive contributions are important at low P T and near the exclusive limit, i.e. when the J/ψ takes over the entirety of the photon momentum. A cut on P QT is usually sufficient to get rid of them. However, one can further cut on a variable called elasticity, defined as z = P Q ·P p P γ ·P p . z corresponds to the fraction of the photon energy taken by the J/ψ in the proton rest frame, with the proton momentum defining the z axis. It can equivalently be expressed as z = 2 E p m T W 2 γp e y in terms of the J/ψ rapidity y (with y and E p defined in the same frame) and the quarkonium transverse mass, m T = m 2 Q + P 2 QT . Diffractive contributions usually lie at z close to unity. On the contrary, the resolvedphoton contributions can become important at very low z where only a fraction of the photon energy is involved in the scattering.
All our computations are done using HELAC-Onia which is well tested for many colliding systems. Let us illustrate now with a generic 2 → 2 partonic reaction, γ + a → Q + k, the connexion between its amplitude, to be computed as discussed in section 2.1, and the ep (differential) cross sections which are the measurable quantities. The ep doubledifferential cross section in P T and z is then obtained via a convolution with the proton PDF and the photon flux from the electron : where . Following Eq. (5) of Ref. [42] under the Weizsäcker-Williams approximation, we have taken and m e is the electron mass. Besides the specific treatment of the J/ψ + charm yield, which we discuss below, the rest of the procedure remains completely standard. For 2 → 3 channels, the formulae would differ but the treatment is also absolutely standard.

Leading-P T approximation of the NLO corrections to
J/ψ + g: NLO In our evaluation of the NLO corrections to J/ψ + g, we will resort to the NLO approximation [23,24] which consists in considering the leading-P T contributions at NLO. These contributions are obtained by imposing a lower cut on the invariant mass of every pair of massless partons, s min i j . The NLO approximation avoids the (slow) computations of loop corrections and is easily interfaceable with Monte Carlo (MC) generators, which is advantageous for future EIC simulations. It has successfully been tested in hadroproduction of single J/ψ [24] and of J/ψ + J/ψ/γ/Z pairs [25][26][27][28].
The underlying idea of this approximation is that, for the leading-P T topologies associated to real NLO emissions, s i j for any i, j pair grows with P T and the result becomes insensitive to s min i j . Only the subleading-P T contributions -which normally cancel the IR divergences of the loop contributions-exhibit a log(s min i j ) dependence which is de facto P T -power suppressed. Such an approximation has never been applied to photoproduction whereas the exact same reasoning applies. To 3 check it, we have compared it in Fig. 2 (red band), in the HERA Run 2 kinematics for J/ψ photoproduction [12], to a complete NLO computation by Butenschön and Kniehl [43] (dashed black line). We have used the exact same parameters and PDFs as in Ref. [43]: 1 J/ψ = 1.32 GeV 3 and the CTEQ6M proton PDF set [44]. One observes a steady reduction of the band width from the variation of s min i j /m c within the interval [1,3]. At P T = 15 GeV, the corresponding uncertainty is already less than 30 % which is smaller than the other uncertainties that we will detail later. Setting s min i j /m c = 2 remarkably reproduces the full NLO result. In what follows, we will use this value for all our NLO predictions.
One has found out that, for the same s min i j variation at a given P T , the resulting uncertainty is slightly larger than for hadroproduction (J/ψ + X [24] and J/ψ + γ [25]). This can easily be traced back to the fact that the number of leading-P T real-emission graphs is twice smaller for photoproduction since there is only a single initial gluon which can radiate. As a case in point, one sees that the LO/NLO ratio is not small at P T ∼ 5 GeV when s min i j /P T is large. In fact, this likely reduces the s min i j (lower) uncertainty at P T ∼ 5 GeV since the LO does not depend on s min i j . This is why we have used a larger variation of s min i j /m c for our illustration: [1,3] instead of [1,2] in Refs. [24,25]. For the record, the quark channel contributions from γ + q → J/ψ + g + q are at the level of 10 ÷ 20% in this H1 kinematical range.

A Variable-Flavour-Number-Scheme treatment of the J/ψ+charm yield
Depending or not whether one considers the presence of charm quarks within the proton, J/ψ+charm production follows either from γ + g → J/ψ + c +c (3 Flavour Scheme (3FS)) at αα 3 s or γ + {c,c} → J/ψ + {c,c} (4 Flavour Scheme (4FS)) at αα 2 s . The situation is similar to single-charm production and both partonic subprocesses can be properly combined using a procedure referred to the Variable Flavour Number scheme (VFNS), outlined in Ref. [45].
A first application to quarkonium production was recently performed by one of us in Ref. [46] to treat J/ψ+charm hadroproduction. We go much further here in the photoproduction case considering this channel both from its contribution to the inclusive yield and to the J/ψ+charm yield. We will also study whether a change in the charm quark PDF, motivated by a possible IC component, can be visible in the overall combined cross section.
It is well known that both schemes should yield compatible results, yet they remain complementary in different phase-space regions. The 3FS is better suited close to the heavy-quark threshold with an exact kinematical treatment of the heavy-quark mass and the 4FS should give better results when the momenta exchanged are much larger than m c . Indeed, the evolution of the charm PDF in the 4FS allows one to resum logarithms of µ 2 F /m 2 c , whereas the 3FS only contains one of such logarithms. Like for bottom-quark-initiated processes [47], for increasing scales and x, the effect of these initial-state logarithms will grow.
The VFNS allows one to combine both schemes keeping their own virtues without double-counting contributions; it follows from the seminal work of Aivazis, Collins, Olness and Tung which lead to the introduction of the ACOT scheme [45]. This has been the object of many studies [48][49][50][51][52][53][54][55][56][57][58][59][60] but none related to quarkonia. However, at Born order, the matching procedure is straightforward and easy to implement provided that one properly considers the issue of the charm-quark mass in the 4FS part given the nonrelativistic nature of the J/ψ: m c cannot simply be set to zero.
Rather one should use x c s γp P γ . At LO, the double counting between the 3FS and 4FS contributions can be identified between, on the one hand, the initial-gluon splitting into the charm pair in the phase-space region where the charm quark interacting with the photon is nearly on-shell and, on the other, the PDF of this charm quark which also comes from a gluon splitting "inside" the proton. At LO, the latter reads wherẽ with the well-known Altarelli-Parisi splitting function P qg (z) = 1 2 z 2 + (1 − z) 2 . The overlap counter-term to be subtracted from the 3FS contribution is then Its computation can be handled by HELAC-Onia like the corresponding 4FS contribution.

Digression about the feed-down contributions
As mentioned in our introduction, the FD contributions from excited quarkonia or b quarks should be analysed case by case. The ideal situation is when these FD are experimentally determined. For instance, if the χ c and ψ yield are measured, one can deduce more or less directly what fraction of the J/ψ they contribute. In the case of inclusive photoproduction, such information is extremely scarce.
The fraction of non-prompt J/ψ was determined [12] for 60 < W γp < 240 GeV, 3 bins in z and P T > 1 GeV. In the bin 0.3 < z < 0.4 2 , it was found to be 14 ± 7%. No information about its P T dependence was given but a tuned Pythia analysis by H1 pointed at a possibly large value at the highest possible P T where the inclusive data were taken. Strangely enough, this was left unnoticed until a recent review [6] by one of us.
We have found it useful to revisit this by tuning Pythia 8.2 [61] using a b analysis by H1 [62] using di-electrons events which extends to large P T . Details are given in Appendix A. In agreement with the expectations of H1 [12], we have found the b FD to be significant around P T = 10 GeV which probably motivates a dedicated theory analysis at NLO with an updated treatment of the b → J/ψ hadronisation; this is however beyond the scope of our paper. In our comparison with the H1 data, we will show both the inclusive data (in black) and our expectation for the prompt data (in grey) which has a larger uncertainty but tends to significantly deviate from the inclusive one.
As what concerns the χ c FD, there is currently no theory or experimental indication that it could be relevant and it is therefore systematically neglected. As for ψ , after the last H1 analysis, ZEUS released in 2012 [13] the only information which we are aware of, namely the ratio ψ /J/ψ for 3 points in P T (up to 3 GeV, though). It was found to be compatible with a LO CS prediction, which directly follows from the ratio of the wave functions at the origin and from the ψ → J/ψ branching. In what follows, we will consider it to be on the order of 20%.

Results
In this section, we present and discuss our results for the inclusive prompt yield at HERA and the future EIC. At variance with the above comparison betwen the NLO and NLO , we use a more modern proton PDF set with uncertainties, namely CT14NLO [63]. We have used the same photon flux as above (see Eq. (3)). The corresponding theoretical uncertainty is evaluated automatically by HELAC-Onia using LHAPDF [64], like the factorisationand renormalisation-scale uncertainties from an independent variation in the interval 1 2 µ 0 < µ F , µ R < 2µ 0 , with µ 0 = m T . We have calculated the mass uncertainty by varying m c in the range m c = 1.5 ± 0.1 GeV for all but the QED channel. Indeed, for this channel, the invariant mass of the photon, 2m c , should coincide with M J/ψ and we chose m c = 1.55 GeV. We have also taken O We first naturally start by HERA for which we can compare our computations with existing data. Since our study aims at the large(r)-P T range where the QED and J/ψ+ charm contributions can be relevant, we will focus on the latest H1 data [12] extending out to P T = 10 GeV. Note that the latest ZEUS measurements [9] only reach P max T 5.5 GeV. 2 where the resolved-photon contributions are likely negligible. 3 using O  , is plotted in orange. Its ratio to the data is shown in the lower panel along with that of the LO. For the last 3 bins, we also show our expectation of the prompt yield after the subtraction of the b FD yield (see Appendix A).
About 10 years ago, it was advocated in a couple of works [29][30][31] that the CSM at NLO was unable to account for the HERA data at large P T at variance with the interpretation of NLO Krämer's computations [1,22]. We see here that the CSM can indeed describe the data. The introduction of J/ψ+charm contribution and the subtraction of the expected b FD improve (compare the red and orange bands) the agreement but the NLO agrees with the data within the uncertainties. We do not wish here to open the debate on whether there is room for additional contribution from CO transition and whether NRQCD NLO fits should be updated with a full study of the CSM uncertainties with modern PDFs and an estimation of b FD. Rather, in view of our discussion for the EIC kinematics, we would like to advocate that a CSM evaluation is certainly valuable in order to provide predictions for future EIC studies in terms of expected rates. In addition, we wish to stress that we anticipate the J/ψ+charm contribution to play a more significant role at large P T and maybe also the QED one in regions where the quark PDF would be less penalised.
We now move to the future EIC. This collider is expected to reach large integrated luminosities, at different √ s ep .
Here, we focus on two specific configurations, namely those with the following proton-electron beam energies: 100 GeV on 5 GeV and 275 GeV on 18 GeV. These respectively result in √ s ep = 45 GeV and √ s ep = 140 GeV.
To avoid the contamination from both diffractive/exclusive and resolved-photon contributions, we apply the following cuts: P T > 1 GeV and 0.05 < z < 0.9. As also done at HERA, we consider specific ranges in W γp . We stress that other choices could be made or that W γp could simply be integrated over. In fact, for the lower energy configuration, √ s ep = 45 GeV, where the cross section is a little smaller, we take a relatively broader interval 10 GeV < W γp < 40 GeV. For √ s ep = 140 GeV, we consider 20 GeV < W γp < 80 GeV. Finally, we impose a cut on the photon virtuality, Q 2 < 1 GeV 2 to remain in the photoproduction region. Fig. 4 shows the results for both the 45 GeV (a) and 140 GeV (b) EIC configurations. The colour code coincides with that of Fig. 3. On the lower panel of both plots, the ratio of each contribution to their total is shown. Let us start by first examining the lowest energy configuration at √ s ep = 45 GeV. First, we notice that one enters the valence region when P T increases. Indeed, the P −4 T QED contribution is no more suppressed by the quark PDF and becomes the leading one (and up to about 50% of the full yield) at the largest measurable 4 P T 10 GeV with an integrated luminosity of L = 100 fb −1 . Second, we note that the J/ψ+ unidentified charm contribution is comparable to the O(αα 3 s ) γg(q) fusion subprocesses. In such an EIC configuration, both new photoproduction contributions that were so far overlooked are thus absolutely relevant.
Moving now to the largest energy configuration at √ s ep = 140 GeV, one notes that the P T spectrum should be observable up to P T ∼ 15 ÷ 20 GeV with L = 100 fb −1 . Again, the O(α 3 ) QED contribution becomes dominant at large P T , being for instance more than 40 − 50% of the full yield at P T 15 GeV. At intermediate P T , the photongluon fusion prevails and thus J/ψ photoproduction in principle remains a good glue probe there. More generally, the production of J/ψ + 2 hard partons is dominant for 10 P T 15 GeV. Such a process could be identified by looking at J/ψ + 2 jets of moderate P T . Moreover, one of them should be relatively close to the J/ψ while the other should recoil on both. This clearly motivates dedicated studies to assess the feasibility of such J/ψ + 2 jet measurements. For a generic discussion about J/ψ + jets production, we guide the reader to [6]. 140 GeV (b). The calculation is performed adopting the same µ F , µ R and PDFs as Fig. 3 with the same meaning for bands.

J/ψ + (identified) charm
We finally move to the associated production of a J/ψ and a charm quark, leading to the observation of a charmed particle. In the past, unexpectedly large yields for the production of J/ψ and a cc pair in e + e − annihilation were measured by the Belle [67,68] collaboration at KEK. Such large yields could finally be explained within the CSM after the inclusion of α 3 s (NLO) corrections [69,70] with a K factors on the order of 1.4 ÷ 1.7. We guide the reader to Ref. [6] for a discussion of the different existing theoretical studies.
This naturally motivated theoretial studies of the corresponding process in the hadroproduction case [46,66,[71][72][73]. We extend these here to the photoproduction case with the potential of the future EIC in mind and use the same LO VFNS computation as for J/ψ + (unobserved) charm channels contributing to the inclusive J/ψ yield which we just discussed. At variance with the results presented in Section 4, we now consider a charm tagging efficiency, which we dub ε c and which then enters the combination of the 3FS and 4FS contributions (see Section 2.4) meant to avoid any double counting: Based on EIC simulations [74,75], we take ε c = 10%. It is a first working hypothesis which includes possibilities to tag the (anti)charm via charmed hadrons or by tagging a charm jet. In order to also study the possible impact of IC, we adopt the CT14NNLO proton PDF set [76], which includes PDF eigensets with some IC effects, namely a "sea-like" one (in green) and a valence-like one, also called "BHPS", (in red). The central eigenset include "no IC" (in blue). We stress that this only affects the combined cross section via the 4FS contribution through Eq. (8) 5 . Fig. 5 shows the results for the J/ψ+charm associated production at the EIC at √ s ep = 45 GeV (a) and √ s ep = 140 GeV (b). On the lower panel of both plots, the ratio to the "no IC" case is shown and compared to the "no IC" relative uncertainty. One can first note that at √ s ep = 45 GeV, even with L = 100 fb −1 , the yield is limited to P T below 5 GeV. On the contrary, the P T reach extends up to P T ∼ 10 − 12 GeV at √ s ep = 140 GeV. In the former configuration, the BHPS valence-like (in red) enhancement is visible, being as large as 4 − 5 times the "no IC" yield at P T ∼ 4 − 5 GeV. More importantly, this deviation is larger than the expected uncertainty (also in blue in the lower panel) of the "no IC" yield. At √ s ep = 140 GeV, one does not probe the valence region and such a deviation becomes too small compared to the "no IC" uncertainty. The EIC at √ s ep = 45 GeV will thus be the place to probe the non-perturbative charm content of the proton with J/ψ+charm production. Besides, it would also be interesting to know how precisely the charm and anticharm contributions to this process could be experimentally separated out in order to look for an IC induced c(x) vsc(x) asymmetry as recently found from lattice QCD studies [77]. A further motivation to study this process is to investigate on colour transfers [78], that can enhance the cross section at low invariant masses of the J/ψ+charm pair, M J/ψ+c . Such a phenomenon should be both measured in photoproduction and hadroproduction.

Conclusions
In this Letter, we have analysed the inclusive, large-P T J/ψ photoproduction at ep colliders. Such a process is a key 5 except for some possible effects on the gluon PDF via the momentum sum rule. with "no IC" (blue), sea-like (green) and "BHPS" valence-like (red) ICs. The lower panel shows the ratio to the "no IC" curves and its relative uncertainty. The kinematical cuts are the same as in Fig. 4. A 10% charm detection efficiency is considered when drawing the horizontal observability lines. measurement to improve our understanding of the production mechanism of quarkonia and a potential good probe of gluons. Our analysis incorporates the leading-P T leading-v next-to-leading α s corrections, i.e. the NLO CSM. In addition, we have studied for the first time the O(α 3 ) QED and the J/ψ+charm contributions, which both also contain leading-P T topologies. We have computed the J/ψ+charm cross section at LO in the VFNS, which properly combines the O(αα 3 s ) 3FS and O(αα 2 s ) 4FS contributions with specific counter-terms, avoiding double counting.
By comparing our NLO cross section to a complete NLO calculation at HERA, we have further validated this leading-P T approximation, already used in hadroproduction for a few processes. We have then reassessed the latest H1 data. It turns out that, contrary to earlier claims, the CSM incorporating NLO leading-P T contributions can describe J/ψ photoproduction at HERA for P T > 5 GeV. The agreement with the data slightly improves when the J/ψ+charm contribution is included, while the LO QED contribution would only play a role at larger P T (on the order of 20 GeV). We have also confirmed with a simple data-driven analysis that the b FD should probably not be ignored as it has so far been the case. When approximately subtracted from the data, the data-CSM agreement indeed further improves.
We have then studied the case of a future EIC which will be able to run at high luminosity (up to 100 fb −1 for one year of data taking) with different √ s ep energies. As such, the EIC should have a larger P T reach. We have provided predictions for the inclusive cross section for prompt J/ψ photoproduction at both the lowest and highest √ s ep configurations, namely 45 GeV and 140 GeV. At variance with HERA, the P −4 T QED contribution turns out to be dominant at large P T , becoming about half of the prompt yield at P T 10 (15) GeV for the √ s ep = 45 (140) GeV configuration. This is due to the prominent role of the quark in the valence region. Furthermore, we have found that the J/ψ+c contribution is on the same order as the O(αα 3 s ) γg(q) subprocesses at √ s ep = 45 GeV. It is thus important to note that J/ψ production at the lower energy configuration of the EIC is no longer a gluon probe at mid and large P T .
The O(αα 3 s ) γg(q) processes can however be disentangled by looking at a J/ψ with 2 jets of moderate P T .
Applying the VFNS to the J/ψ+charm production, we have also demonstrated its observability at the EIC. At √ s ep = 45 GeV, a "BHPS" valence-like behaviour can be observed at P T > 5 GeV, with a possible enhancement as large as 4 − 5 compared to a "no IC" hypothesis . Indeed, in this energy configuration this processes is sensitive to the charm quark in the valence region. This is not so at √ s ep = 140 GeV, where one cannot a priori discriminate between the "IC" and "no IC" hypotheses. The EIC at √ s ep = 45 GeV thus offers a unique opportunity to better constrain the charm quark content via the associated production of a J/ψ with an identified charmed particle. In the future, such studies could surely be extended to other planned colliders, such as the LHeC and FCCeh [79,80]. In addition, dedicated experimential studies for the associated production of J/ψ+jets and FD from heavier quarks or quarkonia at the EIC would be surely beneficial.
In this appendix, we briefly explain the data-driven procedure we have used to estimate the b → J/ψ FD at HERA. We have considered the latest H1 analysis on bb production [62] via di-electron events. The b photoproduction cross section was measured in four P T (b) bins, with P T (b) = P 2 T,b + P 2 T,b /2, in the range [0, 30] GeV. It is important to note that H1 reported on a b, not on a B cross section. We have thus computed the corresponding LO plus parton shower (LO+PS) cross section using Pythia 8.2, adopting the CT14NLO PDF set with m b = 4.75 GeV and µ F = µ R = m T . Note that we have not considered any QCD uncertainty since we have assumed that they are fully correlated between both the H1 kinematics of the b measurements and that of the J/ψ measurements. The corresponding kinematical ranges are indeed very similar. Such theory uncertainties will thus be absorbed in the tuning factor which we discuss now. We have indeed performed a χ 2minimisation to compute a tuning factor, denoted N FD , such that the obtained LO+PS Pythia spectrum reproduce best the H1 b data. We have found N FD = 3.8 ± 0.5. We note at this stage that N FD partially takes into account higher-order effects missing in the LO+PS computation. The resulting cross section is shown in Fig. 6 (a) with the H1 values.
This tuning of the b cross section being made, we have again used Pythia 8.2 to compute the b → J/ψ cross section in the H1 kinematics [12]. Indeed, we expect Pythia 8.2 to account reasonably well for the fragmentation of a b quark into a J/ψ. The obtained cross section is plot in red in Fig. 6 (b). Its uncertainty is purely from that of N FD , thus from the H1 uncertainty of the reported H1 b cross section. As explained above, we have assumed that the theory uncertainties are absorbed in the tuning factor. We have then subtracted this b → J/ψ yield from the inclusive one in order to obtain a data-driven estimate of the prompt cross section. These are respectively shown in Fig. 6 Fig. 3. (black and grey points). To conclude, we note that our data-driven procedure yields a fraction of non-prompt J/ψ slightly smaller, yet compatible with the P T -integrated value reported by H1 [12] for 60 < W γp < 240 GeV for 0.3 < z < 0.4. This means that our data-driven estimate can be seen as a conservative low value of the b FD to be subtracted at large P T .  [62] and our LO+PS computation. The blue (resp. red) histograms depicts the "untuned" (resp. tuned) results. The dotted histograms indicates the uncertainty from the χ 2 -minimisation via N FD = 3.8 ± 0.5. (b) The measured inclusive (black points) and derived prompt (grey points) H1 J/ψ cross sections [12] along with the derived b → J/ψ one (red) obtained after tuning Pythia (red).