Fermionic bound states in Minkowski space: light-cone singularities and structure

The Bethe–Salpeter equation for two-body bound system with spin 1/2 constituent is addressed directly in the Minkowski space. In order to accomplish this aim we use the Nakanishi integral representation of the Bethe–Salpeter amplitude and exploit the formal tool represented by the exact projection onto the null-plane. This formal step allows one (i) to deal with end-point singularities one meets and (ii) to find stable results, up to strongly relativistic regimes, which settle in strongly bound systems. We apply this technique to obtain the numerical dependence of the binding energies upon the coupling constants and the light-front amplitudes for a fermion–fermion 0+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0^+$$\end{document} state with interaction kernels, in ladder approximation, corresponding to scalar-, pseudoscalar- and vector-boson exchanges, respectively. After completing the numerical survey of the previous cases, we extend our approach to a quark–antiquark system in 0-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0^-$$\end{document} state, taking both constituent-fermion and exchanged-boson masses, from lattice calculations. Interestingly, the calculated light-front amplitudes for such a mock pion show peculiar signatures of the spin degrees of freedom.


Introduction
The standard approach to the relativistic bound-state problem in quantum-field theory was formulated, more than a half century ago, in a seminal work by Salpeter and Bethe [1]. In principle, the Bethe-Salpeter equation (BSE) (see also the review [2]) allows one to access the non-perturbative regime of the dynamics inside a relativistic interacting system, as As is well known, apart the celebrated Wick-Cutkosky model [3,4], composed of two scalars exchanging a massless scalar, solving BSE is very difficult when one adopts the variables of the space where the physical processes take place, namely the Minkowski space. Furthermore the irreducible kernel itself cannot be written in a closed form. Nonetheless, in hadron physics it could be highly desirable to develop nonperturbative tools in Minkowski space suitable for supporting, e.g., experimental efforts that aim at unraveling the 3D structure of hadrons. It should be pointed out that the leading laboratories, like CERN (see Ref. [5] for recent COMPASS results) and JLAB (see, e.g., Ref. [6]), as well as the future electron-ion collider, have a dedicated program focused on the investigations of semi-inclusive DIS processes, i.e. the main source of information on the above issue.
Moreover, on the theory side we mention the present attempts of getting parton distributions as a limiting procedure applied to imaginary-time lattice calculations, following the suggestion of Ref. [7] (see Ref. [8], for some recent lattice calculations), as well as the strong caveat in Ref. [9]. All that motivates a detailed presentation of our novel method of solving BSE with spin degree of freedom in Minkowski space, which in perspective could give some reliable contribution to the coherent efforts toward an investigation of the 3D tomography of hadrons.
Our method is based on the so-called Nakanishi integral representation (NIR) of the Bethe-Salpeter (BS) amplitude (see, e.g., Ref. [10] for a recent introduction to the issue, and references quoted therein). This representation is given by a suitable integral of the Nakanishi weight function (a real function) divided by a denominator depending upon both the external four-momenta and the integration variables. In this way, one has an explicit expression of the analytic structure of the BS amplitude, and proceeds to formal elaborations. We anticipate that the validity of NIR for obtaining actual solu-tions of the ladder BSE, i.e. the one we have investigated, is achieved a posteriori, once an equivalent generalized eigenvalue problem is shown to admit solutions. Within the NIR approach that will be illustrated in detail for the two-fermion case in what follows, several studies have been carried out. Among them, one has to mention the work devoted to the investigation of: (i) two-scalar bound and zero-energy states, in ladder approximation with a massive exchange [11][12][13][14][15], as well as two-fermion ground states [16,17]; (ii) two scalars interacting via a cross-ladder kernel [18,19]. A major difference, which separates the above mentioned studies in two groups, is the technique to deal with the BSE analytic structure in momentum space. In Refs. [13][14][15]17,19], it has been exploited the light-front (LF) projection, which amounts to eliminate the relative LF time, between the two particles, by integrating over the component k − = k 0 − k 3 of the constituent relative momentum (see Refs. [13,[20][21][22][23] for details). Such an elegant and physically motivated procedure, based on the non-explicitly covariant LF quantum-field theory (see Ref. [24]), perfectly combines with NIR, and, as discussed in the next sections, it allows one to successfully deal with singularities (see Ref. [25] for an early discussion of those singularities) that stem from the spin degrees of freedom acting in the problem. In general, the LF projection is able to exactly transform BSE in Minkowski space into a numerically affordable integral equation for the Nakanishi weight function, without resorting to the so-called Wick rotation [3]. Specifically, the ladder BSE is transformed into a generalized eigenvalue problem, where the Nakanishi weight functions play the role of eigenvectors.
Differently, Refs. [12,16,18] adopted a formal elaboration based on the covariant version of the LF quantum-field theory [26]. This approach has not allowed one to formally identify the singularities that plague BSE with spin degrees of freedom, and therefore, in this case, the eigenvalue equations to be solved are different from the ones we get. In particular, in Ref. [16] the 0 + two-fermion case is studied and a smoothing function is introduced for achieving stable eigenvalues (indeed the eigenvalues are the coupling constants, as shown in what follows). It should be pointed out that in the range of the binding energies explored in Ref. [16], B/m ∈ [0.01 − 0.5] (m is the mass of the constituents), the eigenvalues fully agree with the outcomes of our elaboration for all the three interaction kernels considered (see [17] and the next sections). As to the eigenvectors, the only case discussed in [16] is in overall agreement with our results (cf. Sect. 6).
The aim of this work is to extend our previous investigations of the ladder BSE [10,[13][14][15]17,27] for a 0 + twofermion system, interacting through the exchange of massive scalar, pseudoscalar or vector bosons. Since Ref. [17] was basically devoted to a validation of our method through the comparison of the obtained eigenvalues and the ones found in the literature, in the present paper we first provide the nontrivial details of the formal approach, which can be adopted for future studies of systems with higher spins. Then we illustrate: (i) for each interaction above mentioned, physical motivations for the numerical dependence of the binding energy upon the coupling constant g 2 we have got; and (ii) the peculiar outcomes of the NIR+LF framework, represented by the eigenvectors of our coupled integral system. Let us recall that the eigenvectors are the Nakanishi weight functions, namely the fundamental ingredient for recovering both the full BS amplitude and the LF amplitudes. Furthermore, we extend our analysis to a fermion-antifermion pseudoscalar system with a large binding energy, i.e. in a strongly relativistic regime.
In particular, after tuning both the fermion mass and the mass of the exchanged vector boson to the values suggested by lattice calculations, we show the LF amplitudes for such a mock pion. They feature the effects due to the spin degrees of freedom, and show the peculiarity of an approach addressing BSE directly in Minkowski space. This preliminary study, once it will be enriched with suitable phenomenological features, could be relevant in providing the initial scale for evolving pion transverse-momentum distributions (TMDs) [28].
The present paper is organized as follows. In Sect. 2, we introduce: (1) the basic equation for the homogeneous two-fermion Bethe-Salpeter equation, with a kernel in ladder approximation, based on three different kinds of massive exchanges, i.e. scalar, pseudoscalar and vector bosons; and (2) the general notation for the BS amplitude of a twofermion 0 + state. In Sect. 3, we illustrate the Nakanishi integral representation for the 0 + bound state of two fermions and review the LF-projection technique, which allows one to formally infer an integral equation fulfilled by the Nakanishi weight function. In Sect. 4, we introduce our distinctive method for formally obtaining the kernel of the integral equation fulfilled by the Nakanishi weight function, and we separate out the light-cone non-singular and singular contributions, by carefully analyzing the end-point singularities, related to the spin degree of freedom of the problem we are coping with. In Sect. 5, we provide our numerical tools for solving the integral equation for the Nakanishi weight function, which is formally equivalent to getting a solution of the BSE in Minkowski space. In Sect. 6, we discuss several numerical results for a two-fermion 0 + state: from the dependence of the binding energy upon the coupling constant, being peculiar for the three exchanges we consider, to compute the LF amplitudes, building blocks of both LF distributions and TMDs. In a forthcoming paper [29] we aim to address a phenomenological, but realistic 4D kernel to study TMDs [30]. In Sect. 7, we present an initial investigation of a fermion-antifermion 0 − state, featuring a mock pion, with input parameters inspired by standard lattice calculations. In Sect. 8, we provide a summary and make concluding remarks to close our work.

General formalism for the two-fermion homogeneous BSE
In this Section, the general formalism adopted for obtaining actual solutions of the BSE for a bound system composed of two spin-1/2 constituents is presented. Though the approach based on NIR is quite general, and it can be extended at least to BSE with analytic kernels, given our present knowledge, the ladder approximation is suitable to start our novel investigation on fermionic BSE, since it allows us to cope with some fundamental subtleties without considering additional, but irrelevant for our most urgent aim, complications. We can anticipate that the mentioned issues are related to the singularities onto the light-cone [25], and the efforts for elucidating them are an unavoidable formal step in order to extend the NIR approach to higher spins (e.g. vector constituents). Among the two-fermion bound systems, the simplest one to be addressed is given by a 0 + bound state, which after taking into account the intrinsic parity can be trivially converted into a 0 − fermion-antifermion composite state, once the charge conjugation is applied.
In what follows, we adopt (i) the ladder approximation for the interaction kernels, modeling scalar, pseudoscalar or vector exchanges, and (ii) no self-energy and vertex corrections, apart a scalar form factor to be attached at the interaction vertices (see below). With those assumptions, the fermionantifermion BS amplitude, Φ(k, p) fulfills the following integral equation [16]: where the off-mass-shell constituents have four-momenta given by p 1(2) = p/2 ± k, with p 2 1(2) = m 2 , p = p 1 + p 2 is the total momentum, with M 2 = p 2 the bound-state square mass, and k = ( p 1 − p 2 )/2 the relative four-momentum. The Dirac propagator is given by and the Γ i are the Dirac structures of the interaction vertex we will consider in what follows, namely Γ i ≡ I, γ 5 , γ μ , for scalar, pseudoscalar and vector interactions, respectively. Moreover, using the charge-conjugation 4 × 4 matrix C = iγ 2 γ 0 , we define Γ 2 = C Γ T 2 C, and is a suitable interaction-vertex form factor. Besides the Dirac structure, the interaction vertex contains also a momentum dependence (due to the exchanged-boson propagation) as well as a coupling constant. In particular, depending on the interaction, one has the following expression for iK in Eq. (1): -for the scalar case -for the pseudoscalar one, -and finally for a vector exchange, in the Feynman gauge, The BS amplitude Φ(k, p) can be decomposed as follows [16]: where φ i are suitable scalar functions of (k 2 , p 2 , k · p) with well-defined properties under the exchange k → −k, namely they have to be even for i = 1, 2, 4 and odd for i = 3. The allowed Dirac structures are given by the 4 × 4 matrices S i , viz They satisfy the following orthogonality relation: Multiplying both sides of Eq. (1) by S i and carrying out the traces, there is a reduction to a system of four coupled integral equations, written as with i, j = 1, 2, 3, 4. The coefficients c i j are obtained by performing the convenient traces: and are explicitly given in Appendix A (see also Ref. [16]). Notice that the coefficients for a pseudoscalar exchange and a vector can easily be obtained from the coefficients for a scalar interaction, as explained in Appendix A.

Nakanishi integral representation and the light-front projection
In the 1960s, for a bosonic case, Nakanishi (see Ref. [31] for all the details) proposed and elaborated a new integral representation for perturbative transition amplitudes, relying on the parametric expression of the Feynman diagrams. The key point in his formal approach is the possibility to express a n-leg transition amplitude as a proper folding of a weight function and a denominator that contains all the independent scalar products of the n external four-momenta. Noteworthy, by using NIR, the analytic properties of the transition amplitudes are dictated by the above mentioned denominator. As a final remark, one should recall that the weight function is unique, as demonstrated by Nakanishi exploiting the analyticity of the transition amplitude, expressed through NIR (see [31]). A step forward of topical interest was carried out in Ref. [16], where the generalization to the fermionic ground state was presented. It should be recalled that originally NIR was established only for the bosonic case, with some caveat about a straightforward application to the fermions, as recognized by Nakanishi himself, who was aware of the possible tricky role of the numerator in the Dirac propagator.
In the spirit of Ref. [16], one can apply NIR to each scalar function in Eq. (7), tentatively generalizing the Nakanishi approach to the fermionic case. Let us recall that the denominator of a generic Feynman diagram contributing to the fermionic transition amplitudes has the same expression as in the boson case analyzed by Nakanishi [2], and this is the main feature leading to NIR. In conclusion, one can write for each φ i (k, p) in Eq. (7) where κ 2 = m 2 − M 2 /4. For each scalar function of the BS amplitude it is associated one weight function or Nakanishi amplitude g i (γ , z ; κ 2 ), which is conjectured to be unique and encodes all the non-perturbative dynamical information.
The power of the denominator in Eq. (12) can be chosen as any convenient integer. Actually, the power 3 is adopted following Ref. [16]. The scalar functions φ i (k, p) must have well-defined properties under the exchange k → −k: even for i = 1, 2, 4 and odd for i = 3. Those properties can be straightforwardly translated to the corresponding properties of the Nakanishi weight function 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. By inserting Eq. (12) in Eq. (10), one can write the fermionic BSE as a system of coupled integral equations, given by where the kernel that includes also the Nakanishi denominator of the BS amplitudes on the r.h.s. is It is necessary to stress that the validity of the NIR for the BS amplitude is verified a posteriori. Namely, if the generalized eigen-equation in Eq. (13) admits eigen-solutions then NIR can be certainly applied to the scalar function φ i (k, p). Let us recall that Eq. (10) formally follows from the BSE in Eq.
(1). One can perform the four-dimensional integration on k in Eq. (14), obtaining where P (1) i j (k, p; γ , z ) P (2) i j (k, p; γ , z ) and the coefficients a k i j , b k i j and d k i j are given in Appendix A. In the above equations, one has where the LF variables k ± = k 0 ± k 3 have been used, and the following notation has been adopted: γ = |k ⊥ | 2 and z = −2k + /M ∈ [−1, 1] (see Ref. [10] for details on this range). Moreover, in Eqs. (16), (17) and (18), one has with By direct inspection of Eq. (20), one can realize that several poles affect the needed integration on the relative fourmomentum. In order to treat the poles in a proper way, we adopt the so-called LF projection onto the null-plane (see, e.g. [20,22,23]), which amounts to integrate over the LF component k − both sides of Eq. (13). It should be pointed out that such a LF projection is specific of the non-explicitly covariant version of the LF framework, and it is a key ingredient in our approach, as will be clear in the next sections. As to the present stage of the elaboration, the LF projection marks one of the differences with the explicitly covariant approach in Ref. [16].
To conclude this Section, it should be emphasized that the announced difficulties in the fermionic case are related to the possible light-cone singularities generated by powers of k − contained in the coefficients c i j (k, k , p) (see also Ref. [25] and Refs. [32][33][34] for more details on the light-cone singularities). In the next sections, the singularities will be identified and formally integrated out.

Light-front projection of the fermionic Bethe-Salpeter equation
The integration of the l.h.s. is straightforward, since it is analogous to the scalar case. After the LF projection one gets what we shortly call light-front amplitudes Notice that ψ i are scalar functions rotationally invariant into the transverse plane and they will be used for constructing [29] the socalled valence component of the two-fermion state, once a Fock expansion of this state is introduced. In Eq. (22), i can be removed, since we are dealing with a bound state and κ 2 > 0. Differently, the integration on k − of the r.h.s. of Eq. (13) has to be done very carefully and it represents the core of our approach. A first, but short, presentation can be found in Ref. [17], together with some important numerical outcomes.
Then performing the LF projection as in Eq. (22), one can transform Eq. (13) into a new coupled integral-equation system for the Nakanishi weight functions, viz., where This coupled system has the attractive feature that one has to deal with a single non-compact variable, γ , and a compact one, z (like for γ and z, respectively).
In the next subsection the integration on k − in Eq. (24) will be discussed in detail.

The kernel operator of the coupled integral-equation system
One can write the kernel L i j as follows: where the non-vanishing F 0;i j , F 1;i j , F 2;i j and F 3;i j are listed in Table 1. Moreover, in Eq. (25) the functions C n are integrals over k − defined by with n = 0, 1, 2, 3.
It is fundamental to notice that the powers of k − in the numerator of Eq. (26) are dictated by the coefficients , which in turn are determined by the Dirac structures of both the BS amplitude (cf. Eq. (7)) and the interaction vertices [cf., e.g., the matrix Γ i in the ladder BSE in Eq. (1)]. Since the power of k − in the numerator lead to singular integrals, as discussed below, one has to bring in mind that the most severe singularities can be determined in advance, once the Dirac structures above mentioned are known. Therefore, if one has to deal with the BS amplitude of a fermion-boson system or a vector-vector one, the highest power of k − can easily be singled out after modeling (in agreement to the required symmetries and the interaction) the needed Dirac structures.
As it is clearly shown in Eq. (26), the k − integration of the kernel L i j could become tricky if the powers of k − in the numerator are not suitably balanced by the ones in the denominator for some values of k + D . The values z = ±1, which in principle could generate troubles, are made harmless by the vanishing values of the Nakanishi weight functions at those values (see the discussion of the boundaries for the scattering case in Ref. [15], which can be adapted to the present discussion). In conclusion, to perform the projection of the kernel K i j onto the LF hyper-plane, one has to carefully study the integrals C n . In particular, we should analyze the behavior of the integrand when the contour arc goes to infinity in order to use the Cauchy residue theorem in Eq. (26). The critical situation occurs for k + D = 0, which was not recognized as the source of the instabilities met in Ref. [16] (they were fixed only through a numerical procedure). Indeed, to single out the actual sources of instabilities makes the NIR approach for systems with spin degrees of freedom sounder.

The integration on k − of the kernel L i j and the LF singularities
In this Section, we present the detailed results for the integration on k − in Eq. (26), considering in the next subsections two different cases: (i) k + D = 0, which will lead us to obtain a result corresponding to Eq. (15) of Ref. [16]; and (ii) any k + D , i.e. including both the previous case and the tricky k + D = 0. The analysis of the latter case represents our main contribution to the topic, since we are able both to achieve a formal result, which explains why in Ref. [16] instabilities were found, and to obtain numerical results for the eigenvalues in full agreement with Ref. [16], where the issue was fixed through the introduction of a smoothing function. Also the eigenvalues obtained in Ref. [35] by using a Euclidean BSE are recovered within our approach.

Non-singular kernel
Let us first consider the case k + D = 0. Then, in the denominator of Eq. (26), for n = 0, 1, 2, 3, one has powers of k − sufficiently large for the arc contribution to vanish at infinity and enabling one to safely apply the Cauchy residue theorem. For k + D > 0 one can close the path in the upper semi-plane and pick up the contribution of the residue at the pole k − u +i . It is understood that for z → −1 and k + D > 0 one has a vanishing result, since the remaining poles belong to the same semi-plane. In this case, one gets where the small, positive quantity η has been put in the theta function in order to strictly exclude the case k + D = 0, and D is given in Eq. (21). Moreover, it is noteworthy that the denominator γ + z 2 m 2 + (1 − z 2 )κ 2 is always non-vanishing for a bound state, since κ 2 ≥ 0.
Analogously, for k + D < 0 we can close the path in the lower semi-plane and pick up the residue contribution at k − d − i , obtaining If the condition k + D = 0 is satisfied, by defining one can get what we call non-singular contribution to the kernel L i j i.e.

4.2
The k − integration of the kernel L i j for any k + D For vanishing k + D the powers of k − in the numerator of L i j make the analytic integration tricky, since one cannot naively close the integration path with an arc at infinity. In Appendix B, all the needed details are given, while in this section the results are summarized and commented.
The main formal tool we have exploited is given by the following singular integral studied in detail by Yan in Ref. [25], which, indeed, belongs to a series of papers devoted to the analysis of the field theory in the infinite momentum frame. The mentioned integral is In order to complete our analysis we also exploited the derivative with respect to y, i.e.
with n > 2. Indeed, in what follows, only first (n = 3) and second (n = 4) derivatives are used.
Then the function C n can be decomposed as follows: with C (ns) n given in Eq. (29) and (see Appendix B) According to the above decomposition, the kernel L i j in Eq. (24) can be written where L (ns) i j is given in Eq. (30) and L Due to the factor v 2 (1 − v) 2 in Eq. (39), which helps us to get rid of the end-point issues at v = 0 and v = 1, one can safely write obtaining with D (s) The following physical interpretation should be associated with the presence of singular contributions to the kernel of the integral equation when k + D vanishes, since it is proportional to the plus momentum carried by the exchanged quantum. Note that k + D ∝ (z − z), where the momentum fractions of the external and internal fermions are 0 < ξ = (1−z)/2 < 1 and 0 < ξ = (1 − z )/2 < 1, respectively. Hence, one can identify the contribution of a LF zero-mode to the kernel when the exchanged quanta propagates along the light-like direction (k + D = 0). In some cases it can give rise to a LF singularity, as mentioned above.
Of course, the problem is avoided if the fermions propagate before the exchange of the light-like quanta. The high virtuality of the quanta with k + D = 0 in general suppresses any intermediate propagation, as for example in the twoscalar BSE, unless the exchanged boson couples with higherspin particles. This happens when the so-called instantaneous terms and/or derivative coupling are considered.
In conclusion, the final form of the homogeneous ladder BSE for two fermions reduces to where L (ns) i j (γ , z, γ , z ) and L (s) i j are given by Eqs. (30) and (39), respectively.
In physical terms the essential contribution of the instantaneous terms and the light-like exchanged quanta should be associated with a two-fermion probability amplitude equally distributed along the light-like trajectory contained in the null-plane.

Numerical method
For the 0 + two-fermion bound state, the system in Eq. (43) contains four coupled integral equations. In order to solve such a system, we introduce a basis expansion for the Nakanishi weight function where (i) w i mn are suitable coefficients to be determined by solving the generalized eigenproblem given by the coupled system of integral equations (43), (ii) G λ i 2m+r i (z), with r i = 0 for i = 1, 2, 4 and r 3 = 1 (due to the symmetry under the exchange z → −z) are related to the Gegenbauer polynomials C λ i 2m+r i (z) and (iii) J n (γ ) to the Laguerre L n (aγ ). In particular, one has with q = (2λ − 1)/4. The following orthonormality conditions are fulfilled: Notice that the kernel L (s) in Eq. (43) contains terms proportional to derivatives of the delta-function δ (z − z) and therefore the numerical evaluation needs derivatives of the basis functions G λ i (z). After inserting the above expansion in Eq. (43), one obtains a matrix form of the system, reducing to solve a discrete generalized eigenvalue problem as in the simple case of two scalars (see ref. [13]). At the end, the generalized eigenvalue problem to be solved looks like with (i) C and D square matrices, and (ii) w the eigenvector that allows one to reconstruct the Nakanishi weight functions (cf. Eq. (44)). The eigenvalue is g 2 , as in the standard way for investigating the ladder BSE, since the binding energy, defined by B = 2m − M is a non-linear parameter, and it is assigned in the range 2 ≥ B/m ≥ 0, given the well-known instability that occurs for odd numbers of boson fields (see, e.g., [36] for the φ 3 case).
In Ref. [17], where the eigenvalues were shown for several cases, the calculations have been carried out taking as the biggest basis the one with 44 Laguerre polynomials and 44 Gegenbauer ones with indices for each Nakanishi weight functions 5/2, 7/2, 7/2, 7/2, respectively. To improve the convergence, in those calculations the parameter a = 6.0 has been adopted, and the variable γ has been rescaled according to γ → 2γ /a 0 with a 0 = 12. The two parameters a and a 0 help to take into account the range of relevance of the Laguerre polynomials and the structure of the kernel, respectively. Finally, a small quantity, = 10 −7 , has been added to the diagonal elements of the discretized matrix on the l.h.s. of the system in Eq. (43).
The study of the eigenvectors is more delicate, even if the final goal is the calculation of the LF amplitudes, Eq.

Results for scalar, pseudoscalar and vector exchanges
The solution of the coupled set of integral equations for the Nakanishi weight functions (43) was obtained for three different ladder kernels (see Sect. 2) featuring scalar, pseudoscalar and vector exchanges, besides an interaction vertex smoothed through a form factor [Eq. (3)]. Our study aims at singling out signatures of the dynamics generated by the different kind of exchanges. Indeed, among the physical motivations of such a broad analysis, we have to put the application to the study of pseudoscalar mesons, originating from vector exchanges between quarks, and in what follows a first attempt, which we call a mock pion, will be discussed (see Sect. 7).
The differences between the coefficients present in the kernel of the coupled equations are due to the peculiar Dirac structures entailed by scalar, pseudoscalar and vector exchanges [see Eqs. (A.3), (A.4)]. Such differences are naturally reflected in the corresponding Nakanishi functions, since they are suitably weighted in the integral equations through the above mentioned coefficients. Eventually, the differences show up in the LF amplitudes, Eq. (22).
In this section we focus on the scalar, pseudoscalar and vector cases, leaving the actual application to the mock pion to the next section.

Binding energy vs. coupling constants
We start by showing the binding energy B/m as a function of g 2 for both scalar and pseudoscalar exchanges, fixing the cutoff in the vertex form factor, Eq. (3), to the value Λ/m = 2. The results shown in Fig. 1, partially presented in Table I and II of Ref. [17], have been obtained by choosing the masses  Fig. 1 allow one to illustrate interesting overall behaviors of the binding B/m, when the boson mass μ changes. In general, for fixed values of Λ/m and B/m, the coupling constant g 2 increases as the mass of the exchanged boson increases. This feature is an expected one, if we take into account that the interaction has a range controlled by the inverse of μ. Hence, if the range shrinks, then an increasing of the interaction strength is necessary for allocating a bound state with given B/m, as is well known also in the nonrelativistic case. In particular, for weakly bound systems, at fixed B/m, the correlation between the coupling constant g 2 and μ/m is almost linear (cf. Tables I and II in Ref. [17] for B/m ≤ 0.05). Moreover, the correlation between g 2 and μ is also driven by the nature of the exchanged boson, as shown by the quantitative differences illustrated in the upper and the lower panels in Fig. 1. As is well known in the non-relativistic framework when the one-pion exchange is investigated, the pseudoscalar interaction is weak in a 0 + state. In particular, the tensor force obtained from the non-relativistic reduction does not contribute in a 0 + state, and the remaining interaction gives a repulsive contribution, making impossible to bind the system. Only when the relativistic effects become important, the binding can be established, but for very large g 2 . In conclusion, in the pseudoscalar case, one is confronted with a much weaker attraction that requires, for a given binding energy, larger couplings than the scalar case needs.
When we approach the collapse of the bound state, i.e. when a composite massless particle is created since B/m → 2, the shape of the binding becomes less sensitive to the variation of g 2 (a vertical asymptote is approached). By increasing the binding, the size of the system becomes smaller and smaller and therefore the range of the interaction, dictated by μ/m, becomes less relevant for the functional dependence of the binding upon g 2 close to its critical value for B/m = 2, namely the ultraviolet regime is now starting to govern the dynamics inside the system. This holds for both exchanges.
Let us now focus on the behavior of d(B/m)/dg 2 , shown in Fig. 2. It is quite peculiar for the scalar case with respect to both pseudoscalar and vector cases (for the figure illustrating B/m for the massless vector exchange see Ref. [17]), since it shows the presence of a minimum (less pronounced for μ/m = 0.5) positioned almost in the middle of the range of g 2 relevant for a given μ. For the pseudoscalar case the minimum is close to the threshold, as a consequence of the repulsive contributions at low binding energies above mentioned. This observation suggests that also for the scalar case the appearance of a minimum be related to some repulsive contributions that develop for B/m → 2 (i.e. M → 0). As a matter of fact, from direct inspection of Table 1, one can single out repulsive contributions that depend upon 1/M. The LF amplitudes ψ i (γ , ξ ; κ 2 ), Eq. (22), for the 0 + twofermion system bound by a scalar exchange are presented in the Fig. 4. We compare two cases: weak binding B/m = 0.1 (upper panels) and strong binding B/m = 1.0 (lower panels). The motivation of such a comparison is given by the attempt of widening our intuition, based on non-relativistic physics, to the extreme binding, where the relativistic effects are expected to be relevant, like in the case of light pseudoscalar mesons. Indeed, in the present work an initial analysis of this physical case will be proposed by introducing and investigating a mock pion (see Sect. 7). To begin, it is useful to compare our results for the Nakanishi weight functions with the outcomes obtained in Ref. [16], within the Nakanishi framework, but inserting a smoothing function. In Fig. 3, the amplitudes g i (γ , z; κ 2 ), corresponding to μ/m = 0.5 and B/m = 0.1 (weak binding) as in Ref. [16], are presented. The upper panel shows g i (γ , z; κ 2 ) for a fixed values z 0 = 0.6 and running γ , while the lower panel illustrates the same quantities, but for γ 0 /m 2 = 0.54 and running z. Recall that the results for the eigenvalues obtained in Ref. [16] coincide with ours (cf. [17]) at the level of the published digits, while the eigenvectors, i.e. the Nakanishi weight functions, are slightly different. Indeed, it is extremely encouraging to observe that so very different methods for taking into account the singularities in the BSE are able to achieve an overall agreement, (at least in weakbinding regime) in describing g i (γ , z; κ 2 ), that have very rich structures. In Fig. 4, results obtained with Λ/m = 2 and μ/m = 0.15 are presented. It should be pointed out that the general behavior shown in Fig. 4   All the LF amplitudes are normalized to ψ 1 (γ = 0, ξ = 0.5), to quickly appreciate the relative strengths (in a forthcoming paper [29], it will be adopted the proper normalization through the BS amplitude for eventually obtaining the LF distributions).
The behavior for running γ shows a difference between the weak and strong binding that can be directly ascribed to the size shrinking of the systems when the binding increases. Therefore the overall growing pattern of the amplitudes at higher γ /m 2 when B/m increases (cf. the l.h.s. of Fig. 4 upper and lower panels) is an expected one. Remarkably, the characteristic momentum-width one can infer from the lines in Fig. 4 is of the order of γ /m 2 ∼ B/m. Another interesting feature is the growth of the amplitude ψ 2 when B/m increases. This reflects the importance of singular terms in the ultraviolet region, since in this case one has the dominant role of the coefficient c S 23 [see Eqs. (A.1) and (A.2) for counting the involved power of k]. Already in Ref. [16] it was noted the strong singular contribution of the kernel at the end points when c 23 is nonzero, which is enhanced in the relativistic case of strong binding.
The dependence on ξ of the LF amplitudes ψ i presented in r.h.s. of Fig. 4 shows the expected peaks around 1/2, except for ψ 3 , which is antisymmetric in z. Less trivial is the comparison between the weak-and strong-binding regime, in particular for ψ 2 and ψ 4 . It can be seen that, for B/m = 1, the LF amplitude ψ 2 accumulates toward the end points, a property which seems to be slightly present also in ψ 1 . This feature can be traced to the same effect seen for running γ , namely the relevance of the coefficient c S 23 , associated to the singular behavior at the end points. For the weak-binding case, the effect is damped, as one could naively guess from the observation that be found for both the pseudoscalar and the vector exchange (see below). In particular, this last case leads us to focus the previous discussion on the coefficients by excluding c S 44 , since c V 44 = 0 but the effect is present. Moreover, given the smallness of ψ 3 one can pay attention only on c S 14 and c S 24 . As a final remark, one should notice in the r.h.s. that the LF amplitudes enlarge their-own range when B/m grows, as expected from the size shrinking of the system.

Light-front amplitudes: pseudoscalar case
The general comments, with a particular emphasis on the relevance of the coefficients producing LF singularities, presented for the scalar case are suitable also for the LF amplitudes ψ i (γ , ξ ; κ 2 ) obtained by using a pseudoscalar exchange with μ/m = 0.15. But a crucial difference is generated by the peculiar dynamics entailed by the spinor coupling characterizing the two cases, as already pointed out for Fig. 1. Such a difference plays a role in explaining the leading position of ψ 4 , which is related to the spin-momentum correlations.
The dependence of ψ i on ξ is presented in the right panels of Fig. 5. For the weakly bound case, one notices that ψ 1 and ψ 2 are somewhat different while in the scalar case they almost coincide. We trace back the reason for that by looking at the spinor structure associated to each ψ i . The Dirac operator for ψ 1 is γ 5 and for ψ 2 is / pγ 5 /M. For both cases in the weakbinding limit, namely in the non-relativistic regime, the vector charge and scalar charge densities of the state are expected to be close, and in the scalar coupling case the corresponding operators commutes with the scalar coupling Dirac operator. For the pseudoscalar case, the γ 5 coupling has opposite commutation properties with the Dirac structure associated to ψ 1 and ψ 2 , and this feature can be identified as the source for the observed difference in the upper-left panel of Fig. 5.
In conclusion, the performed analysis of scalar and pseudoscalar exchanges yields a description coherent with the intuition stemming from the general structure of both the BS amplitude and the kernel, and increases the confidence in the novel approach we are pursuing for solving BSE in Minkowski space. Therefore, to address the vector exchange with such a reliable tool becomes a very stimulating issue, given its possible application to hadron physics. Starting from the solution of the four-dimensional BSE, our goal is to construct a phenomenological tool for analyzing the structure of a simple model of the pion, namely a quarkantiquark pair living in Minkowski space and bound through a massive vector exchange. Before considering a quark-antiquark system, we continue our analysis of the changes in the structure of the 0 + state with different couplings and binding energies in order to fully appreciate the subtleties generated by the dynamics inside the bound state. It is worth bearing in mind how accurate are our calculations by presenting a quantitative comparison with analogous calculations shown in Ref. [16]. The massless exchange between two fermions in the 0 + state was considered [i.e. Eq. (6) with μ 2 = 0]. As shown in Table 2, the agreement is excellent, and, moreover, the possibility to have a formal treatment of the LF singularities allows us to extend the range of the calculations. The behavior of B/m as a function of g 2 (see the corresponding figure in [17]) has the same overall structure found for the pseudoscalar case, with a minimum of the derivative of d(B/m)/dg 2 close to the threshold.
Given that, before treating the mock pion problem, we present results for a massive vector with μ/m = 0.15, still using the form-factor parameter equal to Λ/m = 2. As it is easily seen, in Fig. 6 general patterns similar to the ones observed for the pseudoscalar case can be recognized, when one moves from weak-to strong-binding regime.
The dependence of the LF amplitudes ψ 1 and ψ 2 on ξ , in the right panels of the figures, suffers a dramatic change from weak to strong binding, while the amplitudes ψ 3 and ψ 4 just get wider when the binding increases (apart the change of sign of ψ 3 ). As in the scalar case, the amplitudes ψ 1 and ψ 2 almost coincide for running ξ and γ = 0, in the weak-binding case. This feature can be understood by considering the Dirac structure associated with ψ 1 and ψ 2 in the expansion of the BS amplitude, i.e. γ 5 for the first amplitude and γ 0 γ 5 for the second one, in the CM frame. In the non-relativistic limit, the vector interaction largely reduces to its Coulomb component, whose Dirac structure, γ 0 , has equal commutation properties with the ones associated to ψ 1 and ψ 2 . The same happens when the scalar exchange is considered, since one has simply the identity matrix at the interaction vertex. Moreover, in the weak-binding regime, the scalar and charge densities of the fermionic constituents (acting at the interaction vertices) tend to be the same, while for growing B/m this is not the case. Finally, notice that even the coupling constants g 2 are similar for scalar and vector exchange, in the weak-binding regime.
Moving from weak to strong binding the amplitudes ψ 1 , ψ 3 and ψ 4 become wider as a function of ξ . A new feature arises in ψ 1 , it becomes quite flat and sharply decreasing at the end points, while ψ 2 has the characteristic peaks discussed before. The flattening of ψ 1 appears to be a signature of the vector coupling, since it remains present when the vector mass increases, as illustrated in the next section. Indeed, all the patterns shown in Fig. 6 do not qualitatively change when the mass of the exchanged boson grows.

The mock pion
In this section, we present a first investigation of a simplified model for the pion, which, to some extent, can be considered Lattice-QCD inspired. In particular, we have tuned our parameters, like the masses of quark and exchanged vector boson, according to the outcomes of Lattice-QCD calculations.
As above illustrated, within the formal elaboration one adopts for getting solutions of BSE, it is natural to switch to a BS amplitude pertaining to a fermion-antifermion system. If one multiplies the original fermion-fermion BS amplitude by the charge-conjugation operator, then one gets a compact expression of the BSE with fermionic degrees of freedom, eventually using the standard multiplication rule for matrices. Summarizing, the scalar functions [Eq. (7)] that enter both fermion-fermion (0 + ) and fermion-antifermion (0 − ) BS amplitudes are the same, but the Dirac structures are different, as it must be. The ladder kernel we exploit is the one corresponding to a massive vector exchange, in Feynman gauge [cf. Eq. (6)], with the same vertex form factor in Eq. (3).
Lattice-QCD results obtained in Refs [37][38][39], by adopting the Landau gauge, indicate that the gluon propagator obtained in the space-like region is infrared finite and can be roughly described by using a massive propagator with a dynamical gluon mass μ ∼ 500 MeV. For the moment, we assume such a handy approximation for the gluon propa- Finally, it is worth noticing that the assumed value of the constituent mass gives a binding energy B = 1.44 m = 360 MeV, which leads to a typical size ofh c/B = 0.55 fm, quite close to the experimental pion charge radius of 0.672 ± 0.008 fm [42]. The eigenvalues g 2 , obtained with the above set of parameters, and the corresponding α s are summarized in Table 3, while the LF amplitudes of our mock pion, are presented in Fig. 7. It is possible to recognize an overall behavior of the ψ i 's similar to the one in Fig. 6, where a strong-binding case B/m = 1 for a vector exchange with μ/m = 0.15 and Λ/m = 2 is shown. But it is much more interesting to observe the effects produced by varying exchanged mass μ/m and vertex cutoff Λ/m. As μ/m increases the tails of the LF amplitudes do the same (cf. l.h.s. of Figs. 6 and 7): given the overall size of the system (fixed by the binding energy), one has a reduction of the interaction range, which in turn triggers a shrinking of the system itself. The same kind of effect can be seen also decreasing the size of the interaction vertex, i.e. by increasing Λ/m. More striking is the effect shown in the r.h.s. of Fig. 7, where the LF amplitudes are given as a function of the variable ξ . There, one can observe a quite peculiar two-horn pattern, which becomes more enhanced when Λ/m grows. The striking two-horn patterns are due to the presence of spin degrees of freedom (see, the differences with the case of a two-scalar system in Ref. [13]). Bearing in mind that we will investigate such structures in a forthcoming paper [29], we expect non-trivial impacts of the above feature on the evaluation of both TMDs (see e.g. [43]) and valence distributions (which contain a multiplicative factor ξ(1 − ξ)). In view of this, it is suggestive to recall the wellknown pion distribution amplitude introduced in Ref. [44], which displays a two-horn pattern.

Conclusions and summary
In the present paper, we have described in great detail our approach to getting actual solutions, directly in Minkowski space, of the ladder Bethe-Salpeter equation for a system of two interacting fermions, in a state 0 + , as well as a fermionantifermion system in 0 − channel. This effort makes complete the presentation of our approach, shortly illustrated in Ref. [17], where the most urgent aim was the quantitative comparison with the eigenvalues obtained in Ref. [16] and in the Euclidean space [35]. The large amount of information provided allows one to assess in depth the approach based on the Nakanishi integral representation of the BS amplitude. Hence, one can easily recognize that the approach is a very effective tool for exploring the non-perturbative dynamics inside a relativistic system, since it combines flexibility (one is able to address higher-spin system) and feasibility of affordable calculations. The main ingredients for proceeding toward the numerical solution of the coupled system of integral equations, formally equivalent to the initial BSE, are (i) the Nakanishi integral representation of the BS amplitude and (ii) the LF projection of the BS amplitude, which allows one to get LF amplitudes depending upon real variables. In order to embed the approach based on NIR into a well-defined mathematical framework, it is fundamental to notice that those LF distributions have exactly the form of a general Stieltjes transform [45], and therefore the NIR approach appears to be more general than one could suspect from the Nakanishi work focused on the diagrammatic analysis of the transition amplitudes [31,46]. The analysis of the Stieltjes transform method for investigating the fermionic BSE will be left to future work. Moreover, the LF projection plays also a keyrole in determining the LF singularities that, in principle, plague the method (as shown in Ref. [16], where a function was introduced for smoothing the singular behavior of the involved kernel). In particular, by exploiting previous analyses of the LF singularities [25], we succeeded in formally treating them. As a consequence, we have obtained calculable matrix elements, after introducing an orthonormal basis for expanding the four Nakanishi weight functions needed for achieving the description of the fermionic states (for a 1 + states the weights become eight).
We have shown and discussed both eigenvalues and eigenvectors of the generalized eigenvalue problem one gets after inserting in the BSE three different kernels, featuring massive scalar, pseudoscalar and vector exchanges. Such interactions produce very peculiar form for the correlation between the binding energy and the coupling constant g 2 . From the eigenvectors, i.e. the Nakanishi weight functions, we have evaluated the corresponding LF amplitudes, which show the effects of the exchange in action. It should be recalled that the LF amplitudes are the building blocks for constructing transverse-momentum distributions and valence wave functions (see, e.g. [10,13]). Last but not least, we have presented the application to a simplified model of the pion, by taking the mass of the constituents and the mass of the exchanged vector boson (with a Feynman-gauge propagator) from some typical lattice calculations. The intriguing feature shown in the r.h.s. of Fig. 7 will be further discussed in Ref. [29], where also the transverse-momentum distributions of a 0 − state will be investigated.
In conclusion, we would like to emphasize that the present approach can be certainly enriched with new features impacting both kernel and self-energies of both constituents and intermediate boson, but already at the present exploratory stage it appears to have a great potentiality, as a phenomenological tool for facing with strongly relativistic systems, in Minkowski space.
with c S 13 = c S 31 = 0, B and B defined by The other two sets of coefficients can be obtained from c S i j by exploiting the properties of the Dirac matrices. In particular, one gets for the pseudoscalar exchange and for the vector exchange For carrying out the discussion of the LF singularities one meets, it is useful to decompose the coefficients c i j (k, k , p) so that terms with the same power of k μ are gathered. Namely where the non-vanishing coefficients a n i j , b n i j and d n i j are given by  25), for any value of k + D . Namely, we obtain the general expression of the coefficients C n defined in Eq. (26). They are necessary for determining the kernel L i j . One can write where n = 0, 1, 2, 3 and where the Feynman trick has been applied to the denominators containing k + D for obtaining the final expression. The issue is the possibility or impossibility to apply the Cauchy residue theorem for getting an analytic result. Hence, one has to carefully check if the arc at infinity gives or not a vanishing contribution. As is shown in what follows, a positive answer depends upon the values of n and k + D . The n = 0 case is not affected by any difficulty due to values of k + D . One can always close the arc at infinity, even for k + D = 0. Indeed, if the last case happens, one could be concerned about the end points z = ±1, but, fortunately, one is left with a standard LF integration and gets a finite result (see Ref. [33]). The evaluation of B 0 proceeds by exploiting the residue theorem, leading to , and we have an analogous expression for I(k − d ). Then C 0 is given by For n = 1, the integral B 1 is . (A.14) Then C 1 is given by (A.15) For n = 2, the integral B 2 is , (A.16) where Then C 2 is given by the sum of a singular term and the non-singular one already obtained, i.e.
For n = 3, the integral B 3 is Then C 3 contains both singular and non-singular contributions as in the case of C 2 . One has (A.20)

Appendix C: Singular terms for scalar, pseudoscalar and vector interactions
In this appendix, the final expressions of the singular contributions obtained in Appendix B are explicitly given, in order to facilitate a quick usage of the main results of our work. For the scalar case, the non-vanishing contributions for the singular part of L (A.23) In the above expressions Notice that even for Λ 2 < μ 2 the factor D does not vanish, due to the factor (1 − v)μ 2 in˜ D , which is canceled by a corresponding factor in the square bracket of D .