Bound state equation for the Nakanishi weight function

The bound state Bethe-Salpeter amplitude was expressed by Nakanishi using a two-dimensional integral representation, in terms of a smooth weight function $g$, which carries the detailed dynamical information. A similar, but one-dimensional, integral representation can be obtained for the Light-Front wave function in terms of the same weight function $g$. By using the generalized Stieltjes transform, we first obtain $g$ in terms of the Light-Front wave function in the complex plane of its arguments. Next, a new integral equation for the Nakanishi weight function $g$ is derived for a bound state case. It has the standard form $g= N g$, where $N$ is a two-dimensional integral operator. We give the prescription for obtaining the kernel $ N$ starting with the kernel $K$ of the Bethe-Salpeter equation. The derivation is valid for any kernel given by an irreducible Feynman amplitude.


Introduction
The Bethe-Salpeter (BS) equation [1] for a few-body system takes into account the relativity and, implicitly, all the possible intermediate states. As demonstrated in [2] for the case of the ladder kernel -which displays the strongest (pole) singularities -it is possible to find the bound states and the scattering solutions in Minkowski space by a direct numerical calculations, in spite of the many singularities appearing in the integral and inhomogeneous terms.
There exist however another efficient approach to solve this equation, proposed in [3] and developed in a series papers [4,5,6,7,8,9], which relies on the Nakanishi representation [10] of the BS amplitude Φ(k, p) in terms of a non-singular weight function g(γ, z) as follows: where κ 2 = m 2 − M 2 4 , m is the constituent mass and M is the total mass of the bound state (p 2 = M 2 ). Once g is known, one can compute the BS amplitude by means of (1) and calculate other observables, like the electromagnetic form factors (see [11]).
The weight function g determines also, via another integral representation, the light-front (LF) wave function [4]: In this one-dimensional representation, variable z plays the role of a parameter, varying in the limits −1 ≤ z ≤ 1. The relations between variables (γ, z) and the standard LF variables (k ⊥ , x) are given by γ = k 2 ⊥ and z = 2x − 1. By Eq. (2) one can obtain the LF wave function and calculate its contribution to the form factors as well as the momentum distributions.
In a series of recent papers [4,5,6,7,12,13], two independent methods for finding g have been developed. In the first one [4,5,6,7] the Nakanishi weight function g is determined by solving the double integral equation, previously established in [4]: where the kernel V is expressed via the BS kernel. The solution g of equation (3) was found in the case of an OBE ladder BS kernel in [4] and for the ladder plus cross-ladder one in [5,15]. Equation (3) contains integral terms in both sides, though the integral in the l.h. side is one-dimensional. This fact generates some instability of its numerical solution, specially taking into account that when discretizing the left-hand side it results into an ill-conditioned matrix [4]. An equation for g, getting rid of the left-hand side integral term of (3), was first derived in [6]. This derivation is based on the uniqueness of the Nakanishi representation and it was fulfilled for the ladder kernel only. To clarify this point, let us represent equation (3) in a symbolic form as: In ref. [6], the equation (4) was first transformed into an equivalent one having the form L·(g − V·g) = 0. This transformation used essentially the properties of the ladder kernelV . Its straightforward application to the same equation with a non-ladder kernel would meet some algebraic difficulties. The uniqueness theorem [10] claims that if L·f = 0, it follows that f = 0.
In this way, one obtains the equation written in terms of a unique integral operator V. The main result of [6] consisted in finding equation (5) and the corresponding kernel V. This equation was also numerically solved in [6] and the binding energies were in very good agreement with the solutions of Eq. (3) found in [4]. In addition, the stability of results relative to method and rank of discretization was much better.
In the second method [12,13], the Nakanishi weight function g is extracted from the "Euclidean" 1 BS amplitude Φ E . Since, as mentioned, the equation for Φ E can be solved rather easily and, in addition, Φ E is accessible by the lattice calculations in a full QFT framework, one can try to invert the representation (1) and find in this way, the Nakanishi weight function g. Once known, it can be used to calculate the Minkowski amplitude Φ M and the observables in Minkowski space. Notice that the Nakanishi weight function g is the same for both the Minkowski and Euclidean BS amplitudes but depends on the total mass of the system M through the parameter κ. All the scalar products, resulting in the singularities, are contained in the denominator of (1).
Alternatively, since the LF wave function ψ LF can be also found relatively easy from the corresponding LF equation (given e.g. in [14]), one can try to find g by inverting the representation (2). Here again, the procedure is complicated by the instabilities of the discretized equations (1) and (2) relative to g(γ, z) since both of them belong to the class of so called ill-conditioned problems. We would like to emphasize however that both equations are invertible in the mathematical sense, and so they accept a unique solution g. Only its numerical inversion, when required with some accuracy, is far from being straightforward and requires some inspired methods.
The aim of this paper is to derive from (3) an equivalent equation written in the normal form g = N·g. To this aim, we will not rely on the uniqueness of the Nakanishi representation neither on the ladder kernel, but we will find the explicit form of the inverse operator L −1 in (4) in a general case (Section 2). We will first apply this inverse operation to (2) and express g in terms of the LF wave function (Section 3). Applying the same result to equation (4), we will derive in Section 4 an explicit representation of N = L −1 V . in terms of a one dimension integral on V . Our derivation is valid for any kernel K in the original BS equation or for any kernel V in equation (3). In the case of a ladder kernel we will check the identity N ≡ V.

Analytical inversion of the Stieltjes kernel
The results that follow are based on the observation, that the left-hand-side of (3), and the right-hand-side of (2), can be trivially related to the Stieltjes transform: 2 This integral transform can be analytically inverted [16,17] for an arbitrary ρ. In the particular case ρ = 2 its inversion is given by: where C is a contour in the complex plane of w that starts and ends at the point w = −1 enclosing the origin in the counter-clockwise sense. Two possible choices are illustrated in Fig.  1. We will take hereafter (left panel) the unit circle centered at the origin C(0, 1) although the contour displayed in the right panel could be well adapted for further calculations. The variable w can then be parametrized as w = exp(iφ) and the integration over the contour is reduced to the integration over φ in the limits −π < φ < π. That is Let us now consider an integral relation of the form denoted symbolically as f =L g. By introducing the shifted variable γ ′′ = γ + b, a straightforward application of the inverse Stieltjes transform (8) applied to (9) gives Though the proof of Eq. (7) is given [17] in the general case (6), we will give below a simple demonstration for the particular case of the operator (9) we are interested in.
To this aim, let us write (10) as kernel of the operatorÛ ≡L −1L is written for convenience in the form: by introducing the parameter ǫ, being understood that once the integration (11) is performed, the limit ǫ → 0 should be taken. The integral (11) can be analytically computed and gives: In the limit ǫ → 0, the kernel U(γ, γ ′ ; ǫ) is a sharp function of γ ′ with a maximum at γ ′ = γ given by Therefore, U(γ, γ ′ ; ǫ → 0) is a representative of a delta-function: that isÛ is the unit operatorÛ =L −1L = 1 and so (10) is the inverse of (9). It is interesting to illustrate the inverse equation (10) by means of an anaytically solvable model. This is provided, for simplicity in the case b = 1, by the following pair (f, g) satisfying Eq. (9): By inserting f given by (14) into Eq. (10), the integral can be performed analitically and we indeed reproduce g from (15). Care must be taken however when using a numerical method to perform an integration which involves complex logarithms (see the remark at the end of next section).

Computing the Nakanishi weight function in terms of the light-front wave function
A first straigthforward application of this result can be done by inverting Eq. (2), which express the LF wave function ψ LF in terms of the Nakanishi weight function g. By setting in (10) and introducing explicitly the z-dependence (which plays the role of a parameter) on g, one has g(γ, z) = γ 2π The Nakanishi weight function is thus formally expressed in terms of the LF wave function provided it is known in the complex plane of its first argument.
By inserting g computed by (16) in Eq. (1) one can obtain the BS amplitude Φ, both in Euclidean (Φ E ) and in Minkowski (Φ M ) space. Notice that the possibility to obtain the valence LF wave function from the Minkowski BS amplitude was already established some years ago by a projection technique [14,19]. The contents of Eq. (16) can be considered as a way to realize the reciprocal procedure.
This result is remarkable and highly non trivial since it provides a link betwen two relativistic frameworks -Light Front Dynamics and Bethe-Salpeter equation -which have a very different dynamical contents. The key stone of this link is provided by the Nakanishi weight function g.
The interest of such a possibility in different fields of hadronic physics was emphasized in our previous works [12,13] where we presented an attempt to determine g by numerically inverting a discretized version of Eqs. (1) or (2). The corresponding linear systems turned to be ill-conditioned and, apart from some analytical models like the ones described in previous section -Eqs. (14) and (15) -its determination is not free from uncertainties. Equation (16) constitutes a formal prove of such an inversion.
The practical interest of (16) deserves however further investigations. It relies on the knowledge of the LF wave function on the complex plane. This is indeed technically accessible, as it has been done in [20] for the Euclidean BS equation, but it is not an easy task and it is not clear that could be more efficient that the numerical inversions described in [12,13]. It is worth noticing that g is real for bound states, a property which is far form being explicit in (16). This suggests that other integration contours than C(0, 1) in (7) taking benefit of this property, like the one on the right panel of Fig. 1 with R → ∞, could be more suitable.

Deriving the integral equation for g
By applying the integral transform (10) to both sides of the equation (3) and using eq. (13), we obtain the following equation for the Nakanishi weight function g in the canonical form More explicitly where The operatorN =L −1V in (17) exists if the corresponding integral (19) converges. Some divergence could eventually appear if the integration contour crosses a singularity of V in the complex plane. In this case, the infinitesimal imaginary "−iǫ" term added to masses in the Feynman amplitudes, as usual, is equivalent to a contour deformation allowing to avoid crossing singularities. Since V is related to a Feynman amplitude, it contains −iǫ, the kernel N(γ, z; γ ′ , z ′ ) is finite and the operatorN =L −1V in (17) exists.
The equation (18) determines also a relation between the coupling constant α in the kernel N and the total mass M (or the binding energy B = 2m − M) of the system or its spectrum (if any). As usual, the solution g is found for a given M, and so although no explicitly written g depends on the binding energy.
The relation between the original kernel K appearing in the BS equation and the kernel V in (19) and (3) was derived in Ref. [4]. In the case of a one boson exchange kernel V takes the form: where: and The integral (22) was calculated analytically in [4] and reads: By inserting (21) in the integral (19), the kernel N of Eq. (18) will be determined. Notice however that the benefit of using the analytic expression (23) with respect to the numerical integration (22) is not always granted, neither from the computational point of view nor from the ambiguities arising when using multivalued complex integrands. As mentioned in the Introduction, the canonical form of Eq. (3) was first derived in [6] under the hypothesis of the ladder kernel. These authors wrote this equation -Eq. (24) of [6] -in the symbolical form (5) with a kernel V that, when written in a form close to (21), reads 3 Function h, obtained from h ′ given in [6] by Eqs. (27) and (28), takes the form: where The integrand on C and Q reads and the integration limits are given by y ± = −B±∆ 2A . The argument η of the θ-function in the first term of (26) is: In order that Eqs. (17) and (5) could represent the very same equation, the kernels N and V must be identical to each other. However the complexity of their analytical forms, specially taking into account that N results from the integral (19), makes any direct comparison beyond the scope of the present paper.
To prove this identity we must entirely rely on numerical calculations. More precisely, we will substitute the kernel V defined by (21) in Eq. (19), calculate the integral over φ for a selected set of external variables (γ, z, γ ′ , z ′ ), determine N(γ, z, γ ′ , z ′ ) and compare the result with V(γ, z, γ ′ , z ′ ) given by (25).
We have considered two of these sets, corresponding to negative and positive values of η, the argument of the θ function appearing in Eq. (26). The kernel OBE parameters were fixed to α = 1, m = 1, µ = 0.15, M = 1.9.
For the first set of variables one has η = −1.24 < 0, and the θ function in (26) does not contribute to the result. The numerical calculation gives: As one can see that the N and V kernels coincide with each other within 12 digits, the computer accuracy.
For the second set of variables Here again, the kernels coincide with each other within the same accuracy. Note that these calculations were performed using the analytical expression for the W function (23).
The following warning concerns numerical calculations involving complex multivalued functions. When using Mathematica or fortran packages to evaluate the integral (19) with analytical expressions like Eq. (23) in the integrand, one should take care with the implicit cuts on the multivalued functions (e.g. log(z), √ z, etc.) implied in these tools. For example, the result depends on whether or not we apply to r.h. side of (23) the Mathematica command PowerExpand. The kernel (33) for set (32) was obtained after applying PowerExpand; without applying it, the two kernels are different. 4 Whereas the result for set (30) is independent on using PowerExpand 4 Examples: or not. Note that the equation (22) for W in the integral form does not contain this ambiguity. After numerical calculation of the double integral (over φ and v), N coincides with V given by (25), though with less precision. The results described in this Section confirm the identity of the two kernels N and V, despite their very different formulation and analytical expressions. This identity between these two kernels, proved in the case of a ladder approximations, suggests to use the integral form (19) over V as an alternative way to compute V in a more involved dynamical context or for scattering case. Indeed, the validity of Eq. (18) does not rely on the uniqueness ansatz and has been established under the same conditions than Eq. (3), that is for any kernel K inserted in the BS equation.
The same remark concerning the real character of g discussed at the end of the previous Section holds here in the definition of N itself through the complex integral (19). The solution g as well as the kernel are, in the bound state case, both real and it would be interesting to take benefit on that and make it explicit with an appropriate choice of the contour in Fig. 1.

Conclusion
We have considered the Nakanishi integral representation of the Bethe weight function g. These integral transforms are difficult to invert since they result into ill-conditioned linear systems.
By using an analytic inversion of the Stieltjes integral transform we have first obtained a closed form expressing the Nakanishi weight function in terms of the Light-Front wave function in the complex plane.
Starting with a Bethe-Salpeter equation for the weight function g in the form Lg = V gsee Eq. (3) -previously established in [4,5] we have used the same inversion procedure to write it in the canonical form g = Ng -see Eq. (18) -with the kernel N given by an integral of kernel V .
This equation generalizes, to the case of a kernel given by any set of irreducible Feynman graphs, the equation found in [6] valid only for the ladder kernel and using the uniqueness anzatz. This generalization strongly expands the applicability of the Nakanishi representation to find the solution of the Bethe-Salpeter equation, as well as its applications.
It is worth noticing that the validity of applying the inverse Stieltjes transform for deriving Eq. (18) and solving in this way the Bethe-Salpeter equation relies on the fact that the Bethe-Salpeter kernel is generated by irreducible Feynman diagrams. This allows to use the Nakanishi representation for the Bethe-Salpeter amplitude and obtain in this way the equation (3). When The result depends on the order of factors in the argument of log, on using or not the PowerExpand (!) command and on how log is defined. Chose e.g. z = exp −i π 4 = 0.707(1 − i). The value of log(z) depends on the branch one takes: either log(z) = log exp −i π 4 = − iπ 4 = −0.785i, or log(z) = log exp i − π 4 + 2π = 7iπ 4 = 5.498i. The same happens with √ z in (24) when using Eq. (23). This reflects the fact that we are dealing with multivalued function. Since the numerical calculation of an integral, like (22), is unambiguous we take it as the correct value.
using an arbitrary kernel -e.g. constructed by hand or inspired on phenomenological grounds, however providing converged results -the Nakanishi representation for the corresponding amplitude could not be defined and therefore the applicability of the Stieltjes transform would be not guaranteed.
The first evident application of the method developed here, could be to find the solution g for the ladder plus cross-ladder kernel, previously obtained from Eq. (3) in [5,15]. The equation (18) can also be extended straightforwardly to the two-fermion system starting, for example, from [8].
The possibility to apply an analytic inversion of the Stieltjes transform can be useful in other fields of nuclear and hadronic physics, where the use of integral transforms was pioneered by [21,22,23], in order to avoid the instabilities of the numerical inversion.