Charmonium spectrum with a Dirac potential model in the momentum space

TL;DR

This paper models the charmonium spectrum using a relativistic potential approach in momentum space, fitting known states and predicting higher resonances, including the X(3872), with quantified uncertainties.

Contribution

It introduces a relativistic Dirac potential model with screening factors in momentum space, avoiding nonrelativistic approximations and providing uncertainty estimates for resonance predictions.

Findings

Successfully reproduces the charmonium spectrum structure

Predicts higher-lying resonances including X(3872)

Provides uncertainties and parameter correlations for the model

Abstract

We study the charmonium spectrum using a complete one gluon exchange approach based on a phenomenological relativistic $q\bar{q}$ potential model with Dirac spinors in momentum space. We use phenomenological screening factors to include vacuum quantum effects. Our formulation does not rely on nonrelativistic approximations. We fit the lowest-lying charmonia (below the $D\bar{D}$ threshold) and predict the higher-lying resonances of the spectrum. In general, we reproduce the overall structure of the charmonium spectrum and, in particular, we can reasonably describe the $X(3872)$ resonance mass as (mostly) a $c\bar{c}$ state. The numerical values of the free parameters of the model are determined taking into account also the experimental uncertainties of the resonance energies. In this way, we are able to obtain the uncertainties of the theoretical resonance masses and the correlation among the free parameters of the model.

Tables

Table 1. Table 1: Fitted charmonia for potentials I and II compared to the experimental masses from PDG. n 𝑛 n stands for the principal quantum number, L 𝐿 L for the orbital angular momentum, J 𝐽 J for the total angular momentum and S 𝑆 S for the spin.

Name	$n^{2S+1}{L}_J$	Masses (MeV)
		Potential I	Potential II	Experiment
${\eta }_c$	$1^1{S}_0$	${2990}_{-12}^{+6}$	${2981}_{-12}^{+8}$	$2983.6\pm 0.6$
$J/\psi$	$1^3{S}_1$	${3089.3}_{-6.5}^{+5.7}$	${3101}_{-10}^{+13}$	$3096.916\pm 0.011$
${\chi }_{c0}$	$1^3{P}_0$	${3417.0}_{-5.9}^{+8.4}$	${3416}_{-13}^{+11}$	$3414.75\pm 0.31$
${\chi }_{c1}$	$1^3{P}_1$	${3500.7}_{-3.4}^{+4.4}$	${3506.3}_{-7.1}^{+7.7}$	$3510.66\pm 0.07$
$h_c$	$1^1{P}_1$	${3514.9}_{-4.6}^{+5.9}$	${3521.0}_{-6.7}^{+6.5}$	$3525.38\pm 0.11$
${\chi }_{c2}$	$1^3{P}_2$	${3579}_{-12}^{+5}$	${3563}_{-14}^{+8}$	$3556.20\pm 0.09$
${\eta }_c^{\prime }$	$2^1{S}_0$	${3646}_{-12}^{+7}$	${3647}_{-12}^{+10}$	$3639.2\pm 1.2$
${\psi }^{\prime }$	$2^3{S}_1$	${3674.9}_{-5.7}^{+5.4}$	${3679.2}_{-9.4}^{+9.2}$	${3686.109}_{-0.014}^{+0.012}$

Table 2. Table 2: Fit parameters for the two potentials considered in this work. Error bars are reported at the 1 σ 𝜎 \sigma (68%) CL and take into account all the correlations among parameters (see Section  III.2 ).

Parameter	Potential I	Potential II
$m$ (MeV)	${1455}_{-26}^{+30}$	${1364}_{-63}^{+46}$
$M_0$ (MeV)	${108}_{-25}^{+21}$	$270\pm 45$
${\alpha }_{st}$	${0.63}_{-0.06}^{+0.08}$	${0.544}_{-0.060}^{+0.027}$
${\beta }_v$ (GeV2)	${0.0134}_{-0.0014}^{+0.0009}$	${0.004}_{-0.002}^{+0.001}$
$k_s$	${59}_{-9}^{+11}$	${96}_{-10}^{+14}$
$p_s$ (GeV)	${0.698}_{-0.067}^{+0.065}$	${0.91}_{-0.06}^{+0.11}$
$A$ (GeV-2)	$0.015\pm 0.004$	$-0.006\pm 0.002$
${\beta }_s$ (GeV2)	$0$ (fixed)	${0.013}_{-0.003}^{+0.004}$

Table 3. Table 3: Correlation matrix for the parameters of potential I.

	$m$	$M_0$	${\alpha }_{st}$	${\beta }_v$	$k_s$	$p_s$	$A$
$m$	$1$
$M_0$	$-0.78$	$1$
${\alpha }_{st}$	$0.86$	$-0.40$	$1$
${\beta }_v$	$-0.88$	$0.39$	$-0.93$	$1$
$k_s$	$0.11$	$-0.02$	$0.13$	$-0.23$	$1$
$p_s$	$0.33$	$-0.34$	$0.05$	$-0.16$	$-0.62$	$1$
$A$	$-0.76$	$0.59$	$-0.53$	$0.62$	$0.03$	$-0.70$	$1$

Table 4. Table 4: Correlation matrix for the parameters of potential II.

	$m$	$M_0$	${\alpha }_{st}$	${\beta }_v$	$k_s$	$p_s$	$A$	${\beta }_s$
$m$	$1$
$M_0$	$-0.61$	$1$
${\alpha }_{st}$	$0.40$	$0.31$	$1$
${\beta }_v$	$-0.12$	$0.04$	$0.06$	$1$
$k_s$	$0.31$	$-0.10$	$0.14$	$-0.23$	$1$
$p_s$	$0.02$	$-0.20$	$-0.42$	$-0.17$	$-0.35$	$1$
$A$	$0.22$	$-0.18$	$0.09$	$-0.21$	$0.24$	$0.58$	$1$
${\beta }_s$	$-0.29$	$-0.20$	$-0.52$	$-0.66$	$0.02$	$0.33$	$0.29$	$1$

Table 5. Table 5: Predicted charmonia for potentials I and II compared to the experimental masses and uncertainties from PDG. Notation as in Table  1 . For the X ​ ( 4160 ) 𝑋 4160 X(4160) experimental mass the first error is statistical and the second systematic  Abe:2007sya .

Name	$n^{2S+1}{L}_J$	Mass (MeV)
		Potential I	Potential II	Experiment
—	$1^1{D}_2$	${3788}_{-12}^{+20}$	${3789}_{-13}^{+9}$	—
$\psi \left(3770\right)$	$1^3{D}_1$	${3723}_{-8}^{+41}$	${3749}_{-10}^{+7}$	$3773.13\pm 0.35$
$\psi \left(3823\right)$	$1^3{D}_2$	${3782}_{-11}^{+30}$	${3788}_{-12}^{+9}$	$3822.2\pm 1.2$
—	$1^3{D}_3$	${3849}_{-17}^{+13}$	${3817}_{-20}^{+14}$	—
$X\left(3915\right)$ [$Y\left(3915\right)$]	$2^3{P}_0$	${3867.9}_{-8.6}^{+5.5}$	${3884}_{-13}^{+14}$	$3918.4\pm 1.9$
—	$2^1{P}_1$	${3908}_{-13}^{+11}$	${3923}_{-13}^{+11}$	—
$X\left(3872\right)$	$2^3{P}_1$	${3903.5}_{-9.4}^{+9.3}$	${3918}_{-12}^{+13}$	$3871.69\pm 0.17$
${\chi }_{c2}\left(2P\right)$	$2^3{P}_2$	${3941}_{-16}^{+8}$	${3945}_{-14}^{+17}$	$3927.2\pm 2.6$
$\psi \left(4160\right)$	$2^3{D}_1$	${4059}_{-20}^{+13}$	${4094}_{-20}^{+19}$	$4191\pm 5$
$X\left(3940\right)$	$3^1{S}_0$	${4030}_{-71}^{+9}$	${4043}_{-19}^{+18}$	$3942\pm 6$
$\psi \left(4040\right)$	$3^3{S}_1$	${4046}_{-13}^{+9}$	${4060}_{-11}^{+17}$	$4039\pm 1$
—	$3^3{P}_0$	${4182}_{-47}^{+17}$	${4219}_{-21}^{+26}$	—
$X$(4140) [$Y$(4140)]	$3^3{P}_1$	${4200}_{-31}^{+24}$	${4238}_{-22}^{+29}$	$4146.9\pm 3.1$
$X\left(4160\right)$	$4^1{S}_0$	${4300}_{-100}^{+30}$	${4349}_{-26}^{+31}$	${4156}_{-20}^{+25}\pm 15$

Table 6. Table 6: Predicted states that can be interpreted differently depending on the potential model. Notation as in Table  1 .

Name	Potential	$n^{2S+1}{L}_J$	Mass (MeV)
			Theory	Experiment
$X\left(4260\right)$	I	$4^3{S}_1$	${4311}_{-44}^{+33}$	$4251\pm 9$
$\left[Y\left(4260\right)\right]$	II	—–
$X\left(4360\right)$	I	$3^3{D}_1$	${4312}_{-53}^{+62}$	$4346\pm 6$
$\left[Y\left(4360\right)\right]$	II	$4^3{S}_1$	${4360}_{-28}^{+34}$
$\psi \left(4415\right)$	I	$4^3{D}_1$	${4489}_{-78}^{+65}$	$4421\pm 4$
	II	$3^3{D}_1$	${4378}_{-29}^{+36}$
$X\left(4660\right)$	I	——		$4643\pm 9$
$\left[Y\left(4660\right)\right]$	II	$4^3{D}_1$	${4620}_{-39}^{+54}$

Equations

$$\begin{array}{rl}{\vec{p}}_{1a}& =-{\vec{p}}_{2a}={\vec{p}}_a,\\ {\vec{p}}_{1b}& =-{\vec{p}}_{2b}={\vec{p}}_b,\end{array}$$

$$\vec{q}={\vec{p}}_b-{\vec{p}}_a,$$

$$x=\cos \theta ={\hat{p}}_a\cdot {\hat{p}}_b,$$

$${\vec{q}}^2={\vec{p}}_a^2+{\vec{p}}_b^2-2p_a{p}_{b}x,$$

$$\begin{array}{r}E({\vec{p}}_{ia})=E({\vec{p}}_a)=\sqrt{{\vec{p}}_a^2+m^2},\\ E({\vec{p}}_{ib})=E({\vec{p}}_b)=\sqrt{{\vec{p}}_b^2+m^2},\end{array}$$

$$\Delta E=E({\vec{p}}_b)-E({\vec{p}}_a),$$

$$Q^2={\vec{q}}^2-(\Delta E{)}^2.$$

$$u({\vec{p}}_i,{\vec{\sigma }}_i)=\sqrt{\frac{E({\vec{p}}_i)+m}{2E({\vec{p}}_i)}}\left(\begin{array}{c}1\\ \frac{{\vec{p}}_i\cdot {\vec{\sigma }}_i}{E({\vec{p}}_i)+m}\end{array}\right)\chi ,$$

$${\bar{\mathcal{H}}}_{\text{int}}({\vec{p}}_b,{\vec{p}}_a)=⟨\overrightarrow{p_b}|H^{(v)}|\overrightarrow{p_a}⟩+⟨\overrightarrow{p_b}|H^{(s)}|\overrightarrow{p_a}⟩,$$

$${\bar{\mathcal{H}}}_{\text{int}}({\vec{p}}_b,{\vec{p}}_a)\to {\mathcal{H}}_{\text{int}}({\vec{p}}_b,{\vec{p}}_a)=F_s(p_b){\bar{\mathcal{H}}}_{\text{int}}({\vec{p}}_b,{\vec{p}}_a)F_s(p_a),$$

$$F_s(p)=\frac{1+k_s}{k_s+\exp \left(p^2/p_s^2\right)},$$

$$⟨{\vec{p}}_b|H^{(v)}|{\vec{p}}_a⟩=\alpha \left(J_1^0{J}_2^0{D}_{00}+J_1^{\alpha }{J}_2^{\beta }{D}_{\alpha \beta }\right).$$

$$D_{00}$$

$$D_{\alpha \beta }$$

$$(\Delta E)J_{1,2}^0=\pm \vec{q}\cdot {\vec{J}}_{1,2},$$

$$\begin{array}{rl}⟨\overrightarrow{p_b}|H^{(v)}|\overrightarrow{p_a}⟩=& V^{(v)}(\vec{q})[J_1^0{J}_2^0\left(1-\frac{(\Delta E{)}^2}{Q^2}\right)\\ & -{\vec{J}}_1\cdot {\vec{J}}_2\left(1+\frac{(\Delta E{)}^2}{Q^2}\right)].\end{array}$$

$$J_i^{\mu }=J_i^{\mu }({\vec{\sigma }}_i;{\vec{p}}_b,{\vec{p}}_a)=\bar{u}({\vec{p}}_{ib},{\vec{\sigma }}_i){\gamma }_i^{\mu }u({\vec{p}}_{ia},{\vec{\sigma }}_i).$$

$$⟨{\vec{p}}_b|H^{(s)}|{\vec{p}}_a⟩=V^{(s)}(\vec{q})I_1{I}_2,$$

$$I_i=I_i({\vec{\sigma }}_i;{\vec{p}}_b,{\vec{p}}_a)=\bar{u}({\vec{p}}_{ib},{\vec{\sigma }}_i)u({\vec{p}}_{ia},{\vec{\sigma }}_i).$$

$$V^{(v)}(\vec{q})=-\frac{4}{3}\frac{{\alpha }_{st}}{{\vec{q}}^2}+{\beta }_v\frac{3b^2-{\vec{q}}^2}{({\vec{q}}^2+b^2{)}^3},$$

$$V^{(s)}(\vec{q})=A+{\beta }_s\frac{3b^2-{\vec{q}}^2}{({\vec{q}}^2+b^2{)}^3}.$$

$$V^{(v)}(r)=-\frac{4}{3}\frac{{\alpha }_{st}}{r}+{\beta }_v{\pi }^{2}r\exp (-rb),$$

$$\{\begin{array}{l}\text{potential II}\to \text{model using Eqs. (refvtotq) with }{\beta }_s=0,\\ \text{potential II}\to \text{model using Eqs. (refvtotq) with }{\beta }_s\ne 0.\end{array}$$

$$K(\vec{p})=2\sqrt{{\vec{p}}^2+m^2}.$$

$$[K(\overrightarrow{p_b})+M_0]\Psi ({\vec{p}}_b)+\int d^3{p}_a{\mathcal{H}}_{\text{int}}({\vec{p}}_b,{\vec{p}}_a)\Psi ({\vec{p}}_a)=M\Psi ({\vec{p}}_b),$$

$${\Psi }_{n,\{\nu \}}(\vec{p})=R_{n,L}(p;\bar{p}){\left[Y_L(\hat{p})\otimes {\chi }_S\right]}_{J,M_J},$$

$$R_{n,L}(p;\bar{p})={\left[\frac{2n!}{\Gamma (n+L+\frac{3}{2})}\right]}^{\frac{1}{2}}\frac{s^L}{{\bar{p}}^{\frac{3}{2}}}{\mathcal{L}}_n^{L+\frac{1}{2}}(s^2)\exp \left[\frac{-s^2}{2}\right],$$

$$M_{\{\nu \},n_b,n_a}=M_0{\delta }_{n_b,n_a}+\int d^{3}p{\Psi }_{n_b,\{\nu \}}^{†}(\vec{p})K(\vec{p}){\Psi }_{n_a,\{\nu \}}(\vec{p})+\int d^3{p}_b{d}^3{p}_a{\Psi }_{n_b,\{\nu \}}^{†}(\overrightarrow{p_b}){\mathcal{H}}_{\text{int}}({\vec{p}}_b,{\vec{p}}_a){\Psi }_{n_a,\{\nu \}}(\overrightarrow{p_a}),$$

$$d^2=\sum _i{\left(E_i-M_i\right)}^2,$$

$${\bar{E}}^{\text{j}}=2\sqrt{m^2+{\left(p_{1/2}^{\text{j}}\right)}^2}$$

$$F(\vec{p})=\sqrt{\frac{E(\vec{p})+m}{2E(\vec{p})}},R(\vec{p})=\sqrt{\frac{E(\vec{p})-m}{2E(\vec{p})}},$$

Full text

I Introduction

In the last years the interest in charmonia has been renewed by the experimental observation of new charmoniumlike states Olsen1 ; brambilla ; Godfrey:2008nc ; Esposito:2016noz ; Swanson1 ; Esposito:2014rxa ; Lebed . Some of these states have been easily identified as charmonia, e.g. $h_{c}(1P)$ Rosner ; Rubin and $\chi_{c2}(2P)$ Uehara , as reported by the Particle Data Group (PDG) PDG2016 . However there exist other important resonances that are still under examination. These states, denoted as $X$, $Y$ and $Z$, do not match the predictions of the nonrelativistic or semirelativistic $q\bar{q}$ potential models and some of them have been described as tetraquarks, molecules, glueballs and hybrids Olsen1 ; brambilla ; Esposito:2016itg ; Brodsky:2014xia ; Godfrey:2008nc ; Pilloni:2015doa ; Esposito:2016noz ; Lebed ; Guo:2008yz ; Guo-Galata ; Ber-Brambi ; Faccini:2013lda . Another possible explanation suggests that some of these states may be the result of kinematic effects Szczepaniak:2015eza ; Guo:2014iya that appear in the nearby of the thresholds of the open charm channels. The interpretation of these resonances depends, to a certain extent, on the theoretical approach that is followed to perform the amplitude analysis of the data Pilloni:2016obd ; Kang:2016jxw and to study charmonium spectroscopy.

Quantum chromodynamics (QCD) must be considered, in any case, the underlying field theory for the study of these systems. However, due to the difficulty of finding analytically (approximated) solutions for QCD, alternative approaches have been followed. For example, numerical lattice QCD (LQCD) Liu:2012ze ; Christopher ; Lang ; Gunnar , potential quark models Isgur ; Capstick ; Faustov ; DeSanctis2 ; Cao2012 ; DeSanctis3 ; klink ; poly ; Bram-Gro ; Kawanai1 ; Laschka ; PhysRevD.85.091503 ; Kawanai2 , unquenched quark models (UQM) Santo ; Santo2 ; Ferre-Galata ; FeSa1 ; FeSa2 ; FeSa3 , and Bethe-Salpeter approaches Hilger:2015hka ; Fischer ; Popovici .

LQCD represents a very promising but, at the same time, highly difficult tool. LQCD computations reproduce the structure of the charmonium spectrum and predict many undiscovered states Liu:2012ze ; Christopher ; Gunnar . In particular, LQCD is computationally taxing, especially when open channels and the resonant nature of the states are considered, and still needs to be improved in order to reproduce accurately the hadron spectrum Luscher ; Rummukainen ; Briceno1 ; Briceno2 ; Dudek:2012xn ; mesonsector ; Wilson ; Dudek . Exploratory LQCD simulations of the $D\bar{D}$ scattering in S and P waves have been recently computed Lang but they still need to be improved to achieve the precision of the lowest-energy meson sector Briceno:2016mjc . However, the construction of potential quark models whose interaction has been derived from LQCD Kawanai1 ; Laschka ; PhysRevD.85.091503 ; Kawanai2 , building up models that are not purely phenomenological Bali:2000gf , has proved to be a very useful tool for the study of these physical systems. Historically, the Cornell potential Eichen represents the first interaction that was introduced, with partial success, to study the charmonium spectrum. This model consists of a Coulomb-like potential and a linear term inserted in a standard nonrelativistic Schrödinger equation. The Coulomblike term takes into account reduction of the one gluon exchange (OGE) in the short range region. The linear term introduces quark confinement in the model, in approximate agreement with the spectroscopy of the first excited states. LQCD studies have confirmed the presence of the linear term in the long range effective interaction PhysRevD.85.091503 . However, in this model, the nonrelativistic approximation cannot be really justified and the fine structure of the spectrum, related to the spin-spin, spin-orbit and tensor interactions, has to be introduced by hand in a purely phenomenological way. Furthermore, higher order effects of the field theory reduction are also absent.

Throughout the literature several relativized quark models have been introduced in order to incorporate some relativistic effects Isgur ; Capstick ; Faustov ; DeSanctis2 ; Cao2012 ; DeSanctis3 ; klink ; poly ; Bram-Gro . A possible approach to include some of the relativistic effects in a consistent way is given by the generalization of the Fermi-Breit expansion (originally studied for the electromagnetic interaction) followed in DeSanctis2 ; Cao2012 ; DeSanctis3 to study the resonant states of strongly interacting quarks. We point out that this approach is not fully relativistic because the interaction terms are expanded in powers of $p/m$, where $p$ and $m$ are the quark momentum and mass, respectively. These models also assume that the mean values of $p/m$ are small and, consequently, only the first terms of the expansion are retained. The obtained effective interaction, by means of a Fourier transform, can be expressed in the coordinate space in terms of local operators. As a consequence, the eigenvalue equation for the effective Hamiltonian of these models can be solved with the standard techniques. For the specific case of relativized quark models for hadronic systems, the Bakamjian-Thomas construction may be used in the context of relativistic Hamiltonian dynamics klink ; poly . This approach allows us to perform exact Lorentz boosts on the wave functions of the model, but no connection between the interaction Hamiltonian and the Dirac structure of the underlying theory can be easily established. In other words, one cannot express the interaction of the model in terms of the Dirac covariants that are the building blocks of any relativistic theory, i.e. the scalar, vector, pseudoscalar, tensor and pseudotensor terms. Further details can be found in klink ; poly .

Another approach to improve the dynamical description of the quarkonium spectrum consists in employing the UQM formalism. In the UQM, starting from a relativized quark model, the coupling of the valence quarks to the continuum is taken into account through a self energy term Santo ; Santo2 ; Ferre-Galata . This coupling causes a shift in the energy, generally leading to a reduction in the masses of the quarkonia FeSa1 ; FeSa2 ; FeSa3 in the high-lying states. This effect can also be incorporated, in a more phenomenological way, through the screening of the interaction potential Li ; Gonzalez which produces a similar global result.

It is also possible to perform a four-dimensional fully relativistic computation of the heavy meson spectra solving the Bethe-Salpeter equation. This approach has been taken in the context of the rainbow ladder approximation Hilger:2015hka . This approach presents relevant theoretical advantages; however, the quality of the theoretical spectrum is not completely satisfactory and, in particular, spurious states are obtained Fischer ; Popovici as a consequence of the adopted dynamical approximations Ahlig ; Bijtebier1 .

In the present work we develop a relativistic model inspired by the OGE interaction. In our calculation we do not perform any nonrelativistic expansion. Hence, the well-known spin-orbit, spin-spin and tensor effects are automatically given by the relativistic calculation. Due to the Dirac structure of the interaction, the model is necessarily constructed in the momentum space. We include vector and scalar interaction terms as well as phenomenological screening factors Li ; Gonzalez ; Bai . We solve an integral eigenvalue Hamiltonian equation in the momentum space through a variational method and we perform the diagonalization of the Hamiltonian matrix in a selected orthonormal basis. We determine the parameters of the model by fitting the low-lying charmonia states taking into account the experimental uncertainties. These uncertainties are exactly propagated to the parameters and to the computed spectrum using the bootstrap technique NumericalRecipes ; Landay:2016cjw . This method is computationally expensive but allows a rigorous determination, from the statistical point of view, of the parameters of the model and their uncertainties, as well as the propagation of such uncertainties to the computation of the spectrum. The rigorous study of the uncertainties allows to identify the statistically significant differences between the theoretical and the experimental spectra, determining which states are in need of physics beyond the $c\bar{c}$ picture to be described. Moreover, the method allows to study the correlations among the parameters of the model, which provides insight on how physically independent are the different pieces that build the phenomenological model. As far as we know, no previous works on phenomenological models of charmonia have addressed the problem of rigorously establishing the correlations among parameters, the associated uncertainties of the parameters and their impact in the predicted spectrum.

The paper is organized as follows. In Section II we develop the general structure of our theoretical model. In particular, we show the details of the vector and scalar interactions that we use to obtain the charmonium spectrum. In Section III we discuss the methods used to solve the eigenvalue equation and to determine the free parameters of the model and their uncertainties. We also show how these uncertainties are propagated to the energy levels of the theoretical charmonium spectrum. In Section IV we provide our results compared to the available experimental data and we discuss different phenomenological interpretations. Finally, in Section V we summarize the results of the study. The details on the spin-angle matrix elements are left for Appendix A.

II Theoretical model

In this section we provide a detailed description of our theoretical model for the study of charmonium spectroscopy. The model consists in a three-dimensional reduction of the Bethe-Salpeter equation in which only the contributions of the positive energy Dirac spinors of the quarks are retained. Dirac spinors guarantee that relativity is satisfied and the Dirac structure of the interaction is highlighted. Due to the relativistic nature of the model, the well-known spin-orbit, spin-spin and tensor effects are all automatically included in our calculation. If we expand our model in powers of $p/m$ we recover the standard semirelativistic Fermi-Breit interaction Lifshitz .

II.1 Definitions

As it is customary, we work in the center of mass (CM) reference frame of the $c\bar{c}$ system and we use the system of natural units ($\hbar=c=1$). The quantities of the ket and bra states are labeled with the indices $a$ and $b$, respectively. Hence, the three-momentum of the quark (particle 1) and the three-momentum of the antiquark (particle 2) are written as

[TABLE]

where the relative three-momenta $\vec{p}_{b}$ and $\vec{p}_{a}$ are the actual kinematic variables employed in the calculations. We also introduce the three-momentum transfer

[TABLE]

and the cosine of the angle between $\vec{p}_{a}$ and $\vec{p}_{b}$

[TABLE]

so that

[TABLE]

where $p_{a}$ and $p_{b}$ are the moduli of the respective three-momenta, $\hat{p}_{a}=\vec{p}_{a}/p_{a}$ and $\hat{p}_{b}=\vec{p}_{b}/p_{b}$. The on shell energies of both particles ($i=1,2$) are

[TABLE]

where $m$ stands for the quark and antiquark masses. We introduce the energy difference (for both particles) between the bra and ket states

[TABLE]

and the invariant four-momentum transfer squared ($Q^{2}>0$)

[TABLE]

Finally, for completeness, we give the standard definition of the positive energy Dirac spinors that will be used in the calculation

[TABLE]

where $i=1,2$ is the particle index, $\vec{\sigma}_{i}$ are the Pauli matrices and $\chi$ represents a generic two component spin wave function. Using Eq. (1), we take $\vec{p}_{i}=\vec{p}_{ia}$ and $\vec{p}_{i}=\vec{p}_{ib}$, for the ket and bra states, respectively. The spinors are normalized to ${u}^{{\dagger}}(\vec{p}_{i},\vec{\sigma}_{i})u(\vec{p}_{i},\vec{\sigma}_{i})=1$.

II.2 Total interaction Hamiltonian and screening factors

The matrix elements of the total Hamiltonian in momentum space are obtained adding two terms that represent the vector and the scalar interactions

[TABLE]

where, for brevity, the spin matrices have not been explicitly written. However, the standard phenomenological potentials do not allow for a good reproduction of the spectrum, in particular above the open charm threshold Li ; Vija-Gon-Gar . A more complete dynamical treatment is necessary. Some of the vacuum quantum effects predicted by QCD can be taken into account by unquenching the quark model as in bijker ; ref ; Ferretti:2015fba . In particular, the virtual pair creation effects are incorporated to the Hamiltonian, i.e, the contribution given by the coupling of the $c\bar{c}$ states to the meson-meson continuum Santo ; Santo2 . Below threshold, this coupling to the continuum states gives rise to virtual $q\bar{q}-q\bar{q}$ components in the meson wave function that shift the charmonia Santo2 and bottomonia FeSa1 ; FeSa2 ; FeSa3 masses through a self-energy term. Part of these virtual pair production effects can be taken into account, in a more phenomenological way, introducing a screening factor in the confinement interaction in the coordinate space Li ; Gonzalez .

In this work, we introduce this screening effect phenomenologically, incorporating momentum dependent factors in the interaction Hamiltonian as follows

[TABLE]

where the factors $F_{s}(p_{a})$ and $F_{s}(p_{b})$ take into account the screening effect. In the present work we use the following stepwise function to describe such screening effects

[TABLE]

where $k_{s}$ and $p_{s}$ will be determined (as the rest of the parameters of the model) by fitting to the lowest energy levels of the charmonium spectrum. The details on the fit are provided in Section III.2.

II.3 Vector interaction

In the Coulomb gauge, the vector interaction matrix element reads

[TABLE]

If we choose $\alpha=\alpha_{em}$ (the electromagnetic fine structure constant) we obtain the standard electromagnetic interaction whose expansion up terms $p/m$ gives the Fermi-Breit interaction Lifshitz . We generalize the previous equation to the strong interaction between a quark and an antiquark DeSanctis2 . In particular, $J^{\mu}_{i}$ represents the quark ($i=1$) and antiquark ($i=2$) Dirac four-current that will be given explicitly in Eq. (17). Furthermore

[TABLE]

Given that we use positive energy Dirac spinors, we can write the continuity equation in the form

[TABLE]

where the signs $+$ and $-$ correspond to the quark and antiquark cases respectively. Using Eq. (15) we can rewrite Eq. (12) for the strong interaction as

[TABLE]

The explicit form of the effective vector potential $V^{(v)}(\vec{q})$ is provided in Eq. (20a) of Section II.5. Notice that, compared to Eq. (12), we take $\alpha=\alpha_{st}$, which represents the phenomenological strong coupling constant of the model.

Finally, the quark and antiquark four-currents are given by the standard form

[TABLE]

In this work we do not include a chromomagnetic current term, which seems to be unnecessary given the quality of the obtained spectrum as shown in Section IV.

II.4 Scalar interaction

In hadron spectroscopy the spin-obit effect is in general rather small and cannot be reproduced if we only use a vector interaction. Hence, a scalar interaction, partially responsible for quark confinement Swanson1 ; Swanson2 , is included in the interaction Hamiltonian to improve the description of the charmonium spectrum Cao2012 ; Ebert . Also in the present work we find that the scalar interaction is needed to obtain a good description of the charmonium spectrum. We write the scalar interaction in the following standard way

[TABLE]

with the explicit form of the effective scalar potential $V^{(s)}(\vec{q})$ provided in Eq. (20b) of Section II.5. The scalar vertex is given by

[TABLE]

II.5 Vector and scalar potentials

In order to include different physical effects in our effective interaction, it is necessary to use phenomenological expressions for the potentials in momentum space. There exist a wide range of models for the scalar and vector potentials DeSanctis2 ; Faustov ; Radford . In this work we use the following effective potentials:

[TABLE]

We note that only by a combination of scalar and vector confining potentials is it possible to eliminate the empirically contraindicated spin-orbit interaction in the interquark potential. The vector interaction, Eq. (20a), is given by a regularized Cornell potential in momentum space where the $-\frac{4}{3}\alpha_{st}/\vec{q}^{\>2}$ term in Eq. (20a) corresponds to the standard pure Coulombian strong potential, being $\alpha_{st}$ the effective strong coupling constant and $\frac{4}{3}$ the color factor. The second term corresponds to a confining interaction where the $\beta_{v}$ parameter is related to the vector confinement intensity. The parameter $b$ has been introduced to avoid the $1/{q}^{4}$ singularity in the momentum space interaction and, in the present model, has no direct physical meaning. The expression of the vector potential in the coordinate space reads

[TABLE]

where it is apparent that the second term represents a linear confinement interaction regularized with a decreasing exponential function. For convenience, we fix $b=10^{-2}~{}\text{GeV}$, that corresponds to $1/b\simeq 20~{}\text{fm}$. The $b$ value is chosen small enough to avoid numerical problems and large enough to approximately recover the Cornell potential for the energy region of interest. Regarding the scalar interaction in Eq. (20b), we take a phenomenological constant term, represented by $A$, plus a linear regularized potential, analogous to the confining term of the vector potential. We note that $\beta_{s}$ is (as $\beta_{v}$) related to the confinement strength, and we also fix $b=10^{-2}~{}\text{GeV}$. In this work we use two different prescriptions for the scalar potentials to fit the experimental data. In the first prescription (potential I), we set $\beta_{s}=0$, i.e. we keep the effective constant scalar confinement term $A$ in Eq. (20b) and we omit the scalar linear regularized potential. In the second one (potential II), the value of $\beta_{s}\neq 0$ is determined fitting the data. We study these two prescriptions separately in order to assess the impact of the confinement term of the scalar interaction. In what follows we refer to these two models as

[TABLE]

II.6 Integral equation for the $c\bar{c}$ system

Now, we can write the total integral Hamiltonian equation (in momentum space) for our system. First, we introduce the standard kinetic energy of the quark and the antiquark in the charmonium CM

[TABLE]

The integral Hamiltonian equation takes the form

[TABLE]

where ${\cal H}_{\text{int}}(\vec{p}_{b},\vec{p}_{a})$ is given by Eq. (9), $M_{0}$ represents a phenomenological zero point energy of the spectrum, $M$ is the resonance mass (that is the eigenvalue of the integral equation) and $\Psi(\vec{p})$ is the resonance wave function. For brevity, in the interaction Hamiltonian, the Pauli spin matrices have not been explicitly written and also, in the wave function, the energy and angular momentum quantum numbers have been omitted.

Equation (24) can be obtained performing a three-dimensional reduction of a Bethe-Salpeter equation with an instantaneous interaction and projecting the result onto the positive energy spinors itzy .

III Solution method and fitting procedure

III.1 Solution method

In order to determine the resonance energies of the spectrum we have to solve Eq. (24) which is a Fredholm integral equation of second order Arfken . There are different available methods to solve this kind of integral equations. In particular, we employ a variational procedure. First, we express the wave function of the system by means of a superposition of the wave functions of a given trial basis. Using the standard $LS$ coupling for the angular momenta the wave functions read

[TABLE]

where the radial function in the momentum space $R_{n,L}(p;\bar{p})$ represents the trial function, $n$ is the principal quantum number, $\bar{p}$ is the variational parameter with the dimension of momentum, $Y_{L,M_{L}}(\hat{p})$ is the corresponding spherical harmonic, and $\chi_{S,M_{S}}$ is the spin function. For convenience we introduce the shorthand notation $\{\nu\}=L,S,J$ for the spin-angle quantum numbers of the basis, i.e., the orbital angular momentum ($L$), the total spin ($S=0,1$) and the total angular momentum ($J$). For the actual calculations we take the projection of the total angular momentum on the $z$-axis as $M_{J}=J$. The matrix elements involving other values of $M_{J}$ are related to the $M_{J}=J$ matrix element through the Wigner-Eckart theorem edmonds . For simplicity we do not consider the possibility that different values of $L$ mix, whose effect has been shown to be negligible in semirelativistic models Radford ; DeSanctis2 . For the radial part we select the orthonormal three dimensional harmonic oscillator wave functions

[TABLE]

where ${\cal L}^{L+\frac{1}{2}}_{n}(s^{2})$ are the generalized Laguerre polynomials, $s\equiv p/\bar{p}$ and $\bar{p}$ is the variational parameter. Now, we can diagonalize the Hamiltonian (mass) matrix in Eq. (24) obtaining

[TABLE]

where, as before, the Pauli spin matrices in ${\cal H}_{\text{int}}(\vec{p}_{b},\vec{p}_{a})$ have not been explicitly written. Notice that a six-dimensional integration over $\vec{p}_{b}$ and $\vec{p}_{a}$ appears in the interaction term. We integrate analytically the spin-angle part of the matrix elements (see details in Appendix A) so that only a three-dimensional integration over $p_{b}$, $p_{a}$ and $\cos\theta=x$ is left.

To solve Eq. (27) we use the variational diagonalization procedure introduced in DeSanctis2 . We obtain good numerical convergence for the energy eigenvalues taking the first ten trial wave functions of the basis for each state. Hence, the $10\times 10$ Hamiltonian matrix is diagonalized and minimized through the standard variational approach. Once a numerical solution is obtained, we check that the solution represents a good numerical approximation. We substitute the numerical solution in the lhs of Eq. (24) and verify that it fulfills the integral equation.

III.2 Fitting and uncertainties calculation

The two models under study defined by potentials I and II, see Eq. (22), depend on seven and eight parameters respectively. Potential I depends on $m$ (constituent mass of the $c$ and $\bar{c}$ quarks), $M_{0}$ (phenomenological zero point energy of the spectrum), $\alpha_{st}$ (effective strong coupling constant), $\beta_{v}$ (vector confinement intensity), $A$ (constant scalar interaction), and $k_{s}$ and $p_{s}$ (screening parameters), while the model with potential II depends on the same parameters plus $\beta_{s}$ that represents the scalar confinement intensity.

To determine the values of the parameters, the corresponding uncertainties and the charmonium spectrum we fit the experimental resonance masses given in Table 1, i.e. the ground states of each $J^{PC}$ band plus the $\eta_{c}^{\prime}$ and $\psi^{\prime}$ resonances, whose energies are below the open charm threshold $D\bar{D}$. In doing so, we proceed as follows Fernandez-Ramirez:2015fbq :

We randomly choose values for the masses of the resonances by sampling a Gaussian distribution according to the uncertainties given in Table 1, obtaining a resampled lowest-lying charmonia spectrum; 2. 2.

We use the least-squares method to minimize the squared distance

[TABLE]

where $M_{i}$ represent these resampled ground state experimental masses, i.e. the lowest-lying states with quantum numbers: $i\equiv J^{PC}=$ $0^{-+}$ ($\eta_{c}$ and $\eta^{\prime}_{c}$), $1^{--}$ ($J/\psi$ and $\psi^{\prime}$), $0^{++}$ ($\chi_{c0}$), $1^{+-}$ ($h_{c}$), $1^{++}$ ($\chi_{c1}$) and $2^{++}$ ($\chi_{c2}$), and $E_{i}$ are the theoretical CM energies of the lowest-lying states obtained from the potential models I and II solving Eq. (27). The fits are performed using MINUIT MINUIT .

This procedure is repeated $1000$ times providing enough statistics to compute the expected value of the parameters and the spectrum energies as well as their uncertainties at a $1\sigma$ (68%) confidence level (CL). The expected values of the parameters are computed as the mean values of the $1000$ samples. To compute the uncertainties at a $1\sigma$ CL we take the best 68% of the fits and we compute the differences between the mean value and the highest and lowest masses. In this way our uncertainties are not symmetric. Table 2 provides the fitted parameters for both potentials and Tables 3 and 4 the correlation matrices of the parameters for potentials I and II, respectively.

Table 3 (potential I), shows a strong correlation among the parameters of the vector interaction in Eq. (20a). Furthermore, we note a marked correlation between $A$ and all the parameters of the vector interaction, confirming that, in this model, vector and scalar interactions are physically correlated. Besides, it is worth highlighting the correlation of the parameters in the vector and scalar interactions with the screening parameters $k_{s}$ and $p_{s}$ of Eq. (11). We notice that the parameter $p_{s}$ turns out to be weakly correlated with the parameters of the vector interaction but strongly correlated with the parameters of the scalar interaction. For potential II, see Eq. (22), we have the additional parameter $\beta_{s}$. In this case the parameters are less correlated among them than for potential I, see Table 4. In particular, we point out that the correlations between the effective quark mass ($m$) and other parameters are significantly reduced in potential II compared to potential I, with the exception of the phenomenological zero point energy $M_{0}$. Also, the correlation between $\alpha_{st}$ and $\beta_{v}$ is largely reduced. In this way, the parameters in the vector interaction are uncorrelated among them. Regarding the scalar interaction, we observe that the parameters $A$ and $\beta_{s}$ are almost uncorrelated. However, we noticed a considerable correlation between parameters $\beta_{v}$ and $\beta_{s}$, suggesting that, in any case, a strong physical correlation exists between the confinement terms of the vector and the scalar interactions. The parameter $p_{s}$ of the screening factors is weakly correlated with the other parameters of the interactions except with the phenomenological parameter $A$ in the scalar interaction. The high correlation highlights how the screening impacts the scalar interaction. Finally, comparing the correlation matrices for both potentials, we find that the addition of the confining term in the scalar interaction reduces the correlation among some of the parameters, which allows for a more straightforward interpretation of the physics associated to each term in the potential.

For completeness, we plot the screening function $F_{s}(p)$ in Fig. 1 for the two potentials. We introduce the screening momentum $p^{\>\text{j}}_{1/2}$ ($\text{j}=\text{I},\text{II}$ labels potentials I and II) that are defined by $F_{s}(p^{\>\text{j}}_{1/2})=1/2$ (recall that $F_{s}(0)=1$). We find the values of $p^{\text{I}}_{1/2}=1.4~{}\text{GeV}$ and $p^{\text{II}}_{1/2}=1.9~{}\text{GeV}$. These values correspond to the screening kinetic energy

[TABLE]

obtaining $\bar{E}^{\text{I}}=4.1~{}\text{GeV}$ and $\bar{E}^{\text{II}}=4.7~{}\text{GeV}$. It shows that the screening effect is active above the open charm threshold where a sort of saturation energy must be introduced for the interaction Xiang .

We now consider the specific calculation of the charmonium spectrum and its theoretical uncertainties. Once that $1000$ fits have been performed, the uncertainties of the model parameters can be propagated to the charmonium spectrum using bootstrap NumericalRecipes ; Landay:2016cjw . This method allows us to carry to the spectrum computation the experimental uncertainties together with the correlations among the parameters. This method is more accurate than those using the covariance matrix, although it is computationally more expensive. In more detail, the computation of the charmonium spectrum is performed as follows: For each one of the $1000$ sets of parameters obtained from the fits we compute the spectrum. Hence, for each state of the spectrum we have a set of $1000$ values of the mass. The expected value of each state as well as the uncertainty is computed in the same way it was done for the parameters.

IV Results for the charmonium spectrum

Tables 1 (fitted states) and 5 (predicted states) report the obtained spectra for both models. Figures 2 (potential I) and 3 (potential II) provide a graphical representation of our results compared to the available experimental data from PDG PDG2016 . We can see that the obtained fits (below open charm threshold) are excellent and that, in general, the whole spectrum is well described with both potentials. We note that the statistical uncertainties grow up with the excitation energy of the states.

Table 1 shows that potential II reproduces all the fitted states within the $1\sigma$ CL while for potential I $\chi_{c1}$ and $\chi_{c2}$ fall marginally outside of the $1\sigma$ CL. Hence, potential II produces a slightly better result for the $J^{PC}$ ground states. If we look into the predicted excited states (Table 5) we see larger differences between both potentials, although the produced overall structure of the spectrum is the same.

In this Section we discuss in more detail some of the states as given by potentials I and II, according to the $n^{2S+1}L_{J}$ spectroscopical assignment of the quantum numbers.

IV.1 $X(3915)$, $X(3872)$ and $\chi_{c2}(2P)$

For these states we assign the quantum numbers $2^{3}P_{J}$, with $J=0$, $1$ and $2$ for $X(3915)$, $X(3872)$ delAmoSanchez:2010jr and $\chi_{c2}(2P)$ respectively. Potentials I and II are able to reproduce the hyperfine splitting of the fitted $1P$ states ($\chi_{c0}$, $h_{c}$, $\chi_{c1}$ and $\chi_{c2}$) but the structure of the $2P$ states cannot be reproduced with the same accuracy. In particular, the fact that the $X(3915)$ state has a higher mass than the $X(3872)$. However, we accurately predict the $\chi_{c2}(2P)$ state within uncertainties with both models. These states lay very close to the $\bar{D}D^{*}$ threshold, hence, dynamical effects and degrees of freedom not considered in our models such as threshold effects and molecular admixture in the wave functions can be responsible for the deviations from our predictions.

The $X(3872)$ state is of particular interest and requires a careful discussion. This resonance was discovered by the Belle collaboration in $B$ decays Choi and was confirmed by CDF collaboration Acosta . Various interpretations have been proposed for this resonance whose properties cannot be easily explained within a standard $q\bar{q}$ picture of mesons due to the proximity to the $\bar{D}D^{*}$ decay threshold Santo2 . It has been identified as a tetraquark Terasaki ; Faustov ; Esposito:2014rxa ; Maiani:2004vq , as a $\bar{D}D^{*}$ molecule Wong ; Baru ; Ortega:2010qq ; AlFiky:2005jd , as $c\bar{c}$ core plus molecular components Santo2 , and as a hybrid state Takizawa . Even in the context of semirelativistic quark models there is certain controversy on the assignment of the quantum numbers. In particular, besides the usual $2^{3}P_{1}$ assignment (favored by studies of the decay processes Faccini:2012zv ; Aaij13 ), a $1^{1}D_{2}$ has also been proposed DeSanctis ; Radford . Within the Bethe-Salpeter approach this resonance has been identified with a $3^{3}P_{1}$ state at a mass value of $3912~{}\text{MeV}$ Fischer but, at the same time, a spurious $2^{3}P_{1}$ state at $3672~{}\text{MeV}$ is also obtained. Other calculations with a similar formalism obtain approximately the mass value of $3900~{}\text{MeV}$, but incorporating the $X(3872)$ state in the fitting procedure Hilger:2014nma . The semirelativistic calculations in DeSanctis ; DeSanctis2 , that employ a vector and a scalar interaction, obtain a mass of $3912~{}\text{MeV}$ for the $2^{3}P_{1}$ state. The UQM obtains a mass value around $3908~{}\text{MeV}$ ref ; Santo2 ; Ferreti_alone ; Ferre-Galata . A calculation in the coordinate space with a screened potential has also given a good result of $3901~{}\text{MeV}$ Bai . In general, the nonrelativistic and semirelativistic quark models tend to overestimate the mass of this resonance with respect to the experimental value Radford ; Cao2012 ; Barnes . In the present work we obtain $3903$ and $3918~{}\text{MeV}$ with potentials I and II, respectively. These results are in reasonable but not complete agreement with the experimental value. This result together with the large isospin violation Olsen1 which prevents the identification of this state with a standard charmonium, indicates the need for a more accurate model that includes, in more precisely way, the threshold effects due to the opening of $D\bar{D}$ and $\bar{D}D^{*}$ channels. Finally, we note that both models predict a $1^{1}D_{2}$ state with mass value of $\approx 3790~{}\text{MeV}$. This resonance has not been experimentally observed yet. To complete the analysis of this multiplet, we comment on the $X(3915)$ resonance. This state has been described as a $c\bar{c}$ state Gang and also has been suggested to be a tetraquark Voloshin and a molecule Molina:2012bz . In the $c\bar{c}$ context, this resonance is identified with the $2^{3}P_{0}$ state Gang . However, there exist some issues that question the validity of this assignment. For example, it has been proposed that the $X(3915)$ is basically the $\chi_{c2}(2P)$ Zhou:2015uva ; Baru:2017fgv Among them we highlight the lack of signal in the $X(3915)\rightarrow D\bar{D}$ decay channel Olsen1 . Furthermore, the energy difference between the mass of the $X(3915)$ and its hyperfine splitting partner $\chi_{c2}(2P)$, is smaller than expected Olsen1 ; Olsen2 ; Xia . In our model the expected values of the mass for the $2^{3}P_{0}$ state are $\approx 45~{}\text{MeV}$ (potential I) and $\approx 20~{}\text{MeV}$ (potential II) below the experimental data. Hence, our model does not describe this state with high accuracy, but taking into account the general uncertainties in the hadronic models, an interpretation of this state as a $c\bar{c}$ meson cannot be excluded.

IV.2 $\psi(3770)$, $\psi(4040)$, $\psi(4160)$ and $\psi(4415)$

In the overpopulated zone of the $J^{PC}=1^{--}$ resonances above the $D\bar{D}$ energy threshold, we find some states that are reported by the PDG as charmonia: $\psi(3770)$, $\psi(4040)$, $\psi(4160)$ and $\psi(4415)$, and generally interpreted as $1^{3}D_{1}$, $3^{3}S_{1}$, $2^{3}D_{1}$ and $4^{3}S_{1}$ states respectively Olsen1 .

In our model, with both potentials, we can reproduce the mass of the $\psi(3770)$ with acceptable accuracy. The minor difference between the theoretical and experimental value ($50~{}\text{MeV}$ with potential I and $20~{}\text{MeV}$ with potential II) can be related to the proximity of the $D\bar{D}$ threshold.

The $\psi(4040)$, with an experimental mass value of $4039~{}\text{MeV}$, is generally identified as a $3^{3}S_{1}$ state Olsen1 ; Lebed . However, the theoretical mass value obtained with potential models for the $3^{3}S_{1}$ state is greater than the experimental value Barnes ; DeSanctis2 ; Cao2012 . In our case we can describe this state with sufficient accuracy thanks to the inclusion of the screening functions Bai . In particular, potential I provides a theoretical prediction that agrees with the experimental value within uncertainties. Potential II also produces an acceptable mass value for this state.

The $\psi(4160)$ is interpreted as a radial excitation of the $\psi(3770)$. The predicted mass value is lower than the experimental one, a situation similar to what happens in other screened potential models Bai ; Vijande . If we accept that the $\psi(4160)$ can be considered as a pure $2^{3}D_{1}$ state, we have to conclude that our model and the screened potential models are unable to reproduce this state, supporting a more elaborated nature for this state Barnes-Close ; Asner_etal .

Finally, the resonance $\psi(4415)$ has been generally interpreted as a $4^{3}S_{1}$ state Olsen1 but also as a $5^{3}S_{1}$ state in Bai where the masses of the $\psi$ states are lowered by the inclusion of the screening effects. Within our model it may be tentatively described in two different ways depending on which potential is used (see Table 6). Potential I suggests that we can roughly identify this resonance with a $4^{3}D_{1}$ state but potential II suggests a $3^{3}D_{1}$ state. These unusual tentative assignments for $\psi(4415)$ differ from the standard ones of the semirelativistic potential models Isgur ; Radford where this resonance is well described as a $4^{3}S_{1}$ state. We note that in our model the $4^{3}S_{1}$ state has a mass of $4311~{}\text{MeV}$ with potential I and $4360~{}\text{MeV}$ with potential II, that are better matches for the $X(4260)$ [$Y(4260)$] and the $X(4360)$ [$Y(4360)$] resonances, respectively. Hence, our $4^{3}S_{1}$ state cannot be identified as the $\psi(4415)$.

IV.3 $X(4260)$, $X(4360)$ and $X(4660)$

These resonances with $J^{PC}=1^{--}$ were discovered by BaBar B.Aubert1 and Belle C.Yuan ; Wang1 collaborations and confirmed by a combined analysis of the data Z.Liu . There exist different interpretations of the nature of these resonances: hybrid charmonia Zhu ; Guo:2008yz , tetraquarks Wang ; Galkin-Fau , molecules Close , $c\bar{c}$ mesons Ding and some authors suggest that these signals are the product of the interference between different channels and resonances Chen0 ; Chen1 ; Chen2 . In the charmonium context, these vector mesons are atypical. Their quantum numbers are univocally determined by their production mechanism but some of their features, e.g. the absence of open charm production Galkin-Fau , do not match the $c\bar{c}$ picture brambilla ; Olsen1 .

Our two potentials produce different results for these states. We note that, in our model, the predictions in the high energy region of the spectrum depend strongly on the form of the scalar interaction. Potential I predicts a $4^{3}S_{1}$ state, i.e. $\psi(4S)$, that can be roughly identified with the $X(4260)$ [$Y(4260)$] resonance. This assignment is reasonable if we compare with the bottomonium system. In particular, we compare the experimental mass difference $\Upsilon(4S)-\Upsilon(3S)\approx 224~{}\text{MeV}$ with that of $\psi(4S)-\psi(3S)\approx 254~{}\text{MeV}$. In Liu-Chen a mass of $\approx 4260~{}\text{MeV}$ was obtained for the $\psi(4S)$ state. Also, as in Bai , in our model we can match the $3^{3}D_{1}$ state to the $X(4360)$ resonance. Finally, we note that potential I cannot reproduce the $X(4660)$ [$Y(4660)$] resonance. Potential II depicts a completely different situation. We identify the resonances $X(4360)$ [$Y(4360)$] and $X(4660)$ with the $4^{3}S_{1}$ and $4^{3}D_{1}$ states respectively but we cannot describe the $X(4260)$ state.

IV.4 $X(3940)$ and $X(4160)$

These states were observed by Belle collaboration in double charm production processes Abe:2007sya ; Abe ; Shen . Their decays suggest the $J^{PC}=0^{-+}$ quantum numbers. They have been described in many different ways: hybrids Petrov , molecules Molina:2009ct ; Fernandez , and as the second ($\eta_{c}^{\prime\prime}$) and the third ($\eta_{c}^{\prime\prime\prime}$) radial excitations of the $\eta_{c}$ meson in a standard $c\bar{c}$ picture KTChao ; Estia ; Ping . Both semirelativistic and nonrelativistic potential models Barnes obtain masses for $\eta_{c}^{\prime\prime}$ and $\eta_{c}^{\prime\prime\prime}$ above the experimental values of $X(3940)$ and $X(4160)$, which casts doubts on these assignments.

In Liu-Chen it is assumed that the energy gap between $\psi(3S)-\psi(2S)\approx 353~{}\text{MeV}$ resonances is similar to the $\eta_{c}(3S)-\eta_{c}(2S)$ gap, making it possible to estimate the energy of $\eta(3S)$. The authors found that $\eta(3S)\approx 3992~{}\text{MeV}$, roughly close to the experimental mass of the $X(3940)$. This allows us to identify it as a $3^{1}S_{0}$ state. On the other hand, through a decay analysis in Liu-Chen it is concluded that it is not possible to identify the $X(4160)$ resonance with a $4^{3}S_{0}$ state.

In the screened potential model Bai the correspondence between the $X(3940)$ resonance and the $3^{1}S_{0}$ state is slightly improved. They found $\eta_{c}(3S)=3991~{}\text{MeV}$, but the difference between the experimental mass of the $X(4160)$ and the theoretical mass value of the $4^{1}S_{0}$ state ($4250~{}\text{MeV}$) is large enough to call into question this assignment. In our model, with potential I and the assignments $3^{1}S_{0}$ and $4^{1}S_{0}$, we obtain the values of $4030~{}\text{MeV}$ and $4302~{}\text{MeV}$ for $X(3940)$ and $X(4160)$, respectively. Using potential II, we obtain for the same states the values of $4043~{}\text{MeV}$ and $4349~{}\text{MeV}$, respectively. We see in Table 5 and in Fig. 2 that with potential I we can obtain a rough description of $X(3940)$. Nevertheless, in general, we conclude that our model is not able to describe appropriately these two resonances.

IV.5 $X(4140)$

This resonance was first observed in Aaltonen1 and later confirmed Aaltonen2 . It has been interpreted as a tetraquark Fang ; Stancu ; Wang:2015pea , as a molecule Wang:2009ue ; Wang:2014gwa ; Xiang , as a hybrid charmonium Mahajan and as a $c\bar{c}$ state Gonzalez . Recently, the LHCb collaboration Aaij reported the $1^{++}$ quantum numbers for this resonance, which leads to a $3^{3}P_{1}$ state interpretation, i.e. $\chi_{c1}(3P)$, within the $c\bar{c}$ context Chen:2016iua . The mass value of the $3^{3}P_{1}$ state in the relativistic potential Barnes and screened potential Bai models is above of the mass of the $X(4140)$ [$Y(4140)$] resonance, but, for the case of the latter the difference is just $40~{}\text{MeV}$, so this interpretation is possible. In our case, the $3^{3}P_{1}$ state is obtained with higher values than the experimental data found by LHCb for both potentials. More precisely, potentials I and II predict this state approximately $30~{}\text{MeV}$ and $70~{}\text{MeV}$ above the experimental value, respectively.

IV.6 $\psi(3823)$

This resonance was observed for the first time more than 20 years ago Antoni . Recently, it was observed again by Belle collaboration Bhardwaj and confirmed by BESIII collaboration Ablikim . In agreement with semirelativistic potential models Barnes ; Radford , our model describes reasonably well this state with both potentials, $\approx 3785~{}\text{MeV}$ versus an experimental value of $3822~{}\text{MeV}$, with the standard assignment $1^{3}D_{2}$ ($J^{PC}=2^{--}$). In Bai , it was found a mass value of $3798~{}\text{MeV}$, similar to ours. In the UQM Santo2 a mass of $3736~{}\text{MeV}$ was found, significantly lower than the experimental value.

IV.7 Missing resonances

Our model, as the majority of the $c\bar{c}$ models do, predicts (with both potentials) some states that have not been observed in the experiments yet. In particular, we note the states $1^{1}D_{2}$, $2^{1}P_{1}$, $3^{3}P_{0}$ and $1^{3}D_{3}$. The $1^{1}D_{2}$ state, as already discussed, does not represent a good candidate for the $X(3872)$ and constitutes a missing resonance produced by our model. The $2^{1}P_{1}$ state belongs to an energy region where only the charged mesons $Z_{c}(3900)$ and $Z_{c}(4010)$ have been found. This region is currently empty of charmonialike states, while the radial excitations of the $h_{c}$ state should lie there. The $3^{3}P_{0}$ state that was identified as $X(4140)$ in the past becomes a missing resonance (see Section IV.5) if we consider the latest results from LHCb Aaij . The last missing resonance, the $1^{3}D_{3}$ state, also remains undetected despite that its energy value is close to the open charm threshold.

V Summary and conclusions

We have studied the charmonium spectrum employing a $c\bar{c}$ relativistic Dirac potential model in momentum space with a vector and a scalar interaction. We note that only by a combination of scalar and vector confining potentials is it possible to eliminate the empirically contraindicated spin-orbit interaction in the interquark potential. Positive energy Dirac spinors have been used without performing any nonrelativistic expansion.

In our model, we have employed for the vector interaction a one gluon exchange term plus a confining term. For the latter term we used a standard regularized linear confining interaction transformed to momentum space. For the scalar interaction we have considered two possibilities: (i) potential I, which is given by a constant term only, namely $A$, that accounts for an effective confinement term; and (ii) potential II, that together with the constant term $A$ it also incorporates a standard regularized linear confining term. We have also incorporated phenomenological screening factors that take into account the coupling of the $c\bar{c}$ system with virtual $q\bar{q}$ states. Potential I depends on seven parameters: the $c$ and $\bar{c}$ quark mass $m$, the zero point energy $M_{0}$, the strong coupling constant $\alpha_{st}$, the strength of the confining term in the vector interaction $\beta_{v}$, the constant term of the scalar interaction $A$, and the parameters of the screening factors $k_{s}$ and $p_{s}$. Potential II depends on eight parameters: the same ones as potential I plus the strength of the confining term of the scalar interaction $\beta_{s}$.

Both potential models have been fitted to the experimental masses of eight resonances below the open charm threshold (see Table 1): $\eta_{c}$, $J/\psi$, $\chi_{c0}$, $\chi_{c1}$, $h_{c}$, $\chi_{c2}$, $\eta^{\prime}_{c}$ and $\psi^{\prime}$. In both cases we obtain a good description of these states, including the unusual splitting between $h_{c}$, $\chi_{c1}$ and $\chi_{c2}$ resonances. The uncertainties in the parameters have been computed using bootstrap technique. In this way we can propagate exactly the statistical errors in the data to both the parameters and the resonance masses taking into account all the correlations. We predict the energies of the resonances above the open charm threshold and compute the associated uncertainties (see Tables 5 and 6 and Figs. 2 and 3). We predict correctly the structure of the charmonium spectrum and we obtain a reasonable agreement with most of the available experimental data. The screening effect turns out to be relevant in the description of the spectrum with both potentials, especially at energy values above the open charm threshold. The correlation matrices of the parameters for potential I (Table 3) and II (Table 4) show that there is a small correlation between the screening parameters and the vector interaction parameters. However, the correlation is strong with the parameters of the scalar interaction. Hence, we confirm that the screening factors in our model take phenomenologically (and only partially) into account the excitation of $q\bar{q}$ degrees of freedom. As the excitation energy grows, it becomes more difficult to reproduce the data, indicating that the presence of new degrees of freedom that should be introduced in the theoretical model. This point requires a new and more comprehensive investigation.

We have performed a full statistical error analysis, determining the uncertainties of the parameters and their correlations. We have exactly propagated those uncertainties and correlations to the predicted spectrum. As far as we know, previous error analysis within phenomenological models have been very limited and incomplete. To perform the error analysis we use the bootstrap technique. This method is computationally taxing but provides a rigorous treatment of the statistical uncertainties. Rigorous error estimations allow us to assess if the inclusion of a new effect in the phenomenological model is necessary or not. Moreover, the study of the correlation among the parameters that arises from the error analysis allows us to identify how independent are the different pieces of the model among them. A full error analysis is mandatory to identify which deviations from experimental data can be absorbed into the statistical uncertainties of the models and which can be related to physics beyond the $c\bar{c}$ picture, guiding future research.

In more detail, the $X(3915)$, $X(3872)$ and $\chi_{c2}(2P)$ are identified as $2^{3}P_{J}$ states with $J=0$, $1$, and $2$, respectively. The $\chi_{c2}(2P)$ is accurately reproduced. The $X(3872)$ mass is slightly overestimated ($3904~{}\text{MeV}$ with potential I and $3918~{}\text{MeV}$ with potential II) signaling that this state is mostly a $c\bar{c}$ state with its mass dynamically modified by effects not taken into account in our models, e.g. molecular components or open channel effects. The $X(3915)$ is not well described, and, hence, the $2P$ splitting among these states is not properly accounted for.

We also recover the general structure of the $\psi$ states above the $D\bar{D}$ threshold, namely $\psi(3770)$ ($1^{3}D_{1}$), $\psi(4040)$ ($3^{3}S_{1}$), $\psi(4160)$ ($2^{3}D_{1}$) and $\psi(4415)$ ($4^{3}D_{1}$ with potential I and $3^{3}D_{1}$ with potential II). However, while $\psi(3770)$ and $\psi(4040)$ are sufficiently well described, $\psi(4160)$ and $\psi(4415)$ are not.

Our two models produce different results for the $X(4260)$ [$Y(4260)$], $X(4360)$ [$Y(4360)$] and $X(4660)$ [$Y(4660)$] resonances. Potential I produces a $4^{3}S_{1}$ state that can roughly be identified with the $X(4260)$ resonance. The higher-lying $X(4660)$ resonance is reproduced within the uncertainties of the model with potential II but not with potential I, and the $X(4360)$ resonance is reproduced with both potentials, but within very large theoretical uncertainties. We consider that we are not able to provide a good description of $X(3940)$ and $X(4160)$ resonances with our model, an indication that these states might be dominated by dynamical effects that go beyond the $c\bar{c}$ picture. The $X(4140)$ resonance is approximately reproduced with potential I, suggesting a correspondence with a $c\bar{c}$ $3^{3}P_{1}$ state; however, it is not properly reproduced by potential II despite the large uncertainties of the theoretical prediction. The $\psi(3823)$ state is also well reproduced, in agreement with semirelativistic potential models.

Further investigation should be devoted to understand confinement in a more complete way. We note that potential II, which includes a confining term in both the vector and the scalar interactions, provides slightly better results than potential I (which lacks the scalar confining term), in agreement with the standard argument from spin-orbit reduction that the vector and the scalar potentials have approximately the same form. In addition, the lower correlation among the parameters in potential II, Table 4 compared to the correlation among the parameters in potential I, Table 3, allows usto make a more straightforward connection to the physics associated to each term in potential II than in potential I, in particular the confining interaction term.

The inclusion of relativistic effects in $q\bar{q}$ potential models constitutes an important step forward in the description of charmonia, especially the higher-lying states. However, it is apparent that a deeper dynamical study is also needed. One should use interactions more straightforwardly related to QCD and include other effects such as molecular admixture, tetraquark components, as well as open channels and threshold effects have to be taken into account to obtain a complete description of the experimental states, in particular those whose difference with the quark model prediction cannot be accounted for by the statistical uncertainties.

Acknowledgements.

The authors thank R. Bijker for useful discussions and comments. This paper is part of the PhD thesis of D.M., supported by Departamento Administrativo de Ciencia, Tecnología e Innovación (Colciencias) by means of the loan-Scholarship Program No.6172. D.M. thanks Instituto de Ciencias Nucleares (ICN) at Universidad Nacional Autonóma de Mexico (Mexico City, Mexico) for the hospitality extended to him during the first semester of 2016 and, in particular, for the access to the Tochtli-ICN cluster that was used to perform part of the numerical calculations of this work. The work of C.F.-R. is supported in part by PAPIIT-DGAPA (UNAM, Mexico) Grant No. IA101717, by CONACYT (Mexico) Grant No. 251817 and by Red Temática CONACYT de Física en Altas Energías (Red FAE, Mexico). M.D.S. gratefully thanks INFN Sezione di Roma 3 and Department of Physics of Roma Tre University (Rome, Italy) for the hospitality extended to him during his sabbatical year in $2015--16$ and, in particular, for the use of the computational facilities of the Cluster Roma Tre.

Appendix A Spin-angle matrix elements

To calculate the spin-angle matrix elements we need to factorize the corresponding operators. To this aim we introduce the following quantities

[TABLE]

and

[TABLE]

As an example, we first consider the vector term of the interaction. Then, we generalize the calculation to the case of the vector and scalar interactions of the model. Substituting the Dirac four-currents of Eq. (17) in Eq. (16) and using the definitions in Eqs. (30) and (A), we can rewrite the vector Hamiltonian as follows

[TABLE]

We use the notation from Section II.2 and rename

[TABLE]

where $\bar{{\cal H}}^{(v)}(\vec{p}_{b},\vec{p}_{a})$ is defined in Eq. (10). Recalling Eq. (4) for $\vec{q}^{\>2}$, we can express the rhs of Eq. (33) as the product of the functions $\bar{{\cal H}}^{(v),k}(p_{b},p_{a},x)$ times the operators depending on $\hat{p}_{a},~{}\hat{p}_{b},~{}\vec{\sigma}_{1}$ and $\vec{\sigma}_{2}$. These operators are denoted as $\Lambda_{k}(\hat{p}_{a},\hat{p}_{b},\vec{\sigma}_{1},\vec{\sigma}_{2})$. We obtain a sum of eight terms

[TABLE]

The explicit expressions of the first four terms of this sum are

[TABLE]

Analogously, the other terms of the vector interaction and those of the scalar interaction can be found through straightforward calculations. Next, we perform a double spherical harmonic expansion on both the scalar and vector terms with respect to the angles of $\hat{p}_{b}$ and of $\hat{p}_{a}$ for each $\bar{\cal H}^{(\tau),k}(p_{b},p_{a},x)$ term, where $\tau$ labels the vector ($v$) and scalar ($s$) interactions. We obtain

[TABLE]

where we have expanded in Legendre polynomials $P_{\ell}(x)$

[TABLE]

The second line of Eq. (43) represents the spin-angle operators that can be calculated analytically from the angular part of the wave functions of our basis given in Eq. (25). In doing so, we write

[TABLE]

where $d\Omega_{a}$ and $d\Omega_{b}$ represent the integration elements over the angles of $\hat{p}_{a}$ and $\hat{p}_{b}$, respectively. The matrix elements of the previous equation can be calculated in a standard way by using angular momentum algebra.

Finally, the interaction matrix elements, that is the third term in the rhs of Eq. (27) can be written as

[TABLE]

where two integrations over $p_{a}$, $p_{b}$ and one integration over $x=\cos\theta$ from Eq. (44) remain.