Spectra of heavy–light mesons in a relativistic model

The spectra and wave functions of heavy–light mesons are calculated within a relativistic quark model which is based on a heavy-quark expansion of the instantaneous Bethe–Salpeter equation by applying the Foldy–Wouthuysen transformation. The kernel we choose is the standard combination of linear scalar and Coulombic vector. The effective Hamiltonian for heavy–light quark–antiquark system is calculated up to order 1/mQ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1/m_Q^2$$\end{document}. Our results are in good agreement with available experimental data except for the anomalous Ds0∗(2317)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{s0}^*(2317)$$\end{document} and Ds1(2460)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{s1}(2460)$$\end{document} states. The newly observed heavy–light meson states can be accommodated successfully in the relativistic quark model with their assignments presented. The DsJ∗(2860)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{sJ}^*(2860)$$\end{document} can be interpreted as the |13/2D1⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|1^{3/2}D_1\rangle $$\end{document} and |15/2D3⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|1^{5/2}D_3\rangle $$\end{document} states being members of the 1D family with JP=1-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J^P=1^-$$\end{document} and 3-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3^-$$\end{document}.


Introduction
Great experimental progress has been achieved in studying the spectroscopy of heavy-light mesons in the last decades [1][2][3][4][5][6][7][8]. In the charm sector several new excited charmed meson states were discovered in addition to the low-lying states. For D J mesons, the excited resonances D(2740) 0 , D * (2760) [1], D J (2580) 0 , D * (2650) and D * (3000) [2] were found in the D ( * ) π invariant mass spectrum by the BaBar and LHCb Collaborations. For D s J mesons, besides the wellestablished 1S and 1P charmed-strange states, the excited resonances D s J (2632) [3], D s J (2860) [4], D s J (2700) [5] and D s J (3040) [6] were observed in the D ( * ) K invariant mass distribution by the two collaborations. In the b-flavored meson sector, several excited states were studied in experiment as well as the ground B and B s meson states [9]. The strangeless resonances B J (5840) 0 and B(5970) 0 were found in the Bπ invariant mass spectrum by the LHCb a e-mail: liujb@ihep.ac.cn b e-mail: lucd@ihep.ac.cn and CDF Collaborations, respectively [10,11]. The stranged B * s J (5850) were observed in the B ( * ) K invariant mass distribution by the OPAL Collaboration [12].
The heavy-light meson spectroscopy plays an important role in understanding the strong interactions between quark and antiquark. Meanwhile, it provides a powerful test of the various phenomenological quark models inspired by QCD. Heavy-light mesons have been investigated extensively in relativistic quark models [13][14][15][16][17][18][19], where many relativistic potential models are constructed by modifying or relativizing nonrelativistic quark potential models and additional phenomenological parameters are employed. For the heavylight system one needs a model that can retain the relativistic effects of the light quark. In this work we resort to the originally relativistic Bethe-Salpeter equation [20]. The Bethe-Salpeter approach was widely used in studying mesons so as to embody the relativistic dynamics [21][22][23][24][25][26]. It is rather difficult to solve the Bethe-Salpeter equation for meson states, especially when considering states with large angular momentum quantum number. In order to study the spectrum of heavy-light mesons systematically, we choose to reduce the Bethe-Salpeter equation in the first place.
In our previous work [27], we apply the instantaneous approximation and obtain an equation equivalent to the Bethe-Salpeter equation. The Hamiltonian for the heavylight quark-antiquark system is expanded to order 1/m Q by applying the Foldy-Wouthuysen transformation to the equivalent equation. We find that the leading Hamiltonian is actually not Dirac-like. The interaction we derive is essentially different from the Breit interaction [28][29][30]. In this paper we extend and improve our study of the spectrum of the heavylight mesons D, D s , B and B s . The running of the coupling constant is considered. Moreover, the 1/m 2 Q correction is calculated. Many papers have only considered the leading 1/m Q term in the heavy-quark expansion [27,[31][32][33][34]. Our calculation shows that the 1/m 2 Q corrections to the masses of the mesons are around 50 MeV, which is too large to be neglected. The parameters in the equations are determined by fitting the masses of the 1S and 1P meson states presented by particle data group (PDG) [9], while the states beyond 1P are calculated as a prediction. We find that in the Bethe-Salpeter formalism the linear confining parameter, i.e. the string tension, actually depends on the masses of the constituent quark and antiquark in mesons. The large discrepancy between experimental data and our previous work is decreased in this work. The newly observed heavy-light meson states can be accommodated successfully in our predicted spectra.
This paper is organized as follows. In the next section, we have a brief review of the relativistic quark model. Section 3 is for the solution of the wave equation and the perturbative corrections. In Sect. 4 we have numerical results and discussions. The last section is for a brief summary.

The model
According to the conventional constituent-quark model, the mesons can be seen as a composition of a quark and an antiquark. In the Bethe-Salpeter formalism, the eigenequation for quark-antiquark systems has the general form [20] (1) where p 1 and p 2 relate to the total momentum P and the relative momentum p, as follows: Using the energy-momentum conservation, i.e. p 1 + p 2 = p 1 + p 2 , Eq. (1) can be simplified as Here we choose the interaction kernel as the standard Coulomb-plus-linear form, which is one-gluon-exchange (OGE) dominant at short distances with linear confinement at long distances. If one applies the instantaneous approximation, i.e. neglecting the frequency dependence, the kernel can be written as where the transferred momentum k is defined as Since the interaction kernel K ( p, p , P) is no longer dependent on p 0 , we can perform the integration over p 0 in Eq. (4). After transforming the instantaneous Bethe-Salpeter equation into coordinate space, the wave function of the eigenequation decouples from the time coordinate [35][36][37].
In our previous work [27], with the help of projection operators for the wave function we found that the instantaneous Bethe-Salpeter equation is equivalent to the following equation: where the superscript "1" and "2" stand for the heavy quark Q and the light antiquarkq in the Qq meson, respectively. The operators in the above equation are defined as with the free Dirac Hamiltonians Inserting Eqs. (10) and (11) into (9), we can verify the relation The interaction potential U (r) in Eq. (7) is directly derived from the instantaneous Bethe-Salpeter equation and closely related to the interaction form we assumed in the kernel. It can be written as with where V (r ) and V (−k 2 ) are related to each other according to Fourier transformation. The instantaneous Bethe-Salpeter equation as an integral equation is equivalent to a less complicated differential equation shown in Eq. (7) but it is still difficult to solve. For heavylight systems, the heavy-quark effective theory is applied. It is reasonable to consider the heavy-quark expansion, i.e. the 1/m Q expansion. One can reduce the equivalent eigenequation by calculating the interactions of the heavy-light quarkantiquark meson order by order.
Our goal can be achieved by employing the Foldy-Wouthuysen transformation [38]. The operators involved in Eq. (7) can be divided into two sets: the "odd" O and the "even" E. The name "odd" denotes the operators couple the large and small components of the Dirac spinor, while the "even" operators are diagonal with respect to the large and small components. The main idea of the Foldy-Wouthuysen transformation is to apply a unitary transformation U which retains the "even" operators and eliminates the "odd" operators. If one writes the original Hamiltonian as according to Foldy and Wouthuysen, one obtains the transformed Hamiltonian: The reduction by performing the Foldy-Wouthuysen transformation on Eq. (7) has been detailed in our previous work [27]. Instead of βm being the main term in the common Dirac Hamiltonian shown in Eq. (16), the dominant term is β E in our case: The reduction result is calculated to order 1/m Q in our previous work [27]. With the similar procedure, here we extend the result to order 1/m 2 Q . By inserting the "odd" and "even" operators of Eq. (7) into Eq. (17), we obtain the Hamiltonian expansion. After the Foldy-Wouthuysen transformation, we havẽ with The perturbative termH consists of various terms of order 1/m Q and 1/m 2 Q . We divide it into three parts: wherẽ We can simplify the above equations by inserting an identity matrix (γ (1) 5 ) 2 = 1 between two odd operators of the heavy quark, with the help of the relations {γ 5 , β} = 0, [γ 5 , α] = 0, γ 5 α = and the substitutions β (1) → 1, (1) → σ (1) . Moreover, we can take the substitution h 2 → 1 if h 2 appears at the ends of the expression ofH as in Eqs. (22)(23)(24), since the corrections ofH are calculated as a perturbation toH 0 .
With the considerations above, we obtain our final Hamiltonian where the leading order Hamiltonian H 0 has the form and the subleading Hamiltonian H to order 1/m 2 Q can be written as with the interaction potentials U 1 (r), U 1 (r) and U 2 (r) in the above equations are defined as The leading order Hamiltonian H 0 we obtain for the heavy-light quark-antiquark system in Eq. (26) is not Diraclike as in Refs. [33,39]. We have Its form is more like the form used in relativized quark models [13,40,41]. As for the double-heavy system, we have h 2 → 1 and β (2) → 1, then Eq. (26) can be reduced to which is the Schrödinger formalism extensively used in nonrelativistic or semirelativistic quark models.

Solution of the wave equation
In this section we solve the eigenequation of the leading order Hamiltonian H 0 in Eq. (26). Before doing this we would like to discuss the properties of the solution of the eigenequation associated with H 0 . The eigenequation of H 0 can be written as the above equation is equivalent to which is equivalent to From Eqs. (36) and (38), we have and since (h 2 ) 2 = 1, When we take c = −1, Eq. (36) is transformed to which is not the correct eigenequation for the bound systems we are interested in. Thus we only have c = +1. This is the reason for the substitution h 2 → 1 we use in the last section. If all the eigenfunctions of H 0 for bound states satisfy the relation the eigenfunction set of H 0 is not complete. A complete set is needed to construct the identity operator 1 = i |ψ i >< ψ i | in order to calculate the perturbative correction of H b . Thus we construct a new Hamiltonian. Inspired by the relation h 2 ψ = ψ we transform the potential term in Eq. (36) as then the new Hamiltonian we construct can be written as It is easy to verify that the eigenfunction set of the new Hamiltonian includes both subsets: where the subset {ψ + } is identical to the eigenfunction set associated with the original Hamiltonian H 0 in Eq. (36). Now we turn to solving the eigenequation associated with the new Hamiltonian H 0 , that is, In the heavy-light quark-antiquark system, we treat the heavy quark as a static source, while the light one is described relativistically by a Dirac spinor. It is easy to verify that H 0 commutes with all the elements of the standard operator set { j 2 , j z , K , S z } associated with the free Dirac Hamiltonian. Then the eigenstates of H 0 can be labeled by the quantum number set {n, j, m j , k, s} corresponding to the operator set.
The quantum number k can have two opposite values for an eigenstate with quantum number j: The leading order invariant mass E (0) can be determined by quantum numbers n, j and k, or equivalently by n, j and l. The parity of the bound states is determined by P = (−1) l+1 . The Dirac spinor with quantum numbers j m j and l can be written as where the subscripts l A and l B stand for l and 2 j − l, respectively. The complete expression of y m j j,l (θ, ϕ) can be found in Ref. [34].
For a bound state of a quark and an antiquark the wave function will effectively vanish when the distance between them is large enough. We designate such a large typical distance as L. Then the heavy quark and light antiquark bounded in the meson can be viewed as restricted in a limited space, 0 < r < L. Thus we can expand the radial functions f (r ) and g(r ) by spherical Bessel functions associated with the distance L: where N n and a n are the module and the nth root of the spherical Bessel function j l (r ), respectively. Inserting Eqs. (9) and (11) into (43), we can rewrite H 0 in the matrix form where h.c. stands for Hermitian conjugate and the operator elements are According to Eqs. (48) and (49) we can rewrite the eigenequation of H 0 in the representation of the state basis constructed from spherical Bessel functions. In this representation the operator H 0 can be written in its matrix form: The matrix elements of H 0 in the above equation can be calculated by applying the relation and the eigenequation [27] ( where ( p) is a pseudo-differential operator function and (k) is a normal function, p and k stand for the modules of momentum operator p and momentum k, respectively.
With the normalization condition we easily obtain the matrix elements of H 0 . For the operators regarding the energy of the motion, we have In order to write down the expressions of the elements associated with the interaction potential in a compact form we introduce a symbolic notation: Then we have After calculating every element of the Hermitian matrix of H 0 we diagonalize the Hermitian matrix and obtain eigenvalues and eigenvectors. The eigenvalue of the matrix is the eigenenergy of H 0 , while the eigenvector is associated with the coefficients g i , f α , which are defined in Eqs. (48) and (49). That is to say, the eigenequation shown in Eq. (45) is solved and the corresponding eigenenergy and eigenfunction are obtained.
Here we turn to discussing the perturbative corrections of H defined in Eq. (27). The perturbative term H does not commute with the standard operators introduced for the free Dirac Hamiltonian, but it still commutes with the total angular momentum operator J = j + S and the parity operator P of the bound state.
Thus the quantum number set associated with the total Hamiltonian H = H 0 + H can be denoted by {n, J, M J , P}. By using Clebsch-Gordan coefficients the total wave function of the heavy-light quark-antiquark bound state can be decomposed as follows: with which the corrections and mixings caused by H can be calculated perturbatively. The 1/m Q and 1/m 2 Q perturbative terms are given in Eqs. (28)-(30).
The properties of the eigenfunctions of H 0 are of great help in calculating the perturbative corrections. We have already used h 2 → 1 to get rid of the h 2 at the ends of the perturbative terms. As for the h 2 sandwiched in H b , h 2 → ±1 can be applied due to Eq. (44). H b can be rewritten as Here we define two operators: Then we have As discussed at the beginning of this section the eigenfunction set of H 0 can be divided into two parts {ψ + , ψ − }, where ψ + and ψ − represent the physical and unphysical states, respectively. Inserting the identity operator consisting of the complete set of H 0 in Eq. (68) we obtain The correction in first order perturbation can be written as where δ E (2) n,l, j, From Eq. (69), the correction of H b for ψ + n can be written as With the eigenfunctions we obtain the 1/m Q and 1/m 2 Q corrections can be calculated. Then the masses of all the different J P states are determined.

Numerical results and discussions
The vector and scalar potentials are chosen to have a Coulombic behavior at short distance and a linear confining behavior at long distance. They can be written in a simple form: The running coupling constant α s (r ) in the vector potential is derived from the coupling constant α s (Q 2 ) in momentum space via Fourier transformation. It can be parametrized in a more convenient form [13]: where α i and γ i are parameters which can be fitted according to the behavior of the running coupling constant at short distance predicted by QCD. The behavior of α s (r ) is depicted in Fig. 1. In this work we use the same α s (r ) given in Ref. [13], where the α i and γ i parameters have the values α 1 = 0.25, α 2 = 0.15, α 3 = 0.20, and γ 1 = 1/2, γ 2 = √ 10/2, There are two free parameters in the scalar potential. One is the string tensor constant b which characterizes the confinement of the quark-antiquark system. The other is a phenomenological constant c which is adjusted to give the correct ground state energy level of the heavy-light meson state. The behavior of the parameters b in our model is quite different from those in usual quark models. The slope parameter b is the most essential parameter in all the different kinds of phenomenological quark models where the linear confinement is assumed. The parameter b determines the structure of the calculated spectrum, more specifically, it determines the Regge trajectories and the energy gaps between radial excitations. But unlike the Dirac Hamiltonian in Eq. (34) or the Schrödinger Hamiltonian in Eq. (35), the Hamiltonian for the heavy-light meson states in the Bethe-Salpeter formalism has a different form for the interaction potential, that is, where h 2 has the form In the above equation the diagonal part, i.e. the "even" operator, is the chief contributor to the eigenvalue of the eigenequation. If m 2 tends to infinity the Hamiltonian degenerates into the nonrelativistic case. But if m 2 tends to 0 an additional factor 1/4 will appear and it weakens the ability of the confinement parameter b in V s (r ) to elevate the energy levels of the excitations. That is to say, in the Bethe-Salpeter formalism the energy level is also sensitive to the light-quark mass m q . The experimental data shows that the mass splitting is similar in the D and D s mesons. For example, the energy gaps between 1 3/2 P 2 and 1 1/2 S 0 states are Thus in the Bethe-Salpeter formalism different slope parameters are required to coordinate with different constituentquark masses in order to recover the structure of the heavylight meson spectra. This is also true for the radial excitations. In Fig. 2, the energy gap E between the first radial excitation and the ground state is depicted as a function of m q . We take the D meson as an example to illustrate the dependence on the quark mass. The values of the parameters, which are fitted for the D meson spectrum, are fixed except for the light-quark mass of the D meson. From the different shapes of the dashed, dotted and solid lines according to the three schemes, i.e. the Schrödinger, Dirac and Bethe-Salpeter formalisms, one can find: • When m q is taken large enough the three schemes tend to give the same value for the energy gap E. It indicates the equivalence of the three schemes when dealing with double-heavy mesons. • In the region m q < 1 GeV, which is the case for heavylight mesons, the three schemes give quite different values for the energy gap. It has the pattern E Schr > E Dirac > E B−S . In order to give the same energy gap for a specific meson the confinement parameter should be chosen as b Schr < b Dirac < b B−S . The literature supports this sequence. For instance, b Schr is taken as 0.175 GeV 2 [41], 0.180 GeV 2 [13], b Dirac is taken as 0.257 GeV 2 [34], 0.309 GeV 2 [28], while b B−S can be taken up to 0.400 GeV 2 in this work.  55) is 300 × 300. Since a lot of integrals are involved the Gauss-Legendre quadrature is widely used in our numerical calculation. The spectra of the heavy-light D, D s , B, B s mesons are fitted based on the data given by PDG [9]. In this work we take the masses of the well-established 1S and 1P heavy-light meson states as our input for fitting the parameters. After the fitting the highly excited states beyond 1P are also calculated in the spectra and we identify the newly observed highly excited meson states in our model. The numerical results for the spectra of D, D s mesons are presented in Table 1, while the B and B s mesons are presented in Table 2. The calculated spectra are in good agreement with the experimental data. Our results are compared with the results of two other relativistic models [34,39]. One is derived by a quasipotential approach and the other is obtained by reducing the Bethe-Salpeter vertex function. The result in this work is improved compared with our previous work [27]. Taking the mass difference between the pseudoscalar state and the vector state for example, as shown in Table 1 the discrepancy from experimental data is decreased for D, D * , D s and D * s states as well as other states. Theoretical deviations from experimental data mainly occur in the D s meson sector, specifically, the D * s0 (2317) and D s1 (2460) resonances. Our calculations for the two resonances are about 100 GeV higher than their masses measured in the experiment. The discrepancy may be ascribed to the instantaneous approximation, the naive assumption of the kernel or the α 2 s (r ) contributions, i.e. the loop corrections. However, it is more likely to find an explanation beyond the naive quark model [42]. The masses of the two resonances predicted by the constituent-quark model are generally 100-200 MeV higher than experiments [33,34,[43][44][45]. The mass of D * 0 , 2318±29 MeV is almost identical to the mass of D * s0 , 2317.8±0.6 MeV. It cannot be explained in the conventional quark model if the difference between the two anomalous resonances in the model is merely their light-quark masses m s and m u,d . In this work the confinement parameter b takes different values for different systems but it still is not capable to explain the small mass difference of the two resonances.
As D * s0 and D s1 lie just below the DK and D * K threshold, respectively, the authors in Ref. [46] have suggested that the two resonances may be D * s0 (DK ) and D s1 (D * K ) molecular states, while in Refs. [47][48][49], D * s0 and D s1 are considered as cs states which are significantly affected by mixing with the DK and D * K continua. In Ref. [50], the authors suggest that the discrepancy of the calculated masses in quark models can be qualitatively understood as a consequence of self-energy effects due to strong coupled channels. In Refs. [51][52][53] the interpretation of the heavy J P (0 + , 1 + ) spin multiplet as the parity partner of the ground-state (0 − , 1 − ) multiplet is proposed. Both theoretical and experimental efforts are required in order to fully understand the nature of the anomalous D * s0 (2317) and D s1 (2460) states. As for the D mesons, in the mass region 2500-3000 MeV several resonances are measured by the LHCb Collaboration [2]. The assignments of these states are listed in the upper part of Table 1 [2,9] This work Previous work [27] [2,9] This work Previous work [27] Ref. [39] R e f. [  From Table 1 one can see our results favor the measurements.
As for the D s mesons, several states beyond the 1P state have been observed. Their masses and identifications are presented in the lower part of Table 1. Recently, the LHCb Collaboration identified D * s J (2860) as an admixture of two resonances: D * s3 (2860) − and D * s1 (2860) − [7,8], with their masses measured as 2859 ± 12 ± 6 ± 23 MeV and 2860.5 ± 2.6 ± 2.5 ± 6.0 MeV, respectively. In Refs. [34,39] cited in Table 1, their predictions do not favor this identification, with their calculations generally 60 MeV higher than the measured masses. While our results for both |1 3/2 D 1 and |1 5/2 D 3 are around 2860 MeV, the two resonances can be interpreted as members of the 1D family with J P = 1 − and 3 − . The resonances D s J (2632) , D * s1 (2710) and D s J (3040) are identified as radially exited states with n = 2 in our model. The D s J (2632) was firstly observed by SELEX Collaboration at a mass of 2632.5 ± 1.7 MeV, it can be assigned as the |2 1/2 S 0 . The assignment for D * s1 (2710) is proposed as J P = 1 − in Refs. [54,55], which agree with our prediction as our calculated mass for |2 1/2 S 1 is close to its experimental mass 2708 ± 9 +11 −10 MeV [5]. The D s J (3040) resonance is observed in the D * K mass spectrum at a mass of 3044±8 +30 −5 MeV by the BABAR Collaboration [6]. Here we assign it as |2 1/2 P 1 in our predicted D s meson spectrum.
In the b-flavored meson sector, experimental data for excited B meson states are limited for now. But still several b-flavored mesons are observed [57]. The strangeless resonances B J (5840) 0 and B(5970) 0 were measured by the LHCb and CDF Collaborations, respectively [10,11]. The stranged B * s J (5850) was observed by the OPAL Collaboration [12]. Their masses were measured as M(B J (5840)) = 5862.9 ± 5.0 ± 6.7 ± 0. In Table 2, we can identify B J (5840) 0 and B(5970) 0 as |1 1/2 P 1 and |2 1/2 S 1 , respectively, in the spectrum of B meson, while B * s J (5850) can be assigned as |1 1/2 P 1 in the spectrum of B s meson.
Finally, after solving the wave equation one can obtain not only the eigenenergy of each bound state but also their wave functions. The radial wave functions g n,l, j (r ) and f n,l, j (r ) for physical and unphysical D meson states are depicted as an example in Figs. 3 and 4, respectively. We stress that the solution of the eigenequation associated with the original H 0 in Eq. (26) gives only the wave functions of the physical states depicted in Fig. 3. In Sect. 3 we construct a new H 0 for the heavy-light systems in Eq. (43), the unphysical states depicted in Fig. 4 are due to the new H 0 for which the original one is substituted.

Summary
The spectra of heavy-light mesons are restudied in a relativistic model, which is derived by reducing the instantaneous Bethe-Salpeter equation. The kernel is chosen to be the standard combination of linear scalar and Coulombic vector. By applying the Foldy-Wouthuysen transformation on the heavy quark, the Hamiltonian for the heavy-light quarkantiquark system is calculated up to order 1/m 2 Q . We find that in the framework of an instantaneous Bethe-Salpeter equation the string tension b in the confinement potential is sensitive to the masses of the constituent quarks in the meson. The spectra of the D, D s , B and B s mesons are calculated in the relativistic model. Most of the heavy-light meson states can be accommodated successfully in our model except for the anomalous D * s0 (2317) and D s1 (2460) resonances. In the Bethe-Salpeter formalism, the assumption of the interaction kernel for mesons is rather a priori; kernels with other spin structures can also be studied. In this work, we only restrict our calculations to the spectra of heavy-light  [9] This work Previous work [27] Ref. [39] Ref.  Table 2 continued n j L J Meson E expt. [9] This work Previous work [27] Ref. [39] R e f.