Observing the Minkowskian dynamics of the pion on the null-plane

A dynamical model is applied to the study of the pion valence light-front wave function, obtained from the actual solution of the Bethe-Salpeter equation in Minkowski space, resorting to the Nakanishi integral representation. The kernel is simplified to a ladder approximation containing constituent quarks, an effective massive gluon exchange, and the scale of the extended quark-gluon interaction vertex. These three input parameters carry the infrared scale {\Lambda}QCD and are fine-tuned to reproduce the pion weak decay constant, within a range suggested by lattice calculations. Besides f{\pi}, we present and discuss other interesting quantities on the null-plane, like: (i) the valence probability, (ii) the dynamical functions depending upon the longitudinal or the transverse components of the light-front (LF) momentum, represented by LF-momentum distributions and distribution amplitudes, and (iii) the probability densities both in the LF-momentum space and the 3D space given by the Cartesian product of the covariant Ioffe-time and transverse coordinates, in order to perform an analysis of the dynamical features in a complementary way. The proposed analysis of the Minkowskian dynamics inside the pion, though carried out at the initial stage, qualifies the Nakanishi integral representation as an appealing effective tool, with still unexplored potentialities to be exploited for addressing correlations between dynamics and observable properties.

The pion plays a pivotal role within Quantum Chromodynamics (QCD), since its Goldstone boson nature is associated with the dynamical generation of the mass of the hadrons and nuclei constituting the visible universe. The pion has a rich structure that stems from the spin degrees of freedom of its constituents, which are necessarily associated with the covariant Minkowski space formulation of QCD. Its dynamics entangles in a conspicuous way the space and spin distributions of the fundamental degrees of freedom within a hadron. The pion is still puzzling by its Goldstone boson nature and its composition in terms of quarks and gluons (see, e.g., Ref. [1]), so that its momentum distributions, the most typical dynamical quantities, have been a target of intense investigation in recent years [2][3][4][5][6][7][8][9][10][11], as well as of planned experimental research at the future Electron Ion Collider.
Hadron imaging is driving experimental [12,13] and theoretical [14] research efforts towards the exploration of the Fock-space structure of the light-front (LF) wave functions, even beyond the valence component. By properly selecting the imaging space, the Fock-space structure of the hadron can be revealed by looking at single-parton distributions with or without spin polarization, double-parton distributions (see, e.g., Refs. [15,16]), triple-parton distribu-tions and in general n-parton distributions. It is clear that such a program, when developed, would provide the ideal framework to study in great detail the Fock-space components of the hadron wave function. Therefore, research efforts to explore the Fock-space content of the hadron LF wave function are necessary, either by using Euclidean Lattice discretization, i.e. lattice QCD (LQCD) (see, e.g., Ref. [17]) or continuous QCD techniques (see, e.g., Ref. [18]).
The challenge on the theory side is to extract from Euclidean calculations the relevant observables defined in the Minkowski space. Vigorous research is pursued to extend LQCD calculations, carried out in Euclidean space, and eventually attain the parton distribution functions (PDFs). With such an aim, several strategies have been proposed, like the one based on (i) the quasi-parton distribution functions (QPDFs) [19] (see, e.g., Ref. [20] for early results), (ii) the pseudo-parton distribution functions (PPDFs) [21] (see, e.g., Ref. [22] for the pion case) and (iii) the so-called lattice cross-section method, as applied in Refs. [10,23]. Another approach is the analytical continuation of the solution of the Euclidean Bethe-Salpeter (BS) equation to the Minkowski space. This method uses the Nakanishi integral representation (NIR) [24] of the Euclidean BS amplitude, allowing to perform the analytical extension to the Minkowski space. Despite that the extraction of the Nakanishi weight-function from the Euclidean calcula-tions constitutes an ill-posed numerical problem and is highly challenging, the NIR method has been applied for obtaining the pion valence PDF from the Euclidean amplitude [25] and the idea was further explored in Refs. [26,27]. Although delicate issues on non perturbative renormalization were pointed out in Ref. [28], there is a possibility of using the NIR to compute QPDFs, allowing to bridge the continuum Minkowski-space QCD to LQCD calculations.
As is well-known [29], the phenomenological description of hadron states on the null-plane can take advantage of a meaningful Fock expansion, once a tiny mass is assigned to the exchanged bosons. Hence the LF approach can usefully exploit the powerful physical intuition based on the Fock space, without the difficulties present in the covariant phenomenological models. In order to formally obtain the LF valence wave function from the BS amplitude, one has to perform its projection onto the null-plane. As a matter of fact, by applying the LF projection to the correlator, built with the minimum number of field operators and providing a non-zero matrix element between the vacuum and the hadron state [30,31], one eventually gets the valence wave function (see also, e.g., Ref. [32] and Appendix A). Indeed, one could generalize the procedure, since to a given hadron state one can associate an infinite number of BS amplitudes, with any number of legs, i.e. quarks and gluons, compatible with the hadron quantum numbers. In turn, each BS amplitude, when projected onto the null-plane, gives the corresponding amplitude of the Fock state with the number of constituents equal to the fields present in the BS amplitude itself. In principle, images of the probability densities obtained from those LF amplitudes shed light on the dynamics inside the hadron with an unprecedented level of detail on the Fock content of the hadron state. However, even the BS amplitude with the minimal number of legs, has information on the full LF Fock-space composition of the hadron [30,31]. It should be recalled that gauge links are always required between the quark field operators in an observable [33], like e.g., a photon absorption amplitude, to keep color gauge invariance, while the BS amplitude by itself is not gauge invariant.
On the theory side, the NIR of the BS amplitude can be a useful tool to solve BS equations in Minkowski space [34][35][36][37][38][39][40][41], and provide the parton distribution amplitudes as well as the elementary fragmentation functions. We have no proof yet that QCD, which embodies confinement, allows such integral representation in the non perturbative domain, although bound state solutions of the BS equation (without confinement) can be actually achieved by using NIR (more precisely, by formally converting the BSE into a generalized eigenvalue problem). In order to use the integral representation to solve the BS equation in its full glory, it is necessary to write the quantities entering the kernel, e.g. the quark-gluon vertex and the propagators, in terms of the NIR. As we mentioned, the nice feature of the NIR is that one can analytically extend it from Minkowski to Euclidean space by performing the Wick-rotation, so that a direct comparison with LQCD results can be feasible. It is worthwhile to point out that within the NIR approach one can prove that the form of the valence wave function for asymptotically large transverse-momentum presents the factorization of the dependences upon ξ (the fraction of the longitudinal momentum) and k ⊥ , naturally recovering the power-law fall-off in the UV region [39,42,43]. Such a property can be extended to higher Fock components of the wave function. Furthermore, at the initial scale, even the simplified calculation of the PDF, based on the Mandelstam formula involving the BS amplitude, will include partons from Fock-state components beyond the valence one. This is an immediate consequence of the solution of the BS equation in Minkowski space.
In the perspective of exploring dynamical models, incorporating as much as possible non perturbative features of QCD in Minkowski space, and to take advantage of the results for building useful hadron imaging, we study the response of the pion valence momentum distribution to the variation of (i) the effective masses of both quark and gluon and (ii) the scale governing the size of the interaction quark-gluon vertex. The variation range of the three input parameters in the ladder BSE of a qq bound system in Minkowski space, is suggested by the corresponding quantities suggested by LQCD calculations (see, e.g., Refs. [44][45][46]).
Our approach relies on the use of the ladder onegluon exchange kernel, assuming that the effect of non-planar diagrams are N c suppressed. As recently shown in a study for bosonic bound states with color degrees of freedom [47], N c = 3 is already large enough to reduce the nonplanar contributions in the structure observables of the bound state to at most 5%, even for strongly bound systems. We solve the BS equation adopting the technique based on (i) the Nakanishi integral representation of the BS amplitude and (ii) the LF projection of the BSE (following the initial elaboration of Ref. [48] and the further developments of Refs. [39,40]). The dynamical inputs in the model are the constituent quark mass, the ef-fective gluon mass, ranging between ∼ Λ QCD /10 and 2Λ QCD , and the size of the quark-gluon vertex, of the order ∼ Λ QCD (see, e.g., Refs. [45,46]). In our covariant model the pion state contains an infinite set of LF components, that are built by a qq and any arbitrary number of effective gluons (as needed to obtain the bound-state pole in the four-points Greenfunction). In order to better understand the influence of those higher-Fock states, we analyze (i) the decay constant, (ii) the valence probability and its spin decompositions, (iii) the longitudinal-and transversemomentum distributions, with a particular attention to the end-point behavior of the first one, (iv) the distribution amplitudes and (v) the 3D image on the nullplane, using both the LF-momentum space and the 3D space, described by the covariant Ioffe-time and the transverse coordinates. It has to be emphasized that apart from the obvious exception of the decay constant, all the other quantities are investigated also in terms of their spin decomposition, opening a window on genuinely relativistic effects inside the pion. Such an extensive analysis represents a distinctive feature of our approach carried out directly in Minkowski space, where the physical processes take place. This is a fundamental step toward a future goal of constructing a framework where Euclidean and Minkowskian studies of the dynamics inside hadrons can be made complementary.
The paper is organized as follows. In Sec. II we outline the formalism for solving the BS equation within the Nakanishi Integral Representation framework, and we briefly illustrate its application to the 0 − bound state. In Sec. III, the valence component of the pion is analyzed in terms of its spin decomposition and an analogous study is extended to the LF-momentum distributions. Section IV is devoted to the formulation of the pion decay constant in terms of the NIR, showing its direct relation with the anti-aligned spin component of the valence wave function. In Sec. V the valence wave function of the pion is investigated in the 3D configuration space associated with the null-plane, in parallel with the 3D representation of the pion obtained in momentum space. In Sec. VI, our wide numerical exploration is presented and discussed, ranging from the pion decay constant to the valence probability (with its spin decomposition), and from the LFmomentum distributions to the 3D pion imaging. We close the work in Sec. VII, drawing our conclusions and presenting perspectives of future developments.

II. THE BETHE-SALPETER EQUATION AND THE NAKANISHI INTEGRAL REPRESENTATION
We briefly summarize the formalism for solving the BSE in Minkowski space for the 0 − quark-antiquark bound state, within the approach based on both the LF-projection technique and the NIR [24]. More details can be attained from e.g. Refs. [39,40], (cf. the equivalent treatment within the covariant LF framework of Refs. [36,48]).
For instance, in the case of a positively charged pion, the BS amplitude and its conjugate are given in the coordinate space by (the translation invariance has been applied) where (i) U and D are fields with quantum numbers corresponding to u and d quarks, respectively; (ii) X = η 1 x 1 + η 2 x 2 , with η 1 + η 2 = 1 (in the present case η i = 1/2); (iii) x = x 1 − x 2 and (iv) p is the total momentum, with M 2 = p 2 the bound-state squared mass. It should be pointed out that since one has also to fulfill the Feynman prescription, as encoded into the chronological operator (if one innocently applies the Dirac conjugation then one gets an anti-chronological ordering). In the momentum space, the intrinsic components are given by Notice that the chronological operator acting in the coordinate space generates the presence of +i in the momentum space. The BS amplitude, Φ(k, p), fulfills the BS equation, that in the ladder approximation reads: where (i) the off-mass-shell constituents have fourmomenta given by p 1(2) = p/2 ± k, with p 2 1(2) = m 2 with m the constituent mass, (ii) p = p 1 + p 2 is the total momentum, (iii) k = (p 1 − p 2 )/2 is the relative four-momentum and (iv) q = k − k the momentum transfer. The Dirac and gluon free propagators are given by Notice that the effective gluon propagator is chosen in the Feynman gauge. Moreover, Γ µ is the interaction vertex and Γ ν (q) = C Γ ν (q) C −1 with the charge operator given by C = iγ 2 γ 0 . In the present model, the quark-gluon extended vertex is described by where g is the coupling constant and Λ a suitable scale for taking into account the size of the color distribution of the interaction vertex. A quantitative estimate of this parameter is suggested by Refs. [45,46].
A. Solving the BSE for the 0 − bound state The most general expressions for Φ(k, p) andΦ(k, p) allowed by the parity and the four-momenta at disposal is [48] (see also Ref. [49]) where φ i are suitable scalar functions that depend upon: {k 2 , p 2 , k · p}, and contain the analytical behavior imposed by the Feynman prescription, i.e. +i . They have to fulfill well-defined properties under the exchange k → −k, according to the anti-commutation rule for the involved fermionic fields. As a result, one has even scalar functions for i = 1, 2, 4 and an odd one for i = 3, when k → −k. The allowed Dirac structures are represented by the 4 × 4 matrices S i , given by The above matrices satisfy orthogonality relations that allow one to reduce the BS equation (4) to a system of four coupled integral equations for φ i (k, p).
The scalar functions can be conveniently written in terms of the NIR as follows: where and g i (γ , z ; κ 2 ) are called Nakanishi weight functions (NWFs) of the scalar function φ i (k, p). Noteworthy, g i (γ , z ; κ 2 ) are real functions which are conjectured to be unique (cf. the theorem on the uniqueness by Nakanishi in Ref. [24]). Those functions encode all the non perturbative dynamical information. The power of the denominator in Eq. (9) can be chosen as any convenient integer ≥ 3, since we are considering the BS amplitude (see also Ref. [48]). The properties of the scalar functions φ i (k, p) under the exchange k → −k can be straightforwardly translated into the corresponding properties of the NWFs g i (γ , z ; κ 2 ) under the exchange z → −z , i.e. they must be even for i = 1, 2, 4 and odd for i = 3.
As is well-known, Eq. (9) allows us to perform the LF projection of the BS amplitude, leading to the valence component of the state (see, e.g., Refs. [32,37] and Sec. III). This motivates the application of the same projection to both sides of BSE, with the aim of determining the NWFs, and eventually the full BS amplitude in Minkowski space. In particular, following Refs. [39,40], one starts from the coupled system of integral equations for the scalar functions φ i (k, p) and arrives at a coupled system for the NWFs, viz where the kernel L ij (γ, z; γ , z ), in the ladder approximation, can be found in full detail in Ref. [40]. It is worth noticing that the two-scalar case was also studied by using an ordinary linear integral equation, where on the lhs there is directly the NWF and on the rhs the folding of the NWF with a suitable kernel. This integral equation was obtained exploiting a uniqueness theorem by Nakanishi [24], assumed to be valid for the non-perturbative case and applied within the LF framework in Ref. [37]. Unfortunately, an analogous treatment for the two-fermion case is hindered by the presence of singularities in the interaction kernel (see below and Refs. [39,40]), making it unclear whether or not the Nakanishi theorem can be formally applied. Hence a more careful analysis is necessary and it will be presented elsewhere. We strongly emphasize that the kernel receives contributions from LF singularities originated by the treatment of the spin degrees of freedom acting in the problem. They were successfully taken into account by means of the methods developed in [39] (see Ref. [50] for a previous discussion of those singularities). The above set of integral equations is solved numerically by matrix manipulation algorithms, after expanding the NWFs onto the Cartesian product of Laguerre polynomials, for the γ dependence, and Gegenbauer ones, C λi n with suitable λ i , for the z dependence.

B. Normalization
In order to calculate hadronic properties, in our case the valence probability and momentum distributions, the BS amplitude has to be properly normalized. In the ladder approximation, the normalization reads [51] T r It is worth noting that such a normalization can be easily reverted into the charge normalization, within the adopted ladder approximation.
By using Eq. (7) and performing the Dirac traces, the normalization condition turns to be: where N c is the number of colors and b = (k · p) 2 − k 2 M 2 /M 4 . By introducing the amplitudes φ i given in terms of the NIR, Eq. (9), one can straightforwardly perform the analytical integration of the momentum-loop using Feynman parametrization in Eq. (13), and finally get where Even in the ladder approximation, the normalization, Eq. (13), contains the contributions beyond the valence one from the Fock expansion of the pion state, i.e. it takes into account the infinite sum of states with a quark-antiquark pair and any number of gluons (see e.g. [31,52]).
We should call the attention of the reader to the fact that from the normalization condition Eq. (13), one can arrive at the normalization fulfilled by the amplitudes of the pion-state Fock components (see, e.g., Refs. [32,37]), so that a probabilistic framework can be restored. In particular, the normalization of the valence component is nothing else than the probability to find the component with the lowest number of constituents inside the pion.

III. VALENCE PROBABILITY AND LF MOMENTUM DISTRIBUTIONS
The valence probability and momentum distributions can be derived resorting to the LF quantum-field theory methods (see, e.g., Ref. [29]), where one defines the creation and annihilation operators for particles and antiparticles with arbitrary spin on the null-plane, in order to construct the generic LF Fock state. Actually, the valence component is defined by where b(q 2 , σ 2 ) is the quark annihilation operator and d(q 1 , σ 1 ) is the antiquark one, (16) is related to the BS amplitude by means of the LF projection (see the details in Appendix A) as follows From Eq. (17), one ultimately recognizes that the evaluation of the valence wave function comes from the elimination of the relative LF time between the quark operators entering the BS amplitude (see also Ref. [32]). Alternatively, the valence wave function can be obtained using the quasi-potential expansion method adapted to perform the LF projection of the BS equation and amplitude (see Refs. [30,31,52,53] for details).
As above mentioned, the Fock expansion of the interacting system allows one to restore a probabilistic framework that the BS amplitude is not able to yield. In particular, one can get the valence probability given by where z = 1 − 2ξ = −2k + /M . As it is shown in detail in Appendix A, the valence component can be decomposed into two spin contributions, given by the configurations we indicate as anti-aligned and aligned. The first configuration corresponds to a total spin of the quark-antiquark pair S = 0, while the second one pertains to the S = 1 spin-state. It is important to emphasize that the anti-aligned configuration, yielding the largest contribution to the pion state, is coupled to the eigenstate of the operator L z with eigenvalue z = 0. Within the LF framework, the z-component of the orbital angular momentum is diagonal in the Fock space, being a kinematical generator. Differently the aligned configuration necessarily invokes the eigenvalues z = ±1, given the pseudoscalar nature of the pion. Interestingly, the presence of the aligned-spin contribution is an unavoidable and clear signature of the relativistic dynamical regime inside the pion, since the z = ±1 orbital component is related to the small component of the spinors in the Dirac theory. It is clear that a quantitative study of the aligned component has a pivotal role in understanding the relativistic features of the light hadrons (and the possible relativistic corrections for the heavier ones). Finally, it should be pointed out that the z = ±1 contribution is perfectly compatible with the requested parity conservation (recall the minus signs in the matrix γ 0 entering the representation of the parity operator). The obtained expression for ϕ 2 is where and In the above equations, ψ i are given by Notice that for the two spin configurations, S = 0 and 1, one has the suitable dependence upon the angle in the transverse plane, as dictated by the eigenvalue z , 0 and ±1, respectively. After inserting Eq. (19) into Eq. (18), one can write the valence probability in terms of the valence momentum-distribution density, P val (γ, z), i.e. where with the anti-aligned and aligned probability densities defined by and Recall that |k L(R) √ 2| 2 = |k ⊥ | 2 = γ. The valence longitudinal and transverse LF momentum distribution densities are obtained by properly integrating the valence probability density P val (γ, z). The longitudinal-momentum distribution, with its spin decomposition, is given by with For the transverse-momentum distribution one has with It should be pointed out that φ(ξ) is the unpolarized structure function, which one can access in the deep inelastic limit of the leading order virtual photon absorption process, as illustrated by the diagram on the left side of Fig. 1. In this case, where the final states are given by qq plane waves, the description of the inclusive process that allows one to extract the longitudinal distribution can be obtained by integrating on the final states, so that its pictorial representation is yielded by the box diagram for on-mass-shell quarks, as discussed in Ref. [54] (for a description of the deepinelastic scattering in an exactly solvable LF model with a confining interaction, see, e.g., Ref. [55]). The right side of Fig. 1 represents the contribution from the one-gluon exchange in the final state, that comes from the expansion of the Wilson line, needed for assuring the color gauge invariance of the quark correlator entering in the description of the deep-inelastic process (see, e.g., Ref. [33]). It has to be pointed out that such a contribution in the final state is also necessary for obtaining non-vanishing T-odd transverse momentum distributions (TMDs) [56].

IV. DECAY CONSTANT
A basic observable that one has to reproduce for assessing a given approach is surely the pion decay constant, f π . It is defined in terms of the BS amplitude by (see, e.g., Ref. [57] for details) Contracting with p µ and using the decomposition of the BS amplitude given by Eq. (7), one can perform the trace and obtain It is worth noting that the decay constant is determined only by one component (even under the exchange 1 → 2) of the BS amplitude. By using LF variables, one can manipulate Eq. (22) and get Hence, within the NIR approach the final expression for the decay constant reads (34) It should be recalled that g 2 is properly normalized through Eq. (14).
Equivalently, one can proceed from Eq. (32), by carrying out the 4D integration on the NIR of the scalar function φ 2 (k, p) (see Eq. (9)) and eventually obtaining the result in Eq. (34). As a matter of fact, one has where, after applying the change of variable q = k + z p/2, a Euclidean 4D integration has been performed, since the analytic structure is known. Hence, inserting in Eq. (32) first Eq. (9) and then the result given in Eq. (35), Eq. (34) is re-obtained. Interestingly, one can re-express f π in terms of the anti-aligned component, i.e. S = 0, given in Eq. (20), once an integration by part of the third term is performed, viz.
It is worth noticing that this expression is expected, since the famous argument applied for explaining the prevalence of the muonic channel in the pion decay is based on both the helicity conservation applied to a S = 0 state and the phase space constraint.

V. THE PION IMAGE ON THE NULL-PLANE
The analysis of the deep-inelastic scattering (DIS) is usually carried out in the infinite-momentum frame. In addition, one can study the DIS processes in the target frame, adopting the configuration space. Then, one is able to establish a framework where a more detailed investigation of the space-time structure of the hadrons can be performed. In particular, a fruitful and far-reaching example of the above mentioned program is the introduction of the so-called Ioffe-time [58][59][60]. In the laboratory frame, it has a relevant role in addressing the issue of a quantitative description of the interplay between dynamical regimes governed by short-and long-range interactions. Indeed, one can also introduce a covariant realization (see, e.g., Ref. [61]) of the Ioffe-time, with the same aim of studying the relative importance of short and long light-like distances, e.g. probed in DIS processes. This covariant expression is particularly useful for the analysis of the valence wave function, which, like all LF amplitudes in the Fock expansion of the bound state, is a LF-boost invariant quantity and properly encodes information on the dynamics inside the hadron. It should be pointed out that we use the Ioffe-time in the context of the BS amplitude, i.e. the transition amplitude from a hadronic state to the vacuum, while, e.g., in Refs. [21,62] the focus is on the generalized parton distributions and/or the transverse-momentum distributions.
As is well-known, the null-plane x + = t + z = 0 defines a three-dimensional hypersurface, with space-like distances, where the LF wave functions live. The configuration space associated with this hypersurface has longitudinal coordinate x − = t−z and transverse position x ⊥ . The coordinate x − is called (c = 1) Ioffe-time [58][59][60], with the covariant version given by x·p that in the target frame becomes x − p + /2, on the null-plane. Notably, in the target frame, when the DIS regime is reached, the Ioffe-time measures the light-like distance between the production of a qq pair by the virtual photon and its interaction inside the hadron. Moreover, one can quickly realize that the Björken x Bj , or more precisely, the longitudinal fraction ξ = p + 1 /p + , is the variable conjugated to the Ioffe-time, once the Fourier transform of the matrix elements of the electromagnetic current correlator is analyzed for obtaining DIS structure functions. As suggested by phenomenological analyses (see, e.g., Refs. [63,64]), one realizes that the Ioffe-time (also indicated as the coherence length of the qq pair) is ∝ 1/M x Bj . This follows from the energy-time uncertainty, involving the qq pair off-shell energy and the time interval between its production and interaction with the hadronic medium, in the target frame. Hence, longer and longer coherence lengths pertain to the values of x Bj , where the QCD dynamics is dominated by the IR regime. This enforces the relevance of quantities that depends upon the Ioffe-time, when we aim at disentangling the different light-like distances that the virtual photon probes, and eventually shedding light on, e.g., higher Fock states production, onset of the confinement, valence structure, etc.
Obviously, one can study the valence amplitude also in the configuration space, where the dependence results to be upon the coordinates {z = x − p + /2, b} [61] (notice that the LF-boost invariant definition of the Ioffe-time has been used and b ≡ x ⊥ ). Rather than the Fourier transform of ϕ 2 (ξ, k ⊥ , σ i ; M, J π , J z ), it is physically more interesting to address the distribution probability in the coordinate space, in analogy with the study carried out in the LF-momentum space, where the distribution is given in Eq. (24) and fulfills the sum rule (23). In view of this, it is better to consider the Fourier transform of the two spin components ψ ↑↓ and ψ ↑↑ , Eqs. (20) and (21), that reads where (i) z = 0, 1, for S = 0, 1 respectively, and (ii) for x + = 0 the scalar product x · k reduces to x · k = x − k + /2 − b ⊥ · k ⊥ (in our convention). The amplitudes ψ ↑↓(↑↑) (|k ⊥ | 2 , z) vanish for z outside the interval [−1, 1].
Collecting the results presented in Appendix B, one can write the Fourier transform of the two spin components, Eqs. (20) and (21), in terms of auxiliary amplitudes, where the leading asymptotic behavior for large b is factorized out, i.e.
with κ given in Eq. (10) and (recall that the product z g 3 (γ , z; κ 2 ) i s even in z) Notice that χ ↑↓ and χ ↑↑ are symmetric byz → −z, i.e. the inversion of the light-like axis. We will provide results for the amplitudes χ ↑↓ (z, b) and χ ↑↑ (z, b), where the exponential fall-off in b is factorized out, for the sake of presentation.
From the above elaboration, one can obtain the probability density in the 3D spacez ⊗ b ⊥ , that reads and it fulfills by construction the following relation Finally, we observe that the null-plane components in Eq. (38) atz = 0 can be directly obtained from Euclidean-space calculations, once the spin components of the transverse amplitude are defined as follows The above quantities come from the integration over 1 2 dk + dk − of the BS amplitude (leading to x + = x − = 0). Notably, given the analytic properties of the NIR, an equal result can be obtained if one integrates on ı dk 0 E dk 3 , where k 0 E is the Euclidean momentum component (see Ref. [42] for the analytical details). Thus, it could be of interest to compare the transverse functions obtained by direct calculations in the Euclidean space for the pion BS amplitude with the ones evaluated by solving the BSE in Minkowski space through the NIR.

VI. QUANTITATIVE STUDY
This section is devoted to a wide presentation of a quantitative study of the charged pion structure, carried out within the approach previously outlined. In particular, we discuss results for (i) the decay constant, (ii) the valence probability and its spin decomposition, (iii) the valence longitudinal-and transversemomentum distributions, and finally (iv) the 3D image of the pion, in the space described by the Ioffetime and the transverse coordinates, i.e. {z, b}.
It is important to remind that the pion BSE amounts to a set of four coupled integral equations for the scalar functions φ i (k, p) present in Eq. (7) (see Refs. [39,40] for details), and it is formally transformed into a set of coupled integral equations for the four NWFs g i (see Eq. (11)), within the NIR approach. In turn, this set of equations can be expressed as a generalized eigenvalue problem. In order to accomplish this step, we exploit the property of the NWFs to be real functions, depending upon real variables (one compact and the other non compact) by expanding them onto a bi-orthonormal basis (for details see [40]).
As discussed in Sec. II, while solving the BSE in ladder approximation, one has three input parameters: (i) the constituent-quark mass and the gluon effective one, indicated by m and µ respectively; (ii) the scale Λ of the extended interaction vertex. The values of the adopted parameters cover a fairly broad spectrum, including the values inspired by LQCD results. In practice, we have considered eleven sets of values, where m ranges from 187 to 255 MeV, µ from 0.15 m to 2.5 m (i.e. from about 30 to 600 MeV) and Λ from m to 2m, to be of the same magnitude Λ QCD as suggested in [45,46].
For completeness, we specify the basis and truncation scheme adopted in our calculations and our numerical accuracy. We use for each NWF an expansion of the form where A i kn are the coefficients to be determined and the functions G λ n , and L n are defined by G λ n (z) = (1 − z 2 ) (2λ−1)/4 Γ(λ) n!(n + λ) 2 1−2λ πΓ(n + 2λ) C λ n (z), where C λ n denotes a Gegenbauer polynomial and L n is a Laguerre polynomial.
Moreover, in order to speed up the convergence of the Gaussian integration, it has been chosen a = 6/m 2 . It should be noticed that because of the symmetry under z → −z one has The basis functions defined by (44) obey the orthogonality relations In our calculations, each index λ i is a half-integer, i.e. λ i = i + 1/2, with { 1 , 2 , 3 , 4 } = {2, 4, 6, 6}, representing the best choice after checking the numerical convergence, also varying the number of polynomials in the basis. For obtaining the results presented in this Section, it was used up to N γ = N z = 60 basis functions for each g i (γ, z; κ 2 ). The checked accuracy in the coupling constant is at the shown significant digits; for the valence probability the accuracy is three significant digits; for f π the relative accuracy is better than 0.1%; for φ ↑↓(↑↑) (ξ), the point-wise accuracy for 0.95 > ξ > 0.05 is about 3 significant digits and decreasing towards the end-points, i.e. within the width of the lines shown in the following figures.

A. Static properties: the decay constant and the valence probability
In our calculations of the pion structure, as it is customary in solving BSE, one assigns a value to the quark mass m and using m π = 140 MeV, one gets the binding B, defined by B = 2m − m π . Hence, once we select a value for both the gluon mass µ and the interaction-vertex scale Λ, we can proceed in solving the generalized eigenvalue problem where the quarkgluon coupling constant, g 2 , is the eigenvalue (see Eq. (11), and notice that beyond the ladder approximation one has to cope with a more general eigenvalue problem [43,65]). For the eleven sets of the three input parameters, as shown in detail in Table I, we have evaluated the following main quantities: (i) the adimensional coupling that combines the bare coupling g 2 /4π and the factors from the two extended interaction vertexes (see the definition in Eq. (6)), (ii) the valence probability, and its spin decompositions and (iii) the charged-pion decay constant.
In Table I we also inserted two other quantities, useful for sharpening our physical insights. The first one is the effective strengthᾱ, given bȳ where the value of the average transverse-momentum k 2 ⊥ ∼ √ 0.2 m has a close correspondence to the characteristic scale of the decreasing behavior shown by the transverse-momentum distribution in the model, as it will be clear when presenting results for this quantity in the next subsection. Loosely speaking, the denominator µ 2 + k 2 ⊥ plays the role of an effective mass carried by the gluon in the interaction region relevant for binding the qq system. The second quantity is the adimensional ratio f π /m, that yields an estimate of the strength of the effective quark-pion coupling introduced in low-energy effective approaches, like the Nambu-Jona-Lasinio one (see, e.g., Ref. [66]). Table I is organized according to increasing values of the valence probability, corresponding also to decreasing values of the effective couplingᾱ, since both of them point to the same physical mechanism, as discussed below.
First of all, after slightly tuning the parameters in the range suggested also by LQCD (see the eighth set), one is able to reproduce the experimental value of the pion decay constant, i.e. f exp π ± =130.50(1)(3)(13) MeV [67], as well as the LQCD average value f LQCD π ± = 130.2 (0.8) MeV, as given in Ref. [68]. By varying the sets of parameters, the valence probability P val , and f π , run in the interval [0.64 − 0.71] and [77 − 134] MeV, respectively.
The valence probabilities for the anti-aligned and aligned spin configurations are also shown in Table I, where one recognizes the expected prevalence of the spin S = 0 component, ranging from 0.55 to 0.58, and the minor role of the S = 1 configuration, from 0.09 to 0.14. In any case, it is worth noticing that the S = 1 contribution, which is exclusively relativistic in nature, is by no means negligible, since the relative weight increases up to 30%. From the Fock-expansion standpoint, the size of such a relative weight indicates that the higher-Fock states are quite relevant in the description of the pion state on the null-plane. In fact, the ladder kernel of the BSE when projected onto the LF hyperplane [30,31,52,53] takes into account an infinite number of Fock-components beyond the valence state, built as a qq pair and any number of gluons.
The remaining probability, 1 − P val , is distributed among the first Fock components beyond the valence one, and one should notice the consistent picture that emerges from observing the correlation between the decreasing values of the valence probability and the increasing values ofᾱ. This behavior is rather natural to be expected, asᾱ weights in an effective way the coupling to the higher Fock states present in the dynamical model. Consequently, the largerᾱ, the smaller P val , since more gluons can be present in the intermediate states, and the valence state becomes less likely. As to the pion decay constant, while f π does not show a regular behavior whenᾱ is decreasing, the ratio f π /m has less pronounced variations, since this ratio combines the effect of the higher Fock states through two different quantities. On one side, f π is associated with the anti-aligned component of the valence wave function at the origin, i.e.z = b x = b y = 0 (see Eq. (36)), and its increasing indicates a major role of more compact configurations, necessarily related to the higher Fock states. On the other side, the quark mass determines the binding energy, and larger values of B = 2m−m π are related to smaller size of the pion, i.e. more compact configurations take place (cf. the values of m andᾱ for the sets III, IV and V, as well as IX, X and XI, where µ and Λ do not vary). Hence, the values of f π show a dependence upon quark mass, since both quantities are driven by the relevance of states beyond the valence one.
In conclusion, Table I provides two interesting insights, which can be briefly highlighted: (i) the decay constant is influenced by the compact configurations and the pion size, thus UV and IR properties are reflected inits actual value, with the conspicuous relation to the constituent quark mass, and (ii) also the valence probability encodes signatures of the UV and IR regimes. This will become more clear once the behavior at the end-points of the longitudinal distribution and the high-momentum tail of the transverse momentum distribution are analyzed, and recalling that f π and P val are obtained from both distributions after properly integrating..

B. The valence LF-momentum distributions
In this subsection, we present the dynamical quantities predicted for the pion on the null-plane: (i) the longitudinal-momentum distribution, φ(ξ), Eq. (27), and the transverse one, P(γ), Eq. (29), (ii) the valence LF-momentum density, P(γ, z), Eqs. (24), (25) and (26), (iii) the distribution amplitude, with its spin decomposition, and its transverse counterpart. The calculation of such quantities represents the needed first step, followed by the suitable evolution, for meaningful comparisons with experimental results, e.g., obtained in DIS and semi-inclusive DIS processes, like the pion parton distribution function and the transverse-momentum distributions, as well as extracted from the deeply-virtual Compton scattering, like the generalized parton distributions.  Let us first consider the longitudinal-momentum fraction distribution for the valence LF component, as a function of the quark longitudinal fraction ξ, and its anti-aligned and aligned components. It has to be stressed that on the null-plane, the longitudinalmomentum fraction distribution has the proper support, and therefore it fulfills both normalization and momentum sum rule when we take into account the whole set of the Fock components. This is a remarkable benefit for correctly analyzing DIS processes.
The calculations are shown in Fig. 2, for some illustrative cases, with the normalization equal to 1 (i.e. each distribution is divided by the respective probability) and the model parameters given in Table I Figure 2 shows that the decrease of the effective dimensionless strength of the kernel,ᾱ, broadens φ(ξ) with a regular pattern, which is not too sensitive to the wide variations of both gluon mass and vertex parameter. Notice that even in presence of a very light gluon, the valence longitudinal distribution does not show dramatic changes, although they are visible. As expected, given its genuinely relativistic nature, the aligned distribution has a wider shape than the antialigned one, namely larger values of the quark mo-mentum prevail, but the overall impact on the total distribution is mitigated by the smaller probability associated with the S = 1 configuration. This observation is quantitatively corroborated by the analysis summarized in Table II, where it is shown, for a few examples, the correlation betweenᾱ and the exponent of the function (1 − ξ) η , which we use, as customary, for fitting the valence longitudinal-momentum distributions close to the end-points.
Notice that the distributions are symmetric at the initial scale, while after evolving they cumulate at values ξ ≤ 0.5 (as it will be discussed elsewhere). As pointed out in the previous subsection, the decreasing values of the valence probability and the corresponding increasing value ofᾱ were put in relation with the increasing role of compact configurations, generated by higher Fock states. Plainly, also the increasing values of the exponent η in the fitting functions can be explained through an analogous mechanism. In fact, the growth of the powers, shown in Table II, depletes the valence quark distribution in the UV region, corresponding to ξ → 1. The function (1 − ξ) η becomes smaller and smaller for increasing values of η, and the probability of large longitudinal momenta decreases. This implies a suppression of configurations with the valence qq pair close together.
As already noticed, the aligned longitudinalmomentum distribution is broader than the antialigned one (cf. Fig. 2). This corresponds to a  softer end-point behavior, with exponents systematically smaller than the anti-aligned ones, about η ↑↓ − η ↑↑ ∼ 0.2, as shown in Table II.
The results for the Mellin moments of the valence longitudinal-momentum distributions are shown in TABLE II. An excerpt from the set of exponents of the fitting function (1 − ξ) η for ξ → 1, while varying the set of input parameters. The three columns, labeled by η ↑↓ , η ↑↑ and η, correspond to the anti-aligned, aligned and total valence longitudinal-momentum distributions (cf. Fig. 2 Table III, also with the decomposition in spin contributions. It should be recalled that the present approach predicts a quite robust longitudinal fraction carried by the valence qq pair, according to P val that ranges from 60% to 70%. This can be assessed from the first Mellin moment, after noticing that if the va-lence φ(ξ) is normalized to 1, and not to P val , the value of x is always 0.5, due to the symmetry of the valence wave function around ξ = 0.5. It is interesting to observe that the moments do not change too much for the set of input parameters we have chosen, and the increasing values correspond to the increasing relevance of the momentum distribution close to the end-points. Such a behavior is clearly more evident for Mellin moments associated to higher powers. To carry out a detailed comparison with LQCD results and the analogous outcomes from the continuum QCD approach (see, e.g. Ref. [18]), one has to determine first the initial scale of the valence longitudinalmomentum distribution and then apply the evolution (see, e.g., Ref. [69]).
The second dynamical quantity we have analyzed is the transverse-momentum distribution, P (γ) (Eq. (29)), and its decomposition in spin configurations. The results, with normalization equal to 1, are presented in Fig. 3 for the same set of parameters adopted for the longitudinal-momentum distribution in Fig. 2. The same general behavior found for the longitudinal-momentum distribution can be also recovered for the transverse-momentum one, namely, the size of the tail at large values of γ becomes smaller for larger values ofᾱ, and vice-versa. As we learned, this is correlated to the role played by the typical size of the spatial correlation of the valence qq pair.
The anti-aligned and aligned distributions are also shown in Fig. 3. It is worth noticing that the last one has a slower decay compared to the anti-aligned case, since it carries a factor of γ, stemming from the relativistic spin-orbit coupling, which eventually adds a factor k ⊥ to the momentum dependence of the aligned component in the valence wave function. The typical momentum scale that determines the effective strength,ᾱ, i.e. γ/m 2 ∼ 0.2, can be inferred from the behavior of P (γ)/P (0) close to γ = 0, as shown in Fig. 3. More explicitly, such value of γ roughly indicates the region where P (γ) is relevant for obtaining the actual value of the valence probability. Hence, the above mentioned range of γ can be associated with the size of the region where the interaction is effective in building the qq bound state (in momentum space). In other words, such a low-energy or IR scale should govern the kernel of the BSE, so that it is effective in giving the strength necessary to create the strongly bound qq system, resulting in the pion.
Summarizing the first part of the analysis, one should point out that the end-point behavior of φ(ξ) is strongly correlated to the UV properties of the adopted kernel in the BSE, while the range of γ, where P (γ) is large, is governed by the size of the bound state, i.e. by the IR behavior of the interaction.
An overall view of the pion in the 3D LF-momentum space can be obtained from Fig. 4, where the LFmomentum density P(γ, ξ) (cf. Eqs. (24), (25) and (26)), is shown. One should recall that the longitudinal and transverse distributions are obtained through the suitable integration of the density P(γ, ξ).
Finally, we show the distribution amplitude (DA) [29,[70][71][72]. This quantity is introduced through the factorization of the amplitudes associated with exclusive processes and can be expressed as an integral on the transverse-momentum dependence of the valence wave function. In particular, we have evaluated the following spin decompositions Analogously, we have introduced the transverse distribution amplitude (TDA) by integrating the valence wave function over the fraction of longitudinal momentum carried by the valence quark (recall z = 1 − 2ξ), viz It has to be pointed out that the TDA is the Fourier transform of Eq. (42), namely the transverse amplitude in the transverse-coordinates space. The TDA can be also obtained from Euclidean-space calculations (see, e.g,. Ref. [42]). The results for the two spin configurations of both DA and TDA, obtained by using the parameters of the set VIII, are shown in Fig. 5. It is interesting to observe that the aligned component of the DA is wider and decreases slower at the end-points than the antialigned component, as it happens for the longitudinalmomentum distributions (c.f. Table II  region is governed by the one-gluon exchange, i.e. the short-range interaction, the IR region incorporates the features dictated by the long-range correlations in the transverse-coordinate space. In the 3D space described by the Ioffe-time and the transverse coordinates, one can obtain an image of the pion in terms of the two spin components of the valence wave function, i.e., χ ↑↓(↑↑) (z, b), given in Eq. (39). Such a picture of the pion allows one to understand better the interplay between short and long light-like distances in the description of the hadron structure (see Sec. V). In view of this, one notices that the region with small values of {z, b} is the place where the UV effects should manifest. Beside the 3D image, we also present the transverse amplitudes, ϕ T ↑↓ (b) andφ T ↑↑ (b), given in Eq. (42), since they could be the target of LQCD studies. To perform the calculations shown in this subsection, we have used the parameter set VIII (see Table I) that fits the pion decay constant.
The 3D image of the pion spin components on the null-plane is provided in Fig. 6. Notice that the exponential factor exp(−κ b), present in the Fourier transform of ϕ 2 (ξ, k ⊥ , σ i ; M, J π , J z ) (cf. Eqs. (B2), (B3) and (B4)) is factorized out in both χ ↑↓(↑↑) (z, b), al- lowing to use a linear scale in the 3D plot. For the purpose of the figure, each component is multiplied by the transverse coordinate b, canceling a log-type singularity at b = 0, generated by the Bessel function K 0 . A general feature of both densities is the sharp enhancement forz = 0, i.e., at vanishing light-like distances. Inspired by such an enhancement, and in order to better analyze the physically significant dependence uponz of the valence wave function, we have also studied the absolute value squared of the integrals withψ ↑↓(↑↑) (z, b) given in Eq. (38). These amplitudes have a more direct link to the spin components of the valence LF wave function since they also contain the original exponential factor e −κ b . The integrated amplitudes in Eq. (51) describe each spin configuration where the constituents are at light-like distancez and have relative transversemomentum k ⊥ = 0. As shown in Fig. 7, the amplitudes Ψ ↑↓(↑↑) (z) 2 , have a nice diffraction pattern, that represents a peculiar feature of such quantities, emphasizing the interference content due to the entangled qq pair. Finally, in Fig. 8, we present the Fourier transform of the transverse amplitudes, ϕ T ↑↓(↑↑) (γ) (Eq. (50)), in the transverse-coordinate space, i.e., .
The log-plot allows one to recognize the characteristic exponential decay at large distances, clearly dominated by exp(−κb), which is the familiar fall-off of the wave function of a bound state (κ 2 = m 2 −M 2 /4 > 0). Since the amplitudes depend upon transverse coordinates, they should be accessible to LQCD. Therefore, it could be interesting to compare its predictions with phenomenological calculations in order to get information about the highly nonlinear behavior of QCD at large transverse distances.

VII. CONCLUSIONS AND PERSPECTIVES
Within the light-front framework, where the physical intuition based on the Fock space expansion of the hadron wave function can be used at large extent, we have studied the strongly bound qq system that generates the pion. In our approach, the Bethe-Salpeter equation in Minkowski space is solved by using the Nakanishi integral representation of the BS amplitude. The ladder approximation of the BS interaction kernel in the Feynman gauge has been considered, with the constituent quarks interacting by an effective massive gluon exchange. An extended effective quark-gluon vertex function is introduced through a form factor, that contains a new scale parameter, beside the gluon mass. The range where (i) the form-factor parameter Λ, (ii) the quark effective mass and (iii) the gluon one can vary has been chosen according to LQCD results, with the guideline given by the value of the IR scale Λ QCD . We fine-tuned the parameters around such a scale in order to reproduce f π . Once the BSE has been solved, we formally obtained the valence LF wave function, that is well-defined in terms of the lowest number of fields associated to the BS amplitude, and therefore it contains unambiguous information on the dynamics that one can convey into the interaction kernel of the BSE itself.
The detailed study of the valence Fock component has been carried out by first addressing static (or integral) quantities, like: (i) the pion decay constant, f π , and (ii) the valence probability, whose value is around 70%, according to our calculations. A relevant and peculiar feature of our approach is the pos-sibility of decomposing the valence component in the two allowed spin configuration: the dominant S = 0 component and the purely relativistic S = 1, that yields a ∼ 20% contribution to the valence probability. Such a decomposition also has been applied to the second set of investigated quantities, bringing a considerable wealth of dynamical information. In particular, we have analyzed: (i) the longitudinal and transverse LF-momentum distributions; (ii) the distribution amplitudes that depend upon the longitudinal and transverse LF-momentum components; (iii) the probability densities in momentum and configuration spaces (the last one is given by the Cartesian product of the covariant Ioffe-time and the Euclidean transverse coordinates). For each quantity, we have highlighted signatures of the dynamics governed by the one-gluon exchange and have also emphasized the relevance of the transverse degrees of freedom, more accessible to the LQCD studies.
Future developments of our approach are primarily related to the implementation of both quark and gluon dressed propagators, more realistic dressed quarkgluon vertex, in order to extend to the Minkowski space the successful studies of the spontaneously broken chiral symmetry performed in the Euclidean space. In this way, the kernel of the BSE can be improved systematically, adding step by step new dynamical contents from QCD, that one should explore in the elaboration of hadron models. Furthermore, observables like the electromagnetic form factor and generalized transverse-momentum dependent parton distributions (see, e.g., Ref. [73]) of the pion are within our plans for future studies.  (q 1 ,q 2 ) = T r P σ1σ2 (q 1 ,q 2 ) γ 5 Λ + D σ1σ2 2 (q 1 ,q 2 ) = T r P σ1σ2 (q 1 ,q 2 ) k ⊥ · γ ⊥ γ 5 Λ + , with P σ1σ2 (q 1 ,q 2 ) = v(q 1 , σ 1 )u † (q 2 , σ 2 ) . (A10) The spinors u and v in terms of the LF variables can be found in Appendix B of Ref. [29], namely u(q, σ) = q + + γ 0 m + γ 0 q ⊥ · γ ⊥ 2q + χ σ σχ σ (A11) and v(q, −σ) = q + − γ 0 m + γ 0 q ⊥ · γ ⊥ 2q + χ σ σχ σ , with σ = ±1. Notice that in v there is an opposite sign for the helicity (in the LF framework helicity and third component of the spin coincide, see, e.g., [74]). For a quick evaluation of the above matrix elements in Eqs. (A9), it is useful to introduce the following dyadic products χ σ σχ σ χ σ † σχ σ † = (1 + σγ 5 ) Λ + , (A13) where in Eq. (A14) for σ = ±1 one has to use γ L(R) , with The actual expressions of D 1 and D 2 can be obtained through suitable traces. In particular, for D −σσ 1 one gets where we can properly move the projector Λ + inside the trace for simplifying the whole expression, e.g.
Also D σσ 1 is vanishing, once the following relation is adopted and the standard results for the traces of the Dirac matrixes. Namely, one gets D σσ 1 (q 1 ,q 2 ) = T r The last matrix element D σσ 2 is given by where where J m (x) is the Bessel function of integer order. In particular, the Fourier transform of the spin antialigned and aligned components are associated with integrals of J 0 and J 1 respectively, after performing the angular integration in Eq. (37). By adopting the expression of ψ ↑↓ (γ, z) and ψ ↑↓ (γ, z) in terms of the respective NIRs (cf. Eqs. (20) and (21)), one has to evaluate the following integrals. The first one is where K n (x) is the modified Bessel function of the second kind. The other two integrals are and Notice that F 0 , F 0 and F 1 depend upon z 2 , and this allows one to eliminate odd functions when integrating on z in Eq. (37). The driving exponential fall-off of F 0 , F 0 and F 1 in the asymptotic limit b → ∞ comes from K m (x), which reads: Hence, the leading exponential behavior in the integrals (B2), (B3) and (B4) comes from values of e −b γ +κ 2 +z 2 M 2 4 (as seen from Eq. (B5)) with γ close to 0 and z ∼ 0, namely F 0 (z ∼ 0, γ ∼ 0, b)| b→∞ = = b −1/2 F 0 ((z ∼ 0, γ ∼ 0, b)| b→∞ = = b −1 F 1 ((z ∼ 0, γ ∼ 0, b)| b→∞ ∼ e −b κ , (B6) with κ given in Eq. (10).