Bethe-Salpeter bound-state structure in Minkowski space

The quantitative investigation of the scalar Bethe-Salpeter equation in Minkowski space, within the ladder-approximation framework, is extended to include the excited states. This study has been carried out for an interacting system composed by two massive bosons exchanging a massive scalar, by adopting (i) the Nakanishi integral representation of the Bethe-Salpeter amplitude, and (ii) the formally exact projection onto the null plane. Our analysis, on one hand, confirms the reliability of the method already applied to the ground state and, on the other one, extends the investigation from the valence distribution in momentum space to the corresponding quantity in the impact-parameter space, pointing out some relevant features, like (i) the equivalence between Minkowski and Euclidean transverse-momentum amplitudes, and (ii) the leading exponential fall-off of the valence wave function in the impact-parameter space.


Introduction
In the last decade and a half, a quite effective tool is emerging for solving the Bethe-Salpeter equation (BSE) [1] directly in Minkowski space, i.e. avoiding to look for solutions in the Euclidean space by exploiting the Wick rotation [2]. The novel approach is based on the Nakanishi integral representation (NIR) of the n-leg transition amplitudes, that was proposed long time ago [3,4,5]. In particular, the Bethe-Salpeter (BS) amplitude can be formally written like the NIR for the 3-leg amplitude, namely a proper folding of an unknown Nakanishi weight function and a denominator that contains the analytic structure [6,7,8,9,10,11,12,13]. To be practical, let us quickly mention the main features of the NIR for the BS amplitude: (i) the weight function is real and smooth for bound states (see [14] for the zero-energy scattering states), and (ii) it depends upon real variables, of which one is non-compact and the others are compact; (iii) the denominator must depend only upon the independent scalars that can be constructed from the external momenta. Assuming the validity of the NIR for the bound-state case, and taking advantage of the above mentioned features, one can exactly project onto the null-plane (see, e.g. [15,16]) the BS amplitude integrating over the Light-front (LF) variable k − = k 0 − k 3 (k + = k 0 +k 3 and k ⊥ ≡ {k 1 , k 2 }) and formally obtain the so-called LF valence wave function (cf. [17,18,12]), i.e. the amplitude of the component with the lowest number of constituents when the LF Fock expansion of the interactingsystem state is considered. Remarkably, within the NIR approach, the LF valence wave function is given by a non-singular integral involving the Nakanishi weight function. This suggests to integrate on k − both sides of the BSE, getting an integral equation for the Nakanishi weight function. If there exist solutions for this integral equation (this validates a posteriori the previous assumption), then the BS amplitudes of bound states can be reconstructed. In particular, when the above procedure is applied to the BSE with the irreducible kernel in ladder approximation, a generalized eigen-equation for the Nakanishi weight function is obtained (see, e.g., Refs. [9,13] for the LF case), while for the cross-ladder case one has to deal with a non-linear eigen-problem (see Ref. [10]). Notice that in Refs. [6,7] solutions of the scalar BSE in ladder approximation have been obtained by using (i) standard variables (and not the LF ones), and (ii) exploiting the uniqueness of the Nakanishi weight function.
Aim of the present work is to carefully study both spectrum and 3D structure of the bound states, obtained by solving the ladder BSE for a system composed by two massive scalars interacting through a massive scalar. Such an investigation is a natural extension of the previous analysis of only the ground state [13]. In particular, the structure is studied by means of the 3D representation of the LF valence component, both in momentum and impactparameter (IP) spaces. One of the motivations for starting a detailed analysis of the non-perturbative features of an interacting system in momentum and IP spaces (see, e.g. [19] for an introduction) is given by the increasing interest on this topic in hadronic physics, where the valence component plays an important role in determining the dynamical properties of hadrons. For instance, the valence component is an important dynamical ingredient for evaluating parton transverse-momentum distributions, which depend upon both the Bjorken momentum fraction x and the transverse components of parton momentum [20,21], or parton density distributions in IP space, that can be related to the generalized parton distributions (see, e.g., Ref. [20]).
Our paper is organized as follows. In Sec. 2, we quickly introduce the general formalism (see, e.g., Refs. [12,13,14] for more details) and we present a comparison between Minkowski and Euclidean results for the eigenvalues of the relevant integral equation. In Sec. 3, the valence LF wave function and the corresponding density distributions, evaluated both in transverse-momentum space and impact-parameter one, are discussed, showing our numerical results for the available spectrum together with some interesting formal outcomes of our analysis. In Sec. 4, conclusions are drawn and some perspectives presented.

Minkowski space solutions of the Bethe-Salpeter equation
Let us recall the general formalism we have adopted to solve the BSE in Minkowski space. As it is well known the BSE in momentum space for a relativistic bound state is given by the following homogeneous integral equation where i K(k, k ; p) is the interaction kernel that contains all two-body irreducible diagrams, p µ is the total momentum with the bound state mass given by M = √ p 2 . In the present approach we do not consider the self-energy contribution, so that G (12) 0 (k, p) is the product of two free propagators, with m the constituent mass. The BS amplitude for an s−wave state solution of Eq. (1) can be written in terms of NIR as [9,12,13] Φ(k, p) where κ 2 ≡ m 2 − M 2 /4. By substituting (3) into (1) and integrating over k − on both sides, one can obtain the following generalized integral equation for the Nakanishi weight function (for details see Refs. [9,12,13]): with γ = |k ⊥ | 2 , ξ = (1 − z)/2 and ψ(ξ, k ⊥ ) is the valence light-front wave function (the factor √ 2 comes from the symmetry of the problem; see, for example, [12]). In Eq. (4) V LF is the Nakanishi kernel given in terms of the BS 4D kernel, by In this work we adopt the ladder approximation for the BS kernel: where α = g 2 /(16πm 2 ) and µ is the exchanged-scalar mass. According to [13], we have solved Eq. (4) by using a basis function expansion of the Nakanishi weight function, composed by Laguerre polynomials L j (aγ) (with j = 0, 1, N g ) for describing the γ-dependence (where a is an appropriate parameter, as discussed in [13]) and even Gegenbauer polynomials C (5/2) 2 (z) for the z one (with 2 = 0, 2, ..., 2N z ). More specifically, for the γ-dependence we use an expansion in terms of the functions L j (γ) ≡ √ aL j (aγ)e −aγ/2 , where ∞ 0 dγL i (γ)L j (γ) = δ ij . The expansion in Gegenbauer polynomials is given in terms of the functions G (z) ≡ 3 (2 )! 2 + 5 2 /Γ(2 This last choice is dictated by the symmetry property of the Nakanishi weight function g(γ, z; κ 2 ) = g(γ, −z; κ 2 ), that is requested by the bosonic nature of the adopted constituents. It should be recalled that a definite statistical property of the BS amplitude avoids the socalled abnormal solutions of BSE, namely the ones with negative norm [3,4,7], that are associated with excitations in relative time of the bound states (see Refs. [22,23] for a more recent discussion of the issue). Finally, the z 2 dependence of g(γ, z; κ 2 ) entails a symmetry of the valence wave function, namely In our numerical approach, accurate convergence was achieved for the ground state by using 14 Laguerre (N g = 13) and 10 Gegenbauer (N z = 9) polynomials. For the excited states, the convergence was reached with 26 Laguerre and 10 Gegenbauer polynomials. After introducing the basis function expansion and the ladder approximation Eq. (7), Eq. (4) turns into the matrix form of a generalized eigenvalue problem. In particular, one can symbolically write where g is the eigenvector. Differently from the familiar non-relativistic case, in the eigen-equation (8) the role of eigenvalue is played by the coupling constant α, while the mass of the system M is a parameter that can be assigned, after fixing the exchanged-scalar mass µ. In the standard way of analyzing the BSE within the NIR framework [6,7,8,9,10,11,12,13], one introduces the binding energy as which constrains the range of B/m to the interval between 0 and 2, i.e. 0 ≤ M/m ≤ 2, avoiding in this way the well-known instability of the φ 3 model (see Ref. [24]). In our NIR studies of the ladder BSE, Eq. (8) has a pivotal role. First, after fixing m, µ and B gr , it yields the coupling constant of the ground state, i.e. the smallest value of the coupling constant that we call α gr ; secondly it allows one to calculate the spectrum. Indeed, once we have found the coupling constant α gr , one can find the excited state with respect to B gr by slightly changing Eq. (8), as follows In other words, after fixing m, µ and α gr , we search for values M = 2m − B > M gr that produce eigenvalues λ = 1 (as trivially seen, for M = M gr one has λ = 1).

Comparing Minkowski and Euclidean eigenvalues
In order to check the reliability of the computed masses for the excited states, we provide a comparison between the results of our calculations, obtained in the Minkowski space within the NIR, with those one can evaluate in the Euclidean space. In Table 1, we show the binding energies, in unit of the constituent mass m, for the first, B(1)/m, and the second, B(2)/m, excited states, corresponding to a ground state B(0)/m = 1.9 and different values of µ. The choice of such a large binding energy is motivated by the fact that strongly-bound states should be affected by large relativistic effects. First, we have verified that the values of α gr for the binding energy of the ground state of B(0)/m = 1.9, obtained with Euclidean-and Minkowski-space calculations are the same within our numerical accuracy. Then, we have computed the excited state energies given in Table 1, achieving a very satisfactory agreement between the results evaluated in the two spaces. As a remark on the numerical procedure, it should be pointed out that for values of µ/m smaller than 0.05 or B/m < 0.01 the convergence is quite slow, and it is needed an extrapolation of the results with respect to N g , in order to accurately determine the eigenvalues.
It is also important to show the behavior of energy ratios, B(n)/B(0) with n ≥ 0, for small B/m and µ → 0. For a bound state composed for two spinless bosons exchanging a massless scalar boson, the corresponding relativistic expression in lower orders of α, as derived in [25], is The first term is the non-relativistic limit. As verified in Fig. 1, where B(1)/B(0) and B(2)/B(0) are shown for small values of µ/m, the energy ratios are consistent with the non-relativistic limit. Moreover, the agreement for µ → 0 between the relativistic and non-relativistic eigenvalues is observed only for small values of n. Indeed, as binding energies increase relativistic effects become larger and larger.

Valence light-front wave function and momentum distributions
It is attractive to perform numerical comparisons of dynamical quantities that in perspective could be useful for an experimental investigation of actual interacting systems. In view of this, from the valence LF wave function introduced in Eq. (5) (see the next subsection for the numerical results), one can define both the probability distribution to find a constituent with LF longitudinal Table 1 Comparison of the spectra obtained in the Euclidean space in the Minkowski one, by varying µ/m and, consequently, α gr , but taken fixed the value of the ground-state binding energy to B(0)/m = 1.9. fraction ξ = p + i /P + , given by and the probability distribution to find a constituent with LF transverse momentum k ⊥ = |k ⊥ |, that reads Both LF distributions are normalized to the probability of the valence component, once the BS amplitude itself is properly normalized (see Ref. [26] for a general discussion and Ref. [13] for the application within the NIR). Such a probability yields the probability to find the valence contribution in the LF Fock expansion of the interacting two-scalar state (see, e.g., [18,12,27]). As a matter of fact, one has Notice that k ⊥ can be associated with the intrinsic transverse momentum, in the frame where p ⊥ = 0, which is allowed by the covariance of our description.
Although we have discussed the issue of the proper normalization of the valence state, in what follows we are interested on the overall 3D structure of the valence wave function, and therefore we have simply adopted an arbitrary normalization.

Momentum space valence wave function for excited states
In Figs. 2 and 3, we present the LF wave function of the first (left panels) and second (right panels) excited states, corresponding to the case µ/m=0.1, α gr =6.437, B(1)/m=0. 22 and B(2)/m=0.05 (see Table 1). As clearly shown, the wave function displays the typical feature of the first and second excited states, i.e. one and two nodes, respectively.
By a direct inspection of the corresponding panels for the first excited state in Fig. 2 and for the second excited state in Fig. 3, one observes that, in the plane (ξ , k ⊥ /m), the node structure is present for (k ⊥ /m) 2 < 1 and ξ < 0.75, and it is symmetric with respect to ξ = 1/2. In particular, the node structure moves toward ξ = 1/2 as k ⊥ increases. Such a behavior can be naively expected when Cartesian three-momenta are adopted. As a matter of fact, the relation between Cartesian and LF components is If we assume a dependence upon k 2 for the excited-state valence wave function (instead of the actual dependence upon ξ and k ⊥ , separately), i.e. the same dependence found in phenomenological valence wave functions widely adopted for describing ground states (as the one exploited in the discussion of the nucleon form factors in Ref. [18]), then the behaviors shown in Figs. 2 and 3, for the node structures and asymptotic behaviors of the states, become quite reasonable. Indeed, the assumed excited-state valence wave functions have to display a node at a fixed value for k 2 , and therefore according to Fig. (2), for increasing k ⊥ the variable ξ is constrained to approach 1/2 (i.e. the maximal value of ξ(1 − ξ)), in order to take (almost) constant k 2 . In conclusion, the correlation between the LF components, ξ and k ⊥ , in determining the node position can be largely explained by the rotational invariance of the phenomenological wave functions, if they depend upon k 2 . Notably, our calculations, genuinely in Minkowski space, actually confirm the overall expectation, based on a simple phenomenological Ansatz, that takes into account the rotational invariance. It should be reminded that, within a LF framework, the rotational invariance can be fully recovered only if the whole Fock expansion is considered [18]. We just add that the second node present in the right panel of Fig. 2 is hard to be seen given the scale of the plot.    Fig. 4. The asymptotic k ⊥ behaviors of the first (left frame) and second (right frame) excited states are shown, using the same label convention as given in Fig. 3.
From Eq. (5), one can obtain the asymptotic behaviors, k ⊥ → ∞, of LF valence wave function for bound states, which are given by Such behavior is explicitly shown in Fig. 4 for the first two excited states (for ground-state, see Ref. [13]). It should be emphasized that independently of the value of ξ the wave function is damped as ∼ k 4 ⊥ .
We close this subsection by discussing the equivalence of the transverse-momentum amplitudes in Minkowski and Euclidean spaces [15], respectively defined as where Φ E (k E , p) is obtained from Φ(k, p) after applying the Wick rotation with k 0 → ik 0 E .
Notably, within NIR, one can easily prove that φ T M (k ⊥ ) = φ T E (k ⊥ ), since one can exploit the explicit expression of the analytic dependence of the BS amplitude, as given in Eq. (3). As a matter of fact, choosing the rest frame, one can straightforwardly see that the zeros of the Nakanishi denominator in the complex plane of k 0 are given by: with z ∈ [−1, 1] and γ ∈ [0, ∞]. Therefore, Eq. (3), as a function of complex k 0 , has two cuts, with branch-points at Recalling that κ 2 is positive for bound states, one can show that at the branch-point k b 0+ a cut starts in the upper half-plane for Re k 0 < 0, while at the branch-point k b 0− the cut is placed in the lower half-plane for positive values of Re k 0 , as shown in Fig. 5. If, in Eq. (17), where the Minkowski space is adopted, one considers the integration variable k 0 as a complex one, i.e. k 0 = |k 0 |e iθ , and rotates the angle θ up to 90 o , no singularities are crossed (cf. Fig. 5). Furthermore, assuming that the BS amplitude drops out fast enough for large complex |k 0 |, the Cauchy theorem holds and the Wick rotation [2] can be applied for computing the transverse amplitude. Namely, one can adopt a new integration path, along a purely imaginary k 0 , without dealing with any singular integrals. Consequently, the Minkowski and Euclidean transverse amplitudes, given by Eqs. (17) and (18) are formally equivalent.
The quantitative comparison for the cases µ/m=0.1 and 0.5, with α gr taken from Table 1 (recall that one has always B(0)/m = 1.9), is presented in Fig. 6, showing a very good agreement between the transverse amplitudes, within the accuracy of our numerical approaches. It is worth noticing that such an equivalence gives an additional confidence in NIR, since it should be emphasized that the Euclidean solutions of BSE are not obtained by assuming the NIR for BS amplitudes. Therefore the comparison in Fig. 6 should be considered as a further check of the reliability of the NIR itself, at least at the ladder level, besides the passed tests for both eigenvalues [6,7,9,13] and scattering lengths [14]. Moreover, Fig. 6 illustrates nice and general features of the transverse amplitudes, that appear when the binding energies change. As a matter of fact, the position of the node in the first excited state moves toward smaller values of k ⊥ as the binding energies decreases, i.e. from the left panel (B(1)/m = 0.22) to the right panel (B(1)/m = 0.0082). Analogously, the amplitudes themselves decrease more quickly in momentum space. Both features can be explained by the increase of size of the bound state when the binding energy decreases.

Valence LF wave function in the impact-parameter space
The transverse charge densities have been thoroughly discussed by Miller in Ref. [28], in close relation to the elastic electromagnetic form factor. Indeed, the transverse charge density allows one to properly generalize the well-known non-relativstic relation between form factor and density to a relativistic framework. As a matter of fact, it turns out that for a composite bosonic state, the form factor F (Q 2 = −q 2 ) can be written as where (i) the momentum transfer q µ is evaluated in the Breit frame with q + = 0, (ii) Q 2 = q 2 ⊥ , (iii) b belongs to the transverse plane, called IP space [19], and (iv) ρ(b) is the IP density. It has to be pointed out that the IP density is the sum of contributions from all the LF amplitudes of the Fock expansion of the interacting-system state, such that ρ(b) = ρ val (b) + higher Fock states densities · · · .
The valence term is defined through the valence wave function in the IP space, φ(ξ, b), as follows with normalization (cf. Eq. (14)) d 2 b ρ val (b) = P val . In Eq. (23), the IP-space valence wave function is the 2D Fourier transform of ψ(ξ, k ⊥ ), given by where φ(ξ, b) results to be symmetric with respect to 1 − 2ξ, as a consequence of the already discussed symmetry of g(γ, z; κ 2 ) under the transformation z → −z. Moreover, one can deduce the general behavior for large transverse separations, b = |b|, as illustrated in what follows. By performing the 2D Fourier transformation, the IP-space valence wave function can be written within NIR for the s−wave state as follows where with J n (x) the integer-order Bessel function for n = 0. Also the integration over γ can be analytically carried out, leading to where K 1 (x) is the modified Bessel function of the second kind. The function F (ξ, b) exponentially drops out for b → ∞. Such a behavior can be understood by a close analysis of Eq. (27). First, from the physically-motivated request [28] that φ(ξ, b) is finite for b → 0 (see also [13]), such that one can deduce that g(γ, 1 − 2ξ; κ 2 ) must vanish for γ → ∞. Therefore, the relevant interval of γ in the integral (27) can be taken effectively finite. Exploiting such an observation, one can extract the driving exponential fall-off of F (ξ, b) in the asymptotic limit b → ∞. In this limit K 1 (x) reads: The leading exponential behavior in the integral (27) comes from values of e −b √ γ+κ 2 +(ξ−1/2) 2 M 2 [as seen from Eq. 29)] with γ close to 0. Therefore, where the exponential fall-off is singled out and the reduced function f (ξ, b) should decrease more smoothly for large values of b. It has to be pointed out that an exponential fall-off is expected for bound states, since it is generated by short range interactions, in analogy with the behavior found for the nonrelativistic two-dimensional case. The above feature has been investigated by actually calculating F (ξ, b), and in turn f (ξ, b), for both ground and firstexcited states. In Fig. 7, f (ξ, b) is presented for the ground (left) and firstexcited (right) states. In both cases, we have µ/m = 0.1 and α gr =6.437. It is worth noting in the right panel the nice node structure of the excited state in the whole {1 − 2ξ, b} plane. The tail of the function F (ξ, b) with respect to b is studied in more detail through f (ξ, b), putting in evidence that the steep fall-off of the valence wave function, which is largely taken into account by the leading exponential term, included in the definition (30). This suggests at most a polynomial behavior in f (ξ, b), which is clearly seen in the figure for b < 15/m.

Conclusions and Perspectives
We have investigated both spectrum and excited states of the scalar Bethe-Salpeter equation, in ladder approximation, by getting, for the first time, solutions directly in Minkowski space, within the Nakanishi integral representation of the BS amplitude. A basic ingredient of our approach is the exact projection of the BSE onto the null-plane (see, e.g. [15,16]), that allows one to master in a simple and very effective way the singularities typical of the BS formalism. We have considered an s-wave interacting system composed by two massive scalars and interacting through a massive scalar, extending the study of the ground state performed in Ref. [13] (see Ref. [9] for an analogous study within the explicitly-covariant LF approach), and carefully analyzing the valence wave function both in Minkowski and impact-parameter spaces.
Within the numerical accuracy of our approach, we have found a finite number of excited states for non-zero exchanged-scalar mass, and we have successfully compared our results with the corresponding ones obtained in the Euclidean space, where obviously NIR is not assumed. A detailed study of the valence wave function structure has been carried out in the plane (ξ, k ⊥ ), showing the expected node structure of the first and second excited states. Furthermore, our investigation, both analytic and quantitative, of the transverse-momentum amplitude has allowed to remarkably show the equivalence of the quantity evaluated both in Euclidean and Minkowski spaces. This further strengthens the reliability of the approach based on NIR for solving the BSE, already applied to fermionic systems [11], kernels beyond the ladder one [10] and in the zero-energy scattering case [14]. Finally, we also explored the asymptotic properties of the impact-parameter space valence wave function for large transverse distances, where an exponential fall-off was singled out (similar to the nonrelativistic case in the 3D Euclidean space) and quantitatively tested for the excited states.
In perspective, the present study encourages the extension of the approach based on the NIR to excited states of actual physical systems, as well as to explore results obtained for other dynamical quantities within a wider and deeper comparison between Minkowski and Euclidean calculations.