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/m_Q^2$. Our results are in good agreement with available experimental data except for the anomalous $D_{s0}^*(2317)$ and $D_{s1}(2460)$ states. The newly observed heavy-light meson states can be accommodated successfully in the relativistic quark model with their assignments presented. The $D_{sJ}^*(2860)$ can be interpreted as the $|1^{3/2}D_1\rangle$ and $|1^{5/2}D_3\rangle$ states being members of the 1D family with $J^P=1^-$ and $3^-$.


I 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. While for D sJ mesons, besides the wellestablished 1S and 1P charmed-strange states, the excited resonances D sJ (2632) [3], D sJ (2860) [4], D sJ (2700) [5] and D sJ (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 and CDF Collaborations, respectively [10,11]. While the stranged B * sJ (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. In the 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 heavy-light 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 were widely used in studying mesons so as to embody the relativistic dynamics [21][22][23][24][25][26]. While 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 spectra 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 heavylight quark-antiquark system is expanded to order 1/m Q by applying 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 on the spectra of the heavy-light 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 only considered the leading 1/m Q term in the heavyquark 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 depend 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 III is for the solution of the wave equation and the peturbative corrections. In Section IV, we have numerical results and discussions. The last section is for a brief summary.

II 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]: 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 , the Eq.(1) can be simplified as: (4) 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 time coordinate [35][36][37].
In our previous work [27], with the help of projection operators for the wave function we find 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 Eq.(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 now is equivalent to a less complicated differential equation shown in Eq. (7), but still it is difficult to solve. For heavy-light 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 quark-antiquark 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 that 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) is 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 have: 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 of H ′ as in Eqs. (22∼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 (26) 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 heavylight quark-antiquark system in Eq. (26) is not Dirac-like as in Refs. [33,39]. Its form is more like the form used in relativized quark models [13,40,41]. As for doubleheavy 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.

III 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: (37) From Eqs. (35) and (37), we have: and since (h 2 ) 2 = 1, When we take c = −1, Eq. (35) is transformed to: which is not the correct eigenequation for bound systems we are interested in. Thus we have only c = +1. This is the reason for the substitution h 2 → 1 we use in the last section.
While if all the eigenfunctions of H 0 for bound states satisfy the relation the eigenfunction set of H 0 is NOT complete. An 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.(35) 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. (35). Now we turn to solve 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 commute 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 a 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 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 2j − l, respectively. The complete expression of y mj 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 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 n-th root of the spherical Bessel function j l (r), respectively.
Inserting Eq. (9) and (11) into Eq. (42), we can rewrite H 0 in the matrix form where h.c. stands for Hermitian conjugate and the oper-ator elements are: According to Eqs. (47) and (48), we can rewrite the eigenequation of H 0 in the representation of the state basis constructed by 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] Ω(p) j l (kr)Y lm (r) = Ω(k) j l (kr)Y lm (r), 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 can obtain the matrix elements of H 0 easily. For the operators with respect to energy of motion, we have In order to write down the expressions of the elements associated with the interaction potential in a compact form, here we introduce a symbolic notation: Ô m,lA;n,lB 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 are associated with the coefficients g i , f α , which are defined in Eqs. (47) and (48). That is to say, the eigenequation shown in Eq. (44) is solved and the corresponding eigenenergy and eigenfunction are obtained.
Here we turn to discuss 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 still it 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: 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 composed by the complete set of H 0 in Eq. (67), we obtain: The correction in first order perturbation can be written as: n,l,j + δE (1) n,l,j,J + δE (2) n,l,j,J , where δE (2) n,l,j,J = δE From Eq. (68), 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.

IV 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, α 1 = 0.25, α 2 = 0.15, α 3 = 0.20, and γ 1 = 1/2, γ 2 = √ 10/2, γ 3 = √ 1000/2, the α i and γ i parameters have the same values as in Ref. [13]. 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 confinement parameters b in this work is quite different from those in usual quark models. As discussed in our previous paper, the parameter b is responsible for elevating the energy levels of states with higher quantum numbers n or l, i.e. the calculated energy gaps between the energy levels of the radial or angular excitations and the ground state depend on the parameter b. But unlike the Dirac Hamiltonian: or the Schrödinger Hamiltonian in Eq. (47), the influence of the confinement potential V s (r) in Eq. (26) weakens as the mass of light quark decreases. That is to say, in the Bethe-Salpeter formalism the energy level is also sensitive to the light quark mass m q , which is shown in Fig. 2 by the energy gap ∆E between the first radial excitation and the ground state as a function of m q . We take B meson as an example to illustrate the dependence on the quark mass. The values of the parameters, which are fitted for B meson spectrum, are fixed except for the light quark mass of B 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 < 1GeV, which is the case for heavy-light 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 Our analysis suggests that in the Bethe-Salpeter formalism the string tension b depends on the masses of the quark and antiquark, especially the light one of them. Besides the potential parameters, four quark mass parameters are employed to fit the heavy-light meson spectra. With all the considerations above, our best fitting of the parameters gives the following values m u,d = 0.398 GeV, m s = 0.598 GeV, 0.390 GeV 2 for cq system, 0.421 GeV 2 for bq system, 0.300 GeV 2 for cs system, 0.316 GeV 2 for bs system, c = −0.320 GeV.
In Section III, two numerical parameters L and N are introduced in our calculation. In principle, if the distance L and the size of the expansion basis N are taken to ∞, we can obtain the exact solution of the wave equation. Our calculation shows the solution is stable when L > 5 fm, N > 50. In this work they are taken as L = 10 fm, N = 150, and the size of the matrix of H 0 in Eq. (54) is 300 × 300. The Gauss-Legendre quadrature is widely used in our numerical calculation since lots of integrals are involved. After strenuous computation, the meson spectra is obtained. The spectra of the heavy-light D, D s , B, B s mesons are fitted based on the data in PDG [9]. The numerical results for the spectra of D, D s mesons are presented in Table I, while B, B s mesons are presented in Table II. The calculated spectra are in good agreement with the experimental data. Our results are compared to the results of two other relativistic models [34,39], one is derived by quasipotential approach, and the other is obtained by reducing the Bethe-Salpeter vertex function.
The result in this work is improved comparing to our previous work [27]. Taking the mass difference between the pseudoscalar state and the vector state for example, as shown in Table I,  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 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, but it is more likely to be explained beyond naive quark model [42]. The masses of the two resonances predicted by constituent quark model are generally 100 ∼ 200 MeV higher than experiments [33,34,[43][44][45]. The mass of D * 0 , 2318±29 GeV is almost identical to the mass of D * s0 , 2317.8±0.6 GeV. It cannot be explained in 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 still it is not capable to explain the small mass difference of the two resonances.
As the 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 resonance may be D * s0 (DK) and D s1 (D * K) molecular states, while in Refs. [47][48][49], the 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 groundstate (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.
Besides the well-established 1S and 1P heavy-light meson states, The highly excited states are also calculated in the spectra and the newly observed highly excited meson states beyond 1P are identified in our model.
From Table I, one can see our results favor the measurement.
As for D s mesons, several states beyond 1P state have been observed, their masses and identifications are presented in the lower part of Table I. Recently, LHCb collaboration identifies D * sJ (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 I, 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 sJ (2632) , D * s1 (2710) and D sJ (3040) are identified as radially exited states with n = 2 in our model. The D sJ (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 sJ (3040) resonance is observed in the D * K mass spectrum at a mass of 3044 ± 8 +30 −5 MeV by 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 * sJ (5850) were observed by the OPAL Collaboration [12]. Their masses were measured as: In table II, 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 * sJ (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 Fig. 3 and 4, respectively.  We stress that the solution of the eigenequation associated with the original H 0 in Eq. (26) gives only the wavefunctions of physical states depicted in Fig. 3. In Section III we construct a new H 0 for heavy-light systems in Eq. (42), the unphysical states depicted in Fig. 4 is due to the new H 0 for which the original one is substituted.

V 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 heavy-light quark-antiquark system is calculated up to order 1/m 2 Q . We find that in the framework of 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 D, D s , B and B s mesons 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 in calculating the spectra of heavy light mesons. With the wave functions obtained when solving the wave equation, B and D decays can be studied in further researches.