Effects of next-to-leading order DGLAP evolution on generalized parton distributions of the proton and deeply virtual Compton scattering at high energy

We studied the effects of NLO Q2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q^2$$\end{document} evolution of generalized parton distributions (GPDs) using the aligned-jet model for the singlet quark and gluon GPDs at an initial evolution scale. We found that the skewness ratio for quarks is a slow logarithmic function of Q2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q^2$$\end{document}, reaching rS=1.5-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^\mathrm{S}=1.5{-}2$$\end{document} at Q2=100\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q^2=100$$\end{document} GeV2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} and rg≈1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^\mathrm{g} \approx 1$$\end{document} for gluons in a wide range of Q2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q^2$$\end{document}. Using the resulting GPDs, we calculated the DVCS cross section on the proton in NLO pQCD and found that this model in conjunction with modern parameterizations of proton PDFs (CJ15 and CT14) provides a good description of the available H1 and ZEUS data in a wide kinematic range.


Introduction
Generalized parton distributions (GPDs) have become a familiar and standard tool of Quantum Chromodynamics (QCD) describing the response of hadronic targets in various hard exclusive processes [1][2][3][4][5][6][7][8]. GPDs can be rigorously defined in the framework of QCD collinear factorization for hard exclusive processes [9,10], which allows one to access universal, i.e., process-independent, GPDs in such processes as deeply virtual Compton scattering (DVCS) γ * + T → γ + T , timelike Compton scattering (TCS) γ + T → γ * + T , exclusive meson production by longitudinally polarized photons γ * L + T → M + T , and, recently, photoproduction of heavy (J/ψ, Υ ) vector mesons γ + T → V + T [11,12]. GPDs contain information on the hadron structure in QCD, which is hybrid of that encoded in usual parton distributions and elastic form factors. In particular, GPDs describe the distributions of quarks and gluons in hadrons in terms of two light-cone momentum fractions and the position in the transverse plane. Also, GPDs are involved a e-mail: guzey_va@pnpi.nrcki.ru in the hadron spin decomposition in terms of the helicity and orbital motion contributions of quarks and gluons [4][5][6][7][8], and carry information on the spatial distribution of forces experienced by partons inside hadrons [13].
GPDs are essentially non-perturbative quantities, which cannot be calculated from the first principles apart from first Mellin moments in special cases in lattice QCD [14,15]. At the same time, evolution of GPDs with an increase of the resolution scale Q 2 is predicted by the QCD Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations modified to the case of GPDs, which are presently known to the next-to-leading order (NLO) accuracy [16][17][18]. Therefore, one of directions of phenomenological studies of GPDs is to determine the non-perturbative input for these evolution equations. After early studies of GPDs using various dynamical models of the nucleon structure [19][20][21][22][23][24][25][26][27], one currently focuses on parameterizations of GPDs, which are determined from fitting the available data. The two main contemporary approaches include the flexible parameterization based on the conformal expansion of GPDs [28][29][30][31] and global fits of GPDs [32][33][34][35], which use the double distribution (DD) model [36][37][38][39][40] in the Vanderhaeghen-Guichon-Guidal (VGG) framework; see for details [34]. One should also mention a pioneering study of global QCD fits of GPDs within the neural network approach [41].
The mentioned above analyses present only a partial, model-dependent picture of GPDs in a limited kinematic range. For further progress, it is important to perform a systematic QCD analysis of evolution of GPDs and cross sections of hard exclusive processes involving them. It will enable one to separate the effects of non-perturbative input GPDs from the perturbative DGLAP evolution and help to explore the possibility to use the data on hard exclusive reactions at high energies for constraining GPDs; see, e.g., [42].
In this paper, we calculate the effect of next-to-leading (NLO) QCD evolution on quark and gluon GPDs of the proton using the brute-force evolution method of [16][17][18] and the physical model for input GPDs, which is motivated by the aligned-jet model [27]. Using the obtained results, we calculate the DVCS cross section on the proton in NLO QCD and compare it to the available HERA data. We find that our approach provides a good description of the DVCS data over a wide kinematic range, including most of the data from H1 and ZEUS collaborations for the unpolarized proton target.
2 Aligned-jet model for GPDs and QCD evolution effects

Input GPDs
The aligned-jet model (AJM) [43,44] for photon-hadron interactions at high energies is based on the general observation that in the target rest frame, the incoming photon first fluctuates into quark-antiquark configurations, which then interact with the target. For the photon virtualities Q 2 = O(few) GeV 2 , the qq pair (dipole) is characterized by a small relative transverse momentum (hence the name aligned-jet), the invariant mass of the order of Q 2 , the asymmetric sharing of the photon's light-cone momentum, and the dipole-nucleon cross section, which has the magnitude typical for hadron-nucleon cross sections. Note that in QCD, this parton picture is complemented by the gluon emission and the contribution of quark-antiquark dipoles with large relative transverse momenta, which become progressively important as Q 2 is increased; see the discussion in Ref. [45]. In the AJM model, one obtains for the ratio of the imaginary parts of DVCS and DIS amplitudes at Q 2 = 1−3 GeV 2 , R = T DVCS / T DIS = 2.5−3.5 [27,46], which agrees nicely with the values of R extracted from the HERA data [47]. This in turn means that the effect of skewness of the singlet quark GPDs in the DGLAP region of X ≥ ζ can be neglected (X is the light-cone momentum fraction of the target in the initial state carried by the interacting parton; ζ is the momentum fraction difference between the two interacting partons). This observation is also supported by the analysis of Ref. [6], which showed that the good description of the high-energy HERA data on the DVCS cross section on the proton can be achieved with the forward parton distribution model for the singlet quark GPDs [48,49], i.e., with the δ-function-like profile in the DD model for sea quark GPDs.
In general, modeling and parametrization of GPDs is a non-trivial task since GPDs should satisfy several general constraints: GPDs reduce to usual parton distributions functions (PDFs) in the forward limit; integration of GPDs over the momentum fraction gives the corresponding elastic form factors; as a consequence of Lorentz invariance, Mellin moments of GPDs are finite-order polynomials in even pow-ers of the skewness η = ζ /(2−ζ ) (the property of polynomiality); GPDs obey positivity bounds expressed an inequalities involving GPDs and usual PDFs. While the first three properties can be naturally implemented in momentum representation of GPDs, positivity is most naturally derived in coordinate representation. Hence, it is an outstanding challenge to propose a practical model of GPDs satisfying all these constraints. (Naturally, field-theoretical approaches based on perturbative diagrams will automatically lead to GPDs satisfying all the constraints [26], but they have little usefulness for GPD phenomenology.) Starting from a model for GPDs in the DGLAP region of X ≥ ζ , there is no unique and simple way to reconstruct GPDs in the entire range of X . For instance, the method proposed in [27,47] does not guarantee polynomiality for higher moments of GPDs and conflicts with dispersion relations (DR) for the real and imaginary parts of the DVCS amplitude [50]. In principle, GPDs with the correct forward limit and satisfying the property of polynomiality can be constructed using the so-called Shuvaev transform [51][52][53][54]. However, this method is usually associated with the leading order (LO) phenomenology and brings about a certain skewness dependence of GPDs in the DGLAP region. Similarly, the flexible parameterization of GPDs based on the conformal expansion [28][29][30][31] contains the skewness effect of GPDs in the DGLAP region and corresponds to model-dependent parton distributions in the forward limit.
In this work, to simultaneously have the forward-like GPDs in the DGLAP region and circumvent the aforementioned problem with polynomiality, we take forward-like GPDs for all X and add the so-called D-term [55], which has support only in the Efremov-Raduyshkin-Brodsky-Lepage (ERBL) region of |X | ≤ ζ . Specifically, we use the following model for the singlet quark (one sums over quark flavors q) and gluon GPDs at t = 0 at the initial scale of μ 0 : q(x, μ) and g(x, μ) are the quark and gluon parton distribution functions (PDFs), respectively. Note that, since we explicitly introduced antiquark GPDs, it is sufficient to consider only non-negative X ≥ 0. Also, we assume that similar relations hold for separate quark flavors q, i.e., where n f is the number of active quark flavors. Thus, our model does not assume the flavor symmetry of quark GPDs. We should stress here that GPDs have to be continuous at X = ζ . In addition, GPDs have to satisfy the correct symmetries around the midpoint of the ERBL region, X = ζ /2. As follows from the general properties of GPDs, the singlet quark GPDs H S X, ζ, t = 0, μ 0 ) is antisymmetric in the ERBL region around the X = ζ /2 point, while the gluon GPD H g (X, ζ, t = 0, μ 0 ) is symmetric in the ERBL region; these constraints are implemented in Eqs. (1) and (2).
The function D S (x/η) is the singlet quark D-term [55], which can be expanded in terms of odd Gegenbauer polynomials C 3/2 n in the following form [56]: The coefficients d 1 , d 3 and d 5 were estimated in the chiral quark soliton model at μ 0 = 0.6 GeV in Ref. [20]: Note that due to the lack of numerical estimates, we neglected the possible gluon D-term in Eq. (1). In this case, D S (z, μ) evolves in μ 2 autonomously (without mixing) and its value for μ > μ 0 can be readily calculated.
In summary, our GPD model in Eqs. (1) and (2) corresponds to the correct forward limit, satisfies polynomiality (the D-term satisfies polynomiality by construction), and obeys positivity bounds in the DGLAP region in the small-ξ and t = 0 limit (all positivity bounds discussed in the literature are for the DGLAP region; see Ref. [7]). Indeed, neglecting the kinematically suppressed contribution of the GPD E, the positivity bound for the quark GPDs reads [7] ( where (4) is trivially satisfied with our GPD model of Eq. (2). By construction, see Eq.
(1), in the middle of the ERBL region at x = X − ζ /2 = 0, our singlet quark GPDs become singular and the gluon GPD vanishes. Being a natural artifact of our model imposing the correct GPD symmetry in the ERBL problem, it does not violate general principles of GPDs, does not conflict with factorization for amplitudes of hard exclusive processes, and does not lead to singularities of the DVCS amplitude. Since the main goal of our work is to study the effects of NLO Q 2 evolution of GPDs in conjunction with different baseline PDFs, the simple model of Eq. (1) should suffice.
Note that in this work, we focus on the quark singlet q (q +q) and gluon GPDs: valence (non-singlet) quark GPDs do not mix with singlet quark and gluon GPDs under the DGLAP evolution and do not appreciably contribute to the DVCS amplitude at high energies.

NLO Q 2 evolution of GPDs and error analysis
The determination of parton distribution functions (PDFs) has always been one of the important ingredients for theory predictions. In this respect, more accurate PDFs play an important role in understanding of hadronic properties and the structure of the nucleon [57][58][59][60]. Our knowledge of PDFs has developed both theoretically and computationally. However, results of various groups lead to different predictions of physical observables. As we know, GPDs are quantities that are related to the PDFs in the forward limit and in many phenomenological approaches. To investigate the impact of different PDFs on the GPDs and their evolution, we calculate the effect of next-to-leading order (NLO) DGLAP evolution equations modified to the case of GPDs using the formalism of [16][17][18] and the input GPDs of Eqs. (1). (The early results on leading order (LO) Q 2 evolution of GPDs were presented in Refs. [51,61]). Perturbation theory predicts the evolution of GPDs and, hence, they depend on the factorization scale, μ 2 . Anomalous dimensions and the kernels at NLO accuracy in pQCD can be found in Refs. [62][63][64][65][66].
For the forward PDFs, we used CT14 [67] and the new CTEQ-Jefferson Lab (CJ15) analysis [68]. To study the impact of PDF uncertainties on the GPD evolutions and DVCS cross sections, we include the uncertainties of CT14 and CJ15 PDFs in the calculations of the evolution and in the DVCS cross sections. In this respect, we note that both CT14 and CJ15 are PDF sets with Hessian PDF eigenvector error sets. In this situation, the theoretical uncertainties of PDFs themselves and any physical quantity related to them, such as the GPDs and DVCS cross sections considered here, can be obtained as usual using the 56 and 48 error sets of the CT14 and CJ15 parametrizations, respectively. To this aim, we must first calculate our desired quantity with various error sets. Then we can compute the deviations from the central result and so the contribution to the size of the upper and lower errors through the following relations [69,70]: The last point that should be noted here is the confidence region considered for estimating error bands since different PDF analyses typically utilize different criteria for estimating PDF errors. The CJ15 PDF sets have been provided with 90% C.L. uncertainties considering standard tolerance criterion Δχ 2 = 2.71, while CT14 use a tolerance criterion as Δχ 2 = 100 with the same confidence level [71]. In this work, we display the CT14 and CJ15 errors on GPDs and DVCS cross sections for 90% C.L. region, so that the tolerance used for CJ15 PDFs be matched with CT14, in order to have a reasonable comparison. Figures 1 and 2 show the results for the singlet quark GPD H S (X, ζ, t = 0, Q 2 ) and the gluon GPD H g (X, ζ, t = 0, Q 2 ), respectively, as a function of X at ζ = 0.001 and Q 2 = 1.69, 4, 10, and 100 GeV 2 . Note that Q 2 = 1.69 GeV 2 is the input scale for CT14 and CJ15. As can be seen from these figures, the Q 2 evolution pushes GPDs into the ERBL region of X < ζ as it should be. The discontinuity of the quark singlet GPD at X = ζ /2 is an artifact of our model (see the discussion in Sect. 2.1), which does not affect the physical observables.
In the quark singlet case, the difference between the predictions based on CT14 and CJ15 PDF is small, especially at lower values of the Q 2 resolution scale. At the same time, in the gluon channel the differences between the CT14 and CJ15 predictions are sizable and exceed the associated uncertainties for large values of Q 2 . One should also note that the uncertainties of the resulting GPDs based on CT14 are larger than those for CJ15, which is related to the large uncertainties of CT14 singlet distributions at small x. Generally speaking, our results indicate that the GPD model of Eq. (1) is sensitive to the input PDFs. Therefore, more accurate PDFs are very important for physical observables involving GPDs such as, e.g., DVCS cross sections. Conversely and optimistically, data on the DVCS cross section may provide new constraints for global QCD analysis of PDFs. Our study makes it clear that using more recent version of PDFs and proper scale dependence in our GPDs model describes DVCS data over a large kinematical range.

Effect of skewness
For phenomenological applications of GPDs, it is important to discuss the so-called skewness factor, which describes the connection between GPDs and PDFs and parametrizes the deviation of GPDs from PDFs. To quantify this effect, it is convenient to introduce the following ratios of quark and gluon GPDs and PDFs [29]: Our results for r S (ζ, μ) and r g (ζ, μ) as functions Q 2 = μ 2 at ζ = 0.001 are shown in Fig. 3. One can see from the figure that both r S and r g are slow logarithmic functions of Q 2 . By construction, r S = r g = 1 at the initial evolution scale of Q 2 = 1.69 GeV 2 . As Q 2 is increased, r S slowly increases up to r S ≈ 1.5-2 at Q 2 = 100 GeV 2 , while r g stays at the level of unity for the studied range of Q 2 . These results agree with the predictions of the flexible GPD parameterization based on the conformal expansion, see Fig. 7 of Ref. [29], except for r S at the input Q 2 = 1.69 GeV 2 , where our result lies lower than that of [29].

Evaluation of Compton form factors and DVCS amplitudes
The standard and well-tested way to access GPDs is the process of leptoproduction of a real photon, ep → eγ p, or deeply virtual Compton scattering (DVCS). At the photon level, the γ * p → γ p DVCS differential cross section, dσ DVCS (W, t, Q 2 )/dt, is expressed in terms of the so-called Compton form factors (CFFs), which in the collinear factorization approach [9] are given as convolution of the perturbatively calculable hard scattering coefficient functions with the non-perturbative GPDs. In particular, at high energies the DVCS cross section is by far dominated by the GPD H and its corresponding CFFs. For the flavor singlet contribution (for the quark singlet and gluon CFFs), one has in the symmetric notation: where ξ = ζ /(2 − ζ ); μ is the factorization scale which is usually set equal to the photon virtuality μ 2 = Q 2 . The explicit form for the coefficient functions C can be found in Refs. [28,72,73] for the non-singlet and singlet cases. For instance, for the quark singlet case, the QCD perturbation series reads Hence, to the LO accuracy of pQCD and in leading-twist approximation, the DVCS scattering amplitude (CFF) can be written as The CFFs depend on ξ (or equivalently on Bjorken x B or the invariant energy W ), the momentum transfer t, and Q 2 and, hence, can be extracted from DVCS experiments. Note that our model for the GPD initial conditions does not imply the flavor symmetry of quark GPDs. Note that the purely electromagnetic Bethe-Heitler (BH) bremsstrahlung process leads to the same final state and interferes with DVCS. However, at high energies (small values of x B ), the DVCS process dominates, which allows one to extract the DVCS cross section by subtracting of the BH contribution. In addition to this, cuts on Q 2 and W have been applied by H1 and ZEUS collaborations to enhance the contribution from DVCS process (see the following subsection).
Detailed analytic expressions for the DVCS and BH amplitudes squared and their interference are well known and can be found in Refs. [74,75]. In our analysis we assume an exponential and factorized t-dependence of the DVCS cross section, with a = 8 GeV −2 and c = 0.15 [27]. This simple parametrization agrees with the measurements of the t dependence of the differential γ * p → γ p cross section at HERA (see the following subsection).

HERA DVCS data
Unpolarized DVCS on the proton has been measured in e ± p collisions at HERA by the H1 [76][77][78][79] and ZEUS [80,81] experiments. The list of the DVCS experiments at HERA along with the measured observables, kinematic ranges and corresponding references is given in Table 1.  [76], where the statistical and systematical errors are added in quadrature, is compared to our NLO pQCD results based on the input of Eq. (1) and CT14 [67] and CJ15 [68] PDFs. The shadowed bands represent the uncertainty of the corresponding PDFs    [76][77][78][79] and ZEUS [80,81] measurements (see Table 1). The error bars the statistical and systematic uncertainties added in quadrature. The bands associated with CJ15 and CT14 prediction correspond to the uncertainty of the respective PDFs.
One can see from these figures that within experimental and theoretical uncertainties, the input GPD model based on the CJ15 fit provides a good description of the H1-2001, H1-2005, H1-2007 and H1-2009 data (Q 2 dependence only for the two latter data sets), while the model based on the CT14 fit tends to somewhat overestimate the cross section normalization (it describes well the W dependence of the H1-2005, H1-2007 and H1-2009 data). At the same time, the CT14 parametrization leads to a very good description of the ZEUS data. These results clearly show that, for some selected PDF sets, such as, e.g., the CJ15 and CT14 fits, the AJM GPD model of [27] together with NLO pQCD calculations describes well the high-energy DVCS cross section.
In order to study effects of the NLO DGLAP evolution on GPDs, a detailed comparison of our obtained results with the DVCS γ * p → γ p cross section is shown in Fig. 10.  The DVCS γ * p → γ p cross section as a function of Q 2 (left) and W (right). Our NLO pQCD results are compared to the 2008 ZEUS data [81]; see for details Fig. 4 In this figure, we show the DVCS cross section as a function of W for some selected values of Q 2 = 2.4, 6.2, 9.9 and 18 GeV 2 . Our NLO pQCD predictions are based on the CT14 [67] PDFs; the experimental points are the 2003 and 2008 ZEUS data [80,81]. The inner error bars represent the statistical, and the full error bars the quadratic sum of the statistical and systematic uncertainties. One can see that a very good agreement between our predictions and ZEUS data is achieved for a wide range of Q 2 and W . It illustrates an important role of the Q 2 -dependence of the quark and gluon GPD H for the successful description of the HERA data, which spans a wide range of Q 2 and W .
In summary, we observe very good overall agreement between our NLO pQCD predictions and most of the H1 and ZEUS data. It warrants the application of our framework to forthcoming and planned DVCS measurements at high energies, such as, e.g., at COMPASS at CERN [82], an Electron-Ion Collider (EIC) [42], and the Large Hadron-Electron Collider (LHeC) [83] or Future Circular Collider (FCC-he).
Moreover, using our GPD model as a baseline, one can perform a global fit of all available H1 and ZEUS DVCS data (cross section and its asymmetries) as well as the data from other DVCS experiments with fixed proton targets (HER-MES, JLab) (the latter will require extension of our model Fig. 10 The DVCS γ * p → γ p cross section as a function of W for selected values of Q 2 = 2.4, 6.2, 9.9, and 18 GeV 2 . The NLO pQCD predictions based on the GPD model of Eq. (1) along with the CT14 [67] PDFs are compared to the 2003 and 2008 ZEUS data [80,81] to the remaining GPDs). It is worth mentioning here that all known constraints on GPDs presented in Sect. 2.1 cause the reduction of flexibility of choosing a proper GPDs functional from. The success of global fits existing in the literature (see the brief discussion of these in Introduction) as well as any future attempts to global fitting procedures strongly depends on the choice of data sets and the functional form of the GPDs. Therefore, any advances both in theory and experiments in these regards are most welcome.

Conclusions
The DVCS process is the golden channel to access GPDs and potentially extract them from the experimental observables. Taking advantage of the high invariant energy available in lepton-proton collisions at HERA, the H1 and ZEUS measured the DVCS cross section in a wide kinematic range and studied precisely its dependence on Q 2 , W , and t. These measurements covered the Bjorken x range of 10 −4 < x B < 10 −2 , where sea quarks and gluons dominate. These data sets provide valuable information for GPDs phenomenology and several groups have attempted to extract CFFs and GPDs using them.
In this work, we studied the effects of NLO Q 2 evolution of GPDs using a model for the singlet quark and gluon GPDs at an initial evolution scale motivated by the aligned-jet model of photon-hadron interactions at high energies. Quantifying the skewness and evolution effects by the GPD-to-PDF ratios r S and r g , we found that r S increases logarithmically slowly from r S = 1 at the input scale of Q 2 = 1.69 GeV 2 to r S = 1.5-2 at Q 2 = 100 GeV 2 ; in the gluon channel, r g ≈ 1 for the studied range of Q 2 . This observation agrees with the results of the more sophisticated model of GPDs based on conformal expansion [29].
Using the resulting GPDs, we calculated the DVCS cross section on the proton in NLO pQCD and compared it to the available HERA data. We found that our simple physical model of input GPDs used in conjunction with two modern parameterizations of proton PDFs (CJ15 and CT14) provides good description of the H1 and ZEUS data. It demonstrates that our GPDs model is reliable and flexible enough to use in fitting procedures using a variety of data sets.