Study of charmonium(-like) mesons under a diabatic approach

Studies on charmonium(-like) states play a significant role in understanding how quarks constitute hadrons, which is closely related to the non-perturbative behavior of strong interactions. Since the discovery of the first charmonium meson $ J/\psi $[1, 2] in 1974, a series of charmonium mesons, such as $ \psi(3686) $ [3], $ \chi_{c0}(1P) $ [4], $ \chi_{c2}(1P) $ [5], and others [6, 7], have been successively discovered experimentally. To understand the properties of these states, various potential models have been proposed to describe the interactions between quarks [8, 9]. Among them, the Cornell potential [10] is widely adopted to predict and interpret the charmonium spectrum. While the Cornell potential is rather successful in describing the low lying charmonium states, its predictions for many highly excited states are still controversial and not well verified experimentally. Meanwhile, in the past two decades, a series of charmonium-like states called XYZ states have been discovered experimentally[11], whose masses differ significantly from the predictions of the potential model [12]. Some of them clearly do not fit with the conventional $ q\bar q $ picture of meson state (see Refs. [13−16] for recent reviews).

For the newly discovered charmoinum-like states, a special feature is that the masses of many states lie very close to some meson pair thresholds, which have raised questions over their molecular nature. In general, for higher excited states, the closeness to the open charm thresholds may significantly affect the charmonium spectrum through the virtual meson loop. To describe these states, it is then necessary to take into account such coupled-channel effects. Various approaches have been developed along this line [8, 17−20]. One solution is to introduce the adiabatic approximation.

The adiabatic approximation, also known as the Born-Oppenheimer (B-O) approximation, was initially developed to describe molecular systems [21] and has since found extensive applications in molecular and atomic physics as well as in the meson states in quantum chromodynamics (QCD) [22, 23]. However, the B-O approximation relies on two key assumptions: the adiabatic approximation and single-channel approximation [24−26]. As mentioned above, for charmonium states that are close to the open flavor meson-meson threshold, the mixing between the $ Q\bar Q $ and meson-meson configurations cannot be ignored, which challenges the validity of the single-channel approximation assumed in the B-O approximation. In this case, the system should be described by a set of coupled-channel Schrödinger equations containing non-adiabatic coupling terms (NACTs) [27], which are difficult to solve in practice [24]. To overcome these problems, it is necessary to generalize the B-O approximation and adopt the so-called diabatic formalism.

The diabatic formalism has many advantages. Firstly, in this approach, the dynamics are described by a diabatic potential matrix, which in principle can be calculated directly in lattice QCD and provides a non-perturbative method to consider the mixing of different configurations. Therefore, it offers a way to study the heavy quarkonium(-like) states based on first principles or models inspired by the lattice QCD results. In practice, once the diabatic potential matrix is determined, one may use numerical methods to solve the coupled-channel Schrödinger equation to obtain the solution of the system. In addition, by choosing an appropriate expansion basis, it is possible to make the wave function components correspond to a pure $ Q\overline{Q} $ state or the meson-meson state. In this way, the physical meaning of each channel is intuitive, which makes it more convenient to analyze the results.

Previous studies [24, 28] have applied the diabatic approach to study the quarkonium spectrum with masses below 4.1 GeV. Using Cornell potential to describe the $ c\overline{c} $ interaction, they calculated the charmonium meson spectra of the $ 0^{++} $, $ 1^{++} $, $ 2^{++} $, and $ 1^{--} $ states. However, in these works, the spin-dependent interaction between heavy quarks was not considered. The absence of spin-dependent terms prevented the authors from exploring the fine and hyperfine structures in the charmonium spectrum, which is important for a complete and accurate understanding on the charmonium spectrum.

In this work, we extend the works in Refs. [24, 28] by further considering the spin-dependent interactions between the $ c\overline{c} $, which cause the fine and hyperfine splittings of the chamonium mesons. By solving the coupled-channel Schrödinger equations, we obtain the masses of the charmonium spectrum below 4.1 GeV and the probabilities of finding $ c\bar c $ or various meson pair components in those states. Compared to the results without considering the spin-dependent interactions, we can now obtain more states for comparison with experimental data. Furthermore, this leads to a better description of lower lying states in the charmonium spectrum. Because this formalism provides a unified description of both conventional and unconventional heavy meson states, we can calculate and describe the states like $ \chi_{c0}(3860) $, whose nature is still disputed [29, 30]. By analyzing the probability of different configurations, we can explore whether a state is conventional or not. Hence, our results may provide new insights on some controversial charmonium(-like) states.

The remainder of this paper is organized as follows. In Sec. II, we review the formalism of the diabatic approach, which was developed from the B-O approximation. In Sec. III, we construct the potential matrix of the system while considering spin-dependent terms. In Sec. IV, we apply the diabatic approach to study the charmonium spectrum and discuss the results. Finally, we summarize our conclusions in Sec. V.

The Hamiltonian for a system of heavy quarkonium state can be written as [24]

where $ \mu_{Q\overline{Q}} $ is the reduced mass of the $ Q\bar Q $ system, and $ {\boldsymbol{p}}({\boldsymbol{P}}) $ is the $ Q\bar Q $ relative(total) momentum. Here, we use $ H^\mathrm{light} $ to represent the Hamiltonian describing the energy of light fields like light quarks and gluons and their interaction with $ Q\bar Q $. The Schrödinger equation of the system in the center-of-mass frame of $ Q\bar Q $ can be written as

In the heavy quarkonium system, the mass of the heavy quark is significantly higher than the energy scale of the light fields, which justifies the neglect of the kinetic energy of the heavy field, i.e., using the static limit, when considering the dynamics of the light fields. As a result, the separation of the heavy quarks $ {\boldsymbol{r}} $ can be regarded as a c-number parameter rather than a dynamical operator. Hence, the Hamiltonian $ H^{\mathbf{light}} $ becomes an operator, $H_{\rm static}^{\mathrm{light}}({\boldsymbol{r}})$, solely pertaining to the light field, and now $ {\boldsymbol{r}} $ is a parameter denoting the separation of the heavy quarks. Thus, in the static limit, the dynamics of light fields in the $ Q\bar Q $ rest frame are described by

where $ V_i({\boldsymbol{r}}) $ is the eigenvalue of the reduced Hamiltonian $ H_{\mathrm{static}}^{\mathrm{light}}({\boldsymbol{r}}) $.

To solve Eq. (2), we apply the diabatic expansion to the eigenstates $ |\psi\rangle $, as in Ref. [24] (also see Refs. [25, 26]). In brief, we use $ |\zeta_i({\boldsymbol{r_0}})\rangle $, which is the i-th eigenstate of $H^{\rm light}_{\rm static}({\boldsymbol{r_0}})$ with $ {\boldsymbol{r_0}} $ being a free parameter and taking some fixed value, to expand the state $ |\psi\rangle $ and get

where $ \tilde{\psi}_i({\boldsymbol{r}}^{\prime},{\boldsymbol{r}}_0) $ is the wave function corresponding to the i-th light field state, and $ {\boldsymbol{r}}^{\prime} $ represents the separation between the heavy quarks. Substituting Eq. (4) into Eq. (2) and multiplying $ \langle{\boldsymbol{r}}| $ and $ \langle\zeta_j({\boldsymbol{r_0}})| $ on the left side gives

where the diabatic potential matrix is defined as

The introduction of the diabatic potential matrix is a crucial step in the diabatic approach. Because $ r_0 $ is a free parameter, it is convenient to choose an appropriate $ {\boldsymbol{r_0}} $ such that the light field configuration, i.e., each state $ |\zeta_i({\boldsymbol{r}}_0)\rangle $, corresponds to a pure $ Q\overline{Q} $ state or the meson-meson state. In this way, the physical meaning of each state is clear, which makes the analysis more intuitive. For later convenience, we relabel the $ Q\overline{Q} $ state as the 0 state and use n with its value starting from 1 to label the n-th meson-meson pair state considered in this work. With these substitutions, we then have

to denote various states. The corresponding wave functions are denoted as

The matrix element of the $ Q\bar Q $ interaction is

The matrix elements describing the interactions of the meson-meson pairs are

where $ 1\leq i,j \leq N $, and N is the total number of $ M\bar M $ states considered in this work. The matrix element of the mixing potential is denoted as

where $ 1\leq j \leq N $ and $ V_{0j}({\boldsymbol{r}},{\boldsymbol{r}}_0)= V_{j0}({\boldsymbol{r}},{\boldsymbol{r}}_0) $.

With the above notation, Eq. (5) can be rewritten in matrix form as

where $ \mathbf{K} $ is the matrix composed of kinetic energy terms, $ \mathbf{V}(r) $ is the potential matrix, and $ \mathbf{\Psi}(r) $ is a column vector of the wave functions. The explicit form of $ \mathbf{K} $ is

where $ \mu_{c\bar c} $ is the reduced mass of $ c\overline{c} $, and $ \mu_{M\overline{M}}^{(i)} $ ($ 1\le i\le N $) is the reduced mass of the i-th meson-meson pair.

Analogous to lattice QCD studies [31], we neglect interactions between different meson-meson components, which yields $ V_{ij}({\boldsymbol{r}},{\boldsymbol{r}}_0)=0 $ for $ i\ne j $. The explicit form of $ \mathbf{V}(r) $ can be given as

Then, $ \mathbf{\Psi}(r) $ is defined as

where $ \psi_{c\overline{c}}(r) $ is the wave function of $ c\overline{c} $, and $ \Psi_{M\overline{M}}^{(i)}(r) $ ($ 1\le i\le N $) is the wave function of the i-th meson-meson component. The normalization condition satisfied by the wave function is

where we have defined the probility of finding the $ c\bar{c} $ and $ M\bar{M}_i $ components in the state as

and

In the previous section, we presented the main formalism for the diabatic approach. In this section, we give the specific form of the matrix elements of $ V_{ij}({\boldsymbol{r}},{\boldsymbol{r}}_0) $, which are needed to solve the Schrödinger equation of the system.

In this part, we discuss the interactions of charm quarks. The element $ V_{c\bar c}(r) $, as can be seen from Eq. (9), describes the static energy of the light field state corresponding to a pure $ c\bar c $ state. Therefore, we can take $ V_{c\bar c}(r) $ as the conventional $ c\overline{c} $ potential [22].

There are many models for the effective potential of $ c\overline{c} $, among which the non-relativistic potential model is relatively simple and widely used. The main part of this model contains a color Coulomb potential and confinement potential [32]. In this paper, we use the linear potential as the confinement potential. Hence, the central potential can be written as

where $ \alpha_s $ and b are model parameters.

In addition, we will further consider spin-dependent interactions in the potential in this work, which were not considered in Refs. [24, 28]. Following the method in Ref. [33], we introduce three spin-dependent terms in the Hamiltonian.

Firstly, we consider the spin-spin contact hyperfine potential. This interaction is one of the spin-dependent terms predicted by one gluon exchange (OGE) potential. In this work, we take it as the Gaussian-smeared form

where $ \vec{S}_{c} $ and $ \vec{S}_{\bar{c}} $ are spin operators acting on the spin of quarks and antiquarks, respectively. We take $\tilde{\delta}_\sigma(r)\; = \; (\sigma/\sqrt{\pi})^3{\rm e}^{-\sigma^2r^2}$ [12], and σ and $ m_c $ are model parameters.

The remaining spin-dependent interactions are from the spin-orbit and tensor couplings between the $ c\bar c $, which can also be deduced from the OGE potential. Using the leading-order perturbation theory, the spin-orbit potential $ V_{SO} $ and tensor potential are obtained as

where $ \vec{L} $ is the orbit angular momentum operator, $ \vec{S} $ is the total spin operator, and T is the tensor operator.

In the $ |J,L,S\rangle $ basis, the elements of the spin-orbit and spin-spin operators are diagonal. Furthermore, the off-diagonal elements of the tensor term are very small and can be neglected [34]. Thus, the matrix elements involving spin operators can be given as [12]

After including the spin-dependent terms, the total interaction potential becomes

The values of the parameters ($ m_c,\alpha_s,b,\sigma $) in the potential will be discussed in the next section. Their adopted values are listed in Table 3.

Parameter	Value	Parameter	Value
$ m_c $/GeV	1.4830	$ r_{\rm cut} $/fm	0.202
$ \alpha_s $	0.5461	ρ/fm	0.3
b/Ge$ \mathrm{V^2} $	0.1425	Δ /GeV	0.116
σ/GeV	1.1384

Table 3.  Parameters adopted in the potential matrix.

Next, we consider the interaction potential within meson-meson pairs in Eq. (10). As mentioned in Section II, one advantage of the diabatic approach is that taking an appropriate value for $ r_0 $, each state $ |\zeta_i({\boldsymbol{r_0}})\rangle $($ 1\le i\le N $) corresponds to a concrete pure meson-meson configuration. Thus, the matrix elements of $ V_{M\overline{M}} $ have clear physical meanings. Its diagonal elements represent the internal interaction of each meson-meson pair, and the off-diagonal elements represent the couplings between different meson-meson pairs.

For the off-diagonal elements of $ V_{ij}^{M\bar M} $, as previously mentioned, we set them to 0, i.e., neglecting the couplings between different meson-meson pairs. That is,

For the diagonal elements, we use the potential used in Ref. [24], i.e.,

where $ m_{M_1} $ and $ m_{\overline M_2} $ represent the masses of the mesons in the meson-meson pair. This means that we also do not consider the interactions between mesons here.

Next, we discuss the mixing potential $ V^{\rm mix}(r) $, which couples $ c\overline{c} $ with meson-meson pairs. The mixing potential used in this work is extracted from lattice QCD results [31] following the approach of Ref. [24]. Let us consider the coupling of $ c\overline{c} $ and the ith meson pair. According to Ref. [24], there exists a crossing radius $ r^c $such that the interaction potential of $ c\overline{c} $ equals the threshold mass of the ith meson pair. Thus, at $ r=r_i^c $, we have

Furthermore, the mixing potential $V^{\rm mix}$ satisfies the following conditions. Firstly, when $ r=r^c $, the mixing potential $V^{\rm mix}(r)$ reaches its maximum value. Secondly, when r is far from $ r^c $, the mixing potential quickly approaches 0. Therefore, a reasonable choice for $V^{\rm mix}$ is a Gaussian function, and it can be taken as

where Δ and ρ are parameters representing the maximum value and width of the Gaussian function, respectively. Here, b is the same parameter as in Eq. (19). The values for the parameters ρ and Δ are discussed and given below. Therefore, we get the matrix elements of mixing potential in Eq. (6) as follows:

As mentioned above, the mixing of $ c\bar{c} $ and $ M\bar{M} $ primarily occurs near the crossing radius $ r^c $. In this region, the long-range part of $ V_{c\overline{c}}(r) $ (Eq. (27)) will play a major role in the mixing potential $ V^{\mathrm{mix}}_i(r) $. In our model, the spin-dependent terms are mainly short-range, and the mixing of different $ c\bar{c} $ states due to the tensor terms is also neglected. Therefore, the inclusion of the spin interactions within the $ c\bar{c} $ system shall not alter the form of the mixing potential $ V^{\mathrm{mix}}_i(r) $ in Eq. (32).

In this section, we present the numerical results and discuss them. Firstly, in part A, we discuss some ingredients needed to describe the charmonium states and the values of parameters in this work. Then, we present the calculated results and discussion in part B.

By using the method and formalisms discussed above, we can now study the charmonium(-like) mesons while considering the coupling of $ c\bar{c} $ with meson-meson pairs. Solving the Schrödinger equation by numerical methods, we can obtain the mass spectrum of the charmonium(-like) mesons. In Table 1, we list the various meson-meson pairs considered in this work and their corresponding threshold masses.

$ M\bar{M} $	$T^{M\bar{M} }/{\rm MeV}$
$ D\bar{D} $	3730
$ D\bar{D}^*(2007) $	3872
$ D_s^+D_s^- $	3937
$ D^*(2007)\bar{D}^*(2007) $	4014
$ D_s^+\bar{D}_s^{*-} $	4080
$ {D}_s^{*+}\bar{D}_s^{*-} $	4224

Table 1.  Low-lying open charm meson-meson channels considered in this work and their thresholds [35].

Heavy quark meson states are classified by their quantum numbers, which include isospin (I), G-parity (G), total angular momentum (J), parity (P), and charge conjugation (C). In this paper, we focus on the isoscalar meson states, i.e., $ I=0 $, which also means $ G=C $.

Note that the considered diabatic potential is centrally symmetric, and there are no off-diagonal elements in the $ |J,L,S\rangle $ basis for the spin-dependent operators (here, $ J,L,S $ represent the total angular momentum, orbital angular momentum, and total spin of the system, respectively). Therefore, for a system of $ c\bar{c} $ or meson-meson states, as long as we know the $ J, L, S $ of the system, we can deduce its $ J^{PC} $ quantum numbers. The coupling of $ c\bar{c} $ with meson-meson states can only happen when they share the same $ J^{PC} $. In Table 2, we list the $ J^{PC} $ quantum numbers and corresponding orbital angular momenta of the $ c\bar c $ and meson-meson states that can couple with each other.

$ J^{PC} $	$ l_{c\overline{c}} $	$l_{D_{(S)}\bar{D}_{(S)} }$	$l_{D_{(S)}\bar{D}_{(S)}^*}$	$l_{D_{(S)}^*\bar{D}_{(S)}^*}$
$1^{--}$	0,2	1	1	1,3
$ 2^{++} $	1,3	2	2	0,2,4
$ 1^{++} $	1		0,2	2
$ 0^{++} $	1	0		0,2
$ 0^{-+} $	0		1	1
$ 1^{+-} $	1		0,2	0,2

Table 2.  Quantum numbers of the charmonium(-like) mesons considered in this work. The corresponding orbital angular momenta l of various components are also given, and a blank space means that no possible orbital angular momentum exists.

For a meson state with quantum numbers $ J^{PC} $, by reading the possible configurations from Table 1, we can construct the Schrödinger equation with Eq. (12) for the system. The obtained coupled-channel Schrödinger equation describes the corresponding heavy meson system.

Next, we need to discuss the values of the parameters used in this work. At present, we have six parameters in the potential matrix: $ m_c $, $ \alpha_s $, b, σ, ρ, and Δ. Besides these six parameters, we still need one more new parameter $ r_{\rm cut} $ to cure the singularity problem due to the ${1}/{r^3}$ term in the potential as $ r\to 0 $. Here, we follow the approach of Ref. [36] by taking ${1}/{r^3}={1}/{r_{\rm cut}^3}$ in the region of $ 0<r<r_{\rm cut} $ to resolve this problem. Therefore, we have seven parameters in total. For $r_{\rm cut}$ and the four parameters ($ m_c,\alpha_s,b,\sigma $) associated with the $ c\overline{c} $ interaction, we take their values from Ref. [36], where the same $ c\bar c $ potential was adopted and the parameters were determined by fitting the masses of 12 well-established $ c\bar c $ states.

For the parameters ρ and Δ, which appear in the mixing potential, we essentially follow the approach of Ref. [24]. For ρ, as argued in Ref. [24], the unquenched lattice QCD calculations rule out a large radial scale for the mixing, and the authors adopted $ \rho=0.3 $ fm in their calcultations; here, we follow their arguments and adopt the same value for ρ. However, to account for the possible effects due to a different $ c\bar c $ potential used in this work, we set Δ as a free parameter and obtain its value by fitting the mass of $ \chi_{c1}(3872) $. All the values of the parameters are then determined and collected in Table 3.

Now, we are ready to solve the coupled-channel Schrödinger equation of the system by numerical methods and obtain the spectrum of the heavy meson system. In this work, we use the renormalized Numerov algorithm to perform the numerical calculations, and the details of this method can be found in Ref. [37].

By solving the coupled-channel Schrödinger equation (Eq. (12)), the mass spectrum of the charmonium(-like) mesons and probabilities of each component in them are obtained. The results are presented in Tables 4 and 5, respectively.

$ J^{PC} $	name	Exp.	Ref. [36]	Ref. [24]	Ref. [28]	This work
$1^{--}$	$ J/\psi $	3096.9	3097	3082.4	3082.4	3097.2
	$ \psi(2S) $	3686.1	3679	3664.2	3658.8	3669.9
	$ \psi(3770) $	3778.1	3787	3790.2	3785.8	3782.8
	$ \psi(4040) $	4039.6	4078	4071		4060.1
$ 2^{++} $	$ \chi_{c2}(1P) $	3556.2	3552	3509.6	3508.7	3550.4
	$ \chi_{c2}(3930) $	3922.5	3967	3933.5	3909.0	3934.0
	$ \chi_{c2}(3P) $		4310		4006.6	4012.8
$ 1^{++} $	$ \chi_{c1}(1P) $	3510.7	3516	3510.0	3509.8	3515.0
	$ \chi_{c1}(3872) $	3871.6	3937	3871.7	3871.5	3871.7
$ 0^{++} $	$ \chi_{c0}(1P) $	3414.7	3415	3509.1	3508.8	3416.8
	$ \chi_{c0}(3860) $	3862				3862.0
	$ \chi_{c0}(3915) $	3922.1	3869	3920.4	3918.9
$ 0^{-+} $	$ \eta_c(1S) $	2984.1	2983			2984.6
	$ \eta_c(2S) $	3637.7	3635			3633.7
	$ \eta_c(3S) $		4048			4039.4
$ 1^{+-} $	$ h_c(1P) $	3525.4	3522			3521.5
	$ h_c(2P) $		3940			3927.9

Table 4.  Calculated masses for charmonium(-like) states. The experimental results (Exp.) are taken from the PDG book [35] in units of MeV. The results from a conventional quark model [36] and the results using the diabatic approaches in Refs. [24, 28] are also presented for comparison.

$ J^{PC} $	Mass/MeV	$ c\overline{c} $	$D\bar{D}$	$D\bar{D}^*$	$D_s\bar{D}_s$	$D^*\bar{D}^*$	$D_s\bar{D}^*_s$	$D^*_s\bar{D}^*_s$
$1^{--}$	3097.2	(100, 0)%
	3669.9	(95, 1)%	4%	1%
	3782.8	(0, 98)%		2%
	4060.1	(74, 5)%					20%	2%
$ 2^{++} $	3550.4	(100, 0)%
	3934.0	(69, 6)%			13%	11%	1%
	4012.8	(12, 36)%				51%
$ 1^{++} $	3515.0	100%
	3871.7	6%		94%
$ 0^{++} $	3416.8	100%
	3862.0	95%			4%	1%
$ 0^{-+} $	2984.6	100%
	3633.7	100%
	4039.4	92%					7%	1%
$ 1^{+-} $	3521.5	100%
	3927.9	93%				4%	2%

Table 5.  Masses and probabilities of various components for the charmonium(-like) states with different $ J^{PC} $. The $ c\overline{c} $ probabilities from different values of $ l_{c\bar c} $(see Table 2) are presented in parentheses. A missing entry under a meson-meson configuration means that the corresponding component gives negligible (probability inferior to 1%) or no contribution at all to the state.

In Table 4, we show the calculated results of charmonium mass spectrum together with the corresponding experimentally observed states taken from the Particle Data Group (PDG) book [35]. Furthermore, the results of three other theoretical works are also presented for comparison. The results of Ref. [36] are based on a nonrelativistic quark model without considering the coupled-channel effects from the meson-meson pairs. The other two are taken from works using the same diabatic approach but without considering the spin-dependent interactions between quarks [24, 28].

From Table 4, one can see that our calculated masses can describe the experimental data relatively well. Let us first discuss the lowest lying state of each $ J^{PC} $ quantum number. It can be seen that our results are basically the same as those of Ref. [36]. This is expected because we use the same $ c\bar c $ potential and the coupled-channel effects are very small due to their masses being far away from the open charm threshold. However, when compared to the results obtained by the diabatic approach [24, 28], we find that the results are improved significantly after including the spin-dependent interactions.

Next, let us consider higher lying states. For the $ 1^{++} $ state $ \chi_{c1}(3872) $, it should be noted that, as mentioned above, we utilize this state to determine the parameter Δ in the mixing potential. From Table 5, one can see that the probability of finding $ D\overline{D}^* $ in this state is approximately $ 94\% $. Consequently, we can interpret this state as a predominantly molecular state. Its radial wave functions of various components are plotted in Fig. 1. The root mean square (rms) radius can be obtained as 6.55 fm, which is notably large compared to its $ c\overline{c} $ component with an rms radius of 1.15 fm.

Figure 1.  (color online) Radial wave function of the calculated $ 0^+(1^{++}) $ state with a mass of 3871.7 MeV.

For the $ 0^{++} $ state $ \chi_{c0}(3860) $ [30, 38], our calculated mass is 3862 MeV, which is exactly the same as the central value suggested by the PDG. Its radial wave functions of various components are plotted in Fig. 2. The resulting rms radius is 1.05 fm. The probabilities of various compenents in this state indicate that $ \chi_{c0}(3860) $ is predominantly composed of $ c\overline{c} $ (95%). Thus, it is not surprising that our calculated mass for $ \chi_{c0}(3860) $ is similar to the quark model result [36]. This also means that the conclusion on the $ \chi_{c0}(2P) $ state still depends on the quark model employed in our work. According to the PDG [35], there are two mesons appearing at around 3900 MeV, i.e., $ \chi_{c0}(3860) $ and $ \chi_{c0}(3915) $. Meanwhile, according to Refs. [39−41], it remains unclear whether $ \chi_{c0}(3915) $ should be assigned to the $ 0^{++} $ or $ 2^{++} $ state. Our results suggest that $ \chi_{c0}(3860) $ can take the place of the $ \chi_{c0}(2P) $ state. At the same time, as we do not see another state with $ J^{PC}=0^{++} $ and a mass around 3915 MeV, this may be considered as the support for the argument that $ \chi_{c0}(3915) $ is not a $ 0^{++} $ state. Meanwhile, it is worth mentioning that recent LHCb results [42] indicate the simultaneous existence of the $ \chi_{c0}(3930) $ and $ \chi_{c2}(3930) $ states. Therefore, for a thorough understanding the $ \chi_{c0} $ states around 3900 MeV, further experimental and theoretical efforts are required.

Figure 2.  (color online) Radial wave function of the calculated $ 0^+(0^{++}) $ state with a mass of 3862.0 MeV.

For the first excited $ 2^{++} $ state, as indicated in Table 5, both $D^*\bar{D}^*$ and $D_s\bar{D}_s$ contribute approximately $ 10\% $ in this state. Experimentally, the mass of $ \chi_{c2}(3930) $ is $ 3922.5\pm 1.0 $ MeV. In comparison to the pure $ c\bar c $ results, the inclusion of the coupled-channel effects significantly improves the results. The radial wave functions of its components are plotted in Fig. 3, and an rms radius of 1.35 fm can be obtained.

Figure 3.  (color online) Radial wave function of the calculated $ 0^+(2^{++}) $ state with a mass of 3934.0 MeV.

Regarding the second excited state with a mass of 4012.8 MeV for the $ 2^{++} $ spectrum, no corresponding meson state has been observed experimentally. In our model, this is a typical unconventional state and is composed of the $ D^*\bar D^* $(51%) and $ c\overline{c} $(48%) components. Its radial wave functions are shown in Fig. 4. The corresponding rms radius is 2.33 fm. Compared to the results obtained using the quark model in [36], the threshold effects of the $ D^*\bar D^* $ channel reduce its mass by approximately 300 MeV. Consequently, the experimental search for this state can offer further verification of the approach used in this work.

Figure 4.  (color online) Radial wave function of the calculated $ 0^+(2^{++}) $ state with a mass of 4012.8 MeV.

For the states of the $ 1^{--} $ spectrum, our results are generally closer to the experimental values than those obtained without considering spin-dependent interactions. For the first excited state, the $ c\overline{c} $ component plays a dominant role (96%), while the role of the $ D\overline{D} $ component is minor (4%). The resulting mass is 3669.9 MeV, which is smaller than that of the conventional $ c\bar c $ model, and the experimental value for this state ($ \psi(2S) $) is 3686.097 $ \pm $0.010 MeV. For the second excited state, the $ c\overline{c} $ component remains dominant and is almost a pure D-wave state (98%), which is similar to the results of the quark model [36]. The computed mass for this state is 3782.8 MeV, which is in good agreement with the experimental value for the meson $ \psi(3770) $.

For the third excited state in the $ 1^{--} $ spectrum, the calculated mass is 4060.1 MeV, and the contribution of the $D_s\bar{D}^*_s$ component in this state is significant (20%). Its mass is close to the experimental value for the meson $ \psi(4040) $, with its mass being $ 4039.6\pm 4.3 $ MeV. The radial wave functions of this state are plotted in Fig. 5, and its rms radius is obtained as 1.42 fm.

Figure 5.  (color online) Radial wave function of the calculated $0^-(1^{--})$ state with a mass of 4060.1 MeV.

The low lying $ \eta_c(nS)(0^{-+}) $ states with $ n=1,2 $ and the state $ h_c(1P)(1^{+-}) $ are well described by the conventional quark model. Our results are consistent with the quark model [36] and experimental results. Also, similar to Ref. [36], we have predicted two new states whose masses are 4039.4 MeV with $ J^{PC}=0^{-+} $ and 3927.9 MeV with $ J^{PC}=1^{+-} $. The meson-meson threshold effect has a small but non-negligible impact (<10%) on the results, which makes the predicted masses of these two mesons slightly smaller than the predictions of the conventional quark model [36]. For these states, further experimental information is still needed to verify the predictions.

We end with some discussions on the approximations adopted in this work. It should be noted that, following Ref. [24], our model can only treat bound state problems, and at the energies above certain meson pair thresholds, the mixing with these meson pairs is simply ignored. Due to this simplification, in the present model, we cannot study the coupling or mixing of $ c\bar c $ states with meson pairs at energies above their mass thresholds. We remind the readers that possible states originating from such dynamics are not considered in this work.

It is also worth mentioning that, in the present work, we use the same coupling strength Δ for both the charmed meson and charmed-strange meson channels in the mixing potential. To estimate the uncertainties due to this approximation, we tried to further consider a factor of ${m_{q}}/{m_s}$ with $ m_q=0.33 $ GeV and $ m_s=0.5 $ GeV, which was used in the $ ^3P_0 $ model to account for the difference of the couplings with different flavors, for the coupling strength Δ of the charmed-strange meson channels. The calculated results show that including this factor may cause significant effects on the states lying very close to the meson pair threshold. For example, the $ 2^{++} $ state with a mass of 3934 MeV is very sensitive to this change because it lies very close to the $ D_s^+D_s^- $ threshold. In fact, when using the new coupling strength, this bound state does not exist any more. Meanwhile, for other states, this change only causes moderate effects on the results.

In summary, by employing the diabatic approach and solving the coupled-channel Schrödinger equations while considering spin-dependent interactions in the $ c\bar c $ system, we obtained the masses of the charmonium spectrum below 4.1 GeV. We also studied the probabilities of finding various components, i.e., $ c\bar c $ or meson-meson pairs, in those states. Compared to previous studies utilizing the diabatic approach [24, 28], our results can describe the charmonium spectrum better after including spin-dependent interactions. From our results, $ \chi_{c0}(3860) $ and $ \psi(3770) $ can be considered as candidates for the $ \chi_{c0}(2P) $ and $ \psi(1D) $ states, respectively. In addition, our calculations show that $ \psi(4040) $ and $ \chi_{c2}(3930) $ may have significant molecular components. It is also interesting to note that, in our work, we only get the $ \chi_{c0}(3860) $ state, which differs from the $ \chi_{c0}(3915) $ state reported in Refs. [24, 28]. As there are still disputes on the nature of $ \chi_{c0}(3915) $ [39−42], further experimental and theoretical efforts are still needed to understand its nature.