Exclusive heavy vector meson electroproduction to NLO in collinear factorisation

We compute the exclusive electroproduction, $\gamma^* p \rightarrow V p$, of heavy quarkonia $V$ to NLO in the collinear factorisation scheme, which has been formally proven for this process. The inclusion of an off-shell virtuality $Q^2$ carried by the photon extends the photoproduction phase space of the exclusive heavy quarkonia observable to electroproduction kinematics. This process is relevant for diffractive scattering at HERA and the upcoming EIC, as well as at the proposed LHeC and FCC.


Introduction
The exclusive production of vector mesons has long been an interesting and attractive observable to study. First measured in the fixed target mode and then in diffractive deep inelastic scattering (DIS) events at the ep linear HERA collider more than 25 years ago, they constitute ∼ 10% of the total inclusive DIS cross-section and are characterised by the presence of a large rapidity gap. They provide a means to investigate the phenomenology of quarkonium production and function as more sensitive probes of the low-x and low scale input gluon parton distribution than any other known high-energy physics phenomenon.
In 1993, around the same time as the first measurements of such diffractive activity in a collider environment, the exclusive electroproduction of a heavy vector meson (HVM), V = J/ψ, Υ, ψ(2S), . . . via γ * p → V p, was showcased to be proportional to the square of the gluon parton distribution function (PDF) [1] in the leading logarithmic approximation (LLA) of perturbative QCD (pQCD), within the framework of k T -factorisation. See [2] for a review. This coincides with the leading-order (LO) term in the collinear factorisation scheme, which we use here. This process is on solid ground thanks to a proven factorisation theorem [3] and the pQCD treatment is justified for large virtualities Q 2 of the photon, γ * .
On the experimental side, the HERMES collaboration [4] reported leptoproduction measurements for the lightest vector mesons in the range 1 GeV 2 < Q 2 < 7 GeV 2 in fixedtarget kinematics. Exclusive electroproduction data for the J/ψ HVM via dimuon and dielectron decays has been measured in the collider mode at HERA in a narrow range of photon virtualities at both ZEUS and H1 experiments, extending up to the largest bin of Q 2 = 22.4 GeV 2 [5,6]. As in photoproduction, the cross section exhibits a steep rise with increasing centre of mass energies of the γ * p → J/ψp subprocess. Today, in the LHC era of collider physics, central exclusive photoproduction of vector mesons V have been measured in the forward rapidity interval 2.0 < Y < 4.5 by the LHCb collaboration via ultraperipheral pp → p + V + p collisions instead [7][8][9]. These are driven by the hard scattering subprocesses γp → V p, measured directly at HERA. Here, the photoproduction reaction (Q 2 ≈ 0) is initiated by a real on-shell photon, γ. Despite the near vanishing of this scale, the factorisation theorems are still assumed to hold for photoproduction since the masses of the produced final state heavy mesons are above the perturbative scale.
Various theoretical models within pQCD exist in the literature that provide a description of the exclusive heavy vector meson photo and electroproduction processes [2]. In the colour dipole approach, the exclusive HVM formation is dominated by scatterings in which the photon fluctuates into a qq pair with a transverse separation r ≈ 0, carrying fractions z ≈ 1/2 and 1 − z ≈ 1/2 of the incoming photon momentum. The dipole model formulation is also able to describe light meson production and photon hadron scattering and is equivalent to the k T -factorisation formalism in the leading ln(1/x) approximation. Following earlier work, in [10], the explicit k T integral had been performed in the last cell of evolution, in effect leading to a description beyond the LLA, to mimic a subset of the full next-to-leading order (NLO) contribution. This accounts for those terms in the conventional DGLAP evolution that are enhanced at small x and that in an axial gauge correspond strictly to gluon-ladder-rung Feynman diagrams. However, we emphasise this does not encompass the complete NLO contribution that one would obtain in the conventional systematic evaluation of Feynman diagrams within the MS collinear factorisation framework.
In this work, we remain entirely within the collinear factorisation set-up at NLO to extract the electroproduction renormalised transverse and longitudinal coefficient functions to NLO in the MS scheme. Previously in the literature [11,12], next-to-leading order MS coefficient functions were calculated in the case of photoproduction of HVMs. The authors of [11] constructed the imaginary part of Feynman diagrams via Cutkosky-cuts in the schannel and then restored their real parts using the corresponding u-channel contributions, via a dispersion relation. In our approach for electroproduction, we directly compute the real and imaginary parts of the amplitude using an integral reduction procedure. As will be discussed and explained, the zero photon virtuality limit of our electroproduction coefficient functions should coincide with these photoproduction results. Note that a subset of the authors in [11] also computed the electroproduction of light neutral vector mesonsṼ = ρ 0 , ω and φ [13]. In our computation, both the virtuality of the photon and the mass of the heavy quark constitute massive scales, which add complexity.
The work presented here is the first step in driving the prediction of exclusive heavy vector meson electroproduction to NLO. As such it represents the technical details of the underlying field theory construction. In the case of photoproduction further steps, including scale fixing and a Q 0 subtraction, were required to improve the convergence of the perturbative prediction before the results could be applied in phenomenological studies. For exclusive J/ψ photoproduction these steps were carried out in [14][15][16][17] and used to present predictions which described data from HERA and LHCb and to determine the gluon distribution at low x and scale. Going forward, our electroproduction results will be useful for studying exclusive processes at the upcoming EIC [18] and at the proposed LHeC [19] and FCC [20], which we leave for future work. Our paper is organised as follows. In Section 2, we give our set-up and model assumptions within collinear factorisation at NLO. In Section 3, we outline the workflow of our calculation. In Section 4 and Section 5, we give analytic expressions for the LO and NLO coefficient functions for the quark and gluon initiated subprocesses. We finish, in Section 6, by checking the explicit cancellation of collinear divergences to NLO within a consistent UV and IR subtraction scheme, before making a comparison with literature in Section 7 and presenting our conclusions in Section 8. We describe the matrix element for exclusive HVM electroproduction as the fluctuation of a hard incoming photon with momentum q µ + ∆ µ /2 into a heavy QQ pair, which then interacts with the proton (or nuclei) carrying momentum P µ − ∆ µ /2 via a two-parton colour singlet exchange mechanism, see Fig. 1. The proton recoils slightly with momentum P µ + ∆ µ /2. The modelling of the open quark-antiquark recombination into the observed exclusive final state HVM with momentum q µ − ∆ µ /2 is made, as in [11,12], via LO Non-Relativistic-QCD (NRQCD). In this approach, the amplitude for the production of two onshell heavy quarks is calculated and projected onto the outgoing HVM quarkonium state.

Kinematics and set-up
The amplitude can be expanded in powers of α s and the heavy quark relative velocity. Here, we compute the α s corrections. The relativistic corrections have been studied elsewhere, see [21]. To LO in the NRQCD relative velocity expansion, the momenta of the quark and anti-quark are equal such that their sum equals the momentum of the HVM.
Following the set-up in [22], the three independent momenta defined above are decomposed in terms of high energy light-like Sudakov basis vectors {p, n, ∆ ⊥ } , satisfying p · p = n · n = 0 and p · n = 1. The mean of the incoming and outgoing proton momenta, P µ , defines the collinear direction. In the Bjorken limit to leading-twist accuracy, i.e. neglecting the masses of the nuclei, the kinematics of the process simplify and can be expressed in terms of the virtuality of the photon, Q 2 , the mass, M , of the HVM, and the skewedness parameter, ξ. Here, P µ ≈ p µ and ∆ µ ≈ −2ξp µ , where 2ξ is the 'kick' which the active quark or gluon receives along the collinear direction so that the t-channel momentum exchange, t = ∆ 2 = 0. The probed partons (gluons or quarks) carry momenta p 1 = (x+ξ) p and p 2 = −(x − ξ) p, the momentum fraction x is integrated over in the convolution with the Generalised Parton Distribution functions (GPDs). Inserting the replacements allow us to reduce the number of dimensionless variables appearing in our description by one. Note that these transformed basis vectors still respectp ·p =n ·n = 0 andp ·n = 1 and so all possible scalars formed from our momenta are unaffected by this change. Therefore, to leading order in the relative heavy quark velocity and in the Bjorken limit, the Sudakov decomposed momenta of Figure 2 where Our convention is that all momenta are incoming. Moreover, in accordance with the leading term in the NRQCD expansion, we make the approximation that the on-shell pole mass of the heavy quarks is m = M/2. At leading order only the gluon induced process contributes; At NLO the QQ pair may scatter from a light quark via a gluon exchange and so a new partonic channel opens. We compute in addition the quark induced process: (q = u,ū, d,d, s,s). (2.5)

HVM and GPD spin projections
The S-wave, spin-triplet projection may be written to leading order in the heavy quark relative velocity as [23][24][25] v where, in the non-relativistic limit, we take the vector meson momenta K = 2p 3 = 2p 5 and mass M = 2m. Hereū j β is the outgoing heavy quark spinor and v i α is the outgoing heavy anti-quark spinor. The indices i and j label the colour of the heavy quarks whilst α and β label their spin. O 1 V represents the non-perturbative NRQCD matrix element. The vector S describes the polarisation of the HVM; it satisfies S · * S = −1 and K · * S = 0. Here relative to [25] we have an overall minus sign. It multiplies the overall amplitude and so has no effect on the cross-section.
The quark GPD contraction is implemented as a spin projection of the on-shell quark scattering matrix. We replace the spinors of the quark and anti-quark connecting to the quark GPD, F q (x, ξ), by The factor / p αβ results in a trace over the spin line of the quarks connecting to the GPD at the amplitude level. This can be understood by considering the numerator of a quark diagram. Representing the product of quark propagator numerators as [A] and applying the spin projection we obtain Similarly, for the gluon induced partonic channel, the contraction with the gluon GPD, F g (x, ξ), is implemented as a projection of the on-shell gluon scattering matrix. Throughout we will use dimensional regularisation in d = 4 − 2 space-time dimensions and make the replacement where µ 1 and * 2 ν are the polarisations of the incoming and outgoing gluons respectively. Here, The correct iδ prescription for the poles has been discussed extensively in the literature, see for example [26]. The prescription given here is valid for both DVCS and HVM production [11,27]. The indices a, b are gluon colour indices in the adjoint representation. The factor 1/(N 2 c −1) averages over the gluon colours. The 1/(d − 2) factor is our analytic continuation of the number of transverse polarisations of a gluon in d dimensions. It appears due to the average over the gluon polarisations. The factor of 1/2 is required to prevent double counting when both s and u channel gluon diagrams are computed (as done here) and the momentum fraction x is integrated over from −1 to 1 (see later). The object g µν ⊥ ≡ g µν − p µ n ν − p ν n µ carrying the gluons' Lorentz indices is the perpendicular metric tensor which projects onto the plane perpendicular to p and n.

Lorentz-invariant tensor decomposition
We consider only the vector part of the amplitude at leading twist and at t = 0. Highertwist terms are formally suppressed and axial-vector contributions are neglected, as they are not needed here with an unpolarised nucleon in the initial state. As shown in Section 2.1, in the Bjorken limit at leading twist, all of our external kinematics can be expressed in the Sudakov basis {p, n} with ∆ ⊥ = 0. We decompose the part of the amplitude insensitive to the helicities of the incoming partons in the nucleon target in terms of the available Lorentz structure in this basis. Explicitly, we factor off the polarisation vectors for the incoming photon and outgoing HVM and work with the amputated amplitude T (µν) . 2 It follows that where A, B, C, D and E are the arbitrary coefficients of the decomposition. Imposing local current conservation at the photon vertex, p 4,µ T (µν) = 0, together with the identity K ν T (µν) = 0, constrains the coefficients. The former is the familiar Ward-identity while the latter is, strictly speaking, not but holds at the Feynman diagram level and is true due to our choice of HVM spin projection. With these equations, we obtain the decomposition (2.14) An explicit exposition of the available Lorentz structure has therefore produced a manifest decoupling of the system into two overarching degrees of freedom, parametrised by T ⊥ and T L . Contractions of (2.12) with explicit realisations of the physical transverse and longitudinal polarisation vectors of the photon ε γ µ and HVM ε V * ν pick out one of the two scalar coefficients in each case. This may be seen as follows. The transverse polarisation vectors have only a transverse component in their Sudakov-basis decomposition. By construction, their contraction with p µ , p ν , n µ , n ν in µν vanishes and the only contribution is due to g µν ⊥ : where ε γ ±,µ and ε V * ±,ν are the transverse polarisation vectors for the photon and HVM in the helicity basis, respectively. The corresponding longitudinal polarisation vectors are satisfying ε γ L · ε γ L = −1 and p 4 · ε γ L = 0, with similar relations for the HVM. Then, where there is now no contribution from g µν ⊥ and where we have fixed N in (2.13) such that the contraction of µν with the physical longitudinal polarisation vectors is unity.
The decomposition (2.11,2.12) is therefore readily identifiable as a separation into transverse and longitudinal components, with T ⊥ and T L having the physical interpretation now as the process's transverse and longitudinal form factors, respectively. With this choice 3 , the transverse and longitudinal helicity amplitudes, A ±± = ε γ ±,µ ε V * ±,ν T (µν) and A 00 = ε γ L,µ ε V * L,ν T (µν) , are equal to T ⊥ and T L , respectively. Note also that the introduction of the N factor into µν in (2.13) allows for both two-tensors in multiplication with the scalar coefficients T ⊥ and T L to have mass dimension zero. In this way, one may extract these coefficients in turn through suitable projections onto the T (µν) structure.
We remark that (2.11) coincides with the leading-twist tensor decomposition found in Generalised-Deeply-Virtual-Compton-Scattering (GDVCS), see e.g. [28], upon neglecting the axial-vector and helicity flip contributions. This is to be expected as the only distinction in the kinematical set-up is the final state production of a heavy photon instead of a heavy vector meson that we have here, however this remains indifferent to the construction of the underlying tensorial structure and the applicability of our Ward and Ward-like identities.
The vector part of the amplitude may be written, using collinear factorisation, as where the ellipses represent contributions that only appear beyond the leading twist term and outwith the chiral-even theory with t = 0, but which would appear in e.g. polarised scattering or if the nucleon mass would not be neglected. The renormalised quark and gluon GPDs are denoted F q and F g respectively. C ⊥,q and C ⊥,g are the renormalised quark and gluon vector transverse coefficient functions while C L,q and C L,g are the renormalised quark and gluon vector longitudinal coefficient functions. The dependence of the GPDs F q and F g on the factorisation scale µ F and on t has been suppressed. The dependence of the renormalised coefficient functions C ⊥,q , C ⊥,g , C L,q and C L,g on m, µ F and the renormalisation scale µ R has been suppressed, too. We recall that the Lorentz indices µ and ν are those of the incoming photon and outgoing HVM respectively. Measurements of HVM production from unpolarised targets probe only the charge conjugation even quark GPD; we may therefore replace q F q with the singlet quark GPD F S /2. There is also an additional photon helicity flip term which we neglect.

Overview of the calculation
We generate all LO and NLO Feynman diagrams using QGRAF [29] and select those diagrams that are compatible with our external colour and kinematical constraints. Example Feynman diagrams for the LO and NLO quark and gluon initiated subprocesses are shown in Fig 3. Each selected diagram is converted into an expression through insertion of Feynman rules, derived in an arbitrary linear covariant gauge. The appropriate GPD quark or gluon projector is applied to each diagram, together with the HVM spin projection, and the resulting d-dimensional Dirac traces are computed and handled in FORM4.2 [30]. As our amplitude is written with two free Lorentz indices, µ and ν, we need at most a basis decomposition for the rank-two tensor k µ k ν , where k is the loop momentum. Due to the external colour and kinematical constraints, the integral structures obtained contain in general linearly dependent propagators. They therefore cannot be reduced directly using standard integral reduction tools. We express them first as structures containing only linearly independent propagators by applying a partial fractioning routine in line with the Leinartas algorithm [31]. In this way, we eliminate the linear dependence amongst the propagators and proceed with the integral reduction via REDUZE2 [32] which encodes Laporta's integration by parts algorithm [33].
All results presented below are given in the ERBL region, |ξ| > 1 (i.e. |x| < ξ), where we expect an absence of imaginary parts. The results in the physical region, |ξ| < 1 (i.e. |x| > ξ), may be obtained by restoring the correct analytic continuation. Explicitly, this iŝ The coefficient functions are expanded in the bare strong coupling a 0 = α 0 /(4π) as, with i =⊥, L and j = q, g.

LO results
At leading-order (LO) in α s , there is only a gluon initiated subprocess. There is no contribution from an initial state light quark as it only couples to the outgoing heavy quark system via a loop-induced gluon insertion at next-to-leading-order (NLO). A single gluon exchange is ∝ f abc , which vanishes due to our GPD spin projection ∝ δ ab and antisymmetry of the SU(3) structure constants. The tree level contribution is therefore dominated by a gluon-gluon (pomeron-like) exchange: With the expansion in the bare strong coupling defined in (3.2), the bare (denoted by hat script) transverse and longitudinal coefficient functions at LO arê Here, where g e (g s ) is the electromagnetic (strong) coupling, e q is the photon-quark charge and O 1 V is the NRQCD matrix element. The colour factor T F = 1/2, and There is a cancellation of terms between those appearing in the 1/(d − 2) factor of eqn. (2.9) and those in the numerators of the Feynman diagrams. This observation is important as the tree level results enter into the renormalisation of the NLO result. We will therefore not generate finite surplus contributions of the form / in the computation of the NLO counterterms, see Section 6.
The Q 2 → 0 (i.e. w → ∞ ) limit of these results coincide with the photoproduction result, [11,12], where only the amplitude to produce a transversely polarised HVM is non-vanishing. 4 Note that the limiting point w → 0 is not kinematically attainable. The production of a time-like vector meson (with invariant mass squared, M 2 = 4m 2 > 0) initiated by a space-like photon, Q 2 < 0, does not, as expected within the analytic structure of the S-matrix, produce a pole at LO.

NLO results
At NLO, there are quark and gluon initiated subprocesses. We compute and express our bare quark and gluon NLO coefficient functions in terms of sets of coefficients c i and a universal basis set {f i , i = 1, . . . , 12} of logarithms and dilogarithms that arise from the Feynman loop integrations. They are Here, v =ξ/w 2 and so, with this choice of variables, the analytic continuation to the DGLAP regime, as specified in eqn.

Quark subprocess
The quark q subprocess amplitude, A q , to NLO can be written in the form a .

(5.5)
Here, χ 1,2 are the HVM and GPD projectors and Γ (1,2) a encode the Dirac structure for each spin line in diagram a. The appearance of a trace at amplitude level is particularly noteworthy and is a reflection of our external colour and kinematical constraints.
Each of the four quark diagrams where the external photon couples to the open heavy quark or antiquark lines may be decomposed into a tadpole and two bubble master integrals that emerge from the Laporta algorithm. When the photon instead attaches to the heavy quark mass propagator, there are in addition two triangle master integrals in the diagram decomposition. In this latter case, we have decomposed a pentagon integral, containing linearly dependent propagators, into a sum of triangles, bubbles and the tadpole integral which, by construction, contain only linearly independent propagators.
We normalise our quark NLO coefficient function with the factor arises from our loop integrals and the scale µ 0 is a mass parameter introduced in dimensional regularisation to maintain a dimensionless bare coupling. Here, γ E is the Euler-Mascheroni constant. The bare transverse quark NLO coefficient function can be written in the form where the f i are given in eqn. (5.3) and the c ⊥,i are vw 2 + 32v, c ⊥,11 = 0, c ⊥,12 = 0. Similarly, the bare longitudinal quark NLO coefficient function can be expressed aŝ where the f i are given in eqn. (5.3) and the c L,i are

Gluon subprocess
We normalise our gluon NLO coefficient function with the factor where C is defined in eqn. (5.7). The bare transverse gluon NLO coefficient function can be written in the form where the f i are given in eqn. (5.3) and the c ⊥,i are We have checked that the Q 2 → 0 (w → ∞) limit of the above agrees with the expression given in [11,12] for the photoproduction set-up. We recall that v =ξ 14) The bare longitudinal gluon NLO coefficient function can be written in the form where the f i are given in eqn. (5.3) and the c L,i are Here, C F and C A are the fundamental and adjoint Casimirs of SU(N c ). With N c = 3, C F = 4/3 and C A = 3. The diagram group theory factors recur in three combinations: The coefficients of C F and C A are polynomially reduced and presented above using the MultivariateApart package [34]. The physical symmetryξ → −ξ is again made apparent in our presentation of the NLO gluon coefficient functions.

UV renormalisation and mass factorisation
We renormalise the gluon, heavy quark field and heavy quark mass in the on-shell (OS) scheme. The strong coupling constant is renormalised with light flavours treated in the MS scheme and with the heavy quark loop of the gluon self-energy subtracted at zero momentum. The UV renormalised amplitude A UV may be written in terms of the bare amplitude A using the relation 5 , Here a 0 = α 0 /(4π) and a s = α s (µ 2 R )/(4π) are the bare and renormalised strong couplings, respectively, and we have introduced S = (4π) e − γ E . The heavy quark bare mass parameter is denoted m 0 and the renormalised quark mass is denoted m. The gluon and heavy quark renormalisation constants are denoted Z A and Z 2 and the number of external gluons and heavy quarks are denoted n g and n q , respectively. In the second line the bare amplitude and renormalisation constants are expanded using We have also introduced the mass counterterm amplitude A mct which is computed by inserting mass counter vertices and cross-checked by computing the derivative of the bare amplitude with respect to the heavy quark mass, m. The derivative must be taken prior to equating the HVM mass, M , to twice the heavy quark mass in the HVM projector.
The explicit expressions for the renormalisation constants at one-loop are Although the above renormalisation procedure is stated for the bare amplitudes, it can equally be applied to each of the bare coefficient functionsĈ i=⊥,L, j=g,q . The UV renormalised coefficient functions are thus given bŷ The relevant mass counterterm coefficient functions arê Since the quark amplitudes enter only at a 2 s they are trivially modified by the above procedure.
After UV renormalisation, poles in still remain in both the quark and gluon oneloop coefficient functions and must be absorbed into the definition of the GPDs via mass factorisation. This procedure generates counterterms which render the coefficient functions finite. Concretely, we replace the bare quark singlet (F S ) and gluon (F g ) GPDs with the mass factorised, factorisation scale (µ 2 F ) dependent, GPDs, (6.10) Here V (1) is the coefficient of α s /(4π) in the generalised splitting function V . We take the generalised splitting functions from Ref. [35]. Inserting these relations into the bare version of the factorisation formula (2.19) and absorbing the divergent terms into the coefficient functions we obtain where the mass factorisation counterterms are given bŷ The authors of [36] have also computed the transverse and longitudinal renormalised amplitudes at NLO for the exclusive electroproduction of heavy quarkonia, within the same calculational framework of collinear factorisation and LO NRQCD. Our work therefore serves as an independent, simultaneous analysis of the process and, in addition, a check on the results already available. We have obtained agreement numerically with both their quark and gluon transverse and longitudinal renormalised amplitudes after publication of their erratum, in which typographical errors were corrected and an inconsistency resolved in their high energy limit (ξ → 0) that we had pointed out. Interestingly, in the high energy limit the dominant logarithm ln(1/ξ) is multiplied by ln Q 2 /µ 2 F whereQ 2 = (Q 2 +4m 2 )/4. Therefore, adopting the scale choice µ 2 F =Q 2 would automatically resum these doublelogarithmic terms ∼ α 2 s ln Q 2 /µ 2 F ln(1/ξ) and effectively set the NLO corrections to zero in the high energy limit.
Note that our results are expressed in terms of a basis set comprising a smaller number of logarithms and dilogarithms. The set we have selected is such that all basis functions appearing are manifestly real if the coefficient function is real. Throughout the entire ERBL phase space, all of our basis functions evaluate to a real number and, in addition, only one square root argument is not rationalised. The substitution z 2 = (1 + 2v − vw 2 )/(1 + vw 2 ) would allow all arguments of the functions appearing in our basis set to be rationalised, however this would upset our explicit v → −v symmetry which we deem to be a more important feature to maintain. We have also decided to keep track of all the group theory factors that arise, allowing us to generalise the expressions for N c = 3 and be valid for an arbitrary SU(N c ) gauge theory instead.

Conclusions
In this work, we have computed renormalised coefficient functions for exclusive electroproduction of heavy vector mesons to NLO in the collinear factorisation framework. The description of the QQ → V transition vertex was made within LO NRQCD and the mass of the heavy vector meson set equal to twice the on-shell mass of the heavy quark. Our results respect the physical ξ → −ξ symmetry and coincide with the renormalised coefficient functions for exclusive photoproduction of heavy vector mesons in the Q 2 → 0 limit. They can be used in a phenomenological analysis of the exclusive electroproduction data already measured at HERA and, in time, with the data from LHC, the upcoming EIC and the proposed LHeC and FCC. These results may also help in studying the onset of saturation physics. Importantly, our predictions and these data can collectively provide further constraints on the gluon distribution in nuclei at moderate to low values of the scale Q 2 and Bjorken x.