Heavy-light mesons in a relativistic model

We study the heavy-light mesons in a relativistic model, which is derived from the Bethe-Salpeter equation by applying the Foldy-Wouthuysen transformation to the heavy quark. The kernel we choose is based on scalar confinement 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, Ds, B, Bs meson states are obtained. The spectra are calculated up to the order of 1/mQ, and wave functions are treated to leading order.


Introduction
Recently, great progress has been made measuring the spectra of heavy-light quark-antiquark systems, 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 have been propounded to identify the new discovered states and explain their properties [9][10][11][12][13][14]. Spectroscopy provides a powerful test of the theoretical predictions based on the quark model in the standard model.
Heavy-light quark-antiquark systems Qq play an important role in understanding the strong interactions between quarks and antiquarks. The double-heavy meson sector is rather mature, and nonrelativistic potential models have been proven extremely successful for the description for heavy quarkonia [15]. For heavy-light meson systems, however, 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 have been 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. In Ref. [21], an effective approach is developed by reducing the Bethe-Salpeter equation. Here we apply the Foldy-Wouthuysen transformation to 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 making further approximations. The Breit interaction is widely used to study quark-antiquark systems as a nonrelativistic or semirelativistic model [22][23][24]. 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 leading order. The parameters in the equations are determined by fitting the spectra of the heavy-light meson states presented by the Particle Data Group (PDG) [25]. The identifications of the newly observed heavylight meson states are given based on the theoretical prediction.
The organization of this paper is as follows. Section 2 is for the relativistic model and the effective Hamiltonian of the heavy-light quark-antiquark system. In Section 3, we solve the relativistic wave equation. Section 4 is forthe numerical results and discussion. In Section 5 we have a brief summary.

The model
The Bethe-Salpeter equation for the quark-antiquark system can be written as [26]: (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 [27], which can be reduced to that used in Ref. [28]. 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 in front of the 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. [29,30].
The projection operators are defined as with Applying the projection operators, the four coupled equations 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. [31], that is where If one makes the approximation

073101-2
Eq. (23) can be bought into the Breit form. The Breit equation can be used as a 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 heavy-light 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 If we perform the Foldy-Wouthuysen transformation directly on Eq. (34), the result can be awkward. What we want to achieve is a reduced and simple form of the Hamiltonian for the heavy-light quark-antiquark system.
Here we show that it is possible. Equation (27) is not equivalent to the instantaneous Bethe-Salpeter equation, i.e. the four coupled equations Eqs. (16)- (19), as pointed out in Refs. [32,33]. We find that our goal can be achieved by employing the constraint of Eqs. (18) and (19) which are dropped in Eq. (27). Subtracting Eq. (19) from Eq. (18), we get 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 φ, Applying the equivalent form of Eq. (37), Eq. (36) can be transformed to a new equivalent equation and it is easy to verify that the above equation is equivalent to the two coupled equations (36) and (37). Furthermore, it is equivalent to

073101-3
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 In the preceding 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 . For the case of the heavy-light system we focus on in this paper, however, the reduction results of the two schemes 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). Its left-hand side can be rewritten as: where 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) we can expand Eq. (51) to order 1/E The last termH 2 is of order V 2 v . Comparing to the other terms, its correction is suppressed by the coupling constant, so here we omit it. Considering the substitution β (1) → 1 in the Foldy-Wouthuysen transformation,H 1 can be reduced tõ thenH is transformed to 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.

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. [34,35], we prove first a relation.
Let Ω(p) be the pseudo-differential operator function and Ω(k) be the normal function, here p and k stand for the moduli of the 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 i k· r can be decomposed into a series of spherical harmonics where j l is the l-th order spherical Bessel function, and Y lm (ˆ r) is the spherical harmonic, which satisfies the normalization condition withˆ r 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:

073101-5
The above equation is true for an arbitrary function φ(r), which gives the orthogonality condition for the spherical Bessel function 2 π r 2 dk k 2 j l (kr)j l (kr ) = δ(r − r ).
Furthermore, if we take φ(r) = j l (k r), we can transform Eq. (73) to a simple form by using Eq. (75), that is 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 be of great help in solving the wave equation. Now we turn to solving the eigenequation of H 0 . It is easy to verify that the operators { j 2 , j z , K, S z } are a set of mutually commuting operators which com- . 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 with Thus the normalization condition of Eq. (79) 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 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 spinor part of Eq. (79) can be rewritten as by using the formula (here σ stands for σ (2) ) we have where the relations are used. Inserting Eq. (12) and (15) into Eq. (64), we can rewrite H 0 in the matrix form

073101-6
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, so 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. (86), 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. (83). 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. (95), 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), (88) and (89) with the normalization condition Eq. (82), we can get the matrix elements of H 0 easily.
Here we define a symbolic notation (108) Then we have 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. (104). Then the eigenequation associated with H 0 is solved and the eigenfunction is obtained. Now we discuss the perturbative corrections of H . The perturbative term H does not commute with any of the operators introduced in solving the eigenequation of H 0 , but still commutes with where J = j + S and P is the parity operator. Thus the eigenstates of the total Hamiltonian H = H 0 +H can be labeled by the set of quantum numbers {n, J, M J , P }.
With the help of the Clebsch-Gordan coefficients, we can compose the basis states with the quantum number set {n, k, j, J, M J } by combining the eigenstates of H 0 , then calculating the corrections and mixings in the obtained basis.
The correction in first order perturbation can be written as E n,l,j,J = E (0) n,l,j + 1 m 1 δE (1) n,l,j,J .
Mixing can happen between states with the 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 , with i, j = 1, 2, and the Ψ (0) i,j 's denote the basis states which are defined in Eq. (113).
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.

Numerical results 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 of giving a reasonable spectral structure for each meson system. In the framework of the 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 heavylight D, D s , B, B s mesons mainly based on the meson states presented in PDG [25]. The numerical results for the spectra of heavy-light mesons are presented in two tables. Table 1 is for D, D s mesons and Table 2 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 deviation of the mass of the D s0 * (2317) 0 resonance from experiment can be as large as 100 MeV. In the two tables, one can see that except for 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, but the dropped α 2 s term in the Hamiltonian can also make the result incomplete. In order to refine the present model and improve the calculation, the potentials in the kernel should be considered in more detail [36,37]. 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 have been observed by the LHCb collaboration in the mass region between 2500 and 2800 MeV as well as the region around 3000 MeV. Their masses are measured as [8]:  Fig. 1. The radial wave functions g n,l,j (r) and f n,l,j (r) for the D meson as an example. The wave functions are defined in Eq. (79). They are the radial part of the solution of the eigenequation associated with H0.
In Table 1, 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 ob-tained simultaneously when solving the wave equation. The radial wave functions g n,l,j (r) and f n,l,j (r) for D mesons are depicted in Fig. 1 as an example. [g 2 n,l,j (r)+ f 2 n,l,j (r)]r 2 is the probability density distributed along the quark-antiquark distance r. In Ref. [40], 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 [27]. Their wave functions have the same pattern as ours, though the Hamiltonian H 0 in our model looks very different from the Dirac-like Hamiltonian used by them.

Summary
We have constructed a relativistic model for heavylight quark-antiquark systems by studying the reduction of the instantaneous Bethe-Salpeter equation. The kernel we chose is based on scalar confinement and vector Coulomb potentials, and it shows a Coulombic behav-ior at short distance and a linear confining behavior at long distance. The bound states of the D, D s , B, B s mesons were studied using this model. The predictions of the spectra and the assignments for the newly observed charmed mesons were presented, and our results are in reasonable agreement with experimental measurements. The spectral structure was obtained in the relativistic model with only three potential parameters. We find that the parameter b in the confinement potential depends on the masses of the quarks in the framework of the Bethe-Salpeter equation. Theoretical deviations from experimental data appear mainly in the D s meson sector, especially for the D * s0 (2317) 0 and D s1 (2460) resonances. The discrepancy may be ascribed to the naive assumption for the kernel; the influence of other kernels with different spin structures can be studied in further research. The wave functions of each bound state can be obtained by solving the wave equation and used in the study of B and D decays.