{'fraction_non_alphanumeric': 0.043970935258522, 'fraction_numerical': 0.01616784654564146, 'mean_word_length': 4.210225197808886, 'pattern_counts': {'":': 1, '<': 14, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Study of relativistic bound states for scalar theories in Bethe-Salpeter and Dyson-Schwinger formalism

The Bethe-Salpeter equation for Wick-Cutkosky like models is solved in dressed ladder approximation. The bare vertex truncation of the Dyson-Schwinger equations for propagators is combined with the dressed ladder Bethe-Salpeter equation for the scalar S-wave bound state amplitudes. With the help of spectral representation the results are obtained directly in Minkowski space. We give a new analytic formula for the resulting equation simplifying the numerical treatment. The bare ladder approximation of Bethe-Salpeter equation is compared with the one with dressed ladder. The elastic electromagnetic form factors is calculated within the relativistic impulse approximation.


I. INTRODUCTION
In the Quantum Field Theory the two body bound state is described by the threepoint bound state vertex function or, equivalently, by Bethe-Salpeter (BS) amplitudes, both of them are solutions of the corresponding (see Fig. 1) covariant four-dimensional Bethe-Salpeter equations (BSE) [1]. Up to now most of the studies were restricted to the case when irreducible interaction kernel is approximated by (sum of) single particle exchanges. In this so called ladder approximation the scattering matrix is given by the sum of the generated ladders. It is known, that such approximation is not sufficient when more realistic models are considered [2], [3], [4]. To move beyond this approximation one is in practice confined to the use of some phenomenological Ansatze. (In hadronic physic, these Ansatze are very often made already at the level of two point correlators. For modeling of the gluon propagator in the context of BSE and Dyson-Schwinger equations (SDE), see for instance [5]).
Here we are considering some simple scalar models. An extensive review of BS studies in scalar theories (with at most cubic and non-derivative interaction) can be found in Ref. [6], [7]. The various improvements of simple ladder kernel have been considered, in particular, including the self-energy effects [8], [9] or contributions from crossed box diagrams [10]. The study of the influence on the bound state spectrum following from the infinite resummation of certain ladder and crossed-box diagrams can be found in the paper [11]. Furthermore, there is a number of interesting papers on the solution of Wick-Cutkosky models (with zero mass of the exchanged particle). These solutions employ various effective techniques like the point form of relativistic quantum mechanics [12], variational calculations [13] or the light front dynamics [14].
The standard approach to determine the spectrum and the BS vertex makes use of partial wave decomposition which reduces four-dimensional integral equation into the twodimensional one. The alternative more recently exploited treatment is based on the O(4) hyperspherical expansion [15]. In this approach the BSE is transformed into an infinite set of one-dimensional integral equations. Notable advantage of this approach is a good numerical convergence and easy identification of excited bound states spectra.
Very often, the ladder BSE's are solved with the help of the so-called Wick rotation [16]. However, the backward analytical continuation is quite difficult even for the ladder approximation, while for more complicated cases its proper implementation is unclear or at least highly non-trivial.
In this work we follow the method of solving the BSE directly in Minkowski space [8], [17], [18], in which the problems associated with Wick rotation do not arise. The method is based on utilization of generalized spectral representation for n-point Green functions in quantum field theory [19]. In this treatment the BSE written in momentum space is converted into a real integral equation for a real weight function with number of independent variables dependent on details of the model. We extend the earlier work [8] first to the case with unequal masses of constituents. This then allows us to treat the ladder BSE in which all propagators (of constituents and of the exchanged particle) are fully dressed. This is achieved by the implementation of Lehmann representation of the propagator: G(p 2 ) = dωσ(ω)D(p; ω) , D(p; ω) = 1 p 2 − ω + iǫ ;σ(ω) = R δ(ω − m 2 pole ) + σ(ω) . (1.1) where σ(ω), which is a smooth function, nonzero above a threshold, is determined by the Dyson-Schwinger equations (DSE). In the pole term we can take R = 1, a choice corresponding to the conventional on-shell renormalization scheme (G has a unit residuum when momentum approaches a simple pole at physical mass).
To account for the effect of self-energy we transform the momentum BSE to the form suitable for complementary solution together with the appropriate DSE for propagators. Note here, that the perturbative one loop contribution has been already considered in [8] and certain Euclidian version of this problem has also been investigated [9]. In qualitative agreement with [9] we have found that the critical value of coupling gives the domain of applicability of BSE (at least in its ladder approximation). The couplings below critical one allow only solutions for relatively weakly bound states. It is even more interesting that the effect of the propagator dressing on bound state spectra is rather small. In comparison with the bare ladder approximation the same binding energy is then achieved with the coupling smaller by about several per cents, even for values of the coupling close to the critical one.
Clearly, when we take some or even all particle propagators dressed, the number of spectral integrations increases. Note here, that up to the rather exotic case of massless Wick-Cutkosky model the appropriate solution is not known analytically but must be found numerically. Mainly due to this reason we re-formulate the equation obtained by Kusaka et al [8] and we offer the solution where the appropriate integral kernel is free of any additional numerical integration (see Appendix of Ref. [8] for the original solution). The elimination of this numerical integration then not only improves numerical accuracy but also reasonably decreases the CPU time.
To see explicitly the effect of radiative corrections we compare the derssed BSE results with its bare ladder approximation. We set the parameters of our model to that used in Refs. [15] and [17], to compare the bare ladder solutions to those obtained before in [8], [15], [17].
Having solved equations for spectral functions, one can determine the BS amplitudes in an arbitrary reference frame. This makes this technique suitable for calculations of response to an external fields. In section IV. we briefly introduce the formulas defining the elastic charge form factor G(Q 2 ) in relativistic impulse approximation (RIA). Although the elastic form factor represents a simple dynamical observable, its Minkowskian calculation represents nontrivial task. For this purpose we consider the (massive) Wick-Cutkosky model given by Lagrangian gauged as follows: where covariant derivative is D µ = ∂ µ − ieA µ . In our form factor calculations the effects of scalar dressing was not taken into account, since it would significantly increase computational complexity of the problem. Furthermore, it is assumed that e << g which implies that the interaction of the charged particle field φ 1 with the electromagnetic field can be treated perturbatively. As in Ref. [9], we have chosen the same coupling constant for interaction of the field φ 3 with the fields φ 1 and φ 2 . The form factors were calculated for several sets of masses of constituent and exchanged particles.

II. DRESSED LADDER BETHE-SALPETER EQUATION
The BS amplitude for bound state (φ 1 , φ 2 ) in momentum space is defined through the Fourier transform of Here p 1,2 are the four-momenta of particles corresponding to the fields φ 1,2 that constitute the bound state (φ 1 , φ 2 ). The total and relative momenta are then given as P = p 1 + p 2 and p = (η 2 p 1 − η 1 p 2 ), respectively, and P 2 = M 2 , where M is mass of the bound state. Finally, P · X + p · x = p 1 · x 1 + p 2 · x 2 . From now on we will put η 1 = η 2 = 1/2, which corresponds to the usual separation of center of mass motion for equal mass case, but can be also employed for unequal masses (although X is then not the coordinate of the center of mass).
Introducing the BS vertex function Γ = iG −1 1 G −1 2 Φ, the homogeneous BSE for S-wave bound state reads The bound states appear as poles of the scattering matrix. The normalization condition for the BS vertex function follows from the requirement that the pole appropriate to a given bound state is a simple one: whereΓ(p, P ) is the conjugate of Γ(p, P ) and p 1,2 = ±p + P/2, p ′ 1,2 = ±k + P/2 (for the details see e.g. Ref. [6]).
In this work we do not solve the BSE with most general irreducible scattering kernel V (the most general structure of the kernel V written in terms of its Perturbation Theory Integral Representation (PTIR) can be found in Ref. [8] or [19]). Here we restrict ourselves to the case of dressed ladder approximation with φ 3 -exchange, in which V (p, k; P ) = g 2 G 3 (t), where t denotes the usual Mandelstam variable t = (p − k) 2 . Note that for the bound state of particles (φ 1 , φ 2 ) only the t-channel interaction above is effective, whereas for the bound state of (φ i , φ i ), i = 1, 2 one has to consider also possible u and s channel diagrams. In the present work we study the case (φ 1 , φ 2 ) due to its simplicity and leave the other cases for discussion elsewhere. Let us also recall that in the (dressed) ladder approximation, defined by the exchange of chargeless scalar φ 3 , the photon coupling to φ 1 alone (one particle current, or in other words RIA) is by itself gauge invariant, when taken between the corresponding solutions of BSE. And the normalization condition (2.3) in this approximation reduces to the condition G(Q 2 ) = 1. We will solve the BSE for massive constituents (m 1,2 > 0) and m 3 ≥ 0, the invariant mass of the bound state satisfies 0 ≤ P 2 < (m 1 + m 2 ) 2 .
Taking the dressed kernel and the full propagators of constituent particles into account, the rhs of BSE (2.2) can be written as The interesting unequal-mass ladder case of Ref. [20] is also described by Eq.(2.4), although η 1 = η 2 = 1/2 and p and k are not relative momenta. Since the dependence on momenta in (2.4) is explicit, one might always re-scale Γ(p, P ) to proper relative momentum. The integral representation for the BS vertex function may be written as [19] Γ(p, (2.5) The positive integer n represents a free parameter without clear physical meaning. One can take advantage of this freedom of choice to pick up n so that the numerical solutions of integral equations for spectral functions are made more stable. The spectral functions ρ [n] (α, z) for different n can be related by integration over α by parts. Kusaka et al [8] choose n = 2 for their numerical solution of the BSE, we adopt the same value in this paper. The bare (symmetric) Wick-Cutkosky model corresponds to the choice: α 1 = α 2 = m 2 , the exchanged boson is massless (α 3 = 0) and no radiative corrections are considered. This model is particularly interesting because it is the only example of the nontrivial BSE which is solvable exactly [16]. For this model, there is no freedom in choice of n and (for the S-wave bound state) the expression (2.5) reduces to the one dimensional PTIR: Using a technique similar to the one used in ref. [8], the BSE can be converted to the following real integral equation for the real spectral function: where we denoted λ = g 2 /(4π) 2 . The derivation is presented in the Appendix A, where the explicit expressions for particular choices n = 1, 2 are given. The central results of this work are expressions obtained for V [n] , which are simpler that the ones presented in Ref. [8]. No additional integration is required which decreases the computer time necessary for numerical calculation. Besides, our formulas hold also for unequal masses of the constituents. The extension (together with above mentioned simplification) to the case of more complicated scattering kernel is not so straightforward, but we believe that it is possible. Note also, that due to the property of solid harmonic with respect to the integration over the momentum the presented procedure can easily be generalized for the bound state with non-zero spin (here, the total orbital momentum) [8].
A bound state with equal composite masses is described by the vertex function Γ which is symmetric under the transformation P · p → −P · p. In terms of the weight function this symmetry reads ρ(α, z) = ρ(α, −z). However, there are solutions that do not respect this symmetry even in the case of equal masses. These are usually called ghost solutions and the appropriate amplitudes have a negative norm. Such solutions are often considered to be nonphysical and it is supposed that they point at inner inconsistency in description of relativistic bound states within the BS formalism, at least in the ladder approximation. At this place, it is important to mention that the Lagrangian (1.2) describes the models that are a subset of theories with potentials unbounded from below and in very strict sense they are discarded due to the vacuum instability. On the other hand one can assume, at least for sufficiently small couplings the existence of local minima of the potentials is sufficient to support of the existence of ground state of the theory. While in the large coupling regime, say for g/m >> 1 no reasonable physics can be learned from the perturbation theory and/or from equation like ladder BS by itself. To conclude, we note that such ideas are supported by at least two facts. The ghost BS solutions do appear only for a large value of λ. Furthermore, from the Dyson-Schwinger study we know that the scalar theory studied here makes sense only up to the certain critical coupling, see e.g. [9], [21].
The important question arises, what is the validity of the full theory when renormalization is properly taken into account. Although the quantitative answer lies beyond the approximations used in this article and requires more careful investigation, we make a simple attempt to find the domain in which self-consistent solutions of the BSE and Dyson-Schwinger equations (within the framework of reasonable approximations) exist. Furthermore, we inquire an influence of scalar propagator dressing on the solution of the BSE for the bound states.

III. DRESSING PROPAGATORS BY THE DYSON-SCHWINGER EQUATIONS
The solution of the DSE for the scalar models with the help of spectral decomposition will be discussed in detail in our forthcoming paper [21]. Here we give only brief presentation of the DSE in bare vertex approximation, their renormalization, re-arrangement in terms of the spectral function and some properties of solutions, important for our further discussion of the BSE. In this section we assume m 3 < m 1 + m 2 , so that the propagator of the exchanged particle Φ 3 has an isolated physical pole.
By dressing of the scalar propagators in our study of BSE we mean only the dressing due to the "strong" interaction between scalars, the coupling to the electromagnetic field is neglected. Let us now write the strong interaction part of our Lagrangian (1.2) in terms of bare, unrenormalized quantities (fields and coupling constants), labeled by subscript "0": In the previous section we have chosen the strength of of both couplings to be the same.
Here, we distinguish the bare couplings, anticipating that they are renormalized by different amounts (see below). The kinetic terms are parametrized by the unrenormalized masses m 0i . These masses undergo the infinite mass renormalization To re-scale the residuum of the full propagators to unity, we will complement the infinite mass renormalization by the finite (since the model is superrenormalizable) renormalization of the fields and coupling constants That is, we will employ below the on shell renormalization scheme in which the propagators have unit residua when momentum approaches its mass shell value p 2 → m 2 .
In this paper we consider the Dyson-Schwinger equations in the simplest approximation in which the proper vertices are replaced by the bare ones Γ oi = g oi . Then, the DSE in their unrenormalized form read and making use of (3.3), one gets the re-scaled DSE The renormalization of proper self-energies proceeds by double subtraction: Identifying the appropriate renormalization constants (3.2,3.3) we can immediately write the full propagator in terms of finite physical quantities The DSE for the renormalized self-energies are given by (3.5) and subtraction (3.6).
For the purpose of our BS calculation, we now fix the renormalized couplings and masses as follows where g is coupling constant from (1.2). That is, we will compare the solutions of the BSE for the bare and dressed ladder kernel taken for the same numerical value of unrenormalized and renormalized coupling constant, respectively. The masses are fixed to allow comparison with some of the results of Refs. [8,15,17]. Now, it is straightforward task to evaluate the spectral representation of the renormalized self-energy. Lehmann representation (with unit residuum) for G i reads: Notice that functions σ i , i = 1, 2, have the thresholds at m i,th = (m i + m 3 ) 2 = 2.25m 2 , whereas the function σ 3 has the threshold at m 3,th = (m 1 + m 2 ) 2 = 4m 2 . Analogously, for the for self-energies: The spectral representation for Π R explicitly satisfies Π R (m 2 ) = Π ′ R (m 2 ) = 0 following from (3.6). Rewriting now the relation between G and Π in the form G = D + DΠG (D being the free propagator with the physical mass) and taking its imaginary part, we arrive to the first relation between the spectral functions σ and ρ: 12) where P stands for principal value integration. All the functions in (3.12) are positive and regular above the perturbative thresholds and identically equal to zero elsewhere.
Substituting the spectral representations (3.10) into the DSEs (3.5), making the subtraction as in (3.6) and comparing to the lhs in the form of (3.11), one gets after lengthy algebra: where λ = g 2 /(4π) 2 and the function B(α, β; ω) is related to the Källen function λ as follows: (3.14) Before the numerical treatment the explicit integration-separating the δ-function parts of Lehmann weightsσ-has to be performed. Equations (3.12,3.13) constitute the closed system of integral equations for spectral functions which can be solved numerically by iterations without any additional approximation. So obtained dressed propagators have been used when solving the Bethe-Salpeter equation. The results are discussed in the section V. Before leaving this section we review some important features of our solutions of DSE.
The behavior of the imaginary parts of propagators-the Lehmann functions σ i (α)-for fields Φ 1,2,3 is shown at Fig. 2.
The renormalization constant Z i are calculated from the relation (3.15) From the Fig.3 we can see that the field renormalization constant Z 1,2 changes sign from positive to negative at some critical pointλ crit = g 2 crit /(4πm) 2 ≃ 1.5 ± 0.1, where the error reflects the difficulty of making the numerical estimate of the value for which the solution cannot be found and the dimensionless coupling is defined asλ = λ/m 2 . We did not find any numerical solutions of DSE for couplings larger thanλ crit . It is reasonable to suppose that the quanta associated with the fields φ 1,2 do not describe physical particles whenλ >λ crit .

IV. ELASTIC ELECTROMAGNETIC FORM FACTOR
The electromagnetic form factors parameterize the response of bound systems to external electromagnetic field. The calculation of these observables within the BS framework proceeds along the Mandelstam's formalism [22]. For the elastic scattering on the S-wave bound state, (P 2 i = P 2 f = M 2 ) the current conservation implies the parameterization of the current matrix element G µ in terms of the single real form factor G(Q 2 ) The elastic electromagnetic form factor G(Q 2 ) depends only on the square of photon incoming momentum q and we use the usual SLAC convention Q 2 = −q 2 , so that Q 2 is positive for elastic kinematics.
The matrix element of the current in relativistic impulse approximation (RIA) is diagrammatically depicted in Fig. 4. In this paper we are not taking into account the dressing of the scalar propagators when calculating the charge form factor. Then, the matrix element is given in terms of the BS vertex functions as where we denote P = P i and j µ 1 represents one-body current for particle φ 1 , which for the bare particle reads j µ where p i , p f is initial and final momentum of charged particle inside the loop in Fig. 4, i.e., p i = k + P/2, p f = q + k + P/2.
We have already mentioned in the previous section, that if the vertex functions Γ are solution of the BSE with a kernel corresponding to exchange of single chargeless particle, the RIA defined above is by itself gauge invariant and the normalization condition for the BS amplitudes is equivalent to the normalization G(0) = 1.
The main result of this paper, as far as charge form factor is concerned, is re-writing of the rhs of Eq. (4.2) directly in terms of the spectral weights of the bound state vertex function. It allows the evaluation of the form factor by calculating the integral of nonsingular expression, without having to reconstruct the vertex functions Γ(p, P ) from their spectral representation. The derivation of this integral involves some lengthy algebra and is relegated to Appendix B.

A. Bare Ladder BSE
We have solved the bare ladder BSE for symmetric (m 1 = m 2 = m) scalar theory with bare ladder kernel and bare constituent propagators G i (p i ) = D i (p i , m i ) by iterations of the integral equation for spectral functions. The standard procedure was followed: after fixing the bound state mass (P 2 ) we looked for the solution by iterating spectral function for fixed dimensionless "coupling strength"λ ≡ g 2 /((4π ) 2 m 2 ). If the iterations failed-measure being both the difference of the rhs and lhs of the integral equation and deviation of the auxiliary normalization integral from pre-defined value-we were changingλ (halving intervals of successive guesses) till the solution was found.
In case of Wick-Cutkosky model the one dimensional integral equation (A24) was solved. Although the solution of this one-dimensional integral equation could be found by the inversion of its discretized form, we have tested the iteration procedure (used later also for massive exchange). The equation was discretized by numerical Gauss integration, it appear that it is sufficient to take 40 Gauss point (though the number cited in Table 1 are obtained with 98 points). It is known [6] that P 2 = 0 corresponds toλ = 2 from which our result slightly deviates in the fifth digit. We have also reproduced (up to four published digits) all results for Wick-Cutkosky model from [15]. In the Fig.5 the weight functions are plotted against spectral variable z for several fractions of binding η = √ P 2 /2m . For the massive scalar exchange the two-dimensional integral equations (A31,A35) were solved. We have found (in agreement with [8]) that numerical errors are about one order of magnitude bigger for n = 1, hence n = 2 is preferable and only the results with this choice are discussed below.
For numerical solution we discretize integration variables α and z using Gauss-Legendre quadratures (with tangent mapping from < −1, +1 >→< α min , ∞ > for α ). The equation (A35) is solved on the net of N = N z * N α points which are spread on the rectangle (−1, +1) * (α min , ∞). The value α min is given by the support of the spectral function (see the appendix A). We have not optimized the grid during the iteration procedure as it was done in the study [8]. Instead, we have solved the equation for several different numbers of grid points while keeping fixed the ratio of N α /N z and then extrapolated the results to the "ideal" case with N α = N z = ∞. Examples of numerical convergence for some cases of bound states are presented in the Table 2. In the Table 1 we compare our results for m 3 = m/2 with those of ref. [8].
Below we show the dependence of charge form factor on the parameters of the model: on the range of interaction characterized by the inverse mass of exchanged meson m 3 and on the strength of forces which bind the particles together. To calculate the form factors we first have to solve the BSE for chosen sets of parameters. We vary the parameters as follows: 1. First we solve the BSE for several bound state masses P 2 keeping the ratio m 3 /m fixed 2. Then we vary the mass of the exchanged meson, keeping the masses of all studied bound states fixed (our choice is η = √ P 2 /(2m) = 0.5) and determining the corresponding coupling strengthsλ.
Where the independent numbers were available [15], they agree with our results (see Table  3). If the mass of exchanged particle becomes small (but nonzero) the convergence of our numerical procedure becomes somewhat poorer and more sensitive to the initial guess. For illustration the weight functionsρ [2] is plotted in figs. (6,7) for two different values of exchanged mass.

B. Dressed Ladder BSE
In this subsection we finally discuss numerical solutions of the BSE including the dressing of the propagators. We introduce the dressing by two steps, switching it on first only for the exchanged particle and in the second step also for constituents. We should point out that since the solution of the DSEs in the bare vertex approximation (by which we dress the propagators) breaks down for coupling constants larger thatλ crit = 1.5 we can consider to only rather weakly bound states: for the bare BSEλ = 1.5 corresponds to η = √ P 2 /(2m) ≃ 0.78. The propagator dressing of the exchanged particle in the one loop approximation was already considered in the Ref. [8]. We go beyond the one loop approximation and determine the continuum part of Lehmann weight σ 3 from the DSE (3.12) with the same value of the coupling constant. That is, we insert into Eq. (2.7 ) the dressed kernel with the pole situated at m 3 = m/2. As noted above, the constituent propagators are at this stage left undressed in the BSE, although they all self-energies have been taken fully into account in DSEs. The integration over ω (5.2) in the BSE kernel was performed using Gaussian quadrature with 16 points. Including the kernel self energy slightly decreases (by at most few per cent) the mass of the bound state, even forλ ≃λ crit. . Let us point out, that the kernel of Eq. (2.7) in the dressed ladder kernel approximation is free of any singularities. The accuracy of the numerical solution is comparable to the the bare ladder case. For example, weλ = 0.734 for η = 0.95 for the grid of 32 * 32 points and λ = 0.749(0.752) for the grid of 64 * 64(96 * 96), the convergence is similar to to the case of bare ladder (see the Tab.1). The extrapolated (to very large grid) values of λ ′ s for fractional binding η = 0.999, 0.99, 0.97, 0.95 are showed in Fig.8.
In the next step we have included the self energies of the constituents. As we shall see the effect is relatively small forλ <<λ crit , but increases rapidly asλ → 1. As in the previous case the Lehmann weights have been calculated from the DSEs solved for the same value of the coupling. As the first guess we have used the solution of the BSE linearized inσ(α), i.e. with only one propagator dressed. This guess is rather close to the exact solution for η ≤ 0.9.
The constituent particle in a weakly bound system (η ≃ 1) live near their mass shell. Therefore, one can naively assume that the values of coupling for such weakly bound state should not be strongly affected by dressing of constituent propagators. For deeper bound states we have found that the effect of the dressing of constituents is much larger than that due to the dressing of the kernel (exchanged particle) Fig.8. The couplings for fully dressed BSE are not determined with the same high accuracy as those for the bare ladder BSE, since the grids are not optimalized for very different ratios of α's which appear in the kernel of the BSE.

C. Charge form factor
Various form factors are extensively studied in scalar theories like Wick-Cutkosky model (see for example [23]). In these studies the dependence of the form factor on the binding and on the range of the "strong" interaction has been considered, therefore we perform a similar calculations in our formalism.
In the approach adopted in this paper (employing the spectral representations in the Minkowski space), the bound states masses and corresponding vertex functions can be obtained with good accuracy and in reasonable CPU time. Unfortunately, the calculation of the scalar form factor as outlined in appendix B leads to more complicated results (B11,B12). Even if one would be able to perform analytically all additional Feynman integrations, the formula (B12) still involves the four-dimensional integration over the spectral variables. We are taking also the integrals over four Feynman variables numerically with the help of Gauss-Legendre quadrature, taking the number of points for each of them equal to the one for spectral variable z. Since relatively small number of integration points (from 16 to 40) was taken for each integration, the presented results have to be viewed rather as an estimate of form factor behavior. One can always refine the grids at the expense of longer CPU time.
We have also compared our results to choose obtained in the Gross (spectator) formalism, choosing the "scalar deuteron" parameters (see [24]). In analogy with the real deuteron, the parameters are chosen as: m 3 /m = 138/938.9, η = (2 * 938.9 − 2.3)/2 * 938.9 ≃ 0.9988. The bound state vertex functions were found by solution of the Gross and BS equations. All phenomenological form factors introduced in [24] have been "switched off" (the limit Λ ′ s → ∞ are taken in the "strong" form factors) when calculating the Gross wave function and the bound state current. The e.m form factors were calculated in the spectator RIA and as described in appendix B (using the grid 32 8 ), respectively.
The form factors for several bound states listed in the Tables 1 and 3 are presented in Figs.9 and 10, respectively. In Fig.9 the ratio of the exchanged and constituent mass is kept constant m 3 /m = 0.5 and the mass of the bound state M = √ P 2 is varied. In agreement with physical expectations one sees that as the bound state become more tight the elastic form factor increases. (For the infinitely bound point system is should be equal to unity, even our deepest bound state are of course still only approaching this value.) We put into the same plot for comparison also the scalar deuteron result. How the form factror changes with the mass of the exchanged particle is shown in Fig.10. We included several states for which M = m, bound by the one-boson-exchange potential of the range r = 1/m 3 , which is varied. Two other systems, first with η = 0.4; m 3 = 0.25 and the scalar deuteron, are added for comparison. From both figure we can conclude that the behavior of form factor is determined by the strength of the interaction and its range rather differently for various M. The range of the interaction is more significant for a weaker ones. This agrees with conclusions of ref. [10], which have also compared our results for small exchanged mass η = 0.8, m 3 = 0.15m with Wick-Cutkosky model prediction for η = 0.784, m 3 = 0 (Tables 2  and 1 of [10], respectively), and found only slight difference in the range Q 2 = (0, 100m 2 ).

VI. CONCLUDING REMARKS
The spectral representation was employed for solving the Bethe-Salpeter equation in (3+1) Minkowski space. The new analytical formula for the integral equation kernel has been derived.
The method is efficient solving both the bare and dressed ladder BSE. Solution of Dyson-Schwinger equations for propagators leads to appearance of critical value of the coupling constant, beyond which the solution collapses. This restricts substantially the region in which the effects of dressing can be studied. Since the coupling is rather weak, the dressing leads only to moderate decrease of the bound state masses: even close to the critical value of the coupling the fractional binding of the bound state of the dressed BSE is smaller than the corresponding one for the bare BSE by at most 15 per cents. As an example of application of the obtained vertex functions, we calculated the elastic electromagnetic form factor.
To further develop the method, it would be interesting to extend it to more complicated BS kernel: trying to include the cross boxed contributions, s and u channel interactions etc. It is already known that the "spectral" approach used here is suitable even for more complicated systems, for scalar QED see Ref. [18]. One of our future goals is to manage the complication due to fermionic degrees of freedom.

APPENDIX A: KERNEL FUNCTIONS
In this appendix the real integral equation for the BS vertex weight is derived in detail. The PTIR form for scalar bound-state vertex reads where ρ [n] (α, z) is the real PTIR weight function for the bound state vertex function, and n is a dummy parameter. The function F is given by [19] F (α, z; p, P ) = α − (p 2 + zp The support of ρ [n] (α, z) can be determined in general (see [19]) for arbitrary interaction. In our case of one-boson-exchange interaction kernel one gets a bit higher α min [8]. We discuss our treatment of the lower bounds below. The following procedure is straightforward but a bit exhaustive. The rhs of the BSE (2.2) with the kernel given by the exchange of the single (dressed) particle Φ 3 has to be re-written in the form allowing to extract the integral equation for the spectral function ρ [n] (α, z). The integrand contains the two "constituent" propagators, the denominator F (α, z, k, P ) from the spectral representation of Γ(k, P ) and the propagator of the exchanged particle (all other factors will be skipped for a while for the sake of briefness).
Using the Feynman parameterization technique we first write Next the denominator of (A1) is added: where z ′ = tz + (1 − t)η. Now, we include the propagator of the exchanged particle, the integral over d 4 k and factors ig 2 , defining: where Q = (1 − x)p − xz ′ P/2 and S = 1−z ′2 4 P 2 . Since z ′ lies in the interval < −1, +1 >, 0 ≤ S < (m 1 + m 2 ) 2 /4. Interchanging the integrals over η and t with the help of: and introducing λ = g 2 /(4π) 2 we get Let us separate the t dependence of F (A, z ′ ; p, P ). First we substitute for η into the definition of M 2 , which yields Next, we introduce the notation (indicating explicitly the t-dependence of A and R): in which the t−dependence of F (A(t), z ′ , p, P ) reads Since t-dependence of F (A(t), z ′ ; p, P ) is linear, the integral over t can be taken: and hence Using this result, the BSE can be written as follows: whereσ(α i ) are the Lehmann functions of the dressed propagators, reducing to the δ-function if the dressing is neglected. Obviously, the rhs is still not quite in the desired form, both F (A(0), z ′ , p, P ) and F (A(T s ), z ′ , p, P ) n are functions of momenta p and P . Below we re-write the kernel as a sum of several fractions, which after substitution α ′ = A(T ), T = 0, T ± would allow to use the uniqueness theorem [19] and extract the BSE in the spectral form. This is possible only if the integrals over α ′ on both right and left hand sides of A16 are taken over the same intervals. To show this it is first necessary to prove that R(T ) − S > 0, for T = 0, T ± , since the functions A(T ) → +∞ for x → 0 and x → 1 they have on for 0 < x < 1 the minimum equal to: next, one has to show that these lower bounds are not in conflict with the lower bound for α ′ (and the same lower bound for α). It is simple for the undressed equal mass case, when the condition above taken for T = 0 actually defines the lower bound for α in the following form and the same for α, z → α ′ , z ′ . Since R(T ± ) depend also on α and z and since for the equal mass case R(T ± ) ≥ R(0), the next two constrains clearly conform with the lower bound for α ′ and one can extract from them the upper bound for the integration over α: (compare to eq.(A6) of [8], where one bracket seems to be misplaced). For the unequal mass case or for the dressed propagators this analysis is much more complicated, mostly due to the fact that now for some combination of parameters it can occur R(T ± ) < R(0). The necessary condition R(T ) − S > 0 can again be proven in this case (though after much longer algebra), ensuring that the common lower bound for α's exists. But we could not resolve the conditions (A18) analytically, they are treated numerically. Now, let so go back to eq. (A16) and proceed by considering first the simple case of the symmetric Wick-Cutkosky model.

The Wick-Cutkosky Model
In this model the constituents have the same masses m 1 = m 2 = m and the mass of the exchanged particle is zero (all propagator dressings are neglected). Then for the Swave ground state vertex function, the spectral function depends only on variable z and the denominator enters with power n = 1 (2.6): and in I DDDF we replace Taking the integral over x we find (recall that for equal masses 0 ≤ S < m 2 ): Comparing both sides of the BSE and using the uniqueness theorem [19], we get the wellknown (see e.g. [6]) integral equation for ρ [1] (z): Although its solution is known analytically even for exited states, the energy spectrum still has to be found numerically (up to the P 2 = 0 corresponding to λ = 2). For the purpose of numerical treatment, this equation is usually rewritten in the form: where the temporary auxiliary normalization condition dzρ(z) = 1 was imposed (i.e., if used in further application, the vertex function would have to be re-normalized in accordance to (2.3)).

The BSE for n = 1
We will now bring eqs. (A15,A16) to the desired form for the particular choice n = 1. For the numerical solution this value is not the most suitable one. But since the formal manipulations are in this case the simplest, we treat it first for methodical reasons.
For n = 1 the integrand of (A15) can be decomposed into the sum of simpler fractions with the help of where we have introduced a shorthand notation Notice that the lhs of (A24) is nonsingular for J(α, z) = 0, so when this happens the rhs behaves like 0/0, which calls for some caution in the numerics. So, the integral I DDDF now reads In the last step we introduce the spectral variable α ′ = A(T ) and use the dependence of A(T ) on x to convert the integration over x into integral over α ′ . Picking up explicitly the x-dependence of A(T ) we can write: With the help of these relations we get Using this result in (A16), one gets from the uniqueness theorem the integral equation for the BSE structure function is obtained: Notice, that x ± (T ) depend on α ′ , z ′ (which are fixed from the lhs), α 1 , α 2 , α 3 (which are fixed when the dressing is neglected), P 2 in S (which is given by the binding energy of the system) and for T = 0 also on α, z (through J(α, z) in R(T ± )). This allows to recast the integral equation into the form more convenient for numerical treatment.

BSE for n = 2
Now, we would first describe the necessary modifications for the choice n = 2 which was used in actual numerical calculations in this paper. Going back to (A15) we can for this case write: .
Then, the integral I DDDF can be cast into form: where we have first integrated the second term of the sum over T by parts (to increase the power of F (A(T ), z ′ ; p, P ), the boundary term vanishes when x → 0, 1) and then introduced the integration over α ′ as in the previous section. Now, we can use the uniqueness theorem of PTIR [19] and identify the BS weight function on the right-hand side of BSE: [2] (α ′ , z ′ ; α, z)ρ [2] (α, z) , (A32) Before treating the general case with fully dressed propagator we consider in the next subsection the pure ladder BSE.

APPENDIX B: ELASTIC ELECTROMAGNETIC FORM FACTOR
In this Appendix we derive the expression for the elastic electromagnetic form factor G(Q 2 ). The relevant matrix element is diagrammatically depicted in the Fig.4 and the starting equation for the e.m. matrix element is given by (4.2). The labeling of momenta corresponds to the Fig.4: q is a virtual photon incoming momentum and we put Q 2 = −q 2 > 0, P i = P and P f = P + q are the total momenta of bound state in in-and out-state, respectively. For simplicity, we will consider the equal mass case: m 2 1 = m 2 2 = m 2 . The form factor G(Q 2 ) can easily be extracted from its definition (4.1). After multiplying by P i + P f we get The matrix element is evaluated between on-shell states of appropriate composite scalar, i.e. P 2 = (P + q) 2 = M 2 which implies 2P.q + q 2 = 0 (B2) The one body current j µ reads Using (B2) one can simplify (2P + q) · ((2k + q + P ) = 2k · (2P + q) + 2M 2 + Q 2 2 (2P + q) · ((2k + q + P ) (2P + q) 2 = 1 2 1 + 4k · (2P + q) Taking this into account we can write: Now, we will express G(Q 2 ) in terms of spectral functions of the bound state vertex functions Γ, re-writing first the product of the propagators with the help of the Feynman parametrization.
For the product of the propagators one gets: where the substitution y = s + (1 − s)η has been used and relations P · q = Q 2 /2 and P 2 = M 2 has been replaced used. Next step we combine PTIR (for n = 2) of the product of the bound state vertex functions.
Making use of the Feynman variable t for matching (B7) with the term in large brackets of (B8), the relation for the form factor (B4) can be rewritten as where we have omitted the vertex weight functions and the α, z's integrals (exactly the prefactor in front of the large bracket in (B8)). Integration over the momentum k (with the shift k + (aP + bq)/2 → k then yields   TABLE 1. Dimensioneless couplingλ = g 2 /(4πm) 2 as a function of fraction of binding η = √ P 2 /2m for two cases of exchanged mass m 3 . The case m 3 /m = 0 is the Wick-Cutkosky model. The second case m 3 /m = 0.5 is compared with the result obtained by Kusaka et al [8].     η=0.2 η=0.5 η=0.9 η=0.99 η=0.999 D. Gross D. BSE (16) D. BSE (32) FIG. 9: The behavior of the elastic electromagnetic form factors for various bound states characterized by η. The mass of exchanged particle is fixed to be m 3 = 0.5m, except for the scalar deuteron case (D.), which is calculated for comparison using two different grids.