{'fraction_non_alphanumeric': 0.03917709358451661, 'fraction_numerical': 0.01691738132058672, 'mean_word_length': 4.120455908818236, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Confinement Phenomenology in the Bethe-Salpeter Equation

We consider the solution of the Bethe-Salpeter equation in Euclidean metric for a qbar-q vector meson in the circumstance where the dressed quark propagators have time-like complex conjugate mass poles. This approximates features encountered in recent QCD modeling via the Dyson-Schwinger equations; the absence of real mass poles simulates quark confinement. The analytic continuation in the total momentum necessary to reach the mass shell for a meson sufficiently heavier than 1 GeV leads to the quark poles being within the integration domain for two variables in the standard approach. Through Feynman integral techniques, we show how the analytic continuation can be implemented in a way suitable for a practical numerical solution. We show that the would-be qbar-q width to the meson generated from one quark pole is exactly cancelled by the effect of the conjugate partner pole; the meson mass remains real and there is no spurious qbar-q production threshold. The ladder kernel we employ is consistent with one-loop perturbative QCD and has a two-parameter infrared structure found to be successful in recent studies of the light SU(3) meson sector.


I. INTRODUCTION
QCD models based on the Dyson-Schwinger equations [DSEs] provide an excellent tool in the study of nonperturbative aspects of hadron properties in QCD [1]. Such models can implement quark and gluon confinement [1,2,3,4], generate dynamical chiral symmetry breaking [5,6], and maintain Poincaré covariance. It is straightforward to implement the correct one-loop renormalization group behavior of QCD [7], and obtain agreement with perturbation theory in the ultraviolet region. Provided that the relevant Ward identities are preserved in the truncation of the DSEs, the corresponding currents are conserved. Axial current conservation induces the Goldstone nature of the pions and kaons [8]; electromagnetic current conservation produces the correct hadronic charge without fine-tuning [9].
Previous work [7,10] has shown this model to provide an efficient description of the light-quark pseudoscalar and vector mesons via two infrared parameters within the rainbow truncation of the DSE for solution of the dressed quark propagators coupled to the ladder approximation for the Bethe-Salpeter equation [BSE]. Furthermore, in impulse approximation, the elastic charge form factors of the pseudoscalars [11] and the electroweak transition form factors of the pseudoscalars and vectors [12,13] are in excellent agreement with data. As a collorary, the strong decays of the vector mesons are well-described in impulse approximation without parameter adjustment [14].
Nonperturbative solutions of the DSEs are almost always implemented in the Euclidean metric for technical reasons and this work is no exception. A resulting complication is that the solution of the BSE for meson bound states requires an analytic continuation in the meson momentum to reach the on-mass-shell point P 2 = −M 2 . This causes the quark p 2 in the BSE to vary throughout a domain of the complex plane bounded by a parabola which widens in proportion to the mass M of the meson. For low mass mesons, such as the π and K mesons, this excursion into the complex plane is simple to handle. Although the difficulties encountered in studies of ground state vector mesons are greater, they too may be overcome directly [10]. For mesons with masses greater than about 1 GeV, these excursions into the complex plane are deep enough to encounter singularities in the rainbow DSE solutions for the quark propagators. As a result, little is known about the implementation of the BSE in this situation or even whether a solution is well-defined. For this reason, DSE-BSE model studies of mesons have for the most part been restricted to light mesons. An exception is the study of Ref. [15] in which the Ansatz adopted for the quark propagators in the non-Euclidean domain precluded singular behavior.
It has long been known that the rainbow approximation to the fermion DSE produces complex conjugate singularities on the timelike half of the complex-momentum plane (ℜ(p 2 ) < 0) in QED [16,17]. The phenomenon has been studied in detail within QED [18,19] and in an Abelian-like model of QCD [20]. In these early studies the singularities were attributed to artifacts of the rainbow approximation. A connection between the generation of complex conjugate singularities and the possible absence of real mass poles was noticed in a QED study [16]. Later studies within QCD [2,3,21] raised suggestions that complex singularities may not be artifacts of the approximation at all, but rather may be a property of the full theory, possibly related to confinement and the absence of real quark mass poles. Beyond the rainbow approximation, very little is known about DSE solutions and their singularity structure. These questions and the influence of approximations have yet to be completely resolved. However they are being addressed; for example, in one simple model calculation, the inclusion of the full dressing of the quark-gluon vertex in the DSE leads to the disappearence of all quark singularities in the complex plane save essential singularities at infinity [2].
It has been well-established by Dyson-Schwinger studies that in order to obtain an acceptable amount of dynamical chiral symmetry breaking in agreement with empirical observations, the effective rainbow DSE kernel must have a significant infrared strength. This in turn tends to move mass poles away from the real timelike p 2 axis, at least for a domain characterized by several GeV. This suggests that quark propagators may be successfully modeled in terms of parametrizations based on entire functions [22] or at least based on functions with singularities away from the real axis such as complexconjugate poles.
If such an approach to modeling QCD is to be extended to includeqq andqQ meson spectra beyond 1 or 2 GeV, then it is necessary to determine how solutions of the BSE are defined when complex singularities in the quark propagators may be encountered. This issue is addressed in the following using a simple parmetrization of the quark propagator. Although the method developed herein is first applied to the BSE, it has applications beyond meson spectroscopy. For example, in production processes at moderate momentum transfers, the use of complex mass pole representations for confined quarks has been suggested as a means to handle the sizable timelike momentum carried by the quarks in these processes [23]. One deals with these issues also in applications of the DSE-BSE approach to electromagnetic elastic [11] and transition form factors [12,13] of light mesons at highmomentum transfer.
To simplify the analysis, we adopt a parameterization of the dressed quark propagator in terms of complex conjugate mass poles. In the standard approach to the meson BSE, two of the four integration variables are integrated numerically and it is during the integration of these variables that poles may be encountered as a result of analytic continuation. In this article, a clear implementation, based on Feynman integral techniques, is introduced to map the four-dimensional integration in the BSE onto a one-dimensional domain where numerical methods of evaluating Cauchy integrals are well-established. We find that when a complex conjugate mass pole parametrization of the quark propagators is employed, the would-beqq decay width of the meson generated from a particular mass pole is exactly cancelled by the contribution of the conjugate partner pole. As a result, the meson mass M remains real valued, even above the pseudo-threshold M 2 > ℜ(4m 2 1 ) where m 1 is the lowest complex mass in the propagator representation.
The primary limitation of the method developed here is that it relies on an explicit parameterization of dressed quark propagators in terms of pairs of complex conjugate mass poles. However the limitation is offset by the finding that such a representation can provide excellent fits to both DSE solutions (as is shown in Sec. V) and also to recent lattice-QCD data [24]. Of course, use of such a phenomenological representation will not exactly preserve the dynamical relation between vertices and propagators necessary to satisfy the Ward Takahashi identities of QCD. In particular, without the axial vector Ward Takahashi identity, the chiral-limit pseudoscalar meson solutions will not obey Goldstone's theorem [8]. As a result, the toy model introduced in Sec. V would not be recommended for studies of light pseudoscalar mesons; chiral-symmetry-preserving DSE approaches have already been well-established [7,10] for that sector. Rather, the present work explores techniques for extending the present DSE-based approach to mesons with masses above 1 GeV where one expects to encounter complex conjugate singularities.
The article is organized as follows. In Sec. II, our notation for the ladder BSE for meson states is introduced, as is the complex mass pole representation of the quark propagators which will be employed. In Sec. III, a simplified BSE integral equation is used to illustrate the technique for analytic continuation of the BSE. The fact that the resulting bound state mass will be real is demonstrated explicitly. In Sec. IV, the method is extended to the general BSE for a vector meson. A simple model for the ladder BSE kernel is introduced in Sec. V and numerical results are presented. A discussion of this approach and outlook is presented in Sec. VI.

II. THE BETHE-SALPETER EQUATION
In the following the Euclidean metric is employed. The scalar product of two four vectors is then a · b = 4 µ=1 a µ b µ , and a 2 > 0 for a spacelike 4-vector a µ . The Dirac γ-matrices are Hermitian γ † µ = γ µ and obey the anti-commutation relation {γ µ , γ ν } = 2δ µν . The dressed quark propagator S(p) and meson Bethe-Salpeter (BS) amplitude Γ(p; P ) are solutions of the renormalized DSE, (1) and the BSE, Here D µν (k) is the renormalized dressed-gluon propagator, Γ i ν (q, p) is the renormalized dressed quark-gluon vertex, Γ(p; P ) is the BS amplitude for a quark-antiquark bound-state meson, and K(q, p; P ) is the renormalized two-particle irreducible qq scattering kernel. We consider equal mass quarks with momenta q ± = q ± P/2 where q is the quark-antiquark relative momentum, and P = q + − q − is the meson momentum which satisfies P 2 = −M 2 where M is the meson mass.
In Eq. (2) the double contraction of Dirac indices has been denoted by ⊙. If the Dirac indices of the elements were reinstated in the BSE and written as parenthetic Roman letters (a, b, . . .), then Eq. (2) would appear as Γ (ab) (p; P ) = Λ q K (ab;cd) (p, q; P )S (de) (q + ) Use of the double contraction operator ⊙ allows us to suppress Dirac indices, and simplify the forms of the equations appearing herein. In Eqs. (1) and (2) Λ k ≡ Λ d 4 k/(2π) 4 denotes a translationally invariant ultraviolet regularization of the momentum space integral with mass-scale Λ. Lorentz covariance entails that the solution of Eq. (1) is of the form S(p) −1 = i/ pA(p 2 )+ B(p 2 ). The Lorenz scalars A(p 2 ) and B(p 2 ) are renormalized at spacelike p 2 = µ 2 such that A(µ 2 ) = 1 and B(µ 2 ) = m q (µ), where m q (µ) is the renormalized current quark mass. After renormalization, one is free to remove the regularization by taking the limit Λ → ∞.
Recently, a successful model has been developed that provides an excellent description of the masses, decays and other properties of the light pseudoscalar and vector mesons [7,10]. The model consists of the rainbow-ladder truncations for the quark DSE (1) and meson BSE (2) and the use of an effectiveqq interaction constrained to coincide with perturbative QCD in the ultraviolet domain and containing a phenomenological infrared behavior. In this particular model, the DSE kernel in rainbow truncation is where D free µν (k = p − q) is the free gluon propagator in Landau gauge, and α eff (k 2 ) is the effective running coupling. The ladder truncation of the BSE kernel is where k = p − q and ⊗ denotes the direct (or Cartesian) product of matrices. This truncation scheme is self-consistent in that it ensures that the dressing of the quark-antiquark vector and axial-vector vertices generated by Eqs. (4) and (5) satisfy their respective Ward-Takahashi identities. This feature is important since the axial-vector Ward-Takahashi identity guarantees that, in the chiral limit of zero current quark mass, the ground state pseudoscalar meson bound states are the massless Goldstone bosons arising from dynamical chiral symmetry breaking in accordance with the Goldstone's theorem [7,8]. The vector Ward-Takahashi identity guarantees electromagnetic current conservation for mesons and nucleons if the impulse approximation is used to describe the current in terms of quarks [9,11]. The rainbow-ladder truncation is particularly suitable for the flavor octets of light pseudoscalar and vector mesons where higher-order contributions to the quark-gluon skeleton-graph expansion have significant cancellations [25,26] and so may be safely neglected.
With the ladder truncation from Eq. (5), the BSE for a light-quark vector meson is with k = p − q and Here Γ ρ is a transverse four vector. We will refer to the vector meson under study here as the ρ meson although, due to a number of simplifications made for illustration purposes, the particular model solutions considered here should not be taken as physical representations of the ρ meson.
A linear eigenvalue λ(P 2 ) has been introduced into Eq. (6) so that it will yield solutions over a continuous range of P 2 . A physical bound state of mass M ρ corresponds to an eigenvalue of λ(P 2 = −M 2 ρ ) = 1. To find such a solution, it is clear that Eq. (6) must be analytically continued from Euclidean space where P 2 > 0 to Minkowski space where P 2 < 0 and the on-massshell condition P 2 = −M 2 ρ may be realized. In terms of the Euclidean four momentum P µ , such an analytic continuation would correspond to P 4 = i P 2 + M 2 ρ , from which it follows that quark momenta q 2 ± = q 2 −M 2 ρ /4 ± i q 2 M ρ z that enter the BSE are complex functions of the integration variables q 2 > 0 and the direction cosine −1 ≤ z ≤ +1. Therefore, in order to carryout such an analytic continuation of the BSE (6), the dressed quark propagators S(q 2 ± ) must be known in the complex-momentum plane within a parabolic region having the negative real (timelike) point −M 2 ρ /4 as the apex and extending symmetrically about the real axis.
Studies within the rainbow-ladder truncation [10] find that the non-analytic points of u-and d-quark propagators that occur nearest to the origin q 2 = 0 are located at q 2 = −0.207 ± i0.331 GeV 2 while those of the s-quark propagator are found at q 2 = −0.376 ± i0.602 GeV 2 [14]. This finding also agrees with qualitative observations from an earlier study [27]. From the positions of these non-analytic points one can define a critical mass M c = √ −P 2 > 0 above which these singular points enter the domain of integration for the BSE. For the model of Ref. [10], these critical masses are M c = 1.09 GeV forūu mesons and M c = 1.47 GeV forss mesons. Clearly, the non-analytic points lie outside the domains of integration for ladder BSE calculations of ground state pseudoscalar and vector mesons, but may lie within the integration domain for heavier mesons. This observation provides one of the main reasons that many BSE studies of mesons focus primarily on light meson states.
The appearence of poles in the domain of integration of the BSE raises two questions: 1) does the analytic continuation of the BSE to timelike momenta P 2 < 0 yield physically meaningful quantities, and if so, 2) by what methods does one obtain the solutions of this singular integral equation? The goal of this article is to explore these questions. To facilitate an initial investigation, it is convenient to adopt a parametrization of the quark propagator S(p) as a sum of N pairs of complex conjugate mass poles, where m n are complex-valued mass scales and z n are complex coefficients. Comparison with the general form S(p) −1 = i/ pA(p 2 )+ B(p 2 ), allows identification of the amplitudes A(p 2 ), B(p 2 ), and the quark running-mass allows identification of the scalar and vector amplitudes of the propagator, σ S (p 2 ) and σ V (p 2 ) respectively. These are needed for the BSE. For example, the quark propagator of Eq. (8) has a scalar amplitude of the form and a vector amplitude σ V (p 2 ) which is obtained by removing each of the m n and m * n from the numerators. It is shown in Sec. V that this form of quark propagator with N = 3 provides an excellent fit to DSE solutions that employ realistic effective qq interactions. Furthermore, this form has been used recently to provide a parametric fit of the chiral-quark propagator [24] obtained from lattice-QCD simulations [28,29].

III. ANALYTIC CONTINUATION: A SIMPLE EXAMPLE
A general study of the analytic continuation of the BSE is quite involved and complicated by the Dirac algebra necessary to describe quark-antiquark states. It is prudent to first proceed by simplifying several of these complications and focus on the issue of analytic continuation of the eigenvalue λ(P 2 ) and the associated integrations over the complex-conjugate poles in the quark propagator. These restrictions are lifted in Sec. IV where the approach introduced here is extended to the general BSE.
In the following, assume that the quark propagator is well-represented by a single pair of complex conjugate poles (that is, N = 1 in Eq. (8)). Furthermore, assume that the vector-meson BS amplitude may be represented by the single Dirac amplitude Γ ρ (p; P ) = (δ ρα − P ρ P α /P 2 )γ α V (q; P ), where P ρ is the four momentum of the bound-state vector meson. After these simplifications, the BSE (6) may be converted to a scalar eigenvalue equation for the amplitude V (q; P ) by multiplication from the left by γ ρ followed by a trace over the Dirac indices to obtain . (10) Observe that the kernel of this equation is bilinear in the propagator amplitudes σ S (q 2 ± ) and σ V (q 2 ± ) but only the denominators of these functions are displayed explicity in Eq. (10). The remaining factors from the quark propagators and effective qq interaction have been collected into the scalar function f ab whose detailed form is unimportant to the present discussion but may be derived following the discussion of Sec. IV. In addition, the denominators have been written in terms of µ 2 = m 2 R − m 2 I and ∆ 2 = 2m R m I , where m = m R + im I is the complex mass scale parameter in the quark propagator in Eq. (9). This form allows a parametrization of the square of the quark mass m 2 and its complex conjugate m 2 * in terms of the values a = ±1.
In principle, the q-dependence of the integrand in Eq. (10) is only known once the solution V (q; P ) has been obtained. For the present, consider the artificial situation in which the solution of the BSE is known to be of the form V = g φ(q 2 ) U (z q ) where z q = q·P/ q 2 P 2 and φ(q 2 ) and U (z) are known functions. In this case, the constant g may be eliminated from the BSE and following an integration over d 4 p one obtains, and where F ab (q, P ) is a non-singular, Lorentz invariant function of q and P . The objective is then to evaluate the four integrals λ ab (P 2 ) and show how they may be analytical continued to negative values of P 2 . This is carried out explicitly in the remainder of this section.
When P 2 is continued to negative values, the quark momenta in Eq. (12) satisfy q 2 ± = (q 2 ∓ ) * . The two factors in the denominator of λ a,−a (P 2 ) (that is, a = −b in Eq. (12)) will thus be complex conjugates. As a result, for P 2 sufficiently timelike, there will always be some value of the integration variable q µ for which both factors in the denominator of λ a,−a (P 2 ) vanish simultaneously. The result is a non-analytic point in λ +− (P 2 ) and λ −+ (P 2 ). Such behavior would usually be indentified with the decay threshold for the vector meson into an asymptotically-free quark and antiquark pair. However, in the following, it is shown that when the integrals in Eq. (12) are added together the imaginary parts cancel exactly thereby eliminating the possibility of quark-antiquark production thresholds. In contrast to this, the two factors in the denominator of λ aa (P 2 ) never vanish simultaneously and consequently there are no corresponding non-analytic points in either λ ++ (P 2 ) or λ −− (P 2 ).
To analyze the behavior in detail, it is advantageous to restructure the integral in Eq. (12) so that the analytic continuation from the Euclidean domain P 2 > 0 to the timelike domain P 2 < 0 is straightforward. For values of P 2 > 0, one may combine the two factors in the denominator using the Feynman parametrization, and obtain where y = α − β, and (15) One then eliminates the term linear in q in the denominator by a shift of the origin q ν → q ν − yP ν /2. The resulting denominator depends on P 2 , q 2 , and y but not on q · P . After this shift, the function F ab (q, P ) in the numerator still depends q · P . All other variables save q 2 and y can be integrated out yielding, where q 2 ab (P 2 , y) = − is independent of q 2 , and F ab (q 2 , y; Here z q is the direction cosine defined by q · P = q 2 P 2 z q . The singular behavior due to the propagators is now a single second-order pole in the variable q 2 . Its location q 2 ab (P 2 , y) is given in Eq. (17) and is not in the path of integration for spacelike P 2 > 0; the integral is regular in Euclidean space.
The y integration is performed after the integration over q 2 . Hence, for a given value of y, the q 2 -integration domain is the positive real axis and the integrand has a pole in the complex q 2 plane whose location q 2 ab (P 2 , y) moves with P 2 . The location has a negative real part for P 2 > 0 and remains negative as P 2 passes through zero and becomes moderately timelike. Only when  17) is shown as the dashed curve. The solid curve is the contour required to implement analytic continuation. P 2 < −4µ 2 will the real part of the pole location become positive and have a chance of impinging upon the domain of the q 2 -integration. For that to happen, the imaginary part of q 2 ab (y) would have to vanish. As anticipated above, Eq. (17) shows that this never happens when a = b; for λ ++ (P 2 ) and λ −− (P 2 ) the two factors in the denominator of Eq. (12) are never zero simultaneously. It is therefore trivial to carryout these integrations. However, for a = −b, the imaginary part of q 2 ab (y) can become zero when the Feynman variable y = 0 and care must be taken to properly handle such singular integrals.
As a specific example, the trajectory of the pole location q 2 +− (y) is shown in Fig. 1. It begins at (−µ 2 , ∆ 2 ) when y = −1, moves along the dashed curve as y increases, crossing the real axis when y = 0, and terminates at (−µ 2 , −∆ 2 ) when y = 1. For y < 0, the calculation of the integral is straightforward. However for y > 0, the pole in q 2 has crossed the real axis and the analytic continuation of the integral λ +− (P 2 ) is defined by continuous deformation of the contour in a manner that ensures the pole does not cross the contour of integration. This is illustrated in Fig. 1, for the deformed contour corresponding to y ≈ +1/2. Deformation of a integration contour in this manner is equivalent to adding the pole residue to the result of integration along the real q 2 axis. If the residue term were neglected, one would obtain a discontinuous result due to the crossing of a branch cut onto the wrong sheet.
For the integral λ −+ (P 2 ), the pole trajectory traces out the same path as in Fig. 1 but in the opposite direction; that is, y > 0 corresponds to the upper half q 2plane and the pole residue will be of opposite sign from λ +− (P 2 ). After analytic continuation to P 2 < 0, the results for the four integrals can be expressed as where Θ(x) is the Heaviside step function.
In the timelike region, inspection of Eq. (19) shows that individual integrals λ ab (P 2 ) are complex valued and satisfy The imaginary parts cancel in pairs thereby producing a real result for the eigenvalue λ(P 2 ) = a,b=± λ ab (P 2 ). Notice that this does not mean that one may neglect the addition of residue terms required by the analytic continuation of the integrals. Rather, when ∆ 2 > 0 the residue (second) term of Eq. (19) is not purely imaginary and the sum of the real parts of the residues are required for the final result for λ(P 2 ). Of course, these residue terms appear only for P 2 more timelike than the pseudo-threshold point P 2 = −4µ 2 , that is, for bound state masses M ρ > 2µ. Before extension of this treatment to the general problem, consider the connection between the above discussion of analytic continuation and the usual case encountered above physical two-body thresholds. In particular, consider the limiting case of infinitesimal ∆ for the single integral λ −− (P 2 ). This corresponds to the standard case of employing elementary propagators with the conventional boundary conditions of a negative infinitesimal −iǫ added to the real mass. Above the two-body threshold P 2 < −4µ 2 , there would beqq production described by the residue term in λ −− (P 2 ). It follows that a bound state calculation, based on such an integral alone, would result in a complex mass describing an unstable meson with a finiteqq decay width. However, in the present case the presence of conjugate quark poles ensures that the analytic continuation introduces both positive and negative would-be widths that cancel exactly. The resulting BSE eigenvalue λ(P 2 ) in Eq. (10) will be purely real and therefore bound state masses defined by λ(P 2 = −M 2 ρ ) = 1 will be real and will describe a meson with noqq decay width.

IV. GENERALIZATION TO THE BETHE-SALPETER EQUATION
The technique described in the previous section is generalized to allow solutions to the homogeneous BSE for an arbitrary quark-antiquark bound state. To be specific, we consider an equal quark-antiquark vector meson bound state. One can write the most-general BS amplitude for a massive vector meson as a sum of eight Dirac covariants and Lorentz-scalar amplitudes F α (q; P ), Γ ρ (q; P ) = δ ρµ − P ρ P µ P 2 γ µ F 1 + q µ q F 2 + q µ P F 3 + γ 5 ǫ µανβ γ α q ν P β F 4 + q µ F 5 + σ µν q ν F 6 + σ µν P ν F 7 + q µ σ αβ q α P β F 8 .
In a more concise form this is where Λ α are the Dirac covariants in Eq. (21); for example, Λ α=1 ρ = (δ ρµ − P ρ P µ /P 2 ) γ µ . When all indices of a product of two Dirac covariants are contracted in the manner corresponding to the usual trace operation, one obtains where t βα (q, P ) is a Lorentz scalar function of q 2 , P 2 and q · P , and we have used the double-contraction product ⊙ notation introduced in Eq. (2). Of course, one is free to choose any other set of Dirac covariants Λ α ρ (q, P ) as long as they are complete. In particular, one is free to choose a set of covariants for which t αβ → δ αβ in which case the resulting set of coupled integral equations may be of a simpler form than those obtained herein. However, in practice the behavior of such covariants can be poorly behaved near q = 0, and can lead to numerical difficulties [10]. Therefore, the set of covariants given by Eq. (21) is preferable and the lack of orthogonality is easily accommodated by t αβ (q, P ).
In the previous section, it was shown that when the q 2 and q · P dependence of the BS amplitude Γ µ (q; P ) is known, the analytic continuation of the BSE may be carried out in a straight-forward manner. In practice, the functional form of the BS amplitude is known only after the solution to the BSE has been obtained. The approach outlined in the previous section may be generalized by expansion of the functional dependence of the BS amplitudes F α (q; P ) in terms of complete sets of known functions of q 2 and z q = q · P/ q 2 P 2 . Consider the expansion where U i (z q ) are Chebyshev polynomials of the second kind satisfying the ortho-normality condition and φ s (q 2 ) are a set of ortho-normal functions, that are specified later, and that satisfy The required number of basis terms will be determined by convergence of the BSE solution with respect to the meson observables of interest. For example, one particular choice of basis {φ s (q 2 )} which proves to be efficient for observables dominated by infrared physics may be inefficient for observables sensitive to the ultraviolet. In the present work, we take convergence of the meson mass as the criterion. For brevity, we define the operator O si (p) which projects out the expansion coefficients g α si (P 2 ) from the BS amplitude in Eq. (24). Using Eqs. (22) and (23), one finds The set of coefficients g α si (P 2 ) will be refered to as BS amplitudes since they are equivalent to the BS amplitude Γ ρ (p; P ) in that they contain all the dynamical information regarding the solution of the BSE. The expansion coefficients g α si (P 2 ) are functions of P 2 only. The operator O si (p) is an integral operator that is introduced for brevity and notational convenience.
We take P µ = √ P 2 (0, 0, 0, 1), p µ = p 2 (0, 0, 1 − z 2 p , z p ), and the integration momentum q µ is represented as The integrand in Eq. (28) has no dependence upon angle φ; the only dependence on direction cosine x q is through the p · q dependence of the BS kernel K(p, q; P ). In Eq. (28) we may therefore substitute The integration over x q in Eq. (28) produces a function of p 2 , q 2 , z p , z q and P 2 which can be expanded in the basis functions. That is, we define a quantity V by +1 −1 Upon substituting this back into Eq. (28) above, one obtains a discretized form of the BSE, λ(P 2 ) g α si (P 2 ) = V αA;βB si;wℓ (P 2 ) G AB wℓ,tj (P 2 ) g β tj (P 2 ) , (32) where repeated indices are understood as being summed over. Here we have introduced the projection of the product of propagator amplitudes in the form Thus far, the BSE has been reduced to the discrete eigenvalue problem given in Eq. (32). Physical solutions are at λ(P 2 = −M 2 ρ ) = 1 from which one identifies the existence of a bound state vector meson of mass M ρ . The only task remaining is the determination of the elements in Eq. (32) that make up the kernel of the eigenvalue problem. These are V and G. The generalized interaction V is determined in terms of the kernel by the projection in Eq. (31).
Consider the explicit form of G. Clearly there is a strong similarity between the form of Eq. (33) and the simple integral whose analytic continuation was considered in Sec. III. This is clarified by the expansion of the quark propagator amplitudes in terms of the mass pole terms from Eq. (9) in the form where Z A a=+1 n = z n or z n m n for A = V or S, respectively, and Z A a=−1 n = z * n or z * n m * n , respectively. Then, Eq. (33) can be written as .

(35)
These integrals in Eq. (35) are of the type λ ab (P 2 ) already explored in our discussion of analytic continuation in Sec. III. The only difference is that the number of integrals that have to be performed can be large since it depends on the number of basis terms for representation of the BS amplitudes via Eq. (24), and the number of conjugate pairs of mass poles for representation of the propagators.
The same methods employed in Sec. III may be used here to calculate the integrals G AB st;ij (P 2 ) for all values of P 2 while accounting properly for the analytic continuation. In particular, the Feynman integral method used here will lead to a denominator [q 2 − q 2 nmab (P 2 , y)] 2 , where the pole location is The integration over q 2 is carried out first, followed by the integration over the Feynman parameter y, and finally the integration over z q . Integrations over the two remaining angles in Eq. (35) are trivial. In general, for a quark propagator parametrized in terms of N pairs of complex-conjugate poles, there are (2N ) 2 pole trajectories, half of which impinge on the q 2 integration domain and lead to residue additions. The imaginary parts of the integrals in Eq. (35) cancel in pairs as described in the previous section.
In ladder approximation, the determination of V is straight-forward. There are never any poles encountered when P 2 is analytically continued to negative values. Of course, singularities may be encountered when contributions beyond ladder are maintained in the BS kernel, and these must be treated in a way analogous to G. This possibility is not addressed in the present work.

A. Ladder-Rainbow Kernel
Here we use a particular model for the DSE-BSE ladderainbow kernel to numerically implement the preceeding developments. The model is specified by the "effective coupling" α eff (k 2 ) and we employ the Ansatz [10] 4πα eff (k 2 ) with γ m = 12/(33 − 2N f ) and F (s) = (1 − exp −s )/s. The ultraviolet behavior matches that of the 1-loop QCD running coupling α s (k 2 ). The resulting solution of the ladder-rainbow DSE-BSE system of equations in the UV region generates the correct 1-loop perturbative QCD structure. The first term implements the strong infrared enhancement in the region 0 < k 2 < 1 GeV 2 phenomenologically required [30] to produce a realistic value for the chiral condensate. With m t = 0.5 GeV, τ = e 2 − 1, N f = 4, Λ QCD = 0.234 GeV, and a renormalization scale µ = 19 GeV, it has been found that ω = 0.4 GeV and D = 0.93 GeV 2 give a good description of qq , m π/K and f π with physically acceptable current quark masses, [7,10]. The propagator amplitudes from this rainbow DSE model have recently been shown [14] to have the same qualitative features as lattice QCD simulations [28,29].
The vector meson masses and electroweak decay constants produced by this model are in good agreement with experiments [10]. Without any readjustment of the parameters, this model agrees remarkably well with the most recent Jlab data [31] for the pion charge form factor F π (Q 2 ). Also the kaon charge radii and electromagnetic form factors are well described [11,32]. The strong decays of the vector mesons into a pair of pseudoscalar mesons are also well-described within this model [14,33].

B. Results
The DSE solution for the u/d quark propagator from the model interaction given in Eq. (37) can be well fit with N = 3 pairs of complex conjugate poles in the representation given in Eq. (8). The result displayed in Fig. 2 corresponds to the parameter set This parameterization was chosen so that Z(p 2 ) and M (p 2 ) were reproduced well at both p 2 = 0 and p 2 = 100 GeV 2 while the main features of the momentum dependence are preserved. We imposed the constraint that M (p 2 ) should approach its UV limit from above. The pole locations are p 2 = −m 2 i where m 2 i = µ 2 i + i∆ 2 i . Thus we have µ 1 = 0.455 GeV, ∆ 1 = 0.576 GeV, µ 2 = 1.13 GeV, ∆ 2 = 1.12 GeV, and µ 3 = 1.46 GeV, ∆ 3 = 1.33 GeV.
With use of this propagator representation in the BSE, a u/d quark meson bound state with mass greater than the lowest pseudo-threshold 2µ 1 =0.910 GeV would be needed to test the method we have described. The model-exact ρ meson mass from the kernel Eq. (37) is 0.742 GeV [10] and is therefore not suitable. Heavier u/d mesons have not been carefully studied within this model interaction due to the very singularity issue we are concerned with. In fact the evidence from exploratory studies is that, for example, m a1 in this ladder-rainbow model is only about 200 MeV above m ρ and also does not provide a clear test. To test the present method, we artificially modify the pole parameterization of the quark propagator so that the lowest pseudo-threshold moves below the vector meson mass, and the latter moves to 0.770 GeV. We achieve this with the following parameterization in terms of N = 2 pairs of conjugate mass poles m 1 = −0.400 + i0.200 GeV, z 1 = 0.15, m 2 = 0.550 + i0.350 GeV, z 2 = 0.35 + i0.37 .
The resulting propagator amplitudes are also shown in Fig. 2. This parameterization corresponds to µ 1 = 0.346 GeV, ∆ 1 = 0.4 GeV, and µ 2 = 0.424 GeV, ∆ 2 = 0.621 GeV. With this propagator, the lowest pseudothreshold in the BSE is at 2µ 1 =0.692 GeV. The subsequent BSE solution is not to be taken as a physical represenation of the ρ meson; its purpose is to provide a test calculation where the propagator singularities are within the domain of integation. The basis functions used to expand the q 2 dependence of the BS amplitudes are where w = q 2 /Λ 2 , the typical hadron scale is Λ = 0.800 GeV, L s (w) are Laguere polynomials, and N s is fixed by the normalization condition given in Eq. (26). To obtain better than 1% accuracy for the eigenvalue λ(P 2 ), we found that 7 terms in this basis were sufficient. The number of Chebyshev terms used for the angle basis was one. The first covariant from the set given in Eq. (21) is known to be dominant for the vector meson solution [10] and this was the only one used here.
The resulting BSE eigenvalue λ(P 2 ) is given in Fig. 3 as the solid curve. In the spacelike region (P 2 > 0) the eigenvalue is a monotonically, decreasing function of P 2 . When P 2 decreases below zero and becomes timelike, the eigenvalue λ(P 2 ) increases for a while then crosses unity, at P 2 = −M 2 ρ , signifiying a bound state at M ρ = 0.770 GeV. In Fig. 3 this value of P 2 is shown as a thin, vertical line. On the timelike side of this line, the eigenvalue continues to increase, then begins to decrease and ultimately goes to zero for (P 2 ≪ 0). This is because the eigenvalue is a measure of the magnitude of the BSE kernel and the propagators therein falloff as the difference between the momentum arguments and the pole positions becomes larger than any mass scale in the interaction.
The solid curve in Fig. 3 represents the proper analytic continuation of the BSE from P 2 > 0 to P 2 < 0. The dashed line represents the result from a direct evaluation of the BSE integral without regard to the possibility of singularities and their movement as a function of P 2 in relation to the integration domain. That is, only a transcription from P 2 > 0 to P 2 < 0 has been made in the integrand of Eq. (35); the Feynman technique for combining the denominators and shifting the 4-momentum variable to complete the square has not been made. The integration variables are q 2 and z q =q ·P . It is clear that for the spacelike domain (P 2 > 0), and for the limited timelike domain for which no poles move into the domain of integration (−4µ 2 1 < P 2 < 0, i.e., below the lowest pseudo-threshold), the two approaches give identical results. The dashed line represents the approach to the BSE that has usually been implemented; it has been limited to states light enough to avoid the propagator singularities.
With the parameters given above, the quark poles are first encountered within the BSE integration when the analytic continuation reaches P 2 = −4µ 2 1 , that is, at P 2 = −0.48 GeV 2 . Although the contibutions of the conjugate pair of poles to the imaginary part of the BS eigenvalue λ(P 2 ) cancel, their real contribution produces a discontinuity in the derivative λ ′ (P 2 ) at that point. Careful inspection of the solid curve in Fig. 3 reveals that this cusp occurs at the precise value where the full calculation and the naïve calculation diverge.

VI. DISCUSSION
We have provided a method to obtain meson bound state solutions from the Bethe-Salpeter equation in Euclidean metric when the dressed quark propagators have time-like complex conjugate mass poles within the integration domain. This approximates features encountered in recent QCD modeling via the Dyson-Schwinger equations; the absence of real mass poles simulates quark confinement. For this exploratory study we represent the quark propagators as a sum of complex-conjugate mass poles, and project the BSE on to complete basis sets to represent the momentum and angle dependence. We use Feynman integral techniques to combine the propagator denominators and map the integration to a onedimensional domain. This allows a clear analysis of the analytic continuation in total meson momentum needed to reach the mass shell and which causes the singularities to impinge upon the integration domain. The BSE linear eigenvalue remains real; in other words the eigenmass remains real. The would-be decay width from one pole is exactly cancelled by the effect of the partner pole; the meson is stable against decay into a quark-antiquark pair. This describes the confinement of quark and antiquark within the meson bound state, even though the meson mass is "above threshold".
One of the limitations of the method presented here is that it relies upon the projection of the BSE on to a complete basis set of functions to represent the momentum dependence of the BS amplitudes. In general the behavior of the BS amplitudes is not known until after the BSE is solved. If the chosen basis set is an inefficient representation, one would expect this to result in a lack of convergence. The (vector) meson mass that we have studied here is an integrated quantity dominated by infrared physics; it is not surprizing that convergence is easily achieved with the basis states of Eq. (40) which fall off exponentially in the ultraviolet region. However, we know from studies such as Ref. [10] that the leading ultraviolet power law behavior of the vector meson BS amplitudes is 1/q 2 . Clearly the present basis would be inefficient for applications that depend strongly on the BS amplitude in this domain, such as the asymptotic behavior of the pion charge form factor [34]. For such studies, appropriate basis sets would have to be utilized.
The principal reason that the BS amplitudes are expanded in terms of known basis functions is that the Feynman integral technique for handling the two propagator denominators requires a complex variable shift for the rest of the integrand. In principle, the analytic continuation of the BS amplitude into the complex plane is determined by the dynamics of the BSE and can be determined only through the solution. The implementation of the present method of solution requires a known analytic behavior for the BS amplitudes as expressed through the basis functions. One would expect an inadequacy of the basis in this respect to show up as poor convergence.
Future work using the BSE solution method presented here will include extension of the DSE modeling approach to SU(3) flavor meson states above 1 GeV and meson form factors for Q 2 > 3 GeV 2 .