# Annihilation diagram contribution to charmonium masses

Renqiang Zhang,1,2 Wei Sun,1, Feiyu Chen,1,2 Ying<EMAIL_ADDRESS>Ming Gong,1,2 Xiangyu Jiang 1,2 and Zhaofeng Liu1,2 1 Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049, P.R. China
2 School of Physics, University of Chinese Academy of Sciences, Beijing
100049, P.R. China

###### Abstract

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 $\Delta
M_{\mathrm{hyp}}=M_{J/\psi}-M_{\eta_{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 $\Delta M_{\mathrm{hyp}}$ from both the quenched
approximation and full-QCD simulation with dynamical quarks. A recent full-QCD
calculation gives $\Delta M_{\mathrm{hyp}}=116.2(3.6)$ MeV after the continuum
extrapolation DeTar  _et al._ (2019), which is consistent with the
experimental value $\Delta M_{\mathrm{hyp}}^{\mathrm{exp.}}=113.0(0.5)$ MeV
Zyla  _et al._ (2020). The latest calculation carried out by the HPQCD
collaboration finds $\Delta M_{\mathrm{hyp}}=120.3(1.1)$ MeV at the physical
point after considering the quenched QED effects Hatton  _et al._ (2020).
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 $\Delta M_{\mathrm{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) Okubo (1963); Mandula
_et al._ (1970); ’t Hooft (1976), its contribution to $\Delta
M_{\mathrm{hyp}}$ can be sizable. To be specific, the contribution of
disconnected diagrams to the $J/\psi$ correlation function is of order
$O(\alpha_{s}^{3})$ based on a naive power counting of the strong coupling
constant $\alpha_{s}$, while it is of order $O(\alpha_{s}^{2})$ for the case
of $\eta_{c}$. Especially, since $\eta_{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 $\eta_{c}$ upward in the similar sense of the origin of
$\eta^{\prime}$ mass Schäfer and Shuryak (1998). Furthermore, given the
lattice prediction of the pseudoscalar glueball mass $M_{G}\sim 2.56$ GeV
Morningstar and Peardon (1999); Chen  _et al._ (2006); Sun  _et al._ (2018),
which is close to the $\eta_{c}$ mass, the $\eta_{c}$-glueball mixing may
introduce an additional mass shift of $M_{\eta_{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 $\Delta
M_{\mathrm{hyp}}$.

There are also a few lattice studies on the annihilation diagram contribution
to the masses of $J/\psi$ and $\eta_{c}$ McNeile and Michael (2004); de
Forcrand  _et al._ (2004). It is found the mass of $J/\psi$ is affected little
by this effect, but the $\eta_{c}$ mass can be sizably changed. Thus it is
very possible that the annihilation diagram correction to $\Delta
M_{\mathrm{hyp}}$ comes mainly from the shift of $M_{\eta_{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. Levkova and DeTar
(2011), 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_{\eta_{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
Peardon  _et al._ (2009). 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 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 Lucha  _et al._ (1991) 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

###  II.1 Lattice setup

The gauge configuration ensembles are generated on anisotropic lattices. The
gluonic action is chosen to be the tadpole improved version Morningstar and
Peardon (1997); Chen  _et al._ (2006); Sun  _et al._ (2018); Edwards  _et al._
(2008)

$\displaystyle S_{g}$ $\displaystyle=$
$\displaystyle-\beta\sum\limits_{s<s^{\prime}}\left[\frac{5}{9}\frac{\mathrm{Tr}P_{ss^{\prime}}}{\gamma_{g}u_{s}^{4}}-\frac{1}{36}\frac{\mathrm{Tr}R_{ss^{\prime}}}{\gamma_{g}u_{s}^{6}}-\frac{1}{36}\frac{\mathrm{Tr}R_{s^{\prime}s}}{\gamma_{g}u_{s}^{6}}\right]$
(1)
$\displaystyle-\beta\sum\limits_{i}\left[\frac{4}{9}\frac{\gamma_{g}\mathrm{Tr}P_{st}}{u_{s}^{2}}-\frac{1}{36}\frac{\gamma_{g}\mathrm{Tr}R_{st}}{u_{s}^{4}}\right]$

where $P_{\mu\nu}$ and $R_{\mu\nu}$ are the $1\times 1$ plaquette variable and
the $2\times 1$ rectangular Wilson loop in the $\mu\nu$ plane of the lattice,
respectively, with $s,s^{\prime}$ referring to the spatial direction and
$\nu=t$ referring to the temporal direction. The tadpole parameter $u_{s}$ is
defined through the spatial plaquette
$u_{s}=\langle\frac{1}{3}\mathrm{Re}\mathrm{Tr}P_{ij}\rangle^{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 $\gamma_{g}$ is the bare aspect
ratio parameter $\gamma_{g}\sim 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 Edwards  _et
al._ (2008); Lin  _et al._ (2009); Sun  _et al._ (2018) proposed by the Hadron
Spectrum Collaboration

$\displaystyle S_{f}$ $\displaystyle=$
$\displaystyle\sum\limits_{x}\bar{\psi}(x)\left[m_{0}+\gamma_{t}\hat{W}_{t}+\sum\limits_{s}\frac{1}{\gamma_{f}}\gamma_{s}\hat{W}_{s}\right.$
(2)
$\displaystyle-\frac{1}{4u_{s}^{2}}\left(\frac{\gamma_{g}}{\gamma_{f}}+\frac{1}{\xi}\right)\sum\limits_{s}\sigma_{ts}\hat{F}_{ts}$
$\displaystyle+\left.\frac{1}{u_{s}^{3}}\frac{1}{\gamma_{f}}\sum\limits_{s<s^{\prime}}\sigma_{ss^{\prime}}\hat{F}_{ss^{\prime}}\right]\psi(x)$

where $\hat{F}_{\mu\nu}=\frac{1}{4}\mathrm{Im}(P_{\mu\nu}(x))$ and the
dimensionless Wilson operator reads

$\displaystyle\hat{W}_{\mu}$ $\displaystyle=$
$\displaystyle\nabla_{\mu}-\frac{1}{2}\gamma_{\mu}\Delta_{\mu}$
$\displaystyle\nabla_{\mu}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[U_{\mu}(x)\delta_{x,x+\hat{\mu}}-U_{\mu}^{\dagger}(x-\hat{\mu})\delta_{x,x+\hat{\mu}}\right]$
$\displaystyle\Delta_{\mu}(x)$ $\displaystyle=$ $\displaystyle
U_{\mu}(x)\delta_{x,x+\hat{\mu}}+U_{\mu}^{\dagger}(x-\hat{\mu})\delta_{x,x-\hat{\mu}}-2.$
(3)

For the given $\beta$ and the bare quark mass parameter $m_{0}$, the bare
aspect ratio $\gamma_{f}$ for fermions in Eq. (2) is tuned along with
$\gamma_{g}$ to give the physical aspect ratio $\xi=5$ that is derived from
the dispersion relation of the pseudoscalar ($\eta_{c}$) meson

$\sinh^{2}\frac{\hat{E}}{2}=\sinh^{2}\frac{\hat{m}}{2}+\frac{1}{\xi^{2}}\sum\limits_{i}\sin^{2}\frac{\hat{p}_{i}}{2}$
(4)

where $\hat{E}$ and $\hat{m}$ are the energy and the mass of the hadron in
lattice units, $\hat{p}_{i}$ is the lattice momentum $\hat{p}_{i}=2\pi
n_{i}/L_{s}$.

In practice, the procedure of the tuning of parameters is very complicated,
since the parameters $\\{\beta,\gamma_{g},\gamma_{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 1.

Table 1: The parameters of the gauge ensembles used in this work. The spatial lattice spacing $a_{s}$ is set from the static potential and the Sommer’s parameter $r_{0}=0.491$ fm. We use two bare quark mass parameter $m_{0}$, which give the spin average masses $M_{\mathrm{avg}}$ of $\eta_{c}$ and $J/\psi$ close to the experimental values. Label | $\beta$ | $L_{s}^{3}\times T$ |  $a_{s}$ (fm) | $\xi$ | $m_{0}$ | $N_{\mathrm{cfg}}$ |  $M_{\mathrm{avg}}$ (GeV)
---|---|---|---|---|---|---|---
E1 | $2.8$ | $16^{3}\times 128$ | $0.1026$ | $5$ | 0.0408 | 6000 | $2.672$
E2 | $2.8$ | $16^{3}\times 128$ | $0.1026$ | $5$ | 0.0549 | 6000 | $3.006$

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}\times T=16^{3}\times 128$ anisotropic lattice with the aspect ratio
being set to be $\xi=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~{}\mathrm{fm}$. The bare charm quark mass parameter is set from
the spin average of physical $J/\psi$ and $\eta_{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).

Figure 1:  The profiles of the smearing function $\Psi(r)$ of quark fields.
The data points are the values of $\Psi(r)$ with the number of the
eigenvectors of the spatial Laplacian $N=10,30,50$ and 100. The dashed curves
are naive fits using Gaussian function forms in the interval
$r/a_{s}\in[0,6]$. We take $N=50$ in this study.

###  II.2 Distillation method

Quark-antiquark mesons can be accessed by quark bilinear operators
$\bar{\psi}(\mathbf{x},t)K_{U}(\mathbf{x,y};t)\psi^{\prime}(\mathbf{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 $\psi(x)$ can
be spatially smeared as

$\hat{\psi}_{\alpha}^{a}(\mathbf{x},t)=\Phi^{ab}(\mathbf{x},\mathbf{y})\psi_{\alpha}^{b}(\mathbf{y},t)$
(5)

where $a,b$ are color indices, $\alpha$ refers to the Dirac spin component,
and $\Phi^{ab}(\mathbf{x},\mathbf{y})$ is the smearing function at timeslice
$t$. Hadronic operators constructed out of smeared quark field $\hat{\psi}$
can dramatically reduce their coupling with much higher states. In order to
obtain $\Phi^{ab}$, we adopt the state-of-arts approach, namely, the
distillation method Peardon  _et al._ (2009), 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

$-\nabla_{\mathbf{xy}}^{2}(t)v_{\mathbf{y}}^{(n)}(t)=\lambda_{n}v_{\mathbf{x}}^{(n)}(t)$
(6)

where the spatial Laplacian is expressed explicitly as

$\displaystyle-\nabla_{\mathbf{xy}}^{2}(t)$ $\displaystyle=$ $\displaystyle
6\delta_{\mathbf{xy}}$ (7) $\displaystyle-$
$\displaystyle\sum\limits_{j=1}^{3}\left(\tilde{U}_{j}(\mathbf{x},t)\delta_{\mathbf{x+j,y}}+\tilde{U}_{j}^{\dagger}(\mathbf{x-j},t)\delta_{\mathbf{x-j,y}}\right)$

where $i,j$ are spatial indices and $\tilde{U}$ is the stout-link smeared
gauge link. After that, for a given Dirac index $\alpha$, each eigenvector
$v_{\mathbf{y}}^{(n)}(t)$ is used as the source vector to solve the linear
equation array

$M_{\beta\delta}(\mathbf{z},t^{\prime\prime};\mathbf{x},t^{\prime})S_{\delta\alpha}(\mathbf{x},t^{\prime};\mathbf{y},t)v_{\mathbf{y}}^{(n)}(t)=v_{\mathbf{z}}^{(n)}(t^{\prime\prime})\delta_{\beta\alpha}$
(8)

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 $\alpha=1,2,3,4$, the timeslices $t\in[0,T-1]$ and
the number of the eigenvectors $N$. Finally, multiplying the eigenvectors
$v^{(n^{\prime}),\dagger}_{\mathbf{x}}(t^{\prime})$ to each of the solution
vectors
$S_{\delta\alpha}(\mathbf{x},t^{\prime};\mathbf{y},t)v_{\mathbf{y}}^{(n)}(t)$
one gets for given $t^{\prime},t$ and $\delta,\alpha$ the matrix elements in
the subspace expanded by the $N$ eigenvectors,

$\tau_{\delta\alpha}^{n^{\prime}n}(t^{\prime},t)=v_{\mathbf{x}}^{(n^{\prime}),\dagger}(t^{\prime})S_{\delta\alpha}(\mathbf{x},t^{\prime};\mathbf{y},t)v_{\mathbf{y}}^{(n)}(t),$
(9)

which are usually called perambulators. Based on these perambulators, we can
define a propagating matrix

$\displaystyle P_{\alpha\beta}(\mathbf{x},t;\mathbf{y},t^{\prime})$
$\displaystyle=$
$\displaystyle\sum\limits_{n,n^{\prime}=1}^{N}v_{\mathbf{x}}^{(n)}(t)\tau_{\alpha\beta}^{nn^{\prime}}(t,t^{\prime})v_{\mathbf{y}}^{(n^{\prime}),\dagger}(t^{\prime})$
(10) $\displaystyle\equiv$
$\displaystyle\Phi(\mathbf{x},\mathbf{u};t)S_{\alpha\beta}(\mathbf{u},t;\mathbf{v},t^{\prime})\Phi(\mathbf{v},\mathbf{y};t^{\prime})$

with

$\Phi(\mathbf{x,y};t)=\sum\limits_{n=1}^{N}v_{\mathbf{x}}^{(n)}(t)v_{\mathbf{y}}^{(n),\dagger}(t)\equiv\left[V(t)V^{\dagger}(t)\right]_{\mathbf{xy}}$
(11)

being the desired smearing function, where $V$ is a $3V_{3}\times N$ matrix
with each column being one of the eigenvectors $v^{(n)}(t)$ for
$n=1,2,\ldots,N$ and the row number $3V_{3}$ referring to all the spatial
points and color indices. Obviously the propagating matrix
$P_{\alpha\beta}(\mathbf{x},t;\mathbf{y},t^{\prime})$ can be viewed as the
propagator of the smeared quark field
$\hat{\psi}_{\alpha}(\mathbf{x},t)=\Phi(\mathbf{x,y},t)\psi_{\alpha}(\mathbf{y},t)$
in the background gauge filed $\\{U_{\mu}(x)\\}$, namely,

$P_{\alpha\beta}(\mathbf{x},t;\mathbf{y},t^{\prime})=\langle\hat{\psi}_{\alpha}(\mathbf{x},t)\bar{\hat{\psi}}_{\beta}(\mathbf{y},t^{\prime})\rangle_{U}.$
(12)

Note that if $N=3V_{3}$, then the completeness of $v^{(n)}(t)$ implies that
$\Phi(\mathbf{x,y};t)=\delta_{\mathbf{xy}}\otimes I_{\mathrm{color}}$, such
that $P_{\alpha\beta}(x,y)$ is actually the all-to-all propagator of the
original fermion field $\psi$. One can define the norm of the smearing
function

$\Psi(|\mathbf{x-y}|)=\sum\limits_{t}\left[\mathrm{Tr}\Phi(\mathbf{x,y};t)\Phi(\mathbf{y,x};t)\right]^{1/2}$
(13)

to measure the degree to which the quark field is smeared. Due to the
translational and the rotational invariance, $\Psi(r)$ can be averaged over
all the coordinates $\mathbf{x,y}$ that satisfying $r=|\mathbf{x-y}|$. It is
easy to see that, if all the eigenvectors are involved in defining
$\Phi(\mathbf{x,y};t)$, then the completeness and the orthogonalization
require $\Psi(r)\propto\delta(r)$. If the number of eigenvectors involved is
truncated to be $N<3V_{3}$, then $\Psi(r)$ is smeared from a delta function to
a function falling off rapidly with respect to $r$. Figure 1 shows the
profiles of $\Psi(r)$ in data points at $N=10,30,50$ and 100, where $\Psi(r)$
damps faster when $N$ becomes bigger, and the dashed curves are naive fits
using Gaussian function forms in the interval $r/a_{s}\in[0,6]$. We take
$N=50$ in this study.

Now consider meson operators constructed out of the smeared quark field,

$\displaystyle\mathcal{O}_{A}(t)$ $\displaystyle=$
$\displaystyle\bar{\hat{\psi}}(\mathbf{x},t)K_{U}^{A}(\mathbf{x,y};t)\hat{\psi}(\mathbf{y},t),$
$\displaystyle\mathcal{O}_{B}(t)$ $\displaystyle=$
$\displaystyle\bar{\hat{\psi}}(\mathbf{x},t)K_{U}^{B}(\mathbf{x,y};t)\hat{\psi}(\mathbf{y},t),$
(14)

which have the same quantum number $J^{PC}$. The connected part of their
correlation function can be expressed as

$\displaystyle C_{AB}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-\mathrm{Tr}\left[P(\mathbf{x},t;\mathbf{y^{\prime}},t^{\prime})K_{U}^{B}(\mathbf{y^{\prime},y};t^{\prime})P(\mathbf{y},t^{\prime};\mathbf{x^{\prime}},t)K_{U}^{A}(\mathbf{x^{\prime},x};t)\right]$
(15) $\displaystyle=$
$\displaystyle-\mathrm{Tr}\left[\Phi(\mathbf{x,u};t)S(\mathbf{u},t;\mathbf{v},t^{\prime})\Phi(\mathbf{v,y^{\prime}};t^{\prime})K_{U}^{B}(\mathbf{y^{\prime},y};t^{\prime})\Phi(\mathbf{y,w};t^{\prime})S(\mathbf{w},t^{\prime};\mathbf{z},t)\Phi(\mathbf{z,x^{\prime}};t)K_{U}^{A}(\mathbf{x^{\prime},x};t)\right]$
$\displaystyle=$
$\displaystyle-\mathrm{Tr}\left[V^{\dagger}(\mathbf{x^{\prime}},t)K_{U}^{A}(\mathbf{x^{\prime}},\mathbf{x},t)V(\mathbf{x},t)\tau(t,t^{\prime})V^{\dagger}(\mathbf{y^{\prime}},t^{\prime})K_{U}^{B}(\mathbf{y^{\prime}},\mathbf{y},t^{\prime})V(\mathbf{y},t^{\prime})\tau(t^{\prime},t)\right]$
$\displaystyle\equiv$
$\displaystyle-\mathrm{Tr}\left[\phi^{A}(t)\tau(t,t^{\prime})\phi^{B}(t^{\prime})\tau(t^{\prime},t)\right]$

where

$\phi^{A,B}(t)=V^{\dagger}(\mathbf{x},t)K_{U}^{A,B}(\mathbf{x},\mathbf{y},t)V(\mathbf{y},t)$
(16)

is the characteristic kernel that reflects completely the structure and the
properties of operator $\mathcal{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 $\tau(t,t^{\prime})$ is universal and the
information of the meson operators are completely encoded in $\phi^{A,B}(t)$.
The disconnected part of the correlation function can be similarly derived in
terms of $\phi^{A,B}(t)$ and $\tau(t,t^{\prime})$ as follows

$D_{AB}(t,t^{\prime})=\mathrm{Tr}\left[\phi^{A}(t)\tau(t,t)\right]\mathrm{Tr}\left[\phi^{B}(t^{\prime})\tau(t^{\prime},t^{\prime})\right].$
(17)

Finally, the full correlation function of $\mathcal{O}^{A}(t)$ and
$\mathcal{O}^{B}(t^{\prime})$ is

$G_{AB}(t,t^{\prime})=C_{AB}(t,t^{\prime})+D_{AB}(t,t^{\prime}).$ (18)

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

$C_{ij}(t)=\frac{1}{2}C_{ij}^{I=0}(t)+\frac{1}{2}C_{ij}^{I=1}(t)\equiv
C_{ij}^{I=0}(t)$ (19)

if the flavor wave functions of $I=0,1$ are taken as
$\frac{1}{\sqrt{2}}(\bar{c}c\pm\bar{c}^{\prime}c^{\prime})$. Similarly the
disconnected part should be

$D_{ij}(t)=\frac{1}{2}D_{ij}^{I=0}(t),$ (20)

since gluons are isospin singlets and couple only with $I=0$ states.
Therefore, the full correlation function of the $I=0$ states is

$G_{ij}^{I=0}(t)=C_{ij}(t)+2D_{ij}(t).$ (21)

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 Sun  _et al._ (2018) by
redefining the correlation function as

$\tilde{G}^{0^{++}}(t)=G^{0^{++}}(t)-G^{0^{++}}(t+\delta t)$ (22)

where the vacuum constant term cancels. In practice, we choose $\delta t=5$.
The price we pay for this is a little lose of the signal-to-noise ratio.

Table 2: The gamma matrix combinations $\Gamma$ for the $J^{PC}$ channels corresponding to $1S$ and $1P$ states.  $J^{PC}$ | $0^{-+}$ | $1^{--}$ | $0^{++}$ | $1^{++}$ | $1^{+-}$ | $2^{++}$
---|---|---|---|---|---|---
$\Gamma$ | $\gamma_{5}$ | $\gamma_{i}$ | $I$ | $\gamma_{5}\gamma_{i}$ | $\epsilon_{ijk}\sigma_{jk}$ | $|\epsilon_{ijk}|\gamma_{j}\nabla_{k}$
$1S,1P$ | $\eta_{c}$ | $J/\psi$ | $\chi_{c0}$ | $\chi_{c1}$ | $h_{c}$ | $\chi_{c2}$

Figure 2:  The ratio functions $R(t)$ for the six channels
$J^{PC}=0^{-+},1^{--},(0,1,2)^{++}$ and $1^{+-}$ on gauge ensemble E1. The
shaded curves with error bands are the fit results using the linear part of
Eq. (30) for $1^{--},1^{+\pm},2^{++}$ channels and the results for $0^{\pm+}$
channels using Eq. (32). The darker regions indicate the fit time range.

##  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 $\Gamma$

$R^{\Gamma}(t)\equiv\frac{2D_{11}^{\Gamma}(t)}{C_{11}^{\Gamma}(t)}$ (23)

where the subscript ‘11’ refers to the correlation function of the
$|r|/a_{s}=0$ operator in the operator set of the $\Gamma$ channel. The ratios
$R^{\Gamma}(t)$ at the two quark masses are shown in Fig. 2 and Fig. 3 for
channels with $J^{PC}=0^{-+},1^{--},(0,1,2)^{++}$ and $1^{+-}$. It is seen the
magnitudes of $R^{\Gamma}(t)$ are of order $\mathcal{O}(10^{-3}-10^{-5})$: The
$R^{\Gamma}(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^{\Gamma}(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 ‘$\Gamma$’
temporarily in the following discussion. Usually the connected part of the
charmonium correlation functions and the full correlation functions can be
parameterized as

$\displaystyle C(t)$ $\displaystyle=$
$\displaystyle\sum\limits_{i}\frac{W_{i}}{2m_{i}}e^{-m_{i}t}$ $\displaystyle
G(t)$ $\displaystyle=$
$\displaystyle\sum\limits_{i}\frac{\tilde{W}_{i}}{2\tilde{m}_{i}}e^{-\tilde{m}_{i}t}.$
(24)

Since $G(t)=C(t)+2D(t)$, the ratio of the disconnect part to the connected
part can be expressed as

$R(t)=\frac{2D(t)}{C(t)}\equiv\frac{\sum\limits_{i}\frac{\tilde{W}_{i}}{2\tilde{m}_{i}}e^{-\tilde{m}_{i}t}}{\sum\limits_{i}\frac{W_{i}}{2m_{i}}e^{-m_{i}t}}-1.$
(25)

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_{i}^{2})$. When the
disconnected diagrams are considered, the propagator acquires a self-energy
correction $-\Sigma(p^{2})$, which is independent of the states. Thus the full
propagator contributing to each term of $G(t)$ can be expressed as

$\displaystyle\frac{\tilde{W}_{i}}{p^{2}+\tilde{m}_{i}^{2}}$
$\displaystyle\equiv$
$\displaystyle\frac{W_{i}}{p^{2}+m_{i}^{2}+\Sigma(p^{2})}$ (26)
$\displaystyle=$
$\displaystyle\frac{W_{i}}{p^{2}+m_{i}^{2}}+\frac{\sqrt{W_{i}}}{p^{2}+m_{i}^{2}}(-\Sigma(p^{2}))\frac{\sqrt{W_{i}}}{p^{2}+m_{i}^{2}}+\frac{\sqrt{W_{i}}}{p^{2}+m_{i}^{2}}(-\Sigma(p^{2}))\frac{1}{p^{2}+m_{i}^{2}}(-\Sigma(p^{2}))\frac{\sqrt{W_{i}}}{p^{2}+m_{i}^{2}}+\cdots$

Figure 3:  The ratio functions $R(t)$ for the six channels
$J^{PC}=0^{-+},1^{--},(0,1,2)^{++}$ and $1^{+-}$ on gauge ensemble E2. The
shaded curves with error bands are the fit results using the linear part of
Eq. (30) for $1^{--},1^{+\pm},2^{++}$ channels and the results for $0^{\pm+}$
channels using Eq. (32). The darker regions indicate the fit time range.

with the relations

$\displaystyle\tilde{W}_{i}$ $\displaystyle=$ $\displaystyle
W_{i}\left(1-\frac{d\Sigma(p^{2})}{dp^{2}}|_{p^{2}=-\tilde{m}_{i}^{2}}\right)\equiv
W_{i}(1+\epsilon_{i})$ $\displaystyle\tilde{m}_{i}^{2}$ $\displaystyle=$
$\displaystyle m_{i}^{2}+\Sigma(-\tilde{m}_{i}^{2}).$ (27)

Based on the OZI rule and from the observation in Fig. 2 and Fig. 3, it is
safe to assume $\epsilon_{i}\ll 1$ and $\delta m_{i}=\tilde{m}_{i}-m_{i}\ll
m_{i}$. Thus we have

$\delta
m_{i}=\frac{\Sigma(-\tilde{m}_{i}^{2})}{2m_{i}}\approx\frac{m_{1}}{m_{i}}\left[\delta
m_{1}+\Delta_{i}\left(1+\frac{\Delta_{i}}{2m_{1}}\right)\epsilon_{1}\right],$
(28)

where $\Delta_{i}=m_{i}-m_{1}$ has been defined and $\Sigma(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
$\mathcal{O}(10^{-3})$ or even smaller. If we assume $\epsilon_{i}$ is
insensitive to $p^{2}$ in the energy range of interest, we can take the
approximation $\epsilon_{i}\approx\bar{\epsilon}$ and $\delta
m_{i}\approx\frac{m_{1}}{m_{i}}\delta m_{1}$, which imply

$\tilde{\Delta}_{i}=\tilde{m_{i}}-\tilde{m}_{1}\approx\Delta_{i}\left(1-\frac{\delta
m_{1}}{m_{i}}\right)\approx\Delta_{i}.$ (29)

Thus we have

$\displaystyle R(t)$ $\displaystyle\approx$ $\displaystyle-\delta
m_{1}t+\epsilon_{1}-\frac{\delta m_{1}}{m_{1}}$ (30) $\displaystyle+$
$\displaystyle F(t)\frac{\delta
m_{1}}{m_{1}}\sum\limits_{i>1}\frac{W_{i}m_{1}\Delta_{i}}{W_{1}m_{i}^{2}}e^{-\Delta_{i}t}+\cdots$

where

$F(t)=\left[1+\sum\limits_{i>1}\frac{W_{i}m_{1}}{W_{1}m_{i}}e^{-\Delta_{i}t}\right]^{-1}.$
(31)

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, $\delta 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_{\mathrm{min}},t_{\mathrm{max}}]$ (the
values of $t$ are in the lattice unit $a_{t}$ in the context). The fitted
slopes $\delta m_{1}^{\Gamma}a_{t}$ and the $\chi^{2}/\mathrm{dof}$’s on the
two ensembles (labeled as E1 and E2) are listed in Table 3 along with the
$\chi^{2}/\mathrm{dof}$ in the fitting time window
$[t_{\mathrm{min}},t_{\mathrm{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\sim 10$. In the $1^{--}$ channel,
$R(t)$ shows up good linear behavior in the time interval $t\in[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\in[4,9]$ by less than $2\sigma$. If we include
the data points in the interval $t\in[10,14]$ to the fit, the central value
does not change much by with a mildly larger $\chi^{2}/\mathrm{dof}$. Actually
on E2 ensemble, all the data points in the interval $[4,14]$ converge to a
single line.

Obviously, the $\delta m_{1}$ of the $J/\psi$ 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) Okubo (1963);
Mandula  _et al._ (1970); ’t Hooft (1976), since the two charm quark loops are
mediated by at least three intermediate gluons. For $\chi_{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 Landau (1948);
Yang (1950). It is surprising that this kind of mass shift of $\chi_{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 Morningstar and Peardon (1999); Chen  _et al._
(2006); Sun  _et al._ (2018); Richards  _et al._ (2010); Gregory  _et al._
(2012), which is lower than the $\chi_{c2}$ mass and should result in positive
mass shift if $\chi_{c2}$ and the tensor glueball can mix.

Table 3: The mass shifts $a_{t}\delta m_{1}$ for all the six channels. The fit time ranges $[t_{\mathrm{min}},t_{\mathrm{max}}]$ and the values of the $\chi^{2}$ per degree of freedom ($\chi^{2}/\mathrm{dof}$) are also presented. | $\eta_{c}(0^{-+})$ | $J/\psi(1^{--})$ | $\chi_{c0}(0^{++})$ | $\chi_{c1}(1^{++})$ | $\chi_{c2}(2^{++})$ | $h_{c}(1^{+-})$
---|---|---|---|---|---|---
$[t_{\mathrm{min}},t_{\mathrm{max}}]$ | $[5,20]$ | $[4,9]$ | $[3,14]$ | $[7,14]$ | $[5,10]$ | $[3,13]$
$\delta m_{1}^{\Gamma}a_{t}$ (E1) | 0.00031(1) | $\sim 10^{-6}$ | 0.007(5) | 0.00014(2) | -0.00031(3) | 0.00010(1)
$\chi^{2}/\mathrm{dof}$ | 1.5 | 0.78 | 0.62 | 0.59 | 0.82 | 0.71
$[t_{\mathrm{min}},t_{\mathrm{max}}]$ | $[6,20]$ | $[4,14]$ | $[3,14]$ | $[6,14]$ | $[4,10]$ | $[4,13]$
$\delta m_{1}^{\Gamma}a_{t}$ (E2) | 0.00032(2) | $\sim 10^{-6}$ | 0.003(4) | 0.00006(2) | -0.00029(2) | 0.00012(2)
$\chi^{2}/\mathrm{dof}$ | 1.4 | 0.59 | 1.2 | 0.64 | 0.71 | 1.3

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
Morningstar and Peardon (1999); Chen  _et al._ (2006); Sun  _et al._ (2018);
Richards  _et al._ (2010); Gregory  _et al._ (2012), 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. (III)) to account for
the possible contribution from glueballs. Consequently, Eq. (25) for these two
channels may be modified to

$R^{\prime}(t)\approx R(t)+W_{g}e^{\Delta_{g}t}$ (32)

where $\Delta_{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 and 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 $\delta m_{1}$ of the
ground state pseudoscalar charmonium $\eta_{c}$ is determined to be +3.0(1)
MeV for E1 ensemble and +3.1(2) MeV for E2 ensemble, respectively. $\delta
m_{1}$ is converted into the value in physical units using the temperal
lattice spacing $a_{t}^{-1}=9.62$ GeV. This result is consistent with the
previous result in Ref. Levkova and DeTar (2011) and the result $\delta
m_{1}=3.9(9)$ MeV derived from the glueball-charmonium mixing mechanism Zhang
_et al._ (2022).

The values of the $\chi_{c0}$ mass shift are fitted to be 67(48) MeV for E1
ensemble and 29(38) MeV for E2 ensemble, respectively. The results for
$\chi_{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 $\delta m_{1}$ for $\chi_{c0}$.

Using the lattice spacing $a_{t}^{-1}\approx 9.62$ GeV, we get the mass shifts
$\delta m_{1}$ as follows:

$\displaystyle\mathrm{E1}:~{}~{}~{}$ $\displaystyle\delta
m_{1}(0^{-+})=3.0(1)~{}~{}\mathrm{MeV}$ $\displaystyle\delta
m_{1}(1^{++})=1.3(2)~{}~{}\mathrm{MeV}$ $\displaystyle\delta
m_{1}(2^{++})=-3.0(3)~{}~{}\mathrm{MeV}$ $\displaystyle\delta
m_{1}(1^{+-})=1.0(1)~{}~{}\mathrm{MeV}$ $\displaystyle\mathrm{E2}:~{}~{}~{}$
$\displaystyle\delta m_{1}(0^{-+})=3.1(2)~{}~{}\mathrm{MeV}$ (33)
$\displaystyle\delta m_{1}(1^{++})=0.6(2)~{}~{}\mathrm{MeV}$
$\displaystyle\delta m_{1}(2^{++})=-2.8(2)~{}~{}\mathrm{MeV}$
$\displaystyle\delta m_{1}(1^{+-})=1.2(2)~{}~{}\mathrm{MeV}.$

Since the fitted $\delta 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 $\delta m_{1}$ from E1 and E2 ensembles are
compatible with each other in most channels except for the $1^{++}$ channel,
where the deviation of $\delta m_{1}$ values on the two ensembles is larger
than $3\sigma$. 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^{PC}=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^{PC}$ 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.

###  IV.1 Optimized operators through variational method

As we addressed in Sec. II.2, the distillation method automatically provides
the gauge covariant smearing function $\Phi^{ab}(\mathbf{x},\mathbf{y})$ for
quark fields on each timeslice. To be specific, the smeared quark field
$\hat{c}(\mathbf{x},t)$ can be obtained by
$\hat{c}^{a}(\mathbf{x},t)=\Phi^{ab}(\mathbf{x},\mathbf{y})c^{b}(\mathbf{y},t)$,
where the duplicated spatial coordinate $\mathbf{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 $\hat{c}(x)$. We
introduce the spatially extended operators for each channel as

$\mathcal{O}(r,t)=\frac{1}{N_{r}}\sum\limits_{\mathbf{x},|\mathbf{r}|=r}\bar{\hat{c}}(\mathbf{x},t)K_{U}^{\Gamma}(\mathbf{x,x+r};t)\hat{c}(\mathbf{x+r},t)$
(34)

where $N_{r}$ is the multiplicity of $\mathbf{r}$ with $|\mathbf{r}|=r$, and

$K_{U}^{\Gamma}(\mathbf{x,x+r};t)=\Gamma\mathcal{P}e^{-ig\int_{\mathbf{x}}^{\mathbf{x+r}}\mathbf{A}\cdot
d\mathbf{r}}\equiv\Gamma\mathcal{L}(\mathbf{x,x+r};t)$ (35)

with $\Gamma$ being the specific combination of $\gamma$-matrix that gives the
right quantum number of each of the $1S$ and $1P$ state (The explicit
expressions of $\Gamma$’s are tabulated in Table 2), and
$\mathcal{L}(\mathbf{x,x+r};t)$ being the gauge transportation operator from
$(\mathbf{x},t)$ to $(\mathbf{x+r},t)$ . Obviously, $\mathcal{O}(r,t)$ is
gauge invariant. Thus one can get the kernel function

$\phi^{r}(t)=\frac{1}{N_{r}}\sum\limits_{\mathbf{x},|\mathbf{r}|=r}V^{\dagger}(\mathbf{x},t)K_{U}^{\Gamma}(\mathbf{x,x+r};t)V(\mathbf{x+r},t).$
(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

$\displaystyle G_{ij}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{T}\sum\limits_{t_{0}=0}^{T-1}\langle\mathcal{O}(r_{i},t+t_{0})\mathcal{O}(r_{j},t_{0})\rangle$
(37) $\displaystyle\approx$ $\displaystyle C_{ij}(t)$

where $C_{ij}(t)$ is the connected part

$\displaystyle C_{ij}(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{T}\sum\limits_{t_{0}=0}^{T-1}\left\langle\mathrm{Tr}\left[\phi^{r_{j}}(t+t_{0})\tau(t+t_{0},t_{0})\right.\right.$
(38)
$\displaystyle\times\left.\left.\phi^{r_{i}}(t_{0})\tau(t_{0},t+t_{0})\right]\right\rangle.$

Figure 4:  Effective mass of these correlators with different $r$ (labelled by
$R=r/a_{s}$). $\mathrm{M_{var}}$ labels the effective mass of ground state
optimized operator obtained from $3\times 3$ variational analysis.

In practice, $\mathbf{r}$ are chosen to be spatially on-axis displacements
with $|\mathbf{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_{\mathrm{eff}}(t)=\ln\frac{C(t)}{C(t+1)}.$ (39)

$a_{t}M_{\mathrm{eff}}(t)$ of these operators in the $0^{-+}$ channel are
shown in Fig. 4. The different $t$-behaviors of the effective mass plateaus
show the different coupling of the operators with different $|\mathbf{r}|$ to
different intermediate states. For a given quantum number $J^{PC}$, the
differences among different operators $O_{i}$ in an operator set
$\\{O_{i},i=1,2,\ldots\\}$ can be reflected by their couplings to the
orthogonal complete state set $\\{|n\rangle,n=1,2,\ldots\\}$ with the
normalization condition $\langle n|m\rangle=\delta_{mn}$,

$O_{i}^{\dagger}|0\rangle\equiv\sum\limits_{n}|n\rangle\langle
n|O_{i}^{\dagger}|0\rangle\equiv\sum\limits_{n}W_{in}^{*}|n\rangle$ (40)

where $W_{in}=\langle 0|O_{i}|n\rangle$ 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_{\mathrm{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 $\mathcal{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

$C_{ij}(t)v_{j}=\lambda(t,t_{0})C_{ij}(t_{0})v_{j}$ (41)

to get the eigenvector $v_{i}^{(n)}$ that corresponds to the $n$-th largest
eigenvalue $\lambda^{(n)}(t,t_{0})$. It is expected that the operator
$\mathcal{O}^{(n)}(t)=v_{i}^{(n)}\mathcal{O}(r_{i},t)$ couples most with the
$n$-th lowest state. Thus we have the correlation function $C^{(n)}(t)$ of the
optimized operator $\mathcal{O}^{(n)}$,

$C^{(n)}(t)=\langle\mathcal{O}^{(n)}(t)\mathcal{O}^{(n),\dagger}(0)\rangle\equiv
v_{i}^{(n)}v_{j}^{(n)}C_{ij}(t),$ (42)

whose effective mass function $a_{t}M_{\mathrm{eff}}^{(n)}(t)$ can be defined
similarly to Eq. (39). $a_{t}M_{\mathrm{eff}}^{(n)}(t)$ functions for $n=1,2$
for all the channels are plotted in Fig. 5 (ensemble E1) and Fig. 6 (ensemble
E2). $a_{t}M_{\mathrm{eff}}^{(1)}(t)$ and $a_{t}M_{\mathrm{eff}}^{(2)}(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)$

$\displaystyle C^{(n)}(t)$ $\displaystyle=$ $\displaystyle
W_{1}\left(e^{-M_{1}t}+e^{-M_{1}(T-t)}\right)$ (43) $\displaystyle+$
$\displaystyle W_{2}\left(e^{-M_{2}t}+e^{-M_{2}(T-t)}\right).$

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 $\chi^{2}$ per degree of freedom ($\chi^{2}/\mathrm{dof}$) are
around one and exhibit the fit quality. The fitted masses of $1P$, $2P$
charmonium states are tabulated in Table 4, in which the experimental values
are also given for comparison.

Table 4:  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. | $J^{PC}$ | $\chi_{c0}^{(^{\prime})}(0^{++})$ | $\chi_{c1}^{(^{\prime})}(1^{++})$ | $\chi_{c2}^{(^{\prime})}(2^{++})$ | $h_{c}^{(^{\prime})}(1^{+-})$ | $m_{\mathrm{COG}}$
---|---|---|---|---|---|---
E1 | $M(1P)(\mathrm{MeV})$ | 3039.1(4) | 3083.8(6) | 3109.7(6) | 3094.0(7) | 3093.2(4)
| $M(2P)(\mathrm{MeV})$ | 3410(4) | 3412(7) | 3531(8) | 3526(5) | 3478(5)
E2 | $M(1P)(\mathrm{MeV})$ | 3391.8(6) | 3433.8(9) | 3461.8(7) | 3443.8(9) | 3444.7(5)
| $M(2P)(\mathrm{MeV})$ | 3785(6) | 3792(9) | 3916(11) | 3917(6) | 3860(6)
Exp. | $M(1P)(\mathrm{MeV})$ | 3414.7(3) | 3510.7(1) | 3556.2(1) | 3525.4(1) | 3525.31(7)
| $M(2P)(\mathrm{MeV})$ | 3862(51) | 3871.7(1) | 3922(1) |  | 3898(6)

Figure 5:  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^{PC}=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.

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\neq 0$ has a
’center of gravity’ mass $M_{\mathrm{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_{\mathrm{COG}}$ is defined
as

$\displaystyle M_{\mathrm{COG}}(nP)$ $\displaystyle=$
$\displaystyle\frac{1}{9}\left[5M(n^{3}P_{2})+3M(n^{3}P_{1})+M(n^{3}P_{0})\right],$
(44) $\displaystyle=$ $\displaystyle M(n^{1}P_{1}).$

As shown in Table 4, 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_{\mathrm{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
$\Delta_{\mathrm{HFS}}(nP)=\Delta_{\mathrm{HFS}}=M_{h_{c}}(nP)-M_{\mathrm{COG}}(nP)$
in this paper. On the two ensembles (E1 and E2) in this paper, we obtain
$\Delta_{\mathrm{HFS}}(1P)$ as

$\displaystyle E1:~{}~{}\Delta_{\mathrm{HFS}}(1P)$ $\displaystyle=$
$\displaystyle~{}~{}0.8\pm 0.8~{}\mathrm{MeV}$ $\displaystyle
E2:~{}~{}\Delta_{\mathrm{HFS}}(1P)$ $\displaystyle=$ $\displaystyle-0.9\pm
1.0~{}\mathrm{MeV}.$ (45)

In other words, the relation Eq. (44) is satisfied for $1P$ states with a high
precision. For $2P$ states, we have

$\displaystyle E1:~{}~{}\Delta_{\mathrm{HFS}}(2P)$ $\displaystyle=$
$\displaystyle 49\pm 7~{}\mathrm{MeV}$ $\displaystyle
E2:~{}~{}\Delta_{\mathrm{HFS}}(2P)$ $\displaystyle=$ $\displaystyle 57\pm
9~{}\mathrm{MeV}.$ (46)

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 Zyla  _et al._ (2020). For charmonium
systems, the $1P$ states $(h_{c}(1^{1}P_{1}),\chi_{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_{\mathrm{COG}}=3525.3(1)$ MeV, which is almost
the same as the $h_{c}$ mass $M(h_{c})=3525.4(1)$ MeV. For bottomium systems,
the $1P$ and $2P$ states are below the $B\bar{B}$ threshold and have very
small widths. They are approximately stable particles and have direct
correspondence to the according states predicted by the non-relativistic quark
model. The $M_{\mathrm{COG}}(1P)=9899.9(5)$ MeV and
$M_{\mathrm{COG}}(2P)=10260.3(6)$ also reporduce the $h_{b}(1P)$ mass
$M(h_{b}(1P))=9899.3(8)$ MeV and the $h_{b}(2P)$ mass
$M(h_{b}(2P))=10259.8(1.2)$ MeV.

Figure 6:  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^{PC}=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.

If the relation described by Eq. (44) is somewhat universal, we can use it to
make a prediction of the mass of $h_{c}(2P)$. Experimentally, there are three
particles being assigned to be the candidates of the $\chi_{c0,1,2}$
charmonia, which are $\chi_{c1}^{\prime}(3872)$ (the old $X(3872)$),
$\chi_{c2}^{\prime}(3930)$ and $\chi_{c0}^{\prime}(3860)$, while there is no
experimental evidence of $h_{c}^{\prime}$ yet. The resonance parameters of
these three states are

$\displaystyle\chi_{c0}^{\prime}(3860):$ $\displaystyle
M_{R}=3860(50)\mathrm{MeV},\Gamma_{R}\sim 200~{}~{}\mathrm{MeV}$
$\displaystyle\chi_{c1}^{\prime}(3872):$ $\displaystyle
M_{R}=3871.7(2)\mathrm{MeV},\Gamma_{R}<1.2~{}~{}\mathrm{MeV}$
$\displaystyle\chi_{c2}^{\prime}(3930):$ $\displaystyle
M_{R}=3922(1)\mathrm{MeV},\Gamma_{R}=35(3)~{}~{}\mathrm{MeV}.$ (47)

Except that $\chi_{c0}^{\prime}(3860)$’s width is large, those of the other
two states are relatively small. Using Eq. (44), the mass of $h_{c}(2P)$ can
be tentatively estimated to be

$M(h_{c}(2P))=M_{\mathrm{COG}}(2P)\approx 3898(6)~{}\mathrm{MeV}.$ (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
Edwards  _et al._ (2008); Lin  _et al._ (2009). Our calculation shows that the
hyperfine splitting $\Delta_{\mathrm{HFS}}=M_{J/\psi}-M_{\eta_{c}}$ is
approximately 60 MeV for the two ensembles E1 and E2. This is obviously very
far from the experimental value $\Delta_{\mathbf{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 Liu
_et al._ (2012), where $\Delta_{\mathrm{HFS}}=80(2)$ MeV was obtained. This
discrepancy might be attributed to clover term
$\bar{\psi}\sigma_{\mu\nu}F_{\mu\nu}\psi$ 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 Liu  _et al._
(2012) that a larger coeffecient of the spatial clover term
$\sigma_{ss^{\prime}}\hat{F}_{ss^{\prime}}\propto\mathbf{\sigma}\cdot\mathbf{B}$
can give a larger $\Delta_{\mathrm{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 4), 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 Lucha  _et al._ (1991) (also in Appendix A). The
under-determined splittings in this work may imply that the coefficient of the
temporal clover term $\sigma_{ts}F_{ts}\propto\mathbf{\sigma}\cdot\mathbf{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 Liu  _et al._ (2002); Zhang
and Liu (2001) can produce a larger $\Delta_{\mathrm{HFS}}$ Yang  _et al._
(2013); Gui  _et al._ (2019); Ma  _et al._ (2019). 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 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^{PC}$, 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 $\eta_{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. Levkova and DeTar (2011). These
values are also consistent with mass shift 3.9(9) MeV derived from the
glueball-charmonium mixing mechanism Zhang  _et al._ (2022). For the scalar
charmonium $\chi_{c0}$, we cannot get a reliable result of the mass shift yet.
The contribution of the charm quark annihilation effect to the $J/\psi$ 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 $\chi_{c1}$ and $h_{c}$,
the mass shifts due to the charm annihilation effect are approximately 1 MeV,
which are smaller than that of $\eta_{c}$. This is understandable since the
contribution of two-gluon intermediate states is suppressed to some extent due
to the Landau-Yang theorem Landau (1948); Yang (1950). It is surprising that
this kind of mass shift of $\chi_{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 Morningstar and Peardon (1999); Chen  _et al._
(2006); Sun  _et al._ (2018); Richards  _et al._ (2010); Gregory  _et al._
(2012), which is lower than the $\chi_{c2}$ mass and should results in
positive mass shift if $\chi_{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 Hatton  _et al._
(2020). 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_{\mathrm{COG}}(nP)=M_{h_{c}}(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_{\mathrm{COG}}(1P)-M_{h_{c}}(1P)=0.8\pm 0.8$ MeV and $0.9\pm 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 smaller than those of the non-relativistic quark
model predictions Godfrey and Isgur (1985); Barnes  _et al._ (2005), but
consistent with experimental results, given the assignment that
$\chi^{\prime}_{c0}(3860)$, $\chi_{c1}^{\prime}(3872)$ and
$\chi^{\prime}_{c2}(3930)$ are $2P$ charmonium states.

## acknowlegement

This work is supported by the Strategic Priority Research Program of Chinese
Academy of Sciences (No. XDB34030302), the National Key Research and
Development Program of China (No. 2020YFA0406400), and the National Natural
Science Foundation of China (NNSFC) under Grants No.11935017, No.11775229,
No.12075253, No.12070131001 (CRC 110 by DFG and NNSFC). Y. Chen also
acknowledges the support of the CAS Center for Excellence in Particle Physics
(CCEPP). The Chroma software system Edwards and Joo (2005) and QUDA library
Clark  _et al._ (2010); Babich  _et al._ (2011) 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.

Table 5:  The expectation values of $\langle\mathbf{L}\cdot\mathbf{S}\rangle$ and $\langle S_{12}\rangle$ for $S=1$ and $L\neq 0$ states.  $J$ | $L-1$ | $L$ | $L+1$
---|---|---|---
$\langle\mathbf{L}\cdot\mathbf{S}\rangle_{J}$ | $-(L+1)$ | $-1$ | $L$
$\langle S_{12}\rangle_{J}$ | $-\frac{2(L+1)}{2L-1}$ | $2$ | $-\frac{2L}{2L+3}$

## Appendix

### A. ‘Center of gravity’ mass of heavy quarkonia

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}\approx\mathbf{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 Lucha  _et al._ (1991) is

$\displaystyle H$ $\displaystyle=$ $\displaystyle
2m+\frac{\mathbf{p}^{2}}{m}-\frac{1}{4m^{3}}\mathbf{p}^{4}+V(r)$ (49)
$\displaystyle+H_{\mathrm{SI}}+H_{\mathrm{LS}}+H_{\mathrm{SS}}+H_{\mathbf{T}}$

where $H_{\mathbf{SI}}$ is the spin-independent interaction term whose
explicit form is irrelevant of the discussion here and is omitted,
$H_{\mathbf{LS}}$ is the term for the spin-obital coupling, $H_{\mathrm{SS}}$
is the spin-spin interaction term, and $H_{\mathrm{T}}$ is the tensor
interaction part. If the potential $V(r)$ can be splitted into the vector-like
part $V_{\mathrm{V}}(r)$ and scalar-like part $V_{\mathrm{S}}(r)$, the
explicit expressions of $H_{\mathrm{LS}}$, $H_{\mathrm{SS}}$ and
$H_{\mathrm{T}}$ are

$\displaystyle H_{\mathrm{LS}}$ $\displaystyle=$
$\displaystyle\frac{1}{2m^{2}r}\left(3\frac{d}{dr}V_{\mathrm{V}}-\frac{d}{dr}V_{\mathrm{S}}\right)\mathbf{L}\cdot\mathbf{S},$
$\displaystyle H_{\mathrm{SS}}$ $\displaystyle=$
$\displaystyle\frac{2}{3m^{2}}\mathbf{S}_{1}\cdot\mathbf{S}_{2}\nabla^{2}V_{\mathrm{V}}(r),$
$\displaystyle H_{\mathrm{T}}$ $\displaystyle=$
$\displaystyle\frac{1}{12m^{2}}\left(\frac{1}{r}\frac{d}{dr}V_{\mathrm{V}}-\frac{d^{2}}{dr^{2}}V_{\mathrm{V}}\right)S_{12}$
(50)

with

$S_{12}\equiv
12\left(\frac{(\mathbf{S}_{1}\cdot\mathbf{r})(\mathbf{S}_{2}\cdot\mathbf{r})}{r^{2}}-\frac{1}{3}\mathbf{S}_{1}\cdot\mathbf{S}_{2}\right)$
(51)

where $\mathbf{S}_{1,2}$ are the spins of the charm quark and antiquark,
$\mathbf{S}=\mathbf{S}_{1}+\mathbf{S}_{2}$ is their total spin. For a given
state in terms of the eigenvalues $S,L,J$ of $\mathbf{S}$, the orbital
momentum $\mathbf{L}$ and the total angular momentum
$\mathbf{J}=\mathbf{L}+\mathbf{S}$, the expectation value
$\langle\mathbf{L}\cdot\mathbf{S}\rangle$ reads

$\langle\mathbf{L}\cdot\mathbf{S}\rangle=\frac{1}{2}\left[J(J+1)-L(L+1)-S(S+1)\right].$
(52)

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

$\langle
S_{12}\rangle=\frac{4}{(2L+3)(2L-1)}\left[\langle\mathbf{S}^{2}\mathbf{L}^{2}\rangle-\frac{3}{2}\langle\mathbf{L}\cdot\mathbf{S}\rangle-3\langle\mathbf{L}\cdot\mathbf{S}\rangle^{2}\right].$
(53)

As for the spin-spin interaction, if the $V_{\mathrm{V}}(r)$ is taken to be
Coulomb-like, namely, $V_{\mathrm{V}}\propto 1/r$, then
$H_{\mathrm{SS}}\propto\delta^{3}(\mathbf{r})\mathbf{S}_{1}\cdot\mathbf{S}_{2}$
is actually a contact interaction, whose expectation value vanishes for $L\neq
0$ states, since their radial wave function $\phi_{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 $\Delta M_{\mathrm{hfs}}$
between the spin triplet and spin singlet $S$ states.

The values of $\langle\mathbf{L}\cdot\mathbf{S}\rangle$ and $\langle
S_{12}\rangle$ for $S=1$ and $L\neq 0$ states are listed in Table 5. It is
easy to see that for a $n^{3}L_{J}$ super-multiplet with $J=L-1,L,L+1$, one
has

$\displaystyle\langle\mathbf{L}\cdot\mathbf{S}\rangle_{\mathrm{avg}}$
$\displaystyle\equiv$
$\displaystyle\frac{1}{3(2L+1)}\sum\limits_{J}(2J+1)\langle\mathbf{L}\cdot\mathbf{S}\rangle_{J}=0$
$\displaystyle\langle S_{12}\rangle_{\mathrm{avg}}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{3(2L+1)}\sum\limits_{J}(2J+1)\langle
S_{12}\rangle_{J}=0.$ (54)

Therefore, if one introduces a ‘center of gravity’ mass $M_{\mathrm{COG}}(nL)$
for $n^{3}L_{J}$ states,

$M_{\mathrm{COG}}\equiv\frac{1}{2(2L+1)}\sum\limits_{J}(2J+1)M_{J}(nL),$ (55)

which is the spin averaged mass weighted by the number of polarizations, one
should have

$M_{\mathrm{COG}}=M_{L}(n^{1}L_{L}),$ (56)

where $M_{L}(n^{1}L_{L})$ is the mass of the spin singlet state of $nL$
supermultiplet.

## References

  * DeTar  _et al._ (2019) C. DeTar, A. S. Kronfeld, S.-h. Lee, D. Mohler, and J. N. Simone (Fermilab Lattice, MILC), Phys. Rev. D 99, 034509 (2019), arXiv:1810.09983 [hep-lat] .
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Hatton  _et al._ (2020) D. Hatton, C. T. H. Davies, B. Galloway, J. Koponen, G. P. Lepage, and A. T. Lytle (HPQCD), Phys. Rev. D 102, 054511 (2020), arXiv:2005.01845 [hep-lat] .
  * Okubo (1963) S. Okubo, Phys. Lett. 5, 165 (1963).
  * Mandula  _et al._ (1970) J. Mandula, J. Weyers, and G. Zweig, Ann. Rev. Nucl. Part. Sci. 20, 289 (1970).
  * ’t Hooft (1976) G. ’t Hooft, Phys. Rev. Lett. 37, 8 (1976).
  * Schäfer and Shuryak (1998) T. Schäfer and E. V. Shuryak, Rev. Mod. Phys. 70, 323 (1998), arXiv:hep-ph/9610451 .
  * Morningstar and Peardon (1999) C. J. Morningstar and M. J. Peardon, Phys. Rev. D 60, 034509 (1999), arXiv:hep-lat/9901004 .
  * Chen  _et al._ (2006) Y. Chen _et al._ , Phys. Rev. D 73, 014516 (2006), arXiv:hep-lat/0510074 .
  * Sun  _et al._ (2018) W. Sun, L.-C. Gui, Y. Chen, M. Gong, C. Liu, Y.-B. Liu, Z. Liu, J.-P. Ma, and J.-B. Zhang, Chin. Phys. C 42, 093103 (2018), arXiv:1702.08174 [hep-lat] .
  * McNeile and Michael (2004) C. McNeile and C. Michael (UKQCD), Phys. Rev. D 70, 034506 (2004), arXiv:hep-lat/0402012 .
  * de Forcrand  _et al._ (2004) P. de Forcrand, M. Garcia Perez, H. Matsufuru, A. Nakamura, I. Pushkina, I.-O. Stamatescu, T. Takaishi, and T. Umeda (QCD-TARO), JHEP 08, 004 (2004), arXiv:hep-lat/0404016 .
  * Levkova and DeTar (2011) L. Levkova and C. DeTar, Phys. Rev. D 83, 074504 (2011), arXiv:1012.1837 [hep-lat] .
  * Peardon  _et al._ (2009) M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek, R. G. Edwards, B. Joo, H.-W. Lin, D. G. Richards, and K. J. Juge (Hadron Spectrum), Phys. Rev. D 80, 054506 (2009), arXiv:0905.2160 [hep-lat] .
  * Lucha  _et al._ (1991) W. Lucha, F. F. Schoberl, and D. Gromes, Phys. Rept. 200, 127 (1991).
  * Morningstar and Peardon (1997) C. J. Morningstar and M. J. Peardon, Phys. Rev. D 56, 4043 (1997), arXiv:hep-lat/9704011 .
  * Edwards  _et al._ (2008) R. G. Edwards, B. Joo, and H.-W. Lin, Phys. Rev. D 78, 054501 (2008), arXiv:0803.3960 [hep-lat] .
  * Lin  _et al._ (2009) H.-W. Lin _et al._ (Hadron Spectrum), Phys. Rev. D 79, 034502 (2009), arXiv:0810.3588 [hep-lat] .
  * Landau (1948) L. D. Landau, Dokl. Akad. Nauk SSSR 60, 207 (1948).
  * Yang (1950) C.-N. Yang, Phys. Rev. 77, 242 (1950).
  * Richards  _et al._ (2010) C. M. Richards, A. C. Irving, E. B. Gregory, and C. McNeile (UKQCD), Phys. Rev. D 82, 034501 (2010), arXiv:1005.2473 [hep-lat] .
  * Gregory  _et al._ (2012) E. Gregory, A. Irving, B. Lucini, C. McNeile, A. Rago, C. Richards, and E. Rinaldi, JHEP 10, 170 (2012), arXiv:1208.1858 [hep-lat] .
  * Zhang  _et al._ (2022) R. Zhang, W. Sun, Y. Chen, M. Gong, L.-C. Gui, and Z. Liu, Phys. Lett. B 827, 136960 (2022), arXiv:2107.12749 [hep-lat] .
  * Liu  _et al._ (2012) L. Liu, G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas, P. Vilaseca, J. J. Dudek, R. G. Edwards, B. Joo, and D. G. Richards (Hadron Spectrum), JHEP 07, 126 (2012), arXiv:1204.5425 [hep-ph] .
  * Liu  _et al._ (2002) C. Liu, J.-h. Zhang, Y. Chen, and J. P. Ma, Nucl. Phys. B 624, 360 (2002), arXiv:hep-lat/0109020 .
  * Zhang and Liu (2001) J.-h. Zhang and C. Liu, Mod. Phys. Lett. A 16, 1841 (2001), arXiv:hep-lat/0107005 .
  * Yang  _et al._ (2013) Y.-B. Yang, Y. Chen, L.-C. Gui, C. Liu, Y.-B. Liu, Z. Liu, J.-P. Ma, and J.-B. Zhang (CLQCD), Phys. Rev. D 87, 014501 (2013), arXiv:1206.2086 [hep-lat] .
  * Gui  _et al._ (2019) L.-C. Gui, J.-M. Dong, Y. Chen, and Y.-B. Yang, Phys. Rev. D 100, 054511 (2019), arXiv:1906.03666 [hep-lat] .
  * Ma  _et al._ (2019) Y. Ma, W. Sun, Y. Chen, M. Gong, and Z. Liu, (2019), arXiv:1910.09819 [hep-lat] .
  * Godfrey and Isgur (1985) S. Godfrey and N. Isgur, Phys. Rev. D 32, 189 (1985).
  * Barnes  _et al._ (2005) T. Barnes, S. Godfrey, and E. S. Swanson, Phys. Rev. D 72, 054026 (2005), arXiv:hep-ph/0505002 .
  * Edwards and Joo (2005) R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD), Nucl. Phys. B Proc. Suppl. 140, 832 (2005), arXiv:hep-lat/0409003 .
  * Clark  _et al._ (2010) M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi, Comput. Phys. Commun. 181, 1517 (2010), arXiv:0911.3191 [hep-lat] .
  * Babich  _et al._ (2011) R. Babich, M. A. Clark, B. Joo, G. Shi, R. C. Brower, and S. Gottlieb, in  _SC11 International Conference for High Performance Computing, Networking, Storage and Analysis_ (2011) arXiv:1109.2935 [hep-lat] .