Contact interaction analysis of pion GTMDs

A contact interaction is used to calculate an array of pion twist-two, -three and -four generalised transverse light-front momentum dependent parton distribution functions (GTMDs). Despite the interaction’s simplicity, many of the results are physically relevant, amongst them a statement that GTMD size and shape are largely prescribed by the scale of emergent hadronic mass. Moreover, proceeding from GTMDs to generalised parton distributions, it is found that the pion’s mass distribution form factor is harder than its electromagnetic form factor, which is harder than the gravitational pressure distribution form factor; the pressure in the neighbourhood of the pion’s core is commensurate with that at the centre of a neutron star; the shear pressure is maximal when confinement forces become dominant within the pion; and the spatial distribution of transversely polarised quarks within the pion is asymmetric. Regarding transverse momentum dependent distribution functions, their magnitude and domain of material support decrease with increasing twist. The simplest Wigner distribution associated with the pion’s twist-two dressed-quark GTMD is sharply peaked on the kinematic domain associated with valence-quark dominance; has a domain of negative support; and broadens as the transverse position variable increases in magnitude.


Introduction
It is anticipated that an electron ion collider will be operating in the USA by 2030 [1,2]; construction of a similar machine is being discussed in China [3,4]; new capabilities are expected at Conseil Européen pour la Recherche Nucléaire (CERN) [5]; and the Jefferson Laboratory (JLab) is currently operating at 12 GeV [6]. Each of these facilities has given high a e-mail: jlzhang@njnu.edu.cn b e-mail: phycui@nju.edu.cn c e-mail: jlping@njnu.edu.cn d e-mail: cdroberts@nju.edu.cn (corresponding author) priority to experiments that can yield data that may be used to draw three-dimensional (3D) images of hadrons, i.e. measurements interpretable in terms of generalised or transverse momentum dependent parton distributions: GPDs or TMDs, respectively.
Hadron physics has long focused on one dimensional (1D) imaging of hadrons. It is an ongoing effort, which remains crucial because many puzzles and controversies are unresolved. For instance, even considering what may seem to be the simplest strong interaction system, the pion valencequark distribution has been studied for roughly thirty years, both experimentally [7][8][9][10][11] and theoretically [12][13][14][15][16][17][18][19]; yet, it still attracts vigorous debate. Moreover, the pion's glue and sea distributions are empirically unknown, with theoretical predictions only now becoming available; and kaon distributions are just beginning to receive renewed attention [20][21][22]. The challenge of producing solid predictions for parton distributions within baryons is even greater.
Notwithstanding the need for new, precise data on 1D distributions and associated predictions with a traceable connection to quantum chromodynamics (QCD), the allure of GPDs and TMDs is difficult to resist, given that 3D imaging may enable entirely new aspects of hadron structure to be revealed. Such functions serve as tools with which to probe the multidimensional structure of hadron light-front wave functions (LFWF), thereby providing access to, inter alia: the distributions of mass, pressure and spin within a hadron, both in longitudinal and transverse directions; the sharing of these qualities amongst the various bound-state constituents; and to the spacetime volumes occupied by these constituents, i.e. to their potentially different "confinement" radii.
In order to fully capitalise on 3D imaging data obtained at modern and anticipated facilities, using it to understand the many correlated phenomena which emerge from strong interactions in QCD, methods must be developed that enable GPDs and TMDs to be calculated within frameworks that are mathematically linked to the fundamental theory. To see the importance of this, one need look no further than the thirty year controversy over the pion's valence quark distribution function [21][22][23][24][25][26][27][28][29][30][31].
Herein we explore and illustrate the capacity of generalised parton correlation functions (GPCFs) [32] to serve as a framework for the unified calculation of GPDs and TMDs. As this is a first step, we choose to study the pion and work with a confining, symmetry-preserving treatment of a vector × vector contact interaction (CI) as the foundation for our analysis [33]. A merit of this approach is that, by enabling a largely algebraic treatment of relevant processes and quantities, it provides for an insightful assessment of all results. Moreover, when considered judiciously [33][34][35][36][37][38][39][40], such results may often be interpreted from a QCD perspective because this treatment of the CI preserves many qualities of the leading-order truncation of QCD's Dyson-Schwinger equations (DSEs), itself a sound approach to many hadron observables [41][42][43][44][45][46][47].
Our analysis begins in Sect. 2 with a brief review of the GPCF for a J = 0 hadron. Section 3, augmented by "Appendix A", then describes our CI treatment of the pion and its coupling to photons. The pion GPCF is used in Sect. 4 as the basis for calculating the four twist-two generalised transverse momentum dependent parton distribution functions (GTMDs) associated with dressed-quarks within the CI pion. The discussion highlights the role played by emergent hadronic mass (EHM) in determining the properties of each GTMD. (CI results for all twist-three and twist-four GTMDs are provided in "Appendices B, C", respectively.) In Sect. 5, the twist-two GTMDs are integrated over their light-fronttransverse momentum argument, k 2 ⊥ , to yield results for the pion's vector and tensor GPDs. Features of the derived electromagnetic, gravitational, and transverse-spin distributions are also canvassed. Section 6 shows how one proceeds from GTMDs to TMDs. It provides explicit formulae for all four TMDs supported by the CI in the absence of an adequate model for the Wilson line and highlights their relative sizes and domains of material k 2 ⊥ -support. Section 7 emphasises and illustrates the connection between GPCFs and Wigner distributions by presenting the CI result for a Wigner distribution associated with pion twist-two GPDs and TMDs. A summary and perspective is provided in Sect. 8.

Generalised parton correlation function
We begin by considering the following in-pion quark-quark correlator [32]: where: and k is the relative quark-antiquark momentum. These conventions are illustrated in Fig. 1.
The hitherto undefined quantity in Eq. (1) is the Wilson line, W (− 1 2 z, 1 2 z;n), wheren is a light-like four-vector,n 2 = 0, antiparallel to P,n · P = P − , and the path is chosen as a sequence of line segments [32,48]: The same path is achieved by rescalingn → λn, λ ∈ R, λ > 0; hence, withP 2 = 1, Eq. (1) only depends on The quantity η in Eq. (1) expresses the one remaining degree of freedom, viz. η = sign(n 0 ), in which case η = ±1 describe, respectively, future and past Wilson line trajectories. One passes to generalised transverse-momentum dependent parton distribution functions (GTMDs) by first considering the following partially integrated quantity: where n is a light-like four-vector for which n · P = P + . The object in Eq. (5) is a Dirac-matrix valued function and, as usual, contributions at various orders in a twist expansion can be obtained by appropriate projection operations. Namely, with H being some suitably chosen combination of Dirac matrices, then the scalar functions of interest -the GTMDs -are obtained via Referring to Fig. 1, this operation corresponds to the insertion of H as a connection between the open quark and antiquark lines: ψ(k ∓ /2), respectively. As defined by Eq. (6), a given GTMD is a complex-valued function: the real part is even under the time-reversal operation (T -even), whereas the imaginary part is T -odd. Equally, they are even (odd) under η → −η. (Recall η = ±1 specifies the time-direction of the Wilson line used to define the GTMD.) Beginning with Eq. (6), GPDs are obtained by integration over k ⊥ : only the T -even piece survives, which is independent of η; and the array of TMDs is obtained by setting = 0, which entails ξ = 0.
To proceed, we follow Ref. [40]; namely, retaining m G = 0.5 GeV but setting α 0 /π = 0.36. This combination ensures a good description of π -meson properties. Furthermore, since a momentum-independent interaction cannot support relative momentum between bound-state constituents, we simplify the tensor structure in Eqs. (7), defining the CI RL kernel as follows: When using Eq. (9) in a DSE, it is necessary to impose an ultraviolet regularisation scheme. It should be symmetry preserving so that the results maintain a meaningful connection with the Standard Model. Moreover, since a CI does not produce a renormalisable theory, the associated regularisation mass-scale, uv , is an additional physical parameter. It may be interpreted as an upper bound on the momentum domain within which the properties of the associated system are practically momentum-independent.
As the final step in defining the CI, we include an infrared regularisation scale, ir , when computing all integrals connected with bound-state problems [68]. Since chiral symmetry is dynamically broken by Eq. (9), ensuring the absence of infrared divergences, ir is not a necessary part of the CI's definition. Notwithstanding that, by excising momenta k < ir , one achieves a rudimentary expression of confinement via elimination of quark production thresholds [67,[69][70][71][72][73][74][75]. A natural choice for this scale is ir ∼ QCD . We set ir = 0.24 GeV. Assuming isospin symmetry, it only remains to fix the current-mass, m, of the light quarks. That may be achieved by solving the pion bound state problem specified by the kernel in Eq. (9). In this case, the gap equation for the dressed light-quark propagator is The integral is quadratically divergent. When it is regularised in a Poincaré-invariant manner, the gap equation solution is where M is the dressed-quark mass, momentum-independent in the CI, determined by We define the regularised integral by writing [68] 1 where τ ir,uv = 1/ ir,uv are, respectively, the infrared and ultraviolet regulators described above. Consequently, the gap equation becomes where with (α, y) being the incomplete gamma-function.
In an internally consistent treatment of a vector × vector CI, the Bethe-Salpeter amplitude for the π -meson has the following form [33,35,36]: Here, Q is the pion's total momentum, Q 2 = −m 2 π , m π is the pion mass; M is obtained from the contact-interaction gap equation, Eq. (14); and E π , F π do not depend on the relative quark-antiquark momentum.
The amplitude, π , is obtained from the following homogeneous Bethe-Salpeter equation: Employing the symmetry-preserving regularisation scheme of Refs. [33,36], which emulates dimensional regularisation and requires where C 1 is given in Eqs. (A.1), (A.2) and α = 1 − α, one arrives at the following pair of coupled equations: .4). Evidently, the kernel is only defined after the gap equation has been solved.
Inspection of Eqs. (20), (A.4) reveals that a nonzero value for E π enforces F π = 0, i.e. any theory with a traceable connection to a vector-boson exchange interaction must retain both E π , F π . (When the interaction is momentum dependent, then two other amplitudes are also nonzero [76,77].) If Table 1 With input parameters [35,40]  F π is omitted, then one arrives at a model, which although it may be useful for parametrising data, cannot contribute to the development of insights into characteristics of the Standard Model's Nambu-Goldstone modes [34,36]. Equation (20) is an eigenvalue problem. It has a solution when Q 2 = −m 2 π , at which point the eigenvector is the meson's Bethe-Salpeter amplitude. Working with the on-shell solution, normalised canonically according to Eqs. (A.5), (A.6), the pion's leptonic decay constant is given by (N c = 3): In the chiral limit, i.e. using solutions obtained with m = 0 in Eq. (10), this reduces to [33] Solving Eqs. (10), (20), we obtain the results listed in Table 1, reproducing those reported elsewhere [35,40].
For subsequent use, here we also introduce the dressed photon-quark vertex, γ μ . Using Eq. (9), one has ± = ± /2 (23) Owing to the vector Ward-Green-Takahashi identity (WGTI), preserved in our regularisation of the contact interaction, the solution takes the form [34] where As expected of RL truncation studies of the photon-quark vertex [78,79], the dressing function, P T ( 2 ), exhibits a simple pole at 2 = −m 2 ρ , where m ρ is the mass of the ρ-meson that is generated by the interaction.

Pion twist-two GTMDs
There are three twist-two pion GTMDs. They are obtained with the following choices in Eq. (6): The simplest is that associated with H 1 , which relates to the pion valence-quark distribution function and electromagnetic form factor. We therefore use it to illustrate the computational techniques.
Mapping into Euclidean metric: and since a RL truncation was used to solve the Bethe-Salpeter equation, then internal consistency and preservation of symmetries requires a kindred truncation for the GTMD, in which case where tr D indicates a trace over spinor indices, δ x n (k) = δ(n · k − xn · P), and the "skewness" ξ = [−n · ]/[2n · P], |ξ | ≤ 1.
In proceeding with a WGTI-preserving evaluation of Eq. (28), we first compute the spinor trace; then using the following identities [D(k 2 ) = k 2 + M 2 ]: cancel each common numerator and denominator factor; and finally use Feynman parametrisations to simplify all remaining denominators. In this way, one arrives at : and For later use, we note that one can write Here it is worth recalling a Goldberger-Treiman relation that emerges in a WGTI-preserving treatment of the CI. Namely, in the absence of a Higgs mechanism -so that m = 0 in the gap equation, Eq. (10), and one is dealing with the chiral limit [33]: where the superscript "0" indicates evaluation in the chirallimit. Both M 0 and f 0 π are order parameters for dynamical chiral symmetry breaking (DCSB) [79], which itself is an expression of EHM in the Standard Model [45]. Moreover, Eq. (35) is practically unchanged at physical lightquark current masses. (Similar statements also hold in QCD [93,94].) Consequently, the strength of the pion's canonically normalised Bethe-Salpeter amplitude is a direct measure of EHM; hence, the CI formulae presented above and those to follow reveal that the size and shape of every one of the pion's GTMDs are largely determined by the character of EHM.

Algebraic results
As noted in closing Sect. 2, one proceeds from a GTMD to a GPD by integrating over k ⊥ ; and focusing first on the leading twist GTMDs, one therefrom obtains two GPDs: where H π , E T π may respectively be called the vector (no spin-flip) and tensor (spin-flip) GPDs. The former is directly related to the pion's elastic electromagnetic form factor and gravitational form factors (mass and pressure/stress) [95], whereas the latter provides access to the dependence of the pion's quark distributions on their polarisation perpendicular to the pion's direction of motion (transversity) [96,97]. where with Using the results following Eq. (34) and Eqs. (A.8), it is straightforward to establish that i.e. our CI treatment preserves the time-reversal-invariance property of the GPD. It is nevertheless deficient on the domain E = {x | − ξ < x < ξ} because H π (x, ξ, t) does not satisfy the soft pion theorems [98] A remedy is described elsewhere [81]; to wit, one must include interactions between the two pions in Fig. 1 that would lead to formation of a scalar meson-resonance. Profiting from this understanding, we expand on the Ansatz in Ref. [99] and correct the twist-two vector GPD: where P σ (t) is a quark+antiquark scalar-channel analogue of is consistent with known mathematical GPD constraints and Eqs. (45). 2 Inserting Eq. (39a) into Eq. (40b), the twist-two tensor GPD is obtained: The following remarks are pertinent: ME T π (x, ξ, t) is dimensionless; relative to some other studies, e.g. Refs. [100,101], our normalisation convention in Eq. (39a) entails that E T π (x, ξ, 0) is nonzero in the chiral limit; and once again using the results described in connection with Eqs. (34), (A.8), one finds
at the pion mass in Table 1 (As we have definedH π , this is only approximately true; but if necessary, that is readily corrected following the procedure in footnote 1. the CI dressed-quark distribution amplitude. (d) Using a contact interaction, the 2 The factor ξ 2 in Eq. (46c) should strictly be replaced by θξ ξ / p(ξ 2 ), where p(ξ 2 ) is a simple polynomial, chosen to preserve GPD polynomiality; but that merely complicates numerical analysis without delivering practical improvement. GPD is continuous but not differentiable at x = ±ξ . (This is typical of models whose basis is a separable interaction [81,102].) Beginning with H π , the pion elastic electromagnetic form factor is obtained via It is readily verified by straightforward calculation that the evaluated integral is independent of ξ .
The computed pion form factor is depicted in Fig. 3 upper-panel as the dashed red curve, from which one obtains the associated radius: r em π = 0.44 fm. As discussed in detail elsewhere [33,35,36], the WGTI-preserving treatment of a CI necessarily generates F π = 0 in Eq. (16). Consequently, the CI form factor is hard; 3 namely, it approaches a nonzero constant value as Q 2 → ∞.
It is appropriate now to consider the CI pion vector GPD in impact parameter space [103]: where J 0 is a Bessel function. This density describes the probability of finding a dressed-quark within the lightfront at a transverse position b ⊥ from the pion's centre of transverse momentum (CoTM). Inspecting Eqs. (41) . Consequently, if one omits F π in Eq. (16) so that the pion's elastic electromagnetic form factor is a monopole characterised by a length-scale, where K 0 is a modified Bessel function of the second kind.
Returning to an internally consistent WGTI-preserving CI treatment, so that F π = 0, then the large-Q 2 behaviour of the pion form factor may be characterised via M F → ∞; hence, We have verified these statements numerically. The n = 1 Mellin moment of the twist-two vector GPD delivers the pion's gravitational form factors: where θ 2 relates to the quark mass distribution within the pion and θ 1 is linked to the quark pressure distribution. In a symmetry preserving treatment: θ π 2 (0) = 1; and, following from Eqs. (45), θ π 4 The pion's gravitational form factors are also drawn in Fig. 3. Regarding θ π 1 , the Ansatz used to correct H π on the domain E, Eqs. (46), depends on a representation of the σ -resonance contribution to quark+quark scattering in the scalar channel. To illustrate the associated model-dependent -orange longdashed curve and band. The bands enclose the envelope of curves that fit the lQCD results [105] uncertainty, we used two forms: where m f 0 ≈ (0.48 − i0.28) GeV [106]. The first choice is based on the observation that the CI produces a σ -meson with mass m σ ≈ 2M in the neighbourhood of the chiral limit [35], whereas the second uses instead the pole mass associated with the empirical σ -resonance. Evidently, the uncertainty is noticeable but not large. We find r π θ 1 /r em π = 1.88(13); and a result that is generally softer than the pion's electromagnetic form factor. Turning to θ π 2 , r π θ 2 /r em π = 0.89; and this form factor is generally harder than F em π ( 2 ). The lower panel of Fig. 3 displays a comparison between our CI results and those obtained using lattice-QCD (lQCD), described by [105]: where 2E( ) = 4m 2 π + 2 . The physical interpretation of such distributions is complicated by issues connected with the Poincaré transformation of frame-dependent wave functions in quantum field theory [107]. Nevertheless, they are mathematically well defined; do admit the standard interpretation in systems for which a nonrelativistic approximation can be discussed; and viewed judiciously, can deliver fruitful insights. Moreover, two-dimensional Fourier-transform analogues deliver results of similar magnitude.
Owing to the hardness of CI pion form factors, the integrals in Eqs. (56) do not converge when evaluated using the results for θ π 1,2 ( 2 ) depicted in Fig. 3 -upper panel. We therefore exploit the semiquantitative similarity between CI and lQCD results evident Fig. 3 -lower panel to justify an estimate of the pion's pressure distribution using Eq. (55). The result is depicted in Fig. 4 and the qualitative features are consistent with an intuitive physical interpretation. Namely, the pressure is large and positive in the neighbourhood r 0the dressed-quark+dressed-antiquark are pushing away from each other at small separation; but the pressure changes sign as the separation becomes large, signalling a transition into the domain whereupon the pair experience the effects of confinement forces.
It is important to appreciate that lim r →0 r 2 p π (r ) = 0 in Fig. 4 is an artefact of the simple monopole description of θ 1 ( 2 ) in Eq. (55). In four spacetime dimensions, a quantum field theoretical treatment of form factors always introduces scaling violations, leading to additional ln( 2 /M 2 ) suppression on 2 M 2 . We choose to illustrate the effect of such scaling violation by modifying Eq. (55) as follows: Using this form for θ 1 leads to the blue dot-dashed curve in Fig. 4. In this case, lim r →0 r 2 p π (r ) = 0; yet, the characterising magnitudes are unchanged. An analogue of Eq. (56) has been used to infer the proton's quark pressure distribution from existing data on deeply virtual Compton scattering [108]. Comparing that result with those in Fig. 4 -upper panel, one observes that: (i) the pressure within the pion on the neighbourhood r 0 is roughly five-times larger than that in the proton; and (ii) the two pressure profiles have a similar radial extent. Notwithstanding the issues with Ref. [108] canvassed in Refs. [109,110], profiles analogous to Fig. 4 -upper panel for neutron stars indicate r 0 pressures therein of roughly 0.1 GeV/fm [111]; hence, the core pressures in the pion and neutron stars are commensurate.
A shear pressure distribution can also be defined [97]: whereˆ 2 = 1 =r 2 and j 2 is a spherical Bessel function. Intuitively, r 2 s π (r ) provides an indication of the strength of QCD forces within the pion which act to deform it. Our results are drawn in Fig. 4 -lower panel. Focusing on the more realistic curve, obtained using Eq. (57), these forces peak in the neighbourhood upon which the normal pressure switches sign, i.e. where the forces driving the quark and antiquark away from the core are overwhelmed by attractive confinement pressure. The twist-two tensor GPD expressed in Eqs. (47) is drawn in Fig. 5: it is only nonzero on −ξ < x < 1. Working with this distribution, one obtains the following tensor form factors as the leading Mellin moments: Evaluated using the CI parameters in Table 1, These quantities are subject to QCD evolution; and, as described after Eq. (29), we interpret the results in Eq. (60) as being valid at the hadronic scale, the value of which is discussed in Refs. [21,22]: Using QCD's infrared-finite process-independent effective charge [55],α(k 2 ), to integrate the evolution equations [21,22], one finds Consequently, at ζ = ζ 2 = 2 GeV, This is the scale used in Ref. [100], which reports the following values for these quantities after an extrapolation to the physical pion mass: 0.22 (3), 0.039(10), 5.66 (60), in qualitative agreement with the CI results. Similar conclusions are drawn elsewhere, e.g. Refs. [101,113,114]. The tensor form factors in Eqs. (59) are plotted in Fig. 6, normalised by their 2 = 0 values. Employing this procedure, the depicted form factors are independent of the renormalisation scale [101]. Hence, comparison with the lQCD results in Ref. [100] is meaningful, although quantitative agreement cannot be expected because the lQCD form factors were computed using m 2 π ≈ 20 m 2 empirical π . Bearing this in mind and considering that the CI produces hard pseudoscalar meson form factors, there is reasonable qualitative agreement, e.g.: the radii have the same ordering, r B π 10 /r B π 20 = 1.48(17) (lQCD) vs. 1.14 (herein); and B 10 (t) is generally softer than B 20 (t).
One now has access to the light-front transverse-spin distribution of dressed-quarks within the pion, which is defined in impact-parameter space [100]: with where q π (x, |b ⊥ |) is given in Eq. (50) and J 1 is a Bessel function. For a dressed-quark polarised in the +x direction andŝ ⊥ ·b ⊥ = cos φ ⊥ , ε i j s i ⊥ b j ⊥ = |b ⊥ | sin φ ⊥ . As emphasised above, in an internally consistent CI treatment, all pion form factors are hard; so the integrals that define the transverse densities in Eq. (65) are ill defined. It is nevertheless worth illustrating the character of ρ 1 (b ⊥ , s ⊥ ). We therefore employ the expedient introduced in Eq. (57), choosing the mass-scale "M" to reproduce the CI result for the 2 0 slope of a monopole approximation to the given form factor and setting its 2 = 0 value to match the CI value; to wit, with M F = 1.09 GeV, M B = 1.02 GeV. The result is drawn in Fig. 7. Figure 7 shows that for a dressed valence-quark polarised in the light-front-transverse +x direction, the transverse-spin density is no longer symmetric around b ⊥ = (b x = 0, b y = 0). Instead, the peak is shifted to (b x = 0, b y > 0), with strength transferred from b y < 0 to b y > 0. The average transverse shift is [100]: and the b y profile remains symmetric around the line b x = 0. We interpret these results as being valid at ζ H . The distortion vanishes logarithmically with B π 10 (0) → 0 under QCD evolution, Eq. (62).
Given that E Td π + (x, ξ, t) = −E T u π + (−x, ξ, t), then the three-dimensional profile for a s ⊥ x dressed valenceantiquark is obtained by rotating Fig. 7

Twist-two TMDs
Recall now that one proceeds from a given GTMD to the associated TMD by setting = 0, which also means ξ = 0. At twist-two, our CI treatment (which does not include a Wilson line) produces one nonzero TMD, whose form can be read from Eq. (31) (ς := σ k 2 ⊥ ,0 1 ): This TMD, describing the dressed valence u-quark in the π + , is depicted in Fig. 8. (Note that M 2 f 1 (x, k 2 ⊥ ) is dimensionless.) The root-mean-square value of k 2 ⊥ is defined via Evidently and unsurprisingly, the symmetry-preserving CItreatment produces a hard k 2 ⊥ distribution even at the hadronic sale, ζ H . In contrast, a pion twist-two TMD developed from an interaction with QCD-like momentum dependence yields [115] k 2 ⊥ 1/2 = 0.21 GeV. For additional comparison, we note that a Nambu-Jona-Lasinio (NJL) model has also been used to compute f 1 (x, k 2 ⊥ ) [116]. NJL model results are sensitive to the regularisation scheme employed. Ref. [116] used a Pauli-Villars procedure and obtained results that agree semi-quantitatively with those depicted in Fig. 8.
Owing to gluon radiation and additional fragmentation, the distributions in Fig. 8 which is the π + valence u-quark distribution function.
Since we omit the Wilson line, our result for the pion's Boer-Mulders function is

Twist-three TMDs
In the absence of a Wilson line, the CI supports two nonzero twist-three pion TMDs. The first is obtained from the GTMD E 2 (x, k 2 ⊥ , ξ, t) in "Appendix B.1": This TMD is chiral-odd, viz. it is associated with an interaction-induced quark chirality flip within the target. e(x, k 2 ⊥ ) vanishes in the chiral limit, m π = 0. The upper panel of Fig. 9 depicts the CI result forê(x, k 2 ⊥ ) at the hadronic scale, ζ H . The lower panel highlights the xdependence of its k 2 ⊥ profile: i.e. the |k ⊥ | width ranges from 0.39 GeV at x = 0 to 0.22 GeV at x = 1. Given the hardness of CI form factors, it is most appropriate to make an internally consistent comparison; hence, we observe that Eq. (74) means the width of e(x, k 2 ⊥ ) ranges from 63% → 36% of the width of the chiral-even TMD f 1 (x, k 2 ⊥ ), with mean value 51%. Comparing the images in Fig. 9 with those in Fig. 8, one sees thatê(x, k 2 ⊥ ) is at most two-thirds the size of f 1 (x, k 2 ⊥ ) and typically smaller. In any cross-section, this suppression is compounded by the higher-twist factor m π /n · P.
The second twist-three TMD, which is chiral-even, may be read from "Appendix B.3": (75b) Fig. 10 -upper panel; and the lower panel illustrates the x-dependence of its k 2 ⊥ profile: The |k ⊥ | width varies from 0.32 GeV at x = 0 to 0 at x = 1, owing to the (1 − x) factor in Eq. (75b), i.e. the width of f ⊥ (x, k 2 ⊥ ) ranges from 52% → 0% of the width of the chiraleven TMD f 1 (x, k 2 ⊥ ), with mean value 37%. Comparison of the images in Fig. 10 with those in Fig. 8 reveals thatf ⊥ (x, k 2 ⊥ ) is not more than two-thirds the size of f 1 (x, k 2 ⊥ ) and almost always much smaller. In any crosssection, this suppression is compounded by the higher-twist factor (m π /M)(m π /n · P).

Twist-four TMD
The CI supports a single twist-four pion TMD, which is chiral-even and can be read from "Appendix C.1": f 3 (x, k 2 ⊥ ) is nonzero in the chiral limit so long as the full CI pion Bethe-Salpeter amplitude is used, i.e. F π = 0 in Eq. (16).
We depictf 3 (x, k 2 ⊥ ) in Fig. 11 -upper panel; and in the lower panel sketch the x-dependence of its k 2 ⊥ profile: Here the |k ⊥ | width ranges from 0.34 GeV at x = 0 to 0.31 at x = 1, i.e. the momentum-space breadth of f 3 (x, k 2 ⊥ ) ranges from 56% → 51% of the width of f 1 (x, k 2 ⊥ ), with mean value 53%.
Comparing the images in Fig. 11 with those in Fig. 8, it is plain thatf 3 (x, k 2 ⊥ ) is typically less than one-third the size of f 1 (x, k 2 ⊥ ). This suppression multiplies that introduced into cross-sections by the higher-twist factor (m π /n · P) 2 .
It is also appropriate to remark that whilst the suppression of twist-3 and -4 TMDs relative to the twist-2 TMD is a feature our CI analysis, it is not found in all analyses. The models studied in Ref. [117] highlight this point.

Wigner distribution
Given that (a) GPDs and TMDs can both be obtained directly from Wigner distributions and (b) a given Wigner distribution is obtained by computing a Fourier transform of the associated GTMD at ξ = 0, it is worth presenting a concrete result for the simplest of the Wigner distributions for a dressed-quark in the pion. To this end, recall Eq. (31) and consider Inserting the explicit form of the integrand, one finds where This function has nonzero support on x ∈ [0, 1]. The dimensionless Wigner function in Eq. (80) is plotted in Fig. 12. Each panel shows a different value of |b ⊥ |, viz. 0.1 fm and 0.2 fm. This valence-quark Wigner function is (i) sharply peaked at (x = 1, k 2 ⊥ = 0, b 2 ⊥ = 0); (ii) exhibits power-law suppression as k 2 ⊥ and/or b 2 ⊥ are increased; and (iii) is negative on a neighbourhood (x 1, k 2 ⊥ 0). We anticipate that the analogous Wigner function computed with a realistic interaction will display similar behaviour.

Summary and perspective
We used a vector × vector contact interaction (CI), treated at leading-order in a widely-used symmetry-preserving Dyson-Schwinger equation (DSE) truncation scheme, to calculate an array of twist-two, -three and -four pion GTMDs (Sect. 4, Appendices B, C). Whilst some of the results are particular to the CI, many features are physically relevant, including an observation that the strength and shape of all pion GTMDs are largely set by the scale of emergent hadronic mass (EHM) in the strong interaction. In a few particular cases for which CI limitations were too conspicuous, we augmented the analysis by appealing to continuum-and lattice-QCD results in order to arrive at realistic illustrations of material points.
Concerning GPDs, we found (Sect. 5.2) that the pion's θ 2 mass distribution form factor is harder than its electromagnetic form factor, F em π ; and in turn, F π is harder than the pion's θ 1 gravitational pressure distribution form factor. Concerning the pressure distribution, the peak value, lying in the neighbourhood of the pion's core, is approximately five-times greater than that in the proton; indeed, it is commensurate with the pressure at the core of a neutron star. Moreover, the shear pressure achieves its maximum value when the confinement pressure comes to exceed that generated by the forces driving the quark and antiquark away from the core. (The Ansatz for θ 1 used in the pressure calculations capitalises on lattice-QCD results for this form factor.) The tensor GPD provides information about transversity in the pion; and we found (Sect. 5.3) that polarising a pion's dressed quark in the positive-x direction of the light-fronttransverse plane produces a clear distortion of the transversespin density, shifting its peak in the positive−y direction. This distortion diminishes as the resolving scale is increased.
The pion's GTMDs also provide direct access to its transverse momentum dependent distribution functions (TMDs); and in the absence of a model for the Wilson line, the CI supports four that are nonzero (Sect. 6): one of twist-two, two twist-three, and one twist-four. Our calculations indicate that the twist-two TMD, f 1 (x, k 2 ⊥ ) is largest in magnitude and possesses the greatest domain of material k 2 ⊥ -support. The twist-three distributions, e and f ⊥ , are uniformly smaller; and the twist-four TMD, f 3 , is still smaller. In any crosssection, such suppressions would be compounded by the respective m π /n· P and (m π /n· P) 2 twist-expansion factors.
Wigner distributions are a natural complement to GTMDs, providing an intuitive visual aid to expressing and understanding their physical content. We therefore provided results for a representative example, viz. that associated with the twist-two GTMD that produces the pion's valence-quark distribution function, and electromagnetic and gravitational form factors (Sect. 7). At the hadronic scale, this Wigner function is sharply peaked in the neighbourhood of (x = 1, k 2 ⊥ = 0, b 2 ⊥ = 0) and broadens as the transverse position variable conjugate to the probing momentum transfer, b ⊥ , increases in magnitude. Similar behaviour should be expected of such Wigner distributions calculated with a realistic interaction.
Several extensions of the work described herein immediately suggest themselves. (A) Kindred analyses for the kaon, which would reveal physical effects on GTMDs that arise from constructive interference between Nature's two mass generating mechanisms: EHM and Higgs-boson induced. (B) Development of a practicable realisation of the Wilson line, because it would enable computation of time-reversal-odd GTMDs, whose comparison with the timereversal-even functions calculated herein may yield additional insights that could be exploited in studies using realistic interactions. (C) Repeating this analysis using realistic light-front wave functions for the pion (and kaon), whose profiles are known to explain and predict a diverse array of pseudoscalar meson properties. All these efforts are underway. BK20180323); and Jiangsu Province Hundred Talents Plan for Professionals.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All data generated during this study are represented in this published article.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Useful formulae
Eq. (14) is the first of many integrals appearing herein whose regularised values are expressed in terms of incomplete gamma-functions. In general (n ∈ Z, n ≥ 0): They can usefully be illustrated with simple examples: In general, Such expressions are useful, e.g. in expressing the Bethe-Salpeter kernel in Eq. (20): ω(α, Q 2 )) . (A.4d) We recall here that Eq. (20) is an eigenvalue problem with a solution for Q 2 = −m 2 π , at which point the eigenvector is the pion's Bethe-Salpeter amplitude. When computing observables, one must employ the canonically normalised amplitude, viz. π rescaled such that (A.6) In the chiral limit, viz. using solutions obtained with m = 0 in Eq. (10), Eqs. (A.5), (A.6) impose [33]: (A.7) The function ω(α, Q 2 ) is defined in Eq. (19). Similar arguments appear in the expressions for various pion GTMDs. We list them here.
When describing TMDs, we also use where