Pion Generalized Parton Distributions within a fully covariant constituent quark model

We extend the investigation of the Generalized Parton Distribution for a charged pion within a fully covariant constituent quark model, in two respects: (i) calculating the tensor distribution and (ii) adding the treatment of the evolution, needed for achieving a meaningful comparison with both the experimental parton distribution and the lattice evaluation of the so-called generalized form factors. Distinct features of our phenomenological covariant quark model are: (i) a 4D Ansatz for the pion Bethe-Salpeter amplitude, to be used in the Mandelstam formula for matrix elements of the relevant current operators, and (ii) only two parameters, namely a quark mass assumed to hold $m_q=~220$ MeV and a free parameter fixed through the value of the pion decay constant. The possibility of increasing the dynamical content of our covariant constituent quark model is briefly discussed in the context of the Nakanishi integral representation of the Bethe-Salpeter amplitude.


Introduction
The present theory of strong interaction, the Quantum ChromoDynamics (QCD), should in principle allow one to achieve a complete 3D description of hadrons, in terms of the Bjorken variable x B and the transverse momenta of the constituents. As it is well known, the needed non-perturbative description still represents a challenge, which motivates a large amount of valuable efforts, both on the experimental side (gathering new accurate data, that in turn impose stringent constraints on theoretical investigations) and the theoretical one (performing more and more refined lattice calcua e-mail: salmeg@roma1.infn.it lations and elaborating more and more reliable phenomenological models).
Heuristically, while the short-distance behavior of the hadronic state has been well understood, given the possibility of applying a perturbative approach, entailed by the asymptotic freedom, the long-range part of the hadronic state, which is governed by the confinement, requests non-perturbative tools, suitable for a highly non-linear dynamics. Coping with the difficult task to gain information on the hadronic state, in the whole range of its extension, has been the main motivation for elaborating phenomenological models, which in general play a helpful role in shedding light onto the nonperturbative regime.
Among the phenomenological approaches, covariant constituent quark models (CCQMs) represent an important step forward, since they exploit a quark-hadron vertex fulfilling the fundamental property of covariance with respect to the Poincaré group. Moreover, CCQM's based on the Light-front (LF) framework, introduced by Dirac in 1949 [1], with variables defined by a ± = a 0 ± a 3 and a ⊥ ≡ {a x , a y }, appear to be quite suitable for describing relativistic, interacting systems, like hadrons. Indeed, the LF framework has several appealing features (see, e.g., [2]), quite useful for exploring nowadays issues in hadronic phenomenology. Beyond the well-known fact that the dynamics onto the light-cone is naturally described in terms of LF variables, one should mention: (1) the straightforward separation of the global motion from the intrinsic one (related to the subgroup property of the LF boosts), (2) the largest number of kinematical (i.e. not affected by the interaction) Poincaré generators, (3) the large extent of triviality of the vacuum, within a LF field theory [2] (with the caveat of the zero-mode contributions). In particu-lar, for the pion, one can construct the following meaningful Fock expansion onto the null-plane: |π = |qq + |qq qq + |qq g . . .
where |qq is the valence component. It has to be recalled that an appealing feature of our approach, based on a covariant description of the quark-pion vertex (see [3,4] and references quoted therein), is the possibility of naturally taking into account contributions beyond the valence term.
The experimental efforts are very intense for singling out quantities that are sensitive to the dynamical features of the hadronic states. In particular, in the last decade, it has been recognized that a wealth of information on the 3D partonic structure of hadrons is contained in the Generalized Parton Distributions (GPDs) (see, e.g., Ref. [5] for a general presentation), as well as in the Transverse-Momentum Distributions (TMDs) (see, e.g., Ref. [6] for a detailed discussion). GPDs can be experimentally investigated through the deeply virtual Compton scattering, while TMDs can be studied through semi-inclusive deep inelastic scattering processes, which notably involve polarization degrees of freedom.
Our aim, is to provide a phenomenological model, that has the following main ingredients: (1) a 4D Ansatz for the pion Bethe-Salpeter amplitude, and (2) the generalization of the Mandelstam formula [7] for matrix elements of the relevant current operators (notice that the pion Bethe-Salpeter amplitude is needed in this formula). Remarkably, we introduce only two parameters, namely a constituent quark mass and a free parameter fixed through the value of the pion decay constant. Through our model, we investigate the pion state by thoroughly comparing the results with both experimental and lattice data relevant for the 3D description of the pion. In this paper, we complete the evaluation of the leading-twist pion GPDs, calculating the so-called tensor GPD (see Ref. [3] for the vector GPD and Ref. [8,9] for preliminary calculations of the tensor one). Moreover, in order to accomplish the previously mentioned comparisons, we consider the evolution of quantities that can be extracted from the GPDs, like the parton distribution function (PDF) and the generalized form factors (GFFs). We anticipate that only the leadingorder (LO) evolution has been implemented by using the standard code of Ref. [10]. In particular, the comparison has been performed between our LO results and the experimental pion PDF extracted in Ref. [11] (see Ref. [12] for the NLO extraction) and the available lattice calculations of GFFs as given in Refs. [13][14][15][16].
The paper is organized as follows. In Sect. 2, the general formalism and the definitions are briefly recalled. In Sect. 3, our Covariant Constituent Quark Model is presented. In Sect. 4, the LO evolution of the quantities we want to compare is thoroughly discussed, with a particular care to the determination of the initial scale of our model. In Sect. 5, the comparison of our results with both the experimental PDF and the available lattice calculations is presented. Finally in Sect. 6, the Conclusion are drawn.

Generalities
In this section, the physical quantities, GPDs and TMDs, that allow us to achieve a detailed 3D description of a pion are shortly introduced, since they represent the target of the investigation within our CCQM (for accurate and extensive reviews on GPDs, see, e.g., [5] and on TMDs, see, e.g., [6]). For the pion, given its null total angular momentum, one has two GPDs and two TMDs, at the leading twist.

Generalized Parton distributions
As is well known, GPDs are LF-boost-invariant functions, and they allow one to parametrize matrix elements (between hadronic states) involving quark and gluon fields. In particular, GPDs are off-diagonal (with respect to the hadron fourmomenta, i.e. p f = p i ) matrix elements of the quark-quark (or gluon-gluon) correlator projected onto the Dirac basis (see, e.g., Ref. [17] for a thorough investigation of the pion case). The appealing feature of GPDs is given by the ability of summarizing in a natural way information contained in several observables investigated in different kinematical regimes, like electromagnetic (em) form factors (FFs) or PDFs.
The pion has two leading-twist quark GPDs: (1) the vector, or no spin-flip, GPD, H I π (x, ξ, t), and (2) the tensor, or spinflip, GPD, E I π,T (x, ξ, t) (where I = IS,IV labels isoscalar and isovector GPDs, respectively). In order to avoid Wilsonline contributions, one can choose the light-cone gauge [5] and get and wherez ≡ {z + = z 0 + z 3 , z ⊥ }, ψ q (z) is the quark-field isodoublet and the standard GPD variables are given by with the initial LF momentum of the active quark equal to {k + − + /2, k ⊥ − ⊥ /2}. The factor of two multiplying the vector GPD is chosen for normalization purpose, so that for a charged pion one has where H u π = H IS π + H IV π , and H IS π is odd in x while H IV π is even (see, e.g. [3]). Finally, it is useful for the following to recall the relation with the parton distributions, q(x), viz.
At the present stage, only a few moments of the pion GPDs have been evaluated within lattice QCD, but they represent a valuable test ground for any phenomenological model that aspires to yield meaningful insights into the pion dynamics. In view of the numerical results discussed below, we briefly recall how the Mellin moments can be covariantly parametrized through the GFFs, which are the quantities adopted for comparing lattice calculations and phenomenological results.
The relation between the non-spin flip GPD and the em FF given in Eq. (4) for a charged pion can be in some sense generalized, if one considers Mellin moments of both vector and tensor GPDs. Then one obtains the corresponding GFFs. For instance, one can write the following Mellin moments of both vector and tensor GPDs for the u-quark (see [5,18] for a review) where the symbol [...] indicates the integer part of the argument. In Eqs. (6) and (7), A u n+1,2i (t) is a vector GFF for a u-quark and B u n+1,2i (t) a tensor GFF, respectively. It is worth noting that one can introduce a different decomposition in terms of isoscalar and isovector components instead of a flavor decomposition. In particular, if n + 1 is even (odd) one has an isoscalar (isovector) GFF. A striking feature is shown by the rhs of Eqs. (6) and (7), the so-called polinomiality, i.e. the dependence upon finite powers of the variable ξ . This polynomiality property follows from completely general properties like covariance, parity, and time-reversal invariance; for this reason it can be a good test for any model.
By considering the first vector and tensor moments one gets the following important relations: where B u 1,0 (0) = 0 is the tensor charge for n = 0, also called tensor anomalous magnetic moment (see Ref. [18]). Notably, Eq. (7) leads to the following relation, involving the Mellin moments of the PDF and A u n+1,0 (0), viz.: A physical interpretation of GFFs (see, e.g., [19][20][21]) can be achieved by properly generalizing the standard interpretation of the non-relativistic em FFs to a relativistic framework. Non-relativistically, the em FFs are the 3D Fourier transforms of intrinsic (Galilean-invariant) em distributions in the coordinate space (e.g., for the pion, one has the charge distribution, while, for the nucleon, one has both charge and magnetic distributions). In the relativistic case, one should consider Fourier transforms of GPDs, which depend upon variables invariant under LF boosts. Indeed, only the transverse part of μ can be trivially conjugated to variables in the coordinate space, while for x and ξ (proportional to + ) this is not possible. Therefore, keeping the description invariant for proper boosts (i.e. LF boosts), one can introduce 2D Fourier transforms with respect to ⊥ . Such a Fourier transform allows one to investigate the spatial distributions of the quarks in the so-called impact-parameter space (IPS). In particular, from Eqs. (6) and (7), it straightforwardly follows that, for ξ = 0, only A u n+1,0 ( 2 ) and B u n+1,0 ( 2 ) survive. Due to the LF-invariance of ξ , one has an infinite set of frames (Drell-Yan frames) where ξ = 0. In these frames, where + = 0 and 2 = − 2 ⊥ , one can introduce the above mentioned 2D Fourier transforms in a boost-invariant way (recall that, for a given reaction, the final state or both final and initial states have to be boosted). One can writẽ where b ⊥ = |b ⊥ | is the impact parameter. In general, the Fourier transform of GFFs, for ξ = 0, yield quark densities in the IPS [19][20][21]. In particular,Ã n (b ⊥ ) represents the probability density of finding an unpolarized quark in the pion at a certain distance b ⊥ from the transverse center of momentum.
In addition, if one considers the polarization degrees of freedom, then one introduces the probability density of finding a quark with a given transverse polarization, s ⊥ in a certain Drell-Yan frame. In the IPS, such a probability distribution is It is worth noting that the quark longitudinal (or helicity) distribution density is given only by the first term in Eq. (12), since the pion is a pseudoscalar meson and the term γ 5 / s L in the quark density operator has a vanishing expectation value, due to the parity invariance [16,22]. Equation (12) is quite rich of information and clearly indicates the pivotal role of GPDs for accessing the quark distribution in the IPS. Moreover, as a closing remark, one could exploit the spin-flip GPD E q π T to extract more elusive information on the quasi-particle nature of the constituent quarks, like their possible anomalous magnetic moments, once the vector current that governs the quark-photon coupling is suitably improved (see Sect. 5.4 and Ref. [23] for a discussion within the lattice framework).

Transverse-momentum distributions
TMDs are diagonal (in the pion four-momentum 1 ) matrix elements of the quark-quark (or gluon-gluon) correlator with the proper Wilson-line contributions (see, e.g., Ref. [17]) and suitable Dirac structures. Moreover, TMDs depend upon x and the quark transverse momentum, k ⊥ , that is, not the conjugate of b ⊥ . It should be pointed out that in general the Wilson-line effects must be carefully analyzed, due to the explicit dependence upon k ⊥ (recall that for GPDs such dependence is integrated out). At the leading twist, one has two TMDs, for the pion: the T-even f q 1 (x, |k ⊥ | 2 ), that yields the probability distribution to find an unpolarized quark with LF momentum {x, k ⊥ } in the pion, and the T-odd h q⊥ 1 (x, |k ⊥ | 2 , η), related to a transversely polarized quark and called a Boer-Mulders distribution [24].
The two TMDs allow one to parametrize the distribution of a quark with given LF momentum and transverse polarization, i.e. (see, e.g., Ref. [15,17]) where the dependence upon the variable η in h ⊥ 1 is generated by the Wilson-line effects, whose role is essential for investigating a non-vanishing h ⊥ 1 (see e.g. [24]). At the lowest order, the unpolarized TMD f q 1 , is given by the proper combination of the isoscalar and isovector components, that are defined by After integrating over k ⊥ , one gets the standard unpolarized parton distribution q(x), viz The T-odd TMD, h ⊥ 1 (x, |k ⊥ | 2 , η) needs a more careful analysis, since it vanishes at the lowest order in perturbation theory. As a matter of fact, it becomes proportional to the matrix elements which are equal to zero, due to the time-reversal invariance.
In order to get a non-vanishing Boer-Mulders distribution, one has to evaluate at least a first-order correction, involving Wilson lines (see, e.g., Refs. [17] and [25]). Moreover, by adopting the light-cone gauge and the advanced boundary condition for the gauge field, the effect of the Wilson lines (final state interaction effects) can be shifted into complex phases affecting the initial state (see, e.g., Ref. [26]).

The covariant constituent quark model
The main ingredients of our covariant constituent quark model are two: (1) the extension to the GPDs and TMDs of the Mandelstam formalism [7], originally introduced for calculating matrix elements of the em current operator when a relativistic interacting system is investigated, and (2) a model of the 4D quark-hadron vertex, or equivalently the Bethe-Salpeter amplitude, necessary for applying the Mandelstam approach. In particular, we have assumed a pion Bethe-Salpeter amplitude (BSA) with the following form: where p = p q + pq is the total momentum, t = ( p q − pq )/2 the relative momentum of the qq pair (by using the four-momenta k, and P previously introduced, one has t + p/2 = k − /2, and t − p/2 = k − P). In Eq. (18), S( p q ) = 1/( / p q − m q + ı ) is the fermion propagator and (t, p) the quark-pion vertex. In the present work, only the dominant Dirac structure has been assumed, viz.
(t, p) = γ 5 π (t, p) (19) with (t, p) a suitable momentum-dependent scalar function that contains the dynamical information (see the following subsections for more details). Indeed, Dirac structures contributing to (t, p) beyond γ 5 should be taken into account, but they have a minor impact on the pion BSA, as thoroughly discussed in Ref. [27].
For the sake of completeness, let us recall that the quarkpion vertex fulfills the homogeneous BS equation that reads as follows: where K(t, t ) is the kernel given by the infinite sum of irreducible diagrams (see, e.g., [28]). Finally, it is important to emphasize that our investigation, based on a covariant description of the quark-pion vertex, naturally goes beyond a purely valence description of the pion [3,4].

The Mandelstam formula for the electromagnetic current
The Mandelstam formula allows one to express the matrix elements of the em current of a composite bound system, within a field theoretical approach [7]. It has been applied for evaluating the FFs of both pion [29][30][31] and nucleon [32], obtaining a nice description of both space-and timelike FFs. Furthermore, it has been exploited for calculating the vector GPD of the pion [3,4] and for a preliminary evaluation of the tensor GPD [8,9]. For instance, in the case of the em spacelike FF of the pion, the Mandelstam formula, where the quark-pion vertex given in Eq. (19) is adopted, reads (see, e.g., Ref. [29][30][31]) where R = 2N c m 2 q / f 2 π , f π is the pion decay constant N c = 3 the number of colors, m q the CQ mass and V μ (k, q) the quark-photon vertex, which we have simplified to γ μ in the spacelike region. In presence of a CQ, one could add to the bare vector current a term proportional to an anomalous magnetic moment, namely a term like as in Ref. [23] (where an improved vector current within a lattice framework has been adopted). It should be pointed out that the above mentioned anomalous magnetic moment is not used in the present work. Within CCQM, the expression of the decay constant in terms of π reads (cf Ref. [29]) is the valence wave function. It should be recalled that , properly integrated on κ ⊥ , yields the pion distribution amplitude (DA) (see Eq. (60) and, e.g., Ref. [5] for a general discussion of the DAs and their evolution). The generalization of Eq. (21) to the case of GPDs can be found in Ref. [3,4] for the vector GPD, and in [8,9] for the tensor one, but for the sake of completeness, let us give the expression of both vector and tensor GPDs for the u quark, viz and where j = 1, 2. The δ function allows one to have the correct support for the active quark, i.e. when |ξ | ≤ x ≤ 1. This kinematical region corresponds to the so-called Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) region [33][34][35], or valence region. Moreover, CCQM is able to address also kinematical region beyond the valence one, i.e. −ξ ≤ x ≤ ξ , given the covariance property. This region is called the Efremov-Radyushkin-Brodsky-Lepage (ERBL) region [36,37], or non-valence region. If one adopts a Breit frame with + = − − ≥ 0, then the ERBL region can be investigated.
As a matter of fact, in such a frame one can access the whole range of the variable ξ , i.e. −1 ≤ ξ ≤ 1, and analyze both valence and non-valence regions within the same approach. This allows one to shed light on the interesting topic of the smooth transition from the DGLAP (valence) regime to the ERBL (non-valence) one. The expression of the unpolarized TMD f can easily be obtained from the integrand of the vector GPD, Eq. (24), by recalling the relation in Eq. (16).

The four-momentum dependence of the Bethe-Salpeter amplitude
As above mentioned, in our CCQM we focus on the main contribution to the pion BSA, i.e. the term containing the Dirac matrix γ 5 . This implies that we have to consider only one scalar function for describing the dependence upon the four-momenta present in the problem. Unfortunately, solutions of the homogeneous BSE for hadrons are still lacking in Minkowski space given the extraordinary complexity of QCD, nonetheless very relevant investigations have been carried out in Euclidean space, within the lattice framework [13][14][15][16] or combining BSE and Dyson-Schwinger equation (DSE) (see, e.g., [38] and references quoted therein) or by exploiting a 3D reduction of the BSE itself (see, e.g., [39]). On the other hand, since we would carry on a comparison with a wide set of data, from both experiments and lattice, we resort to adopt a phenomenological Ansatz, which depends remarkably upon only two parameters. This allows us to explore the potentiality of the Mandelstam approach in capturing the main features of the physical quantities under consideration, while having a reasonable predictive power, given the small set of free parameters. The following analytic covariant Ansatz for the momentum dependence of the BSA has been adopted (t, p) where the parameter m R is adjusted to fit f π , while the constants C are fixed through the charge normalization, F π (t = 0) = 1, which amounts to the standard normalization of the BSA, but in impulse approximation. It is worth noting that the expression in Eq. (26) can be cast (see below) in a form suggested by the integral representation of the 4D n-leg transition amplitudes (we are actually interested to the 3-leg amplitude, i.e. the vertex π → qq) elaborated by Nakanishi in the 1960s [40], within a perturbation-theory framework. To quickly illustrate the appealing features of this integral representation, one should consider the n-leg transition amplitude for a many-scalar interacting system, and the infinite set of Feynman diagrams contributing to determine the amplitude itself. In this case, it turns out that the amplitude is given by the folding of a weight function (called the Nakanishi weight function) and a denominator (with some exponent) that contains all the independent scalar products obtained from the n external four-momenta. It has to be pointed out that the analytic behavior of the amplitude is fully determined by such a denominator, and this clearly makes the Nakanishi integral representation a valuable tool for investigating 4D transition amplitudes. For n = 3, one can apply the integral representation to the vertex function for a system composed of two constituents, and we explicitly discuss the analytic structure, i.e. the core of the physical content. Another pivotal motivation that increases the interest on the Nakanishi framework is given by the following computational finding: even if the Nakanishi integral representation has been formally established by considering the whole infinite set of the Feynman diagrams contributing to an amplitude, i.e. a perturbative regime, it has been numerically shown that also in a non-perturbative framework, like the homogeneous BSE (relevant for describing bound systems), the Nakanishi representation plays an essential role for obtaining actual solutions for the vertex function or, equivalently, for the BSA. Applying the Nakanishi representation as an Ansatz for the solution of the BSE one can determine the unknown Nakanishi weight function and achieve a genuine numerical solution of the BSE in Minkowski space. This approach has been applied to the ladder BSE for two-scalar and two-fermion systems (see, e.g., Refs. [41][42][43][44][45] for the Nakanishi approach in Minkowski space and, for the sake of comparison, Refs. [46,47] for two-fermion systems within the Euclidean hyperspherical approach), opening a viable path for phenomenological studies within a non-perturbative regime.
Within the Nakanishi approach, the vertex function (or three-leg amplitude) can be written as follows: where , one obtains Eq. (26). It should be pointed out that while waiting for numerical solutions of the two-fermion system with more refined phenomenological kernels (for the ladder approximation see Ref. [42]), one could perform an intermediate step, still in the realm of Ansatzes, substituting Eq. (26) with Eq. (27), but adopting a different choice of the Nakanishi weight function, e.g. by substituting the simple delta-like form with more realistic functions (see e.g. [44] for the Nakanishi weight functions of a two-scalar system obtained by actually solving the homogeneous BSE in the ladder approximation). In order to set the reference line for the next steps in the elaboration of our CCQM (presented elsewhere), we will adopt the very manageable form given in Eq. (26), in the following comparisons with the experimental and lattice results (see below, Sect. 5).

Evolution of Mellin moments and GFF's
In order to compare our results for PDF and GFFs with experimental data and lattice calculations, it is fundamental to suitably evolve the CCQM outcomes, from the unknown scale μ CCQ to the needed ones, namely μ exp and μ LAT .
Our strategy for determining an acceptable μ CCQ is to study the evolution of the non-singlet PDF Mellin within a LO framework, considering flavor numbers up to N f = 4. It should be pointed out that the choice to adopt the LO framework seems to be well motivated by the phenomenological nature of the CCQM, and by the present uncertainties still affecting both experimental and lattice GFFs.
Let us shortly summarize our procedure for assigning a scale μ to our calculations. The main ingredient to be considered are the Mellin moments of the non-singlet distribution where f NS is related to the unpolarized GPD, as follows: Mellin moments evolve from a scale μ 0 to the scale μ through very simple expressions (see, e.g., [48]), which for the nonsinglet, singlet, and gluon moments read where In Eqs. (31) and (30), the LO anomalous dimensions are indicated by γ (0) ab while the 2 × 2 matrix (0) (n) is given by Let us recall that each anomalous dimension γ (0) ab (n) is obtained from the corresponding LO splitting function. In particular, for the unpolarized case, one has (see, e.g., [48]) where By taking into account the eigenvalues of (0) (n), given by (see [49] for details) one can write the 2 × 2 matrix in terms of projectors and eigenstates as follows: with They fulfill the usual projector properties, i.e. P + + P − = 1, Solutions of Eqs. (31) and (30) are given by Notably, Eq. (44) can be put in a more simple form by using the eigenvalues, γ ± , and the corresponding projectors P ± , viz. [49] − Indeed, we are interested to actually evolve only the moment n = 1 of f NS , since for this moment we can find several lattice calculations (but using different approximations; cf Sect. It should be pointed out that, as explicitly shown in Eqs. (38) and (46), α s (μ, N f ) depends upon the number of flavors N f , at a given scale. Indeed, one has to be particularly careful about the energy scales involved, when one moves from a relatively low μ i = 1 GeV to μ LAT = 2 GeV, large enough to produce a new quark flavor, so that N f increases from 3 to 4. In practice, a two-step procedure has been adopted for moving from α s (μ i , N f = 3) to α s (μ LAT , N f = 4), by properly changing β 0 (N f ), at the threshold μ = m c , i.e. the mass of the charm. In the following, it is also useful to define, at a given energy scale and number of flavors,

QCD evolution of GFFs
Similarly to the more familiar case of PDFs, where the QCD interaction among partons lead to ultraviolet divergences which are factored out and absorbed into a dependence upon the energy scale, also in the case of GPDs one has to deal with the issue of finding and solving evolution equations. As a matter of fact, GPDs do not depend on three variables but on four, namely H (x, ξ, t, μ) and E(x, ξ, t, μ). However, the evolution kernel does not depend on t, so that the relevant variables for the evolution are x, ξ , and μ. One should keep in mind that the evolution of GPDs is produced by the combination of two regimes: (1) the one pertaining to the valence region (|x| > |ξ |) and (2) the one pertaining to the non-valence region (|x| < |ξ |). One could roughly say that the evolution of GPDs interpolates [5] between the two regions and therefore the evolution kernel has to take into account the suitable physical content. In particular, in the valence region a kernel acts with a structure like the one present in the DGLAP equations, while in the non valence region a modified ERBL kernel is involved (see Refs. [51] and [52] for details on the evolution of vector and tensor GPDs, respectively). In our actual comparison, we do not consider the full GPDs, but rather their Mellin moments, since they can be in principle addressed by the lattice calculations. As a matter of fact, GFFs covariantly parametrize the Mellin moments of GPD (see Eqs. (7) and (6)), and evolve through a suitable generalization of Eqs. (43) and (44) (see Refs. [49,[53][54][55][56][57][58]). Let us recall, however, that GFFs are the coefficients of polynomials in ξ that yield the Mellin moments of GPDs and not the Mellin moments themselves: for this reason in general the equations describing GFFs evolution are more complicated than Eqs. (43) and (44). Indeed one can find some notable exceptions where the equations have a simple multiplicative structure.
For the vector GFFs A I ni (t, μ 2 ), one should recall that the evolution of the isoscalar (singlet) GPD, and consequently the evolution of the corresponding Mellin moments, is coupled with the evolution of the gluonic component. This leads one to separate the evolution of GFFs with even and odd n, since for symmetry reasons the even GFFs come from the isoscalar GPDs, while the odd ones come from the isovector GPDs. By repeating the main steps given in Ref. [54] (see also [55,56] where general discussions are presented) for obtaining the evolution equation of both non singlet and singlet vector GFFs, we can express the results in Ref. [54] also as follows: with 0 ≤ ≤ k and . (49) For the singlet vector GFFs we get with 0 ≤ ≤ k + 1, A 00 = 0 and where (0) V is the same 2 × 2 matrix as defined in (33), which depends upon N f through γ (0) GG . In Eqs. (48) and (50), (k) is the usual Euler function. Notice that the dependence upon t in the GFF is not involved in the evolution. Some example of explicit evolution equations are where A 22 (t, μ) is a 2D vector, with quark and gluon components, and L 2 a 2×2 matrix. For A n0 (t, μ) the evolution equations become exactly equal to Eqs. (43) and (44) for the odd and even n's, respectively. Moreover, since γ (0) qq (n = 0) = 0 then L 1 = 1. This is expected since A 10 (t = 0) is the charge (A 10 (t) is the em form factor), namely a measurable quantity and therefore it cannot evolve.
For the tensor GFF B I ni (t, μ), analogous arguments can be carried out, but with a great simplification. In fact, at LO the gluon-quark and quark-gluon transition amplitudes that lead to the corresponding splitting functions are vanishing for the helicity conservation (recall that E T (x, ξ, t) is related to an expectation value with a transversely polarized quark and describes helicity-flip transitions), therefore the anomalous dimension matrix where with a transverse anomalous dimension (notice a factor of 2 difference with [49], due to the different normalization) given by For the sake of completeness, the gluon transverse LO anomalous dimension reads

Numerical results
The reliability of the quark-pion vertex (26), introduced for obtaining the Bethe-Salpeter amplitude (18), has been first checked by comparing our results for a charged pion with the most accurate experimental data not affected by the evolution, i.e. the spacelike em form factor. Theoretically, the em form factor, F π (t), is given by the GFF A 10 (t). In Fig.  1 the results obtained by different models for the em FF are shown as a function of (−t), together with the experimental data. To avoid the use of a log plot, which prevents a detailed analysis, the FF has been divided by the monopole function GeV. Experimental data: as quoted in [3] LF CQM where m q = 0.265 GeV and a dressed quarkphoton vertex were adopted [30,31]; (2) a fit to the lattice data obtained in [23]. It has to be pointed out that the fit to the lattice data was presented in [23] itself, and it has the following monopole expression: with M(m phys π ) = 0.727 GeV. The lattice results were actually obtained for a pion mass m π = 0.600 GeV and then extrapolated to the physical pion mass m phys π = 0.140 GeV, up to t = −4 GeV 2 (see [23]). In Fig. 1, for the sake of presentation, the monopole fit (58) has been arbitrarily extended up to t = −10 GeV 2 . Once the CQ mass is assigned, our CCQM model depends upon only one free parameter, the regulator mass m R in Eq. (26). The value of m R is fixed by calculating the pion decay constant f π (cf Eq. (22)), while the constant C is determined through F π (0), as already mentioned in Sect. 3.2. The PDG experimental value f exp π = 0.0922 GeV [60] has been adopted. In particular, for m q = 0.200, 210, 220 GeV, we got m R = 1.453, 1.320, 1.192 GeV, respectively.
The following comments are in order: (1) a nice agreement between the CCQM results and the experimental FF at low momentum transfer leads to reproduce the experimental value of the charge radius, r 2 exp = 0.67 ± 0.02 fm; (2) beyond −t = 0.5 GeV 2 , CCQM results begin to reveal an interesting sensitivity upon the CQ mass, since if one changes the CQ mass by a 5 % then the corresponding CCQM FF changes by 10-15 %, at high momentum transfer; (3) the CCQM FFs seem to have a similar curvature of the data at high momentum transfer, but in order to draw reliable conclusions, useful for extracting information (and then improving the CCQM), it is necessary to have more accurate data, for −t ≥ 1 GeV 2 .
Another interesting data set to be compared with is given by the photon-pion transition form factor, F * π (−t), measured in the process γ γ * → π 0 [61,62]. Within CCQM, only the LO asymptotic value of such transition FF can be evaluated without adding new ingredients (see for a wide discussion, e.g., Ref. [63] and references quoted therein). As a matter of fact, one gets for high (−t), at LO in pQCD, where φ π (ξ, |t|) is the pion DA evaluated at the scale |t|. The CCQM result (with an undetermined scale for the moment, see the next subsection) is given by where (ξ m π , κ ⊥ ) is defined in Eq. (23). The normalization of φ π follows from Eq. (22). It should be anticipated that CCQM results, both non-evolved and evolved, as shown in Fig. 3, resemble the asymptotic pion DA obtained within the pQCD framework, i.e. φ asy π (ξ ) = 6ξ(1 − ξ), which in turn yields (−t) F * π (−t) → 2 f π ; see Refs. [64,65]. In the following, the values m q = 0.220 geV and m R = 1.192 GeV will be adopted.

Looking for the CCQM energy scale
As is well known, the em FF is not affected by the issue of the evolution, while the other quantities we are interested in, namely the PDF and the GFFs [as well as the DA; see Eq. (60)], have to be properly evolved.
A necessary step for going forward is to assign a resolution scale to CCQM. In order to perform this step, we have taken lattice estimates of the first Mellin moment of f NS (x, μ), whose evolution is determined only by the quark contribution, as normalization of our CCQM (roughly speaking). The starting point is the calculation of both the unpolarized GPD, f NS (x) = 2H I =1 (x, 0, 0), and the corresponding Mellin moments, within our CCQM. In particular, these quantities are shown in Table 1 up to n = 3. To emphasize that there is no direct way to gather information as regards the energy scale μ 0 , a question mark is put in the table.
In the literature there are various lattice results for the first moment at the energy scale of μ = 2 GeV, and we have exploited the ones shown in Table 2. It is worth noticing that the lattice results are not too far from a phenomenological estimate, x phe (μ = 2 GeV), that one can deduce by applying a LO backward-evolution to the value given in Ref. [12], x phe (μ = 5.2 GeV) = 0.217 (11), obtained after a NLO Table 1 Mellin moments of f NS (x) up to n = 3, evaluated within the CCQM with the quark-pion vertex given in Eq. (26), m q = 0.220 GeV, and m R = 1.192 GeV. The energy scale, μ 0 has to be determined (see text) Table 2 Recent lattice results for the first Mellin moment of the nonsinglet f NS (x), at the energy scale μ LAT = 2 GeV. The first and the second lines are the results obtained from unquenched lattice QCD calculations [14,66], while the third result has been obtained in quenched lattice QCD [67] Ref.
After establishing the set of lattice data, we need the value of α LO s at μ = 2 GeV where N f = 4. This value has been obtained starting from α LO s (μ = 1 GeV) = 0.68183 obtained in Ref. [50]. Notice that at the scale μ = 1 GeV only three flavors are active. Then, by using m c = 1. (1, ?). In detail, we calculate first (cf. Eq. (43)) the lattice result at the charm mass scale, viz  Table 2 where M NS (1, μ 0 ) corresponds to our CCQM calculation and β 0 (3) = 9. After determining α LO s (μ 0 , 3), μ 0 is easily found through The results for μ CCQ obtained from the above procedure, applied to the three lattice data, are shown in Table 3 for m q = 0.220 GeV and m R = 1.192 GeV. In particular, the values in the third column of Table 3 are used in the next sections as starting values for the evolution of both the non-singlet PDF and the GFFs. The difference between the three values Table 3 Energy scale of CCQM, μ 0 , as determined from (1) the first Mellin moments calculated within a lattice framework in Refs. [14,66,67] and (2)  Ref.
Lat. 07 [14] 0.271 0.549 1.64 South [66] 0.249 0.506 2.04 χLF [67] 0.243 0.496 2.17 Table 4 Comparison for the second and third Mellin moments of the non-singlet f NS (x), at the energy scale μ LAT = 2 GeV, between the unquenched lattice results of Ref. [14] and the evolved CCQM, where the theoretical uncertainty is generated by the three values for the CCQM initial scale shown in Table 3 x 2 x 3 Lat. 07 [14] 0.128 (18) 0.074 (27) CCQM 0.105 (11) 0.055 (7) of μ 0 in the Table 3 is assumed as a theoretical uncertainty of our results. To complete this subsection, in Table 4, the comparison with the lattice calculation of Ref. [14] for the second and the third Mellin moments is presented.

The evolution of the non-singlet PDF and the comparison with the experimental data
The non-singlet PDF, as already explained, is the simplest to be evolved since one does not need information on the gluon distribution. The evolution has been performed using the FORTRAN code described in [10] that adopts a bruteforce method to solve the LO DGLAP equation for the distribution x f NS (x), and it requests as input the values of (1) μ, the final scale, and (2) the initial N f QCD and μ 0 , as given in Table 3. An important detail in our calculations should be pointed out. For all the values of μ 0 , the evolution has been performed in two steps: first x f CCQM NS (x) has been evolved from μ 0 up to m c = 1.4 GeV and then from m c up to μ = 4 GeV, the energy scale of the experimental data [11]. This is necessary for taking into account the variation of N f , QCD (recall that N f =4 QCD (μ = 2 GeV ) is 0.322 GeV) and consequently α LO s (μ). In Fig. 2, the dashed line is the non-evolved CCQM calculation with m q = 0.220 GeV and m R = 1.192 GeV, while the solid and the dot-dashed lines correspond to our evolved CCQM starting from the initial scales μ 0 = 0.549 GeV and μ 0 = 0.496 GeV, respectively. The differences between the evolved calculations can be interpreted as the theoretical uncertainty of our calculations. However, it is very interesting that our LO-evolved calculations nicely agree with the  Table 3. Full dots experimental data at the energy scale μ = 4 GeV, as given in Ref. [11] experimental data of Ref. [11] for x > 0.5 (see also the same agreement achieved within the chiral quark model of Ref. [53]). On the other hand, it has to be pointed out that refined calculations, like (1) the ones of Refs. [72,73] based on the Euclidean Dyson-Schwinger equation for the selfenergy and (2) the NLO calculation of Ref. [74] based on a soft-gluon resummation, underestimate the PDF tail of the experimental data from Ref. [11], while agree with the analysis of the same experimental data carried out in Ref. [12], within a NLO framework. The reanalysis of the experimental data leads to a tail for large x that has a rather different derivative with respect to the original data from Ref. [11]. For the sake of completeness, in Fig. 3, the CCQM pion DA is presented together with the results at the energy scale μ = 1 GeV and μ = 6 GeV. It is worth noticing that our CCQM evolves toward the pQCD asymptotic pion DA φ π (ξ ) = 6 ξ (1 − ξ) (see, e.g., [64,65]) as the energy scale increases. Analogous results are obtained within the chiral quark model of Ref. [75].

The tensor GPD
We have extended to the tensor GPD our CCQM model already applied to the vector GPD in Refs. [3,4], and in Fig.  4, our final results are shown for some values of the variable ξ and t, but for 0 ≤ x ≤ 1 (preliminary results were presented in Refs. [8,9]). The GPD for negative values of x can be obtained by exploiting the fact that E IS π T (x, ξ, t) is antisymmetric if x → − x, while E IV π T (x, ξ, t) is symmetric (see, e.g., Ref. [5] for details). It has to be pointed out that for ξ → 0 the valence component is dominant (DGLAP regime) while for ξ → 1 the non-valence term is acting (ERBL regime). In view of that, a peak is expected around x ∼ 1 for ξ → 1, as discussed in Refs. [3,4] for the vector GPD.
It is worth mentioning that both isoscalar and isovector tensor GPD calculated within the chiral quark model of Ref. [58] qualitatively show the same pattern (see also Ref. [8] and references therein, quoted for a comparison with results obtained within the LF Hamiltonian dynamics framework).

The evolution of the GFFs and the comparison with lattice data
The first vector GFF A q 10 , i.e. the em FF, is experimentally known, while the other GFFs can be investigated only from the theoretical side. In particular, A q 20 (t, μ), A q 22 (t, μ), B q 10 (t, μ) and B q 20 (t, μ) have been calculated within the lattice framework at the scale μ = 2 GeV [13,15,16]. In this subsection, the comparison between our CCQM predictions and the above mentioned lattice evaluations is presented. It is important to notice that other model calculations of both vector and tensor GFFs are available in the literature (see, e.g., [57][58][59][75][76][77]).
To proceed, we have calculated both vector and tensor GPDs, and then we have extracted the relevant GFFs, by exploiting the polynomiality shown in Eqs. (7) and (6) (see also [3]). The main issue to be addressed in order to perform the mentioned comparison with the lattice data is the evolution of our calculations up to μ LAT = 2 GeV. In the simpler case, represented by the tensor GFFs, the LO evolution of the quark contribution is uncoupled from the gluon one. In particular, the two-step procedure μ CCQ → m c → μ LAT has been Then for μ CCQ ≤ μ < m c one has N f = 3 and gets For m c ≤ μ ≤ μ LAT , the flavor number is N f = 4 and one has In the case of A 20 (t, μ) and A 22 (t, μ) the evolution equation is more complicated, since both GFFs evolve through the following expression: where, for both scales, one has From the definition (51), the exponent in L 2 is a 2 × 2 matrix (see also Eqs. (33), (35), (36), (37) and (34)) that for N f = 3 reads with eigenvalues (see Eq. (39)) At the valence scale, the gluon contribution is vanishing, and therefore one has x q = 1/2. Indeed, the CCQM result amounts to x (μ CCQ ) = 0.47, namely the momentum sum rule is not completely saturated by the valence component at the CCQM scale μ CCQ . This difference originates from the fact that we have a covariant description of the pion vertex, and therefore we have not only a contribution from the valence LF wave function (i.e. the amplitude of the Fock component with the lowest number of constituents), but also from components of the Fock expansion of the pion state beyond the constituent one, like |qq; qq . Without the gluon term at the initial scale (the assumed valence one), A q 2i (t, m c ) is given by (cf. Eqs. (40) and (41)) where and For N f = 4, Eq. (73) changes, since both β 0 and γ (0) GG (1) depend on the flavor number. Therefore with eigenvalues Then the evolution in the second step from m c → μ LAT = 2 GeV reads where .
It should be pointed out that the GFFs A q 2i evolve multiplicatively (recall that the evolution is not influenced by the value of t), given the absence of the gluon contribution at the valence scale, viz From Eq. (80), one realizes that the ratio is also called a charge) can be compared with the same ratio obtained at a different scale, e.g. at μ CCQ . It is understood that the same holds for the tensor GFF. In Figs. 5 and 6, the tensor GFFs B q 10 (t) and B q 20 (t), normalized to their own charges, are shown for both the CCQM model, with m q = 0.220 GeV and m R = 1.192 GeV, and the lattice framework [13,16]. In particular the lattice data are represented by a shaded area, generated by the envelope of curves that fit the lattice data with their uncertainties. In Refs. [13,16], the lattice data have been first extrapolated to the physical pion mass through a simple quadratic (in m π ) expression, and then fitted by the following pole form: where p j and M j are pairs of adjusted parameters, shown in Table 5, for the sake of completeness. In Figs. 7 and 8, the CCQM A 2,0 (t) and A 2,2 (t) (with CCQM parameters different from the ones adopted in Ref. [3]) are presented together with the corresponding lattice results.  Fig. 7 The vector GFF A q 2,0 (t), an isoscalar one, normalized to its own charge. Dashed line CCQM result, corresponding to m q = 0.220 GeV and m R = 1.192 GeV in the vertex (26). The shaded area indicates the lattice data [13] extrapolated to the pion physical mass m π = 0.140 GeV (see text)   If one is interested in a comparison that involves the full GFFs, then it is necessary to specify the scale and, accordingly, to evolve our CCQM results. In particular, since we have a multiplicative evolution, it is sufficient (1) to evolve only the value at t = 0, namely the ones collected in Table  6, through Eqs. (66), (66), (67), (68), (73) and (78) and then (2) to use Eq. (80). As in the case of the evolution of the PDF, we considered the three possible values of μ 0 listed in Table 3. The results are shown in Table 7, together with lattice data [13,16,78] and model calculations, obtained from a chiral quark model [57,76] and an instanton vacuum model [59]. It should be pointed out that within the chiral perturbation theory (see Ref. [79]) one should have the following relation between the so-called gravitational FFs: . This relation is verified by the lattice results, while CCQM does not. Moreover, one should observe that A q 20 (0) slightly differs from x at μ LAT (see Table 2), which contains both quark and gluon contributions.
To have a better understanding of the quality of the comparison between our CCQM results and the lattice data shown in Table 7, we have added our calculation, at t = 0, in Fig. 9, where the lattice results from [13], extrapolated at the physical pion mass, are presented for A q 2,0 (t, μ LAT ) and A q 2,2 (t, μ LAT ). In Fig. 9, the stars at t = 0 represent the CCQM values evolved at the lattice scale (the size of the symbols is roughly proportional to the uncertainties of the initial μ CCQ (cf Sect. 5.1), while the shaded area is the uncertainties produced by the fits to the lattice data, as elaborated in Ref. [13]. It is clear that in order to have a conclusive comparison a more wide lattice data set is necessary, but on the other hand it is impressive that a small quantity, like A q 2,2 (t, μ LAT ), can be extracted with a quite reasonable extent of reliability. In Figs. 10 and 11, analogous comparisons for B q 1,0 (t = 0, μ LAT ) and B q 2,0 (t = 0, μ LAT ) are shown. In particular, Fig. 10 contains both B 1,0 (t = 0, μ LAT ) and B q 2,0 (t = 0, μ LAT ), evaluated within CCQM (stars) and within the chiral quark model of Ref. [57] with different m π . In the figure the lattice data of Ref. [16] are also present. Again, the comparison between values at t = 0 and physical pion mass appears non trivial. In Fig. 11, a recent lattice calculation of B q 1,0 (t = 0, μ LAT ) [78] is compared with our CCQM (triangles). In general, one has an overall agreement, a little bit better for B q 2,0 (0). The knowledge of GFFs allows one to investigate the probability density ρ n (b ⊥ , s ⊥ ) for a transversely polarized u-quark (cf Eq. (12)). In particular, one can address the 3D structure of the pion in the impact parameter space. For instance, one can calculate the average transverse shifts when the quark is polarized along the x-axis, i.e. s ⊥ ≡ {1, 0}. The shift for a given n is given by [16] From the CCQM values evolved at μ LAT , shown in Table 7, one can construct the shifts for n = 1, 2, and then compare with the corresponding lattice results, as given in Ref. [16]. In Table 8, the comparison is shown (recall that A q 1,0 (t = 0) = 1). Obviously, the same observations relevant for Table  7 can be also repeated for Table 8, since it contains the same information but presented in a different context. The values shown in Table 8 indicate that even the simple version of a CCQM is able to reproduce a distortion of the Table 7 GFFs for t = 0 at a scale μ LAT = 2 GeV. The first three rows contain the evolved (see text) CCQM results for m q = 0.220 GeV and m R = 1.192 GeV. The fourth and fifth rows show the lattice extrapolations at the pion physical mass obtained in Refs. [13,16] and in Ref. [78], respectively. The sixth and seventh rows present the calculations from the chiral quark model of Refs. [57,76] and from the instanton vacuum model of Ref. [59], respectively. Notice that the results from [57,76] were not explicitly written in the works, so that they have been extrapolated by the plots presented there Refs. [13,16] 0  . 9 The lattice GFFs A q 2,0 (t, μ LAT ) and A q 2,2 (t, μ LAT ) of Ref. [13] and the CCQM A q 2,0 (t = 0, μ LAT ) and A q 2,2 (t = 0, μ LAT ). Stars CCQM result evolved at μ LAT = 2 GeV (the size of the symbols is roughly proportional to the uncertainties on CCQM initial scale μ CCQ ). Shaded area uncertainties of the fits to the lattice data, as estimated in Ref. [13] (see text) (Adapted from Ref. [13]) transverse density in a direction perpendicular to the quark polarization, and in turn they demonstrate the presence of a non-trivial correlation between the orbital angular momenta and the spin of the constituents inside a pseudoscalar hadron, which attracts a great interest from both experimental and theoretical side (see, e.g., [20]).

Conclusions
A simple, but fully covariant constituent quark model has been exploited for investigating the phenomenology of the leading-order generalized parton distributions of the pion.  1,0 (t = 0, μ LAT ) and B q 2,0 (t = 0, μ LAT ), divided by m π and the corresponding results from the chiral quark model of Ref. [57] and the lattice data of Ref. [16]. Stars CCQM results evolved at μ LAT = 2 for B q 1,0 (0) (upper one) and for B q 2,0 (0) (lower one) (the size of the symbols is roughly proportional to the uncertainties on our initial scale μ CCQ as illustrated in Sect. 5.1). Solid lines results from the chiral quark model of Ref. [57], vs m 2 π . Data points: lattice calculations from Ref. [16]. The vertical line corresponds to the physical pion mass (Adapted from Ref. [57]) The main ingredients of the approach are (1) the generalization of the Mandelstam formula, applied in the seminal work of Ref. [7] to matrix elements of the em current operator between states of a relativistic composite system, and (2) an Ansatz of the Bethe-Salpeter amplitude for describing the quark-pion vertex. Their combination produces a very effective tool that allows a careful phenomenological investigation of the pion, as shown in detail through the evaluation of both vector and tensor pion GPDs. We have also taken into account, at the leading order, the evolution for obtain-  Fig. 11 Comparison between the CCQM B q 1,0 (0) evolved at μ LAT = 2 GeV and the lattice calculation from Ref. [78]. Circles lattice data from [78] for different values of m π . Square extrapolation of the previous lattice data to the physical pion mass, m phys π , as it has been carried out in Ref. [78]. Triangle CCQM result evolved at μ LAT = 2 GeV (the size of the symbols is roughly proportional to the uncertainties on our initial scale μ CCQ , as illustrated in Sect. 5.1) (Adapted from Ref. [78]) Table 8 Mean shifts along the direction perpendicular to the u-quark transverse polarization, s ⊥ ≡ {1, 0}, for n = 1, 2 (cf Eq. (82)). The CCQM results have been constructed from the values in Table 7 (notice  that the uncertainties originate with the three values listed there) CCQMm π = 140 MeV Lattice [16] Lattice [78] b y 1 0.0901 ± 0.0015 fm 0.151 ± 0.024 fm 0.137 ± 0.007 fm b y 2 0.0796 ± 0.001 fm 0.106 ± 0.028 fm ing a meaningful comparison with both experimental data (see Fig. 2 for the comparison with the PDF extracted from the Drell-Yan data in Ref. [11]) and lattice calculations of generalized form factors. Summarizing, the CCQM proves to be quite satisfactory in describing the pion phenomenology, especially considering that the model involves relatively simple calculations and actually admits only one really free parameter (the mass m q of the constituent quark, since m R is constrained by f π ). It is worth noting that the CCQM is elaborated in Minkowski space, and the overall agreement we have shown with the lattice data, obtained in Euclidean space, could be an interesting source of information on the interplay of calculations performed in the two spaces, with a particular attention to the issue of the analytic behavior. In the future the present model could be substantially improved by enriching the analytic structure of the pion BS amplitude through a dynamical approach based on the solution of the BSE via the Nakanishi integral representation [40] (cf. Sect. 3.2), supplemented with a phenomenological kernel. In perspective, given the simplicity and the effectiveness of the approach, one could aim at applying the same model to more complex hadrons than the pion, e.g. to the nucleon within a quark-diquark framework.