{'fraction_non_alphanumeric': 0.05287419399094526, 'fraction_numerical': 0.012566881602414598, 'mean_word_length': 4.035368886432717, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving the Bethe-Salpeter Equation for Scalar Theories in Minkowski Space

The Bethe-Salpeter (BS) equation for scalar-scalar bound states in scalar theories without derivative coupling is formulated and solved in Minkowski space. This is achieved using the perturbation theory integral representation (PTIR), which allows these amplitudes to be expressed as integrals over weight functions and known singularity structures and hence allows us to convert the BS equation into an integral equation involving weight functions. We obtain numerical solutions using this formalism for a number of scattering kernels to illustrate the generality of the approach. It applies even when the na\"{\i}ve Wick rotation is invalid. As a check we verify, for example, that this method applied to the special case of the massive ladder exchange kernel reproduces the same results as are obtained by Wick rotation.


I. INTRODUCTION
There has been considerable recent interest in covariant descriptions of bound states, for example, in conjunction with model calculations of high-energy processes such as deep inelastic scattering. A fully covariant description of composite bound states is essential for the understanding of hadronic structure over the full range of available momentum transfer. In a relativistic field theory the two-body component of a bound state is described by the appropriate proper (i.e., one-particle irreducible) three-point vertex function or, equivalently, by the Bethe-Salpeter (BS) amplitude. An extensive review of the BS equation has been given by Nakanishi [1]. Other reviews of BS equation studies and many further references can be be found in Refs. [2][3][4][5]. The BS amplitude satisfies an integral equation whose kernel has singularities due to the Minkowski metric. The resultant solutions are not functions but mathematical distributions. The singularity structure of these distributions makes it difficult to handle the BS equation numerically in Minkowski space.
We illustrate the BS equation for a scalar theory in Fig. 1, where Φ(p, P ) is the BS amplitude, P ≡ p 1 + p 2 is the total four-momentum of the bound state and p ≡ η 2 p 1 − η 1 p 2 is the relative four-momentum for the two scalar constituents. We have then M = √ P 2 for the bound state mass and also η 1 + η 2 = 1, but otherwise the choice of the two positive real numbers η 1 and η 2 is arbitrary. Note that p 1 = η 1 P + p and p 2 = η 2 P − p. In the nonrelativistic limit the natural choice is η 1,2 = m 1,2 /(m 1 + m 2 ), e.g., η 1 = η 2 = 1/2 for the equal mass case. This is the choice that we will make here.
where the kernel function defined by I(p, q; P ) ≡ −iK(p, q; P ) is the form typically used by Nakanishi [1]. In ladder approximation for a φ 2 σ model we see for example that I(p, q; P ) = g 2 /(m 2 σ − p 2 − iǫ). For ease of comparison with earlier work we will use the latter form of the BS equation here, i.e., Eq. (3).
One approach to dealing with the difficulties presented by the Minkowski-space metric is to perform an analytic continuation in the relative-energy variable p 0 , which is the so-called "Wick rotation" [9]. This has been widely used as a means of solving model BS equations, e.g., see Ref. [8] and references therein. The special case of the ladder BS equation is solved as a function of Euclidean relative momentum in the standard treatment. This is possible since in the ladder approximation the kernel is independent of the the total four-momentum P µ . The difficulties associated with this approach in the general case arise from the fact that since the total four-momentum must remain timelike then P · p becomes complex. In addition, when one uses a "dressed" propagator for the constituent particles or more complicated kernels in the BS equation, the validity of the Wick rotation becomes highly nontrivial. For example, essentially all of the dressed propagators studied previously in Wick-rotated Dyson-Schwinger equation approach contain complex "ghost" poles [8]. Hence, the naïve Wick rotation obtained by the simple transcription of Minkowski metric to the Euclidean metric and vice versa is not valid in general. So while the Euclidean-based approaches certainly play a very important role and can be useful in model calculations, it is preferable to formulate and solve the BS equation directly in Minkowski space. Here we present a method to solve the BS equation without Wick rotation by making use of the perturbation theory integral representation (PTIR) for the BS amplitude [10].
The PTIR method was first used in conjunction with a Wick rotation by Wick and Cutkosky for a scalar-scalar bound state with a massless scalar exchange in the ladder approximation [9,11]. This is now commonly referred to as the Wick-Cutkosky model. They solved the BS equation in terms of a single variable integral representation, which is a special case of the PTIR. Wanders [12] first introduced the more general two-variable integral representation to solve the BS equation for a scalar-scalar bound state with a massive scalar exchange in the ladder approximation. It was shown that all invariant scalar-scalar BS amplitudes have this integral representation in the ladder approximation [13,14] and the weight function for the corresponding BS amplitude was solved formally by means of Fredholm theory [14,15]. This two-variable integral representation was independently proposed for the vertex function on the basis of axiomatic field theory [16][17][18]. Following this proposal of an integral representation for the vertex function, Nakanishi made a detailed and systematic study of the PTIR [10].
The PTIR is a natural extension of the spectral representation for a two-point Green's function to an n-point function in a relativistic field theory. Since the Feynman parametric integral always exists for any Feynman diagram in perturbative calculations, one can always define the integral representation such that the number of independent integration parameters is equal to that of invariant squares of external momenta. Each Feynman diagram contributes to the weight distribution of the parametric integral, so that the complete weight function for the renormalized n-point function is identical to the infinite sum of Feynman diagrams for the renormalized Lagrangian of the theory. Hence, we see that the PTIR of a particular renormalized n-point function is an integral representation of the corresponding infinite sum of Feynman diagrams for the renormalized theory with n fixed external lines. This paper is organized as follows. In Sec. II, we introduce the PTIR of the renormalized scattering kernel for scalar models without derivative coupling. We also introduce the PTIR for the BS amplitude itself. We discuss the structure of the kernel weight function and the BS amplitude. In Sec. III an integral equation for the BS amplitude weight function is derived and expressed in terms of a new type of kernel function, resulting from the kernel weight function and the constituent scalar propagators. The structure of this new kernel function are discussed in Sec. IV. Numerical results are presented in Sec. V and we present our conclusions and directions for future work in Sec. VI.

II. PTIR FOR SCALAR THEORIES
In this treatment we will limit ourselves to studies of bound states of scalar particles interacting through a scalar kernel without derivative coupling. Let φ(x) be the field operator for the constituent scalar particles having a renormalized mass m. We will define µ such that m + µ is the threshold for particle production in the single φ channel. For example, for a φ 2 σ model where the renormalized σ-particle mass satisfies m σ < 2m or m σ > 2m, then we would have µ = m σ or µ = 2m respectively. The renormalized propagator for the φ-particle can be written in the following spectral form; where ρ φ (α) is the renormalized spectral function. Note that ρ φ (α) ≥ 0, (see, e.g., [6]). Following the conventions of Ref. [7], (e.g., pp. 481-487), we can also define centre-ofmomentum and relative coordinates X ≡ η 1 x 1 +η 2 x 2 and x ≡ x 1 −x 2 such that x 1 = X +η 2 x, x 2 = X − η 1 x, and P · X + p · x = p 1 · x 1 + p 2 · x 2 . Hence the Bethe-Salpeter amplitude Φ(p, P ) for the bound state of two φ-particles having the total momentum P ≡ p 1 + p 2 and relative momentum p ≡ (η 2 p 1 − η 1 p 2 ) can be defined as where we have made use of the translational invariance of the BS amplitude. Equivalently, we can write Note that the bound states are normalized such that P |P ′ = 2ω P (2π) 3 δ 3 ( P ′ − P ), where ω P ≡ ( P 2 + M 2 ) 1/2 with M the bound state mass. For a positive energy bound state we must have P 2 = M 2 , 0 < P 2 ≤ (2m) 2 , and P 0 > 0. The normalization condition for the BS amplitude is given by where the conjugate BS amplitudeΦ(p, P ) is defined bȳ

A. PTIR for Scattering Kernel
The scattering kernel I(p, q; P ) ≡ −iK(p, q; P ) describes the process φφ → φφ, where p and q are the initial and final relative momenta respectively. It is given by the infinite series of Feynman diagrams which are two-particle irreducible with respect to the initial and final pairs of constituent φ particles. For purely scalar theories without derivative coupling we have the formal expression for the full renormalized scattering kernel where q 2 i is the 4-momentum squared carried by φ i and s,t and u are the usual Mandelstam variables. This expression has been derived by Nakanishi (PTIR) [10]. Since only six of these are independent due to the relation q 2 1 + q 2 2 + q 2 3 + q 2 4 = s + t + u, this is also the number of independent ξ-parameters. Hence, we need only introduce the "mass" parameter γ and six dimensionless Feynman parameters ξ i with a constraint. The symbol Ω denotes the integral region of ξ i such that Ω ≡ {ξ i | 0 ≤ ξ i ≤ 1, ξ i = 1 (i = 1, . . . , 6)}. The "mass" parameter γ represents a spectrum of the scattering kernel. The function ρ ch (γ, ξ) gives the weight of the spectrum arising from three different channels which can be denoted ch = {st}, {tu}, {us}. This PTIR expression follows since any perturbative Feynman diagram for the scattering kernel can be written in this form, and hence this must also be true of their sum. In a perturbative calculation the weight function ρ ch is calculable as a power series of the coupling constant for a given interaction Lagrangian. For more general theories involving, e.g., fermions and/or derivative couplings, the numerator of Eq. (9) will also contain momenta in general. This would considerably complicate the procedure to be presented here and the problem of how to extend this approach to those cases is as yet unresolved. One can impose additional support properties on the kernel weight functions ρ ch (γ, ξ) by analyzing general Feynman integrals for the theory of interest [10]. However, we do not impose any further conditions at this stage since we wish to incorporate more general cases such as separable kernels, which cannot be written as combinations of ordinary Feynman diagrams.
To illustrate the approach we will discuss two explicit examples here: (a) Scalar-scalar ladder model with massive scalar exchange: The simple t-channel one-σexchange kernel is given by The BS equation with this kernel together with perturbative propagator D 0 is often referred to as the "scalar-scalar ladder model" [1]. The Wick rotated BS equation for this kernel has been studied numerically [19]. We use this kernel as a check of our calculations.
(b) Generalized ladder kernel: A sum of the one-σ-exchange kernel Eq. (10) and a generalized kernel with fixed kernel parameter sets {γ (i) , ξ (i) }. After the Wick rotation this kernel becomes complex due to the p · P and q · P terms, so that solving the BS amplitude as a function of Euclidean relative momentum would be very difficult in this case.

B. PTIR for BS Amplitude
Let us now consider the PTIR of the BS amplitude itself. Since the BS amplitude in the center-of-momentum rest frame [i.e., with P µ = (M, 0)] forms an irreducible representation of the O(3) rotation group, we can label all of the bound states with the usual 3-dimensional angular momentum quantum numbers ℓ and ℓ z . We can thus construct the integral representation of the partial wave BS amplitude in the rest frame and the PTIR automatically allows us to then boost to an arbitrary frame.
For simplicity, let us begin by considering the s-wave amplitude. In order to apply the PTIR to the BS amplitude, it is convenient to first consider the equivalent vertex function Γ(p, P ) defined such that iΦ(p, P ) = D((P/2) + p)Γ(p, P )D((P/2) − p). The vertex function Γ(p, P ) is a one-particle-irreducible 3-point Green's function, which can be expressed as an infinite sum of Feynman diagrams and hence in terms of PTIR has the form derived by Nakanishi [10] In contrast to the case for the scattering kernel, all invariant squares of momenta are independent and so a single weight function ρ 3 (κ, z) is sufficient to describe the sum of all allowed Feynman diagrams. From the perturbative analysis of Feynman graphs for the vertex function one obtains the following support property for the vertex weight function where "max" is taken separately for each fixed set of z. For the vertex function with a timelike total momentum satisfying 0 < P 2 < 4m 2 , it follows that Eq. (11) reduces to the following two-variable representation with the support propertỹ This two-variable representation was first proposed in axiomatic field theory [17,18]. For simplicity, we will now assume that the constituent φ-propagators have their free form, D(p 2 ) → D 0 (p 2 ) = 1/(p 2 − m 2 + iǫ), or equivalently ρ φ (α) = 0 in Eq. (4). The arguments can be extended in a straightforward way to the case ρ φ (α) = 0, but we will not pursue this possibility here. Then attaching two φ-propagators as external lines yields the PTIR for the BS amplitude where we have introduced the non-negative integer dummy parameter n. In order for the expression to be meaningful, we see that we must have the boundary conditions We see that a partial integration of Eq. (15) with respect to α connects weight functions with different n, i.e., and hence, as stated above, the integer n is an artificial or dummy parameter. While the non-negative integer n is arbitrary, a judicious choice can be advantageous in numerical calculations, since the larger we take n, the smoother the weight function becomes for a given Φ(p, P ). From the support property ofρ 3 (β, z) in Eq. (14), it follows that the weight function ϕ n (α, z) = 0 when This support property can be understood from the relation betweenρ 3 (β, z) and ϕ n (α, z), i.e., where . Thus the region of support for ϕ n (α, z), which is specified by the step-functions in Eq. (19), is due to the singularities of the φ-propagators attached as external lines to the proper bound-state vertex Γ. Note that the integral representation in Eq. (15) and the support property Eq. (18) are valid even if we keep the full nonperturbative φ-propagator D(p 2 ) by including the spectral function ρ φ (α). In this case Eq. (19) for the BS amplitude weight function ϕ n (α, z) must be suitably generalized and will include an integration over each of the spectral functions from the two φ-propagators. Now let us consider the BS amplitudes for higher partial waves, i.e, those for bound states with non-zero angular momentum (ℓ > 0). To define these, we consider the rest frame of the bound state, i.e., P µ = (M, 0). The momentum-dependent structures, i.e., the denominators, of the PTIR's for the bound state proper vertex Γ and the BS amplitude Φ are independent of transformations under the little group belonging to the 3-momentum p in this frame. Hence, for higher partial waves with angular momentum quantum number ℓ and third component ℓ z , Γ and Φ can be written as the product of the ℓ-th order solid harmonic Y ℓz ℓ ( p) = | p| ℓ Y ℓz ℓ ( p) and the corresponding PTIR for the scalar (i.e., s-wave) bound state in this frame [1]. A simple way to understand this is that since the system consists only of scalar particles and since in this frame there is only one available 3-vector (i.e., p), then higher partial wave amplitudes must be proportional to Y ℓz ℓ ( p). The solid harmonics are polynomials of their three arguments and since they have the self-reproducing property Eq. (B3) [14], we know that in fact the amplitudes must be proportional to Y ℓz ℓ ( p). Thus we have in the rest frame for a bound state with angular momentum quantum numbers ℓ and ℓ z where the total momentum P is P = ( √ P 2 , 0) and where we have introduced the lower bound α th for the integral range of α, corresponding to the minimal value of the right hand side of Eq. (18). The boundary conditions for ϕ [ℓ] n (α, z) are identical to those for the s-wave Eq. (16). The dummy parameter n can always be taken such that the loop-momentum integral of the BS equation Eq. (3) converges.
By performing a straightforward Lorentz boost we have the integral representation of the partial wave BS amplitude in an arbitrary frame where P is an arbitrary timelike 4-vector with P 2 = M 2 and p ′ = Λ −1 (P )p. The Lorentz transformation Λ(P ) connects P and the bound-state rest frame 4-vector P ′ = (M, 0), i.e, P = Λ(P )P ′ . In the following sections we will study the BS equation Eq. (3) in an arbitrary frame in terms of this integral representation.

III. BS EQUATION FOR THE WEIGHT FUNCTION
In this section we will reformulate the BS equation Eq. (3) as an integral equation in terms of the weight functions. This is the central result of this paper. Substituting the integral representation of the partial wave BS amplitude Eq. (20) together with the PTIR for the scattering kernel Eq. (9) into the BS equation Eq. (3), the right hand side of the BS equation becomes where {a ch , b ch , c ch , . . . , f ch } are different linear combinations of ξ i in each of the three channels, ch= {st}, {tu}, {us}. They are listed in Appendix A. Applying the Feynman parameterization one can perform the q-integral by setting the dummy parameter n such that n + 1 > l/2. This condition ensures that the q-integral is finite. Due to the self-reproducing property of the solid harmonics this integral is again proportional to the solid harmonics Y ℓz ℓ (Λ −1 (P )p) (see Appendix B). Now we multiply the propagators D(P/2 + p)D(P/2 − p) into both sides and absorb them into the RHS expression using Feynman parameterization. The BS equation then becomes The kernel function K [ℓ] n (ᾱ,z; α, z) is defined with a fixed parameter set of the scattering kernel {a ch , b ch , c ch , . . . , f ch } and its definition is given in Appendix C. Comparing Eq. (23) with Eq. (21), and using the uniqueness theorem of PTIR [10], we obtain the following integral equation for ϕ [ℓ] n (α, z): Here we have introduced the total kernel function K [ℓ] n (ᾱ,z; α, z), which is the superposition of K [ℓ] n (ᾱ,z; α, z) with the kernel weight functions ρ ch (γ, ξ) such that Since the weight functions ρ ch (γ, ξ) for the scattering kernel are real functions by their construction, the total kernel function K [ℓ] n (ᾱ,z; α, z) is real, so that the Eq. (24) is a real integral equation with two variables α and z. Thus we have transformed the BS equation, which is a singular integral equation of complex distributions, into a real integral equation which is frame-independent. Once one solves for the BS amplitude weight function, the BS amplitude can be written down in an arbitrary frame. This is clearly advantageous for applications of the BS amplitude to relativistic problems. As we will see the real weight functions can be in fact be real distributions. However, as is evident from Eq. (17) if n is chosen sufficiently large these can be always transformed into arbitrarily smooth functions suitable for numerical solution.
The functions h [ℓ] n (α, z), k [ℓ] n (ᾱ,z; α, z) and g [ℓ] n (ᾱ,z; α, z) are defined as Feynman parameter integrals in Appendix C. The terms containing the function h [ℓ] n (α, z) are independent of z. Since this feature is independent of the kernel parameters, the total kernel function K [ℓ] n (ᾱ,z; α, z) also has this structure for any input scattering kernel. For example, the total kernel function with n = 0 can be written as where H [ℓ] 0 (α, z) and G n (α, z) and g [ℓ] n (ᾱ,z; α, z) with the kernel weight function in Eq. (25), respectively. This structure suggests that the weight function ϕ where c ′ is a constant andφ [ℓ] n (α, z) is a function determined by the kernel. If one substitutes this form into the PTIR for the BS amplitude Eq. (15), the δ-function term gives the product of two propagators with the weight c ′ /2. As discussed in Section II B, this is just a consequence of the fact that the BS amplitude contains the kinematical singularity due to the two constituent propagators. Thus this part of the structure is totally independent of the details of the scattering kernel, (i.e, independent of the scattering kernel weight functions). Now, if the function G 1 λφ where we have introduced an "eigenvalue" parameter λ. Instead of solving for the mass of the bound state for a given scattering kernel, it is more convenient to solve the BS equation as an "eigenvalue" problem for a fixed bound-state mass parameter P 2 . We thus solve for the eigenvalue λ as a function of P 2 and the points P 2 which give λ(P 2 ) = 1 correspond to the masses of the bound state. For the weight function ϕ 0 (α, z) containing the δ-function, (i.e., the weight function Eq. (28) with non-zero c ′ ), one can rescale the weight function such that the strength of the δ-function term is unity. Hence, one can study the bound state problem numerically by iteration even for the n = 0 case. Consider the following formal expression obtained by iterating the kernel There is similarly an expansion for the eigenvalue λ which is the Fredholm series. The scalar-scalar-ladder model has been formally solved by Sato by means of the Fredholm solution. He showed that the second iterated kernel is bounded and that the Fredholm theorem is applicable for the resulting regular equation [15]. Nakanishi extended his results to arbitrary partial wave solutions [14].
To apply rigorous mathematical theorems, such as Fredholm theory, it is necessary to know in detail the singularity structure of the kernel function K [ℓ] n (ᾱ,z; α, z) for given weight functions ρ ch (γ, ξ). Hence, it is very difficult to discuss their applicability for the most general form of the scattering kernel. Instead of considering such theorems, we shall take the more pragmatic path of establishing whether or not the integral equation Eq. (24) is numerically tractable.

IV. STRUCTURE OF THE KERNEL FUNCTION
In this section we consider the structure of the kernel function given in Eq. (25). We first consider the possible singularities of the kernel function K [ℓ] n (ᾱ,z; α, z) for arbitrary n with a fixed kernel parameter set (γ, ξ), i.e., for constant {γ, a ch , b ch , c ch , . . . , f ch }. We will in this section omit for brevity the subscript ch. Apart from the trivial δ-function for n = 0, which becomes a step-function for n = 1, the possible singularities of K [ℓ] n (ᾱ,z; α, z) are those of the functions h [ℓ] n (α, z), k [ℓ] n (ᾱ,z; α, z), and g [ℓ] n (ᾱ,z; α, z). Let us first consider the function g [ℓ] n (ᾱ,z; α, z). The integration over the Feynman parameter is easily performed due to the δ-function. Reflecting the integral range of the Feynman parameter y ∈ [0, ∞), we see that g [ℓ] n (ᾱ,z; α, z) has finite support as a function ofᾱ, α,z and z. The function g [ℓ] n (ᾱ,z; α, z) becomes after the Feynman parameter integral n (ᾱ,z; α, z) where the symbol Pf · stands for the Hadamard finite part. Any singularities that arise from the Feynman parameter integrals should be regularized using the Hadamard finite part prescription. This is consistent with the ordinary iǫ prescription for a calculation of Feynman diagrams in momentum space [10,14]. A discussion of the finite parts of singular integrals and the detailed calculation of g [ℓ] n are given in Appendix D. The variables y κ i for i = 1, 2 are the roots of the equationÃ(α, z;z)y 2 +B κ (α, z;ᾱ,z)y +C κ (ᾱ,z) = 0. See Appendix C for the definition of the functionsÃ(α, z;z),B κ (α, z;ᾱ,z),C κ (ᾱ,z), ζ κ (z, z) and η κ (z). Since the factor ζ κ (z, z)y κ i + η κ (z) in the denominator of Eq. (33) is positive definite for any positive y κ i , this factor does not cause any singularity. Thus forC κ (ᾱ,z) < 0, g [ℓ] n (ᾱ,z; α, z) is regular everywhere except the pole atᾱ = 0, sinceÃ(α, z;z) is non-negative for any α, z andz. ForC κ (ᾱ,z) ≥ 0 the square root of the following factor in the denominator can vanishB whereD κ (ᾱ,z, z) is the regular function defined in Appendix C. This singularity occurs at the boundary of the step function. Since a square root singularity is integrable, this singularity is numerically tractable forC κ (ᾱ,z) > 0. In the general case one can apply the Hadamard finite part prescription (D5) in Appendix D for a regularization. However, when two roots coincide, namely at the point (ᾱ,z) on whichC κ (ᾱ,z) = 0, special care is necessary to perform the α-integral. One may handle this singularity by the regularized expression Eq. (D7) with the finite regularization parameter ǫ. Thus ϕ [ℓ] n (α, z) may have singularities even for a constant kernel parameter set in the general case. Note that this singular structure of g [ℓ] n (ᾱ,z; α, z) is independent of n and l. Now let us turn to the functions h [ℓ] n (α, z) and k [ℓ] n (ᾱ,z; α, z). Since the step function in k [ℓ] n (ᾱ,z; α, z) restricts the integral range of y for a given parameter set, it is enough to consider the integral n (y min , y max ) ≡ ymax y min dy y n+1 (a + y) n−1−l [A(α, z)y 2 + B(α, z)y + C] n+1 (35) in order to discuss the singularity structure. It is easy to show that Eq. (35) can be written as with the functioñ The boundary of the integral y min and y max should be understood as constant, when taking the derivative with respect to α andã, although they are fixed as functions of variables by the step function of Eq. (C16). Since the y-integral should be performed using the Hadamard finite part prescription, the above operation is valid even though I [ℓ] n (y min , y max ) contains singularities. The integral range [y min , y max ] for k [ℓ] n (ᾱ,z; α, z) depends onᾱ,z, α and z and for h [ℓ] n (α, z) is fixed to [0, ∞). Since A(α, z) > 0 for any z and α that is consistent with the support property of the weight function Eq. (18), the integrand always vanishes sufficiently fast for y → ∞. Thus a singularity occurs only if the set of variables;ᾱ,z, α and z, satisfy the condition where the denominator of the integrand Ay 2 + By + C vanishes for some y, namely B 2 − 4AC ≥ 0 1 . For B 2 − 4AC ≥ 0, the integration over y yields and Y = ã y+ã 2 (Ay 2 + By + C). As shown in Appendix C the factor a 2 A −ãB + C becomes a positive definite function D(z) depending only on z after setting a = a. The limit B 2 − 4AC → 0 is also regular, since X is proportional to the factor √ B 2 − 4AC. Thus a singularity occurs only when 1 ± X or Y vanish at the end point of the y-integration. Furthermore, neither the derivative with respect to α nor that with respect tõ a generates any new singular point. They simply change the power of the singular behavior. Since X is independent ofã and ∂Y ∂ã = 2 y y+ã Y , the derivative with respect toã acting oñ I(y min , y max ) does not cause stronger singularities. On the other hand, the derivative with respect to α increases the inverse power of Y and 1 ± X. Thus the possible singularity of the integral I [ℓ] n (y min , y max ) is of the form Pf · 1 (Ay 2 min + By min + C) m or Pf · 1 (Ay 2 max + By max + C) m , where m is a non-negative integer which does not exceed n and where we have used the fact that 1/(1 ± X) ∝ 1/(1 − X 2 ) ∝ 1/(Ay 2 + By + C). It is not difficult to see that the fixed end points (y = 0 and y = ∞), which are independent of the choice of parameters, cause no further singularities provided B(α, z) or C does not vanish. Thus the function h [ℓ] n (α, z) has no singularity, if C = 0. On the other hand, k [ℓ] n (ᾱ,z; α, z) may have singularities. A possible end point of the y-integral, which might give the singularity, is; the positive root y κ i , (i = 1, 2) of the equatioñ A(α, z;z)y 2 +B κ (α, z;ᾱ,z)y +C κ (ᾱ,z) = 0, and the boundary value of y at which κ changes sign: y bd = −a + b 2z b/2−e+az cz−f +b/2 when y bd ≥ 0. Since the Eq. (40) can be written as follows and the factor ζ κ (z, z)y+η κ (z) is positive definite, the end point singularity at y = y κ i is given by Pf · 1/ᾱ n . Thus the possible singuralities of k [ℓ] n (ᾱ,z; α, z) are Pf · 1/ᾱ n and Pf · 1/(Ay 2 bd + 1 I

[ℓ]
n (y min = 0, y max ) at a = 0 can be singular. Such a singular point is independent of the choice of the variables. One can always avoid it by choosing the dummy parameter to be some n ≥ l + 1.
To summarize, the possible singularities of the kernel function K [ℓ] n (ᾱ,z; α, z) for a fixed kernel parameter set {a ch , b ch , c ch , . . . , f ch } are: a trivial δ-function singularity atᾱ = 0 for the choice of dummy parameter n = 0; Pf · 1/ᾱ for n ≥ 1; Pf · 1/(Ay 2 bd + By bd + C) m with 1 ≤ m ≤ n for n ≥ 1; the singularity of Eq. (34) at the boundary of the support for all n and l.

V. NUMERICAL RESULTS
In this section we present numerical solutions for the BS amplitude for bound states in scalar theories using Eq.(24) for two simple choices of scattering kernels: (a) pure ladder kernel with massive scalar exchange, and (b) a generalized kernel combined with the above pure ladder kernel.
The scattering kernel (a), i.e., the one-σ-exchange kernel depending only on t = (p − q) 2 , is given by Eq. (10). This corresponds to choosing for the kernel in Eq.  22)]. This corresponds to choosing ρ st to be some appropriate product of δ-functions multiplied by g 2 . In the pure ladder case it is convenient (and traditional) to factorize out the coupling constant g 2 and a factor of (4π) 2 , by defining the "eigenvalue" λ = g 2 /(4π) 2 [15,14,19]. Thus it is usual to fix the bound state mass P 2 and then to solve for the coupling g 2 , which is what we have done here. In the numerical calculations it is necessary to regularize the integrable square root singularity discussed in the previous section (see Eq.(34)), i.e., we kept the regularization parameter ǫ in the Hadamard finite part Pf·1/x n small but finite, typically ∼ 10 −5 or less. Of course, we verified that our solutions were independent of ǫ provided it was chosen sufficiently small. We began by discretizing the α and z axes and then solved the integral equation by iteration from some initial assumed weight function. It was also, of course, confirmed that the solution was robustly independent of the choice of starting weight function, the number of grid points, and the maximum grid value of α (when the latter two of these were chosen suitably large). We have solved over a range of P 2 and m σ and all solutions reproduced the eigenvalues obtained in Euclidean space after a Wick rotation by Linden and Mitter [19]. We found that this simple method produces convergence in typically 10 cycles or so and for the relatively modest grid choices used gave results which agreed to within approximately 1% of the Wick-rotated ones. In Fig. (2) we plot the weight function φ n=1 (α, z) for the bound state with mass given by P 2 /m 2 = (M/m) 2 = 3.24, where m is the mass of the constituent particles and where the mass of the exchanged particle is m σ /m = 0.5. This solution was generated by solving for n = 0 and then using Eq. (17) to convert to the n = 1 form of the solution. This has the advantage of filtering out numerical noise and the resulting solution is more readily graphically represented. The coupling constant for this case was found to be g 2 /(4π) 2 = 1.02, c.f., the Linden and Mitter result of 1.0349.
For a nontrivial example of the scattering kernel, where the naïve Wick-rotation is not possible, we studied the following simple generalization of the pure ladder case: a sum of the one-σ-exchange kernel (a) and the generalized kernel term with the following fixed kernel parameter set weighted with a factor (1/4)[g 2 /(4π) 2 ]. In addition, the resulting total kernel was symmetrized inz → −z so that we obtain appropriate symmetry for s-wave normal solutions. We obtained these parameters by generating ξ i in the st-channel with a random number generator. We then solved for λ = g 2 /(4π) 2 as for the pure ladder case but for this case the n = 1 was solved for directly, without first obtaining the n = 0 solution. In Fig. (3) we plot the weight function φ n=1 (α, z) with the mass parameters again chosen as P 2 /m 2 = 3.24 and m σ /m = 0.5. Some additional structure is apparent, as expected, and some residual numerical noise is also present in this figure. Methods of reducing this numerical noise by finding more sophisticated choices of grid and integration method are being investigated [20]. The resulting eigenvalue is λ = 0.889, which shows that the addition to the pure ladder kernel has enhanced the binding and so is attractive. We found that the convergence and the stability of the eigenvalues with varying the number of grid points, ǫ, and α max are excellent as for the pure ladder case above.

VI. SUMMARY AND CONCLUSIONS
We have derived a real integral equation for the weight function of the scalar-scalar Bethe-Salpeter (BS) amplitude from the BS equation for scalar theories without derivative coupling. This was achieved using the perturbation theory integral representation (PTIR), which is an extension of the spectral representation for two-point Green's functions, for both the scattering kernel [Eq. (9)] and the BS amplitude itself [Eqs. (20,21)]. The uniqueness theorem of the PTIR and the appropriate application of Feynman parameterization then led to the central result of the paper given in Eq. (24). We demonstrated that this integral equation is numerically tractable for both the pure ladder case and an arbitrary generalization of this. We have verified that our numerical solutions agree with those previously obtained from a Euclidean treatment of the pure ladder limit.
This represents a potentially powerful new approach to obtaining solutions of the BS equation and additional results and applications are currently being investigated [20]. The separable kernel case should be studied, since it has exact solutions and so is a further independent test of the approach developed here. While our detailed discussions were limited to the equal mass case, it is worthwhile to investigate generalizations including nonperturbative constituent propagators and the heavy-light bound state limit to see under which conditions an approximate Klein-Gordon equation can result. It is also important to find a means to generalize this approach to include derivative couplings and fermions.

APPENDIX B: SOLID HARMONICS
Let us consider the integral Applying the Feynman parameterization we have I(n, m, ℓ; p, P ) where F , H and J are F = A(α, z) and A(α, z), B(α, z) and C are defined in Appendix C. Introducing a new variable and recalling the fact Λ −1 (P )P = ( √ P 2 , 0), one can factorize the solid harmonics outside the integral by making use of the following property of the solid harmonics; where F is a sufficiently rapidly decreasing function which gives a finite integral [14]. The integration over q then yields I(n, m, ℓ; p, P ) where we have introduced a new integration variable y such that y = x/(1−x). Thus the integration over the loop momentum reproduces the same solid harmonic. This self-reproducing property of the solid harmonics ensures that the integral representation, Eq. (21), forms an irreducible representation of the Poincare group.
Finally we list some useful linear combinations of the functionsÃ(α, z) and B(α, z); D(z) ≡ a 2 A(α, z) − aB(α, z) + C Changing the integration variable such that Y = ã y +ã 2 (Ay 2 + By + C), (D10) the integral for B 2 − 4AC ≥ 0 can be written as follows; where Y min , X min , . . . are the values of Y and X at y = y min , y max . It is easy to derive the Eq. (38) by performing the integral over X and Y with the definition of the Hadamard finite part (D1).     Figure 3