Double charmonium productions in electron-positron annihilation using Bethe-Salpeter approach

We calculate the double charmonium production cross section within the framework of $4\times 4$ Bethe-Salpeter Equation in the electron-positron annihilation, at center of mass energy $\sqrt{s}= 10.6$GeV, that proceeds through the exchange of a single virtual photon. In this calculation, we make use of the full Dirac structure of 4D BS wave functions of these charmonia, with the incorporation of all the Dirac covariants (both leading and sub-leading). The calculated cross-sections for the double charmonium productions for final states, ($J/\Psi, \eta_c$); ($\Psi',\eta_c$) ; ($J/\Psi, \eta_c'$); and ($\Psi',\eta_c'$) are close to experimental data, and in broad agreement with results of other theoretical models.

Among them, the study of charmonium production in e + e − annihilation is particularly interesting in testing the quarkonium production mechanisms at center of mass energies √ s = 10.6GeV . The investigation of double charmonium production is very important since these charmonium can be easily produced in experiments and hence their theoretical prediction can verify the discrepancy between different theoretical models and experimental data. As regards the dynamical framework, to investigate the double charmonium production is concerned, many approaches have been proposed to deal with the cross -section of double charmonium production [5,6,7,9,8] and the pseudoscalar and vector charmonium production process has been recently studied in a Bethe-Salpeter formalism [10,11]. However, in these studies the complete Dirac structure of P (pseudoscalar)and V (vector) quarkonia was not taken into account, and calculations were performed by taking only the leading Dirac structures, γ 5 , and iγ.ε in the BS wave functions of pseudoscalar and vector charmonia respectively.
In these calculations, we make use of the Bethe-Salpeter equation (BSE) approach [12,13,14,15,16], which is a conventional approach in dealing with relativistic bound state problems. Due to its firm base in quantum field theory, and being a dynamical equation based approach, it provides a realistic description for analyzing hadrons as composite objects, and can be applied to study not only the low energy hadronic processes, but also the high energy production processes involving quarkonia as well.
In our recent works [17,18,19,20], we employed BSE under Covariant Instantaneous Ansatz, which is a Lorentz-invariant generalization of Instantaneous Approximation, to investigate the mass spectra and the transition amplitudes for various processes involving charmonium and bottomonium. The BSE framework using phenomenological potentials can give consistent theoretical predictions as more and more data are being accumulated. In our studies on 4 × 4 BSE, in all processes except in [19], the quark -anti quark loop involved a single hadron-quark vertex, which was simple to handle. However for the transitions such as V → P + γ (where V vector and P pseudoscalar quarkonium), which we have studied in [19] and for the process of double charmonium production, e + e − → V +P , we will study in present work, the process requires calculation involving two hadron-quark vertices, due to which the calculation becomes more and more difficult to handle.
However in the present work, we demonstrate an explicit mathematical procedure for handling this problem using the formulation of 4 × 4 Bethe -Salpeter Equation under Covariant Instantaneous Ansatz. We will use this framework for the calculation of cross section for the production of double charmonia, (J/ψ, η c ), (J/ψ, η c ), (ψ , η c ), and (ψ , η c ) in electron-positron annihilation that proceeds through a single virtual photon at center of mass energies, √ s = 10.6 GeV., where such problems do not enter in our previous papers on 4 × 4 BSE. [17,18,19,20,21,22].
This paper is organized as follows. In section 2, We introduce the detailed formulation of the transition amplitudes for the process of double charmonium production and the numerical results of cross -sections. Finally, we give the discussions and conclusions in section 3.

Formulation of double charmonium production amplitude
We start from the lowest order Feynman diagrams for the process, e + e − → V + P , as given in  The relativistic amplitude M 1 f i for double charmonium production, corresponding to Fig.1, is given by the one-loop momentum integral as: where, s is the Mendelstam variable defined as, s = −(p 1 + p 2 ) 2 , α em = e 2 4π is called the electromagnetic coupling constant and α s = g 2 s 4π is the strong coupling strength, Ψ(P a , q a ) and Ψ(P b , q b ) are the conjugations of the BS wave function of vector and pseudoscalar charmonium respectively.
From the figure, we can relate the momenta of the quark and anti-quark respectively as: and the momenta for the propagators of gluon and quark are given respectively by: As the quark and gluon propagators depend upon the internal hadron momenta q a and q b , the calculation of amplitude is going to involve integrations over these internal momenta and will be quite complex. Hence following Ref. [17,18,19], we simplify the calculation, by reducing the 4dimensional expression of BS amplitude, M 1 f i into 3-dimensional expression of BS amplitude, M 1 f i and employing the heavy quark approximation on the propagators, the momenta k and q 1 can be written as, k ≈ 1 2 (P a +P b ) and q 1 ≈ P a + 1 2 P b , which leads to, k 2 ≈ s 4 and q 2 1 ≈ s 2 +m 2 c [10,11]. Then, with the above approximation and applying the definition of 3D BS wave function in [17,18,19,20] f i can be written in the instantaneous Bethe-Salpeter amplitude form as: where ψ(q a ) and ψ(q b ) are the relativistic BS wave function of pseudoscalar and vector charmonium respectively. We now give details of calculation of double charmonium production for the process, e + e − → V + P in the next section.
• For the production process, e + + e − → V + P The relativistic BS wave function of pseudoscalar and vector charmonium are taken from our recent papers [17,18,19] respectively as: for P b and M b , is the momentum and mass of the pseudoscalar charmonium respectively and N P is the BS normalization of the pseudoscalar charmonium, which is given in a simple form as: and were, ε is the polarization vector of the vector quarkonia, P a and M a , is the momentum and mass of the vector quarkonia respectively and N V is the BS normalization of the vector quarkonia, which is given in a simple form as: where φ P (q b ) and φ V (q a ) are the radial wave functions of P and V , which are solutions of the 3D BSE for pseudoscalar and vector quarkonia respectively (see Ref. [18]). The adjoint BS wave functions for pseudoscalar and vector charmonia are obtained from ψ P,V (q b,a ) = γ 0 (ψ P,V (q b,a )) + γ 0 . With the substitution of adjoint BS wave functions of P and V charmoni into the 3D amplitude Eq. (4), M 1 f i becomes: were Applying the trace theorem and evaluating trace over the gamma matrices, one can obtain the expression: After some mathematical steps, we can get the following expression: Thus we can express the amplitude M 1 f i for the vector and pseudoscalar charmonium production as: The full amplitude for the process e + + e − → V + P , can be obtained by summing over the amplitudes of all the possibilities in Fig. 1. Then, the unpolarized total cross-section is obtained by summing over various V + P , spin-states and averaging over those of the initial state e + e − , which is given as: The total amplitude, |M total | 2 is given as: After integration over the phase space, the total cross-section for vector and pseudoscalar charmonium production is: The algebraic expressions of the wave functions of pseudoscalar and vector charmonium, a) ), for ground (1S) and first excited (2S) states respectively, that are obtained as analytic solutions of the corresponding mass spectral equations of these quarkonia in an approximate harmonic oscillator basis are [18]: where the inverse range parameter β P for pseudoscalar charmonium is β P = (4

Numerical results
We had calculated recently the mass spectrum and various decays of ground and excited states of pseudoscalar and vector quarkonia in [18,19,20]. The same input parameters are employed in this calculation as in our recent works, given in table 1.  Table 1: Input parameters for this study Using these input parameters listed in Table 1, we calculate the cross sections of pseudoscalar and vector charmonium production in our framework are listed in Table 2.
We wish to mention here that the bbbb as well as bbcc production has not been observed so far, but cross sections have been predicted for e − e + → Υ + η b at center of mass energy, √ s = 25 − 30 GeV. We thus wished to check our calculations for double bottomonium (Υη b ) production in electron-positron annihilation for sake of completeness. The same can be studied with our framework used in this work with little modifications. Taking the mass of the b−quark that is fixed from our recent work on spectroscopy of bb states as m b = 5.07GeV [18], and other input  parameters the same, the cross section for production of ground states of η b , and Υ for energy range ( √ s = 28 − 30)GeV is given in Table 3.

Discussions and conclusion
Our main aim in this paper was to study the cross section for double charmonium production in electron-positron collisions at center of mass energies, √ s = 10.6GeV., having successfully studied the mass spectrum and a range of low energy processes using an analytic treatment of 4 × 4 Bethe-Salpeter Equation (BSE) [18,19,20,21,22]. This is due to the fact any quark model model should successfully describe a range of processes-not only the mass spectra, and low energy hadronic decay constants/decay widths, but also the high energy production processes involving these hadrons, and all within a common dynamical framework, and with a single set of input parameters, that are calibrated to the mass spectrum.
Further, the exclusive production of double heavy charmonia in e + e − annihilation has been a challenge to understand in heavy-quark physics, and has received considerable attention in recent years, due to the fact that there is a significant discrepancy in the data of Babar [2] and Belle [1], and the calculations performed in NRQCD [5,6], of this process at √ s = 10.6GeV. Using Leading order (LO) QCD diagrams alone, lead to cross sections that are an order of magnitude smaller than data [2,1,4]. These discrepancies were then resolved by taking into account the Next-to-Leading Order (NLO) QCD corrections combined with relativistic corrections [25,26]), though the NLO contributions were larger than the LO contributions [10].
This process has also been studied in Relativistic Quark Model [7,8], Light Cone formalism [9], and Bethe-Salpeter Equation [10,11]. Relativistic corrections to cross section in Light-Cone formalism (in [9]) were considered to eliminate the discrepancy between theory and experiment.
Further in an attempt to explain a part of this discrepancy, [27] even suggested that processes proceeding through two virtual photons may be important. In Relativistic quark model [7], which incorporated relativistic treatment of internal motion of quarks, and bound states, improvements in the results were obtained.
We wish to mention that in our framework of 4 × 4 Bethe-Salpeter Equation (BSE) using the Covariant Instantaneous Ansatz, where we treat the internal motion of quarks and the bound states in a relativistically covariant manner, we obtained results on cross sections (in Table 2) for production of opposite charge parity cc states, such as, J/Ψ, η c ; Ψ , η c ; J/Ψ, η c ; and Ψ , η c , that are close to central values of data [2,1,4] using leading order (LO) QCD diagrams alone. This validates the fact that relativistic quark models such as BSE are strong candidates for treating not only the low energy processes, but also the high energy production processes involving double heavy quarkonia, due to their consistent relativistic treatment of internal motion of quarks in the hadrons, where our BS wave functions that take into account all the Dirac structures in pseudoscalar and vector mesons in a mathematically consistent manner, play a vital role in the dynamics of the process. We also give our predictions on cross section for double bb production at energies [28][29][30] GeV. in Table 3 for future experiments at colliders.
The main objective of this study was to test the validation of our approach, which provides a much deeper insight than the purely numerical calculations in 4 × 4 BSE, that are very common in the literature. We wish to mention that, we have not encountered any work in 4 × 4 representation of BSE, that treats this problem analytically. On the contrary all the other 4 × 4 BSE approaches adopt a purely numerical approach of solving the BSE. We are also not aware of any other BSE framework, involving 4 × 4 BS amplitude, and with all the Dirac structures incorporated in the 4D hadronic BS wave functions (in fact many works used only the leading Dirac structures, for instance see Ref. [10,11]) for calculations of this production process. We treat this problem analytically by making use of the algebraic forms of wave functions, φ (P,V ) derived analytically from mass spectral equations [18], for calculation of cross section of this double charmonium production process.
This calculation involving production of double charmonia in electron-positron annihilation can be easily be extended to studies on other processes (involving the exchange of a single virtual photon) observed at B-factories such as, e − e + − > Ψ(2S)γ, e − e + − > J/Ψχ c0 , and e − e + − > Ψ(2S)χ c0 . We further wish to extend this study to processes involving two virtual photons, such as the production of double charmonia in final states with C = +1 such as, e − e + → J/ΨJ/Ψ), recently observed at BABAR.
These processes that we intend to study next, are quite involved, and the dynamical equation based approaches such as BSE is a promising approach. And since these processes involve quarktriangle diagrams, we make use of the techniques we used for handling such diagrams recently done in [19,28,29]. Such analytic approaches not only lead to better insight into the mass spectra, and low energy decay processes involving charmonia, but also their high energy production processes such as the one studied here.
Note added: At proof correction stage, we come to know of about the recent experimental observation by CMS collaboration [30] of two excited B + c (2S), and B * + c (2S) states in pp collisions at √ s = 13 TeV. at Large Hadron Collider (LHC). The theoretical study of this process will be a challenge to hadronic physics.