Border and skewness functions from a leading order fit to DVCS data

We propose new parameterizations for the border and skewness functions appearing in the description of 3D nucleon structure in the language of Generalized Parton Distributions (GPDs). These parameterizations are constructed in a way to fulfill the basic properties of GPDs, like their reduction to Parton Density Functions and Elastic Form Factors. They also rely on the power behavior of GPDs in the $x \to 1$ limit and the propounded analyticity property of Mellin moments of GPDs. We evaluate Compton Form Factors (CFFs), the sub-amplitudes of the Deeply Virtual Compton Scattering (DVCS) process, at the leading order and leading twist accuracy. We constrain the restricted number of free parameters of these new parameterizations in a global CFF analysis of almost all existing proton DVCS measurements. The fit is performed within the PARTONS framework, being the modern tool for generic GPD studies. A distinctive feature of this CFF fit is the careful propagation of uncertainties based on the replica method. The fit results genuinely permit nucleon tomography and may give some insight into the distribution of forces acting on partons.


Introduction
Fifty years after the discovery of quarks at SLAC [1,2], understanding of how partons form a complex object such as the nucleon still remains among the main challenges of nuclear and high energy physics. In the last twenty years we have witnessed a new liveliness in the field of QCD approaches to this problem due to the discovery of Generalized Parton Distributions (GPDs) [3][4][5][6][7]. GPDs draw so much attention because of the wealth of new information they contain. Namely, GPDs allow for the so-called nucleon tomography [8][9][10], which is used to study a spacial distribution of partons in the plane perpendicular to the nucleon motion as a function of parton longitudinal momenta. Before, positions and longitudinal momenta of partons were studied without any connection through other yet less complex non-perturbative QCD objects: Elastic Form Factors (EFFs) and Parton Distribution Functions (PDFs). In addition, GPDs have another unique property, namely they are connected to the QCD energymomentum tensor of the nucleon. This allows for an evaluation of the contribution of orbital angular momentum of quarks to the nucleon spin through the so-called Ji's sum rule [4,5]. This energy-momentum tensor may also help to define "mechanical properties" and describe the distribution of forces inside the nucleon [11,12].
It was recognized from the beginning that Deeply Virtual Compton Scattering (DVCS) is one of the cleanest probes of GPDs. The first measurements of DVCS by HERMES [13] at DESY and by CLAS [14] at JLab have proved the usability of GPD formalism to interpret existing measurements, and have established a global experimental campaign for GPDs. Indeed, nowadays measurements of exclusive processes are among the main goals of experimental programs carried out worldwide by a new generation of experiments -those already running, like Hall A and CLAS at JLab upgraded to 12 GeV and COMPASS-II at CERN, and those foreseen in the future, like Electron Ion Collider (EIC) and Large Hadron Electron Collider (LHeC). Such a vivid experimental status is complemented by a significant progress in the theoretical description of DVCS. In particular, such new developments like NLO [15][16][17][18][19][20][21][22], finite-t and mass [23] corrections are now available. Except DVCS, a variety of other exclusive processes has been described to provide access to GPDs, in particular: Timelike Compton Scattering [24], Deeply Virtual Meson Production [25], Heavy Vector Meson Production [26], Double Deeply Virtual Compton Scattering [27,28], two particles [29,30] and neutrino induced exclusive reactions [31][32][33]. For some of those processes experimental data have been already collected, while other processes are expected to be probed in the future.
The phenomenology of GPDs is much more involved than that of EFFs and PDFs. It comes from the fact that GPDs are functions of three variables, entering observables in nontrivial convolutions with coefficient functions. In addition, GPDs are separately defined for each possible combination of parton and nucleon helicities, resulting in a plenitude of objects to be constrained at the same time. This fully justifies the need for a global analysis, where a variety of observables coming from experiments covering complementary kinematic ranges is simultaneously analyzed. So far, such analyzes have been done mainly for Compton Form Factors (CFFs), being DVCS sub-amplitudes and the most basic GPDsensitive quantities as one can unambiguously extract from the experimental data. Recent analyzes include local fits [34,35], where CFFs are independently extracted in each available bin of data, and global fits [36], where CFFs parameterizations are constrained in the whole available phase-space. For a review of DVCS phenomenology we direct the reader to Ref. [37].
The aim of this analysis is the global extraction of CFFs from the available proton DVCS data obtained by Hall A, CLAS, HERMES and COMPASS experiments. We use the fixed-t dispersion relation technique [38] for the evaluation of CFFs at the Leading Order (LO) and Leading Twist (LT) accuracy. For a given CFF, the dispersion relation together with the analytical regularization techniques requires two components: i ) the GPD at ξ = 0, and ii ) the skewness ratio at x = ξ. Ansätze for those two quantities proposed in our analysis accumulate information encoded in available PDF and EFF parameterizations, and use theory developments like the x → 1 behavior of GPDs [39]. They allow to determine a border function [40,41], being a GPD of reduced kinematic dependency x = ξ, and the subtraction constant, directly related to the energy-momentum tensor of the nucleon.
Our original approach allows to utilize many basic properties of GPDs at the level of CFFs fits. We analyze PDFs, but also EFF and DVCS data, that is, we combine information coming from (semi-)inclusive, elastic and exclusive measurements. The analysis is characterized by a careful propagation of uncertainties coming from all those sources, which we achieved with the replica method. Obtained results allow for nucleon tomography, while the extracted subtraction constant may give some insight into the distribution of forces acting on partons inside the nucleon.
This work is done with PARTONS [42] that is the open-source software framework for the phenomenology of GPDs. It serves not only as the main component of the fit machinery, but it is also utilized to handle multithreading computations and MySQL databases to store and retrieve experimental data. PARTONS is also used for the purpose of comparing existing models with the results of this analysis. This paper is organized as follows. Section 2 is a brief introduction of GPDs, DVCS and related observables, with details on the evaluation of CFFs given in Sec. 3. Ansätze for the border and skewness functions are introduced in Sec. 4. Sections 5 and 6 summarize our analyses of PDFs and EFFs, respectively. DVCS data used in this work are specified in Sec. 7. In Sec. 8 the propagation of uncertainties is discussed, while the results are given in Sec. 9. In Sec. 10 we summarize the content of this paper.

Theoretical framework
In this section a brief introduction to the GPD formalism is given. We emphasize the role of quark GPDs, as only those contribute to DVCS at LO. A deep understanding of the basic features of the contributing GPDs is crucial for constructing parameterizations of CFFs. More involved tools, like nucleon tomography, are important for the exploration of the partonic structure of the proton. This section also provides a foundation to DVCS description and illustrates the construction of observables used in our fits.
For brevity, we suppress in the following the dependence on the factorization and renormalization scales, µ 2 R and µ 2 F , which in this analysis are identified with the hard scale of the process Q 2 . A detailed preface to the GPD formalism may be found in one of available reviews [43][44][45][46].

Generalized Parton Distributions
In the following we use the convention for the light cone vectors as in Ref. [44]. In the light cone gauge, quark GPDs for a spin-1 /2 hadron are defined by the following matrix elements: Here, x is the average longitudinal momentum of the active quark, ξ = −∆ + /(2P + ) is the skewness variable and t = ∆ 2 is the square of four-momentum transfer to the hadron target, with the average hadron momentum P obeying P 2 = m 2 − t/4, where m is the hadron mass. In this definition the usual convention is used, where the plus-component refers to the projection of any four-vector on a light-like vector n.
With the help of the Dirac spinor bilinears: which are normalized so thatū(p)γ µ u(p) = 2p µ , one can decompose F q and F q into two pairs of chiraleven GPDs: recognized as two "unpolarized" GPDs H q and E q , and two "polarized" GPDs H q and E q . The relation with one-dimensional PDFs and EFFs is essential for the phenomenology of GPDs. In the forward limit of ξ = t = 0, when both the hadron and the active quark are untouched, certain GPDs reduce to (one-dimensional) PDFs: where q(x) and ∆q(x) are the unpolarized and polarized PDFs, respectively. No similar relations exist for the GPDs E q and E q that decouple from the forward limit. The relation to EFFs can be obtained by integrating GPDs over the partonic variable x: where F q 1 (t), F q 2 (t), g q A (t) and g q P (t) are the contribution of the quark flavor q to the Dirac, Pauli, axial and pseudoscalar EFFs, respectively.
The integrals in Eqs. (11)- (14) do not depend on ξ as a consequence of the Lorentz covariance of GPDs. This feature is generally expressed by a nontrivial property of GPDs known as polynomiality. The property states, that any n-th Mellin moment of a given GPD is always an even polynomial in ξ, of order n + 1 for the unpolarized GPDs and of order n for the polarized GPDs. In particular: where for n = 0 one has the relations given by Eqs.
(11)- (14). The correspondence of GPDs to PDFs and EFFs presages a possibility of studying a spatial distribution of partons inside the nucleon. Indeed, the subfield of hadron structure studies known as nucleon tomography allows one to extract the density of partons carrying a given fraction of the nucleon longitudinal momentum x as a function of the position b ⊥ in the plane perpendicular to the nucleon motion. For unpolarized partons inside an unpolarized nucleon this density is expressed by: where we stress the condition ξ = 0, meaning no change of the longitudinal momentum of the active parton. This density gets distorted when the nucleon is polarized. This effect is described by adding extra terms related to the GPDs H and E. The longitudinal polarization of partons distributed in a longitudinally polarized nucleon according to q(x, b ⊥ ) can be studied with the Fourier transform of GPD H: A representation in the impact parameter space is also possible for the GPD E: A probabilistic interpretation of that result is possible if one changes the basis from longitudinal to transverse polarization states of the nucleon [9]. In such a case, e q (x, b ⊥ ) can be related to a shift of parton density generated in a transversely polarized nucleon.
As indicated in Refs. [9,47], b ⊥ is the distance between the active parton and the point determined by the positions of individual partons weighted by their momenta, so that x i b ⊥,i = 0, where the sum runs over all partons (the struck parton and all spectators as well). This distance is different than that between the struck parton and the spectator system, which is given by: and the parton distribution given as a function of d ⊥ provides a better estimation of the transverse proton size than q(x, b ⊥ ). Another useful quantity appearing in the context of the nucleon tomography is the normalized second moment of q(x, b ⊥ ) distribution given by: A similar quantity can be also defined for the distribution of longitudinal polarization to check how broad this distribution is and how it corresponds to The need of having the proton size finite requires to keep the mean squared distance between the active parton and the spectator system, finite as well, also in the limit of x → 1. We will impose it for the valance quarks as an extra constraint on our Ansatz introduced in Sec. 4.
To avoid a violation of the positivity of parton densities in the impact parameter space, inequalities studied in a series of papers [48][49][50][51][52][53][54][55] must hold. In particular one has the following inequalities, which are proved to be useful to constrain parameterizations of GPDs: For completeness we also show Ji's sum rule, allowing for the evaluation of total angular momentum carried by partons: This feature can be used to investigate the nucleon spin decomposition. We note however, that this analysis concentrates on GPDs H and H, and therefore we will not attempt to give any estimation on J q .

Deeply Virtual Compton Scattering
A prominent role in the GPD phenomenology is played by Deeply Virtual Compton Scattering: where l, N and γ denote lepton, nucleon and produced photon, respectively; the four-vectors of these states appear between parenthesis. Under specific kinematic conditions, the factorization theorem allows one to express the DVCS amplitude as a convolution of the hard scattering part, being calculable within the perturbative QCD approach, and GPDs, describing an emission of parton from the nucleon and its subsequent reabsorption, see Fig. 1. The factorization applies in the Bjorken limit and for −t/Q 2 1, where Q 2 = −(k − k ) 2 is the virtuality of the virtual-photon mediating the exchange of four-momentum between lepton and proton at Born order.
The scattering is described by two angles, see The angle between the leptonic plane (spanned by the incoming and outgoing lepton momenta) and the production plane (spanned by the virtual and outgoing photon momenta) is denoted by φ. The angle between the leptonic plane and the nucleon polarization vector is denoted by φ S . angle between the lepton scattering plane and the direction of transversely polarized target.
An electromagnetic process called Bethe-Heitler has the same initial and final states as DVCS. The total amplitude for single photon production, T , is then expressed by a sum of amplitudes for BH and DVCS processes, which leads to: where I denotes the interference term. The cross section for BH is calculable to a high degree of precision and therefore can be easily taken into account in analyses of experimental data. The interference term provides a complementary information to the pure DVCS cross section, and in a certain kinematic domain allows to access GPDs, even if |T DVCS | 2 is small.
The amplitudes T DVCS and I may be expressed by combinations of CFFs, which are convolutions of GPDs with the hard scattering part of the interaction. CFFs are the most basic quantities that one can unambiguously extract from the experimental data. The way of how CFFs enter the final amplitudes depends on the beam and target helicity states, which provides a welcome experimental filter to distinguish between many possible CFFs and justifies the need of measuring many observables. For brevity we skip the formulas showing how T DVCS and I depend on CFFs. They can be found in Ref. [56]. The evaluation of CFFs is discussed in Sec. 3.

Observables
Let us denote a four-fold differential cross section for a single photon production by d 4 σ b,c t (x Bj , t, Q 2 , φ), where t ∈ {←, →} and b ∈ {←, →} stand for the target and beam helicities, respectively, and c ∈ {+, −} stands for the beam charge. Here, x Bj = Q 2 /(2p · p − p) is the usual Bjorken variable. The cross sections can be used to construct many observables, like cross sections itself, but also differences of cross sections and asymmetries. For instance: Here, the capital letters in the subscripts of observables names denote beam and target polarizations, respectively, with U standing for "Unpolarized" and L standing for "Longitudinally polarized". We also analyze data for "Transversely polarized targets", which are distinguished by the subscript T . These data are provided for two moments, sin(φ − φ S ) and cos(φ − φ S ), which are distinguished by the corresponding labels in the superscripts, as for instance in A −,sin(φ−φ S ) (x Bj , t, Q 2 , φ). Furthermore there are observables probing only the beam charge dependency (subscript C), and those combining cross sections measured with various beam charges to drop either the DVCS or interference contribution (subscripts I and DVCS, respectively). A different group consists of Fourier-like observables related to specific modulations of the φ angle. For instance: Another observable used in this analysis is the slope b(x Bj , Q 2 ) of the t-distribution of the DVCS cross section integrated over φ. Within the LO formalism one can relate this observable to the transverse extension of partons in the proton. However, this requires the following assumptions, which are expected to hold at small x Bj : i ) dominance of the imaginary part of CFF related to GPD H, ii ) negligible skewness effect at ξ = x, iii ) exponential t dependence of the GPD H at fixed x. Since the DVCS t-distribution is usually not exactly exponential, in particular because GPDs for valence and sea quarks may have different t-dependencies, we evaluate b(x Bj , Q 2 ) by probing DVCS t-distribution in several equidistant points ranging from |t| = 0.1 GeV 2 to |t| = 0.5 GeV 2 and perform a linear regression on the logarithmized results. The chosen range of t is typical for the existing measurements of b(x Bj , Q 2 ) by COMPASS [57] and HERA experiments [58][59][60].

Imaginary part
At LO the imaginary part of a given CFF G ∈ {H, E, H, E}, is proportional to the combination of corresponding GPDs, G q ∈ {H q , E q , H q , E q }, probed at x = ξ: Here, the sum runs over all quark flavors (we remind that at LO gluons do not contribute to the DVCS amplitude), e q is the electric charge of the quark flavor q in units of the positron charge e and G q(+) is the singlet (C-even) combination of GPDs: for G ∈ {H, E} and: for G ∈ { H, E}.

Real part
At LO the real part of a given CFF G can be evaluated by probing the corresponding GPD G in two ways, that is, by integrating over one of two lines laying in the (x, ξ)-plane. This duality is a consequence of the polynomiality property required by the Lorentz invariance of GPDs, see Sec. 2. The first evaluating method is the "standard" one, where x values of the involved GPDs are probed at fixed ξ: Here, the quark propagators, 1/(ξ −x) and 1/(ξ +x), enter in a combination given by the type of probed GPDs. One has the difference (sum) for the unpolarized (polarized) GPDs. The second evaluating method is known as fixedt dispersion relation [38] and it involves the integral probing GPDs at ξ = x: Again, the combination of quark propagators depends here on the type of probed GPDs, exactly as for Eq. (37). The additional term in Eq. (38), C G (t), is the so-called subtraction constant. It has the same magnitude but the opposite sign for the CFFs H and E, and it vanishes for the CFFs H and E: After a quick examination of Eqs. (34) and (38), one may notice that the dispersion relation provides a welcome relationship between the real and imaginary parts of the same CFF. As a consequence however, at the LO approximation only GPDs in the limited case of x = ξ and the subtraction constant can be probed.

Subtraction constant
The subtraction constant introduced in Eq. (38) can be related to D-term form factor, D q (t), in the following way: Here, z = x/ξ and D q (z, t) is the D-term [61]. It was originally introduced to restore the polynomiality property in the first models based on double distributions [62], but later it has been recognized as an important element of the GPD phenomenology. Because the D-term vanishes outside the ERBL region |x| < |ξ|, it is not observed in the limit of ξ = 0, and it can be only studied in the "skewed" case of ξ = 0. By expanding the D-term in terms of Gegenbauer polynomials, one can obtain the following series: The first term of this expansion, d q 1 (t), is of a special importance, as it enters the quark part of the QCD energy momentum tensor and it provides an important information on how strong forces are distributed in the nucleon [63].
The subtraction constant can be evaluated by comparing Eqs. (37) and (38): which can be evaluated without principal value prescription, because the singularity at x = ξ is integrable in the expression: Unfortunately, naively setting ξ = 0 in the above formula results in a divergent integral. However, the following moments: are well defined for odd positive j and can be an- has a proper analytic behavior, as described in [64]. Such an analytical continuation can be written as: where we have introduced the analytic regularization technique [40,[64][65][66], given by the following prescription: Here, one subtracts as many terms from the Taylor expansion of f around zero as needed to make the integral convergent, and one treats the compensating terms to be convergent as well.
The analytic properties of the Mellin moments of GPDs has been never proved to be a consequence of general principles (neither it has never been proved to contradict general principles) and because of that can be only treated as a model assumption to be a posteriori confronted with experimental data 1 . We will make such an assumption, and calculate the subtraction constant as: The self-consistency of this approach will lead us to relations (62) and (63) among the otherwise nonrelated parameters of the fitting model.

Ansatz
We present in this study a global extraction of CFFs. According to the terminology used within the GPD community, "global" refers to constraining parameters of an assumed CFF functional forms from various measurements on a wide kinematic range. On the contrary in local extractions, CFFs are extracted as a set of disconnected values in bins of ξ and t (see for instance Ref. [34,35]). We restrict our analysis to the LO approximation and we neglect any contribution coming from higher-twist effects and kinematic (target mass and finite-t) corrections. We adopt the description of cross sections in terms of DVCS and BH amplitudes by Guichon and Vanderhaeghen, used for phenomenology for instance in Refs. [71,72] and publicly available in the open-source PARTONS framework [42]. We point out, that the limited phase space covered by available data, the precision of those data and the plenitude of involved dependencies force us to keep the parameterizations as simple as possible. Otherwise significant correlations appear between fitted parameters, which somehow obscures the interpretation of obtained results. The Ansatz introduced in this section is explicitly given for a factorization scale that one may recognize as the reference scale Q 2 0 at which the model is defined. To include the factorization scale dependence in our fit, that is for the comparison with experimental data of Q 2 = Q 2 0 , we consider the so-called forward evolution, i.e. the one followed by PDFs. The usage of the genuine GPD evolution equations requires the knowledge of GPDs in the full range of x independently on ξ, while in this analysis only the GPDs at x = ξ are considered. It was checked however with the GK GPD model [73][74][75], that the difference between the two evolution schemes is small

Decomposition into valence and sea contributions
In this work we use the decomposition scheme into valence and sea contributions inspired by the double distribution modeling of GPDs [62,76,77]. It gives us: for x = −ξ and G ∈ {H, E}, and for x = −ξ and G ∈ { H, E}. Here, G q val and G qsea are GPDs for valence and sea quarks, respectively. With this decomposition one can replace Eqs. (35) and (36) by one equivalent expression: For the GPDs H q and H q at ξ = 0 we use an Ansatz that is commonly used in phenomenological analyses of GPDs: Here, pdf q G (x) is either a parameterization of the unpolarized PDF, q(x), for the GPD H q or a parameterization of the polarized PDF, ∆q(x), for the GPD H q . The profile function, f q G (x), fixes the interplay between the x and t variables, and it is given by: where A q G , B q G and C q G are free parameters to be constrained by experimental data. This form of f q G (x) allows to modify the classical A q G log(1/x) term by B q G (1−x) 2 in the small x region and by C q G (1−x)x in the high x region. The terms proportional to B q G and C q G were found to work best in the analysis of EFF data (see Sec. 6 analytic at x = 0. To keep the distance between the active quark and the spectator system finite, see Eq. (23), we require to have C q val H = −A q val H . The profile function given by Eq. (55) is more flexible than that used in the GK model [73][74][75], and in the VGG model [76,[78][79][80] . In particular, it should be flexible enough to take into account a different interplay between the x and t variables in the valence and sea regions if required by experimental data. We note that f q G, [47,81] to fit EFF data can not be used in this analysis because of the aforementioned issue with the analyticity caused by For GPDs H q and H q at ξ = x we utilize the concept of skewness function: In our case: where G q (x, 0, t) is given by Eq. (54). We assume the following form of the skewness function: where a q G is a free parameter to be constrained by experimental data. Two other parameters, b q G and c q G , which govern the t-dependence of the skewness function, are fixed in a way to avoid singularities in the evaluation of the subtraction constant. Namely, to use the analytic regularization prescription at fixed t one has: where a and f (x) were introduced in Eq. (48) and δ describes the behavior of PDFs at x → 0: The singularities appear in the two first compensating terms at a = 0 and 1 − a = 0, that is for t ≡ t ∞ 0 = −δ/A q G and t ≡ t ∞ 1 = (1 − δ)/A q G , respectively. The problem does not emerge for higher compensating terms as one typically has 0 < δ < 1 for valence quarks and 1 < δ < 2 for sea quarks. The singularities are regularized by requiring f (0) = 0 at t ∞ 0 and f (0) = 0 at t ∞ 1 , which is achieved by setting b q G and c q G to: where p 0 , p 1 and α are parameters of PDF parameterizations introduced in Sec. 5. We stress that our Ansatz for the skewness function is explicitly defined at x = ξ and it can not be generalized to the case of x = ξ without a nontrivial modification. The form of the skewness function has been selected because of the following reasons: i ) for sufficiently small x and t, the skewness function coincides with a constant value given by a q G . Such a behavior was predicted for HERA kinematics [82] and it was used in one of the first extractions of GPD information [83] from H1 data [84]. These data suggest a q H ≈ 1. ii ) In the limit of x → 1 the skewness function is driven by 1/(1 − x 2 ) 2 . This form has been deduced from Ref. [39], where the power behavior of GPDs in the limit of x → 1 was studied within the pQCD approach. iii ) The subdominant t-dependence in Eq. (58) has been inspired by the skewness function evaluated from GK [73][74][75] and VGG [76,[78][79][80] GPD models, both being based on the one component double distribution modeling scheme [62]. In those models the t-dependence of the skewness function is dominated by b q G +c q G log(1 + x) term. In our Ansatz we multiply it by (1 − x) to avoid any t-dependence at x → 1, which is imposed by Ref. [39], where it was shown that GPDs should not depend on t in this limit, regardless the value of ξ.

CFFs E and E
For E and E we use a simplified treatment justified by the poor sensitivity of the existing measurements on the corresponding CFFs. Moreover, the forward limit of those GPDs has not been measured, resulting in a need of fixing more parameters than for the GPDs H and H, if an analog modeling was to be adopted.
The modeling of the CFF E is similar to that of H and H, i.e. it is based on the dispersion relation with the subtraction constant of the opposite sign as for H. We only consider the valence sector with the forward limit of GPD E q val taken from Ref. [81]: where κ u = 1.67 and κ d = −2.03 are the magnetic anomalous moments for up and down quarks, respectively, and where β u val = 4.65, β d val = 5.25, γ u val = 4 and γ d val = 0. The parameter α u val = α d val ≡ α is fitted to EFF data as described in Sec. 6. The normalization parameter: ensures that: At ξ = 0 the t-dependence is introduced analogously as for the GPDs H and H, i.e.: with f q val E (x) given by Eq. (55), where: are free parameters being fitted to EFF data as described in Sec. 6. It is assumed that A q val E = A q val H . The skewness function used to evaluate E q val (x, x, t) from E q val (x, 0, t) is of the following form: which again is deduced from Ref. [39]. Here, a u val E = a d val E ≡ a q val E = 1 and any t-dependence is neglected. An additional ξ-dependence coming from twist-three and twist-four distribution amplitudes of the nucleon is denoted by f (x). In the present analysis, where only the leading twist is considered, it is assumed For the GPD E the shape of the corresponding form factor is fixed by the GK model [73][74][75]. Only the normalization parameter, N E , is fitted to the experimental data,

Inequalities
We impose two extra constraints implied by the positivity of parton densities in the impact parameter space. Namely, for the profile functions we require to have: being imposed by Eq. (24). The requirements are implemented in a way to penalize combinations of fitted parameters that do not hold the inequalities. It is achieved by introducing a penalty to the χ 2 function, which value is proportional to the maximum violation of a given inequality. Such an implementation allows us to keep the χ 2 function smooth. We report that the inequalities given by Eqs. (70) and (71) have proved to be important to constrain parameterizations for the sea contribution to GPD H and for the valence contribution to GPD H. It is due to a low sensitivity to those contributions and a limited phase space covered by the available data. The inequality (25) is not checked during the minimization. Namely, by following Ref. [47] one can obtain: where for f q H (x) = f q H (x) (condition not imposed in this analysis) one can minimize right hand side over b 2 ⊥ , and get the strongest inequality: which is used for instance in Ref. [81]. To effectively use the inequality (72) in the minimization, one should simultaneously fit all free parameters of f q H (x), f q H (x), f q E (x) and e(x), which presently is not the case, as we fit EFF and DVCS data separately. In addition, for the GPDs H and E we consider only the valence contribution, while the sea sector may have a significant impact on (72) in the range of small x Bj . The validity of (72) has however been checked a posteriori, that is after the χ 2 -minimization. The result of this test shows that the inequality is violated for most of the replicas for very small b 2 ⊥ . To correct this one should change the simplified Ansatz used for GPD E and preferably fit its free parameters simultaneously with those for other GPDs, which is beyond the scope of this analysis.

Approximations and summary
For the GPD H and valence quarks, all parameters of the profile function given by Eq. (55) Sea components are not yet fully symmetric in our fit, as we do not impose the flavor symmetry in PDFs, see Sec. 5. Due to the lack of precision of axial EFF data, for the GPD H all parameters for both profile and skewness functions are fitted to DVCS data. As the contribution coming from ∆q sea is subdominant, we neglect it entirely. For valence quarks, similarly to the GPD H, we allow the profile function to be different for u val and d val , while for the skewness function one has H . For E and E only the valence quarks are considered. For the GPD E all free parameters for both profile function and forward limit given by Eq. (64) are fitted to EFF data. For the GPD E the normalization parameter N E is fitted to DVCS data. In total, we fit 9 parameters to EFF data (see Tab. 2) and 13 parameters are constrained by DVCS data (see Tab. 5).

Analysis of PDFs
The analytic regularization prescription introduced in Eq. (48) requires the function q(x)/x −δ to be nonzero and analytic at x = 0. However, the numeric evaluation of such a function and its derivatives near x = 0 is difficult for typical PDF sets, like those published by NNPDF [85] and CTEQ [86] groups, because of several numerical issues. Namely, interpolations in PDF grids and extrapolations outside those grids for small x, and a delicate cancellation of the numerator and denominator of q(x)/x −δ make the evaluation numerically unstable. The problem can be avoided with functional parameterizations of PDFs, which can be used for a straightforward evaluation of q(x)/x −δ and its derivatives. In addition, such parameterizations allow for a significant reduction of computation time and a precise determination of δ.
In this analysis the parameterization used for both unpolarized and polarized PDFs reads: where: describes the evolution in the renormalization scale, recognized as Q 2 , where Q 2 0 = 2 GeV 2 can be identified as the initial scale. The parameterization contains thirteen parameters, constrained for each quark flavor in a fit to "NNPDF30 lo as 0118 nf 3" set [85] for unpolarized PDFs and to "NNPDFpol11 100" set [87] for polarized PDFs. These fits are performed for a grid of nearly 10000 points equidistantly distributed in the ranges of 10 −4 < x < 0.9 and 1 GeV 2 < Q 2 < 20 GeV 2 . We have selected the sets by NNPDF group as: i ) they are provided for both q(x) and ∆q(x), ii ) for q(x) they are provided at LO and three active flavors, iii ) they provide reliable estimation of uncertainties extracted through the neural network approach and the replica method. The PDF sets are handled with LHAPDF interface [88].
A consistent fit to all replicas of used NNPDF sets allows us to reproduce the original uncertainties of PDFs. Figure 3 demonstrates the agreement between the original sets and our fits. The comparison is satisfactory and in particular the central values of fitted PDFs stay within the original uncertainty bands.

Analysis of Dirac and Pauli Form Factors
Equations (11)- (14) link GPDs to EFFs, thus they can be used to constrain GPDs from elastic data. In this analysis we use this feature to constrain the interplay between the x and t variables for the GPDs H and E, however only for the valence quarks. More precisely, we fix the parameters of f q val H/E (x) defined in Eq. (55) by using the Ansatz for H q val (x, 0, t) given in Eq. (54) and that for E q val (x, 0, t) given in Eq. (67). We achieve this by studying proton, F p i (t), and neutron, F n i (t), Dirac and Pauli form factors: related to Sachs form factors as follows: Sachs form factors can be used to define many observables, in particular: • the magnetic form factor normalized to both magnetic moment, µ p or µ n , and dipole form factor: where i = p, n and with M 2 D = 0.71 GeV 2 . • the ratio of electric and normalized magnetic form factors where i = p, n. • the squared charge radius of neutron: Experimental data for those observables and for G n E (t) come from the sources summarized in Table  1. For the selection of observables and related data we follow Ref. [81], which deals with the subject in great detail. In total, the used sample of data consists of 178 data points covering the range of 0.017 GeV 2 ≤ |t| ≤ 31.2 GeV 2 . Our fit to those data ends for the central PDF replica with χ 2 /ndf = 129.6/(178 − 9) ≈ 0.77. The fit is repeated for all other PDF replicas and for 100 instances of replicated data to propagate the uncertainties of EFF data to the analysis of CFFs. A single instance of replicated data is generated with the following prescription: where v i is the measured value associated to the experimental point i. The total uncertainty, which is also used to evaluate the χ 2 value in the fit to EFF data, reads: where ∆ stat i and ∆ sys i are statistical and systematic uncertainties, respectively, both linked to the point i. The generator of random numbers following a specified normal distribution, f (x|µ, σ), is denoted by rnd j (µ, σ), where j is both the identifier of a given replica and a unique random seed. The final result consists of 201 replicas, where each replica represents a possible realization of fitted EFFs. Such a set of replicas can be used to estimate mean values and uncertainties of the fitted parameters, which we summarize in Table 2. Experimental data are superimposed with the results of our fit in Fig. 4, while distributions of quark EFFs are shown in Fig. 5.  Table 3 summarizes DVCS data used in this analysis. Currently, only proton data are used, while those sparse ones for neutron targets are foreseen to be included in the future. We note, that recent Hall A data [117,118] published for unpolarized cross sections, d 4 σ − U U , are not used in the present analysis since it was not possible to correctly describe them with our fitting Ansatz, and their inclusion was causing a bias on the final extraction of GPD information. This specific question will be addressed in a future study. However we keep the recent Hall A data for differences of cross sections, ∆d 4 σ − LU . In addition, we skip all DVCS data published by HERA experiments since they should be dominated by gluons, and would probably lead to misleading conclusions in our analysis which is adapted to the quark sector. We will back to these problems in Sec. 9.

DVCS data sets
We apply two kinematic cuts on experimental data: The purpose of those cuts is to restrict the phasespace covered by experimental data to the deeply virtual region where one can rely on the factorization between GPDs and the hard scattering kernel. The values of those cuts have been selected with the correspondence to KM analysis [36], which may help to compare results of both analyses in the future. Unpol. PDF unc. Unpol. PDF unc.
Unpol. PDF unc. Unpol. PDF unc.  Table 1 and fitted according to the text. The gray bands indicate 68% confidence level for uncertainties coming from EFF data. The corresponding bands for unpolarized PDFs are indicated by the labels. For data points provided with systematic uncertainties, the inner bars represent statistical uncertainties, while the outer ones are for the quadratic sum of statistical and systematic uncertainties. Otherwise, statistical uncertainties are shown.

Uncertainties
In this analysis all the uncertainties are evaluated with the replica method, which for a sufficiently large set of replicas allows one to accurately reproduce the probability distribution of a given problem. We distinguish four types of uncertainties on the extracted parameterizations of CFFs. They origin from: i) DVCS data, ii ) unpolarized PDFs, iii ) polarized PDFs) and iv ) EFFs. Each type of uncertainty is estimated independently, as described in the following.
The uncertainties coming from DVCS data are estimated with a set of 100 instances of replicated data. Similarly to the analysis of EFFs, see Sec. 6, a single instance of replicated data is generated with the following prescription: Here, v i is the measured value associated to the experimental point i, which comes with statistical, ∆ stat i , systematic, ∆ sys i , and normalization, ∆ norm i , uncertainties. The latter one appears whenever the observable is sensitive to either beam or target polarization, or to both of them. In such cases the polarization describes the analyzing power for the extraction of those observable and the normalization uncertainties are related to the measurement or other determination of the involved polarizations. The total uncertainty, which is also used in the fit of CFFs to evaluate the χ 2 value, is evaluated according to Eq. (84).
The uncertainties coming from unpolarized and polarized PDFs are estimated by propagating our  parameterizations of replicas by NNPDF group, see Sec. 5. Namely, we repeat our fit to DVCS data separately for each PDF replica, that is 100 times for unpolarized PDFs and 100 times for polarized PDFs. A similar method is used to evaluate the uncertainties coming from EFF parameterizations. We repeat our fit to DVCS data for each replica obtained in the analysis of EFF data, see Sec. 6.
As a result of this analysis we obtain a set of 401 replicas, where each of them represents a possible realization of CFF parameterizations. For a given kinematic point the mean and uncertainties can be then estimated by taking the mean and the standard deviation of values evaluated from those replicas.

Performance
For the central PDF and EFF replicas the minimum value of the χ 2 function returned by the minimization routine (Minuit [130] ported to ROOT [131]) is 2346.3 for 2600 experimental points and 13 free parameters, which gives the reduced value equal to 2346.3/(2600 − 13) ≈ 0.91. The values of χ 2 function per experimental data set are summarized in Table  4. In general the agreement between our fit and experimental data is quite good. The only exception is for the COMPASS point, on which we will comment in the following.
In this analysis a set of 401 replicas is obtained that can be used to estimate mean values and uncer-tainties of the fitted parameters. We summarize this in Table 5, where we also indicate ranges in which the minimization routine was allowed to vary the fitted parameters. As one can see from this table, the parameter B qsea H is found to be at the lower limit of the allowed range. Without this limit the fit would end with a much smaller value of B qsea H compensated by a higher value of A qsea H , which we consider to be unphysical. The problem with B qsea H is a consequence of the sparse data covering low-x Bj region, but it may be also a sign of the breaking of the LO description in the range dominated by sea quarks and gluons.
In addition to the fitted parameters, in Table 5 we also show typical values of b q G and c q G evaluated from Eqs. (62) and (63), respectively. These values indicate that the t-dependence in the skewness function is subdominant, which we consider to be an expected feature. Figures 6-9 provide a straightforward comparison between the results of our analysis and some selected data sets coming from various experiments. These plots can be used in particular to estimate the effect of PDF and EFF uncertainties, which one can not judge from Table 4.
In Fig. 7 the comparison for both d 4 σ − U U and ∆d 4 σ − LU data coming from Hall A is shown, however one should keep in mind that only those for ∆d 4 σ − LU are used in this analysis. In general the agreement between the unpolarized cross section data published by Hall A [117,118] and our fit is poor, which can be qualitatively expressed be the reduced value of χ 2 function evaluated for those data equal to 5144.4/594 ≈ 8.66. The agreement is better for φ 0, where for 0 < φ < 45 • and 315 • < φ < 360 • one has 232.4/132 ≈ 1.76. We note that the KM model also has a problem to describe Hall A unpolarized cross section data at LT and LO accuracy [36], and it was suggested [118] that including HT and NLO corrections may be needed to improve the quality of the fit. Figure 9 demonstrates the comparison between results of our fit and experimental data for the tslope b coming from the measurements by COM-PASS [57], ZEUS [59] and H1 [58,60] experiments. We remind that for this observable only the COM-PASS point is used in this analysis. As one can judge from the figure, the agreement between our fit and experimental data is lost for small x Bj , i.e. for the range where sea quarks and gluons dominate. We only leave the COMPASS point to have a coverage in the intermediate range of x Bj , where valence quarks may still contribute significantly. While it is tempting to improve the agreement for small x Bj by adding extra terms to f qsea H (x), like those proportional to (1 − x) n with large n, we refrain from doing that, treating the disagreement as a possible manifestation of NLO effects. The included comparison with GK [73][74][75] and VGG [76,[78][79][80]] GPD models shows that GK manages to reproduce the trend dictated by the HERA data reasonably well. This is not a full surprise, as this model was constrained in the low-x Bj region, however by Deeply Virtual Meson Production (DVMP) data. It is worth pointing out, that the lowest possible order of contribution to DVMP starts with α 1 S , while for DVCS one has α 0 S . It is beyond the scope of this paper, but in the future multichannel analyses of exclusive processes the issue of using the same order of pQCD calculations may become important.

Subtraction constant
The subtraction constant obtained in this analysis is shown as a function of t at Q 2 = 2 GeV 2 and as a function of Q 2 at t = 0 in Fig. 10. One may observe a large uncertainty coming from PDF param- Table 5: Values of the parameters fitted to DVCS data together with estimated uncertainties coming from those data, (un-)polarized PDFs and EFFs. Two last columns indicate the limits in which the minimization routine was allowed to vary the corresponding parameters. In addition, exemplary values of b q G and c q G parameters evaluated at Q 2 = 2 GeV 2 from Eqs. (62) and (63)

Compton Form Factors
The CFFs obtained in this analysis are shown as a function of ξ for the exemplary kinematics of t = −0.3 GeV 2 and Q 2 = 2 GeV 2 in Figs. 11-14. In complement of our results, we also show the curves evaluated at the same kinematics for the GK [73][74][75] and VGG [76,[78][79][80] models. In general, VGG is closer to the results of our fit, which is expected, as this model was primarily designed to reproduce DVCS data in the valence region. We conclude that both models overestimate the imaginary part of the CFF H, which for VGG is compatible with the conclusions of Ref. [34], where local fits to recent JLab DVCS data are reported.

Nucleon tomography
Exemplary distributions of u(x, b ⊥ ) and ∆u val (x, b ⊥ ) observables corresponding to Eqs. (17) and (18), respectively, are shown in Fig. 15. These plots illustrate how parton densities (longitudinal polarization of those partons) are distributed inside the unpolarized (longitudinally polarized) nucleon, however without the possibility of showing the corresponding uncertainties.
This difficulty can be overcome by showing the normalized second moments of q(x, b ⊥ ) and ∆q val (x, b ⊥ ) distributions, see Eqs. (21) and (22) and Figs. 16 and 17, respectively. One may note that b 2 ⊥ q (x) → 0 when x → 1, which is expected, as in this limit the position of the active quark is equivalent to the center of the coordinate system in which the impact parameter is defined. The corresponding distributions illustrating the mean squared distance between the active quark and the spectator system, see Eq. (23), are shown in Fig. 18. In our studies d 2 ⊥ q is finite when x → 1, which we consider to be a welcome feature imposed by our Ansatz. We   Fig. 6: Comparison between the results of this analysis, some selected GPD models and experimental data published by CLAS in Refs. [128,129] for d 4 σ − U U at x Bj = 0.244, t = −0.15 GeV 2 and Q 2 = 1.79 GeV 2 (left) and for A − U L at x Bj = 0.2569, t = −0.23 GeV 2 , Q 2 = 2.019 GeV 2 (right). The solid curves and the gray bands surrounding those curves correspond to the results of this analysis and 68% confidence levels for the uncertainties coming from DVCS data, respectively. The magnitudes of the additional uncertainties to the plotted observables and coming from unpolarized PDFs, polarized PDFs and EFFs, can be separately estimated from the bands located below the curves and labeled with "PDF unc.", "Pol. PDF unc." and "EFF unc.", respectively. Each of those four uncertainties is evaluated in the analysis of the corresponding set of replicas. The inner bars on data points are for statistical uncertainties, while those outer ones are for the quadratic sum of statistical and systematic uncertainties. The dotted curve is for the GK GPD model [73][74][75], while the dashed one is for VGG [76,[78][79][80]    : Comparison between the results of this analysis, selected GPD models and experimental data published by COMPASS in Ref. [57] for the slope b at Q 2 = 1.8 GeV 2 . Except COMPASS data also those obtained in ZEUS [59] and H1 [58,60] experiments are shown. Note that the shown data differ by Q 2 values indicated in the legend. For further description see Fig. 6.
also study a possible difference between b 2 ⊥ u val (x) and b 2 ⊥ d val (x), see Fig. 19. This plot suggests that the distribution of u val (x, b ⊥ ) is narrower (broader) than d val (x, b ⊥ ) in the region of high (low) x. A firm conclusion however is not possible at this moment because of the large uncertainties.

Summary
In this paper we proposed new parameterizations for the border and skewness functions. Together with the assumption about the analyticity properties of the Mellin moments of GPDs, those two ingredients allowed the evaluation of DVCS CFFs with the LO and LT accuracy. The evaluation was done with the dispersion relation technique and it included a deter- mination of the DVCS subtraction constant, which is related to the QCD energy-momentum tensor.
In order to build and constrain our parameterizations we utilized many basic properties of GPDs, like their relation to PDFs and EFFs, the positivity bounds, the power behavior in the limit of x → 1 and even the polynomiality property allowing the evaluation of the subtraction constant by comparing two equivalent ways of CFF computation. Our parameterizations provide a genuine access to GPDs at (x, 0, t) kinematics and therefore they can be used for nucleon tomography.
We performed the analysis of PDFs and obtained a set of functional parameterizations allowing the reproduction of original values and uncertainties. A small number of free parameters appearing in our approach was constrained by EFF and DVCS data. We considered all proton DVCS data, however not all of them entered the final analysis because of the used kinematic cuts and the initial data exploration. Our work was done within the PARTONS project that provides a modern platform for the study of GPDs and related topics. The quality of our fits is quite good. The fit to EFF data returns χ 2 = 129.6 for 178 data points and 9 free parameters, while that to DVCS data returns χ 2 = 2346.3 for 2600 data points and 13 free parameters. The good performance proves that our parameterizations, including the assumed analyticity property, are not contradicted by the used experimental data. We consider the unsurprising discrepancy between our results and the HERA data as a possible manifestation of gluon effects, and the surprising discrepancy with unpolarized cross sections published by Hall A as a possible need for higher twist or kinematic corrections.
Our analysis also favors small values of the subtraction constant. We extract the distributions of q(x, b ⊥ ) and ∆q val (x, b ⊥ ) and we study the properties of those distributions. This analysis is characterized by a careful propagation of uncertainties coming from PDFs parameterizations, EFF and DVCS data, which we achieved by using the replica method. The first successful step towards the reduction of model uncertainties has been done already by selecting PDFs parame-EFF unc.
Unpol. PDF unc. terizations based on the neural network technique. We plan to extend the usage of neural networks to other sectors of our future analyses for further reduction of model uncertainties.