Analytic approach to calculations of mass spectra and decay constants of heavy-light quarkonia in the framework of Bethe-Salpeter equation

This work is an extension of the work in [34] to ground and excited states of 0 ++ , 0 − + , and 1 −− of heavy-light ( cu, cs, bu, bs , and bc ) quarkonia in the framework of a QCD motivated Bethe-Salpeter equation (BSE) by making use of the exact treatment of the spin structure ( γ µ (cid:78) γ µ ) in the interaction kernel, in contrast to the approximate treatment of the same in our previous works [32, 34]), which is a substantial improvement over our previous works [32, 34]. In this 4 × 4 BSE framework, the coupled Salpeter equations for Qq (that are more involved than the equal mass ( QQ ) mesons) are ﬁrst shown to decouple for the conﬁning part of interaction, under heavy-quark approximation, and analyically solved, and later the one-gluon-exchange interaction is perturbatively incorporated, leading to their mass spectral equations. The analytic forms of wave functions obtained from these equations are then used for calculation of leptonic decay constants of ground and excited states of 0 − + , and 1 −− as a test of these wave functions and the over all framework. on


Introduction
During the past few years, there is a growing interest in the experimental and theoretical studies of heavy-light mesons. This interest arose from the discovery of large B 0 − B 0 mixing, leading to the hope that CP violation in B-systems may be observed. Further studies on heavy-light mesons are also important for the determination of Cabibo-Kobayashi-Maskawa (CKM) mass matrix elements. These studies on quarkonia need heavy quark dynamics, which can provide a significant test such as χ b0 (3P ), χ c0 (2P ), X(3915), X(4260), X(4360), X(4430), X(4660) [18]. Further, there are also open questions about the quantum number assignments of some of these states such as X(3915) (as to whether it is χ c0 (2P ) or χ c2 (2P ) [24,25]). Currently there is a lot of excitement about XY Z particles, that are new charmonium like states such as, Z c (3900), Z c (4020) /Z c (4025) [26], which were discovered in BESIII, the states, Y (4260) [27,28](discovered at BABAR), and X(3872) [29] (discovered at Belle). These particles show different features than the conventional charmonium states, and might be good candidates for exotic states, and might even be hybrid or tetra-quark states, or loosely bound charmonium molecules, which is one of the predictions of QCD.
Thus charmonium like states offer us intriguing puzzles. However, since the mass spectrum and the decays of all these bound states can be tested experimentally, theoretical studies on them may throw valuable insight about the heavy quark dynamics. Studies on these hadrons is particularly important, as it throws light on the long ranged QQ and Qq potential, that has not been derived from QCD so far.
In our works, we are not only interested in studying the mass spectrum of hadrons, which no doubt is an important element to study dynamics of hadrons, but also the hadronic wave functions that play an important role in the calculation of decay constants, form factors, structure functions etc. for QQ, and Qq hadrons. This is due to the fact that so far, one of the central difficulties in tests of QCD is lack of knowledge of hadronic wave functions. These hadronic Bethe-Salpeter wave functions calculated algebraically in this work can act as a bridge between the long distance nonperturbative physics, and the short distance perturbative physics. This is further due to the fact that, though these quarkonium states appear to be simple, however, their production mechanism is still not properly understood. Thus the wave functions calculated analytically by us can lead to studies on a number of processes involving QQ, and Qq states. Though the ground-state quarkonia have been shown over the past decade to be rather well described in nonrelativistic QCD, and N.R. potential models, the heavy-light mesons are much more complicated, where one really tests the wide range of soft, hard and soft-collinear scales. Our basic aim has been to develop a model using 4 × 4 BSE that can explain both mass spectrum of QQ, and Qq states as well as their decay widths through various processes using the same set of input parameters that are fixed from their mass spectrum.
In this context, in some of the recent works [30][31][32][33][34], we have been involved in working on the mass spectrum and decay properties of equal mass ground and excited states of scalar, pseudoscalar, vector, and axial vector QQ quarkonia in the framework of a 4 × 4 BSE. These include their leptonic decays, two-photon decays, single photon radiative decays and two gluon decays of these charmonium (cc) and bottomonium (bb) states which have been extensively studied by us in the formulation of Bethe-Salpeter equation. However, we had not so far generalized this 4×4 representation for two-body (QQ) BS amplitude framework to incorporate unequal mass dynamics, which we have now done in the present work. However, the price we have to again pay is to analytically solve a coupled set of equations for all quarkonia, which we have again explicitly shown get decoupled, in spite of the fact that we have used the full structure of the BS wave function, ψ( q) in calculation of γ µ ψ( q)γ µ on the right side of the Salpeter equations. Due to these facts, the system of coupled Salpeter equations encountered in the present work are much more involved and complex than the ones encountered in equal mass quarkonia in [34]. We have explicitly shown that they lead to mass spectral equations with analytical solutions for both masses, as well as eigen functions for the ground and excited states for 0 ++ , 0 −+ , and 1 −− for heavy-light hadrons with quark composition, cu, cs, cb, bu, and bs in an approximate harmonic oscillator basis. We then perturbatively incorporate the One-Gluon-Exchange (OGE), and solve the spectrum of these states. We wish to mention that in unequal mass systems such as Qq, the quarks are not very close together, and the confining interaction dominates over the OGE interactions due to which the perturbative incorporation of OGE term is a good approximation.
The analytical forms of eigen functions for ground and excited states obtained as analytic solutions of spectral equations are then used to evaluate the decay constants and decay widths for leptonic decays of these decays as a test of our framework.
The study of these mesons involves unequal mass kinematics. The unequal mass kinematics (that also gives the partitioning of internal momentum of hadron) used by us that rests on Wightman-Garding (W-G) definitions of momenta between individual quarks is relativistic, and has the advantage that P.q = 0 irrespective of whether the individual quarks are on-shell (P.q = 0) or off-shell (P.q = 0). This W-G partitioning of momenta between individual quarks is a natural choice of momentum partitioning, that allocates most of the internal momenta to heavier quark, while a smaller part of momentum to lighter quark, such that m 1 + m 2 = 1, while for equal mass mesons, the momentum is shared equally between the two quarks, which is what one expects.
The main advantage of our approach in comparison to other BSE approaches is that, we follow analytic methods of solutions for heavy-light quarkonia (whose equations are much more involved than QQ) , that provide a much deeper in sight into the mass spectral problem, and are able to obtain the mass spectrum in terms of the principal quantum number N , and also in the process, we get algebraic forms of wave functions that are used for calculations of various transition amplitudes and decay constants of quarkonia, in contrast to the purely numerical approaches followed by the other works. The correctness of our analytic approach can be judged by the fact that the plots of our wave functions (see [34]) are very similar to the plots of wave functions obtained by numerical approaches [14]. This paper is organized as follows: In section 2, we introduce the formulation of the 4 × 4 Bethe-Salpeter equation under the covariant instantaneous ansatz, and derive the hadron-quark vertex. In sections 3, 4, and 5, we derive the mass spectral equation of heavy-light scalar, pseudoscalar, and vector mesons respectively. In Sections 6, and 7, we derive the decay constants f P for pseudoscalar, and f V for vector Qq states respectively. In section 8, we provide the numerical results and discussion.

Formulation of the 4× 4 Bethe-Salpeter equation
We give here the main points about the 4× 4 BSE under the Covariant Instantaneous Ansatz (CIA), which is a Lorentz-invariant generalization of Instantaneous Approximation (IA), which is used to derive the 3D Salpeter equations [32,34,35]. We start with a 4D BSE for quark-anti quark system with quarks of constituent masses, m 1 and m 2 , written in a 4 × 4 representation of 4D BS wave function Ψ(P, q) as: where K(q, q ) is the interaction kernel between the quark and anti-quark, and p 1,2 are the momenta of the quark and anti-quark, which are related to the internal 4-momentum q and total momentum P of hadron of mass M as, p 1,2µ =m 1,2 P µ ±q µ , where,m 1, ], always satisfy,m 1 +m 2 = 1, and is a natural choice that allocates most of the momentum to the heavy quark, while a smaller part of momentum to the lighter quark in a heavy-light meson, but equal momenta to both quarks in cc mesons.
Making use of Covariant Instantaneous Ansatz, where, K(q, q ) = K( q, q ) on the BS kernel, where q µ = q µ − q.P P 2 P µ is the component of internal momentum of the hadron that is orthogonal to the total hadron momentum, i.e. q.P = 0., while σP µ = q.P P 2 P µ is the component of q longitudinal to P , where the 4-dimensional volume element is, d 4 q = d 3 qM dσ, and following a sequence of steps outlined in [32], we get the covariant forms of four Salpeter equations (in 4D variable q), which are effective 3D forms of BSE, and are valid for hadrons in arbitrary motion. The four independent Salpeter equations are [32]: where Λ ± are the projection operators [32] for each of the constituents. Γ( q) is the 4D hadron-quark vertex function, which enters into the 4D BS wave function, Ψ(P, q) = S F (p 1 )Γ(q)S F (−p 2 ), where We wish to emphasize that the present model after 3D reduction is still covariant. This is due to the fact that we have reduced a fully 4D BSE to 3D BSE (which are actually four Salpeter equations) by use of Covariant Instantaneous Ansatz (CIA), which is a Lorentz-invariant generalization of the Instantaneous Approximation (IA). We thus obtain the covariant forms of Salpeter equations, which are effective 3D forms of BSE, and are valid for hadrons in arbitrary motion.
Regarding the interaction kernel K(q ,q) [34], it can be written as, with colour, spin and orbital parts respectively. For a kernel with the above spin dependence, we can rewrite the hadron-quark vertex in Eq.(3) as [34], where, each of the γ µ s sandwich the BS wave function, ψ( q), with the scalar part of the kernel, where, the confinement part with a sequence of steps can be expressed as The present work is a substantial improvement over our previous works [32,34], in the sense that we have taken the full Dirac structure of the 3D BS wave function, ψ( q) given in Eq. (7) to calculate the spin part, γ µ ψ( q)γ µ that enters into the hadron-quark vertex function, Γ(q) as well as the right hand sides of the 3D coupled integral equations, in contrast to [32,34], where we took only the leading Dirac structures in ψ( q) to evaluate γ µ ψ( q)γ µ , in the integrals on the right of the Salpeter equations in [32,34]. In this work, we have obtained the mass spectral equations with this exact treatment of the spin part γ µ Ψ( q)γ µ . What we further find is that the higher order terms of V c that we had ignored in [32,34] due to negligible coefficients, ω 4 qq associated with these terms, get effectively cancelled out when we take the full Dirac structure of the wave function, ψ( q).
Further, in the present work, we notice that with the use of the exact treatment of the spin dependent part of the kernel in the RHS of Salpeter equations, for case m 1 = m 2 , we get the mass spectrum of equal mass quarkonia (χ c0 , η c , and J/Ψ), for both ground and excited states, where the excited states are closer to data [19] than the excited states obtained in our previous works [32,34].
The framework is quite general so far. Thus, to obtain the mass spectral equation, we have to start with the above four Salpeter equations to solve the instantaneous BS equation.

Mass spectral equation for heavy-light scalar 0 ++ quarkonia
We start with the general form of 4D BS wave function for scalar meson (0 ++ ) in [36]. Then, making use of the 3D reduction and making use of the fact that q.P = 0, we can write the general decomposition of the instantaneous BS wave function for scalar mesons (J pc = 0 ++ ), of dimensionality M being composed of various Dirac structures that are multiplied with scalar functions f i (q) and various powers of the meson mass M as [34] Till now these amplitudes f 1 , and f 4 in equation above are all independent, and as per the power counting rule [11,12] proposed by us earlier, the f 1 , and f 2 are the amplitudes associated with the leading Dirac structures, namely M and P , while f 3 and f 4 will be the amplitudes associated with the sub-leading Dirac structures, namely, q, and 2 P q M . We now use the last two Salpeter equations ψ +− (q) = ψ −+ (q) = 0 in Eq. (2), that can be used to obtain the constraint relations between the scalar functions for unequal mass mesons as The BS-wave function for scalar mesons in Eq. (7) with the help of these constraint relations can be rewritten in terms of only two independent scalar functions (f 1 (or f 3 ) and f 4 ) as We wish to mention that due to these two above equations, the scalar functions f i ( q)(i = 1, ..., 4) are no longer all independent, but are tied together by these relations, due to which the amplitudes get mixed up [34].
The first two Salpeter equations of Eq.(2) lead to a set of coupled integral equations with the full structure of the wave function ψ S (q) in (9) being used to evaluate γ µ ψ S (q)γ µ on the right hand sides of these equations. We proceed in the same way as, [34], where on the right side of these equations, we first work with the confining interaction, V c ( q). We show that these equations can be decoupled, and reduced to algebraic equations in an approximate harmonic oscillator basis, and solve them analytically. These equations are much more involved than the equal mass case [34]. We then incorporate the one-gluon-exchange (OGE) term perturbatively, and obtain the complete spectrum.
These coupled equations (with use of confining interaction alone) are, To decouple these equations, we follow the same procedure in [32,34], where we first add them.
Then we subtract the second equation from the first equation. For a kernel that can be expressed as , we get two algebraic equations which are still coupled. Then from one of the two equations so obtained, we eliminate f 3 (q) in terms of f 4 (q), and plug this expression for f 3 (q) in the second equation of the coupled set so obtained to get a decoupled equation in f 4 (q). Similarly, we eliminate f 4 (q) from the second equation of the set of coupled algebraic equations in terms of f 3 (q), and plug it into the first equation to get a decoupled equation entirely in f 3 (q), which reduces to two identical decoupled equations, one entirely in f 3 (q), and the other that is entirely in f 4 (q) as: It is to be seen here that on RHS of the above two equations, with the exact treatment of the spin part of the kernel, we get only the terms that are linear in V c , (unlike [32,34], where we also obtained quadratic terms of the type, V 2 c , that were very small in magnitude in comparison to V c ). Since the two equations are of the same form in scalar functions f 3 (q) and f 4 (q), that are the solutions of identical equations, we can take, f 3 (q) ≈ f 4 (q)(= φ S (q)). Using the expression for V c (q) given above, we get the equation, where the inverse range parameter β S can be expressed as, Using the method of power series, this leads to the mass spectral equation for scalar mesons as, with the energy eigen value of the scalar mesons, E S = 2β 2 S (N + 3 2 ), where N = 2n + l, with the principal quantum number taking values n = 0, 1, 2, ..., and the orbital quantum number l = 1 that corresponds to P wave states, and the solutions of Eq.(12) are given by the following normalized wave functions that are similar to the wave functions in [34], except for the inverse range parameter β expression that is different from [34] due to the exact treatment of the spin part of the kernel, and also the unequal mass kinematics. The overall structure of these wave functions is very similar to the wave functions derived in [34], except for the algebraic form of β S . They are: Now, we treat the mass spectral equation in Eq. (12), which is obtained by taking only the confinement part of the kernel, as an unperturbed spectral equation with the unperturbed wave functions in Eq. (15). We then incorporate the one gluon exchange term in the interaction kernel perturbatively (as in [34]) and solve to first order in perturbation theory. The complete mass spectra of ground and excited states of heavy-light scalar mesons is where V S coul is the expectation value of V S coul between the unperturbed states of the scalar mesons with l = 1 and n = 0, 1, 2, ..., and γ = The results of our model for mass spectrum for scalar Qq states along with data [19], and other models is given in Table   M χ c0 (4P0) 4.7514 Table 1: Masss spectra of ground and excited states of scalar 0 ++ quarkonia (in GeV) We now derive the mass spectral equations of unequal mass pseudoscalar mesons in the next section.

Mass spectral equations for heavy-light pseudoscalar 0 −+ quarkonia
The general decomposition for the 3D wave function of pseudoscalar mesons obtained from the general 4D form [36] through 3D reduction as in previous section can be written as [34] ψ We use the last two Salpeter equations in Eq.
(2) to find the constraints on the components of the wave function as Plugging Eq. (19) into Eq. (18), we rewrite the wave function for pseudoscalar mesons as We use the first two Salpeter equations of Eq.(2) to obtain the corresponding coupled integral equations of pseudoscalar mesons (with use of confining interaction alone) as, Using the same procedure as in the case of scalar mesons, these two equations can be decoupled, and reduced to two independent algebraic equations as Here, we again see that the scalar functions φ 1 (q) and φ 2 (q) satisfy identical equations, and can be taken as φ 1 (q) ≈ φ 2 (q)(= φ P (q)). Using the expression for V c (q) after Eq. (6), we obtain the mass spectral equation as, whose solutions give the unperturbed mass spectrum (due to confining interactions alone), with the orbital quantum number l = 0 that corresponds to the S states, and β P = . This unperturbed mass spectral equation of pseudoscalar meson is the same as the corresponding spectral equation of scalar meson in Eq. (12), except that β S is replaced by β P , and φ S (q) replaced by φ P (q).
The normalized unperturbed wave functions of 1S, ..., 4S states of pseudoscalar meson with l = 0 are We again incorporate the Coulomb term V P coul associated with the one gluon exchange interaction perturbatively into the original mass spectral equation of pseudoscalar mesons, giving us the complete mass spectra of ground and excited states of heavy-light pseudoscalar quarkonia with orbital quantum number l = 0 as where again the perturbation parameter γ has the same form as in the case of scalar mesons with β 2 S replaced by β 2 P and γ = ω 4 0 C 0 β 2 P (3l + 1), while, the first order correction to the total energy of the system E P is given by the expectation value of the Coulomb term between the unperturbed states of pseudoscalar mesons φ P (nS,q) as The results of our model for pseudoscalar Qq mesons along with data [19] and other models is given in Table 2.
We now give the derivation of the mass spectral equations of vector mesons in the next section.

Mass spectral equations for heavy-light vector 1 −− quarkonia
We again start with the general 4D decomposition [36]. Using 3D decomposition, the wave function of vector mesons can be written as [32,34]: The constraint equations on the components of the wave functions (χ s) can be obtained using the last two Salpeter equations of (2) as χ 3 (q) = χ 6 (q) = 0 Substituting Eq.(29) into Eq.(28), the wave function for vector mesons can be rewritten as Using the first two Salpeter equations, we obtain the coupled integral equations of vector mesons  (with confining interaction alone) as Now, using the same procedures to decouple these equations as in the case of scalar and pseudoscalar mesons, we obtain two decoupled algebraic equations, Here, we see that the scalar functions χ 1 (q) and χ 2 (q) satisfy identical equations, and can be taken We then obtain a single differential equation, which is nothing but the equation of a simple quantum mechanical 3D-harmonic oscillator with coefficients depending on the hadron mass M , and total quantum number N . The wave function satisfies the 3D BSE: which can be rewritten as where is the inverse range parameter, and the total energy of the system is identified as This mass spectral equation of vector meson is the same as the corresponding equation of scalar meson in Eq. (12), except that β S is replaced by β V , and φ S (q) replaced by φ V (q). Therefore, the normalized wave functions of 1S,...,3D states of vector meson, with S and D states corresponding to l = 0 and l = 2 respectively, are Eqs. (34)(35) would lead to degenerate masses for S and D states of Qq, and QQ systems. To get S− D mass splitting, we make use of degenerate perturbation theory. The Coulomb term V V coul associated with the one gluon exchange interaction is perturbatively incorporated into the mass spectral equation, Eq.(34) (that is treated as the unperturbed equation) for vector mesons, as: The complete mass spectral equation of heavy-light vector quarkonia can be put as where V V coul has been weighted by the perturbation parameter γ = the only non-zero expectation values of V V coul are the ones that connect states of the same quantum numbers, n and l. They are: We now give the plots of these normalized wave functions Vs. q (in Gev.) for different states of heavy pseudoscalar and vector mesons ( such as composite of cu, cs, cc, cb and bb) in Fig.1-2 Fig.3-4, respectively. It can be seen from these plots that the wave functions corresponding to nS and nD states have n − 1 nodes. We now calculate the leptonic decays of pseudoscalar and vector Qq mesons in the next sections.

Leptonic decays of pseudoscalar heavy-light quarkonia
The leptonic decays of pseudoscalar quarkonia proceed through the coupling of the quark-anti quark loop to the axial vector current. The leptonic decay constants, f P are defined as, The decay constants can be expressed through the quark-loop integral as, Here, the full 3D BS wave function, ψ P (q) can be taken from Eq. (20), where φ 1 (q), and φ 2 (q) satisfy two identical decoupled equations, Eq.(22), leading to φ 1 (q) = φ 2 (q)(= φ P (q)), which is expressed as, Putting the above expression for ψ P in Eq. (41), and evaluating trace over the gamma matrices on the right side of the equation, we get, To evaluate f P , we multiply both sides of the above equation by Pµ M 2 , and making use of the fact that q.P = 0, and the fact that, we can express f P in terms of a 3D integral, where the 3D wave functions, φ P (q) for pseudoscalar Qq states are given in Eqs. (25), and N P is the 4D BS normalizer obtained through the current conservation condition, where Ψ P ( q) is the 4D BS wave function, while the adjoint BS wave function, Ψ(P, q) = γ 4 ψ † (P, q)γ 4 .
Making us of the fact that in the inverse propagators, S −1 F (p 1,2 ) of the two quarks, their momenta are expressed as, p 1,2 = m 1,2 P ± q, and the 4D volume element, d 4 q = d 3 qM dσ. Integrating over M dσ, and making us of the 3D form of BS wave function, ψ( q), in Eq. (42), evaluating trace over the gamma matrices, and multiplying both sides of the equation by P µ , we get, where we take into account the contribution of the second term, in Eq. (46) to be the same as the contribution of the first term. Leptonic decay constants of 0 −+ quarkonia are given in Table 4 along with data and results of other models.

Leptonic decays of heavy-light vector quarkonia
Leptonic decays of vector quarkonia are defined through the equation, The decay constant f V can be expressed through the quark loop integral, Here ψ V (q) is the full 3D wave function for vector mesons that can be taken from Eq.(30), where χ 1 (q), and χ 2 (q) satisfy two identical decoupled equations, Eq.(32), leading to  is expressed as, where the 4D BS normalizer, N V can be obtained from the current conservation condition in Eq. (46), and following a similar procedure as in the case of pseudoscalar quarkonia as The leptonic decay constants of heavy-light vector mesons are given in Table 5.

Results and Discussion
We have employed a 3D reduction of BSE (with a 4×4 representation for two-body (qq) BS ampli- This feature in our framework is yet to be explored in detail, since in a Poincare covariant framework, [37,38], the numerical results for the amplitudes and masses are independent of the choice of momentum partitioning parameters. This detailed study we intend to do next.
In this work, we make use of the exact treatment of the spin structure (γ µ γ µ ) in the interaction kernel, in contrast to the approximate treatment of the same in our previous works [32,34]). In so doing we do away with the approximation of taking the leading Dirac structures in the structure of 4D BS wave function, Ψ(P, q), which is a substantial improvement over our previous works. We thus first derive analytically the mass spectral equation using only the confining part of the interaction kernel for Qq systems, and calculate the algebraic forms of the wave functions. Then treating this mass spectral equation as the unperturbed equation, we introduce the One-Gluon-Exchange (OGE) perturbatively, and obtain the mass spectra for various states of 0 ++ , 0 −+ , and 1 −− , treating the wave functions derived above as the unperturbed wave functions.
As mentioned earlier, in our works, we are not only interested in studying the mass spectrum of hadrons, which no doubt is an important element to study dynamics of hadrons, but also the hadronic wave functions that play an important role in the calculation of decay constants, form factors, structure functions etc. for QQ, and Qq hadrons, and so far, one of the central difficulties in tests of QCD is lack of knowledge of hadronic wave functions. These hadronic Bethe-Salpeter wave functions calculated algebraically in this work can act as a bridge between the long distance non-perturbative physics, and the short distance perturbative physics. And since these quarkonia are involved in a number of reactions which are of great importance for study of Cabibbo-Kobayashi-Maskawa (CKM) matrix and CP violation, the wave functions calculated analytically by us can lead to studies on a number of processes involving QQ, and Qq states. In this work, we have used these algebraic forms of wave functions to calculate the leptonic decay constants of pseudoscalar and vector Qq quarkonia in Tables   4 and 5 respectively.
We have first obtained the numerical values of masses for ground and excited states of various heavy-light mesons and made comparison of our results with experimental data and other models. In the process, we were also able to reproduce reasonable results for the mass spectrum of ground and excited states of charmonium (cc) such as η c , χ c0 , and J/ψ, which are found to be closer to data than the mass spectrum we had found in earlier works of equal mass quarkonia [32,34]. We then obtained the numerical values of leptonic decay constants for these pseudoscalar and vector heavy-light quarkonia with the same set of input parameters fixed from the mass spectrum.
All numerical calculations have been done using Matlab. We selected the best set of 8 input parameters that gave good matching with data for masses of ground and excited states of heavy-light scalar, pseudoscalar, and vector quarkonia. This input parameter set was found to be C 0 = 0.139, has been introduced to make the linear and harmonic terms of the interaction kernel dimensionally consistent, is chosen to have the form γ = ω 4 0 C 0 β 2 (3l + 1) (where β 2 = β 2 S,P,V ) in order to produce reasonable non-degenerate masses of 2S and 1D, 3S and 2D, 4S and 3D, etc. states of heavy-light quarkonia. We have thus obtained results of masses for B c , B s , B, D s , D, and η c , χ c0 , J/ψ, and are in reasonable agreement with experimental data and other models. We will be using the analytical forms of eigen functions for ground and excited states of heavy-light quarkonia to evaluate the various transition processes involving heavy-light scalar, pseudoscalar, and vector quarkonia as future work.