Covariant extension of the GPD overlap representation at low Fock states

We present a novel approach to compute generalized parton distributions within the lightfront wave function overlap framework. We show how to systematically extend generalized parton distributions computed within the DGLAP region to the ERBL one, fulfilling at the same time both the polynomiality and positivity conditions. We exemplify our method using pion lightfront wave functions inspired by recent results of non-perturbative continuum techniques and algebraic nucleon lightfront wave functions. We also test the robustness of our algorithm on reggeized phenomenological parameterizations. This approach paves the way to a better understanding of the nucleon structure from non-perturbative techniques and to a unification of generalized parton distributions and transverse momentum dependent parton distribution functions phenomenology through lightfront wave functions.


Introduction
Generalized parton distributions (GPDs) were introduced two decades ago [1][2][3] and have been since then a topic of strong interest in the hadron physics community, both on the theoretical and experimental side [4][5][6][7][8][9][10]. A significant amount of beam time of the upgraded Jefferson Laboratory facility is dedicated to their experimental study and they are one of the core scientific cases for building the U.S. electronion collider (EIC). Not only are they an elegant way to unify both parton distribution functions (PDFs) and form factors into a single object, leading to a three-dimensional picture of hadrons [11], but they also provide genuine information and constraints through their skewness dependence. If these constraints have been used to produce phenomenological models [12][13][14][15][16], none of these models is built within a formalism guaranteeing a priori all the constraints to be fulfilled at the same time. The same problems also appears from nonperturbative computations of GPDs using various techniques (see e.g. Refs. [17,18] and refs therein), and only in very few cases particular models have taken care of both constraints (see e.g. Refs. [19,20]).
In this paper, we focus on two of them, the so-called polynomiality and positivity properties. The former states that the Mellin Moments of the GPDs are polynomials (of definite degree) of the skewness. It comes from the Lorentz structure of the matrix elements of local operators defining these Mellin moments. On the other hand, positivity gives an upper bound on the absolute value of GPDs taken at a given kinematics in terms of PDFs. This property can be derived from Hilbert space norms. Consequently, covariant approaches usually successfully recover polynomiality but have hard time dealing with positivity. On the other hand, quantum mechanical based techniques, like Fock states expansions, have the positivity property "built-in" but usually fail at getting polynomial Mellin moments. Facing this issue, model builders have generally to favor one of these two fundamental properties, at the risk of violating the other.
We present here a general solution to this problem. Section 2 will remind the reader the details and subtleties of modeling GPDs through the ways mentioned above. Section 3 introduces the theory behind our technique, while in Sect. 4 we present our method, based on the numerical inversion of the Radon transform with incomplete data. Section 5 is then dedicated to the applications of our algorithm to different examples of lightfront wave functions (LFWFs). In Sect. 6 we discuss our results and the phenomenological relevance of ambiguities related to the inverse of the Radon transform. Section 7 concludes this paper.

GPD theory and modeling
GPDs are defined as a lightfront projection of a non-diagonal hadronic matrix element of a bi-local operator. For the sake of simplicity, we will consider only the case of a twist-2 chiral-even quark GPD of a spin-0 hadron (e.g. a pion): where P (resp. ) is the momentum average (resp. transfer) of the hadron states, t = 2 and x (resp. ξ = − + 2 P + ) is the longitudinal momentum fraction average (resp. transfer) of the quarks. q classically stands for the quark flavor. Any four-vector v is expressed in light cone coordinates From Definition (2.1), it is straightforward to realize that one gets back the PDFs in the = 0 limit, and that integrating over x yields the quark contribution to the electromagnetic form factor. On top of these properties, the physical GPD domain obeys t < t min = − 4ξ 2 M 2 H 1−ξ 2 , M H being the mass of the considered hadron, and ξ ∈ [−1, +1]. In this domain, the support is such that x ∈ [−1, 1]. Due to time reversal invariance, the GPDs defined through the operator of Eq. (2.1) are even in ξ . We will use this property to restrict to ξ ≥ 0 in this paper (unless explicitly stated otherwise).
In this section, we remind the reader the important frameworks that allow to fulfill either positivity or polynomiality, and their consequences on modeling.

Lightfront wave functions and positivity
Lightfront quantization allows the expansion of a hadron state |P, λ of momentum P and polarization λ on a Fock basis: (2.2) where the |N , β; k 1 , . . . , k N denote the N -particles partonic states with each particle carrying a momentum k i . β stands for the relevant quantum numbers. These states are weighted by the LFWFs λ N ,β , containing the nonperturbative physics, and normalized as follows: (2. 3) The measure element in Eq. (2.2) fulfills momentum conservation by construction: (2.5) where i labels the partons. Using this Fock state expansion, one can express GPDs in terms of LFWFs [21]. However, the partonic picture and therefore the way the GPDs are related to LFWFs depends on the considered kinematics. In the socalled DGLAP region (|x| ≥ ξ ), the GPD is given by an overlap of LFWFs having the same number of constituents. Keeping the example of the pion, in the region x ≥ ξ , we have [6]: . . ,x a ,k ⊥a , . . . , (2.6) wherex andk ⊥ are the average momentum variables of the LFWF in the GPD symmetric frame, and x a denotes the momentum fraction of the struck quark (labeled here with a for active). The "hat" and "tilde" variables are the corresponding momenta boosted from the incoming and outgoing hadron frames respectively, and can be related to the "bar" variables through: and In the other part of the DGLAP region (x ≤ −ξ ), there is a similar result.
If we restrain ourselves to the valence contribution to the pion, i.e. the first Fock sector N = 2, this relation can be further simplified. Let us consider the π + case in which the first Fock sector would be ud. We get the following GPD: (2.9) where the LFWF depends only on one set of momenta by virtue of Eq. (2.4), and the hat and tilde relations used are those of either the active quark in Eqs. (2.7)-(2.8) or the spectator (by symmetry of the LFWF). This two-body truncated GPD will be used in later sections. Equation (2.6) highlights the underlying Hilbert space structure, as the GPD appears now to be an inner product between two vectors. Using the Cauchy-Schwartz inequality, it is possible to show the above-mentioned positivity property [21][22][23][24], relating the GPDs to the PDFs denoted by the quark flavour q. For instance in the case of the pion quark GPD H : (2.10) which can be generalized to any hadron, although in more complicated forms.
In the other kinematic area, called ERBL (ξ ≥ |x|), almost the same overlap structure appears, but with LFWFs involving different numbers of constituents N and N + 2. Indeed, there is no trace on β nor N in this case.
The fact that the N and N + 2 particles LFWFs overlap in the ERBL region has significant consequences. First, a lowest order truncation in the Fock space would result in a vanishing GPD in the ERBL region. Then whatever the dependencies in x, k ⊥ , x and k ⊥ of the wave function are, once expressed in the kinetic variables of the GPD symmetric frame, non-polynomial dependencies in ξ appear (see Eqs. (2.7) and (2.8)). If we consider also the fact that we do not expect the LFWFs to be polynomials of momenta, it becomes very unlikely that one can obtain polynomial Mellin moments independently in both kinematic regions. Therefore, we expect compensations between the DGLAP and ERBL regions to occur, canceling non-polynomial dependencies in the Mellin moments. However, the fact that the overlap representation mixes different particles numbers in the ERBL region makes these cancellations unlikely at any finite truncation order of the Fock space. In fact, as already stressed in Ref. [6], independent descriptions of the DGLAP and ERBL regions will almost certainly break polynomiality and hence indicate a violation of Lorentz covariance. This suggests that polynomiality strongly ties the DGLAP and ERBL regions in a subtle way since it relates N -body LFWFs to N +2-body LFWFs. To the best of our knowledge, there has been so far no GPD models built from a consistent computation of 2-and 4-body LFWFs. More generally, polynomiality cannot be readily observed in a GPD representation which is manifestly positive. This raises an important issue about GPD phenomenology. As emphasized in Ref. [25], GPD extractions from experimental data face the curse of dimensionality: fitters have to extract (potentially many) functions of several variables depending on unknown parameters to be determined from measurements. A general first principles GPD parametrization, obeying a priori all theoretical requirements, is still unknown at the time of writing. Consequently there is a risk of an unphysical GPD parametrization being numerically favored in data fitting by lack of a priori constraints. Beyond its own intrinsic interest, a general model-building procedure satisfying both polynomiality and positivity is therefore of direct phenomenological relevance. This problem has been solved previously in a particular case [26,27] by building a covariant extension of an overlap of LFWFs from the DGLAP to the ERBL region. We provide in this paper a general solution to this problem. But before describing it, we clarify the meaning of polynomiality in the following section.

Double distributions and polynomiality
Since flavor plays no explicit role in the following discussion of the properties of the GPD H q , the subscript q will be systematically dropped.

Connection to the radon transform
The polynomiality property states that the Mellin moments 1 −1 dx x m H (x, ξ, t) of the GPD H are polynomials of ξ of degree m + 1: This results from Lorentz symmetry and the decomposition of the matrix elements of local twist-2 quark operators P + 2 q(0)γ + (i ← → D + ) m q(0) P − 2 (where ← → D stands for the left-right covariant derivative) in terms of the generalized form factors c (m) k (t). The further restriction to even powers of ξ is brought by the time reversal invariance of QCD. To simplify the following discussion, we will drop the explicit t dependence.

Different representations
More generally, according to the discovery of Teryaev [28] and Tiburzi [29], F D and D do not constitute a unique parameterization for the integral representation (2.20) of a given GPD H (x, ξ). A function χ(β, α), vanishing 3 on the boundary of , can be used for the following redefinition:  [3], while G was later discovered by Polyakov and Weiss [30]. They are a natural solution of the polynomiality constraint. Since then, DDs have been used to model GPDs based on extracted PDFs in the framework of the Radyushkin Double Distribution Ansatz (RDDA) [31]. The Goloskokov-Kroll model [13,32,33] is a popular example of the family of models following this approach. As a consequence of the success of these phenomenological models, DDs have been remembered by most as a convenient way to implement polynomiality in GPD models.
It has become a common misconception to believe that DDs appear only in this restricted subset of GPD parameterizations. On the contrary, our reasoning above shows that DDs are the essence of the polynomiality property. Obeying the polynomiality property is exactly equivalent to being constructed from a DD, and this is a key argument of our approach.
Three main representations -or schemes -for DDs, all of them related to each other by the transformations Eqs. (2.21)-(2.22), have been so far employed: PW: One DD F PW and an extra one-variable function D PW [30] such that: which, precisely, correspond to the DDs F D (β, α) and δ(β) D(α) used above in (2.18). 3 See Ref. [29] for a thorough discussion of boundary conditions. The most general case does not change the main line of our argument but brings some technical complexity.

BMKS:
One function h BMKS [34] such that: (2.25) P: One function h P [35] such that: By an abuse of terminology, h BMKS and h P are often called DDs. Explicit formulas for χ are given in the literature: • to convert a general (F, G) scheme to a PW scheme [28,29], • to convert a general (F, G) scheme to a BMKS scheme [7,29], • and to convert a BMKS scheme to a P scheme [10].
This last formula is given without proof in Ref. [10], and to the best of our knowledge, none has been published. In fact, the connection between the BMKS and P schemes has remained unclear until recently. Since both BMKS and P schemes will play a central role in the following, we add a derivation of the last converting formula in App. B, and discuss it in details.
In a general (F, G) scheme, following Ref. [28], we define the D-term [30] by: The PW scheme has been the one mainly used in phenomenology through the aforementioned RDDA, complemented with various assumptions concerning the D-term. Attempts to model DDs in the BMKS scheme have been undertaken in Ref. [36] following a new regularization procedure [37,38]. It is worth highlighting that none of these schemes can guarantee a priori the positivity property. The additional P scheme has been designed in an attempt to fulfill both polynomiality and positivity at the same time. It has not been used in the building of phenomenological GPD models to the best of our knowledge.

Quark and anti-quark GPDs
We will use here a general notation h for the DDs in all schemes (h BMKS , h P and h PW = F PW , when omitting an extra D-term). Classically denoting θ the Heaviside function, the DD h can be decomposed as follows: (2.29) where h > and h < are called respectively "quark" (with support on > = ∩ {β > 0}), and "anti-quark" (with support on < = ∩ {β < 0}) distributions [6]; and which yield the following "quark" and "anti-quark" GPDs (corresponding to Radyushkin's original GPDs [39]): with support x ∈ [−ξ, +1], and with support on x ∈ [−1, +ξ ], the total GPD being of course H = H > + H < . The factors C > and C < can be respectively: • both equal to 1 in the PW scheme when the GPD follows a degree m polynomiality property (2.11) or when we consider the GPD minus its D-term contribution Eq. In the following, we will consider only quark GPDs (i.e. β > 0), unless explicitly stated otherwise. It should be therefore understood that H stands for H > .

Polynomiality at work: a simple example
We give a practical illustration of the main assertions of the previous subsection, in particular about how the polynomiality condition is implemented by the representation of a GPD as the Radon transform of a DD in different schemes. Adopting an overall normalization for later convenience, we consider the following DD in the BMKS scheme: (2.32) The restrictions H Toy |DGLAP and H Toy |ERBL of the associated GPD: to the DGLAP and ERBL region read: with support x ∈ [−ξ, +1].
The computation of its Mellin moments can be readily performed, and the first 11 of them are presented in Table 2 in App. C. The DGLAP and ERBL contributions to the Mellin moments are rational fractions, not polynomials! One needs to integrate over both DGLAP and ERBL regions to obtain a polynomial in the variable ξ . This is true at all order as can be seen exactly from the expressions of +ξ −ξ dx H Toy (x, ξ) and +1 +ξ dx H Toy (x, ξ) in Eqs. (C1)-(C2). The contribution of the DGLAP region is a rational fraction with apparent double poles at ξ = ±1. Since both the numerator and its derivative vanish when ξ = 1, this fraction actually possesses only a double pole at ξ = −1. The residues can be straightforwardly computed, exhibiting the general structure of both contributions to the Mellin moments: (2.39) We observe that even with a very simple GPD model, the expressions of the DGLAP and ERBL contributions to the Mellin moments of the GPD are intricate; their ξdependencies are not trivially related but add up to yield a polynomial of degree m + 1: Such a delicate cancellation, which suppresses order by order the poles of two rational fractions to produce a polynomial, cannot be the result of a coincidence. It is difficult to imagine that more complex GPD models will not exhibit such a feature.

A necessary and sufficient condition for the BMKS and P schemes
As a final but not less important result of this section, we carefully examine the implications of expressing as a Radon transform the formal link between a GPD and the onecomponent DD schemes BMKS and P. A GPD expressed as a Radon transform of a DD in the BMKS scheme (see Eqs. (2.24), (2.25)) reads: The l-order Mellin moment (l ≥ 0) of the l.h.s. is a l-degree polynomial of ξ , as can be checked easily. Specializing to the case l = 0 4 yields: Similarly, a quark GPD expressed as a Radon transform of a DD in the P scheme (see Eqs. (2.26), (2.27)) writes: 5 The integral of the l.h.s. (its l = 0 Mellin moment) fulfills 5 : where C P is independent of ξ . On the technical side, we assumed in the derivation of Eq. (2.42) and Eq. (2.44) that we could swap the order of integrations in +1 −1 dx dβdα. This is justified only if the integrated function is summable over [−1, +1] × . However the polynomiality condition (2.11) implies the existence of all moments dβdαβ m α n F(β, α) and dβdαβ m α n G(β, α) for n, m ≥ 0 but does not say anything about dβdαh BMKS (β, α) or dβdαh P (β, α) which may even not exist. Both Eqs. (2.42) and (2.44) precisely rely on the finiteness of these last two integrals. This summability assumption, supplemented by the polynomiality condition, makes possible the application of the Ludwig-Helgason consistency condition (A3) to specify a sufficient and necessary condition for H (x, ξ)/x or H (x, ξ)/(1− x) to be in the range of the Radon transform of summable functions or compactly-supported distributions.
From a strict mathematical point of view, there seems to be no reason to prefer one DD scheme or the other. From the physical point of view, we know that GPDs are continuous functions of x, with support in [−1, +1]. According to perturbative QCD [40], one should expect H (x, ξ) to exhibit a smooth quadratic 6 decrease at large x: should be as singular as the GPD H (x, ξ) itself, and we may expect to be finite. On the contrary, PDF extractions indicate that H (x, ξ = 0) has a power-law divergence at small x. Consequently we expect the Regge singularity of the PDF to be aggravated in +1 −1 dx H(x, ξ)/x, which may even not exist for phenomenological valence quark GPD models based on the RDDA. Although we have employed the same notation H (x, ξ) for the GPDs represented either by h BMKS in Eq. (2.41) and satisfying (2.42), or by h P in Eq. (2.43) and obeying (2.44), it is not self-evident that the same GPD H can be represented with both DD schemes. We will illustrate this point with an example in the following section. At last, irrespective of their values, the integrals would still have to be independent of ξ . We may already anticipate the P scheme to be more tractable for computing purposes.

One GPD, several DD schemes: a simple example
We consider again the GPD H Toy of Eqs. (2.34), (2.35). This model is built to be algebraically simple, not phenomenologically realistic, and the considerations above about GPD phenomenology are temporarily left aside. This model is defined in the BMKS scheme by the DD h Toy BMKS (2.32). The condition (2.42) writes: which is indeed independent of ξ . The scheme transform (C3) converts this BMKS representation to the P representation Eq. (C4). The condition (2.44) writes: which has a non trivial ξ -dependence. There is no contradiction with the previous statements, because the underlying DD in the P scheme h Toy P (C5) is not summable over as testified by Eq. (C6). The integrals over all lines of h Toy P are well-defined and finite (their values are obtained from the GPD H Toy ) but h Toy P is too singular to be summable over . We can pick up the singularities (see Eq. (C7)) of this insightful model in the P scheme and get rid of them. Doing so we build a regularized DD h Toy P, Reg and a regularized GPD H Toy, Reg which differs from the original one H Toy only by the introduction of a D-term (see Eqs. (C10), (C11)). This D-term modifies the condition (2.44) by adding an extra contribution: The original value of the sum rule (2.46) supplemented by this contribution yields: which is now independent of ξ . What did we learn? Condition (2.44) did not natively hold because the underlying DD was not summable over . In principle a GPD can be expressed in any DD scheme, but there is no reasons for DDs to be smooth simple functions in every scheme. In our example the DD h Toy P in the P scheme can even not be identified with a compactly-supported distribution. Adding a D-term produced a GPD model obeying condition (2.44) and deriving from a smooth DD, summable over . It is a different GPD model but a model as theoretically consistent as the initial one. It fulfills polynomiality and has the required behavior under discrete symmetries. It has the same forward limit and satisfies the same positivity inequalities. Since the D-term is immaterial to the m = 1 Mellin moment, the form factor sum rules stays the same too. They only differ by the realization of the additional sum rule (2.44) relying on the m = 0 Mellin moment, yielding either (2.46) or (2.48). To summarize, adding this D-term produced a consistent GPD model which is indistinguishable from the original one based on the sole knowledge of the DGLAP region.

Covariant extension to the ERBL region
The way Lorentz covariance binds the DGLAP and ERBL regions together has been questioned for many years. To the best of our knowledge, the first systematic discussion of this link is the work of Müller and Schäfer [41]. Under analyticity assumptions, they argued that an extension of the GPD from the ERBL to the DGLAP region exists and is unique. They also mentioned that the converse does not hold because of the D-term ambiguity. However this study does not solve the problem of actually extending the GPD from one region to the other. Indeed analytically continuing a GPD is a daunting task only if we cannot start with GPD expressions in simple closed forms. In the following, we present an original discussion based on Radon transform properties of covariant extensions of GPDs from the DGLAP to the ERBL region. We then deduce a general procedure to construct these covariant extensions starting from the knowledge of the values of the GPDs in the DGLAP region.

Intuitive picture
We explained in Sect. 2.2.1 that the polynomiality property requires the existence of compactly-supported distributions h BMKS (β, α), h P (β, α), F PW (β, α), and D PW (α) such that: where writing the first two relations is subject to the conditions mentioned in Sect. 2.2.5. The geometrical content of both equations is the same: the l.h.s. is the integral over lines in the plane of the function appearing in the r.h.s.. Going back to Eq. (2.20), one realizes that DDs are integrated in the (β, α)-plane along the lines: which cross the α-axis at x/ξ and the β-axis at x (see Fig. 1). Consider now a real z such that |z| > 1, i.e. (0, z) / ∈ . The families of straight lines joining (0, z) to (x, 0) with x ∈]0, 1] cover the whole domain > . Thus any point in > contributes to DGLAP kinematics. If |z| < 1, all points in the cone of apex (0, z) and base delimited by (0, 0) and (1 − |z|, 0) brings a contribution to the ERBL region. Proceeding with the limit |z| → 1 − we see that all points in > contribute to ERBL kinematics. The same discussion can of course be carried for < . Therefore the entire DD support but the β = 0 line is integrated over when (x, ξ) covers the DGLAP or ERBL region. At this stage the question of the covariant extension of GPDs from the DGLAP to the ERBL region becomes to know whether it is possible to unravel an underlying DD h(β, α) or F PW (β, α) on uniquely from the DGLAP region. Fig. 1 The domains < and > of the DD (resp. DGLAP and ERBL of the GPD) on the left (resp. right). The Radon transform , which is an integration of h on a line parameterized by the couple (x, ξ), is the operation that sends one domain to the other. The goal of the inversion of the Radon transform is to rely only on DGLAP information, meaning that we have only access in DD space to integration lines that cross the α axis on x/ξ > 1 (red lines). In this example, both x and ξ are positive

Formalization
At the GPD level, a D-term δ(β)D(α) is visible only in the ERBL region. At the DD level, a D-term is irrelevant in the geometric construction of the previous paragraph. There is a minimal ambiguity in the extension from the DGLAP to the ERBL region. Is it the only one? In other words, assume there are two GPDs H 1 and H 2 which are equal over the DGLAP region. Are they equal over the ERBL region, up to a D-term? Without loss of generality, we can consider H = H 1 − H 2 by linearity of the Radon transform. This H is zero in the DGLAP region. Since this reasoning holds up to a D-term, we can express H in the PW scheme assuming that D PW = 0. We want to show that F PW = 0. We will rely on the theorem of Boman and Todd Quinto mentioned in Sect. A. Switching to Radon transform canonical variables with Eqs. (2.14), (2.15), the membership to the DGLAP and ERBL regions becomes: The complementary range in [0, 2π ] corresponds to the physical domain of Generalized Distribution Amplitudes (GDAs) [1,42,43] which are the crossed-channel analog of GPDs. In the DGLAP region, the GPD is zero if |x| > 1, which is equivalent to |s| > | cos φ|, which is larger than | sin φ| for |ξ | < 1, i.e. in the GPD physical domain. We therefore assume here that: The aforementioned theorem of Boman and Todd Quinto asserts that: Selecting s = x 0 cos φ 0 , this last condition implies: Thus F PW vanishes on all lines contributing to the DGLAP region, i.e. F PW = 0 on > . Therefore the knowledge of a GPD in the DGLAP region fully determines it over the whole physical domain up to terms inherited from DDs with support on the line β = 0. This is completely consistent with the discussion about the D-term at the beginning of this section.
What is the manifestation at the GPD level of DDs with support on the line β = 0? Assume we started from a function H source defined in the DGLAP region, and that we were able to identify one DD F PW in the PW scheme such that H source = RF PW . Let us further assume that F PW is properly normalized, i.e. that it yields the correct valence quark number (from the PDF) or the correct electric charge (from the form factor at vanishing momentum transfer). If F PW is a function, its values along the line β = 0 do not contribute to line integrals over since this subset has measure 0. In this case there is no remaining freedom. If F PW is not a function, we already saw the D-term contribution. We cannot exclude a priori contributions like δ(β)D + (α) where D + is an even function. For ξ > 0 this modifies H by the addition of a term: If ξ = 0, this new term manifests itself in the PDF through: Such a term would violate quark number conservation if 1 0 dα D + (α) = 0, but there does not seem to be any first principle reason to forbid such terms if this integral does indeed vanish. Such terms have already been considered e.g. in Ref. [39]. Generally speaking, we end up with two families of ambiguities at the DD level: In principle, there could also exist ambiguities involving derivatives of the Dirac distribution, i.e. δ (n) (β)D ± n (α). Line integrals of such terms contribute to the GPD H with 1/(|ξ |ξ n )D (n) n (x/ξ ) up to a factor ξ depending on whether this term is attached to the DD F (with D + n ) or to the DD G (with D − n ). The presence of such terms does not change the argument of our discussion. Their phenomenological relevance will be discussed elsewhere, and we will stick to the δ(β)D ± (α) ambiguities in the following.
Let us summarize: the knowledge of a GPD in the DGLAP region is enough to constrain it in the ERBL region up to additional terms like (for ξ ∈ [−1, +1]): These terms contribute only to the ERBL region. In other words, the inverse problem of reconstructing a GPD from its restriction to the DGLAP region admits infinitely many solutions. Two distinct solutions differ by D + and D − terms as in Eq. (3.12). This statement is independent of the choice of an underlying DD scheme.
The D + and D − terms modify the polynomiality relation (2.11) by acting on the two highest degree coefficients: Since D + and D − have opposite parities, the F-and G-terms do not simultaneously modify the polynomiality relation. To the best of our knowledge, the ambiguity linked to F-terms has not been discussed so far in this context.

Problem reduction
We have so far discussed that the overlap representation, as expressed by Eq. (2.6), makes sure the positivity of the GPD; that there is no simple way to exploit the overlap in the ERBL region and, for the same price, respect the polynomiality condition at any finite truncation order in Fock space; and, finally, that the DD representation ensures this polynomiality. Thus, an immediate and natural approach to model a GPD, by fulfilling both positivity and polynomiality conditions and exploiting the physical information encoded in the LFWFs would result from: 1. the computation of the overlap DGLAP GPD, symmetric in Fock space (overlap of LFWFs with the same number of constituents), 2. the derivation of a DD from this DGLAP GPD by an inverse problem, 3. the extension of the GPD to its full kinematic domain by means of the DD representation.
This is the program we will apply in the following.
In particular, given a GPD H (x, ξ) with support x ∈ [−1, +1] and with available non-trivial information only in the DGLAP region 0 ≤ ξ ≤ |x| ≤ 1, our goal is to find a DD h such that: (3.14) where the factors C > and C < were defined in Sect. 2.2.3.
In the next section, we describe a well-established numerical procedure to do it. But before, we make the following remarks: • The quark and anti-quark GPDs are not correlated in the DGLAP region. In positive  Only the salmon red domain is used. The blue one is deduced by parity. And the white one, i.e. < , is not correlated in DGLAP, and therefore can be dealt with separately These two properties together reduce the size of a numerical problem by 4 by comparison of a direct numerical inversion of Eq. (3.14). Indeed, we can separate a general GPD into two distinct problems H > and H < , by virtue of linearity of Eq. (3.14) and the non-correlation in the DGLAP region. This limits us to half the DD domain (either > or < ) without loss of generality. And by parity, we reduce again the problem by half. Figure 2 summarizes this. This significantly decreases the computing cost of the numerical inversion 7 and further constrains the target solution.
As stated in Sect. 2.2.3, we will discuss only the case of quark DDs and GPDs. The treatment of anti-quark DDs and GPDs is essentially the same.

Numerical implementation
This section presents the considered numerical implementation of the covariant extension of a GPD from the DGLAP to the ERBL region with its challenges and results. In Sect. 3 we explained that we can solve our physics problem by inverting the Radon transform. This may seem straightforward since the Radon transform is linear, but this task is in fact much more difficult. Indeed the inverse Radon transform may not be continuous (see App. A) and, in a loose sense, two "close" GPDs may be obtained as Radon transforms of very "different" DDs. Since we are facing an incomplete data problem (we know the GPD only in the DGLAP region), the sensitivity to noise is expected to be even stronger than in the complete data problem, where we search a DD from the complete knowledge of a GPD and a GDA over their whole kinematic domains. In this respect, we note that reconstruction artifacts have already been reported [28,44] for the latter situation.
One key remark is in order here. We do not know any closed-form formula for the inverse of the Radon transform restricted to the DGLAP region. It is not even clear that such a 7 Given an algorithm with polynomial complexity O (N p ) where N is the size of the problem, solving two equal-size independent problems would have a O (2 N p ) complexity, which is much better than a joint problem of complexity O (2N ) p . formula even exists. However we do not need it, and it would be of limited practical interest: the potential amplification of noise is related to the discontinuous nature of the inverse Radon transform. It is not the manifestation of a poor numerical scheme or of a badly-designed computing code. It is the inescapable consequence of a precise and general mathematical statement. Even if we had at our disposal a closed-form expression of the inverse Radon transform, we should expect this phenomenon of noise amplification except in the lucky but rare situations where all computations can be performed analytically. As soon as approximations enter the game, the discontinuous nature of the inverse Radon transform may generate some artifacts in the sought-after DD.
The way to go is well-known in the mathematical literature (see e.g. Ref. [45]). Assuming that the underlying DD is smooth enough, it is possible to numerically invert the Radon transform while keeping noise under control. This is called regularization.
This sections falls into three parts. We first discretize our problem to reduce it to the computation of the pseudo inverse of a rectangular matrix. Then we select adequate linear solver and regularization procedure. At last we validate our computing chain with simple but relevant test case scenarios.

Discretization
The goal now is to obtain a discrete problem from the integral equation (2.30), and we will use the usual notation: where A is a m × n matrix, X a vector of dimension n, and B a vector of dimension m.

Mesh
To obtain this finite-dimensional linear problem, the DD space should first be discretized. In an abstract way, we use a set of basis functions v j for the decomposition: where the index j labels the set of basis functions, and therefore the degrees of freedom. Adopting a formalism close to the one of finite element methods (FEM) [46], these basis functions are in one-to-one correspondence to given nodes in the DD domain. Indexing these nodes means indexing the basis functions. Applying this to a given mesh which is a set of vertices (or corners) and edges defining its elements, we can be more explicit. A basis function is non-zero only on elements adjacent to its corresponding node, and the restriction of a basis function to one such element is the Lagrange interpolation with respect to this node, i.e. the polynomial that is equal to 1 on the said node, and 0 on all others. See Fig. 3 for an example of such a basis function.
Following the conventions of FEM [47], we will consider the following classification: P n -Lagrange Used for triangular meshes, where the restriction of a basis function to a triangular element is an interpolating polynomial of total degree at most n. For example, for P 1 , it would be a polynomial of the form a + b β + c α. Q n -Lagrange Used for quadrilateral meshes, where the restriction of the basis function to a mesh element is an interpolating polynomial of partial degree at most n. For example, for Q 1 , it would be a polynomial of the form In the case of linear piece-wise functions (P 1 or Q 1 ), the considered nodes of the basis functions are the vertices of the mesh. For higher orders, the nodes also include other points (such as the middle of the edges for P 2 and Q 2 ). We will also consider constant piece-wise functions and we will call those elements P 0 (which corresponds to d P 0 in FEM notations). In this case, each basis function corresponds to one element (or any node in the interior of the element, e.g. the center of gravity, to keep the same correspondence between nodes and basis functions). Table 1 summarizes this.
Our unknowns h j of Eq. (4.2) correspond to the values of the DD h on the nodes j, and will be recast into the vector X of the discrete problem (4.1).
For the work presented here, we will consider only a triangular mesh, since the domain is a triangle anyway (see Fig. 2), with P 1 or P 0 elements.
We will always use the index j in the following to label the basis functions, i.e. the nodes, and the index k for labelling the triangular elements. Of course, in the case of P 0 , the indices will be interchangeable, since a basis function is defined by an element.

Basis functions
In the case of a triangular mesh, it is natural to use barycentric coordinates to define the basis functions (instead of the Cartesian coordinates). For a given triangle k, we will denote by {λ 1 k (β, α), λ 2 k (β, α), λ 3 k (β, α)} the barycentric coordinates with respect to the three vertices. Note that the number of degrees of freedom is still 2, since λ 1 k + λ 2 k + λ 3 k = 1. A given point (β, α) belongs to the triangle k if all three barycentric coordinates are positive. Moreover, for P 1 elements, they provide natural restrictions for the basis functions, since they are exactly the linear Lagrange interpolations at the vertices. If we denote by (β i , α i ), i = 1 . . . 3, the three vertices of a triangle, then the barycentric coefficient with respect to the first vertex can be written as: and the others similarly by cycling indices. The matrix of the linear problem is determined by the linear operator that transforms a DD into a GPD in Eq. (2.30).
To build this matrix, we only need to know the Radon transform 8 of a basis function: Let us first express this basis function in the P 0 and P 1 cases (superscripts 0 and 1 respectively): wherej is the vertex j recast to the limited set {1, 2, 3} of vertices of the element k. The P 1 basis function is also represented in Fig. 3. Applying the Radon transform on these basis functions yields: where the bounds of the integration α k min , α k max are determined with the three inequalities given by the Heaviside functions (positive barycentric coordinates). For higher order elements, the idea is the same, only the integrated function will change.

Sampling
The next step is then to discretize the GPD variables (x, ξ), i.e. to sample the set of straight lines intersecting the domain . Given that we have only access to DGLAP kinematics, we will use the couples (x, y) ∈ [−1, +1] 2 with y = ξ/x. The choice of (x, y) will determine a line of the matrix. More precisely, the matrix A will have the coefficients: where 1 ≤ i ≤ m indexes the lines of the matrix, and 1 ≤ j ≤ n indexes the columns, i.e. the nodes in DD space. The C > factor was introduced in Eq. (2.30). The size of the matrix is chosen such that we maximize the information, i.e. we need to integrate over lines that cross all the elements of the DD mesh. A value of m ∼ 4 n is empirically satisfying. The matrix can be therefore built by picking random couples (x, y) until we attain the desired size. The results will of course depend on the matrix used and it is interesting to consider this as a source of "statistical error", whereas the regularization procedure (see the following section) would be the source of "systematic error". The statistical error can be managed quite easily and reduced considerably by picking as many samples as we want, whereas the systematic error remains a challenge to estimate.
Once the matrix is built, we use the set of chosen couples (x, y) to build the vector B r.h.s. of Eq. (4.1) with simply: (4.10) In summary, A is a matrix m × n where n is the number of mesh elements for P 0 (or number of vertices for P 1 ) and m the number of straight lines intersecting . Each line will typically cross O( √ n) mesh elements, which means that only O( √ n) coefficients on a matrix line are non-zero and A is a sparse matrix. We need more constraints than parameters (m > n) and we usually use m = 4n, making the rank of A n (i.e. close to full-rank). B is a vector of dimension m, and X of dimension n.

Linear solver and regularization
An additional complexity arises in the selection of the matrix inversion routine. In Sect. 3.2 we assumed the existence of a covariant extension of a DGLAP-restricted GPD H |DGLAP and showed its uniqueness up to the manifestations of ambiguities on the line β = 0. The key question is to know whether such an extension ever exist? Or, stated differently, given a function defined in the DGLAP region (a putative overlap of LFWFs), is it possible to extend it to the ERBL region in a way satisfying polynomiality? Existing criterions as in Eqs. (2.42) and (2.44), complemented by the Ludwig-Helgason consistency conditions, deal with the GPD known over its whole physical domain, not its restriction to the DGLAP region. A numerical solver may have to handle a linear system as in Eq. (4.1) but without one and only one solution. This is common in computerized tomography, not because the solution does not exist (there was one object inserted inside the scanning device), but because the experimental signal comes with noise which may apparently modify the original situation to an inconsistent data problem. In the framework of the Radon transform, causes may be multiple: the integration lines may not cross the same domain (no solutions), or they may be parallel and close to one another and bring redundant information (infinitely many solutions). One efficient way to ensure that the solution always exists and is unique is to turn to a least-square formulation: Search X ∈ R n such that ||AX − B|| 2 is minimum, (4.11) where ||.|| generically denotes a norm in a finite-dimensional vector space. In the present work, we use a recent iterative conjugategradient type algorithm for sparse least-squares problems: LSMR [48]. For inconsistent problems (where the leastsquare formulation is favored), it is equivalent to a minimum residual algorithm for the problem: but it can also solve directly the problem (4.1) when it is consistent, i.e. when the numerical approximation of the target solution is equal to the exact function. In the P 0 case, it means functions that are already piece-wise constant, whereas a P 1 approximation can reproduce exactly a (piece-wise) linear polynomial. This type of algorithms applies naturally its own regularization process, with the number of iterations being the regularization parameter. To illustrate this, we can use the so-called L-curve [49], which is a curve following a regularization parameter (which is in our case the number of iterations) and shows the compromise between the norm of the solution X (the larger the norm, the larger the impact of noise) and the residual norm r , where r = AX − B (which we desire to be small enough to converge to the real solution). This procedure gives the optimal regularization factor to choose for each problem, as the point of maximum curvature of the "L", as shown in Fig. 4.
In practice, as illustrated on Fig. 4, it is very difficult to determine this optimal regularization parameter for the considered problems. A better way to stop the iterations is to consider the stopping criteria used by these algorithms such as LSMR: • For a consistent problem: r ≤ atol A X + btol B ; • For a least-squares problem: A r ≤ atol A r , where atol and btol are the input tolerances.
An empirical value of 10 −5 for the least-square tolerances gives a good compromise between noise and convergence, and this has the benefit of being valid for all considered cases, assuming that the considered DD is smooth and not one of the problematic cases with singularities such as Sect. 5.3 (for which a specific workaround is needed and explained therein). A 10 −5 value for atol (resp. btol) means that the matrix A (resp. the right hand side B) is known exactly up to the fifth decimal, while the rest is numerical noise. Of course, in practice, we can compute analytically A (if the chosen basis functions and mesh allow us to compute the Radon Transform without numerical integration, as it is the case for the method presented here) and B (if the GPD is known analytically), and they are therefore known exactly, up to machine precision. But in the inconsistent (i.e. leastsquares) case, the considered vector B is different from the one due to a GPD calculated from the discrete numerical DD. This difference is the finite limit of the residual, in contrast with the consistent case where the residual has a zero limit. Even though small, when allied with the ill-posed character of the inversion, it can have a large impact on the solution. This is why we consider in practice that B (or equivalently A) is not known exactly and neglect higher decimals ; we apply a regularization procedure by doing so.

Test and validation of the numerics
The first immediate check we can perform to validate the numerical implementation described above, consists in the following. We first take a simple Ansatz for the DD, irrespective of the considered value of t (e.g. t = 0). We then compute the associated analytical GPD by applying the Radon transform, and use only its DGLAP part to apply our numerical inversion and obtain a numerical estimate of the DD. Finally, we compare this result with the original Ansatz. We will apply this testing procedure to the following three quark GPDs, each one of them deriving from a DD in the P scheme: 14) where q(β) is the associated PDF. In particular, we specialize for the case N = 1 and take 9 q(β) = 30 β 2 (1 − β) 2 . We thus obtain a closed algebraic formula for the GPD which, in turn, can be numerically inverted with the procedure described above.
In all cases, we adopted the P scheme as it appears to be a convenient representation for GPD models derived through a covariant extension of an overlap of LFWFs, as discussed in Sect. 2. Therefore it will systematically be employed for 9 This is a very good practical approximation of the result for the valence-quark pion's PDF obtained in Refs. [50,51], within a Bethe-Salpeter and Dyson-Schwinger approach and by including the appropriate correction to the impulse-approximation. It also results directly from the overlap of the LFWF derived from the same Bethe-Salpeter wave function [17]. the analysis of the more practical examples scrutinized in the next section. In the leftmost plots of Fig. 5, we display the comparison of the exact algebraic and numerically approximated h cst P (β, α) for the case (i), while the rightmost ones stand for the comparison of the GPDs directly obtained from both DDs. The upper and lower plots have been obtained, respectively, with piecewise constant (i.e. P 0 elements) and piecewise linear (i.e. P 1 ) basis functions when discretizing the DD space for the reduction to a finite-dimensional problem, as explained in Sect. 4.1. Analogous plots, and similarly arranged, are displayed in Fig. 6 for case (ii) and Fig. 7 for case (iii).
The mesh was generated using the Triangle software [52] with a requirement of maximal area for the triangular elements equal to 0.001, which produced a mesh of 427 vertices and 780 elements. The linear solver used is described in Sect. 4.2. As justified before, we used a tolerance of 10 −5 as a regularization procedure.
All in all, as an overall conclusion we can assert that the numerical inversion approximates very well the three known GPDs as they appear not to differ significantly in all the cases. Several important points should be however stressed out: • The numerical inversion relies only on the knowledge of the GPD within the DGLAP region and its extension to the ERBL is our main goal, wherefore the examination of algebraic and numeric GPDs over the ERBL region is the main outcome of Figs. 5, 6 and 7. • The numerical reconstruction of the DDs may seem quite noisy or far off in some cases, but these discrepancies do not hinder the reconstruction of the GPD, for which the convolution helps smooth these defaults. 10 The physical object of interest is the GPD, not the DD, therefore these discrepancies are not an issue. • It should be noted however that the constant DD can be reconstructed numerically exactly (up to machine precision). Indeed, the regularization in that case is not needed, since there is no distinction between the analytical DD and its discretized version. We could therefore directly invert the discrete problem and recover the exact DD (which would be equivalent to using a tolerance of 0 instead of 10 −5 ), but for the sake of homogeneity, we decided to employ the same method for all shown examples. • It may seem from the plots that the P 1 discretization does not improve on the P 0 result, but this is not true. We chose here to use the same mesh for both P 0 and P 1 elements. This particular mesh has 427 vertices and 780 triangular elements. Hence, the number of degrees of freedom for P 1 (i.e. the number of vertices) is half the one for P 0 (i.e. number of elements). In other words, we attain with P 1 a similar result to the P 0 one but with half the degrees of freedom, i.e. at a much lower cost. It is therefore a significant improvement. In Sect. 5, we will keep

Examples of applications
The approach described in Sects. 2 and 3 can be either applied to numerous existing LFWF-based GPD models to covariantly extend them from the DGLAP region to the ERBL one, or used to build a covariant GPD model, reliable on both DGLAP and ERBL regions, from the knowledge of the LFWF. Although in some particular cases, an analytical derivation of the DD is possible and a full GPD in both DGLAP and ERBL regions can be thereupon obtained, one can only proceed systematically by applying the numerical technique that has been introduced in Sect. 4. In this section, aiming to illustrate the procedure without the intention of being exhaustive, we provide four examples of GPD mod-els, three of which can be extended to the ERBL region both analytically and numerically, allowing us to benchmark our algorithm.

Algebraic Bethe-Salpeter model
We consider first a specific pion GPD model described in Ref. [17], based on the following helicity-0: and helicity-1: contributions to the LFWF, obtained by integrating and properly projecting the pion Bethe-Salpeter wave function resulting from the algebraic model described in [53] where M is a model mass parameter introduced at the level of the quark propagator. A value of M ∼ 0.318 GeV allows to recover the pion charge radius [54]. Equations (5.1), (5.2) correspond to specializing the helicity-0 and helicity-1 components given, respectively, by Eqs. (154) and (155) in Ref. [17] for the asymptotic case ν = 1 therein. Then, exploiting the GPD overlap representation as described in Sect. 2.1 and extending Eq. (2.9) with the two contributions: one is left with 11 : as a fully algebraic result for the DGLAP region, where: encodes the correlated dependence of the kinematical variables x and t. A careful computation allows for the derivation of the following closed expression for the DD in the P scheme given by Eqs. (2.26), (2.27): Thus, we compare the results given by Eqs. (5.4) and (5.7), respectively derived for the DGLAP and ERBL kinematics, with those obtained by the numerical inversion of the linear problem described in Sect. 4; and display the outcome in the leftmost panels of Fig. 8. As can be seen, both numerical and algebraic results compare strikingly well, not only at a qualitative but also a quantitative level. Albeit the successful comparison we obtain for the case t = 0 is satisfactorily enough, aiming at the practical validation of the numerical technique at other values of t, we have also displayed in the rightmost panel of Fig. 8

Algebraic spectator model
If LFWFs have been widely used in attempts to model the pion, the case of the nucleon has also been treated previously, in particular in the pioneering paper of Hwang and Müller [26]. They have developed an algebraic parametrization for two-body LFWFs of the nucleon (described by a constituent quark and a spectator scalar diquark):

Fig. 9
Comparison for the GPD E between our numerical algorithm and the algebraic result of Ref. [26]. Same conventions as in Fig. 8 ϕ with M, m and λ being respectively the nucleon, quark and spectator masses, p being a free parameter. The authors have shown that with such a model, after calculating the overlap of wave functions (see Eqs. (13)(14) of Ref. [26]), one can write the GPD E in the P scheme: (5.11) and they have been able to extract analytically: where N is a constant determined by the usual PDF normalization (obtained with the GPD H, and not E). Therefore, the authors managed to extend their specific model in the ERBL region. Consequently, we use this model as an additional benchmark for our numerical technique. The comparison between our numerical reconstruction and the algebraic result is shown on Fig. 9. We use the same parameters values as in Ref. [26], i.e. M = 1 GeV, m = 0.45 GeV, λ = 0.75 GeV and p = 1. In this case, the normalization constant for the DD is N ∼ 0.176. We stress that in Fig. 9, the qualifying terms "analytical" and "numerical" refer to the DDs. In other words, "analytical result" means that the GPD is calculated through an integration of the "analytical" DD (5.12) in the ERBL region (or through the integration of Eq. (15) in Ref. [26] in the DGLAP region), whereas "numerical result" means that the GPD is calculated from the "numerically" reconstructed DD (from DGLAP information only), through e.g. Eqs. (4.2), (4.8) and (5.11).

Parametrization with Regge behabior
All the previous examples dealt with GPDs smoothly behaving at x = 0. However, phenomenological models for valence quark GPDs often exhibit an integrable singularity, typically a 1/ √ x behavior. If such a GPD is related to a DD h P in the P scheme, then it can be shown that H (x, 0) ∼ +1 −1 dα h P (x, α) at small x. Therefore h P itself may also exhibit an integrable singularity 1/ √ β. The numerical method presented up to now approximates the target DD by piecewise constant or piecewise linear functions on . In particular all approximations are bounded, even if in principle they can be more and more peaked when the mesh gets thinner. As discussed in Sect. 4, careful choices of the size of the mesh and of the number of iterations are essential for the resolution of the inverse problem with adequate control of the numerical noise. With a naive discretization, it is difficult when dealing with singularities to make sense of a solution ; it would probably require a number of iterations that is not attainable due to the necessary truncation of the regularization. We thus adopt a more educated discretization ; knowing that we are dealing with Regge-type singularities, we can adapt our method accordingly, by discretizing: which will be less singular that h P , possibly even free of singularities. This change of target function only modifies the kernel of Eq. (2.30) (it is not a Radon transform anymore), otherwise everything readily follows the same procedure.
Let us exemplify this technique on a simple parameterization of nucleon DD for which, once again, the expression is  Fig. 8 (left panel) already known. We use the RDDA as presented in Sect. 4.3, Eq. (4.14) for N = 1, but this time with a singular PDF 12 : which would give the following DD (with a 1/ √ β singularity) on > : The comparison between the algebraic GPD and the result obtained through the numerical reconstruction of the DD with DGLAP information is shown on Fig. 10.

Gaussian wave functions
Last, on top of the Bethe-Salpeter wave function already mentioned above in the pion case, one can compute a LFWF with a Gaussian Ansatz (such Ansätze have been used for instance in ADS/QCD computations, see Ref. [55] and refs therein). We choose to work with the following one: This wave function yields in the DGLAP region (using Eq. (2.9)): To the best of our knowledge, an algebraic expression for the associated DD (it it exists) has not been published yet. We anyhow examine this case for illustrative purposes and display in Fig. 11 the results for the GPD in both DGLAP and ERBL derived from the DDs obtained with our numerical inversion. Indeed, some qualitative similarities in the shape can be noticed when compared to the results previously obtained with the Bethe-Salpeter LFWF and depicted in Fig. 8. The parameter value of M ∼ 0.315 GeV has been used here.

Chiral symmetry and the soft pion theorem
Pion GPDs should satisfy one extra first principle constraint, tied to the specific role of the pion with respect to chiral symmetry. The isosinglet and isovector pion GPDs H I =0 and H I =1 are defined in terms of the following matrix elements: where = (ψ u , ψ d ) denotes the doublet of u and d quark fields, I is the identity, τ c the Pauli matrices, τ ± = (τ 1 ± iτ 2 )/ √ 2, and π 1 , π 2 and π 3 is a basis of the adjoint representation of the Lie algebra su (2). Expressing theses vectors in terms of charge eigenstates: 3) yields in particular: The additional charge conjugation constraint, H q π + (x, ξ, t) = −H q π − (−x, ξ, t) for q = u, d dictates: Therefore the information for the whole system of pion GPDs can be equivalently encoded in H I =0 and H I =1 , or H u π + and H d π + , which makes H I =0 (resp. H I =1 ) an odd (resp. even) function of x: x, ξ, t) for I = 0, 1 . (6.9) In the chiral limit, the soft pion theorem [56] states that: where ϕ π is the leading-twist pion DA. This theorem holds in any framework where chiral symmetry is properly implemented, as clearly exemplified in a recent independent deriva-tion in the Dyson-Schwinger-Bethe-Salpeter framework [51]. We saw in Sect. 3.2 that an overlap of LFWFs in the DGLAP region characterizes a set of GPDs of the form: H (x, ξ) , (6.12) with δ H defined in Eq. (3.12) and H LFWF (x, ξ) obtained from a numerical inversion of the Radon transform in a particular DD scheme. By linearity, this relation holds in a flavor basis or in an isospin basis. In particular, the x-parity of D + and D − induces two different contributions to the isoscalar and isovector GPDs, and hence to the soft pion theorem: 14) The requirement dz D + (z) = 0 is manifest in this last equation. The normalization of the LFWFs used in the overlap leading to H LFWF imposes the equality of the form factor at vanishing momentum transfer and of the integral of the leading twist DA. Therefore the remaining freedom in the covariant extension of a GPD from the DGLAP to the ERBL region is flexible enough to satisfy the soft pion theorem. The latter is explicitly shown by a more detailed analysis of the algebraic Bethe-Salpeter model in Ref. [57]. We already emphasized that, for any hadron, our model building strategy systematically produces GPDs satisfying a priori polynomiality, positivity, and the correct consequences of discrete symmetries. We have just shown that, in the pion case, the soft pion theorem can always be fulfilled too. Whatever the input overlap is, we can construct a model obeying a priori all the first principles constraints. The reduction formulas to obtain PDFs or form factors impose further restrictions on the choice of the LFWFs -which can be read from the overlap in the DGLAP region -and not on the covariant extension in itself.
At last, to the best of our knowledge, nothing prevents the soft pion theorem to be satisfied in a model with an inadequate chiral symmetry behavior. Fulfilling the soft pion theorem may not be equivalent to enjoying a proper implementation of chiral symmetry breaking. However it is tantalizing to think that the soft pion theorem can only be met if the underlying LFWFs do possess a correct chiral behavior, in which case we may speculate about the physics content of the D − and D + terms restoring this chiral property. These considerations go beyond the scope of the present paper, and are left to a later study. 3), we may expect that different DD schemes will yield different GPD covariant extensions starting from the same GPD in the DGLAP region. This may give a feeling of arbitrariness of the result. On the contrary we will show below that the set of solutions of the inverse problem does not change. Whatever the choice of the DD scheme underlying the inverse of the Radon transform, given an input function in the DGLAP region, the set of GPDs whose restriction to the DGLAP region is this input function stays the same.
To prove this, we only need to study the influence of a DD scheme transformation on the set of solutions described in Sect. 3.2 in the PW scheme, and in particular exhibit the manifestation of the D + and D − ambiguities in the BMKS and P schemes. Since DD scheme transformations are linear, we only have to adapt the general expressions given in App. B to the particular case: Using Eq. (B18), we see that the D + and D − ambiguities in the PW scheme bring a contribution δh P to the DD h P we would obtain by inverting the Radon transform in the P scheme: . The corresponding exercise in the BMKS scheme yields: (6.18) Direct computations of the associated GPD contribution δ H establish that we actually recover the expression (3.12) derived in the PW scheme. The manifestation of the ambiguities related to D + and D − nevertheless presents different forms in different DD schemes.
We will now explain with the simple example H Toy of Eqs. (2.34), (2.35) the impact of the different DD schemes we have considered so far. This GPD model was defined in the BMKS scheme from Eq. (2.32). As mentioned in Sect. 2.2.6 it satisfies the necessary and sufficient condition (2.42) expressing its membership to the range of the Radon transform of summable functions in the BMKS scheme. We can express this GPD in the P scheme by means of Eq. (B19) to obtain a DD h Toy P which is not summable over and violates the analogous necessary and sufficient condition (2.44) in the P scheme.
At the same time, we can numerically invert the Radon transform in the P scheme with the restriction H Toy |DGLAP to the DGLAP region as an input function. The numerical procedure described in Sect. 4 will approximate the DD in the P scheme by a piecewise linear function, and will naturally pick up a (reasonably) smooth solution among all the solutions of this inverse problem. This is indeed what has been successfully checked and showed in Fig. 6 with the second test case of Sect. 4.3. The DD h Toy P, Reg (C9) in the P scheme is smooth (it is a polynomial over > ), and is associated to a GPD equal in the DGLAP region to H where: A direct inspection of Eq. (6.17) reveals that, in the P scheme, the DDs h Toy P and h Toy P, Reg , respectively underlying H Toy and H Toy, Reg , should differ by a term generated byD = D − (and D + = 0): 13 We add the extra sgn function because all computations in App. C were performed under the assumption ξ > 0. It is the presence of sgn(ξ ) or |ξ | that distinguishes the D − and D + contributions since ξ -parity reflects the time invariance of QCD.
which is exactly h Toy P, Sing (C7), the part of h Toy P which had to be subtracted to yield the summable DD h Toy P, Reg in the P scheme. We are therefore left with a clear picture of the inversion procedure. The input function H |DGLAP does not know anything about DD schemes. In principle we may choose any DD scheme to perform the numerical inversion. However, as pointed out here and in Sect. 2.2.6, a GPD may be related to a smooth DD in one DD scheme, and to a singular one in another scheme. There is a difference between DD schemes when actually inverting the Radon transform, and the algorithm of Sect. 4 will naturally select one of the smoothest DDs among the set of solutions of the inverse problem. Consequently, the DD and the GPD provided by the inversion procedure in one DD scheme have no reason to be related by a simple DD scheme transformation to the DD and GPD obtained from the inversion procedure in another scheme. This is also evident when looking at the necessary and sufficient conditions Eqs. (2.42) and (2.44) of Sect. 2.2.5 which seem barely compatible. Thus for a fixed GPD model, we cannot generically expect h P and h BMKS to be both reasonably smooth. However, by properly keeping track of the D − and D + ambiguities in various DD schemes, we observe that the set of solutions to the inverse problem does not depend on the choice of the DD scheme, as it should be.
As a side remark, let us insist on the non-trivial structure of δh P in Eqs. (6.17) and (6.21). A contribution living on the β = 0 line in the PW scheme spans the whole domain > in the P scheme; the singularity structure in DDs may come under different guises when changing schemes: a δ(β) in the BMKS scheme manifests itself in the P scheme as a term potentially singular when β is close to 1, and conversely. Forgetting that h P may not be summable, or even a compactly-supported distribution, a direct application of the support theorem of Boman and Todd Quinto quoted in App. A would have totally missed the D − contribution. The difficulty arises from the fact that physics dictates hypothesis on the DDs F and G, not on the DDs h P or h BMKS . Regarding this, the PW scheme is a natural scheme to lead the discussion on the inversion of the Radon transform, and the propagation of the remaining freedom to fully characterize the set of solutions. At last, the independence of the set of solutions to the inverse problem on the choice of the DD scheme implies in particular that we can select one DD scheme without loss of generality to invert the Radon transform.
We have been working with two one-component DD schemes (P and BMKS) and by an extension of the analysis of App. C, we may expect that some other one-component DD schemes exist. However, we already advocated for the advantages of the P scheme and the drawbacks of the BMKS scheme from a practical point of view. We can exemplify the difference further by considering the algebraic model of Sect. 5.1, whose DD h p is given in Eq. (5.6). At t = 0, h p is simply a quadratic polynomial, and therefore can be easily handled by a computer. Since we know the results algebraically, we can use the scheme change formulas for the BMKS (see App. B) and PW schemes (see Eq. (C12) or Ref. [29]). The scheme transform yields the following results: (6.22) in the BMKS scheme and: in the PW scheme. Several comments are in order. First, trading the P scheme for the BMKS one introduces singularities, making it significantly harder to be extracted numerically. If not singular, the PW remains nevertheless more complicated. A similar statement holds with our toy model H Toy in the PW scheme Eqs. (C14), (C15). On top of this, the P scheme allows to extract a D-term contribution from the numerical inversion, while an inversion in the PW scheme will structurally miss a D-term and yield an incomplete polynomiality. From the model-building perspective (and apart from the pion case and the soft pion theorem), a GPD model directly obtained from the numerical inversion of the Radon transform in the P scheme satisfies all required first principle constraints. The need for D − and D + contributions may be dictated by phenomenology but not by first principles. Contrarily, a GPD model similarly produced as an output of the same algorithm but in the PW scheme will have to be complemented by a D − contribution. For all the reasons above, we think that there is a marked advantage in proceeding with the inversion of the Radon transform in the P scheme.

Conclusion
GPDs computations, either from modeling or ab-initio techniques, remain a hot topic today in hadron physics, and the experiments at COMPASS, Jefferson Lab and on an EIC will certainly have a deep impact on our understanding of hadron structure. However, to fully exploit forthcoming experiments, the community needs models and ab-initio techniques relying on a firm theoretical basis.
For the first time a systematic way to build GPD models such that all properties are a priori fulfilled is outlined. Using the method and algorithm presented here, all first-principle theoretical constraints are met: not only the support property, but also both positivity and polynomiality properties, the soft pion theorem in the case of pion GPDs, and, of course, the remaining, easier-to-implement constraints.
We showed that using of the Radon transform allows, from the overlap representation of GPDs, to covariantly extend results obtained in the DGLAP region to the ERBL one. We emphasized the key role of the different DD schemes and how the choice of the latter impacts numerical reconstructions. As a consequence, even knowing only a GPD in the DGLAP region, it is possible to obtain an extension to the ERBL region including a non-trivial D-term, irrespective of potential additional "D-term-like" contributions. Such extra terms may be used for phenomenological purposes, but are in principle not required as soon as the polynomiality property is obeyed up to its highest degree.
We described the problems tied to the discontinuous nature of the inverse Radon transform. Our numerical results are in very good agreement with algebraic evaluations when available, validating the method and giving confidence to its use in the vast majority of cases where no covariant extension of LFWF overlap is known. These robustness and flexibility emphasize the strength of our technique.
In terms of interpretation however, one should keep in mind that if the Fock state truncation is manifest in the DGLAP region, it is no more visible in the ERBL one. The results we obtain in the ERBL region correspond to the contribution needed to insure that, together with the DGLAP parts, the GPD fulfills polynomiality. The number of Fock states implicitly used in the ERBL domain for that covariant completion is so far unknown. This lack of interpretation in the ERBL region is compensated by the fact that LFWF model builders and ab-initio specialists will be able to extract information from this region, and in fine to compare to experimental data.
At last, on a longer time scale, the present work opens the path to a phenomenology based on LFWFs, using observables related not only to GPDs, but also to TMDs, PDFs, DAs or even Wigner distributions. and "Non-Perturbative QCD 2016", Sevilla, Spain, where first phenomenologically relevant applications have been presented. This work is partly supported by U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract no. DE-AC02-06CH11357, by the Commissariat à? l'Energie Atomique et aux Energies Alternatives, the GDR QCD "Chromodynamique Quantique", the ANR-12-MONU-0008-01 "PARTONS" and the Spanish ministry Research Project FPA2014-53631-C2-2-P.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .
In particular, a DD in a general (F, G)-scheme cannot be cast in the P scheme by the transformations Eqs. (2.21), (2.22) if it does not vanish on the border of .
The extra function C in Eq. (B17) brings a contribution h P,C to the Pobylitsa DD h P through Eq. (B4): where C denotes the first derivative of the function of one variable C. h P,C (β, α) vanishes for β < 0 if C is zero, and C is itself null (since it is an odd function). The only function h P such that Eq. (B1) holds is thus expressed in Eq. (B18). Following similar steps, it is easy to derive a transformation from a (F, G)-scheme to the BMKS scheme. We look for an α-odd function χ BMKS (β, α), vanishing on the border of , such that: G(β, α) − ∂χ BMKS ∂β (β, α) = αh BMKS (β, α) . (B24) This system can be solved for χ BMKS : which can also be found in the literature, e.g. in Ref. [7].

Appendix C: A simple example exhibiting the polynomiality property
We highlight here some subtleties related to the various DD schemes. To apply the discussion of Sect. B leading to Eq. (B17), we need a DD vanishing on the boundary of . We simplify further by considering a quark GPD (hence the related DD vanishes for negative β) and a DD expressed as polynomials to simplify the computations of Mellin moments or GPDs. We scrutinize the simple GPD model of Eqs. (2.34), (2.35), based on the DD (2.32) to show what can be expected in a simple example. Table 2 displays the first eleven Mellin moments of H , as well as the contributions to the Mellin moments of the DGLAP and ERBL regions. The polynomiality Eq. (2.11) property manifestly hold for this set of Mellin moments.

Stated differently, this h
Toy P is such that: However h Toy P has a singularity at β = 1, which makes it non-summable over since: The related GPD is unchanged in both DGLAP and ERBL regions.