Spatial and Momentum Imaging of the Pion and Kaon

Chao Shi, 2, ∗ Kyle Bednar, Ian C. Cloët, and Adam Freese 1Department of Nuclear Science and Technology, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China 2Collaborative Innovation Center of Radiation Medicine of Jiangsu Higher Education Institutions, Nanjing 211106, China 3Center for Nuclear Research, Department of Physics, Kent State University, Kent OH 44242 USA 4Physics Division, Argonne National Laboratory, Argonne, IL 60439 USA


I. INTRODUCTION
Multi-dimensional images of the partonic structure of hadrons are provided by the generalized parton distribution functions (GPDs) [1][2][3] and transverse momentum dependent parton distributions functions (TMDs) [4]. These images encode abundant structural information about hadrons, e.g., the GPDs provide a unified description of form factors and parton distribution functions (PDFs), where the former is related to a hadron's spatial extent and the latter describes the light-cone momentum distribution of partons [5,6] within a hadron. Through xweighted moments the GPDs are connected with hadron matrix elements of the energy-momentum tensor, and therefore shed-light on the spin, energy, and pressure distributions within hadrons [2,7]. Experimentally, GPDs are accessible through hard exclusive processes like deeply virtual Compton scattering (DVCS) or deeply virtual meson production (DVMP). TMDs illustrate the transverse motion of the partons in 3-dimensional momentum space, and therefore complement GPDs. Hadron TMDs can be extracted from semi-inclusive deep inelastic scattering (SIDIS) or Drell-Yan processes.
Calculating GPDs and TMDs directly from the fundamental theory, quantum chromodynamics (QCD), has proven very challenging. Lattice QCD has typically been limited to certain aspects of GPDs and TMDs, such as low x-weighted moments, together with their k 2 T -dependence for TMDs or t-dependence [8][9][10] for GPDs [11][12][13]. However, new approaches, such as Large-Momentum Effective Theory (LaMET) [14][15][16][17], now enable lattice QCD to reveal much richer information on GPDs and TMDs. Model * cshi@nuaa.edu.cn calculations are also crucial, as they can help provide an intuitive picture of the GPDs and TMDs. For instance, using the Nambu-Jona-Lasinio model or the spectral quark model one can calculate the full pion GPD over the entire kinematic range |x | < 1, |ξ | < 1 [18]. Such a calculation is based on non-perturbative covariant Feynman diagrams. An alternative approach is the light front QCD framework, where the GPDs and TMDs are determined through overlap representations in terms of light front wave functions (LFWFs) [19,20]. The unknown elements are then the non-perturbative LFWFs of hadrons. In this work, we will determine the LFWFs of the pion and kaon from a beyond rainbow-ladder Dyson-Schwinger equations (DSE) calculation, and then study the GPDs and TMDs obtained from these LFWFs using overlap representations. The pion and kaon are of particular interest as they emerge as the Goldstone bosons associated with dynamical chiral symmetry breaking (DCSB) in QCD. Pion and kaon GPDs and TMDs are also experimentally accessible-in principle-via hadron-hadron collisions using pion and kaon beams, or through interactions with the virtual meson cloud around nucleon targets [21,22].
To calculate the LFWFs the standard approach is to diagonalize the light-cone QCD Hamiltonian within the light-cone QCD formalism [23]. However, in practice this is numerically very difficult for exact QCD in four spacetime dimensions, and therefore effective interactions are usually adopted [24][25][26][27][28][29]. An alternative approach is to solve the covariant Bethe-Salpeter equation and project the Bethe-Salpeter wave functions onto the light front. This idea originates from a model calculation by 't Hooft [30] and was also used by authors in Refs. [31,32]. In a recent work the pion's leading Fock state LFWFs were obtained from its Bethe-Salpeter wave function provided by a DSE calculation [33], and used to study the pion TMD. In this paper we extend this work to the kaon LFWFs and TMDs, and also the study of pion and kaon GPDs.
In the past few decades the DSE framework has been applied extensively to hadron physics [34,35]. By solving the quark's gap equation and meson's Bethe-Salpeter equation, one obtains the covariant Bethe-Salpeter wave function, from which hadron properties can be determined. The DSEs respect the (approximate) chiral symmetry in the light quark sector, and its dynamical breaking, as demonstrated by satisfying the axial-vector Ward-Takahashi identity (AV-WTI) [36]. Therefore, the extracted LFWFs encode the effects from DCSB and provide a realistic description of the Goldstone bosons at leading-order in the Fock state expansion. Therefore, DSE predictions for the two-particle LFWFs of pion and kaon can provide important insights into the structure of QCD's Goldstone bosons.
This paper is organized as follows: In Sec. II we determine the LFWFs of pion and kaon from their Bethe-Salpeter wave functions. We then study their GPDs and related form factors in Sec. III and Sec. IV. The unpolarized TMDs are determined in Sec. V and a conclusion is given in Sec. VI.

II. PION AND KAON LIGHT FRONT WAVE FUNCTIONS
In light-front QCD hadron states are generally described by a tower of Fock states in a Fock state expansion [23,32]. For a meson with valence quark f and valence anti-quarkh the minimal (2-particle) Fock-state configuration is given by [26,29] where k T is the transverse momentum of the quark f [in a frame where the meson's transverse momentum vanishes (P T = 0)],k T = −k T , x = k + P + is the light-cone momentum fraction of the active quark, andx = 1 − x. The quark helicity is labelled by λ i = (↑, ↓) and δ i j / √ 3 is a color factor.
Ref. [37] showed that for pseudo-scalar mesons there are two independent light front wave functions for the leading Fock state, labeled by ψ 0 (x, k 2 T ) with l z = 0 and ψ 1 (x, k 2 T ) with |l z | = 1. The 2-particle Fock-state configuration is then given by where where k ± T = k 1 ± ik 2 . The LFWFs are obtained from the Bethe-Salpeter wave function via the light front projections [33,38] where the trace is over Dirac indices. The Bethe-Salpeter wave function is defined by χ fh (k, P) = ∫ d 4 z e −ik ·z 0|T f (z)h(0)|M(P) [39,40], and can be expressed as χ fh (k, P) = S f (k + P/2) Γ fh (k, P) S h (k − P/2), where S(k) is the dressed quark propagator and Γ(k, P) the meson's Bethe-Salpeter amplitude [34,41].
In the framework of the DSEs S(k) and Γ(k, P) are obtained by solving the quark gap equation and Bethe-Salpeter equation, respectively. For non-singlet pseudoscalar mesons the AV-WTI should be preserved by carefully selecting truncation schemes. The simplest symmetry-preserving DSE truncation is rainbow-ladder (RL) and has achieved many successes in the study of hadron properties [42][43][44]. A modern extension known as the DCSB-improved (DB) truncation improves upon the RL truncation and provides more realistic description of the pion, kaon, and other hadrons [35]. In this work we employ the existing DB-kernel solution parameterized in Refs. [45][46][47]. Further details about this DSE truncation are given in App. A.
To obtain the pion and kaon LFWFs we first determine an arbitrary k 2 T -dependent moment defined by These can be directly calculated using Eqs. (5) and (6), that is Since we have an analytical form for χ (k, P) obtained by parametrizing the numerical DSE solution, the twodimensional momentum integrations can be completed with the help of Feynman parametrization. In practice, we transform the integration variables to rewrite the integral in the form Comparison with Eq. (7) then reveals that the LFWFs are identified as ψ l z (x, . We present plots of the leading Fock state LFWFs for the pion and kaon in Fig. 1. For concreteness, we focus our discussion to the case of π − and K − , so the d and s are the valence quarks andū is valence anti-quark. In general we find that all the LFWFs are smooth functions decaying as k 2 T increases or x approaches the end-points. As expected for light mesons, the x-dependence of the LFWFs is broad at low k 2 T and get narrower as k 2 T increases, approaching an asymptotic form for large k 2 T proportional to x(1 − x). Fig. 2 provides an example of how the x-dependence of ψ 0 (x, k 2 T ) changes with k 2 T . The strong support of the LFWFs at infrared k 2 T originates from the strength of the covariant Bethe-Salpeter wave functions at low |k T |, which is closely connected to DCSB, as illustrated model-independently in Ref. [36]. Therefore, our LFWFs faithfully inherit the DCSB property from the covariant DSEs calculation. At large k 2 T , the LFWFs decay as ψ 0 (x, k 2 T ) ∼ 1/k 2 T and ψ 1 (x, k 2 T ) ∼ 1/k 4 T , in line with the perturbative QCD expectations [48]. The effects of SU(3) flavor symmetry breaking are clearly apparent in the kaon, as the heavier s quark gains more support at large x and the LFWFs become skewed. This indicates that the s quark carries more of the kaon's light-cone momentum fraction. However, these SU(3) flavor symmetry breaking effects diminish as k 2 T increases. Further analysis of these effects will be given in later sections when GPD and TMD results are presented.
The LFWFs are normalized so that the quark number sum rule ∫ 1 0 dx f (x; µ 0 ) = 1 is satisfied. Therefore, with only the leading Fock state the valence quark distribution function f (x; µ 0 ) is given by FIG. 1. The top row gives the LFWFs for pion and the bottom row gives the kaon results. The left column is ψ 0 (x, k 2 T ) and the right column is . This approximation to the full valence quark distribution function is best at a low hadronic scale µ 0 , which in Ref. [33] was determined to be µ 0 = 520 MeV. In a nonrelativistic system ψ 1 (x, k 2 T ) would vanish because the quarks are in a relative p-wave, however we find that the contribution to the quark number sum rule from ψ 1 (x, k 2 T ) equals 0.36 for the pion and 0.31 for the kaon. Therefore, we find that the valence quarks in both the pion and kaon are highly relativistic. Importantly, the relative strength between ψ 0 (x, k 2 T ) and ψ 1 (x, k 2 T ) in our approach is completely determined by the Bethe-Salpeter wave function, which itself is governed by the underlying quark-gluon interaction. The significant contribution of ψ 1 (x, k 2 T ) to observables likely also implies that higher Fock states may not be negligible in a more realistic calculation. Nevertheless, the higher Fock states are much more difficult to calculate and are beyond the scope of this work.
To calculate H q M (x, ξ , t) we employ its light front overlap representation. This result can be obtained using light-cone quantization and expanding the quark field in Eq. (12) using the canonical field mode expansion and the hadron state ket using a Fock state expansion. Contracting all the operators, one gets the light front overlap representation of the GPD in terms of LFWFs. However, in the ERBL region this requires the overlap of LFWFs with different numbers of constituents, i.e., N and N + 2. The ERBL region is therefore inaccessible is a leading Fock state expansion due to the lack of a 4-particle Fock state.
In a meson M with active quark f , the GPD H f M (x, ξ , t) in the DGLAP region can be expressed as the overlap of LFWFs [19,20,38,49] For the active anti-quarkh, the GPD can be obtained analogously as [20] In the absence of an accessible ERBL region we limit our study to zero skewness, Since the GPD reduces to the PDF at zero momentum transfer, i.e., H q (x, 0, 0) = f q (x). The initial scale µ 0 is determined as follows (see Ref. [33]): At the scale of Q 2 = 4 GeV 2 , the π N Drell-Yan analysis gives averaged momentum fraction of valence quark distribution in pion as 2 x = 0.47 (2) [50,51] and the lattice QCD gives 2 x = 0.48(4) [52]. To match this result, we determine µ 0 = 0.52 GeV, so that x = 0.5 at µ 0 reduces to x = 0.24 at 2 GeV by NLO DGLAP evolution.
The two-dimensional Fourier transform of H q M (x, 0, ∆ 2 T ) gives the IPDs: The IPDs have the interpretation of parton distributions in the transverse plane [5,6], with x the light-cone momentum fraction and b T the transverse separation between the active parton and the origin of transverse center of momentum R T . In the valence picture with two constituents, Fig. 4 we plot ρ q M (x, b 2 T ) for pion and kaon. An important observation is that as x becomes larger, the width of the curves shrinks and the quark distributions are more spatially localized. When x → 1, the width is vanishingly small and the quark stays near the center of transverse momentum. This can be understood since when one quark carries almost all of the light-cone momentum (as x → 1), then R T → r T ,1 and b T ,1 → 0, namely, this quark defines the transverse center of momentum. Alternatively, if we consider the overlap representation of ρ(x, b 2 T ) in terms of LFWFs in the coordinate space, that is [29,53] ρ whereΦ Flavor symmetry breaking effects are clearly evident in Fig. 4. Typically, at smaller x (x = 0.3) there is moreū quark than s quark in kaon over the whole b T range. At larger x (x = 0.7) the situation is reversed. This suggests the s quark is more likely distributed near the center of kaon while the u quark is more spread out. We can also look at 1 Recall that Φ λ 1 λ 2 (x, k T ) has been defined in Eq. (1) and can be easily be related to ψ 0 (x, k 2 T ) and ψ 1 (x, k 2 T ) via comparison with Eqs. (2)-(4). which characterizes the quark density at transverse separation b T . As shown in the lower panel of Fig. 4, the s quark in kaon favors small b T andū quark has a broader distribution, with d quark in pion lying in between. If we look at their mean-squared b T , i.e., It's worth mentioning here that in our calculation, the current quark mass we used are m MeV. This big mass difference gets weakened by the DCSB, and the difference in the u/d and s quark distributions is no longer so dramatic.
Further, on can define the valence-like distribution Fig. 4. However, it's worth mentioning that ρ (0) (b 2 T ) is independent of the renormalization scale, because DGLAP evolution conserves the quark number density at every slice of b T . Equivalently, H (x, 0, t) evolves independently of t [5]. Thus the lower panel of Fig. 4 can also be viewed as the valence (anti-)quark spatial distribution at any scale.

IV. ELECTROMAGNETIC AND GRAVITATIONAL FORM FACTORS
The electromagnetic form factors of a hadron provide important information about its spatial structure. The pion and kaon have one electromagnetic form factor defined by The pion and kaon electromagnetic form factors are also given by the lowest x-weighted moment of their GPDs which is independent of skewness ξ because of the polynomiality property of the GPDs. The result for the pion's electromagnetic form factor obtained using Eq. (19) is given by the dashed curve in Fig. 5. In general we find that our result overshoots the data for all Q 2 , and also the full DSE calculation that uses the Bethe-Salpeter wave function and a dressed quark-photon vertex to directly calculate the pion's form factor [44]. As we will explain, the origin of these discrepancies is naturally explained by the Fock state truncation and the LFWF normalization condition [see Eq. (11)]. At (very) large Q 2 perturbative QCD predicts that the Electromagnetic form factor F (t) of pion in the space-like region. The data is the from NA7 Collaboration [54] (red empty circle) and Jefferson Lab [55] (green filled square). The dashed (blue) curve is based on the unmodified GPD in Eq. (13), while the solid (black) curve uses a GPD with a dressed operator to simulate higher Fock states, see Eq. (24). The dotted curve is the full rainbow-ladder DSE result from Ref. [44] that including an infinite tower of Fock states.
pion's electromagnetic form factor behaves as [56] where w π is the x −1 moment of parton distribution ampli- , where in this case the PDA ϕ π (x, Q 2 ) is normalized at the scale of Q 2 such that where f π = 92.4 MeV is the pion's electroweak decay constant. The DSE calculation based on Eqs. (20)- (22) has been presented in Ref. [45] and the result is reasonable. However, the LFWF normalized by Eq. (22)  The deviation in the low Q 2 region is also easy to understand. The normalized condition for the LFWFs is such that F π (0) = 1. However, as mentioned higher Fock states will carry some charge, which, if included, would cause a modification to the form factor at low to intermediate Q 2 . In addition, there are important contributions that can dramatically change the charge radius but do not impact the charge. Traditionally, these are associated with vector meson dominance (VMD) contributions. VMD is associated with meson poles in the time-like region, where for the pion electromagnetic form factor the rho pole is the most important. In the LFWF approach these VMD contributions can only be obtained by including an infinite tower of Fock states. This is natural in the complete DSE calculation with a dressed quark-photon vertex, but very challenging in a rigorous light-front approach. It is therefore not possible for a leading Fock state calculation-that is intimately connected to underlying QCD dynamics-to give a good description of the electromagnetic form factor for all Q 2 .
With the pion's (Breit-frame) charge radius defined by we obtain from the leading Fock state calculation r c = 0.41 fm, which is significantly smaller than the experiment value of r c = 0.67 fm [57]. A similar result was also found using a relativistic constituent quark model based on an effective qq Hamiltonian [58], where a pion charge radius of r c = 0.45 fm was found [59]. In this work the authors argue that the discrepancy with experiment can be corrected by taking into account the constituent quark charge radius, which is analogous to dressing the vertex as in a full DSE calculation. In a complete DSE calculation the operator that defines the GPDs would be dressed. Such a calculation from the DSE is very difficult and beyond the scope of this work. However, we can use an analogous calculation for this dressed operator from the NJL model to obtain a qualitative measure of the impact of a dressed vertex, or equivalently higher Fock states. Using a dressed operator that defines the GPD from the NJL model [60], we find in the impulse approximation that our leading Fock state DSE result is modified such that where and the modified GPD at zero skewness is denoted by  H d (x, 0, t). Note, in the pion H u (x, 0, t) can be obtained from H d (x, 0, t) by charge symmetry. The second term on the right hand side of Eq. (24) comes from the dressing of the quark vertex in the impulse approximation and provides an additional contribution (see App. B for details).
Using H (x, 0, t), we get the solid curve in Fig. 5 and a charge radius r c = 0.59 fm, with the low to intermediate −t region also significantly improved.
The modification term in Eq.(24) has many interesting properties. For instance, its dressing functioñ F ρ (t) vanishes at t = 0, so the PDF is unchanged, i.e., H (x, 0, 0) = H (x, 0, 0). While at non-vanishing t, the modification term proportional to δ (x) is infinitely negative. Its integration over x yields a finite suppression to the electromagnetic form factor. In terms of the overlap representation, this correction can only be obtained by including an infinite tower of Fock states containingqq pairs. The modification in Eq. (24) brings no change to ρ(x, b 2 T ) for x > 0 and all the results in last section still hold.
The higher moments of the GPD at ξ = 0, i.e., are not affected by this modification term, as a consequence of the δ (x). Among these moments, A q 2,0 (t) contributes partially to the pion's gravitational form factor Θ 2 (t), defined through the matrix element of energymomentum tensor for one-pion states [61] π + (p )|Θ µν (0)|π . (27) with P = p + p , q = p − p and t = q 2 . The form factor Θ 2 (t) is scale independent, while its individual quark contributions A q 2,0 (t) evolve with scale. At the low model scale, the valence picture gives Θ 2 (t) = q A q 2,0 (t). As the scale increases, A q 2,0 (t; µ) evolves accordingly to the evolution of the GPD.
In Fig. 6 we show pion's A d ;π 2,0 (t) (solid red curve) at the scale of 2 GeV and the curve lies within the lattice simulation data. It is closer to the NJL model result (blue dashed) [62] than to the spectral quark model [62]. We have illustrated the kaon GFFs Aū ;K 2,0 (t) and A s;K 2,0 (t) as well.
A light-cone energy radius can be defined in relation to the gravitational form factor A 2,0 (t), and is given by [64] which can be contrasted with an analogous light-cone charge radius defined by r 2 c,LC = −4 ∂ F (Q 2 )/∂Q 2 Q 2 =0 . For the pion we find r u, π c,LC = 0.331 fm and r u, π E,LC = 0.185, meaning the energy radius is about 56% smaller that the light-cone charge radius. Both these radii will be impacted by higher Fock states, however, based on vector meson dominance the light-cone charge radius will increase more because it is impacted by the ρ meson pole whereas the light-cone energy radius is impacted by spin-2 mesons which are much heavier and further from Q 2 = 0. Therefore, we predict that r u, π E,LC /r u, π c,LC = 0.56 is an upper bound on this ratio. For the kaon we find light-cone charge radii   The unpolarized leading-twist TMD is defined as with the gauge link omitted. In terms of the leading Fock state LFWFs the TMD reads [65] f which should be associated with an initial result at a low scale of µ 0 = 520 MeV, just as in the GPD case. We plot the unpolarized TMD for the pion and kaon (s quark) in Fig. 7. Theū TMD in kaon can be simply obtained from s quark distribution by momentum conservation, i.e., The distribution closely resembles the profile of LFWFs. Which suggests that in a purely valence quark picture the quarks are most likely to carry around half of the parent hadron's light-cone momentum with a small intrinsic transverse momentum.
The transverse momentum dependence of the TMD has long been of great interest, and in Fig. 8 we illustrate our results for fixed values of x. Our results decrease with increasing |k T |, being concave at low |k T | and becoming convex as |k T | increases. The inflection point is around 300 MeV. Phenomenologically, postulating a Gaussian |k T |-dependence is popular and Gaussian-based models successfully describe much of the existing data [66][67][68][69][70][71][72]. The gray dot-dash-dash curve is a Gaussian function employed to fit f d 1;π (x = 0.3, k 2 T ) at low |k T |. One can see the fit is good up to around 300 MeV, i.e., it describes well the intrinsic transverse momentum dependence. For large |k T | the Gaussian form would inevitably fail since f q 1 (x, k 2 T ) ∼ 1/k 4 T with our LFWFs.
The x-dependence and k T -dependence in our TMDs are not factorizable, except for at very large k 2 T . For instance, in Fig. 8, the k 2 T is respectively 0.14 GeV 2 and 0.13 GeV 2 when fitting f d In recent years, phenomenological studies of TMDs have appreciated the xdependence of the |k T | behavior and build it this into their parameterizations at the low initial scale [72,73]. In this sense, our result shows qualitative agreement. Note that TMD evolution also generates significant x-dependence in the |k T | behavior, as has been shown in Refs. [33,74]. Finally, we report that for the s quark in kaon, k 2 T is respectively 0.155 GeV 2 and 0.134 GeV 2 when fitting f s 1;K (x = 0.3, k 2 T ) and f s 1;K (x = 0.5, k 2 T ) to the Gaussian form. This is slightly larger than theū or d quarks in pion, and is a measure of SU(3) flavor symmetry breaking.

VI. CONCLUSION
By projecting the mesons' covariant Bethe-Salpeter wave functions onto the light front, we calculate the leading Fock state LFWFs of the pion and kaon. The kaon's LFWFs based on a DSE approach are given for the first time. These LFWFs are significantly enhanced at low |k T | and exhibit the perturbative QCD power law behavior at large |k T |. SU(3) flavor symmetry breaking is revealed in the kaon's LFWFs. We also observe a sizable contribution from the spin-parallel LFWF, suggesting an important role played by the p-wave component in the pion and kaon as relativistic composite particles.
We employ the light front overlap representation given in Eq. (13) to study the GPDs at zero skewness H (x, ξ = 0, t) for the pion and kaon, and using the IPDs ρ(x, b 2 T ), we determine the spatial distribution of the valence quarks. On the light front the quarks with larger light-cone momentum fraction x are generally less spread out in the spatial impact parameter b T . After integration over x, we find the heavier quarks are more concentrated at the center of meson, e.g., the s quark spatial distribution in K − is narrower than theū quark.
Shortcomings in a leading Fock state truncation are exposed in the pion's electromagnetic form factor. An attempt is made to overcome these issues using a dressing of the operator that defines the GPD obtained from the NJL model. The pion and kaon electromagnetic and gravitational form factor then show reasonable agreement with available experimental and lattice data.
Finally, we give the unpolarized TMDs of pion and kaon. The phenomenologically popular Gaussian-like k Tdependence is observed in our result for intrinsic |k T |, but violated at medium and large |k T |. It is also observed that the |k T | behavior in our TMD is slightly x-dependent, suggesting an unfactorizable x-and k T -dependence. In addition, the valence quarks in kaon have a broader transverse momentum distribution, as a consequence of SU (3) flavor symmetry breaking. Starting with the DSE Bethe-Salpeter wave function we have therefore obtain comprehensive insights into the pion and kaon valence quark imaging in both position and momentum space. To aid the calculation of the moments x m l z we use an accurate parametrization of numerical solutions to the gap and BSEs in the DCSB-improved truncation to the DSEs [45,75]. Solutions of the DSE-BSE with the so-called DCSB-improved kernel are available within the literature, both for the pion and the kaon [45,46]. In this work, we employ these results and their available parameterization, that we remind the reader of and slightly modify. The quark propagator S(k) is written as the sum of pairs of complex conjugate poles:   with n = 2. The pseudo-scalar Bethe-Salpeter amplitude Γ(k; P) can be generally decomposed as Γ(k; P) = γ 5 iE(k; P) + / PF (k; P) We employ the dominant terms E(k; P) and F (k; P), which are parameterized by: where ρ u (α) = 3 4 (1 − α 2 ) and {C (1/2) n , n = 0, 1, ..., ∞} are the Gegenbauer polynomials of order 1/2. The value of the parameters are listed in Tab. I. The outgoing quark and anti-quark in the meson carry momentum k + P/2 and k − P/2 respectively, so F (k; P) is even in k · P due to charge parity.   calculated as with momentum assigned in Fig. 9. Here we consider the GPD of isospin 0 or 1, defined as H I =0 = H u + H d and H I =1 = H u − H d . Flavor matrices are implicitly embedded in the elements S, Γ π and Γ ·n. The notation Γ ·n represents the dash line boxed area. We denote the Γ · n as a violet blob in other diagrams. It satisfies the inhomogeneous Bethe-Salpeter equation: Here the red blob is the bare vertex Note the first Dirac delta ensures k ± ∆/2 = l ± ∆/2 at leading truncation. The matrices τ 0 or τ 3 are for isospin 0 or 1, respectively. To solve for Γ · n, one can formally sum the series Its sub-leading term can be evaluated via Mellin moments, i.e., where Ω denotes any of the five Dirac/isospin structures appearing in the NJL model Lagrangian. 2 At ξ = 0 one finds (P · n) ∫ dx x s B I =0,1 (x, 0, t) = −2G ω, ρ Π V V (t)/ n : s = 0 0 : s ≥ 1 (B8) which uniquely determines the result for B I to be: One can calculate the rest terms analogously and their summation gives the overall dressed quark correlator: Putting Eq. (B10) back into Eq. (B1), the first term in the braces gives the bare vertex contribution The second term has x dependence factored out and one easily finds its contribution is proportional to the lowest moment of Eq. (B11). Finally we have Namely, the dressing of the bare vertex introduces an additional term proportional to δ (x), leading back to the modified GPD in Eq. (24). Used functions in Eq. (24) arẽ Here the proper time regularization is used, with parameters determined by hadron mass spectrum and decay constant from Table 1 in Ref. [76], i.e., Λ I R = 0.24 GeV, Λ U V = 0.645 GeV, M = 0.4 GeV and G ρ = 11.0 GeV −2 . Finally, we note that the above dressing diagram doesn't modify the unpolarized TMD within the NJL model. The easiest way to see this is by realizing that in the NJL model the unpolarized PDF is obtained by integrating out the transverse momentum of TMD, i.e., f (x) = ∫ dk 2 T f 1 (x, k 2 T ). Since f (x) = H (x, 0, 0) receives no contribution from the dressing diagrams, and f 1 (x, k 2 T ) is always positive, one deduces any corrections to f 1 (x, k 2 T ) are zero.