Probing gluon GTMDs through exclusive coherent diffractive processes

We extend a previous GTMD model to improve the description of the HERA data on diffractive dijet production, and include exclusive coherent diffractive J/ψ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J/\psi $$\end{document} production data. We find that within our gluon GTMD model context and assumptions, there is considerable tension between the data for these two types of processes concerning the t dependence. Photo- and electroproduction data for protons and nuclei from EIC and UPC data from LHC and RHIC can help to establish whether a common GTMD description is possible, as one would expect, and to facilitate studies of such data we provide predictions for the various experiments. We point out explicitly in which sense this goes beyond the description in terms of GPDs.


I. INTRODUCTION
Coherent diffractive processes form a promising way to probe Generalized Transverse Momentum Dependent parton distributions (GTMDs).For example, exclusive coherent diffractive dijet production in electron-proton collisions has been suggested as a probe of gluon GTMDs in [1] (see also [2]).By measuring the transverse momenta of the two jets one can access the off-forwardness as well as the transverse momentum distribution of gluons, albeit the latter only in an indirect way through a weighted integral where the weight is a function of the observed momenta (see e.g.[3]).Exclusive coherent diffractive J/ψ production in electron-proton collisions can also be used [4,5], but in this case the weight of the integral cannot be varied substantially (double J/ψ would be more suited for that), in that sense it is closer to processes like DVCS that probe GPDs, which are fixed integrals of GTMDs [6].Nevertheless, if the underlying GTMD description is valid, then the various processes should be describable simultaneously by the same GTMDs.There is data from the H1 and ZEUS experiments of HERA, which can be used to check this to a certain extent.In this paper we will attempt such a combined analysis and reach the conclusion that there is considerable tension between the optimal dijet and J/ψ descriptions.This tension can be due to the theoretical assumptions that go into the GTMD description, such as the selected GTMD model form or the model for the J/ψ wave function, but can also have an experimental origin.This offers an excellent opportunity for the future U.S.-based Electron-Ion Collider (EIC) as it can provide additional and more precise data on both processes in different kinematical regions.Other tests of the underlying GTMD description can come from Ultra-Peripheral Collisions (UPCs) data of RHIC and/or LHC.Here it is important to ensure that the processes are exclusive and coherent diffractive which means the proton or nucleus has to remain intact.
In [3] we have considered a McLerran-Venugopalan (MV) based (small x) model in order to obtain a description of diffractive dijet production data of H1.Although the data is neither fully exclusive, nor fully coherent, the kinematics is such that in leading order a GTMD description is expected to be appropriate.The model was indeed able to describe the t dependence and the transverse momentum of the jets, but since it did not include any x dependence, the Q 2 and y dependence was not described very well.For the present study the first step will thus be to include an x dependence, to obtain a better description of the Q 2 and y dependence, before moving on to the description of J/ψ production.We will show that the x-dependent GTMD model can provide a good description of the t dependence of the exclusive coherent diffractive J/ψ production data of H1 and ZEUS for different Q 2 values, including the photoproduction case.However, it turns out there is no choice of parameters that leads a good description of the dijet and J/ψ data simultaneously.The tension is solely in the slope of the exponential fall-off in t, which in the model is entirely governed by the width of the proton profile.We do no intend to resolve this tension in this paper, because there may be other reasons for the tension beyond the GTMD model and future data will be needed to confirm or refute the tension.We use the optimal model for the J/ψ data to obtain predictions for the EIC and for UPCs, rather than the model that has the minimal tension with the dijet data as it would not yield a satisfactory description of either process.
In Sec.II we collect the key ingredients of the diffractive dijet production cross section description in terms of the small-x gluon GTMD of our previous study [3].We then propose an x-dependent GTMD model to improve the description of the HERA-H1 data in Sec.III and show the improved fit in Sec.IV.Our aim is to use as much as possible those fit parameters to also describe diffractive J/ψ production data at HERA and LHC.The diffractive J/ψ production cross section expression in terms of GTMDs is given in Sec.V, followed by discussions about certain phenomenological corrections often considered.The fit results to H1 and ZEUS γ ( * ) p data are presented in Sec.VI where predictions to the future EIC are also provided.By further fitting ALICE Run 2 data on mid-rapidity UPCs off Pb nuclei (in order to include A dependence in the model), we then give predictions for RHIC, LHC (at Run 1 center of mass energy), and EIC (for Au nuclei) in Sec.VII.We summarize our findings in Sec.VIII.

II. DIFFRACTIVE DIJET PRODUCTION IN TERMS OF GTMDS
We first recall the essential expressions for the dijet case in terms of the small-x gluon GTMD which is the angular independent part of the more general GTMD where θ k∆ denotes the angle between k ⊥ and ∆ ⊥ .In this paper we will only consider the angular independent part F [ ] 0 , under the assumption that the contributions of the angular modulations (that enter the cross sections squared) are much smaller.For a detailed analysis of angular modulations (and shape fluctuations) we refer to [7].We also consider the off-forwardness to be entirely in the transverse direction, i.e. we consider zero skewness ξ.
In the strict x → 0 limit and for zero skewness (ξ = 0) one finds the following expression in terms of the dipole scattering amplitude N = 1 − S [ ] C : The average with subscript C is over the color configurations of the proton or nucleus considered.The Wilson loop operator and U [±] are the standard staple-like gauge links in the forward (+) and backward (−) lightcone directions (see e.g.[8,9]).Furthermore, r ⊥ = y ⊥ − x ⊥ and b ⊥ = (x ⊥ + y ⊥ )/2.The MV-like model used for S [ ] will be discussed in the next section.The γ * p differential cross section for Q 2 > 0 can be expressed in terms of amplitudes A T and A L for transverse and longitudinal photon polarizations, respectively: and where j denotes an outgoing jet, y is the inelasticity, e f is the electric charge of a quark of flavor f in units of the positron charge, z = k + 1 /k + is the outgoing quark (jet) longitudinal momentum fraction with respect to the virtual photon longitudinal momentum, as follows: and These expressions show that exclusive coherent diffractive dijet production allows to obtain information on GTMDs even though the transverse momentum dependence is integrated over.The dependence on the external momenta of the weights inside the integrals can be used to study the transverse momentum dependence of the GTMDs.For photoproduction γp one can set Q 2 = 0 and drop the longitudinal part as it does not give any contribution.The expressions also show that exclusive coherent diffractive dijet production allows to obtain information on GTMDs that goes beyond the GPDs.In [6] it was pointed out that the gluon GPD H g at small x can be expressed in terms the small-x gluon GTMD through the following integral 1 : This relation is derived using the operator definitions of the functions, not taking into account renormalization.However, this relation suffers from the same problems as its forward limit, where the collinear PDF f (x) is viewed as the integral of a Transverse Momentum Dependent parton distribution (TMD): f (x) ⊥ ) [10].Since the tail of the TMD f (x, q 2 ⊥ ) behaves as 1/q 2 ⊥ , the integral will diverge logarithmically and requires some form of regularization.More formally, beyond tree level TMD factorization implies that the TMD depends on two scales, the rapidity scale ζ and the renormalization scale µ, satisfying two coupled evolution equations, whereas the collinear PDF only depends on the scale µ and satisfies just one evolution equation.Clearly, Eq. ( 7) suffers from the same problems.Because of this, here we do not provide curves for GPDs based on the GTMDs we obtain, as they would unavoidably depend on the adopted procedure of how to regulate the large transverse momentum behavior to make the integral in Eq. ( 7) converge.Rather than expecting that the TMD determines the collinear PDF through an integral relation, one can instead consider the unambiguous relation in which the collinear PDF determines the large transverse momentum dependence of the TMD, i.e. f (x, q 2 ⊥ ) ∝ α s q −2 ⊥ dy y P ( x y )f (y) , where P denotes a splitting function (ignoring for simplicity the possibility of mixing among various pdfs).Similar expressions hold for the perturbative large transverse momentum tails of GTMDs in terms of GPDs, as recently studied at the one loop level in [11].So rather than using fits of GTMDs to obtain results for GPDs, it is better to use models, lattice determinations, or fits of GPDs to predict the tails of the GTMDs and compare those to GTMD fits.This we do not attempt here, as we are primarily concerned with the small transverse momentum region in the present paper.

III. AN x-DEPENDENT GTMD MODEL
In our earlier study of diffractive dijet production [3] we used a small-x model without any x dependence, because the expressions in terms of the Wilson loop operator were obtained in the strict x → 0 limit [9,[12][13][14].The model was able to describe the H1 data on the differential cross sections as function of t and K ⊥ quite well, for which all data have the same small average x value.But the model did not describe well the data as a function of Q 2 for which the average x value is different for each data point or each bin.
To improve on this we now incorporate an x dependence in the model by replacing the constant free model parameter χ by the following function of x and Q 2 : with x 0 = 3 × 10 −4 and λ = 0.29 based on the model by Golec-Biernat and Wüsthoff (GBW) for the saturation scale [15,16].This value of λ also turns out to allow for a reasonably good description of the Q 2 dependence of the H1 diffractive dijet electroproduction data (using x = s/(yQ 2 )) for diffractive dijet production, but as we will see later J/ψ production prefers a somewhat smaller value λ = 0.22.Our motivation to include this type of x dependence in our model is the observation of geometric scaling behaviour of DIS ep collisions data in the low x and low Q 2 region at HERA [15,16].This feature of the data was well-described by a saturation scale of the form Q 2 s (x) ∼ x −λ with λ ≈ 0.29.Later, it was shown that the total cross-section of γ * p exhibited geometric scaling over a much wider range of Q 2 values (from 0.045 to 450 GeV 2 ) in the x < 0.01 region of HERA data [17] (see also [19,20]).Similar scaling was also observed in diffractive DIS with a specific parameterization [18] and in inclusive eA processes [21].
We note that apart from this new Q 2 dependence introduced through the kinematic relation between x and Q 2 , we do not introduce any additional x and/or Q 2 dependence from QCD corrections.We also consider a fixed coupling constant.The reason for not including evolution is that the scale evolution of GTMDs has not been studied in full 1 Taking into account that the definition of differs from that of F 0 in [6], Eq. ( 7) differs by a factor 2 from the one given in [6].
yet.Even for TMD evolution the interplay between the x and Q2 evolution is not yet fully clear (the scale evolution of TMDs and Sudakov resummation for TMD processes at small x has been investigated in e.g.[22][23][24][25]).Since the kinematic range of EIC is not too different from the one of H1 and ZEUS, we expect evolution not to be essential for obtaining predictions.Given the large uncertainties in the model parameters we expect logarithmic corrections to be of minor importance at this stage.This also avoids the question of what precisely sets the hard scale in the various processes (dijet versus J/ψ, electroproduction versus photoproduction).But the kinematic relation between x and Q 2 for the J/ψ case is affected by the J/ψ mass, of course.
With all these considerations taken into account we arrive at the following McLerran-Venugopalan (MV) [26][27][28] like model for the dipole scattering amplitude: Here the factor e − r r 2 ⊥ is introduced to ensure that the dipole sizes probed are restricted to the perturbative region.This Gaussian weight factor was also introduced in Wigner, Husimi and GTMD distributions [29,30].In practice this should be enforced by the kinematics of the process, i.e. by the large transverse momentum of the jets or the large mass of the produced quarkonium, but by enforcing it explicitly in the model, we can obtain convergent integrals of the GTMD without reference to a process.In principle, r is a free parameter of the model but we will use the fixed value r = (0.4 fm) −2 for both diffractive dijet and J/ψ production, which corresponds to the gluonic radius of the proton used in [31].The free parameter χ is chosen such as to obtain a reasonable description of the data.
The above corresponds to a leading order description of a dipole interacting with a proton or nucleus through two-gluon exchange in the t-channel.The resulting S [ ] C , N , and are purely real.In principle GTMDs can be complex, with an imaginary part that is referred to as the odderon contribution.In the forward limit for unpolarized protons this contribution has to vanish for the dipole case, but odderon contributions may arise for nuclei or from quadrupoles or higher multipoles.Even when the odderon contribution is considered absent down to a certain small x value, nonlinear QCD evolution would generate a nonzero contribution for even smaller x values [32].Therefore, in principle the imaginary part has to be included, but that has not been done in diffractive dijet production thus far.Later on we comment on the size of the expected correction from the imaginary part.
The saturation scale in the model will be taken as with T (b ⊥ ) is the profile of the proton or nucleus.For the proton we use a Gaussian profile with R p the gluonic radius of the proton [31].For heavy nuclei the profile is described by the thickness function [33] T which is obtained from the Woods-Saxon distribution Here, N A is a normalization factor such that T A (b ⊥ = 0) = 1.For the numerical calculations we choose the nuclear radius R A = 1.12A 1/3 fm with A the mass number of the nucleus 2 .For lead (A = 208) we use a A = 0.546 fm and for gold (A = 197) we use a A = 0.535 fm.For the proton we use Q 2 0s,p = 1 GeV 2 , while (similar to what was done in [31]) for the heavy nuclei Q 2 0s,A can be obtained from the relation to be evaluated at the same value of x.Here η will be considered a free parameter to be fitted to the data.Explicitly, for heavy nuclei we have we obtain the following expression for the saturation scale which will determine the saturation scale of the heavy nuclei.The general expectation is that η should have a value near 1.

IV. DIFFRACTIVE DIJET ELECTROPRODUCTION
The diffractive dijet electroproduction cross section ep → e jjp can be related to the γ * T,L p → jjp cross sections in Eqs. ( 3) and ( 4) as follows: In the calculation, we use Λ = 0.24 GeV, C F = 4/3, N c = 3, α em = 1/137, N f = 4, and fixed α s = 0.3.With the above model and cross section expressions we fit the H1 data of [35].Although the data is for the production of at least two jets and not fully exclusive, a leading order description in the measured kinematic regime will be dominated by the production of two jets in any case.Also the data is not fully for coherent diffraction, but the t values are so small and the kinematic cuts on the final state proton are such that the additional light hadrons that may be produced in the process are not expected to alter the process much compared to the coherent case.The data is consistent within errors with actual exclusive coherent diffractive dijet production data from ZEUS [36] (to be specific, this we checked for dσ/dβ at β = x B /x I P = 0.1).The ZEUS data is less differential however and therefore not used here.Finally, the H1 data do not have zero skewness, in fact, it may typically be larger than the x values (the average value for x I P ≈ ξ is around 0.03-0.04[35], whereas the geometric average of the upper and lower value of the x range is 10 −3 ), meaning that the ERBL region is probed rather than the DGLAP one, cf.e.g.[37].These caveats concerning this data should be kept in mind, when we discuss the tension with the best GTMD model description of the J/ψ production case (for which ξ is also not exactly zero, of course).This is also one of the reasons why we do not attempt to find the best model fit that can describe both processes simultaneously.The purpose here is to demonstrate those features that the model fits have in common and those that are in tension, such that future experimental investigations can focus particular attention on these aspects.
The best fit to the t distribution of H1 data [35] is shown in Fig. 1.The t distribution is well described by an exponential fall-off, dσ/dt ∝ exp(−bt), as also is the case for the t distribution for J/ψ production.The slope b is mainly determined by the gluonic radius of the proton R p and we find that the model gives the best fit for R p = 0.49 ± 0.02 fm which leads to the slope b = 6.0 ± 0.5 GeV −2 (in [35] b = 5.89 ± 0.50 GeV −2 is given for the H1 data).Incorporating the statistical and systematic uncertainties of the data which are simply added in quadrature, we give a band around the best fit (central value) χ = 1.5 ± 0.1.Within errors, our model gives a reasonable description of the Q 2 , K ⊥ and y distribution data, an improvement w.r.t. the x independent model of our previous study [3].

V. EXCLUSIVE COHERENT DIFFRACTIVE J/ψ PRODUCTION
The differential cross section of exclusive coherent diffractive vector meson production as a function of momentum transfer squared t = −∆ 2 ⊥ (for ξ = 0) can be written as [4,38] expressed in terms of the elastic scattering amplitudes for the transverse/longitudinal polarization (T /L) of the virtual photon.Here b ⊥ denotes the impact parameter, S(r ⊥ , b ⊥ ) is (the real part of) the dipole-proton scattering amplitude, and (Ψ * V Ψ) T,L is the overlap between the photon and the vector meson wave functions which depends on the transverse dipole size r ⊥ , momentum fraction z carried by the quark, the quark mass m f , the vector meson mass M V , and the photon virtuality Q 2 .We are working in a frame where the colliding virtual photon and proton have zero transverse momentum.We note that in Eq. ( 19) we use the phase factor e i( 1 2 −z))r ⊥ •∆ ⊥ as suggested in [6], rather than e i(1−z))r ⊥ •∆ ⊥ as used in [4].In practice, it does not make much difference though.In models studied by [39], the cross section amplitude using the phase factor from [6] is only 3.5% larger than using the one from [4].
As in the dijet case, the model parameter R p of the proton profile directly determines the t slope of the differential cross section for J/ψ production.However, the latter is found to be quite different from that of the dijet case, with the J/ψ production slope substantially smaller.They are only compatible at the 3σ level, hence there is considerable tension.Thus far no phenomenological study for both types of processes simultaneously has been performed and we will not attempt that here, but in principle both should be describable within the same GTMD based approach.Due to this aspect of tension, the uncertainty in the model parameters is larger than what is suggested by studying only one of the two types of process.Future data is needed to resolve this issue.
The photon-vector meson wave functions overlaps are expressed as [4] ( where φ T,L are model dependent scalar functions and δ is usually taken to be 0 or 1.Because the wave function overlap of the longitudinal part is linearly dependent on Q, the chosen value of δ does not affect photoproduction (Q 2 ≈ 0).However, the chosen δ value significantly affects the large Q 2 cross section as the longitudinal part becomes a larger portion as Q 2 grows: the second term of (Ψ * V Ψ) L grows faster than the first term when δ = 1.In this study, we will use the so-called Gaus-LC (GLC) and boosted Gaussian (BG) vector meson models [4,44].The GLC model is given by: For V = J/ψ the parameters are listed in [4]: GeV.The BG model is given by [40,41]: with 2 .The parameters are listed in [4] (see also [42]): N T = 0.578, N L = 0.575, R 2 = 2.3 GeV −2 , and M V = 3.097 GeV.For both models the parameters are fixed by requiring that the wave function is normalized and that the decay widths f V = f V,T = f V,L are reproduced.We note that other vector meson wave functions and potential models have been considered in order to describe diffractive production of J/ψ and other heavy quarkonia, see for example [43][44][45][46][47][48][49][50].These descriptions can in principle also be translated into GTMD model expressions, which possibly allows for a more direct comparison of models (rather than a comparison of how they describe the data), but that is not our objective here.
Before moving on to the description of the data, we make some observations about whether this process probes a GPD rather than a GTMD.If we restrict to the angular independent part only, then we can express Eq. ( 19) as with δ ⊥ = ( 1 2 − z)∆ ⊥ .This expression shows that also in diffractive J/ψ production one probes an integral over a GTMD with an integrand that depends on the kinematic variables of the process (in this case z and ∆ ⊥ ) and hence different integrals of the GTMD can be obtained in this way, even though with less possibilities for varying the integrand than in the dijet case.This also means that the expression cannot be given in terms of a GPD (which does not depend on z), only in an approximation, as pointed out by [6].As the weight of the integral over the GTMD depends on ∆ ⊥ through δ ⊥ and ∆ ⊥ is generally small and only relevant in a small kinematic region (and the region around z = 1/2 contributes the most), this dependence may be ignored to good approximation leaving a fixed integral over the GTMD, which however still is not an expression in terms of the gluon GPD H g through Eq. ( 7).This would require small |q ⊥ + δ ⊥ |, for which we can consider the first order expansion of the Bessel function , and use the fact that d 2 q ⊥ F [ ] 0 (x, q ⊥ , ∆ ⊥ ) = 0 [6], to find that where in the last step we used Eq. ( 7).It should be stressed that this is an approximation that depends on the relevant range of the q ⊥ integration.Therefore, we will consider the GTMD expression rather than the approximate one in terms of the GPD.Diffractive scattering in terms of GPDs has been studied in [51].

VI. ANALYSIS OF COHERENT DIFFRACTIVE J/ψ PRODUCTION DATA
The total γ * p cross section for J/ψ production studied by H1 [52] is defined as σ tot = σ T + εσ L with ε = (1 − y)/(1 − y + 1 2 y 2 ) (with ε = 0.99 in [4,53]), while by ZEUS [54] it is defined as σ tot = σ T + σ L .Here we will also use ε = 0.99, which corresponds to y 0.13.This differs from the dijet case where a range of y values is considered and the x value depends on y.The x value used in the J/ψ production case is determined by [38,55] We recall that σ L = 0 for the photoproduction case.Diffractive vector meson production at HERA and LHC have been studied extensively with various model approaches, such as [38,39,42,[56][57][58][59][60][61][62].It has been shown in [4,5,63,64] that MV-like models can describe diffractive vector meson production well.In most model studies the dipole scattering amplitude is real and hence A γ * p→V p T,L imaginary, but as mentioned, in reality there will be an imaginary part.The phenomenological correction 1 + β 2 is introduced to account for that contribution.Using dispersion relations an expression for β has been obtained [55,65,66]: This expression implies that only x dependent elastic scattering amplitudes yield nonzero β.Phenomenological studies of HERA data [55,67,68] find that the real part correction (1 + β 2 ) is anywhere between 10% and 25%.Given the large uncertainty in the model fits that we obtain, this correction will not be of much importance.Therefore, we will first fit the data without this correction and then estimate the size of the correction for that fit afterwards.In this way we find that in our case (1 + β 2 ) is in the 10-15% range.More details on this will be presented below.
Another correction often considered comes from taking into account nonzero skewness.The off-forwardness ∆ = P − P for zero skewness, i.e. ξ = −∆ + /(P + + P + ) = 0, means ∆ = ∆ ⊥ and t = −|∆ 2 ⊥ |.In practice ξ will not be exactly zero.Therefore, in order to correct for this the factor R 2 g is included, where R g for gluons at small x and small ξ is given by [66] where H g (x, ξ) is the standard (helicity non-flip) gluon GPD (now for nonzero skewness, unlike in Eq. ( 7)).This factor R 2 g is a substantial correction for HERA kinematics, found to be in the order of 40-70% [55,67,68].We note though that the value H g (x = ξ, ξ) is at the boundary of the ERBL and DGLAP regions, where the function is continuous but not differentiable, hence changing abruptly, so a slight change in x w.r.t.ξ or vice versa can matter considerably.This introduces an uncertainty regarding the actual correction that is needed, as the data span a range of x and ξ values.Furthermore, the specific value x = 2ξ, which is in the DGLAP region, stems from the GPD analysis of [69] for which the region around this value is found to make the dominant contribution.This is quite different from our GTMD approach for which x = 2ξ plays no dominant role and the ratio H g (x, ξ)/H g (x, 0) (i.e. with the same x value in numerator and denominator) for small nonzero ξ seems more appropriate to consider.Moreover, as mentioned, the data are more likely in the ERBL region (x < ξ), although that is not clearly specified in the experimental papers.Hence, given these considerations, here we do not include the R 2 g correction factor in this paper and also not correct for nonzero skewness in another way.Given the large uncertainties in the model fits (due to the large uncertainties in the data), this is not expected to be essential.
In Fig. 2 we show that our MV-like model with free fit parameters χ, R p and λ and without the mentioned corrections can achieve a good description of the t distribution data of HERA [52,54] and of the data on the total cross section as a function of W from HERA (H1 [70] and ZEUS [71]) and LHC (ALICE [72] and LHCb [73,74]).The t dependence of the model is controlled by the proton profile R p , the W slope by λ, while the amplitude of the cross section is χ dependent where large χ means a larger cross section.Here we give preference to the description of the photoproduction data, i.e. the parameters are obtained from a simultaneous description of the t and W dependence of the photoproduction data, where there is only a very narrow range of R p and χ values that describe those data simultaneously.If one would include electroproduction data, the parameters would change considerably (cf.Fig. 3 (left)) such that the photoproduction data would be described less well.Here we note that the H1 and ZEUS data sometimes differ quite a bit from each other, when comparing data sets at the same or very similar values of Q 2 .Therefore, we prefer to let the Q 2 dependence be determined by the model after the parameters are fixed at Q 2 = 0.05 GeV 2 .It then turns out that the t dependence is well described using GLC for all values of Q 2 considered, while BG overestimates the data for large Q 2 .
We show bands of values of χ and R p which describe the photoproduction data qualitatively equally well within the errors of the data points.However, they should not be interpreted as 1σ error bands, as they are not obtained from a fit to all data points through a minimization w.r.t.all parameters simultaneously.Rather, we determine R p from the t dependence, λ from the W dependence, and subsequently obtain χ.Given the uncertainties in the data and the tension with the dijet data, a more sophisticated determination of the parameters and the error bands does not seem to be called for at this point.To be specific about the bands, for GLC we use χ = 1.45 − 1.50 for R p = 0.40 fm and χ = 1.40 − 1.45 for R p = 0.41 fm, while for BG we use χ = 1.10 − 1.15 for R p = 0.40 fm and χ = 1.05 − 1.10 for R p = 0.41 fm.We choose two different possible R p because of the fact that the photoproduction data prefers a steeper slope (R p = 0.41 fm) than the electroproduction data (R p = 0.40 fm).The parameter λ is determined by the slope of the W dependence data of the J/ψ photoproduction total cross section from HERA and LHC in Fig. 2 (right), to which it is very sensitive.It turns out that unlike the dijet case this data prefers λ = 0.22 < λ GBW .We note that trying to obtain a better fit of the W dependence of the total cross section will lead to a less good description of the photoproduction differential cross section dσ/dt.As the model is probably less appropriate for the integral over all t, we have given preference to the latter.We refer to [34] for a combined description of H1 data of total photoproduction cross section and the differential cross section for both coherent and incoherent diffraction, taking into account proton shape fluctuations.
FIG. 2: Left: GTMD model fit of the t dependence of exclusive coherent diffractive J/ψ production data from H1 [52] and ZEUS [54], for two vector meson wave function models: GLC and BG.The systematic and statistical uncertainties are added in quadrature.The bands correspond to the ranges of values of χ and Rp specified in the text.Right: Fit of the model to the W dependence of the total J/ψ photoproduction cross section at HERA and LHC [70][71][72][73][74].The data are only shown for x 0.01 which corresponds to W > 30 GeV (≈ 10M J/ψ ).
In Fig. 3 (left) the χ 2 per degree of freedom (dof) as a function of χ is shown for the case that a fit is made to both the photo-and electroproduction data.It can be seen that also in this case GLC gives a slightly better minimal χ 2 /dof and prefers a χ close to the one of the dijet case.Therefore, the GLC model seems to be preferred.However, when it comes to the t and W slope, all vector meson wave function models require χ, R p and λ values that are smaller than those obtained from the dijet data.This is clearly visible in Fig. 3 (right) where we show the value of b of d 2 σ/dt dQ 2 ∝ e −bt for ep collisions, where b is solely determined by R p in Eq. (11).The preferred proton profile for the J/ψ case has R p = 0.40 − 0.41 fm, while for the dijet case it is R p = 0.49 fm.Here it should be recalled that the J/ψ data is for fixed y (γ ( * ) p) and the dijet data is y integrated (ep) (cf. the bottom right plot of Fig. 1 for the y dependence of the dijet data).The d 2 σ/dtdQ 2 slope data points in Fig. 3 (right) are only available for J/ψ production [54,71,75], while the Q 2 dependence of the dijet slope is extracted from the fit (see Fig. 1).The dijet slope data is given only for Q 2 and y integrated which is b = 5.89 ± 0.50 GeV −2 .The bands of the t slope (b) for the J/ψ case shown in Fig. 3 (right) reflect the aforementioned R p values, while for the dijet case the bands correspond to the 1σ error in R p , i.e.R p = 0.49 ± 0.02 fm, for which a larger R p gives a larger b.Our model shows that the dijet slope is slowly increasing in Q 2 while for J/ψ it is steadily decreasing.This tension cannot be resolved with the present set of just three free parameters.It is also clear that it cannot be attributed to the vector meson wave function, but stems from the proton profile.Of course, it may be (in part) due to the aforementioned caveats about the dijet data, but without additional future data, that can likely not be clarified.
FIG. 3: Left: χ 2 /dof vs χ for dijet production (combined Q 2 , t, K ⊥ , and y dependence data [35] with dof=15) and J/ψ production (combined photo-and electroproduction data [52,54] with dof=35) for two different possible values of Rp.Right: t slope (b) of the model as a function of Q 2 for the dijet and J/ψ cases.The J/ψ data are taken from [54,71,75].Note that the diffractive dijet production case is y integrated, while the J/ψ case is evaluated at W = 90 GeV and for fixed y.
In Fig. 4 we show the W distribution resulting from the model fits to the t-dependence.Both models can describe well the small Q 2 data (photoproduction).For electroproduction GLC overestimates the data by at most 2σ, while BG shows a larger deviation from the data, as shown in Fig. 4 (right).With the obtained model fits we provide predictions for the same process at EIC, which will cover a different kinematic region, but not so different that evolution will play a big role.The left panel of Fig. 5 shows the predictions without the phenomenological correction 1 + β 2 and the right panel with.In Fig. 6 we show the correction by itself.
In the left panel of Fig. 6 the correction obtained with the model is plotted for γ ( * ) p collisions as a function of t and found to be in the 10-15% range, a bit larger in the BG model than in the GLC model, and slowly increasing in Q 2 .This is in agreement with other phenomenological studies [55,67,68].For longitudinal photon polarization a similar size correction is obtained.The right panel shows similar size corrections for UPCs to which we will turn next.Ultra-peripheral Heavy Ion Collisions at LHC and RHIC can be used to study photon-nucleus collisions.In order to make sure that one is dealing with exclusive coherent diffractive production of a J/ψ there need to be rapidity gaps between the J/ψ and both nuclei.We consider the J/ψ to be produced at mid-rapidity.The central rapidity range of the LHC is |Y | < 0.8 and at RHIC |Y | < 1.Here we simply take Y = 0.In that case the colliding photon and gluon will both have an energy of M V /2, which means that the gluon has a momentum fraction of x g = M V / √ s N N .For the LHC Pb-Pb UPC data [76] at √ s N N = 5.02 TeV (Run 2) this corresponds to x g = 6 • 10 −4 and at √ s N N = 2.76 TeV (Run 1) to x g = 0.001.For the RHIC Au-Au UPC data at √ s N N = 200 GeV this corresponds to x g = 0.015, which is at the edge of the range of applicability of the MV-like model that we are using.We will use the ALICE data to fit η and then obtain predictions for RHIC, keeping in mind this caveat.
In order to compare the UPCs to the photoproduction case at EIC, it is useful to know what is W γN in the UPCs.For A-A UPCs at mid-rapidity (Y = 0), the photon-nucleon center of mass energy squared is determined by In Fig. 7 we provide a fit of our model to the t dependence of the differential cross section d 2 σ/dY dt [76] for coherent diffractive J/ψ production at Y = 0 in ultra-peripheral Pb-Pb collisions at √ s N N = 5.02 TeV.According to [76], various models can describe the ALICE data well.One such model is the leading-twist approximation (LTA) of nuclear shadowing [79,80], which combines the Glauber-Gribov formalism with the phenomenology of photon diffraction from HERA.The lower bound of the GLC fit in our model is close to the LTA (low shadowing) model result with a slightly steeper slope, as shown in Fig. 7 (right).Another model, the b-BK [39,81,82], was proposed based on the solution of the Balitsky-Kovchegov equation [83,84] with an impact parameter dependence.Another study that incorporates nucleon shape fluctuations [34] also provides a good fit to the data, including the coherent W dependence of the photoproduction total cross section and the t distribution of coherent and incoherent J/ψ photoproduction data from HERA.While the former two models utilize the BG wave function, our model provides a better description of the t and W dependence data using the GLC wave function.Another wave function model based on the Buchmüller-Tye potential [50] that uses r-b correlation [85] with two different parameterizations: GBW [15,16] and BGBK [86], also gives good agreement with the data [87].Differences in the magnitude and slope of the cross section between other models and ours could also be due to the use of different nuclear radius parameters.
Extrapolating the fit gives a prediction of the first diffractive minimum (or dip) to be at t 0.016 GeV 2 .We find that the dip position is determined by the target profile R A , such that it will move towards a smaller t value for larger R A .The fit turns out to be very sensitive not only to the value of R A , but also to the power of A. With the fit of the model to the ALICE data, we find that η = 0.96 ± 0.01 for GLC and η = 0.95 ± 0.01 for BG give the best fit.Therefore, our model fit indicates that the saturation scale behaves like Q 2 s ∝ A 0.27−0.30, not A 1/3 .The latter in fact does not provide a good fit.Following Eq. ( 16), the saturation scale Q 2 0s,A for the heavy nucleus A depends on η, R A , and R p , and is in all cases found to be between 1.5 and 1.9 GeV 2 .The choice of proton and nuclear profile, particularly the Gaussian profile used in our study, can affect the fitted η value and may differ with different profiles.We did not attempt to find profiles that would lead to η = 1 and do not exclude that that is possible.With the same parameterization, in Fig. 8 we provide predictions for RHIC and LHC (Run 1) at midrapidity.In general, both wave function models give very similar results but GLC shows a somewhat larger band than BG as expected from the previous analysis on γp (see the χ 2 /dof in Fig. 3).The left panel in Fig. 8 is for Au-Au UPCs at √ s N N = 200 GeV.Fig. 8 shows that the first diffractive dip for Au-Au at RHIC is predicted to be around t 0.017 ± 0.001 GeV 2 which is close to its location in RHIC preliminary data [77] and other studies [57,78], while for Pb-Pb at LHC Run 1 it is predicted to be around t 0.015 ± 0.001 GeV 2 .The middle panel in Fig. 8 is for Pb-Pb UPCs at LHC Run 1 with √ s N N = 2.76 TeV.In the rightmost panel of Fig. 8, we provide predictions of e-Au collisions at the EIC for photoproduction (Q 2 = 0 GeV 2 ) and electroproduction (Q 2 = 10 GeV 2 ) for fixed x g = 0.01 which corresponds to W γN = 31 GeV and W γN = 44 GeV, respectively.The bands for each wave function model reflect the uncertainties on R p (which translates to η) and χ.As shown in Fig. 6, the β corrections are in the order of 6-10% for UPCs, which is small compared to the uncertainties, hence not included here.

VIII. DISCUSSION AND CONCLUSIONS
We have shown that incorporating x dependence in our previous gluon GTMD model, along the lines of the GBW parameterization of the saturation scale, will improve the description of HERA-H1 diffractive dijet production data.However, describing diffractive J/ψ production in the same way leads to tension for the slope of the t dependence, which is quite distinct for the two processes, where dijet needs a steeper slope than J/ψ production.In the model this slope is solely determined by the (Gaussian) proton profile and there appears to be no clear way to resolve the tension.A few other differences between the optimal parameter choices for dijet and J/ψ production are found, but these can be reduced by adjusting the J/ψ wave function or by modifying the x dependence w.r.t. the GBW parameterization.For instance, dijet production can be described well by an x dependence of the saturation scale Q s ∝ x −λ with λ = λ GBW = 0.29, while the W dependence of photoproduction of J/ψ's prefers a smaller λ ≈ 0.22.A non-constant λ may be needed, but we did not explore that option, anticipating that future more precise data may shed new light on the differences between dijet and J/ψ production.We do expect that a common gluon GTMD model description of dijet and J/ψ production may be possible once the slope issue is clarified by new data.
With the best fit of our model to combined H1 and ZEUS data on J/ψ production, we have provided predictions for the diffractive J/ψ production in e-p collisions at the future EIC.We further fit our model to UPC data from ALICE (Run 2) to determine the A dependence of the saturation scale.The fit turns out to be very sensitive to the power of A and the nuclear profile R A , which suggests that also the profile function shape will matter considerably.We find an A dependence that is slightly less than the generally expected A 1/3 , to be specific, Q 2 s ∝ A 0.27−0.30 .With the obtained fit we provide predictions for UPCs at LHC (Run 1) and RHIC, and for e-Au collisions at EIC.We have also investigated a phenomenological correction commonly used for J/ψ production which comes from the odderon contribution, which is in the 10-15% range and thus unimportant given the present large uncertainties in the available data and hence in the model.The larger correction from accounting for non-zero skewness that is often considered in GPD approaches in the DGLAP regime does not seem appropriate for our GTMD model and was thus not included.Moreover, a similar correction has not been applied in dijet production, which would affect the comparison.We have also pointed out that dijet and J/ψ production go beyond probing GPDs, rather they probe weighted integrals of GTMDs with weights that depend on external kinematical variables of the process that can be varied and exploited.
The predictions from the presented x-dependent gluon GTMD model for J/ψ production at EIC, LHC, and RHIC, will hopefully facilitate resolution of the t distribution tension with dijet production, clarify the dependence on the skewness probed in the process, and determine the x and A dependence of the saturation scale of heavy nuclei.

and Q 2
is the photon virtuality.The amplitudes A T,L are given in terms of F [ ] 0

FIG. 1 :
FIG. 1: The best fit of diffractive dijet production to H1 data with Rp = 0.49 fm, λ = 0.29, and r = (0.4 fm) −2 .The central value corresponds to χ = 1.5 while the band correspond to χ = 1.4 and χ = 1.6 for the lower and upper band respectively.The systematic and statistical uncertainties are added in quadrature.

FIG. 4 :
FIG. 4: W dependence of the model fits compared to the H1 data for photoproduction (left panel) and for Q 2 = 8.9 GeV 2 (right panel) for the two different wave function models and the same model parameters as in Fig. 2.
78].For LHC Pb-Pb UPC data at √ s N N = 5.02 TeV this corresponds to W γN = 125 GeV and at √ s N N = 2.76 TeV to W γN = 93 GeV, while for RHIC Au-Au UPC data at √ s N N = 200 GeV this corresponds to W γN = 25 GeV.

FIG. 7 :
FIG. 7: (left) Fit of the model to ALICE (Run 2) data [76] of coherent diffractive J/ψ production at midrapidity in ultraperipheral Pb-Pb collisions at √ sNN = 5.02 TeV.We use the same parameter values as in Fig. 2, but fit the additional parameter η that determines the saturation scale for nuclei.(right) Ratio of the model fit to the ALICE data for lower and upper bounds of each wave function.