Energy dependence of the proton geometry in exclusive vector meson production

The gluon radius of the proton is expected to increase at small gluon momentum fractions $x$, an effect which has hitherto not been considered in the dipole model framework. We investigate the energy dependence of exclusive $J/\psi$, $\phi$, and $\rho$ production by introducing three models for $x$ dependence of the gluon thickness function. We allow the transverse width of the proton to increase as $x$ decreases, using novel parametrisations in the spherical proton and the hotspot model. We compare these with a model where the number of hotspots increases as $x$ decreases and confront the models with HERA data. The models exhibit clear differences in the slope of the $t$-spectra and in the cross section ratios between coherent and incoherent events. Comparisons to $t$-slopes and $W_{\gamma p}$ measurements show a preference for models where the proton's size increases as $x$ decreases.


I. Introduction
Diffractive vector meson production is characterized by a momentum exchange with vacuum quantum numbers in the Mandelstam t-channel between the virtual photon and the target in ep and eA collisions [1][2][3]. In the target rest frame, this is described by the dipole picture [4][5][6][7][8][9][10] where the virtual photon splits up into a quark anti-quark dipole which interacts with the target via the strong interaction and then forms the final state. In a coherent event, the target stays intact, while in a dissociative or incoherent event, the target subsequently breaks into fragments. Coherent and incoherent diffractive events have been extensively studied at the HERA ep experiments H1 and ZEUS [11][12][13]. The coherent events are sensitive to the transverse structure of the target, while the measurements of incoherent events provide information about fluctuations in the target wave function.
The measurements of inclusive deep inelastic scattering (DIS) at HERA and its theoretical description with the collinear factorisation framework with DGLAP evolution equations led to an understanding that gluons exist in the form of quantum fluctuations inside the proton, and their density increases steeply at higher energies. Though this framework has been very successful in dealing with hard processes, it could not explain diffractive events. The gluon density is also expected to saturate at small gluon momentum fractions x so that unitarity of the events is upheld, even if a clear signal for saturation is yet to be seen in experiments. The dipole framework provides a unified description of inclusive and diffractive events, as well as exclusive diffraction, at small x. It also is a natural model for describing saturation physics.
In the dipole picture, the impact parameter is a Fourier conjugate to the momentum transfer at the target vertex, thus one can study the transverse structure of the target through cross sections differential to the Mandelstam t variable, only experimentally available in exclusive diffractive events. Exclusive vector meson production in ep collisions, therefore, serves as a good probe of the gluonic radius of the proton in the transverse plane and for their fluctuations. The dipole amplitude in these models are usually parametrisations fitted to HERA reduced cross section measurements.
In [6] the authors introduced an explicit dependence on the target geometry in the form of a gluon thickness profile in the transverse plane. This profile was taken to be independent of collision energies. This geometrical description was further enhanced in [14] by introducing three hotspots of gluon density inside the proton, which were allowed to fluctuate. They were hence able to provide a geometrical description of the incoherent cross sections at HERA. However, these hotspots were still independent of the photon-proton collision energies in the events.
It is expected that the transverse gluonic radius of the proton will evolve and increase for small Bjorken-x as more low momentum gluons are radiated. The increase in the gluonic radius of the proton is compatible with Regge theory predictions, where the extracted coherent slope of the t-distribution increases logarithmically with increasing W γp [15,16]. This phenomenon is known as Gribov diffusion [17]. The analysis in [18] also supports the increase in hotspot width in hotspot models. Here the authors conclude that the transverse diffusion or growth of the hotspots with increasing collision energy is the primary dynamical process underlying the onset of the hollowness effect in pp interactions. Further, they ascribe the measured growth of the total pp cross section to this mechanism. Recent studies on diffractive vector meson production with the description of the proton using JIMWLK evolved Wilson lines [19,20] also indicate the increase in width of the fluctuating gluonic density hotspots with decreasing x. Further, the investigations with light nuclei in [21] suggest an evolution of the transverse width of the fluctuating gluon distribution, which is incorporated in nuclei as sub-nucleon fluctuations. They also find that the JIMWLK evolution predicts an increased deuteron size at small x. Additionally, the fluctuations are expected to evolve in energy and at high energies, there should be no event-by-event fluctuations in the black-disc limit where dipole cross section saturates to unity. The incoherent cross section, sensitive to the event-by-event fluctuations, is vastly suppressed in this limit as shown in [19]. Recent data on J/ψ photoproduction in ultra-peripheral collisions at low-x from the ALICE collaboration at LHC [22] also indicate the suppressed incoherent cross section compared to coherent at high energies. Though the energy dependence of the incoherent cross section is still not measured at AL-ICE, the p t distribution of the decay products of J/ψ is indicative of the suppressed incoherent cross section.
The dipole model provides a clean phenomenological environment to study the effects of energy evolution in proton geometry. It is, therefore, a fertile testing ground for implementing and comparing different aspects of energy dependence of the proton geometry and confront them with existing data as well as guide our understanding in what to expect from future colliders experiments such as the electron-ion collider (EIC) [23,24] and the large hadron-electron collider (LHeC) [25].
In this study, we investigate the energy dependence of the initial state of the proton in exclusive J/ψ, ρ, and φ meson production using the saturated impact-parameter dependent dipole model bSat (also known as IP-Sat), and its linearised version, the bNonSat model. We introduce a novel parametrisation of evolution effects in the proton geometry with and without subnucleon fluctuations. We do this by introducing an explicit x-dependence in the proton and hotspot widths. We also compare it with our implementation of the approach given in [26], in which the number of hotspots have an x-dependence.
The paper is organised as follows. In the next section, we give an outline of exclusive vector meson production in the dipole formalism, as well as of the hotspot model. We then introduce our modifications to these models to take energy dependence into account. In Section III we present the results and compare our predictions with the available HERA data. Finally, we summarise and discuss the main conclusions of the study.

II. The Color Dipole Models
In the dipole picture, the scattering amplitude for the diffractive vector meson production factorises at high energy and is given by the convolution of three subprocesses, as depicted in Fig. 1. First, the virtual photon splits into a quark anti-quark dipole; then, the dipole interacts with the proton via one or many colourless twogluon exchanges, after which it recombines into a vector meson. The amplitude is given by: where L and T represent the longitudinal and transverse polarisation of the virtual photon, r is the transverse size and direction of the dipole, z is the energy fraction of the photon taken by the quark, b is the impact parameter of the dipole relative to the proton and ∆ is the pomeron's transverse momentum where |∆| = √ −t. Here, (Ψ * Ψ V ) is the wave-function overlap between the virtual photon and the vector meson in the final state. The exponential factor has two terms, the b term comes from the impact parameter space while the r term is due to non-forward wave functions which calculated by Bartels, Golec-Biernat, and Peters in [27]. This BGBK factor is important for studying ρ and φ mesons at low Q 2 . The amplitude is a Fourier transform from the transverse coordinate of the quark in the dipole to the transverse momentum transfer in the interaction. The amplitude thus contains information on the spatial structure and fluctuations of the gluon density inside the proton. The dipole cross section dσ qq /d 2 b describes the strong interaction. We use the boosted Gaussian wave function with the parameter values from [9] for the vector-meson wave function. A recent update on different vector-meson wave functions can be found in [28]. The elastic diffractive cross section for a spherical proton (without geometrical fluctuations) is given by: We include initial state fluctuations in the proton's wave function by employing the Good-Walker formalism [29], where the coherent cross section is given by the first moment of the amplitude with respect to the initial state.
In contrast, the total cross section is given by its second moment. The incoherent cross section thus probes the difference between the second and first moment squared, which for Gaussian distributions is the variance of the initial state. Thus, for an event-by-event variation Ω, we have: We consider two versions of the dipole cross section. The bSat model dipole cross section is given by [30]: Due to the exponential functional form in this case, the dipole cross section saturates for large gluon densities xg(x, µ 2 ) and for large dipole sizes r. The scale at which the strong coupling α s and gluon density are evaluated is µ 2 = µ 2 0 + C r 2 and the gluon density at the initial scale µ 0 is parametrised as: where the parameters A g , λ g , C, m f are determined through fits to inclusive reduced cross section measurements. The bNonSat model is a linearised version of the bSat model where : which does not saturate for large gluon densities and large dipoles. This dipole cross section corresponds to a single two-gluon exchange. We use the fit results from [10], which also includes a Gaussian suppression of large dipoles. The transverse profile of the proton is usually assumed to be Gaussian: and the parameter B G is constrained through a fit to the t-dependence of the exclusive J/ψ production rates at HERA [7], and is found to be B G = 4 ± 0.4 GeV −2 , which is fixed for all x values. It should be noted that inclusive observables are not sensitive to the proton's profile, as this will only give contributions at higher twists. Therefore, as we modify the profile function in different ways in the rest of the paper it will not affect the quality of the fit for the bNonSat model at all, and only slightly for the bSat model. We have checked that the description of the HERA reduced cross section remains good in what follows.

A. The Fixed Hotspot Model
Event-by-event fluctuations of initial state gluon density inside the proton can be taken into account by the hotspot model [14,31,32]. Here, the gluons are assumed to be located in density hotspots. This can be implemented by changing the proton profile in eq.(7) in the following way: with which is Gaussian for S g = 0 and peaks more at the centre for non-zero values of S g [33]. Here N q = 3 is the number of the hotspots located at b i sampled from a Gaussian width B qc and B q is the width of the hotspots. B qc and B q control the amount of fluctuations in the proton geometry and are constrained by the coherent and incoherent data. For bNonSat, S g = 0, while for bSat we use S g = 0.3 for J/ψ-production and S g = 0.4 for ρand φ-production. We refer to [34] for a recent Bayesian analysis of the parameters of the hotspot model. Here in the hotspot model, the number and the size of the hotspots is fixed for all x I P values, thus resulting in a fixed width of the proton. We call this the Fixed Hotspot (FH) model. We also include the fluctuations in the saturation scale in our investigations following [31]. The experimentally observed multiplicity distributions and the rapidity correlations in pp collisions need these fluctuations in order to describe the data [35,36], the saturation scale fluctuations are incorporated by letting the saturation scale of hotspots fluctuate independently as follows: The saturation scale is Q 2 [6] in the dipole amplitude. We can implement these fluctuations by changing the normalisation of the profile function. We use σ = 0.4 for J/ψ-production and σ = 0.6 for ρ-and φproduction. These fluctuations play an important role for small |t| [14,33].
The differential cross section receives significant corrections (discussed in detail in [7,31]). Firstly, the dipole amplitude is approximated to be purely imaginary. However, the real part of the amplitude is taken into account by multiplying the crosssection by a factor of (1+β 2 ) with β = tan λπ/2), and λ = ∂ log(A γ * p→V p T,L )/∂ log(1/x I P ). Secondly, to take into account that the two gluons may have different momentum fractions, a skewedness correction to the amplitude is introduced [37], by a factor R g (λ) = 2 2λ+3 / √ π · Γ(λ g + 5/2)/Γ(λ g + 4) with λ g = ∂ log(x I P g(x I P ))/∂ log(1/x I P ). In our model, we calculate both these corrections using a spherical proton. They have significant contribution at low momentum transfer |t| (around 40-60 %).
Next we consider three separate modifications to the proton's gluon density distribution.

B.
Evolution effects on the proton's transverse size First, we explore these evolution effects in a spherical proton (without initial-state fluctuations). We introduce the evolution effects on the size of the proton by parametrising the width of the proton as a function of x I P . This is implemented by introducing an x I P dependence in the profile function by changing eq.(7) as follows: where B p determines the normalisation and λ p governs the evolution of the width of the proton in the transverse plane. As x I P = (M 2 V +Q 2 +|t|)/(W 2 γp +Q 2 −m 2 p ) (where M V and m p are the vector-meson and proton masses, respectively), the parameter values are determined through a fit to the t-spectra at different W γp values of exclusive J/ψ production rates at HERA. The gluonic radius for the spherical proton is r rms = 2B G (x I P ) which now varies with Bjorken-x.

C. Evolution effects on the hotspots' transverse size
To include the evolution effects of the proton's size in the hotspot model we make the width of hotspots x I P -dependent in two different ways. As x I P decreases we expect more emissions of small x I P gluons inside the hotspots, thus increasing the hotspot size in the impact parameter space. We first implement this by changing the profile function in eq.(8) as follows: having the same profile of hotspots as in eq.(9) with, where B hs determines the normalisation and λ hs governs the evolution of the width of hotspots and hence the width of the proton in the transverse plane. We call this model the Varying Hotspot Width (VHW) model. In the models where the nucleon size grows, for instance [19,20,41], the Froissart bound [42,43] limits the energy dependence of the radius of hadron to be logarithmic or weaker. Hence we also consider a parametrisation for width of hotspots as which results in a logarithmic increase of the radius of the proton. We will refer to this parametrisation as the 'logarithmic model'. A similar parametrisation with a quadratic dependence on the rapidity was considered to implement the growth of hotspots in a recent analysis [44] within the color glass condensate framework [45,46] to study the J/ψ production in pp and pPb collisions. We compare both the logarithmic and power law parametrisations for the growth of the hotspots' width to show that the power law growth with small values of the evolution parameter λ hs also follow the unitarity principle in the kinematic region considered in this study. The gluonic radius for the proton in these cases is given by r rms = 2(B qc + B q (x I P )). Thus the average gluonic radius of the proton also increases with decreasing x in addition to the increase in size of the hotspots. At large energies, as the size of the hotspots increase, the hotspots will begin to overlap, which will result in less fluctuations in the proton profile. Hence, at large energies, we expect a suppression of the incoherent cross section in these models.
The transverse profile of the proton with a fluctuating gluon distribution as a function of the centre of mass energy for the photon-proton system W γp in VHW model is illustrated in Fig. 2 (left).

D. Evolution effects on the number of hotspots
We next investigate another x I P -dependence of the hotspot model, for which the number of hotspots increases with decreasing x I P , while keeping the hotspots' width fixed. We refer to this model as the Varying Hotspot Number (VHN) model. This kind of a model was introduced in [26]. In that analysis the authors considered the GBW parametrisation of the dipole amplitude and assumed a factorised form for the thickness function. Following [26], the number of hotspots in the model is parametrised as a function of x I P . In this model the profile is the same as given in eq.(9) with the number of hotspots increasing stepwise at small x I P as: We implemented this prescription into the dipole models described above. The hypothesis underlying this implementation is that as Bjorken-x decreases, the gluon density increases steeply and instead of gluon accumulating in the original three hotspots of the FH model, the number of hotspots itself grows. The parameters p 0 , p 1 and p 2 are determined through a simultaneous fit to the differential and the total incoherent cross sections for J/ψ production. The VHN model is expected to exhibit less fluctuations at high energy as the numerous hotspots begin to overlap, which hence suppresses the incoherent cross section. The profile of the proton in the VHN model in illustrated in Fig. 2 (right).

III. Results
For a spherical proton with increasing transverse width, the parameters in the thickness function are de-  t-dependence of the differential cross section for exclusive J/ψ photo-production in the non-saturated (left) and saturated (right) FH, VHW, and logarithmic models fitted with the data from [1].
termined through a fit to the t-dependence of exclusive J/ψ cross sections at HERA at different W γp [1]. The optimal values of the parameters are found to be B p = 2.3 ± 0.03 GeV −2 and λ p = −0.062 ± 0.002 with a χ 2 /ndf = 48.06/48. The transverse gluonic radius of the proton now increases at small x I P for these parameter values.
In table I, we provide the values of all the parameters in the FH and VHW models with and without saturation, as well as for the logarithmic model. The parameters values are determined through a simultaneous fit to the t-dependence of coherent and incoherent vector meson HERA measurements. For J/ψ photo production, the fit is done to the available t-dependence data [1,11] in different bins of W in the kinematic region 45 ≤ W ≤ 251 GeV and 0.1 ≤ t ≤ 1.0 GeV 2 while for the ρ electro production the fit is done to the available data [3,13] in bins of Q 2 with 0.1 ≤ t ≤ 1.0 GeV 2 . Note that for rest of the plots the models prediction are compared with the data. For the bSat case we use a modified profile originally introduced in [33] to explain the coherent and incoherent data well at small momentum transfer |t|. The poor description of the ρ and φ data is partly due to tensions between the H1 and ZEUS measurements which makes a good fit to all available data difficult. We implement the VHN model with parameters p 0 = 0.011, p 1 = −0.56 and p 2 = 165 keeping the width of the hotspots the same as in the FH model in table I. Fig. 3 depicts the t-dependence (first row) and the energy dependence (second row) of the exclusive J/ψ photoproduction for a spherical proton, with and without evolution effects for different values of W γp . We obtain a similarly good agreement with the current data for the energy dependence in both the models with and without evolution effects. In the models with evolution effects, the total cross section increases at small W γp and is suppressed at large W γp . For the bSat case, the model with evolution effects seems to underestimate the LHCb data points.
In Fig. 4 we show the t-spectrum of the exclusive J/ψ photo-production in the saturated and non-saturated FH and VHW models as well as the logarithmic model. The t-spectrum description in the VHW and the logarithmic models are very similar. As the hotspot width increases, it leads to an increase in the average gluonic radius of the proton and as a result the slope of the t-spectrum increases at higher energies. We have not shown the VHN model predictions as they are identical to that of the FH model. We see that even though the shapes of the W γp and t spectra change, the overall description of the measurements remains equally good. We see that for small W γp , the t-slope decreases while for large W γp it increases, and the measured data seems to prefer a smaller slope overall. These effects off-set each other in the qual- for J/ψ photo-production in the saturated FH, VHW, VHN, and logarithmic models compared with the data from [11].
ity of the fits of the models. We further note that the bNonSat and bSat models give very similar predictions for the evolution effects on the size of the proton. This is consistent with earlier studies [9,33], which have not been able to separate the two models for describing available ep measurements. However, comparisons with UPC measurements with heavy ions at the LHC and RHIC experiments show a clear preference for the bSat model [10]. Hence, the rest of this study will be restricted to the bSat model.
In Fig. 5, we study the energy dependence of the incoherent cross section as well as the ratio of incoherent to coherent cross sections for J/ψ photo-production in the saturated FH, VHW and VHN models as well as the logarithmic model. Many model-dependent effects get cancelled in such a ratio, such as real-and skewednesscorrections and other model uncertainties, making it a clear measure of the energy-dependence of suppression of initial state fluctuations in the proton. In the VHW model we obtain a higher cross section at small W γp while the cross section is suppressed at large energies. As the hotspot size increases in the VHW model and the logarithmic model at small x I P , the hotspots begin to overlap and as a result the fluctuations are reduced at higher energies. The hotspot similarly begin to overlap at small x I P in the VHN model resulting in the suppressed incoherent cross section. We also plot the ratio of incoherent to coherent cross sections. In the VHW and the logarithmic models, this ratio decreases with x I P , while it remains nearly constant for the FH model. For the VHN model, this ration is first constant and then falls in discreet steps as the number of hotspots are constant, N q = 3 for 20 ≤ W γp ≤ 120 GeV, and then increases in integer steps, which results in the wobbly shape of the VHN curve in Fig.5. The current HERA data does not distinguish between these scenarios but theoretically there is a clear difference at large energies. Unlike reference [26], the incoherent cross section does not fall after reaching a maximum in this kinematic range, which is similar to what was found in the MV model with explicit JIMWLK evolution in [19].
In the left panel of Fig. 6, we plot the extracted coherent slope B D as a function of W γp in all the models. Here, B D is defined by fitting dσ/dt ∝ exp(−B D |t|) to the coherent t-distribution. This can be interpreted to correspond to the proton's effective transverse size. For the spherical proton with an energy-dependent width, and in the VHW model, the extracted slope clearly grows with W γp , faster for the latter. For the fixed sized spherical proton, as well as for the FH and VHN models, the transverse width remains constant even at higher energies. The logarithmic model represents the Froissart bound here and we see that all models lie below this for W γp 10 3 GeV which is where the power law of the VHW model crosses the B D value of the logarithmic model. The models' predictions are compared with the available HERA measurements which show a tendency that the average gluonic radius of the proton increases at higher energies. The H1 data shows a more pronounced increase in the size of the proton than the ZEUS data. In the right panel of Fig. 6 we test the dependence on the result on the choice of vector meson wave-function. We show the result from the logarithmic model using the default Boosted Gaussian wave-function as well as the Gauss-LC wave function [6]. We see that the uncertainty coming from the choice of wave-overlap is around 5%, and that it does not affect the resulting shape of the t-slope with respect to W γp , which is consistent with [7].
In Fig. 7, we study the energy dependence of the incoherent cross section and the ratio of incoherent to coherent cross sections for ρ and φ electro-production in the saturated FH, VHW and VHN models. As for J/ψ production, the incoherent cross sections are suppressed at high energies in the VHW and VHN models as compared to the FH model, but the available measurements are unable to distinguish between the models. The suppression of the incoherent cross section is more pronounced in the VHN model than in the other two. The ratio of incoherent to coherent cross section depicts the evolution of fluctuations with energy more clearly as the ratio decreases in the VHW and VHN models as we decrease x I P , while this ratio remains constant for the FH model. This is similar to what we observed for the energy dependence of the J/ψ meson.

IV. Conclusions and Discussion
We have considered three approaches for including energy dependence of the proton geometry into the bSat (also known as IP-Sat) and bNonSat dipole models for describing exclusive diffractive rates. Firstly, we considered the spherical proton, where we allowed the proton's gluonic width to increase with decreasing x I P , using a simple power parametrisation where the proton width becomes B G (x I P ) = B p x λp I P . We noted that this improved the descriptions of the t-spectrum for coherent J/ψ production at different W γp . However, as this model does not include any initial state fluctuations it cannot describe incoherent diffraction.
We also modified the hotspot model in three different ways. Firstly, in the VHW model, by allowing the hotspots' width to increase as B q (x I P ) = B hs x λ hs I P and in the logarithmic model as B q (x I P ) = b 0 ln 2 ( x0 x I P ). Thirdly, in the VHN model we allow the number of hotspots to increase at smaller x I P following the parametrisation of [26]. We note that the VHW and logarithmic models only need one new parameter, while the VHN model needs three. We show that these approaches give similar (but distinct) suppressions of the incoherent cross section at large W γp , as the hotspots begin to overlap. However, there is a clear difference between these approaches in the slope of the t-spectrum in exclusive diffraction. The slope increases with W γp in the VHW model, as well as in the logarithmic model where it grows similarly but slightly slower, while it remains constant in the VHN model.
Currently, the available measurements are unable to distinguish between these scenarios. As seen in table I, all models give a similar fit result. However, considering the W γp spectrum of incoherent diffraction and the incoherent to coherent ratio, as well as the t-slope show a preference for the VHW and logarithmic models. This situation may improve drastically with new measurements. The EIC will be able to measure with increased precision at the small-to-medium W γp region. These measurements would be able to resolve the tension between the H1 and ZEUS measurements which can be seen in Fig. 6, and help restrict the model parameters of the proton's profile further. The LHeC would be able to significantly extend the current phase space by probing the proton's gluonic radius and the incoherent suppression at large W γp , which will largely improve our understanding of the x-dependence of the proton's profile and the fluctuations spectrum. Similarly, current measurements of exclusive vector mesons in pA collisions in UPC at LHC and RHIC can potentially measure the incoherent to coherent ratio at large W γp in the near future. We note that the VHW model has an exponential dependence on rapidity, which will not be physical if the phase space of the analysis is extended, which makes the logarithmic model the strongest candidate for taking this study further.
One may also consider introducing a hybrid model in which both the hotspot width and number vary with the momentum fraction, or even vary around a mean for fixed x I P . Conceptually, the interpretation of such a model would become ambiguous. As the hotspots begin to overlap, the effective number of hotspots is not well defined, and the description with fewer colder hotspots becomes indistinguishable from having more and hotter hotspots. It is, therefore, phenomenologically clearer to keep the for ρ (left) and φ (right) electro-production in the saturated FH, VHW, VHN, and logarithmic models compared with the data from [13].
model of the proton geometry as simple as the experimental measurements allow. Thus, such a hybrid model would not be desirable.
In the future, it would be interesting to extend these models to AA ultra-peripheral collisions at the RHIC and LHC experiments, and to electron-ion collisions for the EIC, as we expect the hotspots' overlap to become more pronounced at lower W γp due to the heavy nucleus geometry. We also plan to implement these models into the event generator sartre framework [47,48].