Generalized parton distributions through universal moment parameterization: zero skewness case

We present a global analysis program for the generalized parton distributions (GPDs) based on conformal moment expansion. We apply the strategy of universal moment parameterization to fit both the collinear parton distribution functions (PDFs) from phenomenology and generalized form factors from lattice calculations, and show that the parameterization is flexible enough to accommodate these constraints. In addition, we can also fit direct lattice calculations of GPDs from large-momentum effective theory. In this work we focus on the analysis of $t$-dependent PDFs which correspond to GPDs in the $\xi \to 0$ limit. The strategy also applies to the $\xi \not =0$ region with extra parameters, and therefore can be fitted to experimental observables in the future. With a demonstrative example of fitted GPDs, we exhibit the quark transverse angular momentum densities of the proton as well as the impact parameter space distributions of quarks in both unpolarized and transversely polarized protons.

Despite the efforts put into investigating GPDs, limited knowledge of them has been obtained at this point due to their higher-dimension nature. More specifically, GPDs, collectively denoted as F , are functions of (x, ξ, t) with x the parton momentum fraction, ξ the skewness parameter describing the momentum transfer in the longitudinal direction and t the total momentum transfer squared. Typically, we have constraints from: parton distribution functions (PDFs) in the forward limit ξ = t = 0, generalized form factors (GFFs) with t-dependence and Compton form factors (CFFs) with ξ-and t-dependence. None of these place 3D constraints, causing a huge gap in the determination of GPDs.
Due to the lack of knowledge, we are faced with the so-called inverse problem when reconstructing GPDs, which states that GPDs can not be uniquely determined with, for instance, the combination of PDFs and the CFFs [42]. Therefore, we are in need of a global analysis program that can combine different constraints on GPDs to solve the inverse problem. One possible way is to parameterize the GPDs with finite (and extendable) numbers of free parameters, such that one can match the number of free parameters to the number of constraints on the GPDs till they are completely determined with controlled uncertainties. This is in the same spirit of determining proton collinear PDFs from various experimental data [43,44].
With the above motivation, lots of GPD parameterization methods are proposed in the literature [45][46][47][48][49][50][51][52][53][54][55]. They generally fall into three categories: the double distribution representations, the moment space representations and other dynamic GPD models. The Radyushkin double distribution ansatz (DDA) is one of the earliest GPD parameterization method [45,46]. Under the double distribution framework, the Vanderhaeghen-Guichon-Guidal (VGG) model [48,56,57] and the Goloskokov-Kroll (GK) model [49,58] were developed for phenomenological applications. In general, the DDA expresses the GPD in terms of an integral expression of the double distribution. Then by constructing proper double distributions, the GPD will be reduced to the corresponding PDF in the forward limit and satisfies the polynomiality conditions. However, since the original double distributions have only one free parameter to control the overall ξ dependence, such GPDs parameterization methods lack the flexibility when fitting to measurements with different ξ [59].
The limitation of DDA leads to explorations of other possibilities. Another approach to constructing GPDs is to consider the polynomial expansion of GPDs. Two major GPD parameterization methods based on the polynomial expansion are the dual parameterization [47,60,61] and the Kumerički-Müller (KM) model [50,52] commonly used in DVCS analysis, which are shown to be mathematically equivalent in ref. [62]. The basic idea of the polynomial expansion or the moment expansion is to expand the GPD in terms of a complete set of polynomials where the expansion coefficients are the corresponding moments. Then the GPD can be expressed as the sum of these polynomials without loss of generality. The polynomial expansion requires the polynomials to properly imitate the behavior of the GPD and extra care will be needed for practical uses.
In addition, there are various dynamic models of GPDs [53-55, 63, 64]. These models generally approximate the hadron states as expansions in the Fock space of quarks and gluons, usually with the help of light front wave functions. Since the parton distributions are the matrix elements of certain light-cone operators between hadrons, they are calculable when the hadron states are replaced with the bound states of quarks and gluons with some effective interactions. In this sense, the parameterization of GPDs is turned into the task of determining the free parameters in the wave functions and the effective interactions.
However, a general parameterization strategy that combines various physical con-straints on GPDs with the lattice calculations still seems to be lacking in the literature. Therefore, in this work we present a program to parameterize GPDs through universal moment parameterization (GUMP), which can be fitted to lattice calculations, PDFs and experimental measurements of CFFs altogether and help solve the inverse problem. We will continue exploring the moment expansion, especially the conformal moment expansion, for reasons to be explained. Our work here is closely related to the KM model, but the goal of our program is quite different. The KM model aims to just fit to the experiment data, and thus only the sea distributions are parameterized by the general moment expansion, whereas the valence distributions are given by some functional forms on the crossover line only. More efforts are put into the sea rather than the valence distributions in the KM model because experiments such as DVCS can only measure the CFFs, which are mostly from the sea distributions that dominate at small x or ξ. As for this program aiming to combine all the possible constraints on GPDs including both experiments and lattice calculations, both sea and quark distributions are important. Therefore, a general extension of the previous parameterization methods to all GPDs is necessary in order to combine the different constraints. The organization of the paper is as follows. In section 2, we discuss the physical properties of GPDs and introduce the well-established conformal moment expansion of GPDs. In section 3, we make our assumptions based on phenomenological consideration. In section 4 we apply them to build the t-dependent PDFs and study the nucleon 3D structures with the fitting results. In the end, we conclude in section 5.

Review on conformal moment expansion of GPDs
In this section, we will introduce the general technique and review salient features of conformal moment expansion for GPDs [50], which will be our basis for the GUMP program. We start with the physical constraints on GPDs which are important guides to finding proper parameterization strategies.

Physical constraints on GPDs
GPDs are simply functions F with 3 arguments (x, ξ, t), which can be defined as the offforward matrix element of light-cone operators as, with n the light-cone vector (n 2 = 0), O(λn) the light-cone operators of the quark/gluon and P /P the initial/final nucleon momenta. We define the average nucleon momentum P ≡ (P + P )/2, the momentum transfer ∆ ≡ P − P and the momentum transfer square t ≡ ∆ 2 . Besides, the skewness parameter ξ is defined as ξ = −n · ∆/(2n ·P ). In the forward limit, GPDs are reduced to the corresponding PDF f (x), so we have which is the most important and stringent constraint on GPDs. The two arguments x and ξ of GPDs have support x, ξ ∈ [−1, 1], whereas the momentum transfer squared t could take any negative value t < 0. Since GPDs are typically defined symmetric in ξ, we will consider only positive ξ > 0 without loss of generality. GPDs do not exist everywhere in the above 3D subspace. For instance, since the skewness parameter ξ corresponds to the longitudinal projection of the momentum transfer, there exists a minimum |t| min for each ξ as with M the nucleon mass such that it takes infinite momentum transfer |t| → ∞ to reach |ξ| → 1, as shown in figure 1. Accordingly, we restrain ourselves to the kinematically allowed region. Furthermore, there exists non-analyticity on the (x, ξ) plane, specifically on the crossover line |x| = ξ. Since GPDs are essentially parton distributions with non-zero momentum transfer ∆, their physical interpretation varies accordingly [3]: with small momentum transfer in the near forward region ξ < |x|, GPDs can still be interpreted as the amplitudes of emitting and reabsorbing a parton in the nucleon and behave like the PDFs, whereas with large momentum transfer in the region ξ > |x|, they are interpreted as the amplitudes of emitting/absorbing a pair of partons and behave like the distribution amplitudes (DAs). The two regions will be referred to as the PDF region (ξ < |x|) and the DA region (ξ > |x|) DA PDF PDF respectively, and they are joined together at the crossover line ξ = |x| where one of the partons has zero momentum fraction. The collinear factorization requires the GPDs to be continuous at the crossover lines, but their derivatives might not be. Consequently, these crossover lines can be extremely subtle for the parameterization of GPDs on the (x, ξ) plane.
In addition, GPDs are also constrained by the polynomiality condition, which states that the Mellin moments of GPDs must be finite order polynomials of ξ such that [3], The polynomiality condition is the result of Lorentz invariance, and it holds even though the integral in eq. (2.4) involves both the PDF region (|x| > ξ) and the DA region (|x| < ξ) in spite of the non-analyticity at the crossover lines. While other physical constraints exist for GPDs, the above constraints will be the main focus of this work.

Conformal moment expansion and scale evolution of GPDs
With the physical constraints discussed above, we now introduce the conformal moment expansion and show how the constraints can be imposed with such expansion. Consider if we expand the GPD in terms of a complete set of polynomials, denoted as C i (x), such that with F C i (ξ, t) the corresponding coefficient for fixed ξ and t. Since C i (x) form a complete set of basis, we have the completeness condition and the orthonormal condition with ρ C (x) the weight function of the polynomials C i (x). Then the F C i (ξ, t) is simply the moment of F (x, ξ, t): The above is generally true for any complete set of polynomials C i (x) and there are infinite possible choices. Therefore, one naturally would ask: which choice is the most suitable one for GPDs? An answer is motivated by the renormalization group equations of GPDs. But before getting into that, we start with a simpler case -the scale evolution of PDFs, or the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations [65][66][67].
The leading-order DGLAP equations for non-singlet quark distributions reads, with P (z) the splitting function. Generally, one needs to solve the integro-differential equation in order to perform the scale evolution of PDFs. However, if one performs the Mellin transform of the DGLAP equation according to the DGLAP equation will read, in Mellin space, where q n (Q) and P n are the Mellin transformations of q(x, Q) and P (x) according to eq. (2.11) and the DGLAP equation becomes multiplicative in the Mellin moment space which holds true to all order in perturbation theory. Therefore, if one parameterizes the PDFs in the Mellin moment space, the scale evolution equations can be easily solved.
The above moment representation also applies to the construction of GPDs. Consider the scale evolution equations of non-singlet quark GPDs in the form of [1,10,62,[68][69][70][71][72][73], the evolution kernel. It can be shown that in the forward limit, the evolution kernel is reduced to the splitting function of DGLAP equation: while for |x|, |x | < 1 it is related to the Efremov-Radyushkin-Brodsky-Lepage (ERBL) kernel [74][75][76][77] for the DA scale evolution: For this reason, the PDF and DA regions of GPDs have also been called DGLAP and ERBL regions respectively. Analogous to the analysis of PDFs and their Mellin moments, the GPD moments that diagonalize the right-hand side (RHS) of eq. (2.13) will renormalize multiplicatively. To the leading order, such a solution is found to be the Gegenbauer polynomials C with V (0) the leading-order evolution kernel and γ j the corresponding anomalous dimension (eigenvalue) of each Gegenbauer polynomial. For the quark we have λ = 3/2, and for the gluon λ = 5/2 should be used. Accordingly, the conformal moments defined as are multiplicatively renormalizable, which will not hold if the higher order evolution effects are taken into account and more care will be needed then [51]. Here an extra normalization constant is multiplied to ensure that in the forward limit, the conformal moments are reduced to the Mellin moments, which shall be expected as the evolution equation itself will be reduced to the DGLAP equation. Also note that there exists a mismatch of conventions in the literature that the jth conformal moment F j (ξ, t) actually corresponds to the nth Mellin moment with n = j + 1 instead of j, according to eq. (2.11). In this work, the two will be distinguished by having j for the conformal moments and n or s for the Mellin moments respectively, and they are related by n = s = j + 1. Inspired by the leading-order evolution and with eq. (2.6), in general GPDs can be written in terms of a formal summation as [62], where we have the weight function of Gegenbauer polynomials ρ 3 2 C (x) = 1 − x 2 and some extra normalization constants to match normalization of the conformal moments F j (ξ, t) in eq. (2.17). The whole prefactors of the conformal moments can be collected to define where an extra (−1) j factor is extracted for future convenience. Analogous to the partial wave expansion, the p j (x, ξ)s are called the conformal wave functions and F j (ξ, t) the conformal partial wave amplitude respectively [50]. Also note that with the above conformal partial wave expansion, the non-analyticity at |x| = ξ can be taken care of -GPDs might not have well-defined derivatives there, but the moments still exist such that the polynomiality condition can be imposed.

Resummation of conformal partial wave expansion
In moment expansion, the reconstruction of PDFs or GPDs requires summing over the contributions of all moments. Such summation does not always converge, and analytic methods shall be developed to perform the formal summation.
Consider the Mellin moment representation of the PDF f (x) for example. In order to express the PDF in terms of its Mellin moments, one can write f (x) in terms of a formal summation as which is obviously not summable, since the δ(x) function and its nth derivatives δ (n) (x) are singular. Even though the RHS of eq. (2.21) can be shown to be equivalent to the PDF f (x) formally, the expression is of no practical use. Therefore, the concept of formal summation has to be introduced, which implies that we consider the RHS of eq. (2.21) as an abstract mathematical entity regardless of the convergence of the series.
To actually reconstruct the PDF f (x), the inverse Mellin transformation can be employed: where the Mellin moments f s are analytically continued in the sense that s can take complex values rather than just positive integers n. This method has been adopted for the numerical analysis of PDFs in the literatures [78,79]. Note that the two equations (2.21) and (2.22) define the same PDF f (x), though they look totally different. To resum eq. (2.19) for GPDs, one obstacle is that the Gegenbauer polynomials are only complete on the interval [−1, 1], while the argument of the Gegenbauer polynomials x/ξ in the above expansion can take any real value. For this reason, we restrict the expression eq. (2.19) to the region |x| < ξ, while extension to the PDF region |x| > ξ has to be defined through analytic continuation. This can be done with the Schläfli integral, as discussed in ref. [50,62], where the integration contour encircles the line segment u ∈ [−1, 1]. The analytic continued wave functions can then be derived with this integral representation and one gets [62], These wave functions are equivalent to the Gegenbauer polynomials definition of conformal wave functions when j is integer, and they can be used for more general propose, for instance, the resummation of moments. Note that the above conformal wave functions are continuous at |x| = ξ but their derivatives are not, which is consistent with the behavior of GPDs at the crossover lines. Two techniques, the Mellin-Barnes integral [51,52] and the Shuvaev-Noritzsch transform [60,61] have been used in the literature to reconstruct the GPDs with moments. Both methods are shown equivalent mathematically [62], and we will focus on the former one here. The Mellin-Barnes integral is based on the Sommerfeld-Watson transformation, assuming f (z) is analytic. The proof is straightforward, as the sin(πz) in the denominator of RHS leads to poles at all integers z with the residue of π(−1) n f (z) and the equation is simply the residual theorem applied to the RHS. Then, the formal summation in eq. (2. 19) can be written in terms of the following Mellin-Barnes integral in the same manner [51,52] F where the constant c should be carefully chosen to have all the non-negative poles while the poles of conformal wave functions p j (x, ξ) and conformal moments F j (ξ, t) should be  outside the contour, see more discussions in ref. [62].
The comparison of the different methods used for PDFs and GPDs are shown in figure  3. In principle, the Mellin-Barnes integral should reproduce the PDFs in the forward limit ξ = 0. Indeed, since the conformal wave functions p j (x, ξ)s have the limit behavior the formal summation in eq. (2.19) will be reduced to eq. (2.21). However, due to the singular behavior of δ functions, the Mellin-Barnes integral will not converge in the ξ → 0 limit. Instead, one can use a finite but small ξ as the regulator for the δ function, and the Mellin-Barnes integral will reproduce the PDF when x ξ, as shown in figure 4. Following the above framework, for any (well-behaved) given set of moments F j (ξ, t), the corresponding GPDs can be calculated and then used for all kinds of phenomenological applications. The GUMP program for global analysis follows this approach by parametrizing the general forms of the moments or universal moment parametrization. More details concerning the parameterization of moments will be discussed in the next section.
Note that while the GPDs calculated above are ready for all kinds of applications, there exist some shortcuts, such as the evaluation of CFFs and scale evolution, as these calculations are easier to perform in moment space. The detailed discussion of this topic is in ref. [51].

Universal parameterization of GPD moments
In the last section, we show that the GPDs F (x, ξ, t) can be equivalently expressed in terms of their moments F j (ξ, t) without model assumption. In this section, we will discuss how to find out the universal parameterization of moments.
Generally, we are interested in GPDs of (anti)quarks with different flavors as well as the gluons, and conventionally the quarks and antiquark GPDs are collected together by defining the F (q) (x, ξ, t) such that [62] where F q (x, ξ, t) is the quark GPD and Fq(x, ξ, t) the antiquark GPD with support 1 > x > −ξ. It is common to reexpress them in terms of the valence and sea distributions analogous to the forward PDF as, such that for each quark flavor, the different distributions can be expressed in terms of the corresponding valence and sea distributions.

Polynomiality condition and small ξ expansion of GPD moments
To impose the physical constraints on GPDs, we start by noting that the polynomiality condition, originally proved for the Mellin moments of the GPDs, can be equivalently applied to the conformal moments. Since each Gegenbauer polynomial C (λ) j (x) can be expressed in terms of linear combination of x k , and vice versa: (3.7) Thus the polynomiality condition of Mellin moments immediately leads to the polynomiality condition of conformal moments. The reverse can be proved in the same way using eq (3.6), so the two polynomiality conditions are equivalent. Therefore, we can write the conformal moments of GPDs as where the form factors F j,k (t) with double subscripts do not depend on ξ. While the GPDs F (x, ξ, t) are completely determined by F j,k (t), limited information on them is known. Two limit behaviors, the PDF limit ξ = 0 and the DA limit ξ = 1, can be considered separately. In the semi-forward limit ξ = 0, the conformal moments F j (ξ, t) are given by F j,k (t) with k = 0: where we know that the GPD will be reduced to the PDF and the conformal moments correspond to the Mellin moments of the PDF if we further take t = 0. Thus, the F j,0 (t) at zero momentum transfer t = 0 is entirely determined by the forward PDF f (x), which can always be determined well from global fitting, and we have the following requirements Commonly in the PDF analysis, for instance in the CTEQ global analysis [43], the ansatz of the PDF f (x) is chosen, whereas further improvements can be done based on this. This ansatz can be immediately transformed to the moment space, and we have the following ansatz for the Mellin moments of the PDF which will be our building blocks for the GPD at ξ = 0 or small ξ.
On the other hand, in the DA limit ξ = 1, we have which could have completely different behavior than the forward case. Actually, they must behave differently. The reason is that the formal summation in eq. (2.19) with the moments in eq. (3.12) does not converge, and thus it can only be resummed by the Mellin-Barnes integral method or other equivalent methods. Consequently, GPDs corresponding to such moments do not generally vanish at |x| = ξ. While it is expected as GPDs do not generally vanish at the crossover line, for large ξ, especially for ξ close to 1, such moments lead to GPDs that do not vanish at the end point |x| = 1. Therefore, the moments at large ξ must converge faster than the moments in eq. (3.12) in order to have proper end-point behavior.
Due to the lack of constraints from physically observable to the large ξ behavior of GPDs, we will restrain ourselves to the region with small ξ to avoid such complexity, whereas the large ξ behavior (which also corresponds to the large t behavior as shown in figure 1) of GPDs requires separate studies. Then, though the moments of the GPDs are generally unknown, in the small ξ region, the first few terms will dominate: with k cut the cut-off parameter which controls the ξ dependence, and generally larger k cut is needed for more flexible ξ dependence to contain the larger ξ region. In the region with small ξ, it is assumed that the ansatz used for the PDF in eq. (3.12) applies not only to the leading moment with k = 0, but also the first few moments with k = 2, 4... which can certainly break down as k gets large.
Note that the small ξ expansion does not apply to the GPDs themselves, since GPDs are not analytic at the crossover line |x| = ξ. However, the moments of GPDs are always polynomials of ξ, and they are always regular and expandable in the small ξ region, which, on the other hand, need proper resummation in order to recover the GPDs.

Regge trajectory and t-dependence of GPDs
The last piece to complete the universal moment parameterization is the t-dependence of GPDs. In the framework of moment expansion, the t-dependence is encoded in all the form factors F j,k (t) and thus the t-dependence of the form factors F j,k (t) are crucial. While each term can be in principle fitted to the experimental measurements or lattice calculations, it is not practical to do so, because the resummation procedure requires a set of F j,k (t) with j goes to infinity. Neither experimental measurements nor lattice calculations allow us to access all the moments, especially the higher order ones, and thus some extrapolation methods that take the known moments into account and allow us to extrapolate the higher moments will be necessary. Motivated by that, we introduce the Regge trajectory in this subsection and discuss how it can help build the t-dependence of GPDs.
The starting point is the crossing symmetry, with which processes of different channels but with the same Feynman diagrams can be connected, as shown in figure. 5. The scattering amplitude T can be written with these Feynman diagrams regardless of which channel one chooses as a function of the Mandelstam variables s, t, u, ignoring the spin degrees of freedom for simplicity. The Mandelstam variables s, t, u for 2 → 2 processes are defined as s ≡ (P + q) 2 , t ≡ (P − P ) 2 and u ≡ (P − q ) 2 . The 3 Mandelstam variables s, t, u satisfy the condition s + t + u = P 2 + P 2 + q 2 + q 2 , with which we can eliminate the variable u for simplicity.
One key hypothetical property in Regge theory is the duality [80], which states that the scattering amplitude is given by the sum of all t-channel resonances and also equals to the sum of all s-channel resonances. Therefore, the t-dependence of the scattering amplitude can be tentatively written as with g t (s, u) the residuals of the resonance poles and the RHS sums over all the t-channel resonances. In principle, the scattering amplitude can then be written directly with all the resonances. However, there are potentially infinite resonances and summing over their contributions are highly non-trivial. Thus, one will need to collect the poles together in order to sum over them, which leads to the Regge theory [81,82]. Regge theory is, in short, a theory about complex angular momentum J (different from the conformal spin j). The key observation is that the mass of the resonances increases with increasing angular momentum J on the hadron spectrum. Thus, a linear Regge trajectory usually written as α(t) ≡ α + α t with α(t) the angular momentum and t the mass squared, also known as the Reggeon. The major assumption in Regge theory is that the resonances and their analytic continuation to complex J are the only singularities in the scattering amplitudes. Then the scattering amplitude can be expressed in terms of all the Regge trajectories or Reggeons.
Note that the scattering amplitude T (s, t, u) is not explicitly a function of J. To separate the contributions of different J, it can be written in terms of the partial wave expansion of a t-channel scattering as (3.16) Then the partial wave amplitudes T J (t, u) will have poles at J = α(t) for each Regge trajectory according to the assumption made for Regge theory, and they can be written with these poles as which sums over all the Regge trajectory and r(t, u) is the corresponding residual. Note that eq. (3.15) and eq. (3.17) are equivalent descriptions of those resonance poles. Their relation is intuitive -eq. (3.15) approaches those poles with a fixed integer J in the M 2 direction, whereas eq. (3.17) approaches those poles with fixed M 2 in the J direction. Then, we can discuss the t-dependence of the GPDs or equivalently the form factors F j,k (t). Since they are essentially the nucleon matrix elements of certain quark (gluon) operators which are closely related to the quark (gluon)-nucleon scattering amplitude, the above arguments also apply to them, at least qualitatively. Consider the partial wave expansion in eq. (3.16). In the higher energy limit s → ∞ and ξ → 0, the cos(θ t ) can be expressed in terms of ξ as Then with the asymptotic behavior of the Legendre polynomials, 19) one finds that the partial wave amplitude with angular momentum J has the asymptotic behavior of ξ −J in the high energy limit. Then it can be quickly identified that the term F j,k (t) corresponds to the angular momentum J = j + 1 − k, which leads to , (3.20) with r(t) the corresponding residuals of the Regge poles. Again, if we stay in the region where momentum transfer |t| is small, the leading lightest Regge trajectory is sufficiently good to describe the t-dependence of the GPDs, while more terms can always be added to make it more flexible. Therefore, we reach the finial proposed ansatz for the moments of GPDs as with k = 0, 2, 4 which breaks down at large k. Note that the residuals r i,j,k (t) are smooth functions of t. In order to restore the corresponding forward limit in eq. (3.12), it should satisfy r i,j,k (t = 0) = j + 1 − k − α i,k . However, given the present amount of data, it should be safe to set it to be simply constant: which can be further improved once more data are available. Then there are essentially 4 parameters for each term in the summation on the RHS of eq (3.21).
In the end, note that the proposed form in eq. (3.21) is equivalent to the ansatz in KM model [52] and the similar proposed forward-like function in the dual parameterization [62] except for the SO(3) partial wave expansion that we dropped for simplicity. However, the main goal of this work is to focus on both sea and valence distributions and the corresponding constraints from lattice calculation in for instance [32][33][34][38][39][40][41] besides the experimental measurements which only partially constrain the GPDs.

A demonstrative example: t-dependent PDF and proton tomography
Following the above global analysis strategy for GPDs, in this section we show an example to apply it to the ξ = 0, t-dependent PDF and proton tomography from various constraints. Note that once the parameterization is specified, in principle one just needs to put all the constraints together, and try fitting to them with for instance, χ 2 minimization. However, the amount of parameters needed for a global fitting is huge -we needed at least 4 parameters for each different GPDs (H, E,H andẼ) and each flavors (u, d, g and even more) at ξ = 0 and even more for non-zero skewness ξ. Therefore, as a proof of principle we start by considering the H and E GPDs of u and d quarks at ξ = 0, which are the so-called t-dependent PDF H u/d (x, t) and E u/d (x, t). These two quantities are necessary to understand the angular momentum densities of u and d quarks [2] which will be our main focus in this section.

Fitting the valence t-dependent PDF
In this subsection, we use the proposed form in eq. (3.21) for the t-dependent PDF H u/d (x, t) and E u/d (x, t). We take k = 0 in the semi-forward ξ = 0 limit, and set i max = 1 for the simplest ansatz. Since in lattice QCD the flavor singlet u + d and non-singlet u − d combinations are commonly separately calculated, we apply such ansatz for the H/E u+d (x, t) and H/E u−d (x, t) respectively. We will neglect the sea contributions in the fitting to lattice results here since the sea contributions to GFFs are suppressed compared to the valence contributions. According to the CTEQ18 PDFs [43], the moments of sea distributions are about one order of magnitude smaller than the valence ones, where similar suppression is expected for the GFFs. Therefore, one should not expect strong constraints to the sea distributions from lattice results, and the experimental measurements are crucial to get the sea distributions.
For the t-dependent PDF H val u/d (x, t), we can simultaneously fit to the forward PDF q(x), for which we choose the recommended CTEQ18qed results [43,44], the quark GFFs [32] and lattice calculated isovector t-dependent PDF H val u−d (x, t) [40]. We used the ManeParse [83] Mathematica package to read the CTEQ18qed PDF and the MINUIT2 minimizer [84] for χ 2 analysis, and the results are given in figure 6, where the generalized form factors are thus only A n0 (t) and B n0 (t) survive in the semi-forward limit ξ = 0. Note that the χ 2 per degree of freedom (d.o.f.) is 3.0 for our 4 parameters fitting. The major source of χ 2 is the discrepancy between lattice calculated form factors and the global fitted PDFs -the moments of PDF from global fitting are significantly lower than the lattice calculations. Therefore, the χ 2 per d.o.f. indicates the systematic uncertainties including the effects of unphysical pion mass of m π = 496 MeV, omitted disconnected diagrams and other effects in lattice calculation that are not taken into account.
Similar fitting can be done for the isoscalar valence distributions H val u+d (x, t) as shown in figure 7. However, since the isoscalar t-dependent PDF H val u+d (x, t) is not calculated in ref. [40], we fit to the forward PDF H val u+d (x, t = 0) and the lattice isoscalar form factors A u+d 10 (t), A u+d 20 (t) and A u+d 30 (t) [32]. Again, the major source of χ 2 is the discrepancy between lattice calculated form factors and the global fitted PDFs which do not indicate the failure of the parameterization form used and can be mitigated if the systematical uncertainties are under control.
For the t-dependent PDF E val u/d (x, t), the forward constraints from PDFs no longer exist, and we have to rely on the lattice calculated t-dependent GPDs E val u/d (x, t) [40]. Though the isoscalar GPD E u+d (x, t) is not calculated there, there is strong evidence showing that the isoscalar form factors B u+d (t)s are consistent with zero [32,34]. Therefore, we will assume E val u+d (x, t) = 0 with which the flavor separation becomes trivial. The results for the E val u−d (x, t) are shown in figure 8. The χ 2 per d.o.f. is about 0.8, which is really due to the lack of precision in the lattice results, since there are no constraints from better determined forward PDFs. Therefore, the uncertainties of the GPD E val u−d (x, t) are significantly larger.

The angular momentum sum rule and proton tomography
With the fitted t-dependent PDFs H u/d (x, t) and E u/d (x, t), interesting information of the nucleon 3D structure can be revealed. On the one hand, as mentioned at the beginning of the section, these GPDs carry information about the transverse angular momentum densities of quarks. On the other hand, as first pointed out by M. Burkardt in refs. [6,7], the GPDs in the ξ → 0 limit can be Fourier transformed in the transverse space to get the impact parameter space parton distribution which describe the 3D structure of the nucleon. In this section, we will explore the proton 3D structure with the above fitted GPDs.
In figure 9, we plot the transverse angular momentum densities of quark with eq. (4.1). Note again that, this plot needs to be further improved with the sea distributions, whose contributions to the angular momentum are, however, mostly in the small x region and one order of magnitude smaller than the valence ones. The contributions from the d quark are always negative according to the fitted results, which are essentially the consequence of the assumption E u+d (0) = 0 supported by the lattice calculations of quark and gluon generalized form factors.
In addition to the parton transverse angular momentum densities which only make use of the GPD in the forward limit ξ = t = 0, we can also take advantage of the t-dependent of the GPDs and study the impact parameter space distribution of the quarks. For instance, the quark densities in unpolarized protons are given by the Fourier transform of the GPDs Figure 10: Plots of quark densities ρ q (x, b) for both u (the upper row) and d (the lower row) quarks in unpolarized proton with increasing momentum fraction x. The densities typically get more concentrated but suppressed for increasing x.
where H q (x, b) is defined as the 2-dimensional Fourier transformation of H q (x, −∆ 2 ). Then we plot the quark distribution with the above expression and the fitted GPDs in figure 10. The quark densities typically get more concentrated but suppressed for increasing x, since the number densities should decrease for quarks with larger momentum due to momentum conservation, and the less deflected they will be. Another set of densities that is of interest is the quark densities in a polarized proton. It is worth noting that the transverse angular momentum densities defined in eq. (4.1) are interpreted as the parton angular momentum densities in a transversely polarized nucleon only as pointed out in ref. [7], see also [85][86][87]. The original quark density in a transversely Figure 11: Plots of the intrinsic quark densities ρ X q,In (x, b) for both u (the upper row) and d (the lower row) quarks in transversely polarized proton (in the X direction) with increasing momentum fraction x. Both the u and d quark densities are shifted in the y direction which contribute to the transverse angular momentum J X , while the u contributions are positive (+y direction) and the d contributions are negative (−y direction). polarized nucleon defined in ref. [7] reads with E q (x, b) the 2-dimensional Fourier transformation of E q (x, −∆ 2 ) and the nucleon polarized in the X (distinguished from the parton momentum x) direction. However, with such quark densities, one find so the transverse angular momentum densities can not be recovered this way. It is later found that this is due to the nucleon's center of mass motion -the center of the nucleon has transverse displacement when transversely polarized and thus the longitudinal motion of the whole nucleon will contribute to the transverse angular momentum which does NOT correspond to the intrinsic angular momentum (spin) of the nucleon. This has been discussed in ref. [88] and carefully studied in ref. [89] where the intrinsic quark densities are defined as with which one has With the above fitted GPDs H u/d (x, t) and E u/d (x, t), the intrinsic quark densities in transversely polarized nucleon and be studied as well as shown in figure 11. While the quark densities in a polarized proton get more concentrated but suppressed for increasing x just like in the unpolarized case, we find that clearly the quarks have non-zero transverse displacement, which contributes to the transverse angular momentum of the nucleon along with their large longitudinal momenta.

Sea distributions and experimental observables
Here we discuss the effects of the sea distributions, which are important sources of systematic uncertainties in the above analysis, since in principle both sea and valence distributions contribute to the generalized form factors. As mentioned above, the sea distributions can be dropped because they are mostly in the small x regions and their contributions to the form factors are suppressed compared to that of the valence distributions. On the other hand, this indicates that even if we consider the sea distributions in the above analysis, they will not be constrained as much. Therefore, it is crucial to also include the experimental observables, especially the ones with small ξ (since it can not be zero) in order to constrain the sea distributions. Then, we are inevitably faced with the effects of both none-zero skewness ξ and the sea distributions simultaneously. Fortunately, in the KM model [52], they show that the sea distributions with such parameterization method can successfully fit to the experimental data which proves the flexibility of this parameterization when extending to non-zero ξ. In this sense, our results here with the valence distributions are complementary to the sea distributions obtained in the KM models. On the other hand, GPDs at none-zero ξ are faced with the inverse problem [42], as the experimental measurements of CFFs alone cannot determine the xdependence at non-zero ξ. The solution (compromise) in the KM model is to model the x-dependence of the GPDs at non-zero ξ with the same parameters as the PDFs, except for a different overall normalization constant, which might be too constraining. However, there have been results of lattice calculated x-dependence of GPDs with LaMET [38] at none-zero ξ. These results, even just qualitatively, can be very helpful in mitigating or resolving the inverse problem at non-zero ξ.
Therefore, to obtain the GPDs at non-zero ξ reliably, we need to extend our present results with sea distributions and fit to both experimental observables and lattice results at none-zero ξ as well. A global fitting program is then necessary for this purpose, which will be studied in the future coming paper.

Conclusion
We present a flexible parameterization program for GPDs employing the universal moment expansion which can be fitted to the forward PDFs, the lattice calculation and the experimental measurements when the skewness parameter ξ is small. Under such parameterization, the polynomiality conditions of GPDs can be easily imposed and the non-analyticity at x = ξ will be taken care of. Since the leading order QCD evolution is multiplicative in the conformal moment space, the QCD scale evolution can be performed easily.
We briefly discuss the major parameterization methods of GPDs in the literatures and introduce the necessary techniques for building the GPDs in moment space. Following the previous works on parameterizing GPDs with moments, the t-dependence of GPDs is parameterized with the Regge theory which shows reasonable agreement with lattice calculated generalized form factors. Then by expanding the moments as powers of ξ, the ξ-dependence of GPDs can be gradually extended and fitted to observables. The xdependence of GPDs is the least known except in the forward limit, which needs further inputs, most likely from lattice QCD.
As a proof of principle, we apply our fitting program to construct the quark t-dependent PDFs H u/d (x, t) and E u/d (x, t) by fitting to the forward PDFs, lattice calculated t-dependent PDFs and lattice calculated generalized form factors. With these fitted results, we study the quark transverse angular momentum densities as well as the impact parameter space distributions of quarks in both unpolarized and transversely polarized protons, and show that the quark densities in transversely polarized nucleons are shifted in the transverse direction of impact parameter space, such that the transverse angular momentum can be generated. We also comment on the future development of this program to include sea distributions, non-zero skewness ξ and experimental observables.