Table Of ContentTwo Photon Decays of η from Lattice QCD
c
Ting Chen,1 Ying Chen,2 Ming Gong,2 Yu-Hong Lei,1 Ning Li,3 Chuan Liu,4,5,∗ Yu-Bin
Liu,6 Zhaofeng Liu,2 Jian-Ping Ma,7 Wei-Feng Qiu,2 Zhan-Lin Wang,1 and Jian-Bo Zhang8
(CLQCD Collaboration)
1School of Physics, Peking University, Beijing 100871, China
2Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
3School of Science, Xi’an Technological University, Xi’an 710032, China
4School of Physics and Center for High Energy Physics, Peking University, Beijing 100871, China
5Collaborative Innovation Center of Quantum Matter, Beijing 100871, China
6School of Physics, Nankai University, Tianjin 300071, China
7Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
8Department of Physics, Zhejiang University, Hangzhou 311027, China
We present an exploratory lattice study for the two-photon decay of η using N = 2 twisted
c f
6 mass lattice QCD gauge configurations generated by the European Twisted Mass Collaboration.
1 Two different lattice spacings of a = 0.067fm and a = 0.085fm are used in the study, both of
0
which are of physical size of 2fm. The decay widths are found to be 1.025(5)KeV for the coarser
2
lattice and 1.062(5)KeV for the finer lattice respectively where the errors are purely statistical. A
n naive extrapolation towards the continuum limit yields Γ (cid:39) 1.122(14)KeV which is smaller than
u the previous quenched result and most of the current experimental results. Possible reasons are
J discussed.
8
2
I. INTRODUCTION and then to detect the real γγ pairs. Improvements of
] the measurement for two-photon branching fraction of
at Charmonium systems play a major role in the under- charmonia will soon be reached in the future.
l standing of the foundation of quantum chromodynamics On the theoretical side, charmonium electromagnetic
-
p (QCD), the fundamental theory for the strong interac- transitions have been investigated using various theoret-
e tion. Duetoitsintermediateenergyscaleandthespecial ical methods [6–14]. In principle, these processes involve
h
featuresofQCD,bothperturbativeandnon-perturbative both electromagnetic and strong interactions, the for-
[
physics show up within charmonium physics, making it mer being perturbative in nature while the latter being
4 an ideal testing ground for our understanding of QCD non-perturbative. Therefore, the study for charmonium
v from both sides. transitions requires non-perturbative theoretical meth-
6
Two-photon decay width of η has been attracting ods such as lattice QCD. Normal hadronic matrix ele-
7 c
considerable attention over the years from both theory ment computations are standard in lattice QCD, how-
0
0 and experiment sides. For example, it is related to the ever, processes involving initial or final photons are a bit
0 process gg → ηc relevant for charmonia production at moresubtle. SincephotonsarenotQCDeigenstates,one
. Large Hadron Collider (LHC) and the small-x gluon dis- hastorelyonperturbativemethodsto“replace”thepho-
2
tribution function from the inclusive production of η tonstatesbythecorrespondingelectromagneticcurrents
0 c
6 which describes the non-leptonic B mesons decays [1]. that they couple to. The details of this idea was illus-
1 Furthermore, two-photonbranchingfractionforcharmo- trated in Ref. [15, 16]. Using this technique, the first ab
: niumprovidesaprobeforthestrongcouplingconstantat initioquenchedlatticecalculationoftwophotondecayof
v
the charmonium scale via the two-photon decay widths, charmonia was reported in Ref. [17]. They found a rea-
i
X whichcanbeutilizedasasensitivetestforthecorrections sonable agreement with the experimental world-average
r for the non-relativistic approximation in the quark mod- values for ηc and χc0 decay rates. However, an un-
a
els or the effective field theories such as non-relativistic quenched lattice study is still lacking. In this paper, we
QCD (NRQCD). wouldliketofillthisgapbyexploringthetwophotonde-
On the experimental side, considerable progress has cayratesofηc mesoninlatticeQCDwithNf =2flavors
been made in recent years in the physics of charmonia of light quarks in the sea. The gauge configurations uti-
via the investigations from Belle, BaBar, CLEO-c and lizedinthisstudyaregeneratedbytheEuropeanTwisted
BES [2–5]. Two methods can be utilized to measure Mass Collaboration (ETMC) [18–29], where the twisted
thetwo-photonbranchingfractionforcharmonium. One mass fermion parameters are set at the maximal twist.
is reconstructing the charmonium in light hadrons with This ensures the so-called automatic O(a) improvement
two-photon fusion at e+e− machines. The other one is foron-shellobservableswhereaisthelatticespacing[30].
to make pp¯pairs annihilated to charmonium with decay This paper is organized as follows. In Section II, we
brieflyreviewthecalculationstrategiesforthematrixel-
ement for two-photon decay of η . The matrix element
c
is normally parameterized using a form factor, which in
∗ Correspondingauthor. Email: [email protected] turnisdirectlyrelatedtothedoublephotondecayrates.
2
ThefollowingsectionIIIaredividedintothreepartscon- culation in the future and conclude.
taining details of the simulation: In section IIIA, we in-
troduce the twist mass fermion formulation and give the
parameters of the lattices used in our simulation. In sec-
II. STRATEGIES FOR THE COMPUTATION
tion IIIC, the continuum and lattice dispersion relations
for η are checked. In section IIIE, numerical results of
c
the form factor are presented which are then converted In this section, we briefly recapitulate the methods for
tothedecaywidthofη meson. Ourfinalnumbercomes the calculation of two-photon decay rate of η presented
c c
out to be smaller than the world-average experimental in Ref. [17]. The amplitude for two-photon decay of η
c
result and barely agrees with the previous quenched re- canbeexpressedintermsofaphotontwo-pointfunction
sult. Possible reasons are discussed for this discrepancy. inMinkowskispacebymeansoftheLehmann-Symanzik-
In Section IV, we discuss possible extensions of this cal- Zimmermann (LSZ) reduction formula,
(cid:90)
(cid:104)γ(q1,λ1)γ(q2,λ2)|ηc(p)(cid:105) = − lim (cid:15)∗µ(q1,λ1)(cid:15)∗ν(q2,λ2)q1(cid:48)2q2(cid:48)2 d4xd4yeiq1(cid:48).y+iq2(cid:48).x(cid:104)Ω|T(cid:8)Aµ(y)Aν(x)(cid:9)|ηc(pf)(cid:105). (1)
q1(cid:48)→q1
q2(cid:48)→q2
Here |Ω(cid:105) designates the QCD vacuum state, |η (p)(cid:105) is ity, respectively. Then one utilizes the perturbative na-
c
the state with an η meson of four-momentum p and tureofthephoton-quarkcouplingtoapproximatelyinte-
c
|γ(q ,λ )(cid:105) for i=1,2 denotes a single photon state with grate out the photon fields and rewrites the correspond-
i i
corresponding polarization vector (cid:15)(q ,λ ), with q and ing path-integral as,
i i i
λ being the corresponding four-momentum and helic-
i
(cid:90) (cid:90) (cid:90)
DADψ¯DψeiSQED[A,ψ¯,ψ]Aµ(y)Aν(x)= DADψ¯DψeiS0[A,ψ¯,ψ](cid:0)...+ e2 d4zd4w
2
×(cid:2)ψ¯(z)γρψ(z)A (z)(cid:3) (cid:2)ψ¯(w)γσψ(w)A (w)(cid:3)+...(cid:1)Aµ(y)Aν(x). (2)
ρ σ
Theintegrationoverthephotonfieldscanbecarriedout ing the disconnected diagrams, one arrives at the follow-
byWickcontractingthefieldsintopropagators. Neglect- ing equation,
(cid:90)
(cid:104)γ(q1,λ1)γ(q2,λ2)|ηc(p)(cid:105)=(−e2) lim (cid:15)∗µ(q1,λ1)(cid:15)∗ν(q2,λ2)q1(cid:48)2q2(cid:48)2 d4xd4yd4w d4zeiq1(cid:48).y+iq2(cid:48).xDµρ(y,z)Dνσ(x,w)
q1(cid:48)→q1
q2(cid:48)→q2
(cid:8) (cid:9)
×(cid:104)Ω|T j (z)j (w) |η (p )(cid:105). (3)
ρ σ c f
(cid:8) (cid:9)
In this equation, form (cid:104)Ω|T j (z)j (w) |η (p )(cid:105). This quantity is non-
ρ σ c f
perturbative in nature and should be computed using
(cid:90)
Dµν(y,z)=−igµν d4k e−ik.(y−z) , (4) lattice QCD methods.
(2π)4 k2+i(cid:15)
The current operators such as j (x) appearing in
ρ
isthefreephotonpropagator,whichinmomentumspace Eq. (3) are electromagnetic current operators due to all
will cancel out the inverse propagators outside the in- flavors of quarks. However, we will only consider the
tegral in Eq. (1) and Eq. (3) when the limit is taken. charm quark in this preliminary study. Contributions
Effectively, each initial/final photon state in the prob- duetootherquarkflavors,e.g. up,downorstrange,only
lem is replaced by a corresponding electromagnetic cur- come in via disconnected diagrams which are neglected
rent operator which couples to the photon and eventu- in this exploratory study. Another subtlety in the lat-
ally one needs to compute a three-point function of the tice computation is that, with c(x)/c¯(x) being the bare
3
charm/anti-charm quark field on the lattice, composite The resulting expression (3) can then be analytically
operators such as the current j (x) = Z (g2)c¯(x)γ c(x) continuedfromMinkowskitoEuclideanspace. Thiscon-
ρ V 0 ρ
needs an extra multiplicative renormalization factor Z tinuationworksaslongasnoneoftheq2 istootime-like.
V i
which we infer from Ref. [31]. To be specific, for the two Tobeprecise,thecontinuationisfineaslongasthevirtu-
set of lattices used in this study, the values of the renor- alitiesofthetwophotonsQ2 ≡(−q2)>−M2 whereM
i i V V
malization factor Z (g2) are 0.6103(3) and 0.6451(3) is the mass of the lightest vector meson in QCD [15, 17].
V 0
for the lattice size 243 × 48 at β = 3.9 and 323 × 64 For quenched lattice QCD, the lightest vector meson is
at β = 4.05, respectively. Annihilation diagrams of J/ψ. However, for our unquenched study, it is safe to
the charm quark itself are also neglected due to OZI- take M = m , i.e. the mass of the ρ meson. Using
V ρ
suppression. In fact, in our twisted mass lattice setup, suitable interpolating operator (denoted by O (x)) to
ηc
we introduce two different charm quark fields with de- create an η meson from the vacuum and reversing the
c
generate masses, so that this type of diagram is absent, operator time-ordering for later convenience, we finally
see subsection IIIA. obtain,
(cid:15) (q ,λ )(cid:15) (q ,λ ) (cid:90)
(cid:104)η (p )|γ(q ,λ )γ(q ,λ )(cid:105)= lim e2 µ 1 1 ν 2 2 dt e−ω1|ti−t|
c f 1 1 2 2 tf−t→∞ 2ZEηηcc((ppff))e−Eηc(pf)(tf−t) i
(cid:28)Ω(cid:12)(cid:12)(cid:12)(cid:12)T(cid:110)(cid:90) d3xe−ipf·xOηc(x,tf)(cid:90) d3yeiq2·yjν(y,t)jµ(0,ti)(cid:111)(cid:12)(cid:12)(cid:12)(cid:12)Ω(cid:29), (5)
where O (x) is an interpolating operator that will cre- (cid:104)Ω|O (x,t )jν(y,t)jµ(0,t )|Ω(cid:105). Of course, one has to
ηc ηc f i
ate an η meson from the vacuum and ω is the energy compute the above three-point functions for each t and
c 1 i
of the first photon. The kinematics in this equation is perform an integration (summation) over t .
i
such that four-momentum conservation p = q +q is Apartfromtheabovementionedthree-pointfunctions,
f 1 2
valid. This equation serves as the starting point for our wealsoneedinformationfromη two-pointfunction. For
c
subsequent lattice computation. Basically, the current example, in the above equation, Z (p ) is the spectral
ηc f
that couples to the first photon is placed at the source weightfactorwhileE (p )istheenergyforη withfour-
ηc f c
time-slice t , the second current is at t while the final momentum p = (E ,p ). These can be inferred from
i f ηc f
η meson is at the sink time-slice t and we are led to the corresponding two-point functions for η . For this
c f c
the computation of a three-point function of the form purpose, two-point correlation functions for the interpo-
lating operators O are computed in the simulation:
ηc
C(pf;t)≡(cid:88)e−ipf·x(cid:104)Ω|Oηc(x,t)Oη†c(0,0)|Ω(cid:105)−t(cid:29)→1 |ZEηc((ppf))|2e−Eηc(pf)·T2 cosh(cid:20)Eηc(pf)·(cid:18)T2 −t(cid:19)(cid:21) , (6)
x ηc f
where Z (p ) = (cid:104)Ω|O |η (p )(cid:105) is the corresponding The three-point functions, denoted by G (t ,t), that
ηc f ηc c f µν i
overlap matrix element. need to be computed in our simulation are of the form,
Gµν(ti,t)=(cid:28)Ω(cid:12)(cid:12)(cid:12)(cid:12)T(cid:110)(cid:90) d3xe−ipf·xOηc(x,tf)(cid:90) d3yeiq2·yjν(y,t)jµ(0,ti)(cid:111)(cid:12)(cid:12)(cid:12)(cid:12)Ω(cid:29). (7)
Keeping the sink of η fixed at t = T/2, we compute desired matrix element is obtained by using the results
c f
G (t ,t) across the temporal direction for all t and t of G (t ,t) for different combinations of t and t and
µν i i µν i i
on our lattices. For a fixed ti, one has to use sequen- integrate over ti with an exponential weight e−ω1|ti−t|.
tial source technique to obtain the t dependence of the In practice, the integral is replaced by a summation over
three-point function. Then, the same calculation is re- t . To explore the validity of this replacement, we have
i
peatedwithavaryingti. Then,accordingtoEq.(5),the checked the behavior of the integrand some of which are
4
illustrated in Fig. 1. It is seen that these integrand as One can perform a chiral twist to transform the quark
a function of t indeed peak around the corresponding t fields in physical basis to the so-called twisted basis as
i
values. follows:
The matrix element (cid:104)ηc|γ(q1,λ1)γ(q2,λ2)(cid:105) can be pa- (cid:16)u(cid:17) (cid:18)χ (cid:19)
rameterized using the form factor F(Q2,Q2) as follows, =exp(iωγ τ /2) u
1 2 d 5 3 χ
d
(cid:104)ηc|γ(q1,λ1)γ(q2,λ2)(cid:105)=2(23e)2m−ηc1F(Q21,Q22)(cid:15)µνρσ (cid:16)s(cid:17)=exp(iωγ τ /2)(cid:18)χs (cid:19)
×(cid:15)µ(q1,λ1)(cid:15)ν(q2,λ2)q1ρq2σ, (8) s(cid:48) 5 3 χs(cid:48)
(cid:16)c(cid:17) (cid:18)χ (cid:19)
where(cid:15)µ(q ,λ ),(cid:15)ν(q ,λ )arethepolarizationvectorsof =exp(iωγ τ /2) c (11)
1 1 2 2 c(cid:48) 5 3 χ
the photons while q and q are the corresponding four- c(cid:48)
1 2
momenta. The physical on-shell decay width Γ for η to
c where ω =π/2 implements the full twist.
twophotonsisrelatedtotheformfactoratQ2 =Q2 =0,
1 2 Two sets of gauge field ensembles are utilized in this
which will be referred to as the physical point in the
work,eachcontaining200gaugefieldconfigurations. We
following, via
shall call them Ensemble I and II respectively. The ex-
(cid:18) (cid:19) plicit parameters are listed in Table I. The correspond-
16
Γ=πα2 m |F(0,0)|2, (9) ingrenormalizationfactorZ (g2)andthevalencecharm
em 81 ηc V 0
quark mass parameter µ are taken from Ref. [31].
c
where α (cid:39) (1/137) is the fine structure constant.
em
Therefore, to extract the physical decay width, we sim-
TABLE I. Parameters for the gauge ensembles used in this
ply compute the corresponding three-point functions in
work. See Ref. [31] and references therein for notations.
Eq. (5) and then extract the form factors F(Q2,Q2) at
1 1 Ensemble β a[fm] V/a4 aµ m [MeV] aµ Z (g2)
various virtualities close to the physical point. Then, we sea π c V 0
I 3.9 0.085243×480.004 315 0.2150.6103(3)
canextracttheinformationforF(0,0)yieldingthephys-
II 4.050.067323×640.003 300 0.1850.6451(3)
ical decay width. Although the physical decay width
is only related to F(0,0), the behavior of F(Q2,q2) at
1 2 For the meson operators, in the physical basis, we use
non-zero virtualities are also of physical relevance when
simplequarkbi-linearssuchasq¯Γqandthecorresponding
studying processes involving one or two virtual photons.
formintwistedbasiswillbedenotedasχ¯ Γ(cid:48)χ whichcan
q q
bereadilyobtainedfromEq.(11). Forlaterconvenience,
these are tabulated in table II together with the possible
III. SIMULATION DETAILS
JPC quantum numbers in the continuum and the names
of the corresponding particle in the light and the charm
A. Simulation setup
sector. The current operators that appear in Eq. (7) are
also listed.
Inthisstudy,weusetwistedmassfermionsatthemax-
imaltwist. Themostimportantadvantageofthissetupis
TABLE II. Local interpolating operators for vector and
theso-calledautomaticO(a)improvementforthephysi-
pseudo-scalar states and the current operators that appear
calquantities. Tobespecific,weuseN =2(degenerate
f in Eq. (7) in both physical and twisted basis, q¯Γq =χ¯ Γ(cid:48)χ .
u and d quark) twisted mass gauge field configurations q q
ThenamesofthecorrespondingparticleandtheirJPC quan-
generated by the European Twisted Mass Collaboration
tum numbers in the continuumare also listed. The index for
(ETMC). The other quark flavors, namely strange and i, µ and ν are 1,2,3.
charm quarks, are quenched. These quenched flavors are ρ/J/ψ π/η jµ jν
c
introducedasvalencequarksusingtheOsterwalder-Seiler
Γ γ γ γ γ
(OS)typeaction[18,32]. FollowingtheRefs.[18,22,23], i 5 µ ν
inthevalencesectorweintroducethreetwisteddoublets, Γ(cid:48) γi 1 γµ γν
(u,d),(s,s(cid:48))and(c,c(cid:48))withmassesµ ,µ andµ ,respec-
l s c JPC 1−− 0−+ 1−− 1−−
tively. Within each doublet, the two valence quarks are
regularized in the physical basis with Wilson parameters
of opposite signs (r = −r(cid:48) = 1). The fermion action for
the valence sector reads
(cid:18)χ (cid:19) B. Twisted boundary conditions
S =(χ¯ ,χ¯ )(D +m +iµ γ τ ) u
u d W crit l 5 3 χ
d
(cid:18) (cid:19) Inordertoincreasetheresolutioninmomentumspace,
χ
+(χ¯ ,χ¯ )(D +m +iµ γ τ ) s particularly close to the physical point of Q2 = Q2 =
s s(cid:48) W crit s 5 3 χ 1 2
s(cid:48) 0, it is customary to implement the twisted boundary
(cid:18) (cid:19)
χ conditions(TBC)[31,33–35]inrecentlatticeformfactor
+(χ¯ ,χ¯ )(D +m +iµ γ τ ) c . (10)
c c(cid:48) W crit c 5 3 χ computations, see e.g. [36]. We have also adopted the
c(cid:48)
5
L=24pf =(000) L=32pf =(000)
t/a=4
t/a=4
0.09 t/a=8 0.09 t/a=8
0.08 tt//aa==1126 0.08 tt//aa==1126
0.07 t/a=20 0.07 t/a=20
t/a=24
a) 0.06 a) 0.06 t/a=28
/ /
t 0.05 t 0.05
a, a,
/ /
(ti 0.04 (ti 0.04
Gµν 0.03 Gµν 0.03
0.02 0.02
0.01 0.01
0
0
−0.01 −0.010 5 10 15 20 25 30
0 5 10 15 20
ti/a ti/a
n =(0−1−2); n =(0 0 0)
2 f
n =(0−1−2); n =(0 0 0)
2 f
lattice size: 323×64
lattice size: 243×48
FIG. 1. The integrand in Eq. (5) versus t for various insertion points t obtained from our simulation with ensemble I (left
i
panel), and ensemble II (right panel). We take n = (0,−1,−2); n = (0,0,0) in this example. The insertion points are
2 f
t=4, 8, 12, 16, 20 and t=4, 8, 12, 16, 20, 24, 28 for ensemble I and II, respectively.
twistedboundaryconditionsforthevalencequarkfields, Theallowedmomentaonthelatticearethusmodifiedto
also known as partially twisted boundary conditions.
(cid:18) (cid:19)(cid:18) (cid:19)
2π θ
The quark field ψθ(x,t), when it is transported by an q = n + , (16)
amountofLalongthespatialdirectioni(i=1,2,3),will i L i 2π
change by a phase factor eiθi,
for i=1,2 where n ∈Z3 is a three-dimensional integer.
i
ψθ(x+Lei,t)=eiθiψθ(x,t), (12) By choosing different values for θ, we could obtain more
values of q and q than conventional periodic bound-
where θ = (θ ,θ ,θ ) is the twisted angle for the quark 1 2
1 2 3
ary conditions. In this paper, apart from the untwisted
field in spatial directions which can be tuned freely. In
caseofθ =(0,0,0),wehavealsocomputedthefollowing
this calculations, we only twist one of the charm quark
cases: θ =(0,0,π), (0,0,π/2), (0,0,π/4) and (0,0,π/8).
fieldinbothvectorcurrents,theothercharmquarkfields
Thesechoicesofferusmanymoredatapointsinthevicin-
remain un-twisted. If we introduce the new quark fields
ity of the physical kinematic region.
cˆ(cid:48)(x,t)=e−iθ·x/Lc(cid:48) (x,t), (13)
θ
it is easy to verify that cˆ(cid:48)(x,t) satisfy the conventional C. Meson spectrum and the dispersion relations
periodicboundaryconditionsalongallspatialdirections;
i.e, cˆ(cid:48)(x+Le ,t) = cˆ(cid:48)(x,t) with i = 1,2,3 if the orig-
i Beforecalculatingthematrixelementwithtwophoton
inal field c(cid:48) (x,t) satisfies the twisted boundary condi-
θ decay from η , the mass for η and ρ state and the en-
c c
tions (12). For Wilson-type fermions, this transforma- ergy dispersion relations for η must be verified. This is
c
tion is equivalent to the replacement of the gauge link; particularlyimportantforourstudyduetothefollowing
i.e, reasons. Firstly, we use the N = 2 twisted mass con-
f
U (x)⇒Uˆ (x)=eiθµa/LU (x), (14) figurations, the sea quarks contains u and d quark field.
µ µ µ
therefore virtual ρ state can enter the game. Thus, we
for µ = 0,1,2,3 and θµ = (0,θ). In other words, each should calculate the ρ mass so as to ensure the photon
spatialgaugelinkismodifiedbyaU(1)-phase. Thenthe virtualities Q2,Q2 >−m2 in this simulations. Secondly,
1 2 ρ
current vectors that appear in Eq. (7) are constructed wedoneedtheinformationfromη correlationfunctions,
c
using the hatted and the original charm quark field as, thevalueofE (p)andZ (p)inordertoextracttherel-
ηc ηc
(cid:40) jν(y,t)=c¯(y,t)(γ )cˆ(cid:48)(y,t), evant matrix elements. Finally, we should also check the
ν (15) dispersion relation of ηc which is quite heavy in lattice
jµ(0,t )=cˆ(cid:48)(0,t )(γ )c(0,t ). units(around0.95)inoursimulationandthereforesome
i i µ i
6
The Lattice Dispersion Relation of Emsemble I and II The Continuum Dispersion Relation of Emsemble I and II
2
2.2 Emsemble I
Emsemble I
Emsemble II
Emsemble II
Z=1.049(15)
2 I 1.8 ZI=0.819(12)
Z =0.938(28)
II Z =0.796(24)
II
1.8
1.6
)2
=
E1.6
(
2hn 2E1.4
is41.4
1.2
1.2
1 1
0.8
0 0.05 0.1 0.15 0.2 0.25 0.3 0.8
3 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
4 sin2(pi=2) p2
i=1
P
FIG.2. η mesondispersionrelationobtainedfromsimulationsonEnsembleI(crosses)andEnsembleII(opencircles). Inthe
c
left/rightpanel,thehorizontalaxisrepresentsthelattice/continuumthree-momentumsquaredvariable. Straightlinesinboth
panelsarethecorrespondinglinearfitsusinglattice/continuumdispersionrelationsEq.(18)orEq.(19). Thefittedparameters
Z and Z together with their errors are also shown.
I II
kinematic factors (qρ and qσ) that enter Eq.(8) might and its lattice counterpart,
1 2
need modifications accordingly.
FollowingEq.(6),theenergyEηc(pf)forηc statewith 4sinh2 E(p) =4sinh2 m+Z ·4(cid:88)sin2(cid:16)pi(cid:17) . (19)
three-momentum pf can be obtained from the corre- 2 2 latt 2
sponding two-point function via i
cosh(E (p ))= C(pf;t−1)+C(pf;t+1) . (17) Forfreeparticles,theconstantsZcont andZlatt shouldbe
ηc f 2C(p ;t) closetounity. InFig.2,weshowthiscomparisonforthe
f
two dispersion relations of the η states in our simula-
The two point function is symmetric about t = T/2. c
tion. In the left/right panel, the dispersion relations are
In real simulation we average the data from two halves
illustrated using lattice/continuum dispersion relations,
about t = T/2 to improve statistics. We use the effec-
respectively. In both panels, points with errors are from
tive mass plateaus at zero three-momentum for the η
c simulations on 323×64 (open circles) or 243×48 (stars)
and ρ state to obtain the masses which are then listed
lattices. Straight lines are the corresponding linear fits
in Table III. The mass of the η comes out to be lighter
c to the data. It is seen that, although both dispersion re-
than its physical value since these values are still finite
lations can be fitted nicely using linear fits, the slope for
lattice spacing values. When extrapolated towards the
thenaivecontinuumdispersionrelation, i.e. Z isdef-
continuum limit, the mass will become compatible with cont
initely different from unity, see e.g. right panel of Fig. 2,
the experimental value. The mass of the ρ here serves
while its lattice counterpart Z is close. This suggests
to restrict our kinematic regions where analytic continu- latt
that, for the η state, we should use the lattice disper-
ation is justified. c
sion relations instead of the naive continuum dispersion
relation. This is not surprising since η is quite heavy in
c
TABLE III. The meson mass values for η and ρ obtained latticeunits. Thismodificationofthedispersionrelation
c
from the two ensembles in this work. doeshaveconsequencesonourdeterminationoftheform
Ensemble mηc[MeV] mρ[MeV] factor.
I 2678(3) 903(88)
To illustrate this difference further, we plot the quan-
II 2812(2) 1051(50) tity E2(p)/(m2+p2) as a function of p2 in lattice units
(i.e. a2p2, note that at these small values of p2, the
Similarly, we obtain the energies for η at non-
c difference between p2 and the lattice version pˆ2 is negli-
vanishing momenta via Eq. (17) which then can be uti-
gible) for two of our ensembles. This is shown in Fig. 3.
lized to verify the the following two dispersion relations:
It is seen that this quantity deviates from unity by as
the conventional one in the continuum,
much as 10% even at rather small values of p2. This is
E2(p)=m2+Z ·(cid:88)p2 , (18) actually caused by the difference between the rest mass
cont i
and the kinetic mass of the η meson.
i c
7
1.5
Emsemble I
Emsemble II
1.4
1.3
)
2p 1.2
+
2m
(
= 1.1
2
E
1
0.9
0.8
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
p2
FIG. 3. The quantity E2(p)/(m2+p2) is plotted vs. p2 in lattice units for the η meson for two of our ensembles. It is seen
c
that it deviates from unity at rather small values of p2 which is caused by the difference between the rest mass and kinetic
mass of the meson.
D. Kinematics straint Q2 > −m2. 1 Since p = q +q , this means
1 ρ f 1 2
that, for a given p , a choice of q completely specifies
f 1
q and vice versa. We therefore take several choices of
2
In order to fully explore the form factor close to the q = n (2π/L) by changing three-dimensional integer
1 1
physical point Q2 = Q2 = 0, we performed a parameter n . At this stage, we can compute the energy of the
1 2 1
scan in the two virtualities. The following notations will first photon ω , since ω2 = q2 −Q2. It turns out that
1 1 1 1
be utilized. First of all, in the continuum, we will use wecanalsocomputethevirtualityofthesecondphoton,
q todesignatethefour-momentumofthetwophotons. Q2 = |q |2 − ω2, since ω = E − ω and q is also
1,2 2 2 2 2 ηc 1 2
We will also use ω to denote the temporal component known by the choice of q . One has to make sure that
1,2 1
of q , i.e. ω ≡ q0 . When the photons are on-shell, the values of Q2 thus computed do satisfy the constraint
1,2 1,2 1,2 2
we have ω = |q| with q being the corresponding Q2 > −m2 otherwise it is omitted. This procedure is
1,2 1,2 1,2 2 ρ
three-momentum. The so-called virtuality of the pho- summarized as follows:
tons are defined as the corresponding four-momentum
1. Pick p and θ. Obtain E (p ) from dispersion
squared: Q2 ≡(−q2 ). f ηc f
1,2 1,2 relation (19);
On the lattice, however, there are also lattice coun-
2. JudiciouslychooseseveralvaluesofQ2 inasuitable
terparts of the above notations, arising from the lattice 1
range, say Q2 ∈[−0.5,+0.5]GeV2;
dispersion relation (19). For that we simply add a hat 1
on the corresponding variable. For example, we will use 3. Pick values of n such that q = n (2π/L). This
1 1 1
ωˆ =2sinh(ω /2) to denote the lattice version of ω . fixesbothω andQ2,usingenergy-momentumcon-
1 1 1 1 2
servation;
The computation has to cover the physical interesting
kinematic region. For this purpose, we have to scan the 4. Make sure all values of Q2,Q2 > −m2, otherwise
1 2 ρ
corresponding parameter space. We basically follow the the choice is simply ignored;
followingstrategy: Wefirstfixthefour-momentumofη ,
c 5. Foreachvalidatedchoiceabove,computethethree-
p =(E ,p ),andplaceitonagiventime-slicet =T.
f ηc f f pointfunctions(7),thetwo-pointfunctions(6)and
Note that we just have to fix p = n (2π/L) and E
f f ηc eventually obtain the hadronic matrix element us-
can be obtained from the dispersion relation (19). This
ing Eq. (5).
effectively puts η on-shell. Here we also have the free-
c
dom to pick a value for the twist angle θ. Then, we
judiciously choose several values of virtuality Q2 around
1
the physical point Q2 =0. To be specific, we picked the 1 Thisisvalidwiththephysicalρmesonmass. Ourlatticevalues
1
range Q2 ∈ [−0.5,+0.5]GeV2, which satisfies the con- yieldalessstringentconstraint.
1
8
E. Form factors using the twisted boundary conditions together with dif-
ferent combinations of the lattice momenta, we are able
to populate the physical region close to Q2 = Q2 = 0
In order to compute the desired hadronic matrix ele- 1 2
rather effectively. We have tried both correlated and un-
ment (cid:104)η (p )|γ(q ,λ )γ(q ,λ )(cid:105) in Eq. (5), we choose to
c f 1 1 2 2
correlated fits on our data. The central values for the
place η state at a fixed sink position t = T/2. This
c f
fitted parameters are compatible, however, the error es-
sink position is then used as a sequential source for a
timatesaresomewhatdifferent. Weadoptthecorrelated
backward charm propagator inversion. We compute this
fits as our final results. Fits for other set of parame-
with all possible source positions t and insertion point
i
ters are similar and the final results are summarized in
t. This method allows us to freely vary the value of ω ,
1
Q2 (as discussed in previous subsection) and to directly Table IV for reference.
1 Having obtained the results for F(Q2,0), we can fit it
inspect the behavior of the integrand in Eq. (5). 1
again with another one-pole form,
Taking p = 0 for the η state as an example,, we
f c
show the behavior of the integrand in Fig. 1 for in- F(Q2,0)=F(0,0)/(1+Q2/ν2) (21)
sertion positions t = 4,8,12,16,20 for ensemble I and 1 1
t = 4,8,12,16,20,24,28 for ensemble II. It is seen that with F(0,0) and ν2 being the fitting parameters. This
the integrand is peaked around t = t, making the con- is illustrated in Fig. 6 for two of our ensembles. Again,
i
tributions close to this point the dominant part of the correlated fits are adopted here.
matrixelement. Forthelatticetheory, theintegrationof Apart from fitting the data in a two-step procedure as
t in Eq. (5) is replaced by a summation. described above, we have also tried to fit the data in a
i
When passing from the matrix element to the form one-step method. When we plug Eq. 21 into Eq. 20 and
factors, one should be careful about the form of the mo- assuming that we are only interested in the value of the
menta to use. Recall that these momentum factors orig- form factor close to the physical point, we may Taylor
inate from derivatives in the continuum. On the lattice, expand it assuming both Q2 and Q2 are small,
1 2
they should be replaced by the corresponding finite dif-
F(Q2,Q2)=F(0,0)+aQ2+bQ2 , Q2,Q2 ∼0.(22)
ferences,i.e. oneshouldusethelatticeversionofthemo- 1 2 1 2 1 2
mentum: q0 → 2sinh(q0/2) and qi → 2sin(qi/2). Since Thus,wecouldfitthedatainaregionclosetotheorigin
thespatialmomentathatweareusingarerelativelysmall with F(0,0), a and b being the fitting parameters. This
in lattice units, the effect of this replacement might be is illustrated in Fig. 7 for two of our ensembles. In each
optional. However, forthe0-thcomponent, sinceeachof case, 35 data points of (Q2,Q2) close to the origin are
1 2
thephotonisroughlyhalfoftheηc energywhichislarge taken and the corresponding form factors F(Q21,Q22) are
in lattice units as we discussed in subsection IIIC, this obtained. Then using a linear fit in both Q2 and Q2, c.f.
1 2
replacement does make a difference. Eq. (22), the form factors at the origin are obtained for
AccordingtoEq.(5),thematrixelementandtherefore both ensembles. Again, correlated fits are adopted here.
also the form factor F(Q2,Q2) should be independent of The fitting results are summarized in Table VI.
1 2
the insertion point t. We indeed observe this plateau be- When computing the physical double photon decay
havior in our data which is illustrated in Fig. 4 for the width, according to Eq. (9), one has to plug in the mass
case of Q2 = 0 as an example. Other cases are simi- oftheη meson. Whatwereallycomputeonthelatticeis
1 c
lar. Fitting these plateaus then yields the corresponding the combination of correlation functions which is related
values for the matrix element (cid:104)ηc|γ(q1,λ1)γ(q2,λ2)(cid:105) or to the matrix element (cid:104)ηc|γ(q1,λ1)γ(q2,λ2)(cid:105) via Eq. (5).
equivalently the form factor F(Q2,Q2). When we parameterize this particular matrix element in
1 2
To describe the virtuality dependence of the form fac- terms of form factor in Eq. (8), the relation involves m
ηc
tor, we adopt a simple one-pole parametrization to fit as well. Therefore, the decay width turns out to be pro-
our data. portionaltom3 : Γ∝m3 |(cid:104)η |γ(q ,λ )γ(q ,λ )(cid:105)|2. Here
ηc ηc c 1 1 2 2
it is then quite different if one substitutes in the value of
F(Q2,Q2)=F(Q2,0)/(1+Q2/µ2(Q2)), (20) m obtained on the lattice, or the true physical value of
1 2 1 2 1 ηc
mphys. = 2.98GeV, the two differs by about 10% for the
whereF(Q2,0)andµ2(Q2)areregardedasthefittingpa- coηacrser lattice and about 5% for the finer lattice. There-
1 1
rametersatthegivenvalueofQ2. Sincemeasurementsat fore, if one would substitute in the true physical mass, it
1
different values of Q2 or Q2 are all obtained on the same will result in a 15% difference in the value of Γ for the
1 2
setofensembles,weadoptthecorrelatedfits,takinginto finer lattice and about 30% for the coarser one.
account possible correlations among different Q2 values. The reason for the above mentioned difference is the
The covariance matrix among them are estimated using following. We are taking the value of the valence charm
a bootstrap method. quark mass parameter µ from Ref. [31]. There, it is as-
c
As an example, taking Q2 = −0.5GeV2, the fitting sumedthat,whenthecontinuumlimitistaken,thevalue
1
results are shown in Fig. 5. It is seen that this sim- of m will recover its physical value. However, being on
ηc
ple formula describes the data rather well even for quite a finite lattice, the computed value of m comes out to
ηc
large values of Q2. We therefore have taken all available belessthanthecorrespondingphysicalvalue. Thediffer-
2
valuesofQ2 intothefittingprocess. Noticealsothat, by ence of the two is in fact an estimate of the finite lattice
2
9
L =24 Q2 = 0(GeV2)
1 L =32 Q2 = 0(GeV2)
0.4 1
0.4
n =(0 −1 −2)
2 n =(0 −1 −2)
2
n =(−1 −1 −2)
2 n =(−1 −1 −2)
0.3 2
n =(−1 −2 −2) 0.3
2 n =(−1 −2 −2)
2
)
22 )
Q 22
Q
,0.2
21 ,0.2
Q 21
Q
F( F(
0.1 0.1
0 0
0 5 10 15 20 0 10 20 30
t/a t/a
Ensemble I Ensemble II
FIG.4. Theplateauoftheformfactorobtainedbyanintegration(summation)overt forthree-pointfunctionG (t ,t)with
i µν i
ensemble I (left panel) and ensemble II (right panel). We take Q2 = 0; n = (0,0,0) in this particular plot. Different data
1 f
points correspond to different choices of n as indicated.
2
L=24 Q2=−0.5(GeV2) L=32 Q2=−0.5(GeV2)
1 1
0.2 0.2
p =0,θ =0 p =0,θ =0
0.18 pf1=1,θ3=0 0.18 pf1=0,θ3=π
f1 3 f1 3
0.16 p =0,θ =π 0.16 p =0,θ =π/2
f1 3 f1 3
0.14 pf1=0,θ3=π/2 0.14 pf1=0,θ3=π/4
p =0,θ =π/4 p =0,θ =π/8
2Q)20.12 pff11=0,θ33=π/8 2Q)20.12 f1 3
2Q,1 0.1 2Q,1 0.1
( (
F0.08 F0.08
0.06 0.06
0.04 0.04
0.02 0.02
0 0
−1 0 1 2 3 4 −1 0 1 2 3
Q22(GeV2) Q22(GeV2)
Ensemble I Ensemble II
FIG. 5. The fitted results for F(Q2,Q2) = F(Q2,0)/(1+Q2/µ2(Q2)) by one-pole form factor for Ensemble I (left figure)
1 2 1 2 1
and Ensemble II (right figure) at a fixed value of Q2 = −0.5GeV2. Different data points correspond to different parameter
1
combinations as indicated in the Figure. p denotes x-component of the momentum p of η , and θ represents z-component
f1 f c 3
of twisted angle θ.
spacingerror. Infact,m isnottheonlyfactorwhichaf- and using Eq. (21), with the values of m obtained
ηc ηc
fectstheresults. TherenormalizationfactorZ (g2)that from each ensemble substituted in, we obtain for the de-
V 0
wequotedinTableIalsodependsonthelatticespacing. cay width Γ = 1.019(3)KeV for the coarser and Γ =
Therefore, we think it is more consistent to substitute in 1.043(3)KeV for the finer lattice ensembles. These re-
the values of m computed on each ensembles. In the sults for the form factor F(0,0) together with the corre-
ηc
end,ofcourse,oneshouldtrytotakethecontinuumlimit sponding results for the decay width are summarized in
when the lattice computations are performed on a set of Table V.
ensembles with different lattice spacings. As for the one-step fitting procedure using Eq. (22),
If using the two-step fitting procedure using Eq. (20) we obtain Γ = 1.025(5)KeV for the coarser and Γ =
10
TABLE IV. The summary of fitted results for F(Q2,0) and µ2(Q2) using Eq. (20) for Ensemble I (left four columns) and
1 1
EnsembleII(rightfourcolumn). Thetotalχ2 valueandthecorrespondingtotaldegreesoffreedomisalsolistedinthecolumns
labelled by χ2/dof.
Ensemble I Ensemble II
Q2(GeV2) F(Q2,0) µ2(Q2)(GeV2)χ2/dof Q2(GeV2) F(Q2,0) µ2(Q2)(GeV2)χ2/dof
1 1 1 1 1 1
-0.5 0.11521(46) 7.79(38) 2.22/11 -0.5 0.11297(44) 8.59(42) 0.27/8
-0.4 0.11353(42) 7.82(36) 2.33/11 -0.4 0.11163(47) 8.62(49) 0.24/7
-0.3 0.11187(39) 7.83(36) 2.30/11 -0.3 0.11031(49) 8.61(54) 0.20/6
-0.2 0.11038(41) 7.83(37) 2.12/10 -0.2 0.10901(52) 8.69(53) 0.19/6
-0.1 0.10874(39) 7.87(36) 2.30/10 -0.1 0.10771(54) 8.62(64) 0.15/5
0 0.10721(43) 7.90(45) 1.79/8 0 0.10645(59) 8.67(59) 0.15/5
0.1 0.10581(46) 7.82(52) 1.37/7 0.1 0.10523(67) 8.74(60) 0.13/5
0.2 0.10432(44) 7.85(49) 1.54/7 0.2 0.10402(71) 8.70(73) 0.03/3
0.3 0.10283(44) 7.92(47) 1.43/7 0.3 0.10191(56) 7.74(47) 2.9/3
0.4 0.10142(45) 7.96(44) 1.51/7 0.4 0.10056(58) 7.83(45) 2.8/3
0.5 0.10012(51) 7.82(51) 0.95/5 0.5 0.09936(63) 7.64(55) 2.3/2
L =24
0.2
L =32
0.2
0.15
0.15
)
0
)
, 0
2Q1 0.1 2,1 0.1
( Q
F (
F
0.05 0.05
0
0
−1 −0.5 0 0.5 1
−1 −0.5 0 0.5 1
Q2(GeV2)
Q2(GeV2) 1
1
Ensemble II
Ensemble I
FIG.6. F(Q2,0)isagainfittedwithaone-poleform: F(Q2,0)=F(0,0)/(1+Q2/ν2)forEnsembleI(leftfigure)andEnsemble
1 1 1
II (right figure).
1.062(5)KeV for the finer lattice ensembles. The results width to be very sensitive to the pion mass. Also, since
for the form factor F(0,0) and Γ(η → γγ) are consis- bothofourensembleshavem L∼3.3,wedonotexpect
c π
tent with each other for Ensemble I using two different verylargefinitevolumeerrorsaswell. Sincewehaveonly
types of fitting procedure. However, for Ensemble II, a two ensembles, it is not possible to make reliable extrap-
combined fitting using Eq. (22) gives a larger result for olation towards the continuum limit. However, if one
both F(0,0) and Γ(η → γγ). We think the value using would try a naive continuum limit extrapolation, assum-
c
the combined fit is more reliable since it gives a much inganO(a2)error,weobtainΓ=1.082(10)KeVwhichis
less value of χ2/dof. In the combined fitting, the naive also listed in Table V. There are of course other sources
continuum extrapolated result for the decay width reads of systematic errors, e.g. the neglecting of the so-called
Γ=1.122(14)KeV.ThefittedresultsforF(0,0)together disconnected contributions, the quenching of the strange
with the corresponding results for the decay width are quark, etc. Therefore, we decided not to quantify the
summarized in Table VI. systematic errors in this exploratory study. However, as
wediscussedabove, thedifferenceintheη massalready
c
Let us now discuss the possible systematic errors. Al- indicatesthattheremightbeafinitelatticespacingerror
thoughthemassofthepioninthetwoensemblesarerel- at the order of 15% for the finer and 30% for the coarser
ativelyheavy, wedonotexpectthedoublephotondecay
