Detecting flavor content of the vacuum using the Dirac operator spectrum
Jian Liang, Andrei Alexandru, Yu-Jiang Bi, Terrence Draper, Keh-Fei Liu, Yi-Bo Yang
DDetecting flavors of vacuum from the Dirac operator spectrum
Jian Liang,
1, 2, 3, ∗ Andrei Alexandru, Yu-Jiang Bi, Terrence Draper, Keh-Fei Liu, and Yi-Bo Yang
6, 7, 8, † ( χ QCD Collaboration) Guangdong Provincial Key Laboratory of Nuclear Science, Institute of Quantum Matter,South China Normal University, Guangzhou 510006, China Guangdong-Hong Kong Joint Laboratory of Quantum Matter,Southern Nuclear Science Computing Center, South China Normal University, Guangzhou 510006, China Department of Physics and Astronomy, University of Kentucky, Lexington, KY 40506, USA Department of Physics, The George Washington University, Washington, DC 20052, USA Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China School of Fundamental Physics and Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China International Centre for Theoretical Physics Asia-Pacific, Beijing/Hangzhou, China
From simulations on 2+1 flavor domain wall fermion ensembles at three lattice spacings withtwo of them at physical quark masses, we compute the spectrum of the Dirac operator up to theeigenvalue λ ∼
100 MeV using the overlap fermion. The spectrum is close to a constant below λ ≤
20 MeV as predicted by the 2-flavor chiral perturbative theory ( χ PT) up to the finite volumecorrection, and then increases linearly due to the nonvainishing strange quark mass. Furthermore,one can extract the light and strange quark masses with ∼
20% uncertainties from the spectrumdata with sub-percentage statistical uncertainty, using the next to leading order χ PT. Using thenon-perturbative RI/MOM renormalization, we obtain the chiral condensates at MS 2GeV as Σ =(260 . . . . .
8) MeV) in the N f = 2 (keeping the strange quark mass at the physicalpoint) chiral limit and Σ = (232 . . . . .
8) MeV) in the N f = 3 chiral limit, wherethe four uncertainties come from the statistical fluctuation, renormalization constant, continuumextrapolation and lattice spacing determination. Note that Σ / Σ = 1 . Introduction:
The QCD vacuum includes both quarkand gluon. A gluon interacts with itself directly, buta quark can only interact with other quarks indirectlythrough a gluon. So a fundamental issue is how a quark’sproperties depend upon the flavors of virtual quarks inthe vacuum.The quark Dirac operator D/ = ( ∂ µ − igA µ ) γ µ is anefficacious tool to address this question, as the quark fla-vors in the QCD vacuum (the “sea” quarks) just affectthe D/ though the gluon field A µ . Based on the Banks-Casher relation [1], the near-zero density ρ ( λ ) of D/ is pro-portional to the vacuum chiral condensate (cid:104) ¯ ψψ (cid:105) , whichis the order parameter of spontaneous chiral symmetrybreaking. By treating λ as an imaginary quark mass [2],chiral perturbative theory ( χ PT) reproduces the Banks-Casher relation as its leading order approximation, andprovides a quantitative guide on the next-to-leading or-der (NLO) correction of ρ ( λ ) from the masses of differentquark flavors in the vacuum, in the small quark mass re-gion where χ PT is still valid. But such a χ PT calculationbecomes quite complicated beyond the next-to-leadingorder (NLO) region when the finite volume effect is in-volved, and its predictability is limited as there are many ∗ [email protected] † [email protected] parameters in NNLO.On the other hand, nonperturbative lattice QCD canprovide first-principles information about ρ ( λ ) which hasbeen shown to be renormalizable [8, 9], while it is chal-lenging to be precise and accurate enough to provide theflavor information in the vacuum. First of all, one mustuse the overlap fermion [3] or the projected Domain-wallfermion [4] to obtain ρ ( λ ) accurately at λ ∼ O (10) or more computer time thanthat of the standard Wilson-like discretization of D/ whichbreaks chiral symmetry explicitly. At the same time, oneneeds to solve the smallest eigenpairs (eigenvalue andcorresponding eigenvector) of D/ , count n ( λ ) and extract ρ ( λ ) accurately. This is O (100) more expensive than thecost of a standard quark propagator. Thus the stochasticmethod has been proposed to provide some estimates of ρ ( λ ) [7, 10–12].During investigations in the last decade, the near-zero D/ eigenpairs have been found to be extremely benefi-cial in lattice QCD calculations of the correlators to im-prove the signal-to-noise ratio using low-mode substitu-tion (LMS) for the noise grid source [13–15]. These lowmodes tend to saturate the long-range parts of hadroncorrelators and are useful in calculating quark loops withthe noise estimator only for the high modes [16]. Inthis work, we solve the low-lying eigenpairs of the chi-ral fermion directly, and then present the result of the a r X i v : . [ h e p - l a t ] F e b spectrum D/ with the smallest statistical uncertainty todate on the 2+1 flavor physical quark mass ensembles attwo lattice spacings. Methodology : The overlap fermion was proposed inRef. [17, 18] to construct a discretized fermion actionsatisfying the Ginsburg-Wilson relation D ov γ + γ D ov = aM D ov γ D ov [19], D ov = M (cid:16) γ (cid:15) (cid:0) H w ( − M ) (cid:1)(cid:17) , (1)where (cid:15) ( H w ) = H w √ H is the matrix sign function, H w ( − M ) = γ D w ( − M ) and D w is the Wilson Diracoperator with a negative mass such as M = 1 . /D w . (cid:15) ( H ) can be decomposed into acombination of the small and large eigenvalue regions, (cid:15) ( H w ) = (cid:88) λ H i <λ Hc λ H i | λ H i | v H i ( v H i ) † + (1 − (cid:88) λ H i <λ Hc v H i ( v H i ) † ) N (cid:88) i d i H i +1w , (2)in which the eigenpair ( λ H i , v H i ) satisfies the relation H w v H i = λ H v H i , and d i =0 , , ,. . . . . . ,N are the Chebyshevpolynomial coefficients to approximate (cid:15) ( H w ) for all theeigenvalues larger than a given cutoff λ Hc with a given ac-curacy [20]. The computation cost of D ov is proportionalto the polynomial order N and one can increase λ Hc toreduce N until λ H become very dense for λ H c ∼ . /a .Even for optimal λ Hc we still have N ∼ N makes the computational cost of D ov to be twoorders of magnitude higher than that of a D w .The chiral Dirac operator can be defined through D ov , D c = D ov − M D ov = M γ (cid:15) ( γ D w ( M ))1 − γ (cid:15) ( γ D w ( M )) , (3)and satisfies the same commutation relation D c γ + γ D c = 0 as that of D/ in the continuum. Eacheigenvector of D c is exactly the same as that of D ov ,while their eigenvalues λ c and λ ov satisfy the relation λ c = λ ov − M λ ov . Using the Arnoldi factorization algo-rithm [21], one can obtain the low-lying eigenpairs of D ov ,and thus those of D c . Both D ov and D c are similar to D w (0) at the continuum limit while providing a properultraviolet cutoff to preserve chiral symmetry. Then wecan define the eigenvalue density of D c as ρ ( λ c ) ≡ ∂n ( λ c ) V ∂λ c ,where iλ is an eigenvalue of D c , V = L × T is the 4-D volume, and n ( λ ) is the number of eigenvalues in therange (0 , λ ]. Since ρ ( λ c ) is so far the best approximationof the long distance spectrum of a quark on a discretizedlattice, so that we will just treat λ c as λ in the followingdiscussions.Based on the partially quenched chiral perturbativetheory (PQ χ PT) which can accommodate the valencequark mass being different from the sea quark masses, one can derive the formula to describe (cid:104) ¯ ψψ (cid:105) of the va-lence quark as a function of both valence and sea quarkmasses. In PQ χ PT [2], the density ρ ( λ ) of the low-lyingeigenvalues as the chiral condensate (cid:104) ¯ ψψ (cid:105) with virtualquark mass iλ in a finite volume V = L × T has theexpression, ρ ( λ, V ) = lim (cid:15) → π (cid:0) (cid:104) ¯ ψψ (cid:105)| m v = iλ − (cid:15) − (cid:104) ¯ ψψ (cid:105)| m v = iλ + (cid:15) (cid:1) = Σ π Re[ Z v ( iλ, m seaq ) ˆΣ P Q ( iλV, m q V )] , (4)where Σ ≡ −(cid:104) ¯ ψψ (cid:105)| m q → is the chiral condensate in thechiral limit of N f flavors, and m q is the light or strangesea quark mass. The standard NLO correction Z v =1 + β ( L, T ) N f − N f F √ V + N f − N f λ Σ32 πF + O ( m seaq , λ ) (5)can have leading order finite-volume correction which is ∼
4% with N f = 2 for a (5 . lattice [22], and alsothe λ -dependent finite volume correction,ˆΣ P Q ( iλV, m q V ) = 1 − V ( (cid:88) q = u,d,s iλ + m q + 12 λ (cid:80) q = u,d,s m q ) + O (cid:0) ( m q ) , λ , V (cid:1) , (6)is suppressed at large λ . It is important to note that Z v can differ from 1 by ∼
30% or more if we keep the strangequark mass at the physical point, and it makes the Σ inthe N f =2 and 3 chiral limits quite different from eachother. We skip the complete NLO expression of ρ ( λ, V )since it is very lengthy especially in the N f =2+1 case,and the interested reader can find it in Ref. [2]. TABLE I. Information of the RBC ensembles [23, 24] usedin this calculation. The pion and kaon masses are in unit ofMeV.Symbol L × T a (fm) m π m K Z S N cfg
48I 48 ×
96 0.1141(2) 139 499 1.135(1)(16) 30364I 64 ×
128 0.0837(2) 139 508 1.020(1)(15) 30448IF 48 ×
96 0.0711(3) 234 564 0.986(1)(15) 18532IF 32 ×
64 0.0626(4) 371 558 0.951(2)(14) 50
Numerical setup : We list the four ensembles we used inthis work in Table I, which include two ensembles at phys-ical light and strange quark masses. We will use just thefirst three ensembles to extract ρ ( λ ), and the last ensem-ble with smallest lattice spacing will just be used to showthe lattice spacing dependence of the scalar renormaliza-tion constant. With one step of HYP smearing [25], weneed 800 H w eigenpairs and 1000 pairs of D c eigenpairsto reach the upper bands 0.158 and ± . i respectivelyon the largest 64I ensemble, while the same upper bandscan be reached with only ∼ H w eigenpairs and ∼ Q Z S-1 (MS-bar, 2GeV)0.063 fm0.071 fm0.084 fm0.114 fm
FIG. 1. The inverse of the scalar current renormalization con-stant at MS 2GeV and different lattice spacings, as a functionof RI/MOM off-shell scale Q . The a Q dependencies in thefigure at different lattice spacings are similar, which meansthat the residual Q dependence will vanish in the continuumlimit. pairs of D c eigenpairs on a 24 ×
64 ensemble at 0.1105(2)fm.The bare chiral condensate and ρ ( λ ) we obtained aboveare under the lattice regularization and require renormal-ization. To renormalize the scalar quark bi-linear oper-ator, we use the regularization independent momentumsubtraction (RI/MOM) scheme [26, 27] under the Lan-dau gauge,¯ Z S ( Q ) = 12 Z RIq ( Q )Tr[Λ( p, I )] p = − Q = C m q + Z S ( Q ) + O ( m q ) ,Z RIq ( Q ) = Z A
48 Tr[ γ γ µ Λ( p, γ µ γ )] p = − Q ,Z A = 2 m q (cid:104) ¯ ψγ ψ | π (cid:105) m π (cid:104) ¯ ψγ γ ψ | π (cid:105) | m q → , (7)and further convert the result into MS 2GeV as (cid:104) ¯ ψψ (cid:105) r = Z S (cid:104) ¯ ψψ (cid:105) b and ρ ( λ ) r = Z m ρ ( λ ) b = Z − S ρ ( λ ) b . Follow-ing the definition of the RI/MOM scheme, we define thequark self energy through the axial vector normaliza-tion constant Z A whose uncertainty is generally less than0.1% at all the lattice spacings. The relation Z m Z S = 1is satisfied automatically and it avoids the systematic un-certainty from the additional chiral symmetry breakingin the Wilson-like actions. In term of Z − S , the residualRI/MOM scheme Q dependence is proportional to a and will vanish in the continuum, since the slopes of a Q shown in Fig. 1 are roughly the same at all the latticespacings we have. The statistical uncertainty in Fig. 1 isat 0.1% level due to the volume source propagator [28]with a given Q , but the total systematic uncertaintieswill be ∼ ∼
90% of this 1.5%), and also the value ofΛ
QCD , scale running, lattice spacing, and also fit range.Fortunately most of the systematic uncertainties are fullycorrelated at all the lattice spacings, and will not be en- larged in the continuum extrapolation. Our results ofthe scalar current renormalization constants are listed inTable I, with two uncertainties from the statistics andsystematics. l (GeV)( p r ( l )) (GeV) dl =0.5 MeV dl =1.0 MeV dl =2.0 MeV dl =4.0 MeV FIG. 2. The λ dependence of (cid:0) πρ ( λ ) (cid:1) / results at 0.084 fm,and different curves show the results with different bin sizes.The number of eigenvalues in the bin will be smaller than 25when δλ < λ ≥ Results:
In Fig. 2, we plot our (cid:0) πρ ( λ ) (cid:1) / result at0.084 fm with 0.5, 1, 2, 4 MeV bin size. The uncer-tainty and fluctuation of our results become much largerwith smaller bin size since the the number of eigenvaluesin each bin is fewer and violates the central limit the-orem requirement, but the uncertainty is just 1% withthe smallest bin size. With the 0.5 MeV bin size, onlythe first two data points drop significantly and it sug-gests that the finite volume effect is relatively small withthe 5.5 fm box at the physical pion mass. Based on thestandard statistical requirement to have at least 25 sam-ples in each bin ( ∼ ρ ( λ )calculations [5–7, 9, 10, 12]. It is understandable sincethe number of eigenvalues is almost proportional to thephysical volume, and we used the largest volume to datein this work.Next, we show the results at both 0.084 fm and 0.114fm in Fig. 3. If we ignore the difference between thephysical point (the u / d quark masses m l ∼ πρ ( λ ) is roughly the chiral condensate Σ atthe N f =2 chiral limit when λ (cid:28) m s . As shown in Fig. 3,our results of (cid:0) πρ ( λ ) (cid:1) / (red and purple bands) with λ <
50 MeV are consistent with the present lattice QCDaverages Σ / =272(5) MeV [29] (gray band). On theother hand, the results at two lattice spacings have verysignificant difference ( ∼
5% on ρ ( λ ) with λ ∼
30 MeV)given their 0.2% uncertainties. Such a difference can beaccounted for in the range of λ ∈ (10 , λ Σ F a term, resulting in the purple and red points l (GeV)( p r ( l )) (GeV)0.114 fm0.084 fm0.114 fm, a l subtracted0.084 fm, a l subtractedFLAG FIG. 3. The (cid:0) πρ ( λ ) (cid:1) / at 0.084 fm (red band) and 0.114fm (purple band) with physical light and strange quarkmasses. By subtracting a single λ Σ F a term, the dataat two lattice spacings (red and purple data points) canagree with each other within 1% difference. The grayband shows the present lattice QCD average of Σ / = (cid:0) π lim m l → lim λ → lim V →∞ ρ ( λ ) (cid:1) / [1, 29]. in Fig. 3. After the subtraction, we can still see visibledifferences at ρ ∼ χ PT, we always subtract the λ Σ F a term and drop thedata points with λ <
13 MeV.We start from the N f =2+1 PQ χ PT, which takes thestrange quark into account. With the light and strangequark masses obtained in Ref. [23] and the NLO formulaprovided in Ref. [2], the χ /d.o.f. is generally large re-gardless of the range of λ we use. Since the NLO cor-rection between the light quark mass and pion mass is ∼
10% which is much larger than the statistical uncer-tainty of the quark mass itself, we treat the quark mass m l , chiral condensate Σ and pion decay constant F inthe N f = 3 chiral limit as free parameters, to fit the datain the range λ ∈ (15 , O ( a ) correction. The results are listed in Ta-ble II. The strange quark mass is still treated as an inputusing the values determined in Ref. [23]. Our estimateof Σ / , F and m l are 232.6(0.9)(1.2)(0.7)(0.7) MeV,67.8(1.2)(0.1)(3.1)(0.8) MeV and 2.945(13)(44)(48)(8)MeV respectively, with four uncertainties coming fromthe statistical fluctuations, renormalization constants,continuum extrapolation (taking the difference betweenthose at the finest lattice spacing and continuum), andthe uncertainty of the lattice spacing determination.Note that the systematic uncertainty of F from therenormalization constant is small as it just depends onthe accurate Z A but not on the inaccurate Z S . Com-paring to the previous determinations of Σ / and F ,our determinations are much more precise. On the otherhand, the value we obtain for m l is 7 σ (or 12%) lowerthan the value 3.34(5) MeV in Ref. [23], which explains why our attempts to fit the data with that quark massresult in a large χ . Thus the NNLO effect in the ρ ( λ ) ex-pression could be essential to obtain the light quark massaccurately, and we should take this 12% as a systematicuncertainty in our determination of m l .An interesting alternative is to treat the strange quarkmass as a free parameter, in order to check what the datacan tell us in this case. For this choice, we obtain Σ / , F , m l and m s as 231(8)(1)(1)(1) MeV, 66(8)(0)(1)(1)MeV, 2.89(35)(5)(11)(1) MeV and 97(24)(1)(13)(0), re-spectively. The value of m s is in agreement with thelattice average 92.0(1.1) MeV, although it is not preciseafter the uncertainty is doubled during the continuum ex-trapolation, and the other quantities agree with the pre-vious determination with fixed m s , but the uncertaintiesare much larger.For comparison, we also consider the N f =2 PQ χ PTcase, which assumes the strange quark mass to be muchheavier than λ and is thus decoupled from the analysis.As shown in Eq. (5), ρ ( λ ) is independent of λ with N f =2in the leading order up to the finite volume correction,but a linear λ dependence appears in the N f = 3 case.The N f =2 form is consistent with what we see in Fig. 3and the N f =2 approximation only works when λ is muchsmaller than the strange quark mass. Since the available λ range is much smaller and we don’t know the exactrange we should use, we have to fix either the pion de-cay constant F in the N f = 2 chiral limit or the lightquark mass m l in the fitting. With the light quark mass3.34(5) MeV obtained in Ref. [23] and the NLO formulaprovided in Ref. [2], we get the N f =2 chiral limit Σ / ∼
265 MeV and F ∼
120 MeV with χ /d.o.f. ∼
10 or morein the range λ ∈ (15 ,
40) MeV, and F will be even largerif we reduce the maximum λ used in the fit. In contrast,if we use the lattice QCD average F = 86 . λ ∈ (12 ,
30) MeV with χ /d.o.f. around 1. (We reduce the bin size to 3 MeV inorder to have 6 data points per bin.) Then we fit the Σ / and m l at two lattice spacings separately and extrapo-late them to the continuum with linear O ( a ) correc-tion, as listed in Table II. Our estimates of Σ / and m l are 260.3(0.7)(1.3)(0.7)(0.8) MeV and 4.34(35)(6)(43)(1)MeV, respectively. It is understandable that the quarkmass we obtained can be 30(14)% larger than the latticeQCD averages 3.36(4) MeV due to the NLO effect, as wediscussed in the N f =2+1 case. Summary and outlook : Based on the precise calcula-tion of the spectral density ρ ( λ ) of D/ on two latticespacings with a statistical uncertainty at the 0.2% level,we determine the chiral condensate at MS 2 GeV to beΣ = (260 . . . . .
8) MeV) in the N f = 2 chi-ral limit and Σ = (232 . . . . .
7) MeV) inthe N f = 3 chiral limit; and the pion decay constant inthe N f = 3 chiral limit is F = 67 . . . . . TABLE II. Summary of the low energy constants and light quark masses obtained from the fit using the finite volume NLOPQ χ PT forms in Ref. [2]. N f =2 N f =2+1 N f =2+1, fit m s Lattice spacing Σ m l Σ F m l Σ F m l m s tics, renormalization constant, continuum extrapolation,and lattice spacing determination. The ratio Σ / Σ =1.40(2)(2) reflects the difference between the N f =2 chirallimit ( m s ∼
90 MeV) and the N f =3 chiral limit ( m s =0),as in the NLO expression of Z v .Thanks to the enhancement of the chiral log, the phys-ical light quark mass m l ∼ Z v ( λ, m l )differ from 1 (that in the chiral limit) by 10%. Thus wecan use ρ ( λ ) to determine the light quark mass precisely,up to the systematic uncertainty from NNLO corrections.On the other hand, we observe a clear λ dependence of ρ ( λ ) in the range of λ ∈ (20 , N f =2+1 NLO PQ χ PT form. It means that thestrange quark mass has a significant impact on ρ ( λ ).At the same time, determining the D/ spectrum and lowenergy constants are not the only application of solvingthe low lying eigenvectors of the overlap fermion. Thoseeigenvectors dominate in the quark propagator when thequark mass is as light as the physical one. Thus thoseeigenvectors we obtained also pave the way to calculatethe quark mass dependent pseudoscalar mass and decayconstant efficiently, which can be used to determine thelight and strange quark masses precisely and cross checkour determination of the low energy constants Σ, Σ , and F above. ACKNOWLEDGMENT
We thank the RBC and UKQCD collaborations forproviding us their DWF gauge configurations. The cal-culations were performed using the GWU code [30, 31]through HIP programming model [32]. This work is par-tially supported by Guangdong Major Project of Ba-sic and Applied Basic Research No. 2020B0301030008and Science and Technology Program of Guangzhou No.2019050001. A. A is supported in part by U.S. De-partment of Energy, Office of Science, Office of Nu-clear Physics under grant no DE-FG02-95ER40907. Y.B. is supported in part by the National Natural Sci-ence Foundation of China (NNSFC) under Grant No.12075253. T. D and K. L are supported by the U.S.DOE Grant No. DE-SC0013065 and DOE Grant No.DE-AC05-06OR23177 which is within the framework ofthe TMD Topical Collaboration. Y. Y is also supportedby the Strategic Priority Research Program of ChineseAcademy of Sciences, Grant No. XDB34030303. Thenumerical calculation is supported by the Strategic Pri-ority Research Program of Chinese Academy of Sciences,Grant No. XDC01040100, and also the supercomputingsystem in the Southern Nuclear Science Computing Cen-ter (SNSC). [1] T. Banks and A. Casher, Nucl. Phys. B , 103 (1980).[2] P. H. Damgaard and H. Fukaya, JHEP , 052 (2009),arXiv:0812.2797 [hep-lat].[3] T.-W. Chiu and S. V. Zenkin, Phys. Rev. D59 , 074501(1999), arXiv:hep-lat/9806019 [hep-lat].[4] R. C. Brower, H. Neff, and K. Orginos, Comput. Phys.Commun. , 1 (2017), arXiv:1206.5214 [hep-lat].[5] H. Fukaya, S. Aoki, T. Chiu, S. Hashimoto, T. Kaneko,J. Noaki, T. Onogi, and N. Yamada (JLQCD, TWQCD),Phys. Rev. D , 074501 (2011), arXiv:1012.4052 [hep-lat].[6] H. Fukaya, S. Aoki, S. Hashimoto, T. Kaneko, J. Noaki,T. Onogi, and N. Yamada (JLQCD), Phys. Rev. Lett. , 122002 (2010), [Erratum: Phys.Rev.Lett. 105,159901 (2010)], arXiv:0911.5555 [hep-lat].[7] G. Cossu, H. Fukaya, S. Hashimoto, T. Kaneko, and J.-I. Noaki, PTEP , 093B06 (2016), arXiv:1607.01099 [hep-lat].[8] L. Del Debbio, L. Giusti, M. Luscher, R. Petronzio, andN. Tantalo, JHEP , 011 (2006), arXiv:hep-lat/0512021[hep-lat].[9] L. Giusti and M. Luscher, JHEP , 013 (2009),arXiv:0812.3638 [hep-lat].[10] K. Cichy, E. Garcia-Ramos, and K. Jansen, JHEP ,175 (2013), arXiv:1303.1954 [hep-lat].[11] G. P. Engel, L. Giusti, S. Lottini, and R. Sommer, Phys.Rev. Lett. , 112001 (2015), arXiv:1406.4987 [hep-ph].[12] G. P. Engel, L. Giusti, S. Lottini, and R. Sommer, Phys.Rev. D , 054505 (2015), arXiv:1411.6386 [hep-lat].[13] A. Li et al. ( χ QCD), Phys. Rev.
D82 , 114501 (2010),arXiv:1005.5424 [hep-lat].[14] Y.-B. Yang, A. Alexandru, T. Draper, M. Gong, and K.-F. Liu, Phys. Rev.
D93 , 034503 (2016), arXiv:1509.04616[hep-lat]. [15] J. Liang, Y.-B. Yang, K.-F. Liu, A. Alexandru,T. Draper, and R. S. Sufian, Phys. Rev.
D96 , 034519(2017), arXiv:1612.04388 [hep-lat].[16] M. Gong et al. ( χ QCD), Phys. Rev.
D88 , 014503 (2013),arXiv:1304.1194 [hep-ph].[17] T.-W. Chiu, Phys. Rev.
D60 , 034503 (1999), arXiv:hep-lat/9810052 [hep-lat].[18] K.-F. Liu, Int. J. Mod. Phys.
A20 , 7241 (2005),arXiv:hep-lat/0206002 [hep-lat].[19] P. H. Ginsparg and K. G. Wilson, Phys. Rev.
D25 , 2649(1982).[20] L. Giusti, C. Hoelbling, M. Luscher, and H. Wittig,Comput. Phys. Commun. , 31 (2003), arXiv:hep-lat/0212012 [hep-lat].[21] W. E. Arnoldi, Quart. Appl. Math. (1951).[22] P. Hasenfratz and H. Leutwyler, Nucl. Phys. B343 , 241(1990).[23] T. Blum et al. (RBC, UKQCD), Phys. Rev.
D93 , 074505(2016), arXiv:1411.7017 [hep-lat].[24] R. D. Mawhinney (RBC, UKQCD), (2019),arXiv:1912.13150 [hep-lat].[25] A. Hasenfratz, R. Hoffmann, and F. Knechtli,
Contentsof lattice 2001 proceedings , Nucl. Phys. Proc. Suppl. ,418 (2002), [,418(2001)], arXiv:hep-lat/0110168 [hep-lat]. [26] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa, andA. Vladikas, Nucl. Phys.
B445 , 81 (1995), arXiv:hep-lat/9411010 [hep-lat].[27] Y. Bi, H. Cai, Y. Chen, M. Gong, K.-F. Liu, Z. Liu,and Y.-B. Yang, Phys. Rev.
D97 , 094501 (2018),arXiv:1710.08678 [hep-lat].[28] J.-W. Chen, T. Ishikawa, L. Jin, H.-W. Lin, Y.-B. Yang,J.-H. Zhang, and Y. Zhao, Phys. Rev.
D97 , 014505(2018), arXiv:1706.01295 [hep-lat].[29] S. Aoki et al. (Flavour Lattice Averaging Group), Eur.Phys. J.
C80 , 113 (2020), arXiv:1902.08191 [hep-lat].[30] A. Alexandru, C. Pelissier, B. Gamari, and F. Lee, J.Comput. Phys. , 1866 (2012), arXiv:1103.5103 [hep-lat].[31] A. Alexandru, M. Lujan, C. Pelissier, B. Gamari, andF. X. Lee, in
Proceedings, 2011 Symposium on Ap-plication Accelerators in High-Performance Computing(SAAHPC’11): Knoxville, Tennessee, July 19-20, 2011 (2011) pp. 123–130, arXiv:1106.4964 [hep-lat].[32] Y.-J. Bi, Y. Xiao, W.-Y. Guo, M. Gong, P. Sun, S. Xu,and Y.-B. Yang,
Proceedings, 37th International Sympo-sium on Lattice Field Theory (Lattice 2019): Wuhan,China, June 16-22 2019 , PoS