Annihilation diagram contribution to charmonium masses

In this work, we generate gauge configurations with $N_f=2$ dynamical charm quarks on anisotropic lattices. The mass shift of $1S$ and $1P$ charmonia owing to the charm quark annihilation effect can be investigated directly in a manner of unitary theory. The distillation method is adopted to treat the charm quark annihilation diagrams at a very precise level. For $1S$ charmonia, the charm quark annihilation effect almost does not change the $J/\psi$ mass, but lifts the $\eta_c$ mass by approximately 3-4 MeV. For $1P$ charmonia, this effect results in positive mass shifts of approximately 1 MeV for $\chi_{c1}$ and $h_c$, but decreases the $\chi_{c2}$ mass by approximately 3 MeV. We have not obtain a reliable result for the mass shift of $\chi_{c0}$. In addition, it is observed that the spin averaged mass of the spin-triplet $1P$ charmonia is in a good agreement with the $h_c$, as expected by the non-relativistic quark model and measured by experiments.


I. INTRODUCTION
The hyperfine splitting of charmonia, namely ∆M hyp = M J/ψ − M ηc is a good quantity to test the precision of lattice QCD studies when charm quarks are involved. There have been many lattice efforts on ∆M hyp from both the quenched approximation and full-QCD simulation with dynamical quarks. A recent full-QCD calculation gives ∆M hyp = 116.2(3.6) MeV after the continuum extrapolation [1], which is consistent with the experimental value ∆M exp. hyp = 113.0(0.5) MeV [2]. The latest calculation carried out by the HPQCD collaboration finds ∆M hyp = 120.3(1.1) MeV at the physical point after considering the quenched QED effects [3]. Obviously, this result, with a much smaller error, still deviates from the experimental value about +7.3(1.2) MeV. Actually, most of the lattice QCD calculations of ∆M hyp do not consider the contribution of charm quark annihilation diagrams (or disconnected diagrams) to the charmonium correlation functions, which may be the major cause of the systematic uncertainty.
Although this kind of contribution is expected to be highly suppressed due to the Okubo-Zweig-Iizuka rule (OZI rule) [4][5][6], its contribution to ∆M hyp can be sizable. To be specific, the contribution of disconnected diagrams to the J/ψ correlation function is of order O(α 3 s ) based on a naive power counting of the strong coupling constant α s , while it is of order O(α 2 s ) for the case of η c . Especially, since η c is a flavor singlet pseudoscalar, its coupling to gluons can be enhanced due to the U A (1) anomaly of QCD, which may shift the η c upward in the similar sense of the origin of η mass [7]. Furthermore, given the lattice prediction of the pseudoscalar glueball mass M G ∼ 2.56 GeV [8][9][10], which is close to the η c mass, the η c -glueball mixing may introduce an additional mass shift of M ηc . Of course the light hadronic intermediate states also contribute in the presence of the light dynamical quarks. Therefore, the above effects should be taken into account if one wants to obtain a more precise theoretical determination of ∆M hyp .
There are also a few lattice studies on the annihilation diagram contribution to the masses of J/ψ and η c [11,12]. It is found the mass of J/ψ is affected little by this effect, but the η c mass can be sizably changed. Thus it is very possible that the annihilation diagram correction to ∆M hyp comes mainly from the shift of M ηc . However, no quantitative results have been obtained due to the large statistical errors yet. A more sophisticated lattice study can be found in Ref. [13], where the disconnected part of the pseudoscalar correlation function D(r) with respect to the four-dimensional distance r is complicatedly parameterized by intermediate states and quite a few excited states. The major results are that the disconnected diagrams result in an approximate +2 MeV shift of M ηc in the quenched approximation and an approximate +8 MeV shift in the N f = 2 + 1 full QCD case. However, for the lack of sea charm quarks, the mass shift of charmionium states stemming from the disconnected diagrams can only be derived indirectly and the complicated parameterization of D(r) may introduce theoretical uncertainties.
In this work, we investigate directly the contribution of charm sea quark loops to the masses of charmonium states. To do so, gauge configurations with charm sea quarks are necessary, such that the theory is in a unitary manner for charm quarks. As an exploratory study, we generate gauge ensembles on anisotropic lattices with N f = 2 degenerate sea quarks whose mass parameters are tuned to be close to the charm quark mass (The effect of light u,d and s sea quarks is ignored temporarily for simplicity). In the presence of charm sea quarks, the full correlation functions of charmonia (including the connected and disconnected parts) have well-defined spectral expressions, from which the masses can be derived directly. The key task is to perform a precise calculation on the disconnected diagrams involved. Technically, we adopt the distillation method to treat the charm quark annihilation effects [14]. On the other hand, a large statistics is mandatory to get the precise signals of disconnected diagrams. Therefore, for the computational expense to be affordable, we generate arXiv:2110.01755v2 [hep-lat] 7 Mar 2022 large gauge ensembles on anisotropic lattices with the spatial lattice spacing a s being larger than the temporal lattice spacing a t . Based on this numerical prescription, we would like to investigate directly the effect of charm quark loops on the masses of 1S and 1P charmonium states.
As we will have a large statistics for charmonium correlation functions, we can test the 'Center of Gravity' relation for both 1P and 2P states and try to give an estimate of the mass of h c (2P ). In the quark model picture, the 'Center of Gravity' mass is expected to be equal to the mass of the spin singlet [15] for charmonium states. With a large statistics, it may be possible to test this relation precisely on the lattice. This paper is organized as follows: In Sec. II we will give a brief introduction to the lattice setups for the calculation and the distillation method. Section III is devoted to the lattice derivation of mass shifts of 1S and 1P charmonium states owing to the disconnected diagrams. As a byproduct, In Sec. IV we check the agreement (or deviation) between the 'center of gravity' mass and the mass of the spin-singlet state for 1P and 2P charmonium from the connected part of charmonium correlation functions. The summary can be found in Sec. VI.

II. LATTICE SETUP AND THE DISTILLATION METHOD
A. Lattice setup The gauge configuration ensembles are generated on anisotropic lattices. The gluonic action is chosen to be the tadpole improved version [9,10,16,17] S g = −β where P µν and R µν are the 1 × 1 plaquette variable and the 2 × 1 rectangular Wilson loop in the µν plane of the lattice, respectively, with s, s referring to the spatial direction and ν = t referring to the temporal direction. The tadpole parameter u s is defined through the spatial plaquette u s = 1 3 ReTrP ij 1/4 , and is tuned self-consistently (The tadpole parameter of the temporal gauge links is set to be u t = 1 as usual). The parameter γ g is the bare aspect ratio parameter γ g ∼ a s /a t with a s and a t being the lattice spacing in the spatial and temporal direction, respectively. For fermions, say, charm quarks in this work, the action is chosen to be the anisotropic version of tadpole improved tree-level clover action [10,17,18] proposed by the Hadron Spectrum whereF µν = 1 4 Im(P µν (x)) and the dimensionless Wilson operator readŝ (3) For the given β and the bare quark mass parameter m 0 , the bare aspect ratio γ f for fermions in Eq. (2) is tuned along with γ g to give the physical aspect ratio ξ = 5 that is derived from the dispersion relation of the pseudoscalar (η c ) meson whereÊ andm are the energy and the mass of the hadron in lattice units,p i is the lattice momentump i = 2πn i /L s . In practice, the procedure of the tuning of parameters is very complicated, since the parameters {β, γ g , γ f , u s , m 0 } are highly correlated with each other. We omit the details here and just present the final ensemble parameters in Table I. As an exploratory study, we ignore the effect of light quarks and generate gauge configurations with N f = 2 flavors of degenerate charm sea quarks on an L 3 × T = 16 3 × 128 anisotropic lattice with the aspect ratio being set to be ξ = a s /a t = 5. The lattice spacing a s is determined to be a s = 0.1026 fm through the static potential and r 0 = 0.491 fm. The bare charm quark mass parameter is set from the spin average of physical J/ψ and η c masses. This lattice setup is based on two considerations: First, it is known that the signals of gluonic operators are very noisy and demand a fairly large statistics. Secondly, we will use the distillation method to tackle with the quark annihilation diagrams (see below).

B. Distillation method
Quark-antiquark mesons can be accessed by quark bilinear operatorsψ(x, t)K U (x, y; t)ψ (y, t) on the lattice. Theoretically, these operators can couple with all the meson states and possible multi-hadron systems which have the same quantum number. If one is interested in the lowest-lying states, the quark field ψ(x) can be spatially smeared aŝ where a, b are color indices, α refers to the Dirac spin component, and Φ ab (x, y) is the smearing function at timeslice t. Hadronic operators constructed out of smeared quark fieldψ can dramatically reduce their coupling with much higher states. In order to obtain Φ ab , we adopt the state-of-arts approach, namely, the distillation method [14], which is briefed as follows.
The distillation method starts with solving the eigenvalue problem of the gauge covariant second order spatial Laplacian operator on each timeslice of each gauge configuration where the spatial Laplacian is expressed explicitly as where i, j are spatial indices andŨ is the stout-link smeared gauge link. After that, for a given Dirac index α, each eigenvector v (n) y (t) is used as the source vector to solve the linear equation array where the matrix M is the Dirac matrix of fermions on the lattice and S = M −1 is the corresponding fermion propagator. This procedure is repeated for all the Dirac indices α = 1, 2, 3, 4, the timeslices t ∈ [0, T − 1] and the number of the eigenvectors N . Finally, multiplying the eigenvectors v (n ), † x (t ) to each of the solution vectors S δα (x, t ; y, t)v (n) y (t) one gets for given t , t and δ, α the matrix elements in the subspace expanded by the N eigenvectors, which are usually called perambulators. Based on these perambulators, we can define a propagating matrix being the desired smearing function, where V is a 3V 3 ×N matrix with each column being one of the eigenvectors v (n) (t) for n = 1, 2, . . . , N and the row number 3V 3 referring to all the spatial points and color indices. Obviously the propagating matrix P αβ (x, t; y, t ) can be viewed as the propagator of the smeared quark field ψ α (x, t) = Φ(x, y, t)ψ α (y, t) in the background gauge filed {U µ (x)}, namely, Note that if N = 3V 3 , then the completeness of v (n) (t) implies that Φ(x, y; t) = δ xy ⊗ I color , such that P αβ (x, y) is actually the all-to-all propagator of the original fermion field ψ. One can define the norm of the smearing function to measure the degree to which the quark field is smeared. Due to the translational and the rotational invariance, Ψ(r) can be averaged over all the coordinates x, y that satisfying r = |x − y|. It is easy to see that, if all the eigenvectors are involved in defining Φ(x, y; t), then the completeness and the orthogonalization require Ψ(r) ∝ δ(r). If the number of eigenvectors involved is truncated to be N < 3V 3 , then Ψ(r) is smeared from a delta function to a function falling off rapidly with respect to r. Figure 1 shows the profiles of Ψ(r) in data points at N = 10, 30, 50 and 100, where Ψ(r) damps faster when N becomes bigger, and the dashed curves are naive fits using Gaussian function forms in the interval r/a s ∈ [0, 6]. We take N = 50 in this study. Now consider meson operators constructed out of the smeared quark field, which have the same quantum number J P C . The connected part of their correlation function can be expressed as is the characteristic kernel that reflects completely the structure and the properties of operator O A,B (t). The last equation of Eq. (15) is the generic expression for mesonic two-point functions in the distillation method, where the perambulator τ (t, t ) is universal and the information of the meson operators are completely encoded in φ A,B (t). The disconnected part of the correlation function can be similarly derived in terms of φ A,B (t) and τ (t, t ) as follows The equation (18) is the basic expression we apply to calculate the corresponding correlation functions of 1S and 1P charmonia. There are subtleties in considering the contribution from the annihilation diagrams to the charmonium masses since there are two degenerate sea quark flavors in this work. This means there is a SU (2) global flavor symmetry, named as the isospin symmetry in principle. Thus the connected part of the correlation functions are actually if the flavor wave functions of I = 0, 1 are taken as . Similarly the disconnected part should be since gluons are isospin singlets and couple only with I = 0 states. Therefore, the full correlation function of the I = 0 states is A special attention should be paid to the 0 ++ channel. The quark loops in this channel have nonzero vacuum expectation values and should be subtracted. For this we adopt a 'time-shift' subtraction scheme [10] by redefining the correlation function as where the vacuum constant term cancels. In practice, we choose δt = 5. The price we pay for this is a little lose of the signal-to-noise ratio.

III. ANNIHILATION DIAGRAM CONTRIBUTION IN DIFFERENT CHANNELS
As the first step, we compare the relative magnitudes of the contribution of the annihilation diagram to the connected counterpart in each channel, which is measured by the ratio for a specific channel labeled by Γ where the subscript '11' refers to the correlation function of the |r|/a s = 0 operator in the operator set of the Γ channel. The ratios R Γ (t) at the two quark masses are shown in Fig. 2 and Fig. 3 for channels with The R Γ (t)'s of 0 −+ , 0 ++ and 2 ++ channels are much larger than those of the other three channels. This is understandable since the charm quark loops in 0 −+ , 0 ++ and 2 ++ channels are mediated by two gluons to the lowest order of QCD, while the charm quark loops in other three channels are mediated by at least three-gluon intermediate states.
On the other hand, the disconnected part of 0 −+ can be also enhanced by the QCD U A (1) anomaly. Especially, the R Γ (t) of 1 −− is two order of magnitude smaller and turns out to be approximately independent of t.
The contribution of the disconnected diagrams to charmonium masses can be estimated as follows. For convenience, we drop the superscript 'Γ' temporarily in the following discussion. Usually the connected part of the charmonium correlation functions and the full correlation functions can be parameterized as Since G(t) = C(t) + 2D(t), the ratio of the disconnect part to the connected part can be expressed as Intuitively, each mass term in the correlation function C(t) is given by the propagator in the momentum space, namely, W i /(p 2 − m 2 i ). When the disconnected diagrams are considered, the propagator acquires a self-energy correction −Σ(p 2 ), which is independent of the states. Thus the full propagator contributing to each term of G(t) can be expressed as with the relations Based on the OZI rule and from the observation in Fig. 2 and Fig. 3, it is safe to assume i 1 and δm i =m i − where ∆ i = m i − m 1 has been defined and Σ(p 2 ) is assumed to be varying very slowly with respect to p 2 . The latter is obviously justified by the observation that R(t) is usually of order O(10 −3 ) or even smaller. If we assume i is insensitive to p 2 in the energy range of interest, we can take the approximation i ≈¯ and δm i ≈ m1 mi δm 1 , which implỹ Thus we have where The last term in Eq. (30) will die out when t increases. Therefore, if R(t) shows up a linear dependence on t in a time region, the contribution of the disconnected diagrams to the ground state mass, namely, δm 1 , can be extracted by the slope of R(t) in this region. From Fig. 2 and Fig. 3, one can see that, for the (1, 2) ++ , 1 −− and 1 +− channels, R(t) do show up linear behaviors in different time regions. Therefore, in these time windows, we can drop the last term in Eq. (30) and fit the R(t) using linear functions. The fit results are shown in these plots by curves with error bands, where the bands in darker colors illustrate the fit time window [t min , t max ] (the values of t are in the lattice unit a t in the context). The fitted slopes δm Γ 1 a t and the χ 2 /dof's on the two ensembles (labeled as E1 and E2) are listed in Table III along with the χ 2 /dof in the fitting time window [t min , t max ] given in the table. It should be noted that the large errors of the ratio function R(t) come from the disconnected diagram which becomes very noisy beyond t ∼ 10. In the 1 −− channel, R(t) shows up good linear behavior in the time interval t ∈ [4,9] (very close to a constant) on both ensemble E1 and E2, so we attribute the behavior of R(t) in the interval [9,14] to be mostly the statistical fluctuations. Actually, the data points in this interval deviate from the linear behaviors determined from the interval t ∈ [4,9] by less than 2σ. If we include the data points in the interval t ∈ [10,14] to the fit, the central value does not change much by with a mildly larger χ 2 /dof. Actually on E2 ensemble, all the data points in the interval [4,14] converge to a single line.
Obviously, the δm 1 of the J/ψ is tiny and consistent with zero within the errors. This can be understood as the result of the strong suppression based on the Okubo-Zweig-Iizuka rule (OZI rule) [4][5][6], since the two charm quark loops are mediated by at least three intermediate gluons. For χ c1 and h c , the mass shifts due to the charm annihilation effect are small but non-zero. This is understandable since the contribution of two-gluon intermediate states is suppressed to some extent due to the Landau-Yang theorem [19,20]. It is surprising that this kind of mass shift of χ c2 is relatively large and negative. The reason for this is not clear yet.
Intuitively, lattice QCD studies predict that the lowest tensor glueball mass is approximately 2.2-2.4 GeV [8-10, 21, 22], which is lower than the χ c2 mass and should result in positive mass shift if χ c2 and the tensor glueball can mix.
The situation of the scalar and the pseudoscalar channels is a little more complicated. It is seen in Fig. 2 and Fig. 3 that the ratio functions R(t) in these two channels increase rapidly when t is large. This observation signals that there may be lighter states than the ground state charmonium contributing to the correlation function G(t) in the scalar and pseudoscalar channels. Actually, lattice QCD calculations predict that the masses of the lowest scalar and the pseudoscalar glueballs are 1.5-1.7 GeV and 2.4-2.6 GeV [8-10, 21, 22], respectively. When the disconnected diagrams are considered, these glueballs can contribute as intermediate states to the correlation function G(t). Therefore, more mass terms should be added into the expression of G(t) (Eq. (24)) to account for the possible contribution from glueballs. Consequently, Eq. (25) for these two channels may be modified to where ∆ g = m 1 − m G > 0 is a parameter resembling the mass difference between the ground state glueball and the corresponding charmonium state. Of course, one-glueball or multi-glueball states (if they exist) can also appear in other channels, but Fig. 2 and Fig. 3 show that their contributions have no sizable effects and can be ignored in practice.
We tentatively use the functions form of Eq. (32) to fit R(t) of the scalar and the pseudoscalar channels, with R(t) in this equation being approximated by a linear function. The fit results are illustrated in Fig. 2 Fig. 3 by curves with error bands, where one can see that the function form describes the data very well in large time ranges. This manifests the efficacy and the necessity of the second term in Eq. (32). The mass shift δm 1 of the ground state pseudoscalar charmonium η c is determined to be +3.0(1) MeV for E1 ensemble and +3.1(2) MeV for E2 ensemble, respectively. δm 1 is converted into the value in physical units using the temperal lattice spacing a −1 t = 9.62 GeV. This result is consistent with the previous result in Ref. [13] and the result δm 1 = 3.9(9) MeV derived from the glueballcharmonium mixing mechanism [23].

and
The values of the χ c0 mass shift are fitted to be 67(48) MeV for E1 ensemble and 29(38) MeV for E2 ensemble, respectively. The results for χ c0 have large errors due to nonexistence or the shortness of the linear behavior. Due to the large error, we cannot get a reliable final result of the δm 1 for χ c0 .
Using the lattice spacing a −1 t ≈ 9.62 GeV, we get the mass shifts δm 1 as follows: Since the fitted δm 1 (1 −− ) is very small, we do not quote its value here but only say that the charm annihilation effects in 1 −− are negligible. The mass shifts δm 1 from E1 and E2 ensembles are compatible with each other in most channels except for the 1 ++ channel, where the deviation of δm 1 values on the two ensembles is larger than 3σ. The reason for this discrepancy is not clear yet and should be investigated in the future.

IV. THE CENTER OF GRAVITY MASS OF P -WAVE CHARMONIUM
Now that we have calculated the perambulators of charm quarks on the two gauge ensembles of large statistics, we would like to take a look at the charmonium spectrum. In this work, we will focus on the lowest two states in J P C = 0 −+ , 1 −− , (0, 1, 2) ++ and 1 −+ channels. In the previous sections we have shown that the inclusion of the disconnected diagrams does not change much the masses of charmonium states, so in the following discussions we only consider the connected part of the correlation functions. For each J P C quantum number, we will first build an operator set, then calculate the correlation matrix of the operators in this set. After that, we will solve the generalized eigenvalue problem (GEVP) to obtain the optimized operators that couple most to specific states. What follows are the technique details.

A. Optimized operators through variational method
As we addressed in Sec. II B, the distillation method automatically provides the gauge covariant smearing function Φ ab (x, y) for quark fields on each timeslice. To be specific, the smeared quark fieldĉ(x, t) can be obtained byĉ a (x, t) = Φ ab (x, y)c b (y, t), where the duplicated spatial coordinate y means the summation over the space volume. Thus the quark bilinear operators for meson states in this study can be built in terms of the smeared charm quark fieldĉ(x). We introduce the where N r is the multiplicity of r with |r| = r, and with Γ being the specific combination of γ-matrix that gives the right quantum number of each of the 1S and 1P state (The explicit expressions of Γ's are tabulated in Table II), and L(x, x + r; t) being the gauge transportation operator from (x, t) to (x + r, t) . Obviously, O(r, t) is gauge invariant. Thus one can get the kernel function (36) Because the disconnected part is far more smaller than connected part and has little effect on present discussion, the correlation function G ij (t) can be approixmated by their connected part as where C ij (t) is the connected part In practice, r are chosen to be spatially on-axis displacements with |r|/a s = 0, 1, 2, 3, 4, 5, 6, 7. The differences between the different operators can be monitored through the effective masses of the corresponding correlation functions C(t) of the operators a t M eff (t) = ln C(t) C(t + 1) .
where W in = 0|O i |n is defined. Therefore, the differences of W in show how different the operators in this operator set are, and also result in the different t-behaviors of the effective masses M eff (t) of the corresponding correlation functions in the early time region. Since in each channel we only focus on the lowest two states, we select operators O(r, t) with r/a s = 0, 3, 6 to compose an operator set, which are expected to couple to the intermediate states more differently, as manifested by the effective masses of these correlation functions. Based on the operator set in each channel, we first carry out the variational method analysis to obtain the optimized operators that couple most to specific states. That is to say, we calculated the correlation matrix {C ij (t)} of each operator set and solve the generalized eigenvalue problem to get the eigenvector v (n) i that corresponds to the nth largest eigenvalue λ (n) (t, t 0 ). It is expected that the operator O (n) (t) = v (n) i O(r i , t) couples most with the nth lowest state. Thus we have the correlation function C (n) (t) of the optimized operator O (n) , whose effective mass function a t M    eff (t) are clearly separated from each other, but still have mild time dependence in the early time range due to slight contamination from higher states. In order to obtain the mass values more precisely, we perform two-mass-term fits to C (n) (t) where the second mass term is adopt to account for the higher state contamination. W 1 and W 2 are coefficients of the first mass term and the second mass term respectively. The fit results with the fitted parameters are also plotted in Fig. 5 and Fig. 6 as curves with error bands, where the extensions of the curves illustrate the fit windows. For all the fits, the values of the χ 2 per degree of freedom (χ 2 /dof) are around one and exhibit the fit quality. The fitted masses of 1P , 2P charmonium states are tabulated in Table IV, in which the experimental values are also given for comparison. The non-relativistic quark model expects that for the same principal quantum number n, the n 2S+1 L J multiplet with S = 1 and L = 0 has a 'center of gravity' mass M COG , which is the spin average mass of this multiplet and should be equal to the mass of the spin singlet counterpart (the derivation of the relation can be found in the Appendix A). For example, for the L = 1 multiplets n 3 P J multiplet, M COG is defined as M COG (nP ) = 1 9 5M (n 3 P 2 ) + 3M (n 3 P 1 ) + M (n 3 P 0 ) , As shown in Table IV, the masses of the 1P states are determined very precisely with the statistic error being less than 1 MeV, so we can check the relation of Eq. (44). The difference between the M COG (nP ) and the mass of the spin singlet h c (nP ) is sometimes called the hyperfine splitting of nP charmonium states, which is denoted by ∆ HFS (nP ) = ∆ HFS = M hc (nP ) − M COG (nP ) in this paper. On the two ensembles (E1 and E2) in this paper, we obtain ∆ HFS (1P ) as In other words, the relation Eq. (44) is satisfied for 1P states with a high precision. For 2P states, we have E1 : ∆ HFS (2P ) = 49 ± 7 MeV E2 : ∆ HFS (2P ) = 57 ± 9 MeV.
There is a substantial deviation from Eq. (44). We are not sure the reason for this deviation yet. It is possible that 2P charmonium states do have a non-zero hyperfine splitting. It is also possible that the masses of 2P states are not determined precisely enough. This should be explored in future studies. It is interesting to note that the experimental results support the relation of Eq. (44) to a very high precision [2]. For charmonium systems, the 1P states (h c (1 1 P 1 ), χ c0,1,2 (1 3 P 0,1,2 )) have been well established. According to the PDG data, the 'center of gravity' mass of 1P charmonia is M COG = 3525. (48)

V. SOME CAVEATS
We would like to point out the possible caveats of our lattice setup. As we addressed in Sec. II, the tadpole improved tree level clover action proposed by the Hadron Spectrum Collaboration is adopted for charm quarks in this study [17,18]. Our calculation shows that the hyperfine splitting ∆ HFS = M J/ψ − M ηc is approximately 60 MeV for the two ensembles E1 and E2. This is obviously very far from the experimental value ∆ HFS = 113(3) MeV, and signals large discretization errors. The Hadron Spectrum Collaboration also found this large discrepancy in its calculation of the charmonium spectrum [24], where ∆ HFS = 80(2) MeV was IV. The masses of 1P and 2P states derived from ensemble E1 and E2. The masses are converted into the values in physical units. The 'center of gravity' masses are the spin averages of the spin-triplet nP charmonium masses. The experimental results are also presented for comparison.  The mass plateaus of the correlation functions of the optimized operators obtained through GEVP method on gauge ensemble E1. In each of the six channels J P C = 0 −+ , 1 −− , (0, 1, 2) ++ and 1 −+ , the mass plateaus corresponding to the lowest two states are shown. For each correlation function, a two-mass-term fit is performed with the second mass term being introduced to take into account the residual contamination of higher states. The darker colored bands illustrate the fit results using the corresponding time window.
obtained. This discrepancy might be attributed to clover termψσ µν F µν ψ of the fermion action, which results in the spin-orbital and spin-spin interactions between the charm quark and antiquark in the non-relativistic approximation. It is found [24] that a larger coeffecient of the spatial clover term σ ss F ss ∝ σ · B can give a larger ∆ HFS . On the other hand, even the masses of the 1P and 2P charmonium states are determined very precisely in this work (see in Table IV), it is obvious that the fine splittings between the 1P states are also smaller than those of the experimental values. In the non-relativistic approximation, these splittings are due to the spin-orbital and the tensor interactions [15] (also in Appendix A). The under-determined splittings in this work may imply that the coefficient of the temporal clover term σ ts F ts ∝ σ · E of the fermion action also has large discretization uncertainties. It may also because we ignored light u,d,s quarks. Based on these observations, it seems that the fermion action we adopt in this work may not be a very good anisotropic version for charm quarks. It is noted that another version of the clover fermion action on the anisotropic lattices [25,26] can produce a larger ∆ HFS [27][28][29]. Of course, whatever the lattice setup is, the continuum limit should be the same, but the smaller the discretization uncertainties, the smaller systematic errors after the continuum extrapolation.

VI. SUMMARY
In this work, we generate gauge ensembles with two flavor degenerate dynamical charm quarks on anisotropic lattices. This lattice setup enables us to investigate the contribution of charm quark annihilation diagrams to charmonium masses in an unitary theoretical framework The mass plateaus of the correlation functions of the optimized operators obtained through GEVP method on gauge ensemble E2. In each of the six channels J P C = 0 −+ , 1 −− , (0, 1, 2) ++ and 1 −+ , the mass plateaus corresponding to the lowest two states are shown. For each correlation function, a two-mass-term fit is performed with the second mass term being introduced to take into account the residual contamination of higher states. The darker colored bands illustrate the fit results using the corresponding time window.
for charm quarks. The distillation method is adopted to realize both the smearing scheme of quark fields and the disconnected diagrams in calculating meson correlation functions. With large statistics, the effects of disconnected diagrams can be derived with a high precision.
For a given quantum number J P C , the mass shift of charmonium masses due to the disconnected diagrams can be derived from the ratio R(t) of the disconnected part of the charmonium correlation function to the connected part. It is found that the R(t) are relatively large for the scalar and the pseudoscalar channel, and the time dependence of R(t) imply that there may be contribution from lighter states than the corresponding lowest charmonium states to the correlation functions when the disconnected diagrams are considered. This is not surprising since glueballs can appear as intermediate states in this case, whose masses are predicted by lattice QCD to be lower than charmonium states. By considering the contribution of glueballs, the mass shift of the ground state pseudoscalar charmonium η c is +3.0(1) MeV and +3.1(2) MeV, respectively,for the two gauge ensembles involved, which is consistent with the previous result in Ref. [13]. These values are also consistent with mass shift 3.9(9) MeV derived from the glueball-charmonium mixing mechanism [23]. For the scalar charmonium χ c0 , we cannot get a reliable result of the mass shift yet. The contribution of the charm quark annihilation effect to the J/ψ mass is very tiny as expected by the OZI rule, since the two charm quark loops are mediated by at least three intermediate gluons. For χ c1 and h c , the mass shifts due to the charm annihilation effect are approximately 1 MeV, which are smaller than that of η c . This is understandable since the contribution of two-gluon intermediate states is suppressed to some extent due to the Landau-Yang theorem [19,20]. It is surprising that this kind of mass shift of χ c2 is relatively large and negative (-3.0(3) MeV for E1 and -2.8(2) for E2). We don't have an explanation to this yet. Lattice QCD studies predict that the lowest tensor glueball mass is approximately 2.2-2.4 GeV [8-10, 21, 22], which is lower than the χ c2 mass and should results in positive mass shift if χ c2 and the tensor glueball can mix. Note that in this study we only consider the charm sea quark effects, the effects of u, d, s quarks have not been taken into consideration yet, which are expected theoretically to affect the mass spectrum of charmonia also [3]. In the presence of u, d, s quarks, the situation will be much more complicated due to charmoniums decays and appearance of light hadron intermediate states. Therefore, the study of their impact on charmonium masses is still a challenging task in the present stage of lattice QCD.
The relation M COG (nP ) = M hc (nP ) is expected by the non-relativistic quark model. We observe that this relation is satisfied to a very high precision level for 1P charmonia, namely M COG (1P ) − M hc (1P ) = 0.8 ± 0.8 MeV and 0.9 ± 1.0 MeV at the two charm quark masses used in this work. For the 2P charmonia, we observe a large deviation from this relation. There are two possible reasons for this deviation: this relation actually does not hold for 2P states, or the masses of the 2P states derived in this work have contaminations from higher states. This should be explored in future studies. On the other hand, for the two charm quark masses, the mass splittings between the 1P and 2P charmonium we obtained are TABLE V. The expectation values of L · S and S12 for S = 1 and L = 0 states.
smaller than those of the non-relativistic quark model predictions [30,31], but consistent with experimental results, given the assignment that χ c0 (3860), χ c1 (3872) and χ c2 (3930) are 2P charmonium states. Chen also acknowledges the support of the CAS Center for Excellence in Particle Physics (CCEPP). The Chroma software system [32] and QUDA library [33,34] are acknowledged. The computations were performed on the CAS Xiandao-1 computing environment, the HPC clusters at Institute of High Energy Physics (Beijing) and China Spallation Neutron Source (Dongguan), and the GPU cluster at Hunan Normal University. In the quark model picture, the spectroscopy of charmonia can be studied by solving the non-relativistic Schrödinger equation of bound states with QCD-inspired inter-quark potentials V (r) along with the relativistic corrections. To order v 2 ≈ p 2 /m 2 , where v is the velocity of the (anti-)charm quark with mass m within a bound state, the generalized Breit-Fermi Hamiltonian [15] is where H SI is the spin-independent interaction term whose explicit form is irrelevant of the discussion here and is omitted, H LS is the term for the spin-obital coupling, H SS is the spin-spin interaction term, and H T is the tensor interaction part. If the potential V (r) can be splitted into the vector-like part V V (r) and scalar-like part V S (r), the explicit expressions of H LS , H SS and H T are with S 12 ≡ 12 (S 1 · r)(S 2 · r) where S 1,2 are the spins of the charm quark and antiquark, S = S 1 + S 2 is their total spin. For a given state in terms of the eigenvalues S, L, J of S, the orbital momentum L and the total angular momentum J = L + S, the expectation value L · S reads L · S = Obviously this expectation value for L = 0 or S = 0 states vanishes. It is easy to see that the expectation value of S 12 also vanishes for L = 0 or S = 0. Otherwise the expectation value of S 12 reads (53) As for the spin-spin interaction, if the V V (r) is taken to be Coulomb-like, namely, V V ∝ 1/r, then H SS ∝ δ 3 (r)S 1 ·S 2 is actually a contact interaction, whose expectation value vanishes for L = 0 states, since their radial wave function φ nl (r) at the origin r = 0 is zero. For L = 0 states (n 1 S 0 and n 3 S 1 states), the spin-spin interaction results the hyperfine splitting ∆M hfs between the spin triplet and spin singlet S states.
The values of L · S and S 12 for S = 1 and L = 0 states are listed in Table V. It is easy to see that for a n 3 L J super-multiplet with J = L − 1, L, L + 1, one has L · S avg ≡ 1 3(2L + 1) J (2J + 1) L · S J = 0 where M L (n 1 L L ) is the mass of the spin singlet state of nL supermultiplet.