Table Of ContentTopological Charge in Two Flavors QCD with
1
Optimal Domain-Wall Fermion
1
0
2
n
a
J
2
] Tung-HanHsieh∗1,Ting-WaiChiu2,3,Yao-YuanMao2 (for the TWQCDCollaboration)
t
a 1 ResearchCenterforAppliedSciences,AcademiaSinica,Taipei115,Taiwan
l
- 2 DepartmentofPhysics,andCenterforTheoreticalSciences,NationalTaiwanUniversity,
p
e Taipei10617,Taiwan
h 3 CenterforQuantumScienceandEngineering,NationalTaiwanUniversity,Taipei10617,
[
Taiwan
1
v
2
0
We measurethe topologicalchargeofthegaugeconfigurationsgeneratedbylattice simulations
4
0 of2flavorsQCDona163×32lattice,withOptimalDomain-WallFermion(ODWF)atNs=16
1. and plaquette gauge action at b =5.90. Using the Adaptive Thick-Restart Lanczos algorithm,
0
we project the low-lying modes of the 4D effective Dirac operator of ODWF, and obtain the
1
1 topologicalchargeandtopologicalsusceptibility. Ourresultoftopologicalsusceptibilityagrees
v: with thesea-quarkmassdependencepredictedby thechiralperturbationtheory,andprovidesa
i determinationofthechiralcondensate.
X
r
a
TheXXVIIIInternationalSymposiumonLatticeFieldTheory,Lattice2010,June14-19,2010,Villasimius,
Italy
∗Speaker.
(cid:13)c Copyrightownedbytheauthor(s)underthetermsoftheCreativeCommonsAttribution-NonCommercial-ShareAlikeLicence. http://pos.sissa.it/
TopologicalChargeinLatticeQCDwithODWF Tung-HanHsieh
1. Introduction
In Quantum Chromodynamics (QCD), the topological susceptibility (c ) is the most essen-
t
tial quantity to measure the topological charge fluctuation of the QCD vacuum, which plays an
important roleinbreaking theU (1)symmetry. Theoretically, c isdefinedas
A t
Q2
c = d4xhr (x)r (0)i= t (1.1)
t Z (cid:10)W (cid:11)
where
1
r (x)= 32p 2e mnls tr[Fmn (x)Fls (x)], Qt ≡Z d4xr (x), (1.2)
Q isthetopological charge (which isaninteger forQCD),r (x)isthetopological charge density,
t
and W is the space-time volume. In the chiral perturbation theory (ChPT), the quark mass depen-
denceofc wasderivedatthetreelevel[1]in1992,andithasbeenextendedtotheone-looporder
t
recently [2]. ForN =2,theseChPTformulascanbewrittenas
f
S
c /m = , (1.3)
t u
1+m /m
u d
S 3M2 M2 m m
c /m = 1− p ln p +K (m +m )+2(2K +K ) u d ,(1.4)
t u 1+mu/md (cid:20) (cid:18)32p 2Fp2(cid:19) m s2ub 6 u d 7 8 mu+md(cid:21)
where S is the chiral condensate, Fp the pion decay constant, and Ki are related to low energy
constants L. A salient feature of c is that it is suppressed in the chiral limit (m →0) due to the
i t u
internal quarkloops. Mostimportantly, (1.3)and(1.4)provide aviable waytoextractS andother
lowenergy constants.
However, one cannot determine Q reliably using the link variables, due to the rather strong
t
fluctuations atshortdistances. Instead, weconsider theAtiyah-Singer indextheorem [3]
Q =n −n =index(D), (1.5)
t + −
where n± isthe number ofzero modes of themassless Dirac operator D ≡gm (¶ m +igAm )with ±
chirality. Thus,bycoutingthezeromodesofthemasslessDiracoperator,wecandetermineQ and
t
c foranensembleofgaugeconfigurations.
t
Recently, the topological susceptibility has been measured in unquenched lattice QCD with
exact chiral symmetry, for N =2 lattice QCD with overlap fermion in a fixed topology [4], and
f
N =2+1 lattice QCD with domain-wall fermion [5]. These results are in good agreement with
f
theChPTatthetreelevel(1.3).
In this work, we measure the topological charge of the gauge configurations generated by
latticesimulationsof2flavorsQCDona163×32lattice, withtheOptimalDomain-WallFermion
(ODWF) [6] at N = 16, and plaquette gauge action at b = 5.90. Mathematically, ODWF is a
s
theoretical framework which can preserve the chiral symmetry optimally for any finite N . Thus
s
the artifacts due to the chiral symmetry breaking with finite N can be reduced to the minimum,
s
especially inthechiralregime. The4Deffectiveoperator ofmasslessODWFis
1−(cid:213) Ns T 1−w H
D=m [1+g S (H )], S (H )= s=1 s, T = s w, (1.6)
0 5 opt w opt w 1+(cid:213) Ns=s1Ts s 1+w sHw
2
TopologicalChargeinLatticeQCDwithODWF Tung-HanHsieh
which is exactly equal to the Zolotarev optimal rational approximation of the overlap Dirac oper-
ator [7]. That is, S (H )=H R (H ), where R (H ) is the optimal rational approximation of
opt w w Z w Z w
(H2)−1/2 [8].
w
Weoutlineourschemeofprojectingthelow-lyingeigenmodesofDasfollows. Using[DD†,g ]=
5
0,andtherepresentation
D|q i=l (q )|q i, l (q )=m (1+eiq ), (1.7)
0
weobtain
1±g
S |q i ≡P H R (H )P |q i =±cosq |q i , P = 5, (1.8)
± ± ± w Z w ± ± ± ± 2
where |q i=P |q i+P |q i=|q i +|q i . Thus, we can perform the eigenmode projection on
+ − + −
theoperator S insteadofD. Moreover, |q i arerelatedtoeachotherthroughtherelation
± ±
1
|q i= (g S−e−iq )|q i , q 6=0,±p ,±2p ,... (1.9)
isinq 5 ±
Assumingthezeromodesonlyresiding ineither+or−chirality, weproject thelow-lying modes
of D as follows. First, we project the smallest eigenmodes of S . If D has zero modes with
+
positive chirality, then the smallest eigenvalues of S are equal to −1, and the corresponding full
+
eigenvectors can be obtained using Eq. (1.9). Otherwise, we have to check whether D has zero
modes with negative chirality by projecting the largest eigenmodes of S . If D has zero modes
−
with negative chirality, the largest eigenvalues of S are equal to +1, and the corresponding full
−
eigenvectors canbeobtained using Eq. (1.9). Finally, fornon-zero low-lying eigenmodes, wecan
pick either +or−chirality ofEq. (1.8)forprojection, then obtain thefull eigenvectors using Eq.
(1.9).
2. The TRLanalgorithm
Ingeneral, thelow-mode projection ofalarge sparse matrixAcan becarried out viaiterative
process, such as the Lanczos algorithm or the Arnoldi algorithm. The common procedure is to
construct anorthonormal basisfromtheKrylovsubspace startingfromaninitialvectorr ,
0
K (A,r )=span r ,Ar ,A2r ,...,Am−1r . (2.1)
0 0 0 0 0
(cid:8) (cid:9)
Thelinearcombinations of{Air ,i=0,1,...m−1}formtheRitzvectorsofA,whichconverge to
0
theeigenvectors ofAasmbecomesverylarge.
SinceH andS areHermitianmatrices, itisnaturaltousetheLanczosalgorithm. Thebasic
w ±
Lanczosprocess(seeFig.1-(a))istoconstruct thefollowingfactorization
AQ =Q T +b q eT, (2.2)
m m m m m+1 m
where{q ,q ,...,q }formanorthonormal basis, Q =[q ,q ,...,q ],e isthem-thcolumnof
1 2 m m 1 2 m m
them×midentitymatrix,andT isatridiagonalmatrixwithdiagonalelementsa andsub-diagonal
m i
elementsb . ThentheRitzpairs(lˆ ,xˆ)canbeobtained from
i i i
Tˆ =U†T U , X =Q U , (2.3)
m m m m m m m
3
TopologicalChargeinLatticeQCDwithODWF Tung-HanHsieh
The update of the T matrix
Input: r ,b =||r ||,q =0
0 0 0 0
Lanczos
For: i=1,2,... m
m
• q =r /b
i i−1 i−1 The 1st loop m k
• p=Aq
i k
• a =q†p
i i After restart
• r = p−a q −b q k diagonalize truncate
i i i i−1 i−1 Lanczos
• b =||r|| k m k
i i
(a) (b)
Figure1: (a)ThebasicLanczosprocess;(b)TheLanczosalgorithmwithThick-Restart.
H : k′=240,m=340 S : k′=100,m=200,n =240
w + p
method N_restarts #ofA·v time(s) speedup N_restarts #ofA·v time(s) speedup
ARPACK 388 35390 136460 1.00 13 1050 112632 1.00
TRLan 999 100140 572951 0.24 12 1300 105790 1.06
n -TRLan 383 59145 78058 1.75 11 1030 90496 1.24
Table1: Thebenchmarkoflow-modeprojectionofH andS foranon-trivialgaugeconfiguration(Q =3)
w + t
at b =5.9 andm a=0.01. Herek′ is the numberof low-modesprojected, m is the dimensionofKrylov
q
subspace,andn isthenumberoflow-modesofH usedforlow-modepreconditioningintheprojectionof
p w
thelow-modesofS .
+
where Tˆ is diagonal with eigenvalues lˆ ,U is a unitary matrix, and X =[xˆ ,xˆ ,...,xˆ ]. As m
m i m m 1 2 m
becomesverylarge, theRitzpairs(lˆ ,xˆ)converge totheeigenmodes (l ,x)ofA.
i i i i
However, the Lanczos algorithm suffers from the following problems. Some specious Ritz
values mayrepeatedly appearwhenmgoestolargervalues. Thisisduetothefactthatthevectors
{q} lose orthogonality rapidly in the finite precision arithmatics. Hence one has to apply the
i
Gram-Schmidt procedure to re-orthogonalize them during the iteration. Also, it requires a large
number ofq toproject afewRitzpairs, consuming alarge amount of memoryandslowing down
i
thecomputation. Awaytoovercometheseproblemsistoperformarestart.
Thereareseveralrestart schemes, inwhichtheThick-Restart Lanczosalgorithm (TRLan)[9]
turns outto bethemost efficient forour purpose. TheTRLanalgorithm isillustrated inFig.1-(b),
where the squares in the diagram stand for the square matrix T , in which the non-zero matrix
m
elementsaredenoted byblacklines.
Obviously, the performance of the TRLan algorithem depends on the truncated dimension
k and the dimension m of the Krylov subspace. Now the question is what is the optimal value
of k (for a given m) such that for each restart, the reduction factor d of the residual of the j-th
j
(not converged) Ritz pair is maximized, while the number of floating-point operations (FPOs) is
minimized. In other words, at each restart, we try to find the optimal k by maximizing the object
function [10]
d lˆ −lˆ
f(k)= j , d ≃C (1+2g ), g = k+1 j, (2.4)
#ofFPOs j m−k lˆ −lˆ
m k+1
4
TopologicalChargeinLatticeQCDwithODWF Tung-HanHsieh
where C (z) is the Chebyshev polynomial of degree m−k. For the number of FPOs, we
m−k
only count the dominant parts which are directly related to the dimension m, namely, the re-
orthogonalization and the update of Ritz vectors. This modified algorithm is called the Adaptive
Thick-Restart Lanczos algorithm, which can attain substantial performance gain, with faster con-
vergence andfewerFPOs.
300 200
150
200
100
100 50
0
0
The # of converged modes -50 The # of converged modes
The dimension k The dimension k
-100 -100
0 100 200 300 400 0 2 4 6 8 10 12
The number of restarts The number of restarts
(a) (b)
Figure2: Thenumberof convergedRitz pairsandthe dimensionk versusthenumberof restartsforthe
n -TRLanalgorithm:(a)H (m=340);(b)S oftheoverlapoperator(m=240).
w +
Forbenchmarking, weprojectthelow-lyingeigenmodesofH andS ,usinganontrivialcon-
w +
figuration (Q =3) generated by the 2-flavor QCD simulation on a 163×32 lattice, with ODWF
t
fermion at N =16 and m a=0.01, and the plaquette gauge action at b =5.90. In Table 1, we
s sea
present our benchmark results forthree different methods: ARPACK(Arnoldi algorithm withIm-
plicit Restart), TRLan (Thick-Restart Lanczos), and n -TRLan (Adaptive Thick-Restart Lanczos).
For TRLan (by definition k =k′, where k′ is the number of low-modes we want to project, and
k is the truncated dimension at each restart), we see that its rate of convergence tends to be very
slow for H , with a rather large condition number (∼106). On the other hand, for n -TRLan, k is
w
automatically adjusted tobeanoptimal valueineachrestart (k′<k<m,seeFig. 2), thustherate
of convergence isdramatically improved. It turns out that, for both H and S , n -TRLanisfaster
w ±
thanTRLanandARPACK.Therefore,weusen -TRLantoprojectthelow-lyingeigenmodesofH
w
andD.
3. Lattice setup
Simulations are carried out for two flavors (N = 2) QCD on a 163×32 lattice at a lattice
f
spacing a∼ 0.11 fm, for eight sea quark masses in the range m a=0.01,0.02,···,0.08 [11]. For
q
thegluonpart,theplaquetteactionisusedatb =5.90. Forthequarkpart,theoptimaldomain-wall
fermionisusedwithN =16. Afterdiscarding300trajectoriesforthermalization, weaccumulated
s
3000−3200 trajectories intotal foreachseaquark mass. Fromthesaturation ofthebinning error
oftheplaquette, aswellastheevolution ofthetopological charge, weestimatetheautocorrelation
timetobearound10trajectories. Thuswesampleoneconfiguration every10trajectories, thenwe
have270−290configurations foreachseaquarkmass.
5
TopologicalChargeinLatticeQCDwithODWF Tung-HanHsieh
Foreachconfiguration,wecalculatethezeromodesand80conjugatepairsofthelowest-lying
eignmodes of the overlap Dirac operator. We outline our procedures as follows. First, we project
240low-lyingeigenmodesofH2 usingn -TRLanalogorithm,whereeacheigenmodehasaresidual
w
less than 10−12. Then we approximate the sign function of the overlap operator by the Zolotarev
optimal rational approximation with 64 poles, where the coefficents are fixed with l 2 =(6.4)2,
max
and l 2 equal to the maximum of the 240 projected eigenvalues of H2. Then the sign function
min w
error is less than 10−14. Using the 240 low-modes of H2 and the Zolotarev approximation with
w
64 poles, we project the zero modes plus 80 conjugate pairs of of the lowest-lying eignmodes of
the overlap operator with the n -TRLanalgorithm, where each eigenmode has a residual less than
10−12.
4. Results
In Fig. 3, we plot the topological susceptibility c = (Q −hQ i)2 /W as a function of the
t t t
sea quark mass m . For all eight sea quark masses, our d(cid:10)ata are well fi(cid:11)tted by a linear function
q
F+Gm withtheinterceptF=0.47(2.37)×10−5 andtheslopeG=5.50(4)×10−3. Evidently,the
q
intercept isconsistent withzero, in agreement withthe c PTexpectation (1.3). Equating theslope
toS /N ,weobtainS =0.0110(7). InordertoconvertS tothatintheMSscheme,wecalculatethe
f
renormalization factor ZMS(2GeV)using the non-perturbative renormalization technique through
m
the RI/MOM scheme [12], and obtain ZMS(2GeV) = 0.5934(10) [13]. Then the value of S is
m
transcribed to
S MS(2GeV)=[261(5)(7) MeV]3, (4.1)
which is in agreement with our previous results : S MS(2GeV)=[245(5)(12) MeV]3 for N =2
f
latticeQCDwithoverlapfermioninafixedtopology[4],andS MS(2GeV)=[253(4)(6) MeV]3for
N =2+1lattice QCDwithdomain-wall fermion[5]. Thestatistical errorrepresents acombined
f
error (including those due to a−1 and ZMS). The systematic error is estimated from the turncation
m
of higher order effects and the uncertainty in the determination of lattice spacing with r =0.49
0
fm. Since our calculation is done at a single lattice spacing, the discretization error cannot be
quantified reliably, but we do not expect much larger error because our lattice action is free from
O(a)discretization effects.
5. Concluding remarks
We have measured the topological charge of the gauge configurations generated by lattice
simulations of2flavors QCDona163×32lattice, withOptimalDomain-Wall Fermion (ODWF)
at N =16 and plaquette gauge action at b =5.90. We use the Adaptive Thick-Restart Lanczos
s
algorithmtocomputethelow-lyingeigenmodesofH2andtheoverlapDiracoperator. Ourresultof
w
c agreeswiththesea-quarkmassdependencepredictedbytheChPTatthetreelevel,andprovides
t
a good determination of the chiral condensate. Furthermore, our result of c also agrees with the
t
theNLOChPT,Eq. (1.4). Wewillpresentthedetailsinaforthcoming paper[14].
6
TopologicalChargeinLatticeQCDwithODWF Tung-HanHsieh
0.0014
N = 2
f
0.0012
0.0010
4V] 0.0008
e
G
χ[t 0.0006
0.0004
0.0002
0.0000
0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16
m [GeV]
q
Figure3: Thetopologicalsusceptibilityoftwo-flavorsQCDwiththeoptimaldomain-wallfermion.
Acknowledgments
This work is supported in part by the National Science Council (Nos. NSC96-2112-M-002-
020-MY3,NSC99-2112-M-002-012-MY3, NSC96-2112-M-001-017-MY3, NSC99-2112-M-001-
014-MY3, NSC99-2119-M-002-001) and NTU-CQSE (Nos. 99R80869, 99R80873). We also
thankNCHCandNTU-CCforproviding facilities toperformsomeofthecomputations.
References
[1] H.LeutwylerandA.Smilga,Phys.Rev.D46,5607(1992).
[2] Y.Y.MaoandT.W.Chiu[TWQCDCollaboration],Phys.Rev.D80,034502(2009)
[3] M.F.AtiyahandI.M.Singer,AnnalsMath.87,484(1968).
[4] S.Aokietal.[JLQCDandTWQCDCollaborations],Phys.Lett.B665,294(2008)
[5] T.W.Chiu,T.H.HsiehandP.K.Tseng[TWQCDCollaboration],Phys.Lett.B671,135(2009)
[6] T.W.Chiu,Phys.Rev.Lett.90,071601(2003);Nucl.Phys.Proc.Suppl.129,135(2004)
[7] H.Neuberger,Phys.Lett.B417,141(1998);R.NarayananandH.Neuberger,Nucl.Phys.B443,305
(1995)
[8] T.W.Chiu,T.H.Hsieh,C.H.HuangandT.R.Huang,Phys.Rev.D66,114502(2002)
[9] K.WuandH.Simon,SIAMJ.MatrixAnal.Appl.22,602,2000.
[10] I.Yamazaki,Z.Bai,H.Simon,L.W.Wang,andK.Wu,Tech.Rep.LBNL-1059E(2008).
[11] T.W.Chiuetal.[TWQCDCollaboration],PoSLAT2009,034(2009)[arXiv:0911.5029[hep-lat]];
PoSLAT2010,099(2010)
[12] G.Martinelli,C.Pittori,C.T.Sachrajda,M.TestaandA.Vladikas,Nucl.Phys.B445,81(1995)
[13] T.W.Chiuetal.[TWQCDCollaboration],inpreparation.
[14] T.W.Chiu,T.H.Hsieh,Y.Y.Mao[TWQCDCollaboration],inpreparation.
7
