An approach to calculation of mass spectra and two-photon decays of $c\overline{c}$ mesons in the framework of Bethe-Salpeter Equation

In this work we calculate the mass spectrum of charmonium for $1P,...,3P$ states of $0^{++}$ , $1^{++}$, as well as for $1S,...,4S$ states of $0^{-+}$, and $1S, ...,5D$ states of $1^{--}$ along with the two-photon decay widths of ground and first excited state of $0^{++}$ quarkonia for the process, $O^{++}\rightarrow \gamma\gamma$ in the framework of a QCD motivated Bethe-Salpeter Equation. In this $4\times 4$ BSE framework, the coupled Salpeter equations are first shown to decouple for the confining part of interaction, under heavy-quark approximation, and analyically solved, and later the one-gluon-exchange interaction is perturbatively incorporated leading to mass spectral equations for various quarkonia. The analytic forms of wave functions obtained are used for calculation of two-photon decay widths of $\chi_{c0}$. Our results are in reasonable agreement with data (where ever available) and other models.


Introduction
An important role in applications of Quantum Chromodynamics (QCD) to hadronic physics is played by charmonium (cc) and bottomonium (bb), which are built up of a heavy quark and heavy anti-quark. By 1 arXiv:1610.03234v5 [hep-ph] 13 Oct 2017 definition, heavy quark has a mass m, which is large in comparison to the typical hadronic scale, Λ QCD . Quarkonia are characterized by at least three widely separated scales [1]: the hard scale (the mass m of heavy quarks), the soft scale (relative momentum q ∼ mv), and the ultra soft scale (typical kinetic energy, E ∼ mv 2 of heavy quark and anti-quark). The appearance of all these scales in the dynamics of heavy quarkonium makes its quantitative study extremely difficult. Quarkonium systems are crucially important to improve our understanding of QCD. They probe all the energy regions of QCD from hard region, where perturbative QCD dominates, to the low energy region, where non-perturbative effects dominate.
Quarkonium states are thus a unique laboratory where our understanding of non-perturbative QCD may be tested.
There has a been a renewed interest in recent years in spectroscopy of these heavy hadrons in charm and beauty sectors, which was primarily due to experimental facilities the world over such as BABAR, Belle, CLEO, DELPHI, BES etc. [2,3,4,5,6], which have been providing accurate data on cc, and bb hadrons with respect to their masses and decays. In the process many new states have been discovered such as χ b0 (3P ), χ c0 (2P ), X(3915), X(4260), X(4360), X(4430), X(4660) [6]. The data strongly suggests that among new resonances may be exotic four-quark states, or hybrid states with gluonic degrees of freedom in addition to cc pair, or loosely bound states of heavy hadrons, i.e. charmonium molecules. 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 ) [7,8]). Thus charmonium offers us intriguing puzzles.
However, since the mass spectrum and the decays of all these bound states of heavy quarks can be tested experimentally, theoretical studies on them may throw valuable insight about the heavy quark dynamics and lead to a deeper understanding of QCD further. Studies on mass spectrum of these hadrons is particularly important, since it throws light on the QQ potential, since the long range confinement potential can not be derived from QCD alone. Further, though these states appear to be simple, however, their production mechanism is still not properly understood. These mesons are involved in a number of reactions which are of great importance for study of Cabibbo-Kobayashi-Maskawa (CKM) matrix and CP violation. In this paper we also study the two-photon decays of scalar quarkonia, χ c0 . These decays are sensitive probes of quarkonium wave functions.
The non-perturbative approaches, such as Effective field theory [9], Lattice QCD [10,11,12], Chiral perturbation theory [13], QCD sum rules [14,15], N.R.QCD [16], Bethe-Salpeter equation (BSE) [17,18,19,20,21,22,23,24,25,26], and potential models [27,28,29] deal have been employed to study heavy quarkonia. Recent progress in understanding of non-relativistic field theories make it possible to go beyond phenomenological models, and for the first time face the possibility of providing unified description of all aspects of heavy quarkonium physics. This allows us to use quarkonium as a bench mark for understanding of QCD, and for precise determination of Standard Model parameters (e.g. heavy quark masses, QCD coupling constant, α s ), and for new physics searches.
All this opens up new challenges in theoretical understanding of heavy hadrons and also provide an important tool for exploring the structure of these simplest bound states in QCD and for studying the non-perturbative (long distance) behavior of strong interactions.
In present work, we do the full mass spectral problem of charmonium for 1P, ..., 3P states of 0 ++ , 1 ++ , as well as for 1S, ..., 4S states of 0 −+ , and 1 −− along with the two-photon decay widths of ground and first excited state of scalar quarkonia, χ c0 for the process, O ++ → γγ in the framework of a QCD motivated Bethe-Salpeter Equation under Covariant Instantaneous Ansatz (CIA), employing the full BS kernel comprising both, the long-range confinement, and the short range one-gluon-exchange (coulomb) interactions.
We do understand that in QQ quarkonia, the constituents are close enough to each other to warrant a more accurate treatment of the coulomb term. Though for bb systems, the coulomb term will be extremely dominant in comparison to confining term, and should not be treated perturbatively. However, seeing our mass spectral results for cc systems, it may not be so unreasonable to treat the coulomb term perturbatively for cc systems. This is specially so for orbital excitations of these states, where the centrifugal effects [34] ensure that the c − c separation is large enough to feel the effect of confining term more strongly than the coulomb term. We further wish to state that some of earlier works [34,35,36] have treated the OGE(coulomb) term perturbatively for charmed mesons and baryons, while some works [37,38] did not take into account the importance of coulomb term for heavy quarkonium systems.
Thus, in the present paper, the coupled Salpeter equations for scalar (0 ++ ), and axial vector (1 ++ ) quarkonia are first shown to decouple for the confining part of interaction, under heavy-quark approximation, and the analytic forms of mass spectral equations are worked out, which are then solved in approximate harmonic oscillator basis to obtain the unperturbed wave functions for various states of these quarkonia. We then incorporate the one-gluon-exchange perturbatively into the unperturbed spectral equation, and obtain the full spectrum. The wave functions of scalar (0 ++ ) quarkonia, are then used to calculate their two-photon decay widths. We further extend the mass spectral calculations of pseudoscalar (0 −+ ), and vector (1 −− ) quarkonia in [30] with the perturbative inclusion of one-gluon-exchange effects.
The approximations used in this analytic treatment of the confining interaction are shown to be fully under control. This work is an improvement on our earlier work [30] on mass spectral problem for pseudoscalar and vector states of quarkonia on lines of some of the earlier works [31,32,33], where we used only the confining interaction.
A quarkonium state is classified by quantum numbers, J P C , where J = L + S, parity, P = (−1) L+1 , and charge conjugation, C = (−1) L+S . With this classification, while the lowest states, l = 0 are present in 0 −+ (pseudoscalar), the lowest state (l = 0), and the second orbitally excited (l = 2) state are present in 1 −− (vector) quarkonia. However, the first orbitally excited (l = 1) states are present in 0 ++ (scalar) and 1 ++ (axial-vector) quarkonia. The same holds true for their radial excitations. This paper is organized as follows: Section 2, deals with the mass spectral calculations of ground and excited states of 0 ++ quarkonia. Section 3 deals with the mass spectral calculation of ground and excited states of 0 −+ , and 1 −− quarkonia, while Section 4 deals with the mass spectral calculations of ground and excited states of 1 ++ quarkonia. Section 5 deals with the calculation of two-photon decay widths of χ c0 . Section 6 deals with numerical results and discussions.

Mass spectra of scalar quarkonia
In the center of mass frame, where q µ = (q, i0), we can write the general decomposition of the instantaneous BS wave function for scalar mesons (J pc = 0 ++ ), of dimensionality M as where it can be shown by use of a power counting rule proposed in [23,25] that the Dirac structures associated with the amplitudes f 1 and f 2 are leading, and will contribute maximum to calculation of any scalar meson observable, while those associated with f 3 , and f 4 are sub-leading. We now use the two constraint equations ψ +− (q) = ψ −+ = 0 in the 3D Salpeter equations, Eq.(16) of [30], to reduce the four scalar functions into two independent scalar wave functions. For equal mass system, these two equations are reduced to Applying the above constraint conditions in Eq.(2) to wave function in Eq.(1), we rewrite the relativistic wave function of the state (0 ++ ) in the form Here, it is to be noted that by use of the above constraint equations, we have reexpressed f 1 in terms of f 3 . Thus, the Instantaneous BS wave function of (0 ++ ) state is determined by only two independent functions, f 3 and f 4 . Putting the wave function in Eq.(3) above, along with the projection operators defined in Eq.(13) of [30], into the first two Salpeter equations, Eq.(16) of [30], and by taking the trace on both sides, we obtain the equations: Solution of these equations needs information about the BS kernel, K(q,q ) [25,30], which is taken to be one-gluon exchange like as regards the colour ( 1 4 λ 1 . λ 2 ) and spin (γ µ γ µ ) dependence, and has a scalar part V (q,q ), written as, Thus, the scalar part V (q,q ) of the kernel involves both the OGE term V OGE , arising from the onegluon exchange, as well as the confining term, V c . We first ignore the V OGE term, and work only with the confining part, V c in Eq. (5), and write Eqs.(4) as, where the spin dependence of the interaction is contained in the factor, Θ S = γ µ ψ( q)γ µ . The scalar part of the confining potential (that involves the colour factor 1 2 λ 1 . Here V c is the part of V c without the delta function. To handle these equations, we first integrate these equations over d 3 q , performing the delta function integration that arises due to presence of V c (q,q ) in the integrand, and get two coupled algebraic equations with V ( q) on RHS. To decouple them, we first add them. Then we subtract the second equation from the first equation, and 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). Thus, we get two identical decoupled equations, one entirely in f 3 (q), and the other that is entirely in, f 4 (q).
The calculation up to this point is without any approximation. However, we notice that if we employ the approximation, ω ≈ m on RHS of the two algebraic equations 1 (this is justified since in the confinement region, the relative momentum between heavy quarks in the bound state can be considered small, since the heavy-quark is expected to move with non-relativistic speeds (see Ref. [39]), and these quarks can be treated as almost on mass shell), these decoupled equations can be expressed as: We see that we get two identical decoupled equations which resemble the harmonic oscillator equations, but for the term involving V 2 c on the right side of these equations. We wish to mention that in recent studies on mass spectra of pseudoscalar and vector quarkonia [30], it was seen that good agreement with data on masses and various decay constants/decay widths of ground and excited states of η c , η b , J/Ψ, and Υ is obtained for input parameters,  Table 1 below, where it can be seen that ω 4 qq ω 2 qq and hence the second term on RHS of Eqs.(7) contributes < 1%, than the first term on RHS of this equation for cc, so that it can be dropped. Table 1: Numerical values of coefficients, Ω S = mΘ S ω 2 qq , and Ω S = Θ 2 S 4 ω 4 qq associated with the terms involving V c and V 2 c respectively, for scalar mesons χ c in RHS of Eqs. (7), and their percentage ratio for the input parameters of our model mentioned above.
Thus the RHS of both equations in Eq.(7) has only the term, −mΘ s V c f 3,4 ( q). Now, putting the spatial wave functions that can be employed to do analytic calculations of various transition amplitudes for different processes. Our approach may lead to a little loss of numerical accuracy, but it does lead to a much deeper understanding of the mass spectral problem. The approximations used by us have been shown to be totally under control. The plots of our algebraic forms of wave functions for scalar, pseudoscalar, vector and axial vector quarkonia are very much similar to the corresponding plots of wave functions in [20] obtained by purely numerical methods, which validates our approach.
part, V c in Eq. (7), the wave functions f 3 , and f 4 then satisfy identical 3D BSE for equal mass heavy scalar mesons: Thus the solutions of these equations, f 3 (q) ≈ f 4 (q)(= φ s (q)). With the use of the above equality of amplitudes, and reexpressing f 3 in terms of f 1 , the complete wave function Ψ s ( q) can be expressed as: We can reduce this equation into 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: Now, with the use of leading Dirac structures, we can to a good approximation (as in case of pseudoscalar and vector mesons [30]) express Θ S = 4. This is due to the fact that γ µ ψ( q)γ µ ≈ γ µ (M f 1 )γ µ = 4ψ(q), due to M I (I being the unit 4 × 4 matrix) being the most leading Dirac structure in the scalar meson wave function, and the terms with q 2 /m 2 have negligible contributions, in heavy quark limit, and can be dropped. Thus the above equation can be put in the form which can in turn be expressed as, where, β s = ( ) 1/4 , and with total energy of the system expressed as, Putting the expression for the laplacian operator in spherical coordinates, we get where l is the orbital quantum number with the values l = 0, 1, 2, 3..., corresponding to S, P, D, ... wave states respectively. This is a 3D harmonic oscillator equation, whose solutions can be found by using power series method. Assuming the form of the solutions of this equation as, φ(q) = h(q)e −q 2 2β 2 s , the above equation can be expressed as, The eigen values of this equation can be obtained using the power series method as, with l = 1. Thus, to each value of n = 0, 1, 2, 3, ... will correspond a polynomial h(q) of order 2n + 1 in q, that are obtained as solutions of Eq. (15). The odd parity normalized wave functions φ(q) thus derived are: The plots of wave functions for scalar quarkonia, χ c0 are given in Fig. 1 below.
where the orbital quantum number, l = 1. Now, treating the mass spectral equation, Eq.(12) as the unperturbed equation, with the unperturbed wave functions, φ S (nP, q) for scalar mesons in Eq. (17), we now incorporate the OGE (coulomb) term in this equation.
Then the above mass spectral equation can be written as: Treating the coulomb term as a perturbation to the unperturbed mass spectral equation, we can write the complete mass spectra for the ground and excited states for equal mass heavy scalar (0 ++ ) mesons using first order perturbation theory as: with l = 1, where < V S coul > is the expectation value of V S coul between the unperturbed states of given quantum numbers n (with l = 1) for scalar mesons, and has been weighted by a factor of γ = C 2 0 β 4 S (m 1 +m 2 ) 2 to have the coulomb term dimensionally consistent with the harmonic term. Its expectation values for the 1P, 2P , and 3P states are:  Table 2: Mass spectrum of ground and excited states of χ c0 with quantum numbers J P C = 0 ++ in GeV.
units with the above set of parameters.
From the mass spectral equation, one can see that, the mass spectra depends not only on the principal quantum number N , but also the orbital quantum number l. We are now in a position to calculate the numerical values for mass spectral of heavy equal mass scalar meson with the input parameters of our model. The results of mass spectral predictions of heavy equal mass scalar mesons for both ground and excited states with the above set of parameters is given in table 2.
We now derive the mass spectral equations with the incorporation of the OGE (Coulomb) term for pseudoscalar and vector quarkonia, and obtain their solutions in the next section (the preliminary calculations using only the confining part of interaction were done in [30]). of [30] given as, with non-perturbative energy eigen functions for l = 0(S), and for l = 2(D) states obtained as solutions of above spectral equation in approximate harmonic oscillator basis for pseudoscalar quarkonia (for states 1S, ..., 4S) and vector quarkonia (for states 1S, ..., 3D) as Eq. (41) in [30], for the perturbative calculation of the short ranged one-gluon-exchange interaction shown later in this section.
The plots of unperturbed wave functions for pseudoscalar and vector quarkonia are given in [30]. The unperturbed mass spectra is expressed as (see [30]), The mass spectra of vector charmonium and bottomonium states using the above spectral equation was found to have degenerate S, and D states [30]. However with the incorporation of the OGE (Coulomb) term the mass spectral equation can be written as: Now, treating the coulomb term as a perturbation to the unperturbed mass spectral equation, Eq. (22), and treating the wave functions φ P,V ( q) in Eq. (41) of [30], as the unperturbed wave functions, we can write the complete mass spectra of ground (1S) and excited states for equal mass heavy pseudoscalar (0 −+ ) and vector (1 −− ) mesons respectively using the first order degenerate perturbation theory as: (N +   3 2 ); N = 2n + l; n = 0, 1, 2...., (25) where < V P,V coul > (has again been weighted by a factor of γ P,V = (m 1 +m 2 ) 2 , as in scalar (0 ++ ) quarkonia in previous section) is the matrix element of V P,V coul between unperturbed states in Eq.(41) of [30] of given quantum numbers n, and l, (with n = 0, 1, 2, 3, ..., and l = 0 for pseudoscalar mesons, while with l = 0, 2 for vector mesons). It is to be noted that, V P,V coul connects only the equal parity states with the same values of quantum numbers n. The only non-vanishing matrix elements of the perturbation between states with the given quantum numbers n and l are listed below: The non-zero values of < V coul > given above not only lead to the lifting up of the degeneracy between the S and D levels with the same principal quantum number N in vector quarkonia, but also leads to bringing the masses of different states of vector and pseudoscalar quarkonia closer to data, as can be seen from the mass spectral results for pseodoscalar (0 −+ ), and vector (1 −− ) quarkonia, which are compared with the experimental data [6], and other models for each state (where ever available), as given in Tables   3 and 4 Table 4: Masses of ground, radially and orbitally excited states of heavy vector quarkonium, J/ψ in BSE-CIA along with their masses in other models and experimental data (all units are in GeV).

Mass spectral equation for axial vector 1 ++ quarkonia
The general form for the relativistic Salpeter wave function of 3 P 1 state with J P C = 1 ++ can be expressed as in [17,19]. In the center of mass frame, we can then write the general decomposition of the instantaneous BS wave function for axial vector mesons (J pc = 1 ++ ) of dimensionality M as: With use of our power counting rule [23,24], it can be checked that the Dirac structures associated with amplitudes g 1 and g 2 are O(M 1 ), and are leading, and would contribute maximum to any axial vector meson calculation. Following a similar procedure as in the case of scalar mesons, we can write the Salpeter wave function in terms of only two leading Dirac amplitudes g 1 , and g 2 . Plugging this wave function together with the projection operators into the first two Salpeter equations in Eq.(16) of [30], and taking trace on both sides and following the steps as for the scalar meson case, we get the coupled integral equations in the amplitudes g 1 , and g 2 : To decouple these equations, we follow a similar procedure as in the scalar meson case, and get two identical decoupled equations, one entirely in g 1 (q), and the other that is entirely in, g 2 (q). In the limit, ω ≈ m on RHS, and due to the fact that for axial quarkonia again, ω 4 qq ω 2 qq these equations can be expressed as: where, it can be checked that both g 1 and g 2 satisfy the same equation and hence we can approximately where Θ A = γ µ Ψ( q)γ µ . Now, with the use of leading Dirac structures, we can again to a good approximation express Θ A = 2. This is due to the fact that γ µ ψ( q)γ µ ≈ γ µ γ 5 γ ν (iM g 1 )γ µ = 2γ 5 γ ν (iM g 1 ) ≈ 2ψ( q).
This mass spectral equation has the same form as the mass spectral equation for scalar quarkonia, Eq. Eqs. (18) for 0 ++ case, and thus the unperturbed wave functions φ A ( q) for 1 ++ would then have the same algebraic form as φ s ( q) in Eqs. (17), but with β s → β A : The plots of wave functions for axial vector quarkonia are given in Fig.2 below Now, the perturbative inclusion of the coulomb term will reduce mass spectrum for axial vector meson in the same form as Eqs. (19) - (21) for the scalar case, except that the inverse range parameter β s has got to be replaced by β A . This complete mass spectral equation for ground and excited states of 1 ++ is: The solutions of the above spectral equation is: 2 ); N = 2n + l; n = 0, 1, 2...., (33) with l = 1, where < V A coul > is the expectation value of V A coul between the unperturbed states of given quantum numbers n (with l = 1) for axial vector mesons, and has been weighted by a factor of (m 1 +m 2 ) 2 to have the coulomb term dimensionally consistent with the harmonic term. Its expectation values for the 1P, 2P , and 3P states are: From the mass spectral equation, one can see that, the mass spectra again depends not only on the principal quantum number N , but also the orbital quantum number l. We are now in a position to calculate the numerical values for mass spectra of heavy scalar quarkonia with the input parameters of our model as mentioned above. The results of mass spectral predictions of heavy equal mass scalar mesons for both ground and excited states with the above set of parameters is given in table 2. The results on masses of ground (1P ) and excited (2P and 3P ) states of quarkonia, χ c1 are given in Table 5.  Table 5: Mass spectrum of ground and excited states of χ c1 with quantum numbers J P C = 1 ++ in GeV.
units with the above set of parameters.

Two-photon decays of scalar quarkonium
We now study the two-photon decay width of a scalar quarkonium (0 ++ ), which proceeds through the quark-triangle diagrams shown in Fig.3. Let P be the total momentum of the scalar quarkonia, and k 1.2 be the momenta of the two emitted photons with polarizations ε 1,2 respectively. Then we can write P = k 1 + k 2 , and let 2Q = k 1 − k 2 . The invariant amplitude for this process can be written as: where e Q = + 2 3 e for cc, The 3D structure of Dirac wave function, Ψ s ( q) is given in Eqs. (9), and (17). The propagators for the third quark in the two diagrams is expressed as, S F (q ∓ Q) = −i( / q∓ / Q)+m (q∓Q) 2 +m 2 . Now for heavy hadrons, where the system can be regarded as non-relativistic, it is a good approximation to take the internal momentum q << M , and hence q 2 << Q 2 , where it can be seen that Q 2 = M 2 4 . Using the propagator expressions given above, and evaluating trace over the gamma matrices, we can write the invariant amplitude given above as: where F S is the decay constant for scalar quarkonium, χ c0 . The decay width for the process can then be expressed as,  Table 6: Two-photon decay widths of ground and first excited states of χ c0 (1P ) and χ c0 (2P ) in eV with the input set of parameters given after Eq.(7)

Discussions
We have employed a 3D reduction of BSE (with a 4 × 4 representation for two-body (qq) BS amplitude) under Covariant Instantaneous Ansatz (CIA) for deriving the algebraic forms of the mass spectral equations for scalar, pseudoscalar, vector and axial vector quarkonia using the full BS kernel comprising of the onegluon-exchange and the confining part in an approximate harmonic oscillator basis that led to analytic solutions (both eigen functions and eigen values). We thus obtain the mass spectra of cc quarkonia for ground and excited states of 0 ++ , 0 −+ , 1 −− ,and 1 ++ states. The mass spectral results for all these states which are compared with the experimental data [6], and other models for each state (where ever available), are given in Tables 2-5 respectively. The masses, and the algebraic forms of eigen functions for each quarkonium state so obtained will be used for calculating their various transitions.
Here we wish to mention that in principle, in QQ quarkonia, the constituents are close enough to each other to warrant a more accurate treatment of the coulomb term. Though for bb systems, the coulomb term will be extremely dominant in comparison to confining term, however, it may not be so unreasonable to treat coulomb term perturbatively for cc systems. Here, we wish to point out that if we are getting reasonable results for orbital cc excitations, it is mainly due to centrifugal effects [34] which ensure that c − c separation is large enough to feel the effect of confining term more strongly than the coulomb term.
Further, the present approach is on lines of some earlier works [34,35,36] where the OGE(coulomb) term is treated perturbatively for charmed mesons and baryons. Similarly in a recent work [37], the 3D harmonic oscillator wave functions were used as a trial wave functions for obtaining cc spectrum and their decays, where these trial wave functions did not take into account the importance of coulomb potential for heavy quarkonium systems. Similar treatment was earlier followed in [38].
We further wish to point out that this analytic approach (giving explicit dependence of spectra on principal quantum number N ) under heavy quark approximation gives a much deeper insight into the spectral problem, than the purely numerical approaches prevalent in the literature. The plots of the analytical forms of wave functions for various J P C states, obtained as solutions of their mass spectral equations are also are given in Figs.1-2 for scalar and axial vector quarkonia. The correctness of our approach can be gauged by the fact that these plots are very similar to the corresponding plots of amplitudes for these quarkonia in [20] obtained by purely numerical methods. The analytical forms of eigen functions for ground and excited states so obtained can be used to evaluate the various other processes involving scalar and axial vector quarkonia, which we intend to do next.