Large-j Expansion Method for Two-Body Dirac Equation

By using symmetry properties, the two-body Dirac equation in coordinate representation is reduced to the coupled pair of radial second-order differential equations. Then the large-j expansion technique is used to solve a bound state problem. Linear-plus-Coulomb potentials of different spin structure are examined in order to describe the asymptotic degeneracy and fine splitting of light meson spectra.


Introduction
Two-body Dirac equations (2BDE), i.e., the Breit equations [1] and its generalizations [2,3,4,5,6,7,8,9,10], are used frequently for the description of relativistic bound state problem, especially in nuclear [11,12] and hadronic [13,14,15,16] physics. Apart from two free-particle Dirac terms, the 2BDE may include potentials which are local matrix-functions in the coordinate representation. This form provides an intuitive understanding of the interaction and may suggest a proper physical choice of the potential in phenomenological models.
But the 2BDE are pathological if certain interaction terms are not treated perturbationally. The set of radially reduced equations [13,5,7] may possess non-physical energy-dependent poles at finite distance r between particles [17,18]. Correspondingly, an exact boundary-value problem becomes incorrect mathematically or improper for assumed physical treatment.
Here we consider a possibility to avoid pathological peculiarities of 2BDE using a pseudoperturbative technique similar to 1/N [19,20,21] or 1/ℓ expansions [22,23]. These methods are applicable to the case of a strong coupling and are little affected by boundary peculiarities of the boundary-value problem.
In our case natural expansion parameter is 1/j, where j is the conserved total angular momentum. After the radial reduction is performed, the 2BDE takes the form of the set of eight coupled first-order differential equations [5,7]. Using a chain of transformations we reduce it to the pair of coupled second-order equations and apply the 1/j expansion technique. The method is applied to the potential model of meson based on the 2BDE.

2-body Dirac equation and its radial reduction
In the centre-of-mass reference frame the 2BDE has the form: where F (r) is a 16-component wave function, h a (p) = α a · p + m a β a ≡ − i α a · ∇ + m a β a , a = 1, 2, are Dirac Hamiltonians of free fermions of mass m a and U (r) is an interaction potential. If F (r) is presented in 4 × 4-matrix representation, the operators α a and β a act as follows: α 1 F = αF , α 2 F = F α T etc, where α and β are Dirac matrices. The potential U (r) is the Hermitian matrix-function, it is invariant under rotation and space inversion transformations (so that the total Hamiltonian H = h 1 + h 2 + U is too). Its general form is parametrized by 48 scalar function of r = |r| [10]. Of physical meaning are potentials admitting field-theoretical interpretation of interaction. In particular, potentials reflecting a spin structure of vector and scalar relativistic interactions are used frequently in potential quark models of mesons. We will consider such a model in the Section 5 using few examples of scalar and vector potentials known in a literature. In the present section the structure of potential is not essential.
In order to apply a pseudoperturbative expansion method to the 2-body Dirac equation let us transform it to an appropriate form.
First of all we perform a radial reduction. Following [5,7] we put the wave eigenfunction F (r) of the total angular momentum j and the parity P into the 2 × 2 block-matrix form: for the parity P = (−) j±1 states, and into a similar form for the parity P = (−) j states but with superscripts interchanged as follows: (A, 0) ↔ (−, +). Here n = r/r, the bispinor harmonics φ A (n) corresponds to a singlet state with a total spin s = 0 and an orbital momentum ℓ = j, and φ 0 (n), φ − (n), φ + (n) correspond to triplet with s = 1 and ℓ = j, j + 1, j − 1. Then for j > 0 the eigenstate problem (1) reduces to the set of eight first-order differential equations with the functions s 1 (r), . . . , v 2 (r) and the energy E to be found. It is convenient to present this set in the following matrix form. Let us introduce the 8dimensional vector-function: X(r) = {s 1 (r), s 2 (r), t 1 (r), . . . , v 2 (r)}. Then the set of radial equations reads: where the 8×8 real matrices H(j) and V(r, E, j) = G(j)/r +m+U(r, j)−EI possesses properties unit matrix and m ± = m 1 ± m 2 ) and j-and P -dependent matrices H(j), G(j) are constant (i.e., free of r), and matrix-potential U(r, j) comes from interacting term of the equation (1). For the case j = 0 components s 2 = t 2 = u 2 = v 2 = 0 so that the dimension of the problem (3) reduces from 8 to 4. It turns out that rank H = 4 (2 for j = 0). In other words, only four equations of the set (3) are differential while remaining ones are algebraic. They can be split by means of some orthogonal (i.e., of O(8) group) transformation. In new basis we have where J (2) is the symplectic 4 × 4 matrix. Thus we arrive at the set Eliminating X 2 from (4) by means of (5) we get a differential set for the 4-vector X 1 while X 2 then follows from the algebraic relation X 2 = −V −1 22 V 21 X 1 . The elimination of X 2 causes non-physical energy-dependent singular points (apart of r = 0 and physical singularities of potentials) in matrix elements of V ⊥ . Now we present the 4-vector X 1 in 2 + 2 block form, eliminate then Φ 2 and arrive at the second-order differential equations for 2-vector Φ 1 : The matrix V 22 is diagonal for all potentials considered in Section 5 (and many other ones). In these cases we can perform the transformation: providing for the operatorL the form which is as close as possible to 2-term form: here W(r, E, j) is a symmetric 2 × 2 matrix, J (1) is 2 × 2 symplectic matrix and {·, ·} + denotes the anticommutator.
We are going to apply the 1/j expansion method to the equation (7). In many physically interesting cases the function Z vanishes or it is negligible at j large. Thus the wave equation has a 2 × 2 matrix 2-term form which is convenient for application of the method. In other cases the third term of the operator (7) contains a first-order derivative via off-diagonal matrix elements only. This form is tractable too, but with more tedious calculations. We do not consider such equations in this paper. Before proceeding further, we study a simpler example of a single 2-term relativistic equation.

Todorov equation via 1/ℓ method
Here we consider the Todorov-type equation describing the relativistic system of two interacting scalar particles in the centre-of-mass reference frame [24,25,26]: Here p = − i ∇, the quasipotential U (r, E) depends on energy E of the system, and the binding parameter b(E) is the following function of E, The corresponding radial equation takes the form where ℓ is the angular momentum quantum number, and Let us consider motion of the system in the neighbourhood of classical stable circular orbit. Given ℓ > 0, the radius r c = r c (ℓ) of the stable circular orbit and the corresponding energy E c = E c (ℓ) satisfy conditions: and One puts r = r c + ∆r and E = E c + ∆E where ∆r and ∆E are small in some meaning, and expand the function W (r c + ∆r, E c + ∆E, ℓ) in power series with respect ∆r and ∆E. Then due to the conditions (11) the leading terms of this expansion represent the harmonic oscillator potential and other ones are anharmonic terms. If the conditions (11) hold for any large value of ℓ it is possible by renormalization of ∆r and ∆E to single out in the equation (9) the ℓ-independent harmonic oscillator problem and anharmonic perturbations which disappear if ℓ → ∞. This is the idea of 1/ℓ expansion method. Application of pseudoperturbative techniques of this type [19,20,21,22,23] to our case meets two peculiarities: the equation (9) represents a nonlinear spectral problem, and an exact solution of the equations (11) may appear to be unknown or too cumbersome for practical use. Thus we modify the technique.
Let us introduce the parameter λ = 1/ √ ℓ which is small at ℓ large. Since the exact form of the functions r c (ℓ) and E c (ℓ) is unknown in general, we first determine asymptotics r c ∼ r ∞ (λ), b c = b(E c ) ∼ b ∞ (λ) at λ → 0 which may be found much easier. Then the functions r c (ℓ) and E c (ℓ) can be presented in the form: where expansion coefficients ρ (n) , µ (n) , n = 1, 2, . . . (and thus the analytical functions ρ(λ) and µ(λ)) can be found, step by step, from the conditions: and ∂ 2W (ρ, µ, λ)/∂ρ 2 > 0; here the dimensionless functionW (ρ, µ, λ) is constructed by the direct use of (12) in (10) and normalizing in order thatW (ρ, µ, λ) to be regular at λ → 0, Now we go to the dimensionless variable r → ξ and spectral parameter b(E) → ǫ, in terms of which the equation (9) takes the form If the functions ρ(λ) and µ(λ) satisfy the conditions (13), the equation (15) is nonsingular at λ → 0. This is true even if we use the first-order approximate solution to (13) in (14), Indeed, using the notation Singular terms are absent if the following set of equations holds: Besides, zero-order terms which are linear in ξ disappear if Notice that the equations (19) and (20)-(21) represent the conditions (13) in the zeroth and first orders of λ, respectively. Thus the equations (19) hold identically and (20)- (21) are linear equations with ρ (1) and µ (1) to be found. In zero-order approximation the equation (15) reduces to the harmonic oscillator problem The higher-order terms in the expansion (18) can be considered as perturbations of the oscillator problem (22). They depend, in general, on the spectral parameter ǫ and can be taken into account by means of the perturbative procedure [25] is appropriate to this case. Otherwise the treatment is similar to [19,20,21,22,23]. The eigenvalues in zero-order approximation ǫ nr = [ω(2n r + 1) + ν]/κ, where n r = 0, 1, . . . is a radial quantum number, are to be corrected by means of higher orders of perturbative procedure. Then, using of the 2nd equation of (14) in (8) gives us the energy spectrum. 4 Breit-type equation via 1/j method At this point we return to the radial 2BDE in the formL(E)Φ 1 = 0, where the 2 × 2 matrix operatorL(E) is given by equation (7) with the last term neglected. Let us put where Ψ 1 and Ψ 2 are components of Φ 1 . Then the equation (6) can be presented as a pair of coupled ordinary second-order differential equations: We will treat this system perturbationally using the pseudosmall parameter λ = 1/ √ j.
Let us suppose for a moment that the right-hand side of the system (26)- (27) can be ignored, so that these equations decouple. Then we can apply to each of the equations the scheme of the Section 3. We define radii and energies of circular orbits by means of the conditions: Then we single out asymptotics of these functions of λ by means of the relations: i + · · · , µ(λ) = 1 + λµ and, using the relations we reformulate the equation (26) in terms of the dimensionless variable ξ 1 and the spectral parameters ǫ 1 while the equation (27) -in terms of ξ 2 and ǫ 2 . Finally, we perform expansion of the equations into powers of λ and solve them separately. Now we are going to take actual coupling of the equations (26) and (27) into account. First of all, we note that the variables ξ 1 and ξ 2 are not of one another, and the spectral parameters ǫ 1 and ǫ 2 are also not independent. Thus we should choose common variables in both equations.
On the contrary, the function w 2 may have a different behaviour at λ → 0. Here we consider three cases.
1. Let r 2∞ = r 1∞ and b 2∞ = b 1∞ . Then w 2 = O(λ −n ), n ≥ 0 (except perhaps very special examples which we do not consider). In this case one can solve formally the equation (30) in favour of ψ 2 (ξ) as follows: This representation leads to the loss of solutions for ψ 2 which are not analytical in λ and thus have nothing to do with the perturbation procedure. The use of (34) in the r.h.s. of (29) permits us to eliminate ψ 2 from (29) and thus to obtain a close wave equation for ψ 1 . The structure and treatment of this equation are the same as those of the equation (15). Moreover, it is obvious from (34) that at least the zero-and the first-order terms of ψ 2 vanish. Thus the r.h.s. of (29) does not contribute in lower orders of perturbation procedure. In zero-order approximation we have the oscillator problem.  (30) is still valid. The only difference from the case 1 is that the r.h.s. of equation (29) may contribute in the first order of λ.
In both the above cases we used the dimensionless variable ξ 1 and obtained a closed eigenstate equation (which we will reference to as the problem 1) for the wave function ψ 1 (ξ 1 ) and the spectral parameter ǫ 1 . We can proceed with the variable ξ 2 and obtain the problem 2 for the function ψ 2 (ξ 2 ) and the parameter ǫ 2 . One might be inclined to think that both problems 1 and 2 are equivalent and lead to the same spectrum (in terms of energy E). Actually, different problems complement one another. This is evident from equation (28) leading to the relation: Indeed, in both 1 and 2 cases |ǫ 2 − ǫ 1 | → ∞ if λ → 0. It does mean that an arbitrary energy level E calculated by means of eigenvalue ǫ   2 ) are small and do not change qualitatively this picture. Thus different problems generate different branches of the energy spectrum of the original set of equation. In this respect the following special case differs essentially from the previous ones.
In the zero-order approximation we obtain the coupled pair of wave equations (on the contrary to the cases 1 and 2 where we had a single wave equation). In physically meaningful cases (of Section 5, for example) they have the form: where χ = lim λ→0 y = const, and parameters ν i , κ and ω are related to the functions w i (i = 1, 2) by the equation of the type of (23), (24) and (25). The equations (35), (36) can be evidently reduced to the pair of similar equations but with parametersν i = {ν 1 +ν 2 ± (ν 1 − ν 2 ) 2 + 4χ 2 }/2 (i = 1, 2) andχ = 0. Thus they become split equations of the form (22). The eigenvalues ǫ corresponding to the first and second equations are separated by finite constantν 1 −ν 2 . Thus the corresponding states mix in higher orders of perturbation procedure.

Application: Regge trajectories of mesons
Here we apply the pseudoperturbative treatment of 2BDE in meson spectroscopy. It is known [27] that spectra of heavy mesons are described well by the nonrelativistic potential model with QCD-motivated funnel potential u(r) = u l (r) + u C (r), where The Coulomb part (37) of this potential describes a nonrelativistic limit of the vector onegluon exchange interaction while the linear part (38) is suggested by the area law in the lattice approximation of QCD and has presumably scalar or scalar-vector nature. Description of light meson spectroscopy needs application of appropriate relativistic models. Most of them are related to the string theory. From the theoretical viewpoint the most interesting are QCD-motivated relativistic models embracing properties of both heavy and light mesons. Such models should reflect the scalar-vector structure of interaction and should lead to funneltype potential in the nonrelativistic limit.
A natural candidate for the relativistic potential model is the2BDE with a short-range vector potential and a long-range scalar one. At least three general structures of vector potential are used in the literature, with u v (r) = u C (r) or another short-range potential; here u ′ (r) = du(r)/dr. The potential (39) is only a static part of vector interaction (see [5]). The relativistic vector field kinematics is taken into account in the potential (40) (see [16,8]) which, for the Coulomb case, was first proposed by Eddington and Gaunt [28,29]. In the generalization (41) of the Breit potential [1] retardation terms have been added [8]. Two different scalar potentials, come from different couplings of scalar mediating field with fermionic fields. The first one (42) arises from the Yukawa interaction (see [6]) while the second one (43) corresponds to so called "minimal" coupling [15]. The latter and also two following potentials can be treated as static approximation of various QFT-motivated scalar quasipotentials [30,31,32,17]: The perturbative treatment of Breit-type equations has been used for calculating a fine splitting in spectra of heavy mesons [14,16]. Light meson spectra are essentially relativistic and need a nonperturbative statement of the problem which is inconsistent because of non-physical singularities of radial equations. To avoid these difficulties in numerical calculations one is forced to invent sophisticated potentials and impose rather artificial boundary conditions [15].
Using the pseudoperturbative treatment of 2BDE with different combinations of potentials (37)-(45) we obtain analytical expressions for meson mass spectra and estimate a role of general structure and input parameters of potentials in the model. We consider mass spectra of lightest mesons (containing u and d quarks only) and try to reproduce their following general features: i) Mass spectra of light mesons fall into the family of straight lines in the (E 2 , j)-plane known as Regge trajectories.
iv) Spectrum is ℓs-degenerated, i.e., masses are distinguished by ℓ (not by j) and n r .
v) States of different ℓ possess an accidental degeneracy which fact causes a tower structure of the spectrum.
For this purpose we use the nonrelativistic potential function (37) and (38) in vector and scalar potentials of different spin structure (39)-(45) and calculate pseudoperturbative spectrum in zero-order approximation. Classification of states then is done using singlet-triplet properties of large-large component of wave function (2) in the nonrelativistic limit.
If the vector short-range interaction is ignored and scalar potentials (42)-(45) are used with u s (r) = ar the pseudoperturbative mass (i.e., energy) of meson in zero-order approximation has the following form: here m + = m 1 +m 2 , and k, η, ζ, δ 1 , δ 2 , κ are dimensionless constants depending on the potential chosen. Four families of energy levels E i (i = A, 0, −, +) form trajectories in the (E 2 , ℓ)-plane which are nearly straight. Indeed, ζ = 0 ÷ 2 for all the potentials considered and rest masses m a (a = 1, 2) of lightest mesons are small compared to the energy scale √ σ. Thus the parameter ζm + √ a determining a curvature of trajectories is small. Parameters δ 1 and δ 2 determine a common shift of all the trajectories and are not important for the present discussion. Below we discuss the calculated values of parameters k, κ and η determining the slope of trajectories and their degeneracy properties.
In the (42) case k = 4 so that the slope σ = ka matches quite well to that of property ii); η = 2 causes accidental degeneracy typical for the harmonic oscillator; but κ = 4 leads to j-dependence of energy (not ℓ-dependence) so that the ls-degeneracy is absent.
In the (43) case k = 4 and η = 2, so that the slope and the accidental degeneracy are the same as in the (42) case; κ = 4 − 3 √ 2 ≈ −0.243 provides an approximate ls-degeneracy, with accuracy 6 %; the splitting is of order of the actual ss-splitting (see property vi)).
Taking into account the vector short-range interaction (one of potentials (39)-(41) with u v (r) = u C (r)) results in a parallel shift of Regge trajectories. The value of the shift is of the order αa, it depends on the vector potential chosen and is different (in general) for different trajectories E i (i = A, 0, −, +).
It has been proved in the framework of single-particle Dirac equation a possibility of confinement by means of vector and equally mixed vector-scalar long-range interactions [33,34,35]. We examined these cases in 2BDE approach using different vector potentials (39)-(41) with u v (r) = ar. Corresponding zero-order pseudo-perturbative spectra have a form similar to (46). The difference is that k = k i = 8 ÷ 12 is two times or more larger than desired, and is different for i = A, 0, −, + (i.e., trajectories are not parallel).

Summary
The Breit equation and its generalizations (2BDE) possess non-physical singularities. In some cases these points lay far from the physically important domain but they make a boundary problem incorrect or physically improper [17,16]. In order to avoid this difficulty and to use the 2BDE in the relativistic bound state problem, especially for the case of strong coupling, we develop the 1/j expansion method.
The method is based on the large-N or large-ℓ techniques applicable to the radial Schrödinger equation. In our case the 2BDE is reduced to the coupled pair of quasipotential-type equations which structure causes principal modification of known techniques. Other changes are related to the fact that the equations represent a nonlinear spectral problem with cumbersome quasipotentials.
We apply this pseudoperturbative method to the 2BDE with the linear+Coulomb potential of different scalar-vector structure. In all cases in the zero-order approximation it was obtained the Regge trajectories which are linear asymptotically. Linear potentials of two scalar structures (43) and (45) which was discarded in [17] as nonphysical (because of singularities in 2BDE) reproduce well in our case general properties of light meson spectra. In particular, the slope σ = ka of light meson trajectories fit well to the experimental value if the parameter a is taken from the nonrelativistic potential model [27]. The third linear potential of [17] (with no singularities in 2BDE) does not match to experimental data.