Study of the pseudoscalar glueball in J/ψ radiative decays

1,2,3 ∗ 4,5 †
Long-Cheng Gui , Jia-Mei Dong 1 , Ying Chen , and Yi-Bo Yang 6‡
1
Department of Physics , Hunan Normal University, ChangSha, 410081 , China
2
Key Laboratory of Low-Dimensional Quantum Structures and Quantum
Control of Ministry of Education, Changsha 410081, China
3
Synergetic Innovation Center for Quantum Effects and Applications(SICQEA),
Hunan Normal University, Changsha 410081,China
4
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
5
School of Physics, University of Chinese Academy of Sciences, Beijing 100049, China
6
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China

We aim to explore the production rate of the pseudoscalar glueball in J/ψ radiative decay by lattice
QCD in quenched approximation. The calculation is performed on three anisotropic lattices with the
spatial lattice spacing ranging from 0.222(2) fm to 0.110(1) fm. As a calibration of some systematical
arXiv:1906.03666v3 [hep-lat] 29 Sep 2019

uncertainties, we first extract the M 1 form factor V̂ (0) of the process J/ψ → γηc and get the result
V̂ (0) = 1.933(41) in the continuum limit, which gives the partial width Γ(J/ψ → γηc ) = 2.47(11)
keV. These results are in agreement with that of previous lattice studies. As for the pseudoscalar
glueball G0−+ , its mass is derived to be 2.395(14) GeV, and the form factor V ˆ(0) of the process
J/ψ → γG0−+ is determined to be V̂ (0) = 0.0246(43) after continuum extrapolation. Finally,
the production rate of the pseudoscalar glueball is predicted to be 2.31(90) × 10−4 , which is much
smaller than that of conventional light q q̄ η states. After the subtraction of the phase space factor,
the couplings of J/ψXγ are similar where X stands for η states and the pseudoscalar glueball.
Possibly, the UA (1) anomaly plays an important role for the large couplings of gluons to the flavor
singlet η states in J/ψ radiative decays.

PACS numbers: [Link], [Link], [Link], [Link]

I. INTRODUCTION Lattice QCD is the ab initio nonperturbative approach

for solving QCD and plays a key role in the investigation
Quantum chromodynamics (QCD) predicts the ex- of the low energy strong interaction phenomena. In
istence of glueballs, namely, the bound states of the glueball sector, lattice QCD in the quenched
gluons. Last several decades witnessed the intensive and approximation predicts that the masses of the lowest
extensive investigations of glueballs both in experiments lying scalar, tensor and pseudoscalar glueballs are
and theoretical studies[1–12]. Experimentally, there roughly 1.5 − 1.7 GeV, 2.2 − 2.4 GeV, and 2.6 GeV,
are ten scalar mesons observed with approximately respectively[1, 2]. Recent lattice calculations with
degenerated masses around 1.5 GeV, among which dynamical quarks seemingly support these predictions
f0 (1370), f0 (1500), and f0 (1710) are the three isoscalars. and do not observe large unquenched effects[5, 6]. As far
According to the quark model, if these states can be as the pseudoscalar glueball is concerned, its predicted
sorted into the q q̄ SU (3) nonet, then the surplus one mass, say, around 2.6 GeV, is much higher than that of
isoscalar hints the existence of an additional degree of η(1405). This discrepancy cannot be easily mediated by
freedom, possibly a glueball states, which can be either considering the glueball-meson mixing in the presence of
one of the isoscalars mentioned above, or mixes with dynamical quarks. It is interesting to notice that some
the conventional q q̄ states. A similar consideration phenomenological studies advocate η(1405) and η(1495)
applies to the pseudoscalar channel: the three isoscalar be the same state which appears differently in different
pseudoscalar mesons η(1295), η(1405), and η(1495) also final states due to some dynamical mechanism[13]. If
motivate the conjecture of the existence of a pseudoscalar this is the case, then there is no redundant pseudoscalar
glueball in this mass range. Actually, there are many meson mass region and subsequently no need for an
phenomenological studies assigning η(1405) to be the additional degree of freedom such as a glueball in 1.5
most likely candidate for the pseudoscalar glueball due GeV mass region. As such, one may wish to search for
to its large production fraction in the J/ψ radiative the pseudoscalar glueball in the energy range beyond 2
decays. However, this assignment has tensions with the GeV according to the lattice predictions.
prediction of the pseudoscalar glueball mass from lattice J/ψ radiative decays are regarded as an important
QCD. hunting ground for glueballs, owing to its the gluon-rich
environment and cleaner background. Apart from their
masses, the production rates of glueballs in J/ψ radiative
∗ guilongcheng@[Link] serve as additional key criteria for the identification
† cheny@[Link] of glueballs, if they can be derived reliably from the
‡ ybyang@[Link] theoretical calculation. Also some quenched lattice QCD
 2

efforts have been made to calculate these production following formula,

rates of the scalar and tensor glueballs [14, 15]. Since
the pure gauge glueballs are well defined hadron states 1 |~k| X
ΓJ/ψ→γX = 2 pi , p~f )|2 , (1)
|MλJ/ψ ,λγ (~
in the quenched approximation, the electromagnetic form 24π MJ/ψ
λJ/ψ ,λγ
factors of J/ψ radiatively decaying into glueballs can be
extracted directly by calculating the matrix elements of where p~i and p~f are the momenta of J/ψ and X,
the electromagnetic current between J/ψ and glueballs. respectively, ~k = p~i − p~f = −~ pf is the decaying
With these form factors, the branch fraction of J/ψ momentum of the emitted photon, λJ/ψ and λγ stand
radiative decaying into the pure gauge scalar glueball is for the different polarizations of J/ψ and the photon,
predicted to be 3.8(9) × 10−3 . It is interesting to notice and MλJ/ψ ,λγ is the on-shell transition amplitude. The
that the sum of the observed branching fraction of the
processes J/ψ → γf0 (1710) → γ +anything gives a value magnitude of ~k can be defined through the masses of
of roughly 2.0 × 10−3 which is very close to the above J/ψ (denoted by MJ/ψ ) and X ( denoted by MX ), say,
M 2 −M 2
predicted value for the pure gauge glueball, while that of |~k| = J/ψ X . The transition amplitude contains all
2MJ/ψ
J/ψ → γf0 (1500) is an order of magnitude smaller. the dynamics of the decay, and can be expressed to the
Recently, the BESIII collaboration has performed the lowest order of QED as
partial wave analysis to the process J/ψ → γφφ and
observed a new resonance X(2500) with the resonance MλJ/ψ ,λγ = ∗µ (~k, λγ )hX(~ µ
pf )|Jem |J/ψ(~
pi , λJ/ψ )i, (2)
parameters MX = 2470+15+101 +64+56
−19−23 and ΓX = 230−35−33 [16].
In the process J/ψ → γη 0 π + π − , BESIII observed a where ∗µ (~k, λγ ) is the polarization vector of photon
complicated structure in the η 0 π + π − invariant mass µ
ef ψ̄f γ µ ψf is the relevant electromagnetic
P
and Jem =
spectrum and reported two new resonances X(2120) f
and X(2370), whose quantum number are likely J P C = vector current with the summation over all possible
0−+ [17]. These new resonances lie in the mass range quark flavors. If the sea quark contributions through
of the pseudoscalar glueball predicted by lattice QCD disconnected diagrams can be neglected, one can use
µ
and are worthy of further experimental investigations. Jem = ec c̄γµ c as an approximation for the electromag-
However, in order to unravel the nature of these states, netic decays of charmonia. Usually, the matrix elements
µ
more theoretical inputs are desired. hX(~ pf )|Jem (0)|J/ψ(~pi , λJ/ψ )i can be expressed in terms
In this work, we calculate the production rate of of form factor through the multipole decomposition. For
the pseudoscalar glueball from lattice QCD by adopting the pseudoscalar X, the explicit expression is
the similar strategy in the scalar and tensor cases. hX(~ µ
pf )|Jem pi , λJ/ψ )i = M (Q2 )µνρσ pi,ν pf,ρ
(0)|J/ψ(~
Even though the dynamical lattice QCD simulation
is dominant nowadays, the study of glueballs is still × σ (~
pi , λJ/ψ ) (3)
challenging since it requires much higher statistics where Q2 ≡ −(pi −pf )2 and M (Q2 ) is the multipole form
comparing to the usual hadron. Therefore, we would like factor, which is sometimes also expressed in terms of a
to perform an exploratory investigation by generating dimensionless form factor V (Q2 ) as M (Q2 ) = m2V (Q2 )
.
the quenched gauge configurations with large statistics X +mJ/ψ

using the anisotropic lattices. In order to calibrate part Thereby the partial decay width can be written as
of uncertainties, we first calculate the partial width of 1 ~3 2
ΓJ/ψ→γX = |k| |M (0)|
the process J/ψ → γηc and compare with the results of 12π
previous lattice calculations and experiments. We admit 1 |~k|3 2
that glueballs may have strong mixing with conventional ≡ |V (0)| . (4)
3π (mJ/ψ + mX )2
mesons in the real world, and the quenched effects on
our result of glueball production rates cannot be reliably It is clearly seen that if the matrix elements in Eq. (3) are
estimated in the present stage. known, the form factor M (Q2 ) (or equivalently V (Q2 ))
This work is organized as follows: Section II gives can be derived to give the decay width directly. Actually,
an introduction to the formalism for calculating the this goal can be achieved in lattice QCD study by
radiative transition width of J/ψ from lattice QCD. calculating the relevant three-point correlation functions
Section III presents the calculation details and results (3),µi
X
including the parameters of the lattice, the relevant ΓJ/ψ→γX (~ pi , ti = 0;~
pf , tf ; ~q, t) = ei~pi ·~x ei~q·~y e−i~pf ·~z
spectrum, and transition form-factors. We give the ~
x,~
y ,~
z
conclusion and some discussions in Section IV. µ †
hOX (~z, tf )Jem (~y , t)OJ/ψ,i (~x, ti )i,
(5)
II. FORMALISM
where OX and OJ/ψ,i are the interpolating field operators
for X and J/ψ, respectively. After the intermediate
Generally speaking, in the rest frame of J/ψ, the (3),µi
state insertion, the three-point function ΓJ/ψ→γX is
partial decay width of J/ψ radiatively decaying into a
parametrized as
pseudoscalar meson X can be calculated through the
 3

(3),µi e−Em tf e−(En −Em )t

P
ΓJ/ψ→γX (~
pi , ti = 0; p~f , tf ; ~q, t) = 4Em En hΩ | Of (0) | m, p~f i
m,n,λn
µ † (6)
× hm, p~f | Jem (0) | n, p~i , λn ihn, p~i , λn | OJ/ψ,i (0) | Ωi
tf t0 e −EX tf −(EJ/ψ −EX )t
e ∗(J/ψ) (X) µ
−−−−−→ 4EJ/ψ EX Zi Z hX(~ pf ) | Jem (0) | J/ψ(~ pi , λJ/ψ )i

where EX and EJ/ψ are the energies of X and J/ψ,

(J/ψ)
TABLE I. Relevant input parameters in this work. The
respectively, Zi = hJ/ψ(~ pf , λJ/ψ )|OJ/ψ,i (0)|Ωi, and spatial lattice spacing as is determined from r0−1 =
(X)
Z = hX(~ pf )|OX (0)|Ωi, which can be extracted from 410(20) MeV by calculating the static potential.
the relevant two-point functions
β ξ as (fm) Las (fm) L3 × T Nconf.
(2) 2.4 5 0.222(2) 1.78 83 × 96 20000
e−i~p·~x hΩ | Oi (~x, t)Oj† (0, 0) | Ωi
X
Γij (~p, t) = 2.8 5 0.138(1) 1.66 123 × 144 20000
t→∞ Zi Zj∗ −E(~p)t 3.0 5 0.110(1) 1.76 163 × 160 10000
−−−→ e . (7)
2E(~p)
Practically, one can carry out a joint fit to the two- (s) (t)
TABLE II. The renormalization constants ZV and ZV of
point functions and the three-point function to extract
the spatial and temporal components of the vector current
the desired matrix elements in Eq. (3), from which the for β = 2.4, β = 2.8 and β = 3.0 lattices. Two momentum
multipole form factors can be derived at different Q2 . modes, (0, 0, 0) and (1, 0, 0), are used for the derivation.
(t) (t) (s)
β ZV (0, 0, 0) ZV (1, 0, 0) ZV (1, 0, 0)
III. NUMERICAL DETAILS 2.4 1.288(5) 1.299(11) 1.388(15)
2.8 1.155(3) 1.159(3) 1.110(7)
3.0 1.106(4) 1.114(6) 1.062(6)
As addressed in Sec. I, we perform the calculation
in the quenched approximation. Since J/ψ, ηc , and
the pseudoscalar glueball are heavy particles, in order
to obtain good signals with high resolutions in the Following the scheme proposed by Ref. [21], ZV (as ) is
temporal direction, we generate the gauge configurations extracted using the ratio of the ηc two-point function and
on anisotropic lattices with the temporal lattice much the related three-point function evaluated at Q2 = 0,
finer than the spatial lattice, say, ξ = as /at 
(2)
1, where as and at are the spatial and temporal (µ) p; tf = n2t )
pµ 12 Γηc (~
lattice spacings, respectively. In practice, we choose ZV (t) = , (8)
p, p~, n2t , t)
p) Γ(3),µ (~
E(~
ξ = 5. The gauge action we use is the tadpole
improved gauge action [18] whose discretization error where the factor 1/2 accounts for the effect of the tempo-
is expected to be O(a4s , αs a2s ). Three gauge ensembles ral periodic boundary condition, and the superscript µ of
with large statistics are generated at different lattice ZV (as ) is used to differentiate the temporal component
spacings for the continuum extrapolation, and the from the spatial ones, since they are not necessarily
relevant ensemble parameters are listed in TableI, where the same on the anisotropic lattices. Figure 1 plots
the lattice spacings as are determined from r0−1 = (µ)
ZV (t) with respect to t for the three lattices. ZV ’s
410(20) MeV by calculating the static potential. For
are extracted from the plateaus and the values are listed
fermions, we use the tadpole improved clover action for
in Table II. Obviously, the renormalization constant
anisotropic lattices [19]. The parameters in the action are (s)
tuned carefully by requiring that the physical dispersion ZV (a) of the spatial components of the vector current
(t)
relations of vector and pseudoscalar mesons are correctly deviate from that of the temporal component, ZV (a)
reproduced at each bare quark mass [20]. The bare charm by a few percents. This deviation is understandable
quark masses for the two lattices are set by the physical since the space-time interchange symmetry is broken on
mass of J/ψ, mJ/ψ = 3.097 GeV. anisotropic lattices, apart from the imperfect tuning of
In the quenched approximation, since there are no the bare speed of light in the fermion action. In practice,
sea quarks, the electromagnetic current contributing to we perform the calculation in the rest frame of the final
the radiative transitions of charmonia involves only the pseudoscalar state. In this frame, according to Eq. (3),
charm quark, say, Jem µ
(x) = ec J µ (x) with J µ (x) = the matrix element of the temporal component of the
c̄γµ c(x), which is the one we adopt in this study. It is a vector current is zero due to the appearance of the
conserved vector current and need not be renormalized totally antisymmetric tensor µνρσ , and we extract the
in the continuum. However, on a finite lattice, it is form factor M (Q2 ) only from the corresponding matrix
not conserved anymore due to the lattice artifacts and elements of the spatial components of the vector current.
(s)
receives a multiplicative renormalization factor ZV (as ). So we only need the renormalization factor ZV .
 4

2.0 Energy of J/ on = 3.0 lattice

0.9
ZVt, fit(0, 0, 0) = 1.2877(72) n2=0
ZVt, fit(1, 0, 0) = 1.2935(85) n2=1
1.8 ZVs, fit(1, 0, 0) = 1.3434(80) 0.8 n2=2
n2=3
n2=4
1.6 0.7 n2=5
n2=6
Z (t)

n2=8
0.6

E(p)
1.4 n2=9
n2=12
1.2 0.5

0.4
1.0
0 10 20 30 40
t
0.3
0 10 20 30 40 50 60 70 80
t
1.6
ZVt, fit(0, 0, 0) = 1.15577(94)
ZVt, fit(1, 0, 0) = 1.1527(77)
1.4 ZVs, fit(1, 0, 0) = 1.1019(66) FIG. 2. Effective energy plateaus E(~ p) of J/ψ for β = 3.0.
E(~p)’s are averaged over the momentum modes with the same
1.2 p|2 . It is seen that signal-to-noise ratios are good for the
|~
momentum modes we are using.
Z (t)

1.0

0.8 through
(2)
X
0.6 p, t) = −
Γij (~ e−i~p·~x TrhS † (~x, t; ~0, 0)Γi S(~x, t; ~0, 0)Γj i
0 10 20 30 40 50 60 70
t ~
x
(9)
1.6 where S(~x, t, ~0, 0) is the point-source propagator of the
ZVt, fit(0, 0, 0) = 1.1084(38)
ZVt, fit(1, 0, 0) = 1.1142(55) charm quark. The effective energy plateaus E(~ p) of J/ψ
1.4 ZVs, fit(1, 0, 0) = 1.0623(64) are illustrated in Fig. 2 for β = 3.0, where E(~ p)’s are
averaged over the momenta with the same |~ p|2 . We check
1.2 the dispersion relation of J/ψ by calculating the squared
Z (t)

speed of light ĉ2

1.0
E 2 (~
p) − m2J/ψ
0.8 ĉ2 = . (10)
p|2
|~
0.6
0 10 20 30 40 50 60 70 80 It is found that the largest deviation of ĉ2 from one is
t
less than 4% on all the three lattices. For illustration,
(s) (t) we plot ĉ2 with respect to different momentum modes on
FIG. 1. The renormalization factor ZV and ZV with
the β = 3.0 lattice in Fig. 3, where the data points are
respect to t for the three lattices. The panels from top to
bottom correspond to β = 2.4, 2.8, 3.0, respectively the averaged values over the momentum modes n̂ with
the same |n̂|2 by assuming the approximate rotational
symmetry.
(3),µi
A. J/ψ → γηc transition The three-point functions ΓJ/ψ→γηc contributed by
the connected diagrams (disconnected diagrams are
neglected) are calculated through the expression
There have been quite a few lattice studies on the
decay process J/ψ → γηc . We would like to carry the (3),µi
ΓJ/ψ→γηc (~
pi , ti = 0; p~f , tf ; ~q, t)
similar calculation and make a comparison with previous
e−i~pf ·~x+i~q·~y hOΓf (~x, tf )(c̄γµ c)(~y , t)Oi† (~0, 0)i
X
studies as well as the experimental value, which serve as =
a calibration to some discretization uncertainties of our ~
x,~
y
lattice setup. We work in rest frame of the pseudoscalar
ei~q·~y trhHΓ† f (~y , t; ~0, 0; p~f , tf )γ5 γµ S(~y , t; ~0, 0)γi i
X
=−
(such as ηc ) with J/ψ moving with a definite momentum
~
y
p̂i = 2πn̂/Las , where n̂ ranges from (0, 0, 0) to (2, 2, 2).
As mentioned above, both the three-point functions and (11)
two-point functions are required in order to extract the where
µ
desired hadronic matrix elements of the current Jem . We X
choose the quark bilinear operators Oi (x) = (c̄Γi c)(x) for HΓf (y, 0; tf ; p~f ) ≡ ei~pf ·~x S(~y , t; ~x, tf )γ5 S(~x, tf ; ~0, 0)γ5
ηc (Γi = γ5 ) and J/ψ (Γi = γi , i = 1, 2, 3), such that the ~
x
two-point function with momentum p~ can be calculated (12)
 5

TABLE III. The mass of J/ψ, ηc and the form factor V̂ (0) on
c2 = E (p)|p|2mJ/
2 2

1.04 three lattices.

β mJ/ψ (GeV) mηc (GeV) V̂ (0)
1.02 2.4 3.097(1) 2.995(1) 2.152(34)(107)
2.8 3.102(1) 3.007(2) 1.962(14)
1.00 3.0 3.105(1) 2.995(1) 1.971(18)
c2

∞ 1.933(41)
0.98

0.96 three lattices we are using are plotted in terms of Q2 in

2 4 6 8 10 12 Fig. 4. The β = 2.4 lattice is coarse (L3 × T = 83 × 96),
n2 such that results has larger systematic errors. We try
to fit with the data at the smallest three Q2 and get
FIG. 3. The squared speed of light ĉ2 versus the momentum V̂ (Q2 ) = 2.045(36). We regard the difference of the fitted
modes n̂ ranging from (0, 0, 0) to (2, 2, 2) on the β = 3.0 values using different Q2 range as systematic error. The
lattice. Each data point shows the averaged value of ĉ2 over errors on other two lattices are statistical errors.
the momentum modes n̂ with the same |n̂|2 by assuming the
In order to obtain the on-shell form factor V̂ (Q2 =
approximate rotational symmetry. The largest deviation of
ĉ2 from one is less than 4%.
0), we adopt the following function form to do the
extrapolation,
Q 2
can be obtained by the sequential source technique [21]. − 16β
V̂ (Q2 ) = V̂ (0)e 2
(15)
In order to increase the statistics, we repeat the same
calculations T times (where T is the temporal lattice size) which is inspired by the simple quark model with the
by setting a point source on a different time slice each harmonic oscillator wave functions of ηc and J/ψ, as
time. With the related two-point functions calculated addressed in Ref. [21]. The extrapolations are also
accordingly, a straightforward way to extract the inter- illustrated in Fig. 4 by curves with error bands. It is
µ
esting matrix elements hηc (~ pf , λf )|Jem (0)|J/ψ(~ pi , λi )i is interesting to see that this kind of function form describes
to fit the three-point function and two-point function the data very well (On our coarsest β = 2.4 lattice, there
simultaneously according to Eq. (6) and Eq. (7) through is a clear deviation from the curve, which is tentatively
the jackknife analysis. To suppress the contribution of attributed to the relatively large discretization error of
excited states, we use this formula Q2 calculated on this lattice). The results after the
v extrapolation are listed in Table III and shown in Fig. 5.
u 2Ei Γ(2) (~
u
µi (3),µi pi , tf − t) Since we have three lattices with different lattice spacings
R (~q, tf , t) = Γ (~q, tf , t)t (2) i (2) as , we also perform a linear extrapolation
Γi (~ pi , t)Γi (~ pi , tf )
V̂ (0, as ) = V̂ (0)cont. + Aa2s
v
u (2) (16)
u 2Ef Γf (~ pf , t)
× t (2) (2)
(13)
Γf (~ pi , tf − t)Γi (~ pi , tf ) to get the final result of V̂ (0)cont. ,

which gives flatter plateaus. In practice, the energies Ei,f V̂ (0)cont. = 1.933(41), (17)
are derived from two-point functions in the joint fit of the
two-point and three-point functions. We can get the form from which we give the prediction of the partial width of
factors by solving the Eq. (3). Based on the OZI rule, the process J/ψ → γηc ,
we neglect the contribution from the quark annihilation
diagrams and only consider the contribution of connected α|~k|3 64
diagrams. As such we compute the form factor V̂ (Q2 ) ΓJ/ψ→γηc = 2
|V̂ (0)|2
which is related to V (Q2 ) by (mJ/ψ + mηc ) 27
= 2.47(11)keV, (18)
2
V (Q ) = 2 × e × V̂ (Q2 ),
2
(14)
3 where the fine coupling constant takes the value at the
1
charm quark mass scale, α = 134 , and mJ/ψ and mηc
where the factor 2 comes from the insertion of the assume the experimental values.
electromagnetic current to both the quark and antiquark We compare our result with those from previous lattice
lines, 2/3e is the electric charge of charm quark. In QCD studies in Table IV where one can see that all the
the expression above, the renormalization constant of the results reach a consensus within errors. This assures us
µ (s)
spatial component of Jem , say, ZV , has been implicitly that, in our study, the systematic uncertainties are not
incorporated into V̂ (Q2 ). The extracted V̂ (Q2 ) on the important in the charmonium sector.
 6

2.2 V(Q 2 = 0) = 2.152(34) 2.0 V(Q 2 = 0) = 1.962(14) 2.0 V(Q 2 = 0) = 1.971(18)

2.0 V(Q 2) 1.8 V(Q 2) V(Q 2)
1.8
1.8 1.6
1.6
1.6 1.4
1.4
V(Q 2)

V(Q 2)

V(Q 2)
1.4 1.2
1.2 1.2
1.0
1.0 1.0
0.8
0.8 0.8
0.6
0.6
0.6
0 1 2 3 4 0 1 2 3 4 5 0 1 2 3 4 5
Q 2(GeV 2) Q 2(GeV 2) Q 2(GeV 2)

FIG. 4. The form factor V̂ (Q2 ) in Eq. (3) of J/ψ to γηc where Q2 = −(PJ/ψ − Pηc )2 . Exponential functions in Eq. (15) are
adopted for extrapolation. The panels from left to right corresponds to β = 2.4, 2.8, 3.0 lattices, respectively.

glueball operators, which is outlined as follows. First,

2.6 the last four of the ten Wilson loops in Fig. 3 of Ref. [2]
VJ/ c (0) = 1.933(41)
are used as prototypes, and then six smearing schemes
2.4 (different combinations of the singlelike smearing and
the double-link smearing , as addressed in Ref. [1])
2.2 are applies to each of these prototype loops. Thus we
V(Q 2 = 0)

obtain 24 different Wilson loops as the basis operators.

2.0 Since the lattice counterpart of the quantum number
0−+ in the continuum is A−+ 1 , where A1 is one of
1.8 the five irreducible representations A1 , A2 , E, T1 , T2 of
the spatial symmetry group O of the cubic lattice, we
1.6
apply the 24 operations of O group to each of the basis
0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 operators and obtain 24 copies of it, whose proper linear
as2(fm2)
combination gives the representation of A−+ 1 . Thereby
we get 24 different operators with the quantum number
FIG. 5. The continuum extrapolation of form factor V̂ (Q2 )
of J/ψ → γηc on three lattices. The lattice spacing as is
A−+
1 , which compose a operator set {φα , α = 1, 2, ..., 24}.
in physical units, f m. The blue band shows the continuum Based on this operator set, we calculate the matrix of the
extrapolation linear in a2s . correlation functions C̃(t) = {C̃αβ } with

1X
C̃αβ (t) = hφα (t + τ )φβ (τ )i (19)
T τ
TABLE IV. Comparison of V̂ (0) with previous lattice results.

Lattice setup V̂ (0) where we sum over τ to increase the statistics. Finally,
Nf = 2 + 1 [22] 1.90(7)(1) by solving the generalized eigenvalue problem
Nf = 2 [23] 1.92(3)(2)
Nf = 2 [24] 2.01(2) C̃(tD )V = λ(tD )C̃(0)V, (20)
Quenched [21] 1.85(4)
Quenched (this work) 1.933(41) we can get the eigenvector V = {vα , α = 1, 2, . . . , 24}
corresponding to the maximal eigenvalue λmax (tD ) ≡
e−mmin (tD )tD (mmin (tD ) is close to the mass of the ground
state), from which we can obtain the optimal operator
B. The partial decay width of J/ψ radiatively Φ(t) for the ground state P pseudoscalar glueball |Gi by
decaying into the pseudoscalar glueball the combination Φ(t) = α vα φa . In this work, we
set tD = 1. The correlation function of Φ(t) can be
We extend the similar study to the process of J/ψ parametrized as
radiatively decaying into the pseudoscalar glueball. It
1X
is known that the signals of glueballs are always noisy, C(t) = hΦ(t + τ )Φ† (τ )i
such that a large statistics is required. On the T τ
other hand, an optimal interpolation operator, which 2
ZG
couples predominantly to the ground state |Gi of the ' e−mG t = W e−mG t , (21)
2mG V3
pseudoscalar glueball, is mandatory for us to extract the
µ
desired matrix element hG|Jem (0)|J/ψi reliably from the where mG is the mass of the ground state, V3 = L3 a3s is
related three-point functions in Eq. (6). In doing so, we the spatial volume of the lattice, and ZG = |h0|Φ(0)|Gi|.
adopt the strategy used in [1, 2] to construct the optimal The effective mass of C(t) is plotted in Fig. 6, where one
 7

suppressed to some extent. Practically, we fix the time

1.2 interval of glueball operator and vector current operator
= 2.4, Mfit=0.6130(40)
= 2.8, Mfit=0.3633(12) as one time slice which mean tf − t = 1, since that the
1.0
= 3.0, Mfit=0.2735(12) optimal glueball operator projects almost totally on the
= 2.4, data
0.8 = 2.8, data ground state pseudoscalar glueball. After the matrix
= 3.0, data elements hG|J µ |J/ψi are derived, in analogy with the
0.6 J/ψ → γηc case, we can obtain the form factors V̂ (Q2 )
M

at different Q2 , where Q2 = −(pJ/ψ − pG )2 (note that

0.4
pG = (mG , ~0) since the pseudoscalar is at rest). The form
0.2 factor V̂ (Q2 ) derived on the three lattices are plotted in
Fig. 7 with respect to various Q2 .
0.0
1 2 3 4 5 6 7 8 9
Because only the on-shell form factor, say, V̂ (Q2 = 0),
t enters into the formula of the transition width of J/ψ →
γG0−+ ,
FIG. 6. The effect mass of pseudoscalar glueball extracted
16 |~k|3
from the two-point correlation function which constructed by Γ(J/ψ → γG0−+ ) = α |V̂ (0)|2 , (24)
optimal glue operators on three lattices. 27 (mG + mJ/ψ )2
we should perform an extrapolation of V̂ (Q2 ) from Q2 6=
can see that the plateau almost starts from the beginning 0 to Q2 = 0. However, in contrast to the case of
of time t. Obviously, there is still some contribution from J/ψ → γηc where the extrapolation function form of
higher states which manifests by the slight increment of V̂ (Q2 ) (Eq. (15)) can be inspired by the wave functions of
the effective mass toward t = 0. To check the extent J/ψ and ηc in the nonrelativistic quark model, we have no
of the higher state contamination, we fit C(t) through theoretical information for the Q2 dependence of V̂ (Q2 )
a single-exponential function and find that the deviation in J/ψ radiatively decaying into glueballs through the c̄c
of W from one is at a level of few percents (note that annihilation. Anyway, it is seen in Fig. 7 that V̂ (Q2 )
C(t) is normalized as C(0) = 1). This means that C(t) depends mildly on Q2 , therefore a polynomial fit in Q2
is almost totally dominated by the contribution from can be safe here. The extrapolation formula is taken as
the ground
√ state and therefore W ' 1, or equivalently,
V̂ (Q2 ) = V̂ (0) + aQ2 + bQ4 . (25)
ZG ' 2mG V3 can be a good approximation. The mass
of the ground state pseudoscalar glueball on the three which can describe the data satisfactorily for all the
lattices are listed in Table V . We obtain the mass of three lattices, as shown in Fig. 7 with blue bands. The
pseudoscalar glueball as 2.395(14) GeV after continuum extrapolated V̂ (0) on the three lattices are listed in
extrapolation. This value is lower than that in Ref. [2], Table V. In order to get the V̂ (0) in the continuum limit,
but it is consistent within errors. we also carry out a linear extrapolation in a2s
With the optimal operator Φ(t), the relevant three-
point function can be calculated and parametrized as V̂ (0) = V̂ (0)cont. + ca2s .

(3),µi Figure 8 shows V̂ (0) (red data points) at different lattice

ΓJ/ψ→γG −+ (~
pi , ti = 0; p~f , tf ; ~q, t) spacings and its continuum limit (blue point), where
0

1 X i~q·~y the blue band illustrates the linear extrapolation in a2s .

= e hΦ(tf + τ ) × J µ (~y , t + τ )Oi† (~0, τ )i
T Finally, the form factor V̂ (0) in the continuum limit is
τ,~
y
determined to be
X e−mG (tf −t) e−EV t
' hΩ|Φ(0)|Gi V̂ (0)cont. = 0.0246(43), (26)
2mG V3 2EV
V,r
which gives the decay width of J/ψ → γG0−+
×hG|J µ (0)|V (~ pi , r)|Oi† (0)|Ωi.
pi , r)ihV (~ (22) (
0.0215(74) keV mG = 2.395 GeV
Similar to the case of J/ψ to ηc γ, we use this formula Γ(J/ψ → γG0−+ ) =
0.0099(34) keV mG = 2.56 GeV
(27)
p
4V3 mG EJ/ψ (~q) at different pseudoscalar masses according to Eq. (24).
µi (3),µi
R (~q, tf , t) = Γ (~q, tf , t) Consequently, the production fraction of the pseu-
C(tf − t)
v doscalar glueball in the J/ψ radiative decay is estimated
(2) to be
ΓJ/ψ (~q, tf − t)
u
u
×t (23)
(
(2) (2)
ΓJ/ψ (~q, t)ΓJ/ψ (~q, tf ) 2.31(80) × 10−4 mG = 2.395 GeV
Br(J/ψ → γG0−+ ) = .
1.07(37) × 10−4 mG = 2.56 GeV
to extract the matrix elements hG|J µ |J/ψi, through (28)
which the contribution from excited states can be
 8

0.06
TABLE V. The mass of pseudoscalar glueball, mG , and the
V(Q 2 = 0) = 0.0307(59) form-factor V̂ (0) of J/ψ → γGps . The continuum limits of
V(Q 2) mG and V̂ (0) are achieved by linear extrapolations in a2s .
0.05
β mG (GeV) V̂ (0)
0.04 2.4 2.724(18) 0.0307(59)
V(Q 2)

2.8 2.550(13) 0.0294(32)

0.03 3.0 2.464(11) 0.0247(33)
Continuum limit 2.395(14) 0.0246(43)
0.02
2.560(35)(120)[2]

0.01
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
Q 2(GeV 2) 0.08
V c(0) = 0.0246(43)
0.07
0.06
V(Q 2 = 0) = 0.0294(32) 0.06
V(Q 2)
0.05 0.05

V(Q 2 = 0)
0.04
0.04
0.03
V(Q 2)

0.03 0.02

0.01
0.02 0.00
0.00 0.02 0.04 0.06 0.08 0.10
as2(fm2)
0.01
0.0 0.5 1.0 1.5 2.0
Q 2(GeV 2)
FIG. 8. The continuum extrapolation of the form-factor
V (Q2 = 0) of J/ψ → γGps . The result show that the form-
0.06 factor V (0) is very weak dependence to lattice spacing.
V(Q 2 = 0) = 0.0247(33)
V(Q 2)
0.05
decaying into a pseudoscalar X is through the M 1 decay,
0.04 such that the partial width is proportional to |~k|3 with
|~k| = (m2J/ψ − m2X )/(2mJ/ψ ) being the magnitude of the
V(Q 2)

0.03 decaying momentum of the final state photon. Because

the mass 2.4 GeV of the pseudoscalar glueball is close
0.02 to the mass of J/ψ, the partial width is suppressed by
the kinematics. In order to obtained a fair comparison,
0.01 we would like to subtract the phase space factor and
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75
Q 2(GeV 2) introduce an effective coupling gX through the definition

FIG. 7. The form factor V̂ (Q2 ) and the extrapolated value 1 2 |~k|3
Γ(J/ψ → γX) = αg , (29)
V̂ (Q2 = 0) of J/ψ → γG0−+ where Q2 = −(pJ/ψ − pG )2 . We 3 X m2J/ψ
use the function V̂ (Q2 ) = V̂ (0) + aQ2 + bQ4 to perform the
extrapolation. where 1/3 accounts for the spin average of the J/ψ, α is
the fine coupling constant. Obviously, gX describes the
coupling of the gluons generated through cc̄ annihilation
IV. DISCUSSION to the pseudoscalar meson X. Since in experiments
only the branching fractions Br(J/ψ → γX) can be
Obviously, production fraction of the pure gauge measured directly, we express gX explicitly in terms of
pseudoscalar glueball is quite small in the J/ψ radiative these branching fractions along with the total width ΓJ/ψ
decays, especially when comparing with 3.8(9) × 10−3 as
of the pure gauge scalar glueball and roughly 1% of " 5
# 12
the tensor glueball. Furthermore, this value is also 24ΓJ/ψ Br(J/ψ → γX)mJ/ψ
gX = . (30)
much smaller than those of most known pseudoscalars. α (m2J/ψ − m2X )3
For example, the branching fraction of J/ψ → γη 0
is 5.13(17) × 10−3 , which is one orders of magnitude In practice, for η, η 0 , and η(2225), we take the branch
larger. One of the reasons is that J/ψ radiatively fractions directly from the PDG data. For η(1405) and
 9

η(1475), their production rates in the J/ψ decays are not at mπ = 650 MeV. In Ref. [29], the authors use both the
differentiated from each other in PDG, so we also take fermionic operator ψ̄γ5 ψ and q(x), and get compatible
them as a whole and sum over the branching fractions results for mη2 on several gauge ensembles, whose chiral
of final states γK K̄π, γγρ, γηπ + π − , and γρ0 ρ0 . For extrapolation gives mη2 = 772(18) MeV.
η(1760) we add up the branching fractions of γρ0 ρ0 and The discussion above helps to understand why the
γωω, for X(1835) we add use the sum of the branching couplings of gluons to η states are not suppressed in
ratios of the final states γ + (π + π − η 0 , pp̄, ηKS KS ). The comparison with the coupling of the pseudoscalar glueball
branch fractions of J/ψ radiative decay are listed in in the J/ψ radiative decays. This also implies the
Table VI and the derived gX ’s are listed in Table VII. coupling gX cannot be used as a characteristics for
The striking observation is that the coupling gX for the identifying the pseudoscalar glueball. Anyway, since the
pseudoscalar glueball is comparable with or smaller than lattice QCD studies predicts the pseudoscalar glueball
those of the known nonflavored pseudoscalars (note that mass is around 2.4-2.6 GeV, we may wish to check if there
the smallness of gη is due to the dominance of the flavor are any candidates in this mass region. With the world
octet component of η(547)). This is in sharp contrast largest J/ψ event ensemble, the BESIII Collaboration
to the usual expectation based on the naive αs -power is performing a scrutinized partial wave analysis on
counting that the gluons in the J/ψ radiative decay the radiative J/ψ decay processes. In the process
couple more strongly to glueballs than q q̄ mesons. Even Jψ → γη 0 π + π − , BESIII observes resonancelike structure
though this naive αs -power counting is not justified in the X(2120) and X(2370) in the invariant mass spectrum
low energy QCD regime, empirically there is an OZI rule of η 0 π + π − [17, 30] , however, their spin and parity have
observed in many hadronic processes that the processes not been determined yet. X(2370) is confirmed in
mediated through gluons, say, involving the Feynman the invariant mass spectra of η 0 K K̄ by BESIII [31],
diagrams without continuous quark lines connecting the and its production fractions in the processes J/ψ →
initial and final states, tend to be strongly suppressed. γX(2370) → γη 0 K + K − and J/ψ → γX(2370) →
In this sense, the production of η states in the J/ψ γη 0 KS KS have been estimated to be [1.86 ± 0.39(sta.) ±
radiative decays seems to be OZI-violated. The QCD 0.29(syss.)] × 10−5 and [1.19 ± 0.37(sta.) ± 0.18(syss.)] ×
chiral anomaly may play an important role in these 10−5 , respectively. In the partial wave analysis of the
processes. In QCD, the flavor singlet axial vector current process J/ψ → γφφ, BESIII also observes a component
N
Pf X(2500), whose preferred spin-parity assignment is
j5µ (x) = q̄k (x)γ5 γ µ qk (x) is not conserved, pseudoscalar [16]. Its production fraction in the process
k=1
J/ψ → γX(2500) → γφφ is determined to be [1.7 ±
Nf
X 0.2+0.2
−0.8 ]×10
−5
. Whether they are the same object or not,
∂µ j5µ (x) = 2Nf q(x) + 2imk q̄k (x)γ5 qk (x), (31) X(2370) and X(2500) reside in the pseudoscalar glueball
k=1 mass region predicted by lattice QCD, whose production
even in the chiral limit mk → 0, where Nf is the rates are compatible with that of the pseudoscalar
number of the quark flavor, and the anomalous gluonic glueball in this work. Recently, BESIII has finished
g2 µν the data collection of 10 billion J/ψ events, hopefully
term q(x) = 16π 2 trGµν (x)G̃ (x) comes either from the
the properties of these states can be determined more
regularization of the linearly divergent one-loop diagrams
precisely in the near future.
of the vector-vector-axial vector current vertex (the
triangle diagram) in the perturbation theory or the chiral
transformation noninvariance of the fermion measure in
V. SUMMARY
the path integral formalism. As shown in Ref. [26, 27],
the matrix element related to the chiral anomaly can
be sizeable, and obviously violated the OZI rule. In We carry out the first lattice calculation of the produc-
other words, the QCD UA (1) anomaly can enhance the tion rate of the pseudoscalar glueball in J/ψ radiative
coupling of η states to gluons, and this nonperturbative decays in the quenched approximation. We generate
effect results in the violation of the OZI rule when flavor gauge configurations on three anisotropic lattices with
singlet pseudoscalar mesons are involved. different lattice spacings, which facilitate us to perform
the continuum limit extrapolation. In order to calibrate
There is also evidence from lattice QCD study that the some of the systematic uncertainties, we first calculate
topological charge density q(x) couples to η 0 strongly. partial decay width of J/ψ → γηc . The related M1
A Nf = 2 + 1 lattice simulation uses q(x) to study η 0 transition form factor is determined to be V̂ (0) =
and observes a clear contribution of η 0 to the correlation 1.933(41), which gives the partial width Γ(J/ψ → γηc ) =
function hq(x)q(0)i [28]. The authors get the result 2.47(11) keV. These results are consistent with the results
mη0 = 1019(119) MeV at the physical pion mass of previous lattice calculations.
consistent with the physical η 0 mass. Similar lattice By applying the variational method to a large operator
studies have also been carried out in the Nf = 2 set, we obtain an optimal operator which couples
case [9, 29]. In Ref. [9], the authors get the mass of predominantly to the ground state pseudoscalar glueball
the isoscalar pseudoscalar η2 to be mη2 = 890(38) MeV G. In this work, mG is determined to be 2.395(14) GeV,
 10

TABLE VI. The branch fractions of J/ψ radiative decay to pseudoscalar light mesons from the PDG data[25].
pseudoscalar mesons final states branching ratios
η γη (1.104 ± 0.034) × 10−3
η0 γη 0 (5.13 ± 0.17) × 10−3
γη(1405/1475) → γK K̄π (2.8 ± 0.6) × 10−3
γη(1405/1475) → γγρ0 (7.8 ± 2.0) × 10−5
η(1405/1475)
γη(1405/1475) → γηπ + π − (3.0 ± 0.5) × 10−4
γη(1405/1475) → γρ0 ρ0 (1.7 ± 0.4) × 10−3
γη(1760) → γρ0 ρ0 (1.3 ± 0.9) × 10−4
η(1760) −3
γη(1760) → γωω (1.98 ± 0.33)

× 10−4
γX(1835) → γπ + π − η 0 +0.34
2.77−0.40
× 10
−5
X(1835) γX(1835) → γpp 7.7+1.5
−0.9
× 10
−5
γX(1835) → γKS0 KS0 η 3.3+2.0
−1.3
× 10
η(2225) γη(2225) 3.14−0.19 × 10−4
+0.50

QCD UA (1) anomaly which is totally a nonperturbative

TABLE VII. The gX of flavor-singlet pseudoscalar mesons.
effect.
Pseudoscalar (X) gX Even though this study is performed in the quenched
η 0.0108(2) approximation and the uncertainty in the presence
η0 0.0259(8) of dynamical quarks is not controlled, we hope our
η(1405/1475) 0.0313(41) result can provide useful theoretical information for
η(1760) 0.0255(25) experiments to unravel the properties of the possible
X(1835) 0.0123(12)
pseudoscalar glueball.
η(2225) 0.0167(17)

G0−+ 0.0126(22)
ACKNOWLEDGEMENTS

L.C. thanks J. Liang for useful discussions. This

and the on-shell form factor of J/ψ → γG is derived work is supported by the National Key Research and
as V̂ (0) = 0.0246(43), in the continuum limit, from Development Program of China (No.2017YFB0203202).
which we obtain the following partial decay width and The numerical calculations are carried out on Tianhe-
the production rate 1A at the National Supercomputer Center (NSCC)
in Tianjin and the GPU cluster at Hunan Normal
Γ(J/ψ → γG0−+ ) = 0.0215(74) keV University. We acknowledge the support of the National
Science Foundation of China (NSFC) under Grants
Br(J/ψ → γG0−+ ) = 2.31(80) × 10−4 (32) No. 11575196, No. 11405053, No. 11335001, No.
11621131001 (CRC 110 by DFG and NSFC), No.
We introduce an effective coupling gX to describe the U1832173, No. 11705055. Y.C. is also supported by the
interaction between the pseudoscalar X and the gluonic CAS Center for Excellence in Particle Physics (CCEPP)
intermediate states in the processes J/ψ → γX, as and National Basic Research Program of China (973
defined in Eq. (30). It is interesting to see that all the Program) under code number 2015CB856700. Y.Y.
gX ’s are comparable for the pseudoscalar glueball and the is supported by the CAS Pioneer Hundred Talents
nonflavored q q̄ pseudoscalars (η states). We tentatively Program. Our matrix inversion code is based on QUDA
attribute the large production rates of the η states to the libraries[32]and the fitting code is based on lsqfit[33] .

[1] C. J. Morningstar and M. Peardon, Phys. Rev. D 60, [4] S. Nussinov and R. Shrock, Phys. Rev. D80, 054003
034509 (1999), arXiv:9901004 [hep-lat]. (2009), arXiv:0907.1577 [hep-ph].
[2] Y. Chen, A. Alexandru, S. J. Dong, T. Draper, [5] C. M. Richards, A. C. Irving, E. B. Gregory, and C. Mc-
I. Horváth, F. X. Lee, K. F. Liu, N. Mathur, Neile, Phys. Rev. D 82, 34501 (2010), arXiv:1005.2473v2.
C. Morningstar, M. Peardon, S. Tamhankar, B. L. Young, [6] E. Gregory, A. Irving, B. Lucini, C. McNeile, A. Rago,
and J. B. Zhang, Phys. Rev. D 73, 014516 (2006), C. Richards, and E. Rinaldi, J. High Energy Phys. 10,
arXiv:0510074 [hep-lat]. 170 (2012), arXiv:1208.1858.
[3] G. Hao, C.-F. Qiao, and A.-L. Zhang, Phys. Lett. B642, [7] W. I. Eshraim, S. Janowski, F. Giacosa, and D. H.
53 (2006), arXiv:hep-ph/0512214 [hep-ph]. Rischke, Phys. Rev. D 87, 054036 (2013).
 11

[8] C.-F. Qiao and L. Tang, Phys. Rev. Lett. 113, 221601 [22] G. C. Donald, C. T. H. Davies, R. J. Dowdall,
(2014), arXiv:1408.3995 [hep-ph]. E. Follana, K. Hornbostel, J. Koponen, G. P. Lepage,
[9] W. Sun, L.-C. Gui, Y. Chen, M. Gong, C. Liu, Y.-B. Liu, and C. McNeile, Phys. Rev. D 86, 094501 (2012),
Z. Liu, J.-P. Ma, and J.-B. Zhang, Chin. Phys. C 42, arXiv:1208.2855.
093103 (2017), arXiv:1702.08174. [23] D. Becirevic and F. Sanfilippo, J. High Energy Phys. 01,
[10] L. Olbrich, M. Zétényi, F. Giacosa, and D. H. Rischke, 28 (2012), arXiv:1206.1445.
Phys. Rev. D 97, 014007 (2018). [24] Y. Chen, D.-C. Du, B.-Z. Guo, N. Li, C. Liu, H. Liu,
[11] E. Klempt and A. Zaitsev, Phys. Rep. 454, 1 (2007), Y.-B. Liu, J.-P. Ma, X.-F. Meng, Z.-Y. Niu, and J.-B.
arXiv:0708.4016 [hep-ph]. Zhang, Phys. Rev. D 84, 034503 (2011), arXiv:1104.2655.
[12] W. Ochs, J. Phys. G40, 043001 (2013), arXiv:1301.5183 [25] M. Tanabashi et al. (Particle Data Group), Phys. Rev.
[hep-ph]. D 98, 030001 (2018).
[13] W. Qin, Q. Zhao, and X.-H. Zhong, Phys. Rev. D 97, [26] M. Gong, Y.-B. Yang, J. Liang, A. Alexandru, T. Draper,
096002 (2018), arXiv:1712.02550. and K.-F. Liu (χQCD Collaboration), Phys. Rev. D 95,
[14] L. C. Gui, Y. Chen, G. Li, C. Liu, Y. B. Liu, J. P. Ma, 114509 (2017).
Y. B. Yang, and J. B. Zhang, Phys. Rev. Lett. 110, [27] Y.-B. Yang (χQCD), Pos LATTICE2018, 017 (2019)
021601 (2013), arXiv:1206.0125. , arXiv:1904.04138 [hep-lat].
[15] Y. B. Yang, L. C. Gui, Y. Chen, C. Liu, Y. B. Liu, J. P. [28] H. Fukaya, S. Aoki, G. Cossu, S. Hashimoto, T. Kaneko,
Ma, and J. B. Zhang, Phys. Rev. Lett. 111, 091601 and J. Noaki (JLQCD Collaboration), Phys. Rev. D 92,
(2013), arXiv:1304.3807. 111501 (2015).
[16] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. D [29] P. Dimopoulos, C. Helmes, C. Jost, B. Knippschild,
93, 112011 (2016), arXiv:1602.01523 [hep-ex]. B. Kostrzewa, L. Liu, K. Ottnad, M. Petschlies,
[17] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. C. Urbach, U. Wenger, and M. Werner, Phys. Rev. D
Lett. 106, 072002 (2011), arXiv:1012.3510. 99, 34511 (2019).
[18] C. J. Morningstar and M. Peardon, Phys. Rev. D 56, [30] M. Ablikim et al. (BESIII Collaboration), Phys. Rev.
4043 (1997), arXiv:9704011 [hep-lat]. Lett. 117, 042002 (2016).
[19] C. Liu, J. Zhang, Y. Chen, and J. Ma, Nucl. Phys. B [31] V. Prasad, in 16th Conf. Flavor Phys. CP Viol. (2018)
624, 360 (2002). arXiv:1810.00180 [hep-ex].
[20] J. Zhang and C. Liu, Mod. Phys. Lett. A 16, 1841 (2001), [32] M. Clark, R. Babich, K. Barros, R. Brower, and
arXiv:0107005 [hep-lat]. C. Rebbi, Comput. Phys. Commun. 181, 1517 (2010),
[21] J. J. Dudek, R. G. Edwards, and D. G. Richards, Phys. arXiv:0911.3191.
Rev. D 73, 74507 (2006), arXiv:0601137 [hep-ph]. [33] G. P. Lepage, B. Clark, C. T. H. Davies, K. Hornbostel,
P. B. Mackenzie, C. Morningstar, and H. Trottier,
Nucl. Phys. Proc. Suppl. 106, 12 (2002), arXiv:hep-
lat/0110175 [hep-lat].