Heavy-Light Mesons In A Relativistic Model

We study the heavy-light mesons in a relativistic model, which is derived from the instantaneous Bethe-Salpeter equation by applying the Foldy-Wouthuysen transformation on the heavy quark. The kernel we choose is based on scalar confining and vector Coulomb potentials, the transverse interaction of the gluon exchange is also taken into account in this model. The spectra and wave functions of $D$, $D_s$, $B$, $B_s$ meson states are obtained. The spectra are calculated up to the order of $1/m_Q$, and wave functions are treated to leading order.


I Introduction
Recently, great progress has been made in the measurement of the spectrum of heavy-light quark-antiquark system, especially in the charm sector. Several new excited states were established in both D and D s mesons [1][2][3][4][5][6][7][8]. Many theoretical models were propounded to identify the new discovered states and explain their properties [9][10][11][12][13][14]. The spectroscopy provides a powerful test of the theoretical predictions based on the quark model in the standard model.
Heavy-light quark-antiquark system Qq plays an important role in understanding the strong interactions between quark and antiquark. Double-heavy meson sector is rather mature, and nonrelativistic potential models have been proven extremely successful for the description for heavy quarkonia [15]. As for heavy-light meson system, one needs a model that can include the relativistic effects of the light quark.
In order to describe the relativistic effects in the bound states, we resort to the Bethe-Salpeter equation. Many studies of mesons were carried out in the Bethe-Salpeter approach [16][17][18][19][20]. It is very difficult to solve the Bethe-Salpeter equation, especially when we consider states with large total angular momentum number J. We can reduce the Bethe-Salpeter equation to an approximate but easier form, so as to study mesons systematically. We apply the Foldy-Wouthuysen transformation on the heavy quark and expand the interaction terms to order 1/m Q . The instantaneous approximation can only contribute to corrections of order (1/m Q ) 2 if one considers the heavy quark on shell as in the heavy quark effective theory (HQET). In this paper we take the instantaneous approximation for the Bethe-Salpeter equation.
One can bring the instantaneous Bethe-Salpeter equation into the Breit form by taking further approximation. The Breit interaction is widely used to study quarkantiquark systems as an nonrelativistic or semirelativistic model [21][22][23]. While in this paper, we take the instantaneous Bethe-Salpeter equation in its original form, from which one can derive a more reliable model that preserves the relativistic effects of the light quark.
We treat the equation in the 1/m Q expansion to the leading order. The parameters in the equations are determined by fitting the spectra of the heavy-light meson states presented in Particle Data Group (PDG) [24]. The identifications of the newly observed heavy-light meson states are given based on the theoretical prediction.
The organization of this paper is as follows. Section II is for the relativistic model and the effective Hamiltonian of the heavy-light quark-antiquark system. In section III, we solve the relativistic wave equation. Section IV is for the numerical result and discussion. In section V we have a brief summary.

II The model
The Bethe-Salpeter Equation for the quark-antiquark system can be written as [25]: (p / 1 − m 1 )χ(p)(p / 2 + m 2 ) = d 4 p ′ (2π) 4 K(p, p ′ , P )χ(p ′ ), (1) where the total momentum P (in the center-of-mass system, P = (E, 0)) and the relative momentum p are defined as: with It is necessary to assume some form for the kernel in order to reduce the Bethe-Salpeter equation. Here we take the kernel as the simplest form [26], which can be reduced to that used in Ref. [27]. The kernel can be written as: where k is the transferred four-momentum, Inserting Eq. (5) into Eq. (1), the Bethe-Salpeter equation can be written as: The above equation can be transformed to a more convenient form that is symmetric with respect to particles and antiparticles by using the charge conjugation transformation. We rewrite Eq. (7) as: After applying the condition we arrive at where ψ(p) is the charge conjugated wave function, and the minus sign appeared before V s indicates that the confinement force is always attractive. By taking the instantaneous approximation, we can perform the p 0 integration and decompose the instantaneous Bethe-Salpeter equation into four equations. Details of this derivation can be found in Refs. [28,29].
The projection operators are defined as with Applying the projection operators, the four coupled equation can be written as where φ(p) = dp 0 ψ(p 0 , p), One can combine the four coupled equations and obtain the equation of φ instead of its projective components φ ±± as in Ref. [30], that is where Λ(p) = Λ If one makes the approximation Eq. (23) can be bought into the Breit form. The Breit equation can be used as an nonrelativistic or semirelativistic approach to study the interactions of quarks. In order to take the relativistic corrections of the light quark into account, we keep Λ(p) in its original form instead of taking it as 1.
Eq. (23) can be transformed into coordinate space according to The coordinate form of Eq. (23) is with where V v (r) and V s (r) are the Fourier transformations of V v (−k 2 ) and V s (−k 2 ), respectively. The wave function φ(r) has two spinor indices, i.e. 4 × 4 = 16 components, it is difficult to solve an eigenequation with so many components. For the heavylight system, one can reduce the eigenequation by applying the nonrelativistic approximation of the heavy quark Q and expand the interactions in the order of 1/m Q .
In an elegant and systematic way the reduction can be achieved by using the Foldy-Wouthuysen transformation.
Here we designate the heavy quark Q and the light antiquarkq by the index 1 and 2 in Eq. (27), respectively. Thus the transformation is performed on the heavy quark with the superscript "1". We reduce the Hamiltonian in the Dirac representation, where If the original Hamiltonian is written in the form where O is the "odd" operator, typical examples are the matrices α and γ, while E is the "even" operator, examples of this class of operators are 1, β and Σ. According to the Foldy-Wouthuysen transformation, the transformed Hamiltonian reads The Hamiltonian in Eq. (27) is  [31,32]. We find that our goal can be achieved by employing the constraint of Eqs. (18) and (19) which are dropped in Eq.
We write Eq. (27) together with Eq. (35) as: It is easy to verify that the above two equations are equivalent to the instantaneous Bethe-Salpeter equation, i.e. the four coupled equations Eqs. (16)∼ (19). Here we expect to obtain one single equation that is equivalent to the above two equations by combining them.
Eq. (37) is equivalent to If the above equation is true, then on the other hand, if Eq. (37) is true, we take ϕ = h 1 φ, then Applying the equivalent form of Eq.(37), Eqs. (36) can be transformed to a new equivalent equation it is easy to verify that the above equation is equivalent to the two coupled equations (36) and (37). Furthermore, it is equivalent to Noticing the relations we can have two different equivalent forms of Eq. (42), they are and Here we guess a form that is equivalent to Eq. (46) and Eq. (47) by employing the common part of them, that is Left multiplying the above equation by 1 2 (1 + h 1 h 2 ), we get the form of Eq. (46). On the other hand, we can obtain the form of the above equation from Eq. (47) if we take ψ = 1 2 (1 + h 1 h 2 )ϕ. That is to say, Eq. (48) is equivalent to Eq. (46) and Eq. (47), thus equivalent to Eqs. (36) and (37), i.e. the instantaneous Bethe-Salpeter equation.
In the foregoing paragraphs, equivalent equations of the instantaneous Bethe-Salpeter equation are obtained. In Appendix A, the reduction of the instantaneous Bethe-Salpeter equation for nonrelativistic systems is discussed. We find that the reductions of the Hamiltonian with the approximation (25), i.e. the Breit interaction, and without the approximation (25), i.e. the choice in this paper, do not differ until order 1/m i m j . While for the case of heavy-light system we focus in this paper, the reduction results of the two scheme can be essentially different. Now we return to the reduction of the Hamiltonian for the heavy-light system. Eq. (49) can be written as (the wave function in the equation is omitted): where now we consider performing the Foldy-Wouthuysen transformation on Eq. (51), the left of it can be rewritten as: the odd and even operators in the above equation are: where h.c. stands for Hermitian conjugate. By using the new form of Eq. (33) wherẽ The last termH 2 is of the order V 2 v . Comparing to the other terms, its correction is suppressed by the coupling constant, here we omit it. Considering the substitution β (1) → 1 in the Foldy-Wouthuysen transformation,H 1 can be reduced tõ thenH is transformed tõ and we obtain the final Hamiltonian The scalar and vector potentials we choose have the simple form The potentials are chosen to have a Coulombic behavior at short distance and a linear confining behavior at long distance.

III Solution of the Wave Equation
In this section, we solve the eigenequation of H 0 and then deal with the perturbative corrections of H ′ .
In order to generalize the algorithm which was used in Refs. [33,34], we prove a relation in the first place.
Let Ω(p) be the pseudo-differential operator function and Ω(k) be the normal function, here p and k stand for the modules of momentum operator p and momentum k, respectively. Let ψ(r) be a function of r that can be written as with the help of we have The exponential factor e ik·r can be decomposed into series of spherical harmonics where j l is the l-th order spherical Bessel function, Y lm (r) is the spherical harmonics, which satisfies the normalization condition withr the unit vector along the direction of r. Inserting Eqs. (68) and (71) into Eq. (70), after simplification we arrive at If we take Ω(p) = 1, then Ω(k) = 1, the above equation is simplified as: The above equation is true for arbitrary function φ(r), which gives the orthogonality condition for the spherical Bessel function Furthermore, if we take φ(r) = j l (k ′ r), we can transform Eq. (73) to a simple form by using Eq. (75), that is Ω(p)j l (kr)Y lm (r) = Ω(k)j l (kr)Y lm (r).
The above equation is actually an eigenequation of Ω(p). If we expand the unsolved wave function in terms of the spherical Bessel function, the above relation can do a great help in solving the wave equation. Now we turn to solving the eigenequation of H 0 . The spinor part of Eq. (107) can be rewritten as by using the formula we have are used. Inserting Eq. (12) and (15) into Eq. (64), we can rewrite H 0 in the matrix form then H 0 can be written as: where Since the quark and antiquark are bound in the meson, when the distance between them is large enough, the wave function drops dramatically. The wave function will effectively vanish at a typically large distance L, the quark and antiquark can be viewed as being restricted in a limited space, 0 < r < L. The unsolved functions f (r) and g(r) can be expanded in terms of the spherical Bessel function in the limited space, we write the basis as where the superscripts "A" and "B" stand for the upper and lower parts of Eq. (78), and N n is the module of the spherical Bessel function with a n the n-th root of the spherical Bessel function j l (x) = 0. N n is added in accordance with the normalization condition in Eq. (111). In the limited space, the momentum k is discrete, we can have the relevance a n L ⇐⇒ k.
The eigenfunction of H 0 can be expanded in the orthonormalized basis {ψ A i , ψ B α }: In the numerical calculation the above summation can be truncated at a large integer N , According to Eq. (87), we can rewrite the eigenequation of H 0 in the representation of {ψ A i , ψ B α }. In this representation, the operator H 0 has its matrix form Applying Eqs. (77), (80) and (81) with the normalization condition Eq. (110), we can get the matrix elements of H 0 easily.
here we define a symbolic notation Diagonalizing the Hermitian matrix of H 0 , we can get the eigenenergy of H 0 and the coefficients g i , f α , which are defined in Eq. (97), then the eigenequation associated with H 0 is solved and the eigenfunction is obtained. Now we discuss the perturbative corrections of H ′ . It is easy to verify that the operators {j 2 , j z , K, S z } are a set of mutually commuting operators which commute with H 0 , where j = L + S (2) , S (2) = 1 2 Σ (2) , K = β (2) (Σ (2) · L + 1), S = 1 2 σ (1) . Then the eigenstates of H 0 can be labeled by the corresponding set of quantum numbers {n, j, m j , k, s} and the eigenequation associated with H 0 can be written as where Ψ (0) n,k,j,mj ,s (r) = g n,l,j (r)y mj j,l (θ, ϕ) if n,l,j (r)y mj j,2j−l (θ, ϕ) In Eq. (108), Y m±1/2 l is the spherical harmonics. For j = l ± 1/2, we have dΩ (y mj j,l ) † y mj j,l = 1. (110) Thus the normalization condition of Eq. (107) is read as ∞ 0 dr r 2 (f 2 n,l,j + g 2 n,l,j ) = 1.
For a state with quantum number j, the operator K can have two opposite eigenvalues ±(j + 1/2) with where the quantum number k is defined as n,k,j,mj,s (r).
The zeroth order invariant mass E 0 n,k,j is equivalently determined by n, j, l, and the eigenstate has the parity P = (−1) l+1 . The Solution of the eigenequation associated with H 0 is detailed in the next section.
The perturbative term H ′ does not commute with any of the operators introduced in solving the eigenequation of H 0 , but still commutes with Mixing can happen between states with same quantum numbers J, M J , P , but with different quantum number j. In this paper the mixing between states with j = l±1/2 is considered. Mixing can be described by the mass matrix. The mass matrix is calculated perturbatively in the basis, which can be written as where i,j 's denote the basis states which are defined in Eq. (114).
In Eq.(65), there is a term (1 + h 2 )/2 on both left and right sides of the expression. We can let (1 + h 2 )/2 operate on the wave functions of H 0 and obtain new wave functions, then deal with the middle term of Eq.(65) with the new wave functions.
the above equation shows the transformation from the wave functions of H 0 to the new wave functions with the operation of (1 + h 2 )/2. With the new wave functions the perturbative contributions of the H ′ to the eigenvalues of states can be calculated easily, they are given below: (1) J P = 0 − Only state with j = 1/2, l = 0 can construct the J P = 0 − state, we have (2) J P = 1 − Both states with j = 1/2, l = 0 and j = 3/2, l = 2 can construct the J P = 1 − state, we have and (3) J P = 0 + Only state with j = 1/2, l = 1 can construct the J P = 0 + state, we have E n,1, 1 2 ,0 = E n,1, 1 (4) J P = 1 + Both states with j = 1/2, l = 1 and j = 3/2, l = 1 can construct the J P = 1 + state, the matrix elements of the mass matrix are With the matrix elements given above, one can get the eigenvalues of the two mixing states and the mixing angle easily by diagonalizing the mass matrix.
Once the eigenfuncion of H 0 is obtained, the perturbative correction of H ′ can be calculated directly by following the procedures presented above.

IV Numerical result and discussion
The parameters used in this work are the quark masses, m u,d , m s , m c , m b and three potential parameters α s , b, c. In the calculation of the spectra, we find that the model is capable to give reasonable spectral structure for each meson system. In the framework of Bethe-Salpeter equation, the parameter b is responsible for the energy levels of states with higher quantum numbers n or l. We find that if the parameter b is taken as a constant for all the meson systems, the energy levels for the excited meson states decrease as the value of Numerical calculation shows that the solution of the wave equation is stable when L > 5 fm, N > 35. Here we take L = 10 fm, N = 50. We fit the spectra of the heavy-light D, D s , B, B s mesons mainly based on the meson states presented in PDG [24]. The numerical results for the spectra of heavy-light mesons are presented in two tables. Table I is for D, D s mesons and Table II for B, B s mesons. The obtained spectra are in reasonable agreement with the experimental measurements. Theoretical deviations from experimental data mainly appear in the D s meson sector, where the calculated mass of the D * s0 (2317) 0 resonance is the worst. In the two tables, one can see that apart from the D * s0 (2317) 0 and D s1 (2460) resonances, mass calculations for the rest of the resonances are in good agreement with experimental data. The discrepancy may be ascribed to the naive assumption of the kernel, the dropped α 2 s term in the Hamiltonian can also make the result incomplete. The highly excited meson states are also calculated in the spectra and the newly observed charmed meson states are identified in our model.
As for D mesons, several resonances are observed by [4]. D sJ (2632) is measured to be 2632.5±1.7 MeV in Ref.
[1], we assign it as the |2 1/2 S 0 state with J P = 0 − . For D * s1 (2710), it is measured to be 2708 ± 9 +11 −10 Mev in Ref. [3], it can be assigned as the |2 1/2 S 1 state with J P = 1 − , which agrees with Refs. [35,36]. The D sJ (3040) resonance is observed by BABAR Collaboration at a mass of 3044 ± 8 +30 −5 Mev [4]. Several possible assignments of this resonance is discussed in Ref. [9]. It is identified as the |2 1/2 P 1 state with J P = 1 + in our predicted mass spectrum for D s meson. Recently, LHCb collaboration identifies the resonance D * sJ (2860) as an admixture of spin-1 and spin-3 resonances and determines the masses and widths of the two states to be [5,6] In Table I, One can see that our predictions for the masses of the |1 3/2 D 1 and |1 5/2 D 3 states are around 2860 Mev. Our model appears to be able to interpret the two states being the J P = 1 − and 3 − members of the 1D family. We can assign D * s1 (2860) − as the |1 3/2 D 1 state with J P = 1 − and D * s3 (2860) − as the |1 5/2 D 3 state with J P = 3 − .
The wave function of each bound state can be obtained simultaneously when solving the wave equation. The radial wave functions g n,l,j (r) and f n,l,j (r) for D meson are depicted in Fig.1 as an example. [g 2 n,l,j (r) + f 2 n,l,j (r)]r 2 is the possibility density distributed along the quarkantiquark distance r. In Ref. [37], the radial wave functions are also presented, they are calculated in a model derived by reducing the spectator equation in relation to the Bethe-Salpeter equation [26]. One can find that their wave functions share the same pattern with ours, though the Hamiltonian H 0 in our model looks very different from the Dirac-like Hamiltonian used by them.

V Summary
We construct a relativistic model for heavy-light quark-antiquark systems by studying the reduction of the instantaneous Bethe-Salpeter equation. The kernel we All the terms of the transformed Hamiltonian are the same as that of the Breit interaction except for the last term. In our scheme, we have OO = 1 2 (β (1) + β (2) )V 2 1 2 (β (1) + β (2) ), (148) in the heavy quark limit, E 0 ≈ m 1 + m 2 . While with the approximation, the result is Eq.(149) differs from Eq.(151) by a minus sign. In conclusion, for double-heavy system we can take the approximation and have a simpler form of wave equation if we disregard the last term in Eq. (147).