Pseudo-scalar mesons: light front wave functions, GPDs and PDFs

We develop a unified algebraic model which satisfactorily describes the internal structure of pion and kaon as well as heavy quarkonia ($\eta_c$ and $\eta_b$). For each of these mesons, we compute their generalized parton distributions (GPDs), built through the overlap representation of their light-front wave function, tightly constrained by the modern and precise knowledge of their quark distribution amplitudes. From this three-dimensional knowledge of mesons, we deduce parton distribution functions (PDFs) as well as electromagnetic form factors and construct the impact parameter space GPDs. The PDFs for mesons formed with light quarks are then evolved from the hadronic scale of around 0.3 GeV to 5.2 GeV, probed in experiments. We make explicit comparisons with experimental results available and with earlier theoretical predictions.


I. INTRODUCTION
The generalized parton distributions (GPDs) have emerged as a comprehensive tool to describe hadron structure probed in hard scattering processes [1][2][3][4][5][6][7][8]. GPDs connect hadron electromagnetic form factors (FFs) [9], measured in elastic processes, to longitudinal parton distributions (PDFs) which are probed in deep inelastic scattering [10]. They provide a kaleidoscopic view of the three-dimensional spatial structure of hadrons, written as a function of the longitudinal momentum fraction x, the momentum transfer t and the skewness variable ξ (longitudinal momentum transfer); furthermore, Fourier transform of the GPDs yields the transverse spatial distribution of partons correlated with x [11], the so called impact parameter space GPDs (IPS-GPDs).
Most physical observables related to mesons can be calculated through a combined knowledge of their Bethe-Salpeter amplitude (BSA) and the quark propagator [12,13]. While in principle it can accurately be achieved through a cumbersome computation of the quark propagator Schwinger-Dyson equation (SDE) and the Bethe-Salpeter equation (BSE) in close connection with full QCD [14], calculation of a plethora of experimentally interesting quantities such as FFs [15][16][17][18][19][20][21][22], distribution amplitudes (PDAs), PDFs [23][24][25][26][27][28][29], and, specially, GPDs [30][31][32][33][34], remains a highly non-trivial pursuit. However, our understanding of the intricate interplay between the quark propagator and the meson BSA permits us to construct their simple Anstäze sufficiently efficacious to make reliable predictions and amicable enough to offer algebraic manipulations. In this article, we carry out this algebraic model (AM) construction for pseudo-scalar mesons in terms of a form-invariant spectral density. The most attractive feature of this AM is that the spectral density is explicitly written in terms of the leading-twist PDA whose existing reliable information allows us to circumvent the need to construct any ad hoc Ansatz for the spectral density.
We begin with an evidence-based Ansatz for the quark propagator and the BSA in terms of a spectral density function which is form-invariant for all the ground-state pseudo-scalar mesons. The Bethe-Salpeter wave function (BSWF) can then be readily constructed whose subsequent projection on to the light front yields the highly sought-after light-front wave function (LFWF). Its integration over the transverse momentum squared (k 2 ⊥ ) gives us access to the valence-quark PDA. We exploit our current detailed and accurate knowledge of the PDAs of pseudo-scalar mesons [27,35] to determine the parameters of our model. We use the overlap representation of the LFWF [3] to compute the GPDs of pion, kaon, η c and η b . From this three dimensional knowledge of these mesons, different limits/projections lead us to deduce the PDFs, FFs and the IPS-GPDs which are then compared to available experimental extractions of these observables.
This article has been organized as follows. In Sec. II, we present our generalized AM for the quark propagator and the BSA of the pseudo-scalar mesons under consideration in terms of a spectral density function. It allows us to derive the leading-twist LFWF in Sec. III by merely appealing to the definition of its Mellin moments. The resulting LFWF permits establishing a closed and simple algebraic connection with the PDA, so that the need to specify a spectral density is completely avoided if the PDA is known. Sec. IV details the extraction of GPDs through the overlap representation of the LFWFs as suggested in [3], in the so called DGLAP kinematic region, and a series of distributions derived therefrom: PDFs, FFs and IPS-GPDs. In Sec. V, we particularize our AM arXiv:2207.06550v1 [hep-ph] 13 Jul 2022 to produce a collection of distributions, using inputs from previous SDE predictions, and compare (when possible) with available theoretical calculations and experimental data. Finally, in Sec. VI, we present a summary of our work and the scope of our model.

II. ALGEBRAIC MODEL
The internal dynamics of a meson can be described, in a fully quantum field theoretic formalism, via its BSWF. In terms of the associated BSA (Γ M ), and the quark (antiquark) propagators (S q,(h) ), the BSWF reads where k − = k − P/2 and P 2 = −m 2 M is the (negative) mass squared of the meson M. The labels q andh which denote the valence quark and antiquark flavors are in general different but might also be the same. Although the propagators and BSA might be obtained from solutions of the corresponding SDEs and BSEs, useful and relevant insight can be extracted from sensibly constructed simpler models.
Plain expressions for the quark (antiquark) propagator and BSAs that capture QCD's key non-perturbative traits are given by where ∆(s, t) = (s + t) −1 ,∆(s, t) = t∆(s, t), k w = k + (w/2)P . Herein, M q(h) is a mass scale akin to a constituent mass for a given quark flavor and n M is a normalization constant that will be determined later. The function ρ M (w) can be regarded as a spectral density, whose particular form determines the pointwise behavior of the BSA, therefore having a crucial impact on the meson observables. The parameter ν > −1 controls the asymptotic behavior of the BSA; this is discussed in detail below. Finally, Λ 2 w ≡ Λ 2 (w) is defined as follows: Notice that, unlike kindred models [30][31][32][33][34][36][37][38][39][40] which have been employed successfully to compute an array of GPD-related distributions, we have promoted Λ → Λ w to incorporate a w-dependence. Keeping in mind the efficacy of earlier models, we point out some key differences which yield a simplification of relevant integrals and closed algebraic expressions relating different distributions: • We retain the constant term from the original models, setting it to M q .
• There is a term linear in w which is the only term not symmetric under w ↔ −w. This asymmetry allows us to study mesons with different flavored quarks and is hence accompanied with the multiplicative factor of (M h 2 − M 2 q ). When quarkantiquark are of the same flavor, this term ceases to contribute by construction.
• Then we have a quadratic term in w 2 which has the coefficient proportional to m 2 M . We choose the coefficients of each power of w to ensure that the condition: guarantees the positivity of Λ 2 (w). When quark and antiquark have the same flavor, the left part of the inequality is trivially satisfied. We consider isospin symmetry, i.e., M u = M d . In any other case, like that of a kaon or heavy-light mesons, the ratio Mh/M q must be set with care. Realistic solutions of the quark SDE provide useful benchmarks [41]. Note that m M < Mh +M q is satisfied for Nambu-Goldstone bosons. One can find sensible values for the constituent masses to uphold this inequality for the ground state pseudo-scalar mesons.
Combining Eqs. (1)-(3), the BSWF acquires the following Nakanishi integral representation (NIR): where the profile function,ρ ν M (w), has been defined in terms of the spectral density as The function M q,h (k = p+P, P ) has the following tensor structure that also characterizes χ M (k, P ): Due to the trace over Dirac indices, cf. eq. (13), the last two terms containing an even number of γ-matrices in the above eq. (8) do not contribute to the leading-twist light-front wave function (LFWF) and, consequently, the PDA. The function D ν q,h (k, P ) is a product of quadratic denominators, Feynman parametrization enables us to combine the denominators in Eq. (9) into a single one. Then a suitable change of variables and a subsequent rearrangement in the order of integration yields the expression: F M (α, σ ν+2 ) = ν(ν + 1)  [27,35], and parameterized according to Eqs. (54). For comparison, the asymptotic distribution, φasy(x) = 6x(1 − x), is also shown.
where σ = (k − αP ) 2 + Λ 2 1−2α , and α, β are Feynman parameters. Since only (1−β) ν−1 depends on β, integration over dβ can be performed directly, thus yielding F M (α, σ ν+2 ) = 2 ν (ν + 1) As we explain in the Appendix, this extra algebraic integration allows us to completely deriveρ ν M (w) in terms of the PDA. In the next section we shall explicitly see that, when employing this model for the BSWF, many quantities and relations of interest can be obtained in a purely analytical manner.

III. LIGHT FRONT WAVE FUNCTIONS AND PARTON DISTRIBUTION AMPLITUDES
For a quark q within a pseudo-scalar meson M, the leading twist (2 particle) light-front wave function, ψ q M , can be obtained via the light-front projection of the meson's BSA as: where δ x n (k M ) = δ(n · k − x n · P ); n is a light-like fourvector, such that n 2 = 0 and n · P = −m M ; a mentioned before, x corresponds to the light-front momentum fraction carried by the quark. The trace is taken over color and Dirac indices. The notation dk ≡ d 2 k π has been employed and the 4-momentum integral is defined as usual: The moments of the distribution are: From Eqs. (10)- (15), one arrives at where σ ⊥ = k 2 ⊥ + Λ 2 1−2α . Uniqueness of the Mellin moments, Eqs. (15)- (16), implies the connection between the Feynman parameter α and the momentum fraction x; therefore one can identify the LFWF as Notice that the above expression resembles the one derived, for instance, in [30,31,40]. However, the crucial difference is the w-dependent definition of Λ w , Eq. (4). As mentioned before, its particular form enables additional simplicity and allows amicable algebraic manipulation as will be evident shortly.
where f M is the leptonic decay constant of the meson. From Eqs. (12) and (17), it is seen that the only term in the above equation that depends on k ⊥ is 1/σ ν+1 ⊥ , then Combining Eqs. (17)- (19) we arrive at the following algebraic relation between ψ M (x, k 2 ⊥ ) and φ M (x): The compact result above is a merit of the AM we have put forward. Throughout this manuscript, we shall employ dimensionless and unit normalized PDAs, 1 0 dxφ q M (x) = 1. The resulting PDA and LFWF are expressed in a quasiparticle basis at an intrinsic scale, intuitively identified with some hadronic scale, ζ H , for which the valence degrees of freedom fully express the properties of the hadron under study. Most results herein are quoted at ζ H (unless specified otherwise). However, for the sake of simplicity, the label ζ H shall be omitted. It is worth reminding that the quark and antiquark PDA are connected via momentum conservation, a constricted and firm connection that prevails even after evolution [42][43][44]. Some practical corollaries of the AM and Eq. (20): • Given a particular form of φ q M (x), the ψ q M (x, k 2 ⊥ ) can be obtained quite straightforwardly.
• As long as we have reliable access to φ q M (x), there is no actual need to construct the profile functioñ ρ ν (w) (although it can be properly identified, as we explain in the Appendix).
• It also works the other way around. A sensible choice ofρ ν M (w) and model parameters yields algebraic expressions for both φ q M (x) and ψ q M (x, k 2 ⊥ ).
• In fact, the present AM can be reduced to the toy model employed in Refs. [37][38][39] with appropriate substitutions. It also faithfully reproduces the results obtained from the more sophisticated Ansatz in Ref. [30][31][32].
• The degree of factorizability of the LFWF is clearly exposed through Eqs. (4) and (20).
Regarding the last point let us consider the chiral limit The bracketed term no longer depends on x; hence, the x and k ⊥ dependence of ψ M (x, k 2 ⊥ ) has been completely factorized. Conversely, as captured by Eq. (20), a nonzero meson mass and quark/antiquark flavor asymmetry, namely m 2 M = 0 and (M 2 h − M 2 q ) = 0, yield a LFWF which correlates x and k 2 ⊥ . So one should expect an increasingly dominant role of x and k 2 ⊥ correlations in heavy-quarkonia and heavy-light systems. Notably, a soft Q 2 -dependence might also be introduced in the definition of the PDA [45,46], Eq. (18), producing the following compact expression: Clearly, φ(x; Q 2 → ∞) = φ(x), which is the limit we take for the sake of the discussion. In the next section, we shall exploit the virtues of Eq. (20) to compute the pseudo-scalar meson GPDs in the overlap representation.

IV. GENERALIZED PARTON DISTRIBUTIONS
The valence quark GPD can be obtained from the overlap representation of the LFWF [5], namely: If p (p ) denotes the initial (final) meson momentum, then P = (p+p )/2 and −t = ∆ 2 = (p−p ) 2 (the latter defines the momentum transfer); . Both x and ξ have support on [−1, 1], but the overlap representation is only valid in the DGLAP region, |x| > ξ. The kinematical completion (the extension to the ERBL domain), required to fulfill the polinomiality property [5], can be achieved through the covariant extension from Refs. [33,34,36,37]. Notwithstanding, the GPD is even in ξ and only nonzero for the valence quark if x > −ξ (the antiquark GPD is non-zero if x < ξ); hence, in the following, we shall restrain ourselves to ξ ≥ 0. Notice again that Eq. (24) implies that the meson is described as a two-body Fock state. This picture is then valid at the hadronic scale, in which the fully dressed quark/antiquark quasiparticles encode all the properties of the meson.
We now work out the expression for the valence quark GPD in detail by substituting Eq. (20) in Eq. (24) As usual, integration on k ⊥ can be performed by introducing Feynman parametrization and a suitable change of variables, such that the integral in Eq. (25) becomes where the function M 2 (u) depends on the model parameters, as well as the kinematic variables x, ξ, t. It acquires the form M 2 (u) = c 2 u 2 + c 1 u + c 0 , where Thus the GPD can be conveniently expressed as Notice that, in the chiral limit, M 2 (u) reduces to and so the integration on du in Eq. (28) can be carried out algebraically for specific values of ν > − 1. In particular, ν = 1 recovers the results in [33,34,[37][38][39]. Beyond the chiral limit, an algebraic expression is found for t = 0: where pFq (u, v, w, z) is the regularized hypergeometric function. Conversely, taking ξ = 0, an expansion of M 2 (u) around −t ≈ 0 yields an algebraic solution for Eq. (28): In the next section we will focus on the forward limit of the GPD (t = 0, ξ = 0) which defines the valence quark PDF. For the time being, we can make an insightful connection with light-front holographic QCD (LFHQCD) approach Ref. [47,48], recalling the following representation for the zero-skewness valence quark GPD therein: wheref q M is some profile function to be determined. An expansion around −t ≈ 0 of this expression, and a subsequent comparison with Eq. (31), enable us to identifŷ It is also useful in extracting insights concerning the IPS-GPDs, as will be addressed below.
We now proceed to discuss the derivation of PDFs, FFs and IPS-GPDs, as inferred from the knowledge of the GPDs in the DGLAP kinematic region.

A. Parton distribution functions
The first term of the Taylor expansion in Eq. (31) corresponds to the valence quark PDF, namely where q M (x) is unit normalized. Recalling that the distributions have been derived at ζ H , the corresponding antiquark PDF is simply obtained as Furthermore, the factorization properties of the LFWF in the chiral limit yields the simple relation: thus stressing that the degree of factorizability of the AM is manifest via the quantity Λ 2 1−2x . As long as we have m 2 M ≈ 0 and also (M 2 h − M 2 q ) ≈ 0, a factorized LFWF will produce sensible results. This is the case of the SDE results from Refs. [26,27], in which Eq. (36) was employed to compute the kaon PDF from its PDA. For the purpose of this work, factorizability will not be assumed and we shall consider the more general case, Eq. (34). The set of relations described in this Section  [49], is rescaled according to the ASV analysis in [50]. Lower panel-The dotted (cyan) line corresponds to u-in-K valence-quark PDF, the dot-dashed (blue) line is the analogous for thes quark and the solid (purple) line corresponds to the u valencequark in the pion. The error bands account for the variation of the initial scale, ζH = 0.33 (1 ± 0.1) GeV. also shows that if the input PDA behaves like φ(x → 1) ∼ (1−x) (as prescribed by QCD, [42]), the PDF will exhibit the large-x behavior q M (x; ζ H ) ∼ (1 − x) 2 . Finally, it is worth recalling that neither Eq. (36) nor Eq. (35) remain valid for ζ>ζ H , due to the evolution equations obeyed by the PDFs [51][52][53][54].
All distributions described so far have been obtained from the LFWF at the hadron scale, ζ H ; as described before, at this low-energy scale, the fully dressed quasiparticles (valence-quarks) express all hadron properties. This is also the case of the valence-quark PDF which, computed at ζ H , entails that all the hadron's momentum is carried by the fully-dressed valence quarks. From the experimental point of view, the access and interpretation of PDFs and GPDs at ζ H imply certain technical and conceptual complications [10]; only above certain en-ergies, typically the mass of the proton, parton distributions can be properly extracted. In particular, experimental data for the case of the pion is only available at ζ = ζ 5 := 5.2 GeV [50,55] (the same for the u K (x)/u π (x) ratio [56]), whereas ζ = ζ 2 := 2 GeV is a typical scale for lattice QCD and phenomenological fits [57][58][59]. To produce a consistent picture when evolving the hadronic scale PDF, we shall follow the all orders scheme introduced in Refs. [24][25][26][27]60] for pion and kaon PDFs, extended to their GPDs in Ref. [30,31], and employed recently in the calculation of the proton PDFs as well [61]. This scheme is based upon the assumption that an effective chargeα allows all beyond leading-order effects to be absorbed within it, thus arriving at a leading-orderlike DGLAP evolution equation. Notably, if the evolution is performed via the computation of several Mellin moments, it is not necessary to specify the pointwise behavior of the effective charge [31] (assuming its existence would be sufficient). To evolve the distributions directly, the exercise we carry out in this article, we takeα from Ref. [27], which implies setting ζ H = 0.33(1 ± 0.1) GeV. In Section V, we present numerical results for evolved pion and kaon PDFs for specific model inputs described therein.

B. Electromagnetic form factor
The contribution of the q quark to the meson's elastic electromagnetic form factor (EFF) is obtained from the zeroth moment of the GPD: an analogous expression holds for the antiquarkh, such that the complete meson EFF reads where e q,h are the valence-constituent quarks electric charges in units of the positron charge. Due to polynomiality properties of the GPD, the EFF does not depend on ξ, therefore one can simply take ξ → 0: A Taylor expansion around t ≈ 0 yields where r q M denotes the contribution of the quark q to the meson charge radius, r M . Comparing the above equations with the integration of Eq. (31) on x, one obtains a semi-analytical expression for r q M : showing the charge radius is tightly connected with the hadronic scale PDF (and thus with the corresponding PDA). The antiquark result is obtained analogously. This contribution to r M reads: wherefh M (x) is defined in analogy to its quark counterpart in Eq. (33),fh Summing up the quark and antiquark contributions, the meson charge radius reads: Clearly, in the isospin symmetric limit, e q + eh = 1 yields F M (t) = F q M (t) and so r M = r q M . For the neutral pseudoscalars, in the isospin symmetric limit, e q +eh = 0 implies F M would be strictly zero, producing r M = 0; thereby we focus only on the individual flavor contribution (F M → F q M ) in such cases (e.g. heavy quarkonia). Finally, note that if the charge radius is known, then Eqs. (42)(43)(44)(45) can be employed to fix the model parameters.

C. Impact parameter space GPD
The IPS-GPD can be obtained straightforwardly by carrying out the Fourier transform of the zero-skewness GPD, H q M (x, 0, t): where J 0 is the cylindrical Bessel function. This distribution is interpreted as the probability density of finding a parton with momentum fraction x at a transverse distance b ⊥ from the centre of transverse momentum of the meson under study. It is extracted in its totality by the GPD's properties in the DGLAP region. Exploiting the representation of the GPD from Eq. (32), we can obtain an analytic expression: Containing an explicit dependence on the PDF, Eq. (47) reveals a clear interrelation between the momentum and spatial distributions. In fact, the PDF is recovered from Furthermore, considering the mean-squared transverse extent (MSTE),  Table 1. Dashed (black) line is the SDE result for the pion [15]. Right panel-It shows analogous results for the kaon FF. Diamonds, rectangles and circles represent the experimental data from Refs. [62][63][64]. Lower (gray) band is the SDE result for the kaon [65].
FIG. 8. ηc and η b electromagnetic FFs. Left panel-The purple band represents our ηc results with the model parameters described in Section V. The band width accounts for a 5% variation of the benchmark charge radius in Table 1. Right panel-Analogous results for η b . For comparison, we have included lattice QCD results from Refs. [66,67], as well as SDE-driven predictions in the contact interaction (CI) model and a former algebraic model for heavy quarkonia [68,69].
the IPS-GPD defined in Eq. (47) yields the plain relation: Integrating over x, and comparing with Eq. (45), one is left with a compact expression for the expectation value: i.e. the expectation value of the MSTE of the valence quark is directly correlated with the meson charge radius.
In the isospin symmetric limit, the following expected result [30,31] is recovered: Interestingly, in the chiral limit, all the algebraic expressions form this Section, valid only at ζ H , become plainly analogous to those from the factorized Gaussian model in [30,31].
In the following section we shall provide a collection of results for the distributions discussed so far, using SDE predictions as model inputs.

V. COMPUTED DISTRIBUTIONS
Now that we have shown a variety of algebraic results for different distributions of partons (and some other quantities), we will particularize the inputs of the AM. The starting point is Eq. (20), which directly relates the leading-twist LFWF with the PDA such that, with the prior knowledge of φ q M (x), the LFWF is derived straightforwardly; the produced physical picture would be valid at ζ H . Given the robustness of the SDE formalism to compute PDAs, we shall employ predictions obtained within this framework as model inputs [27,35]. The specific set of PDAs we consider is the following (x = 1 − x): The expressions above properly capture our contemporary knowledge of such distributions, namely, the soft endpoint behavior and the dilation/compression with respect to the asymptotic distribution [42]: As can be seen in Fig. 1, pion and kaon PDAs are dilated with respect to φ asy (x), while those containing heavy quarks are narrower. As noted for the kaon, the asymmetry between the s and u-quark masses produces a skewed distribution, while the rest of the PDAs are symmetrical. The remaining ingredients are the parameter ν and the constituent masses M q . Regarding the former, ν = 1 is a natural choice since it yields the correct asymptotic behavior of the BSWF [12]. Concerning the values of the constituent masses, we shall employ available experimental [70], SDE [15,[20][21][22]71] and lattice QCD [66,67] results on the charge radii as benchmarks, and determine M q via Eq. (45). Table I collects the constituent quark masses that define our AM and the corresponding charge radii.   With the AM fully determined, the produced LFWFs are shown in Fig. 2. It is clear that the heavier mesons exhibit a much slower damping as k 2 ⊥ increases. Furthermore, just as the PDAs, the LFWFs as a function of x are found to be more compressed in this case.
The valence-quark GPDs are then obtained appealing to the overlap representation of the LFWF, Eqs. (24,28). Pion and kaon results are shown in the bottom panel of Fig. 3, while those of η c and η b can be found in Fig. 4. The GPDs for the heavier mesons naturally have a narrower profile along the x-axis and are harder along the −t-axis. Moreover, the upper panel of Fig. 3 also displays a comparison between the GPDs obtained directly from Eq. (28) and the approximate representation of (32). The derived valence quark PDFs are found in Fig. 5. As one would expect from Eq. (34), the characteristic features exhibited by the PDAs, of dilation and narrowness, are filtered into PDFs. To emphasize it, we notice that the plots in the above mentioned figure display the scale-free parton-like profile Given our preferred value M s ≈ 1.8 M u , the s-in-K momentum fraction at the hadronic scale is <x; ζ H > s K = 0.55, about 4% larger than typical values [26,27]. The pion and kaon PDFs are then evolved from the hadronic scale, ζ H = 0.33(1 ± 0.1) GeV, to the experimentally accessible scale of ζ 5 := 5.2 GeV. The evolution procedure is detailed, for instance, in Refs. [31,60]. Fig. 6 displays the outcome. In the top panel of this figure, the valence quark as well as gluon and sea quark pion PDFs are shown. At the evolved scale, we find typical values of momentum fraction distribution in pion [26,27]: <x; ζ H > val π = 0.41(4), <x; ζ H > sea π = 0.14(2), <x; ζ H > glue π = 0.45 (3). The bottom panel of Fig. 6 compares the valence quark PDFs in pion and kaon. Then again, our choice of M s produces a slightly larger momentum fraction for the s valence-quark at such scale, <x; ζ 5 > s K = 0.25, and a smaller one for the u quark, <x; ζ 5 > u K = 0.17. Concerning the large-x exponents of the valence quark distributions, we find that where β ef f must be interpreted as an effective exponent rather than that obtained from the known evolution equations of β(ζ H ) [27,72]. Moreover, the x-domain of applicability and interpretation of β(ζ H ) is not without its ambiguities and requires special care [73]. The electromagnetic FFs are displayed in Figs. (7,8). As can be noted therein, pion and kaon FFs agree with the available experimental data [62][63][64] and previous SDE calculations [15,65]. The η c FF is compared with lattice QCD [66,67] and SDE results in the contact interaction (CI) model [68]. Similarly, the η b result is contrasted with CI model results and with previous determinations with an AM for heavy quarkonia [69]. Both η c and η b form factors show a satisfactory compatibility with earlier reliable predictions. The IPS-GPDs are derived from the approximate LFHQCD-inspired parametrization of the GPD, introduced in [48] and quoted in Eq. (32). For illustrative purposes, we have considered the convenient representation of Eq. (50), which produces the pion and kaon results shown in Fig. 9. The quark region is identified with x>0, while the antiquark lies in the x<0 domain. The symmetry in the pion case is a natural consequence of the isospin symmetry, whereas the contraction on the s-in-K distribution is a result of M s being larger than M u . In fact, as the constituent quark mass becomes larger, it is expected that the quark plays an increasingly major role in determining the center of transverse momentum; furthermore, the distributions become narrower and the maximums become larger. Given the compact representation of the IPS-GPDs, the values (x max , b max ⊥ ) q M where b q M (x, b ⊥ ) acquires its global maximum, can be readily identified: and x max is the real-valued solution of It is thus clear that a constant PDF yields the point particle limit (|x max |, b max ⊥ ) q M → (1, 0). The location of the maximum and its value are reported in Table II for different mesons. Finally, according to Eq. (52) and our model inputs, we report the expectation values of the MSTE for the kaon: while for the heavy quarkonia and pion in isospin symmetric case we can infer the result from Eq. (53).  For the pion and kaon, one can visually verify these tabulated results in Fig. (9). This completes the presentation of computed results.

VI. SUMMARY AND SCOPE
In this article, we put forward a fairly general AM for the pseudoscalar meson BSWF, which preserves its primary attractive feature of guaranteeing most calculations continue to be analytic. For systematic and visual clarity, we italicize its main features and our key results as follows: • The key functions of the model are the spectral density ρ M (w) and Λ(w), which play the defining role for the dominant BSA, Eq. (3), of the pseudoscalar mesons we study.
• The function Λ(w), defined through Eq. (4), is quadratic in w which is as high as we can go in the power of this polynomial while still preserving the analytic nature of the calculations involved. In all previous models, Λ(w) was merely taken as a constant mass scale Λ.
• Allowing Λ to become a function of the variable w allows us to connect LFWF with PDA algebraically, Eq. (20), without having the need to rather arbitrarily concoct the spectral density.
• Despite having emphasized the previous point, the fact remains that the spectral density can be extracted unequivocally through the knowledge of the PDA.
• Given the most up to date pseudo-scalar meson PDAs, we merely need to fix ν and M q,h . As ν = 1 is a natural choice, we can safely say that the quark mass is the only free parameter to fix the model.
• Crucially, the measure of factorizability of x and k ⊥ in the LFWF is evident through Eqs. (4) and (20). An immediate consequence is the hadronic scale relation between the PDA and the PDF, Eq. (34). This factorization is completely reinstated in the chiral limit, thus reproducing known results [33,34,[37][38][39] as a particular case.
• With the exception of the leading-twist PDAs (which are external inputs) and the charge radii (used as benchmarks to set the values of the constituent masses, M q ), the rest of the distributions and other quantities derived herein are predictions.
Notably, our ingenuous model faithfully reproduces previously known results concerning pions [30][31][32]. Our findings for the kaon are slightly different from those reported therein but can readily and correctly be attributed to the larger strange quark mass favored by this model. However, the description of pion and kaon is compatible with our experimental understanding of these mesons. It is worth mentioning that our pion valence-quark PDF is also compatible with the results from Ref. [74], in which the authors also obtain a NIR for the BSWF, but through the resolution of the corresponding BSE (modeling some of the ingredients that go into the latter). Novel results employing sophisticated mathematical techniques also validate the NIR approach [75]. The distributions reported for η c and η b , and other related quantities, are a completely novel feature of our study. In general, when a comparison is possible, our results also show agreement with other theoretical treatments such as SDEs, lattice QCD, as well as with experimental results. The structure of π 0 , η c and η b is currently being investigated within this model, via two photon transition form factors. In the future we expect to carry out kindred studies of mesons with heavy-light quarks as well as adapt the procedure to study the entangled η−η system and eventually baryons.  (12,(16)(17)(18), it is possible to derive the relation between the PDA and the spectral density ρ M : where the variable y = 1 − 2x has been introduced and we have used the definitions ϕ(y) ≡ φ q M ( 1 2 (1 − y)) and F N = 4 3 π 2 f M n M . The above integral equation can be inverted to a differential equation by differentiating three times, with respect to y, and summing up the resulting equations. This procedure yields an expression for ρ M in terms of derivatives of ϕ: η N ρ M (y) = λ (2) ν (y)ϕ (y) + λ (1) ν (y)ϕ (y) + λ (0) ν (y)ϕ(y) , where η N is a normalization factor such that λ (0) ν (y) = 2νχ 2 + Λ 2 y χ 2 + − 2 1 + (1 − ν)y 2 + ν Λ 2 y +4y(1 − ν) 2Λ 2 y − νχ 2 + Λ 2 y χ + χ − − ν(1 − ν)χ 4 + + 2νχ 2 + Λ 2 y − 8Λ 4 y χ 2 − /Θ y , with the definitions χ ± = (1 − y)Mh ± (1 + y)M q and Θ y = −4 1 − y 2 χ 3 + Λ 4 y . By setting ν = 1, λ (1,0) ν are reduced to λ (1) Furthermore, in the chiral limit: ensuring that our model recovers known result [15]: ⇐⇒ ρ M (w) = ρ asy (w) := 3 4 (1 − w 2 ) .
Beyond the chiral limit, but still keeping the most natural choice ν = 1, the corresponding pion and kaon spectral densities are plotted in Fig. 10. The input PDAs, parametrized according to Eqs. (54), are displayed in the upper panel of Fig. 1.
Though for the purposes of this work (namely computing LFWFs, GPDs and distributions derived therefrom) the determination of ρ M is not required at all, it is worth stressing that the AM we have introduced enables a straightforward derivation of the spectral density from the prior knowledge of the PDA, thus avoiding the need of assuming a particular ad hoc representation for ρ M . This shall be useful for future explorations that require the explicit knowledge of the BSWF, and hence of the spectral density.  Table I.