Generalized parton distributions through universal moment parameterization: non-zero skewness case

We present the first global analysis of generalized parton distributions (GPDs) combing lattice quantum chromodynamics (QCD) calculations and experiment measurements including global parton distribution functions (PDFs), form factors (FFs) and deeply virtual Compton scattering (DVCS) measurements. Following the previous work where we parameterize GPDs in terms of their moments, we extend the framework to allow for the global analysis at non-zero skewness. Together with the constraints at zero skewness, we fit GPDs to global DVCS measurements from both the recent JLab and the earlier Hadron-Electron Ring Accelerator (HERA) experiments with two active quark flavors and leading order QCD evolution. With certain choices of empirical constraints, both sea and valence quark distributions are extracted with the combined inputs, and we present the quark distributions in the proton correspondingly. We also discuss how to extend the framework to accommodate more off-forward constraints beyond the small $\xi$ expansion, especially the lattice calculated GPDs.


Introduction
Over the past two decades, there has been growing interest in the higher-dimensional structures of the nucleon in the studies of the non-perturbative strong interactions described by quantum chromodynamics (QCD). Consequently, generalized parton distributions (GPDs) [1][2][3] that unify the elastic form factors (FFs) and Parton Distributions Functions (PDFs) into a single set of 3-dimensional (3D) functions and contain important information about the nucleon mass, angular momentum and mechanical properties [2,4,5] have gained increasing attention. While it is believed that GPDs can be experimentally accessed by exclusive productions of particles off nucleons such as deeply virtual Compton scattering (DVCS) [6] and deeply virtual meson production (DVMP) [7,8], obtaining them and the corresponding nucleon 3D structures from experimental data has remained a dream. On the one hand, persisting efforts are needed for an adequate amount of data as required for the extraction of such high-dimensional quantities. On the other hand, the extraction of the critical information from experimental data has also been challenging.
Several recent breakthroughs have made this possible. Large data sets of exclusive measurements are generated at Jefferson Lab (JLab) [9][10][11][12], and more will come from the planned Electron-Ion Collider (EIC) [13,14] at Brookhaven National Lab (BNL). Additionally, novel developments in lattice QCD allow the access to the nucleon structures from first principle calculation [15][16][17][18], providing information almost impossible to get with current and future experiments. This recent progress, together with the perennial efforts on the global extraction of PDFs [19][20][21], FFs [22] as well as lattice calculation of generalized form factors [23][24][25][26], have pushed the study of nucleon 3D structures to a new stage, which, on the other hand, requires extra works to put all the above inputs together through global analysis to obtain the state-of-art nucleon 3D structures with GPDs.
Therefore, we are in need of a global analysis program of GPDs with both experimental data and lattice computation. Lots of efforts are put into parameterizing and accessing GPDs from various inputs in the literature [27][28][29][30][31][32][33][34][35]. However, a framework that could combine experimental data and lattice computation still seems to be lacking. In the previous work [36], we proposed the GPDs through universal moment parameterization (GUMP) program for such a purpose and studied the zero-skewness case as an illustrative example. In this work, we follow the previous setup and extend it to non-zero skewness, which allows us to also include the exclusive measurements in the global analysis. With such a program, we perform the global analysis of GPDs that includes global PDFs [19][20][21], FFs [22] and DVCS measurements from both JLab [9][10][11][12] and previous H1 [37] experiments at Hadron-Electron Ring Accelerator (HERA) along with the relevant lattice calculations related to GPDs [17,23] for the first time.
The organization of the paper is as follows. In sec. 2, we introduce the setup for the GUMP program and discuss how to build the 3D GPDs with a handleable set of parameters. In sec. 3, we discuss more details of the fit and present the extracted Compton form factors (CFFs) and GPDs, where we also discuss how to extend the present framework to accommodate more future inputs. In the end, we conclude in sec. 4.

Conventions and fitting procedure
In the previous work [36], we discussed how GPDs can be parameterized in terms of their moments based on the technique that has been systematically developed in ref. [30] and used for instance in ref. [31]. In this section, we will discuss how we apply such a framework to the GPDs at non-zero skewness, including GPDs of various different species and flavors.

Convention and decomposition of GPD
We start by considering the four leading-twist GPDs H, E, H and E, where each of them has different flavors u, d and g corresponding to the up and down quarks as well as the gluon. Although the strange quark might have sizable effects that should be considered for a more careful treatment, the flavor separation of GPDs is typically challenging due to the lack of a flavor-sensitive probe [38]. Strange quark distributions will not be well constrained with just DVCS measurements, and thus we will leave this for the future work with more flavor-sensitive measurements.
For simplicity, the four different species of GPDs will be treated almost equally in this work. Except that the vector GPDs H and E satisfy slightly different scale evolution equations from the axial-vector ones H and E [39], the same parameterization will be used for all four GPDs. Although GPDs or combinations of GPDs of different species H, E, H and E correspond to different helicity amplitudes and might have different behaviors accordingly, knowledge of different GPDs species is limited and thus implementing a different parameterization for them individually may not lead to any significant difference. Therefore, we collectively denote all GPDs as F q,g (x, ξ, t) where F = {H, E, H, E} and q = {u, d}, with x the parton momentum fraction, ξ the skewness parameter and t the total momentum transfer squared and parameterize them in the same manner. We note that the skewness parameter ξ generally takes the value from −1 to 1, but it will be considered to be positive hereafter, as GPDs are commonly defined to be symmetric in ξ, subject to their parity symmetry.
An intriguing feature of the GPDs is that they consist of two regions with totally different physical interpretations. In the PDF-like region where x > ξ or x < −ξ, GPDs resemble the PDFs which are interpreted as the amplitudes of emitting and reabsorbing a parton (quark, antiquark or gluon). On the other hand, in the distribution amplitude (DA)-like region where −ξ < x < ξ, they resemble the DAs instead and are interpreted as the amplitudes of emitting/absorbing a parton-antiparton pair. Therefore, GPDs do not naturally have distinguishable quark and antiquark components, especially in the DA-like region. This feature indicates that the parameterization of GPDs should be flexible in both the PDF and DA-like regions, since different physics are involved.
Motivated by such property, one can write GPDs as, where the ∓ depends on the parity of the GPDs which takes − for vector GPDs H and E and + for axial vector GPDs H and E. Here the subscriptq stands for the quark distributions excluding antiquark, to be distinguished from the subscript q for the quark flavor. The three terms above describe the amplitudes corresponding to the quark, antiquark, and quark-antiquark pair respectively. Accordingly, Fq(x, ξ, t) and Fq(x, ξ, t) have support x > −ξ, while the F qq (x, ξ, t) has support ξ > x > −ξ. In the semi-forward limit where ξ = 0 and t generally non-zero, Fq(x, ξ, t) and Fq(x, ξ, t) reduce to the quark and antiquark t-dependent PDFs with support x > 0, whereas F qq (x, ξ, t) vanishes with its vanishing support. We note that GPDs can always be decomposed into such three parts that each term represents the behavior of GPDs in different regions respectively. The extra DA terms F qq (x, ξ, t) allow extra flexibility for the parameterization of GPDs particularly in the DA-like region, as we will discuss with more details in the next section. While the quark and antiquark GPDs Fq(x, ξ, t) and Fq(x, ξ, t) are the natural generalization of the quark and antiquark PDFs, the DA terms F qq (x, ξ, t), also known as the D-terms in some other contexts [40], only exists in the off-forward case. Such terms live in the DA-like region only, and they will be the dominant effects in the large ξ limit ξ → 1 or the asymptotic limit where renormalization scale µ → ∞ [41,42]. Since the DA terms vanish at the crossover line x = ±ξ, they can always be expressed with a set of complete polynomials like Gegenbauer polynomials, which is exactly the conformal moment expansion. Therefore, each of the three terms can be expressed in terms of their conformal moments and satisfies the polynomiality condition.
In this work, we will focus on the quark and antiquark GPDs Fq(x, ξ, t) and Fq(x, ξ, t), whereas the extra DA terms F qq (x, ξ, t) will not be put into the global analysis for two reasons. First, as discussed in the previous work, we will consider the high energy limit where ξ → 0 ideally and the DA-like region disappears in such a limit. Second, since the DA terms vanish at the crossover line x = ±ξ, they are hard to extract with measurements of DVCS or similar processes which provide effectively just CFFs -the DA terms contribute to the real part of the CFFs only, corresponding to the subtraction terms in the dispersion relations of CFFs [31]. Consequently, the behaviors of GPDs in the DA-like region given by the quark and antiquark GPDs Fq(x, ξ, t) and Fq(−x, ξ, t) will be ambiguous, since one can add extra DA terms without affecting the CFFs too much in the small ξ limit. This will be discussed with more details in the next section. We note that another notation of quark and antiquark GPDs has been used in the literature, see for instance ref. [43] Fq(x > ξ, ξ, t) = F q (x, ξ, t) , with the same convention for the ∓ sign. Since the quark and antiquark GPDs are defined with the same function F q (x, ξ, t) partitioned into two parts, the full quark GPD F q (x, ξ, t) cannot be written as the sum of them F q (x, ξ, t) = Fq(x, ξ, t) ∓ Fq(−x, ξ, t). Due to the loss of linearity in such a GPD decomposition, which will also break the polynomiality condition for quark and antiquark GPDs respectively, this definition will not be used in this work. 1 To summarize, we will consider quark and antiquark GPDs Fq(x, ξ, t), Fq(x, ξ, t) as well as the gluon GPDs F g (x, ξ, t) in this work with four twist-two GPDs H, E, H and E, forming a basis with 20 GPDs. Conventionally, we also rewrite the basis by defining the valence combination: F q V (x, ξ, t) ≡ Fq(x, ξ, t)−Fq(x, ξ, t), and thus the basis can be written with F q V (x, ξ, t), Fq(x, ξ, t) and F g (x, ξ, t) equivalently, analogous to the basis commonly chosen for PDFs global analysis in the literature [19][20][21].

Parameterization of moments
After introducing the GPDs necessary for consideration, we briefly discuss the parameterization method for GPDs. We have introduced the GUMP parameterization in the previous work [36] and more about the conformal moment representation can be found in refs. [30,31]. Generally, we express all GPDs F (x, ξ, t) in terms of their conformal moments F j (ξ, t) in the form of, with p j (x, ξ) the known conformal wave functions. Therefore, all GPDs F (x, ξ, t) are equivalently given by F j (ξ, t). With the polynomiality condition of GPDs [3], the moments can be expressed as polynomials of ξ of given order:

4)
and then these F j,k (t) can be used to construct the GPD F (x, ξ, t).
In the high energy limit, the skewness parameter is small where the first few terms dominate, and we have [31] where we implicitly assume that the F j,k (t) are only non-zero for j > k according to the polynomiality condition. Practically the moments F j,k (t) can be written in terms of the shape governed by the Euler beta function B, the Regge trajectory for the t-dependence and the extra residual term β(t): with which the GPDs are parameterized in terms of the free parameters in eq. (2.6). In the simplest case, only one set of ansatz is needed, so we set i max to be just 1.
In the previous work, we simply set the residual function β(t) to be 1 for the extraction of valence quark distributions in the zero-skewness case, since their algebraic decaying behaviors in |t| can be parameterized well by the Regge trajectory. However, for the sea distributions, it has been observed that in high energy processes, the different crosssections drop exponentially as the momentum transfer |t| increases e.g., for DVCS [37,44], which differs from the power-law behavior indicated by the Regge theory. Therefore, we incorporate the extra exponential behavior into the β(t) to set β(t) = exp(bt). Two β(t)s of different b slope are embedded in the parameterization of the vector GPDs Hq ,g and the axial vector GPDs Hq ,g respectively, with q = {u, d} for sea quarks. We note that due to the lack of flavor-sensitive probe, we assume the same b slope for the gluon and sea quarks of different flavor and also the same Regge slope α for them which is fixed to be α q,g = 0.15 from the pomeron trajectory [45].
Even with such simplification, the parameters are still too many to be fully determined from measurements. The lack of off-forward constraints makes it almost impossible to get the shape of GPDs at non-zero skewness, which is known as the deconvolution problem of GPDs [46], stating that the shape of GPDs can not be uniquely determined from CFFs. Therefore, extra empirical constraints are still needed for the extraction of GPDs. The simplest choice for the ξ-dependent terms is to set them to be proportional to the leading terms, as has been done in ref. [31]: with R k the ratio between them. In this work, we have just one parameter R 2 for the ξ 2

GPDs species and flavors Fully parameterized GPDs linked to
Proportional constants  terms for each GPD F (x, ξ, t) that accounts for its extra ξ-dependence.
Another difficulty in GPD parameterization is about the flavors and species of GPDs. While the extraction of one 3D function is already challenging, the simultaneous extraction of GPDs with multiple flavors and species is yet more difficult. Two of the four different species of GPDs, E and E, do not have corresponding forward PDFs, unlike H and H which reduce to the PDF f (x) and helicity PDF ∆f (x) respectively, adding extra difficulties to their extraction. Consequently, we have to reduce the number of parameters associated to these two GPDs due to the lack of constraints. To do so, we set the two valence quark distributions to be proportional (E u V ∝ E d V and E u V ∝ E d V ) and the sea quark as well as the gluon distributions to be proportional to the corresponding H and H GPDs.
We note that these extra empirical constraints, both for the ξ-dependent terms and for the E/ E GPDs, are added for purely practical purposes. As mentioned above, they simply represent the lack of information of GPDs in the off-forward region, which can and should be improved in the future with more inputs from both lattice calculations and experimental measurements. In table 1, we collect whether GPDs of different species and flavors are either parameterized independently with eq. (2.6) or linked to the other GPDs with some empirical assumptions in this work.
Finally, we comment on the comparison of the GUMP program to the other GPD global analysis programs with conform moments framework, specifically the Kumerički-Müller (KM) model and its extensions [31,39,47,48]. As clearly stated in the previous work [36], the GUMP program follows the general conformal moment parameterization of GPDs, and adopts many useful observations there to make the GUMP parameterization practical. On the other hand, this program focuses more on the extraction of all four leading-twist GPDs including both valence and sea quarks of different flavors with the combined inputs from both lattice calculations and experiments. Therefore, though in a similar spirit, the GUMP program allows a more comprehensive study of the GPDs, especially the valence parts which are effectively constrained by lattice calculations while less accessible from experiments.

Inputs and fitting strategy
With the above parameterization, the global analysis can be performed with adequate inputs to pin down all the parameters. In this subsection, we will discuss how we select inputs from all the available results for the global analysis.
We start with the forward inputs where GPDs reduce to PDFs. While we could have various inclusive measurements as the inputs rather than the global PDFs extracted from a specific work to avoid the bias, there are several reasons we choose the extracted PDFs for the GPD global analysis. First, the extraction of the PDFs themselves is considerably involved and requires dedicated works, especially for the simultaneous extraction of both unpolarized and polarized PDFs [49,50]. Second, much less is known for the off-forward behaviors of GPDs compared to the forward ones. Consequently, a fine forward analysis can not be matched with an equivalent off-forward analysis, making it less urgent to improve the forward part of the analysis. Therefore, we avoid repeating the forward fittings that have been studied extensively by the PDF global analysis community and take their globally extracted PDFs as the inputs. More specially, we consider one of the recent analyses by the JAM collaboration [49], where both the unpolarized and polarized proton PDFs are extracted simultaneously.
The off-forward inputs, on the other hand, consist of various constraints including globally extracted FFs [22], lattice calculations [17,23] and exclusive measurements, or the DVCS measurements [9][10][11][12]37] more specifically for this work. For the form factors, we take the globally extracted charge FFs [22] rather than fit to the elastic scattering data directly [35] for the same reasons as for the PDFs. Since the charge form factors are averaged over all flavors, we combine the charge form factors of both proton and neutron to obtain the charge form factors of up and down quarks respectively, assuming isospin symmetry and ignoring the contributions of strange and other heavy flavors.
The lattice inputs include both calculations of generalized form factors and x-dependence of GPDs [17,23]. Putting the lattice QCD results together with other experimental measurements can be quite subtle, due to the systematic uncertainties of the lattice results that are hard to fully understand. Progress has been made in getting the systematic uncertainties under control and calculations from different groups are converging for certain quantities like the (axial) charge form factors, see for instance the review in ref. [51], but more efforts are still in demand to get the other results consistent. Therefore, in this work we consider the lattice results from a single group only [17,23] 2 to avoid the potential tension among lattice results from different groups. On the other hand, we adjust their weights in the global analysis by increasing the uncertainties of the results to account for the unknown systematic uncertainties. We note that the calculation of the x-dependence of GPDs at non-zero skewness has been done in the literature [17], although the reliable regions with controlled systematic uncertainties get even more subtle due to the irregular behavior of GPDs at x = ξ. Thus, the inclusion of the x-dependence of GPDs at non-zero skewness in the global analysis will be left to the future work.
Last but not the least, we have the experimental exclusive measurements. In principle, the exclusive measurements can and should include as many processes as possible to get better constraints on the GPDs as well as to test the universality of GPDs. However, the main challenge for putting different processes together is the mismatch in the amount of data available. Two types of Compton scattering processes, DVCS and Time-like Compton Scattering (TCS) [52], have been measured at JLab. Much more DVCS data have been obtained than that of TCS, of which the first measurement was made just recently [53]. On the other hand, DVMP typically requires much larger virtuality than DVCS to suppress the unwanted contributions of transverse polarization, which are mostly accessible with colliders such as HERA [45,54] and the future EIC [13,14]. Therefore, we start by considering the DVCS measurements only, which has the broadest global kinematical coverage among these process that includes both the low x B region covered by H1 [37] at HERA and the medium x B region covered by various experiments at JLab [9][10][11][12]. It is worth noting that the DVCS process is mostly sensitive to the quark distributions whereas the gluon distributions are best obtained from meson production or other gluon-sensitive processes. Obtaining the gluonic distributions from them is indeed of high interest and importance, which will be carried out in a separate work.
With all the inputs above, it seems straightforward to simply put them together and perform the global analysis with the parameterization described before. However, it will not be quite practical with the large number of free parameters, which typically means extremely slow convergence in the search of best-fit parameters. Besides, the evaluation of GPDs, especially with scale evolution, is much more computationally intensive, compared to that of the PDFs for instance, making the global analysis even less handleable. Therefore, we split the global analysis into two steps by first performing the semi-forward fit at zero skewness ξ = 0 (while the momentum transfer squared t can still be non-zero) and then fit to the off-forward inputs with non-zero skewness ξ = 0. Apparently, the twostep fitting introduces bias that favors the semi-forward constraints over the off-forward constraints. However, given the fact that much more semi-forward constraints can be obtained from lattice and experiments than the off-forwards ones, it is a reasonable and practical assumption for such a large system of parameters. Figure 1: A plot showing the fitting procedure of the GUMP program. The semi-forward (ξ = 0) fit that consists of four individual fits for each t-dependent PDF is performed before the off-forward (ξ = 0) fit. The results of the semi-forward fit are then passed to the off-forward fit as fixed parameters, whereas the off-forward parameters will be determined with the off-forward constraints.

Global analysis and extracted GPDs
In the previous section, we discussed the parameterization of GPDs as well as the inputs to constrain the parameters, so we would need to find the set of parameters that fit to the measurements best, for which we employ the iminuit interface of Minuit2 [55,56] as the minimizer. In this section, we will present the results of the global analysis as well as the extracted CFFs and GPDs.

Basics of the fit
As discussed before, the whole fit will be split into the semi-forward (ξ = 0) part and the off-forward (ξ = 0) part to avoiding dealing with the huge parameter set in a single fit. Furthermore, in the semi-forward case constraints on different species H, E, H, and E decouple, unlike the off-forward case where all four of them are involved. This allows one to further decompose the semi-forward fit into four separate fits that do not interfere for each of the four t-dependent PDFs. Therefore, the whole fitting procedure eventually consists of five individual fits as shown in figure 1. Correspondingly, the total χ 2 can be decomposed into five parts: where each term χ 2 F with F = {H, E, H, E} stands for the χ 2 for the semi-forward fits of the corresponding t-dependent PDFs. Then one just needs to minimize these χ 2 s separately in order to perform the fit and the results are summarized in table 2. More details of the fitted parameters are presented in Appendix A.
Several comments are as follows: first, we take the globally fitted unpolarized and polarized PDFs [49] with 31 points sampled in the region x ∈ [0.005, 0.6] for each flavor, indicating 155 points for the t-dependent PDFs H and H respectively. Note that since the parameterization we used has only 3 parameters in the forward limit for each PDF, much less than the PDF global analysis [19], the sample size can not be too large correspondingly. The extra constraints from PDFs on H and H account for their larger N data as shown in table 2, while the other around 50 data for each t-dependent PDF are taken from the globally fitted FFs [22] and lattice calculations [17,23].  Table 2: A table summarizing the total χ 2 versus the number of data points N data as well as the χ 2 per degree of freedom χ 2 ν . Since the JLab and H1 DVCS measurements are fitted simultaneously in the off-forward fit, the χ 2 /ν is replaced with χ 2 /N data instead to estimate their separate contributions.

Sub-fits
Second, in the interest of this work, we select the JLab DVCS data with larger Q 2 (Q > 1.8 GeV) and smaller x B (x B < 0.5). The larger Q 2 is required by the factorization theorem and also to suppress the higher twist effects, while the low x B region is selected in accord with the small ξ(x B ) expansion such that we have the expansion parameter ξ 2 0.1. Even with such selection, it still leaves much more JLab DVCS data than the H1 data (which do not need selection as they are already in the large Q 2 and small x B region), since both the azimuthal φ dependence and beam-polarized cross-sections are measured at JLab. Therefore, the JLab data will form the dominant input in the off-forward fit even though they cover about the same number of kinematical points on the (x B , t) space. In figures 2, we show the kinematical coverage of the DVCS measurements at both JLab and H1 and present some typical examples of the fit to the DVCS measurements.
To summarize, the reduced χ 2 s are all around 1 which indicates generally good fits. For the forward fit, the main issue is that the 3 forward parameters for the H PDF cannot perfectly describe the x-shape ranging from x ∈ [0.005, 0.6]. Although this could be improved with more sets of model ansatz, it would also require more off-forward inputs accordingly. As for the off-forward fit, the main challenge seems to be the potential highertwist effects in the JLab DVCS measurements given that the typical Q 2 of them is around 4 GeV 2 , which still allows sizeable twist-three contributions or even twist-four contributions especially at large momentum transfer [12,58,59]. For instance, the LU cross-section plot at the lower-right of figure 2 still deviates from a perfect sine shape as predicted in the Q 2 → ∞ limit, which could be due to the higher-twist effects.  [22]. The DVCS and interference cross-sections are calculated with the formulas in ref. [57] including the contributions of twist-two CFFs only.

The extraction of Compton form factors
As discussed before, since the off-forward inputs contain DVCS measurements only which are effectively combinations of CFFs, they do not constrain the x shape nor the flavor structures of GPDs well. Thus, the flavor-dependent x shape of GPDs extracted under the empirical constraints in subsection 2.2 will be obviously model-dependent. Therefore, before moving on to the extracted GPDs which relies on the choice of ansatz, we first discuss the extracted CFFs of which the comparison is less model-dependent. More specifically, we will compare the CFFs extracted in this work with the locally extracted CFFs in ref. [12] and the CFFs predicted by the KM15 model [60]   Before commenting on the comparison of the extracted CFFs, we first note that even the local extraction of CFFs suffers from the degeneracy issue. The total DVCS crosssections are generally quadratic equations of all the 8 CFFs: with F i = {ReH, ImH, ReE, ImE, Re H, Im H, Re E, Im E} which are the real or imaginary parts of the CFFs corresponding to the GPDs H, E, H and E. For each combination of beam polarization P b and target polarization P t , the pure DVCS (DVCS), interference (INT) and Bethe-Heitler (BH) contributions are quadratic, linear, and constant in the CFFs respectively. Ideally, one would need all 8 possible combinations of P b and P t in order to disentangle the 8 CFFs, see for instance ref. [61], but there could still be degeneracy in the solutions as the nature of quadratic equations. In this work, the degeneracy will be more severe, since only two polarization configurations, unpolarized or polarized beam with unpolarized target (UU and LU), are considered. For instance, one can show that with UU and LU cross-section, the CFF E only shows in the quadratic terms, multiplied to either itself or H, implying that the quadratic terms are invariant under the transformation and the same for the imaginary part. This degeneracy of CFFs in the DVCS cross-sections leads to the ambiguity in the sign of the extracted Re H and Re E as well as the Im H and Im E shown on the right of figure 3, where the extracted CFFs in this work seem to take opposite sign compared to that of the local extraction. 3 Besides this explicit example, there might be other implicit degeneracies in the extracted CFFs which could affect the reliability of such an extraction. Although such degeneracy makes the comparison of the extracted CFFs more subtle and less intuitive, it can certainly be improved in the future with more polarization configurations taken into consideration. On the other hand, many of the CFFs extracted in this work agree well with the local extraction as shown in figure 3, adding more confidence to the extraction of those CFFs.

The extracted GPDs at non-zero skewness
Compared to the CFFs extraction discussed in the previous subsection, the extraction of GPDs will involve more model dependence. The lack of off-forward constraints except the CFFs is the main challenge in the GPD extraction, and the extracted x-shape of GPDs will depend on the ansatz chosen correspondingly. In addition, as the CFFs are averaged over different flavors, the extracted flavor structures could be ambiguous as well. Though suffering from these ambiguities, we still present the extracted GPDs here as an illustration of how GPDs are constrained by the inputs, with the caveat that the GPDs, especially their off-forward behaviors, are not uniquely determined at this point. These results can certainly be improved with more lattice calculation at non-zero skewness as well as more flavor-sensitive data in the future.
In figure 4, we present the extracted PDFs H, E and H, and E at ξ = 0, t = 0, and the reference scale µ 0 = 2 GeV. The two PDFs H and H are fully parameterized and fitted to the globally extracted unpolarized and polarized PDFs in ref. [49], and therefore they agree well with the reference values there, which are not shown in the plots though. On the other hand, due to the lack of the information, the PDFs E and E are extracted with the empirical constraints summarized in table 1 with globally extracted FFs and lattice calculation of form factors as well as GPDs. The E PDFs turn out to be quite significant due to the contributions of the pion pole, according to the lattice calculated form factors of E [23] that we fit the E PDF to. As for the E PDFs, the shape of the E PDFs are obtained combining the lattice calculation of E GPDs and other relevant form factors from both experiments and lattice results. The u and d quark E GPDs are almost the opposite of each other, in accord with the observation that the flavor isoscalar form factors B u+d and gluonic gravitational form factors B g are consistent with zero according to the lattice results [23,62,63]. We note that the valence part of the above PDFs are obtained in the semi-forward fit of the previous work [36] as well. However, they are quantitatively slightly different from the previous results because of the different constraints used here as well as the effects of the sea quark distributions that were not considered in the previous work. The gluon PDFs are also implemented in this work, and they enter the scale evolutions and mix with the quark GPDs. However, since they are not as constrained by the DVCS measurements, they are not presented here.
In figure 5, we present quark GPDs at ξ = 1/3 and t = −0.69 GeV 2 and reference scale µ 0 = 2 GeV. In general, the extracted GPDs oscillate in the DA-like region (−ξ < x < ξ) with a damping tail in the PDF-like region (x > ξ or x < −ξ), as the result of the conformal partial wave expansion which expands the GPDs in the DA-like region in terms of Gegenbauer polynomials that oscillate. We note that due to the DA terms of GPDs, namely the F qq terms in eq. (2.1), the behavior of GPDs in the DA-like region could look very different from here. With off-forward inputs mainly just CFFs, only the GPDs at the crossover line x = ±ξ are effectively constrained, which are then extrapolated with decaying tails to the PDF-like region, whereas the DA-like regions are not uniquely determined which we will discuss more in the next subsection.
The above extracted PDFs and GPDs are generated with the open-source codes of this program [64] which could be used to generate other observables including DVCS crosssection measurements as well, although we should note that only the central values of them are available at present. While the error propagation is a crucial part of a global analysis, the process is extremely computationally intensive. With 20 functions of three variables where each point must be calculated through a numerical contour integral, sampling them over more than 50 parameters for the statistical uncertainties would be very challenging though not as useful since the main uncertainties are still from the systematics. Therefore, we will leave the error estimation to the future works once the task will be more practical.

GPDs in the DA-like region and DA terms
In the end of this section, we discuss more about the GPDs in the DA-like region. Recall that in the previous section, we discussed that in the small x B expansion the DA-like region gets less relevant, and so we did not consider the DA-terms F qq in the global analysis. However, the behavior of GPDs in this region is of genuine interest as well, which is less known compared to our knowledge of GPDs in the PDF-like region e.g., from PDFs. This region of GPDs is equally important and can be accessed directly from lattice calculations at non-zero skewness [17]. We note that these constraints on GPDs in the DA-like regions were not imposed in the global analysis since only very few of them are available at present. However, in this subsection, we will discuss how the current framework under small x B expansion can be extended to accommodate these constraints with the extra DA-terms F qq and produce smoother GPDs in the DA-like region.
In eq. (2.3), we showed that GPDs can be expressed as the sum of their conformal partial wave, where each term is given by a rescaled Gegenbauer polynomial [30]: that is non-zero only the in DA-like region. Therefore, one can always add finite terms like this to the GPDs freely without changing the GPDs in the PDF-like region. Such terms are called the DA terms in this work. Since the Wilson coefficients are antisymmetric in x for the CFFs H and E, adding terms with even j which are symmetric in x will keep both the CFFs H and E invariant, which correspond to the so-called shadow GPDs [46]. On the other hand, adding terms with odd j which are antisymmetric in x will modify the real part of the CFFs H and E while keeping their imaginary parts the same, which are also known as the subtraction terms [31]. Similar arguments apply to the H and E CFFs too, except that their Wilson coefficients are symmetric in x. These terms are in principle hard to extract from experiments, however, they do affect the generalized form factors and can be constrained by the FFs as well as lattice calculations. Therefore, one should parameterize these extra DA terms, at least part of them, and fit them to these constraints in order to obtain the shape of GPDs in the DA-like region. There have been lattice calculations of the GPD shape [17], which could be used to constrain such terms and determine the shape of GPDs in the DA-regions. However, since the results contain only the isovector u − d combination of H and H GPDs and such calculations will break down at x = ±ξ, they do not pose enough constraints on the GPDs for global analysis. Therefore, we will leave the extra fitting of the DA terms to those constraints in the future work with more information in the DA-like region. In figure 6, we show an example of how GPDs can be tuned with the extra DA terms to fit to the constraints in the DA-like regions, where we take the isovector GPDs H u−d calculated on lattice [17] as the reference value. Figure 6: The isovector GPD H u−d tuned to fit the reference value calculated on lattice [17] at reference scale µ 0 = 2 GeV. The dashed line is the original extracted shape of GPDs. By adding extra terms in the DA-like region, we obtain the tuned GPDs as the red curve which approach the reference value shown as blue blocks.
In figure 7, we also show the tuned isovector GPD H u−d on the (x, ξ) plane at t = −0.69 GeV 2 and reference scale µ 0 = 2 GeV. We again note that such results are obtained under the ansatz and empirical constraints used in this work.

Conclusion and outlook
We extend the previous work [36] of the GUMP program to the non-zero skewness case and perform the global analysis of quark GPDs combing experimental measurements of DVCS, relevant lattice calculations for GPDs and PDFs and FFs from global analysis for the first time, whereas the gluon GPDs will be carried out in a separate work with other gluon-sensitive processes such as DVMP.
We argue that empirical constraints are still needed for the global analysis of GPDs, given the extremely large system of GPDs that one needs to consider and the limited knowledge of them currently. With these empirical constraints, we extract the GPDs from global analysis with the above inputs and present the globally extracted PDFs, CFFs and GPDs, with the caveat that more inputs, including more polarization configurations for the DVCS measurements and more lattice results of GPDs at non-zero skewness, are still needed to improve the reliability of such an extraction.
We also discuss the general framework to extend the current program which focuses on the small ξ region of GPDs to allow the analysis of the DA-like regions of GPDs that will be more relevant at larger ξ. We argue that besides the quark and antiquark GPDs, the extra Figure 7: The isovector GPD H u−d on the (x, ξ) plane tuned to fit the reference value calculated on lattice [17] at t = −0.69 GeV 2 and reference scale µ 0 = 2 GeV. The blue (x = ξ) and yellow (x = −ξ) curves correspond to GPDs on the two crossover lines respectively. A cut at ξ = 0.1 is made since GPDs in the DA-like region −ξ < x < ξ get singular when ξ approaches 0.
DA terms are crucial in describing the GPDs in the DA-like regions, which can improve the parameterization with more flexibility and without damaging the physical constraints like polynomiality conditions. We present an example how the DA terms could modify the GPDs in the DA-like region while keeping the CFFs the same, and therefore can be used to parameterize and fit the shape of GPDs in this region.
In the future works, we will consider the meson production processes in the global analysis to better constrain the gluon GPDs. In addition, we also consider adding the strange quark distributions to the analysis which might have sizable effects. Besides, we will also include the DA terms in the global analysis to fit the GPDs in the DA-like region once enough constraints on the GPDs in the DA-like region are obtained.

A GUMP parameters and their best-fit values
In this appendix we present more details of the fit, especially the best-fit parameters obtained from the global analysis. In  We note that the N, α, β, α , b are the parameters in the moments of GPDs according to eq. (2.6) where α(t) ≡ α+α t corresponds to the linear Regge trajectory and β(t) = exp(bt) corresponds to the extra exponential term. Each of these parameters have a superscript representing its GPD species (H, E, H or E) and a subscript representing its flavor (u V , u, d V ,d or g). The R E sea and R E sea are the ratio of the E( E) GPDs to the H( H) GPDs for the sea quarks and gluons. Besides them, there are also parameters like R H u,2 that are the ratio of the ξ 2 terms to the forward terms as defined in eq. (2.7) for GPDs with different species and flavors. We again note that more details of them and the GUMP program that generate the GUMP GPDs, CFFs and cross-sections are available online [64].