Topological susceptibility in 2+1flavor QCD with chiral fermions
Abstract
We compute the topological susceptibility of 2+1flavor lattice QCD with dynamical Möbius domainwall fermions, whose residual mass is kept at 1 MeV or smaller. In our analysis, we focus on the fluctuation of the topological charge density in a “slab” subvolume of the simulated lattice, as proposed by Bietenholz et al. The quark mass dependence of our results agrees well with the prediction of the chiral perturbation theory, from which the chiral condensate is extracted. Combining the results for the pion mass and decay constant , we obtain = 0.227(02)(11) at the physical point, where the first error is statistical and the second is systematic.
EPJ Web of Conferences \woctitleLattice2017 \selectlanguageenglish
1 Introduction
Computing the topological susceptibility, , has been a challenging task for lattice QCD. One of the problems is that the topological charge is sensitive to the violation of chiral symmetry. The quark mass dependence of comes entirely from sea quarks, or a small quantum effect suppressed by , to which contribution from the discretization systematics is relatively large. Simulating QCD on a sufficiently fine lattice is another challenge, as the global topological charge tends to be frozen along the Monte Carlo history [1]. Due to these difficulties, the study of the quark mass dependence of and its comparison with the ChPT formula [2, 3, 4],
(1) 
has been very limited, and only some pilot works with dynamical chiral fermions on rather small or coarse lattices have been performed [5, 6, 7]. Here denotes the chiral condensate, is a combination of the lowenergy constants (renormalized at ) at nexttoleading order (NLO) [8], and and are the physical values of the pion mass and decay constant, respectively.
In this work we improve the computation of in three ways. The first is to employ the Möbius domainwall fermion [9] for the dynamical quarks. It realizes good chiral symmetry, i.e. the residual mass is kept at the order of MeV or lower. As will be shown below, our results have only a mild dependence on the lattice spacing, up to fm. Unlike the overlap fermion that we employed in the previous studies [5], the use of the domainwall fermion allows us to sample configurations in different topological sectors.
The second improvement comes from the use of subvolumes of the simulated lattices. Since the correlation length of QCD is limited, at most by , there is no fundamental reason to use the global topological charge to compute . The use of subvolumes was tested in our previous simulations with overlap quarks [5, 7], where the signal of was extracted from finite volume effects on the meson correlator [10, 11]. In this work, we utilize a different method, which was originally proposed by Bietenholz et al. [12, 13, 14] (similar methods were proposed in [15] and [16]). The method is based on a correlator of local topological charge density, which is designed to be always non negative, and thus to reduce statisical noise. We confirm that 30%–50% subvolumes of the whole lattice, whose size is fm, are sufficient to extract . Moreover, the new definition shows more frequent fluctuation than that of the global topological charge on our finest lattice.
The third improvement is to take the ratio of the topological susceptibility to , a product of the pion mass and decay constant squared calculated on the same ensemble. This ratio,
(2) 
where is a combination of the NLO low energy constants, is free from the chiral logs, as well as possible finite volume effects at NLO. Moreover, the chiral limit of the ratio, 1/4, is protected from the strange sea quark effects. We can, therefore, precisely estimate the topological susceptibility at the physical point by measuring , , and at each simulation point.
We also employ the YangMills (YM) gradient flow [17, 18, 19] in order to make the global topological charge close to integers, to remove the UV divergences, and to reduce the statistical noise. With these improvements, we achieve good enough statistical precision to investigate the dependence of on the sea quark mass. In fact, the topological susceptibility shows a good agreement with the ChPT prediction (1). We determine the value of chiral condensate, as well as , taking the chiral and continuum limits from formulas (1) and (2). The same set of data was also used to calculate the meson mass [20], which was extracted from the shorter distance region of the correlator of the topological charge density. Further details of this work may be found in [21].
2 Numerical setup
In the numerical simulation of QCD, we use the Symanzik gauge action and the Möbius domainwall fermion action for gauge ensemble generations [22, 23, 24]. We apply three steps of stout smearing of the gauge links for the computation of the Dirac operator. Our main runs of flavor lattice QCD simulations are performed with two different lattice sizes and , for which we set = 4.17 and 4.35, respectively. The inverse lattice spacing is estimated to be 2.453(4) GeV (for ) and 3.610(9) GeV (for ), using the input fm [25]. The two lattices share a similar physical size fm. For the quark masses, we choose two values of the strange quark mass around its physical point, and 3–4 values of the up and down quark masses for each . For our lightest pion mass around 230 MeV ( = 0.0035 at = 4.17) we simulate a larger lattice with the same set of the parameters to control the finite volume effects. We also perform a simulation on a finer lattice (at [ GeV] and 285 MeV). For each ensemble, 5000 molecular dynamics (MD) time is simulated. Numerical works are done with the QCD software package IroIro++ [24].
In this setup, we confirm that the violation of the chiral symmetry in the Möbius domainwall fermion formalism is well under control. The residual mass is MeV [26] by choosing the lattice size in the fifth direction = 12 at = 4.17 and less than 0.2 MeV with = 8 at = 4.35 (and 4.47).
On generated configurations, we perform 500–1640 steps of the YM gradient flow (using the conventional Wilson gauge action) with a stepsize 0.01. At every 200–400 steps (depending on the parameters) we store the configuration of the topological charge density and compute its twopoint correlator.
In the following analysis, we measure the integrated autocorrelation time of every quantity, following the method proposed in [1]. The statistical error is estimated by the jackknife method (without binning) multiplied by the square root of autocorrelation time.
The pion mass and decay constant are computed combining the pseudoscalar correlators with local and smeared source operators. Details of the computation are presented in a separate article [27].
3 Topology fluctuation in a “slab”
We use the conventional gluonic definition of the topological charge density , the socalled clover construction [28]. Since the YM gradient flow smooths the gauge field in the range of fm of the lattice, a simple summation over all sites gives values close to integers.
As is well known, the global topological charge suffers from long autocorrelation time in lattice simulations, especially when the lattice spacing is small. This is true also in our simulations. At the highest , drifts very slowly with autocorrelation time of .
Instead of using the global topological charge , we extract the topological susceptibility from a subvolume of the whole lattice . Since the correlation length of QCD is limited by at most , the subvolume should contain sufficient information to extract , provided that its size is larger than . One can then effectively increase the statistics by , since each piece of sublattices may be considered as an uncorrelated sample. Moreover, there is no potential barrier among topological sectors: the instantons and antiinstantons freely come in and go out of the subvolume, which should make the autocorrelation time of the observable shorter than that of the global topological charge.
There are various ways of cutting the whole lattice into subvolumes and computing the correlation functions in them. After some trial and error, we find that the socalled “slab” method, proposed by Bietenholz et al. [12] is efficient for the purpose of computing . The idea is to sum up the twopoint correlators of the topological charge density, over and in the same subvolume:
(3) 
Here the integration over and in the spatial directions runs in the whole spatial volume while the temporal sum is restricted to the region , which is called a “slab”. Since the YM gradient flow is already performed, there is no divergence from the points of . denotes an arbitrary reference time. Due to the translational invariance, the slabs of the same thickness are equivalent, and one can average over . This method is statistically more stable than the other subvolume method we applied in [5, 7] since is guaranteed to be always positive.
If we sample large statistics on a large enough lattice volume, should be a fraction of . Namely, should be a linear function in . Its leading finite volume correction can be estimated using the formula in [11]:
(4) 
where is an unknown constant, and is the mass of the first excited state, the meson. Note that for , the formula gives a simple linear function in plus a constant. Also, note that in the limit of , reduces to the global topology .
Assuming the linearity in , one can extract the topological susceptibility through
(5) 
with two reference thicknesses and . In our numerical analysis, is averaged over the temporal direction. Since the data at and are not independent, we choose and in a range 1.6 fm . In the numerical analysis, we replace by to cancel a possible bias due to the long autocorrelation of the global topology.
4 Results
At , which corresponds to the lattice spacing fm, both and fluctuate well. The data on this lattice, therefore, provide a good testing ground to examine the validity of the slab subvolume method, comparing with the naive definition of the topological susceptibility, .
In Fig. 1, observed at the lightest sea quark mass , on two different volumes and are plotted as a function of . The data converge to the linear plus constant function (4) at , which corresponds to fm. The slope, or , is consistent with that from global topology shown by solid and dotted lines for and lattices, respectively. We also observe the consistency between the and data, which suggests that the systematics due to the finite volume is well under control. The “linear constant” behavior is also seen in ensembles with heavier quark masses. Moreover, the extracted values of the topological susceptibility from the slope show a good agreement with the ChPT prediction, which confirms the validity of the slab method.
At higher values, we also find a reasonable slope at the lightest quark mass for each and , as shown in Fig. 2. For heavier masses, however, some curvature is seen. We consider this curvature is an effect from the bias of the global topological charge. This observation is consistent with previous works (see [1] for example), which reported that the heavier pion mass ensembles show the longer autocorrelation of the topological charge, and the larger deviation of from zero. We determine the shorter reference fm using the data at the lightest quark mass and always choose fm. In order to estimate the systematic errors due to nonlinear behavior, we compare the results with 1) those obtained from different reference times , and , and 2) those obtained without the subtraction of in the definition of the topological charge density. The larger deviation is treated as a systematic error. For the statistical error estimates we follow the method proposed by the ALPHA collaboration [1] assuming a double exponential structure of the autocorrelation function.
The top panel of Fig. 3 presents our data for from all ensembles plotted in physical units. The errorbar in the plot represents the statistical and systematic errors due to the choice of the references and , added in quadrature. On the horizontal axis, the quark mass defined in the scheme at 2 GeV is
(6) 
where the renormalization factor is nonperturbatively computed in [30]. We find a linear suppression near the chiral limit. We also find that there is no clear dependence on and .
First, we compare our results directly to the ChPT formula (1). We perform a twoparameter ( and ) fit to the data at (solid curve in Fig. 3) and (dashed curve) separately. Here we also perform the same fit but omitting the heaviest two points, and take the difference as an estimate for the systematic error in the chiral extrapolation. Since the heaviest points have the following problems, 1) a strong bias is seen in the global topology, 2) ChPT is less reliable, and 3) there is a mismatch between different , we take the result without them as our central values. Note, however, this inclusion/elimination affects but is stable against the change in the fitrange. Namely, the chiral condensate is effectively determined by the lower quark mass data. We then estimate the continuum limit by a constant fit. Comparing our result from the constant fit with a linear extrapolation of the central values, we take the difference as an estimate for the systematic error in the continuum extrapolation.
Next, using our data for the pion mass and decay constant together with , obtained from each ensemble, we take the ratio (2). By a linear oneparameter fit as shown in the bottom panel of Fig. 3, we determine and the ratio at the physical point. In the same way as the determination of and , we take the chiral and continuum limits of both quantities. Note that the fixed chiral limit at of the ratio helps us to determine these quantities.
Finally let us discuss other possible systematic effects. In our analysis, the ensembles satisfying are used and we do not expect any sizable finite volume effects. In particular, our lightest mass point has . We have used configurations at the YM gradient flowtime around fm. We confirm that the flow–time dependence is negligible in the range fm fm. We conclude that all these systematic effects are negligibly small compared to the statistical and systematic errors given above.
5 Summary and discussion
With dynamical Möbius domainwall fermions and the new method using slab subvolumes of the simulated lattice, we have computed the topological susceptibility of QCD. Its quark mass dependence is consistent with the ChPT prediction, from which we have obtained
(7)  
(8) 
where the first error comes from the statistical uncertainty at each simulation point, including the effect of freezing topology. The second and third represent the systematics in the chiral and continuum limits, respectively. The value of is consistent with our recent determination through the Dirac spectrum [31]. We have also estimated the NLO coefficient
(9)  
(10) 
where is renormalized at the physical pion mass, while is renormalization invariant. It is interesting to note that and include a combination of the coefficients , which is supposed to be unphysical in ChPT unless dependence is considered. These are important for possible couplings of QCD to axions [32].
We thank T. Izubuchi, and other members of JLQCD collaboration for useful discussions. Numerical simulations are performed on IBM System Blue Gene Solution at KEK under a support of its Large Scale Simulation Program (No. 16/1714). This work is supported in part by the Japanese GrantinAid for Scientific Research (Nos. JP25800147, JP26247043, JP26400259, JP16H03978), and by MEXT as “Priority Issue on PostK computer” (Elucidation of the Fundamental Laws and Evolution of the Universe) and by Joint Institute for Computational Fundamental Science (JICFuS). The work of GC is supported by STFC, grant ST/L000458/1.
References
 S. Schaefer et al. [ALPHA Collaboration], Nucl. Phys. B 845, 93 (2011) doi:10.1016/j.nuclphysb.2010.11.020 [arXiv:1009.5228 [heplat]].
 Y. Y. Mao et al. [TWQCD Collaboration], Phys. Rev. D 80, 034502 (2009) doi:10.1103/PhysRevD.80.034502 [arXiv:0903.2146 [heplat]].
 S. Aoki and H. Fukaya, Phys. Rev. D 81, 034022 (2010) doi:10.1103/PhysRevD.81.034022 [arXiv:0906.4852 [heplat]].
 F. K. Guo and U. G. Meißner, Phys. Lett. B 749, 278 (2015) doi:10.1016/j.physletb.2015.07.076 [arXiv:1506.05487 [hepph]].
 S. Aoki et al. [JLQCD and TWQCD Collaborations], Phys. Lett. B 665, 294 (2008) [arXiv:0710.1130 [heplat]].
 T. W. Chiu et al. [JLQCD and TWQCD Collaborations], PoS LATTICE 2008, 072 (2008) [arXiv:0810.0085 [heplat]].
 T. H. Hsieh et al. [JLQCD and TWQCD Collaborations], PoS LAT 2009, 085 (2009).
 J. Gasser and H. Leutwyler, Annals Phys. 158, 142 (1984). doi:10.1016/00034916(84)902422
 R. C. Brower, H. Neff and K. Orginos, Nucl. Phys. Proc. Suppl. 140, 686 (2005) doi:10.1016/j.nuclphysbps.2004.11.180 [heplat/0409118]; arXiv:1206.5214 [heplat].
 R. Brower, S. Chandrasekharan, J. W. Negele and U. J. Wiese, Phys. Lett. B 560, 64 (2003) doi:10.1016/S03702693(03)003691 [heplat/0302005].
 S. Aoki, H. Fukaya, S. Hashimoto and T. Onogi, Phys. Rev. D 76, 054508 (2007) [arXiv:0707.0396 [heplat]].
 W. Bietenholz, P. de Forcrand and U. Gerber, JHEP 1512, 070 (2015) doi:10.1007/JHEP12(2015)070 [arXiv:1509.06433 [heplat]].
 W. Bietenholz, K. Cichy, P. de Forcrand, A. Dromard and U. Gerber, PoS LATTICE 2016, 321 (2016) [arXiv:1610.00685 [heplat]].
 H. MejíaDíaz, W. Bietenholz, P. de Forcrand, U. Gerber and I. O. Sandoval, arXiv:1712.01395 [heplat].
 E. V. Shuryak and J. J. M. Verbaarschot, Phys. Rev. D 52, 295 (1995) doi:10.1103/PhysRevD.52.295 [heplat/9409020].
 P. de Forcrand, M. Garcia Perez, J. E. Hetrick, E. Laermann, J. F. Lagae and I. O. Stamatescu, Nucl. Phys. Proc. Suppl. 73, 578 (1999) doi:10.1016/S09205632(99)851433 [heplat/9810033].
 M. Lüscher, JHEP 1008, 071 (2010) [Erratumibid. 1403, 092 (2014)] [arXiv:1006.4518 [heplat]].
 M. Lüscher and P. Weisz, JHEP 1102, 051 (2011) [arXiv:1101.0963 [hepth]].
 C. Bonati and M. D’Elia, Phys. Rev. D 89, no. 10, 105005 (2014) [arXiv:1401.2441 [heplat]].
 H. Fukaya et al. [JLQCD Collaboration], arXiv:1509.00944 [heplat].
 S. Aoki et al. [JLQCD Collaboration], arXiv:1705.10906 [heplat].
 T. Kaneko et al. [JLQCD Collaboration], PoS LATTICE 2013, 125 (2014) [arXiv:1311.6941 [heplat]].
 J. Noaki et al. [JLQCD Collaboration], PoS LATTICE 2013, 263 (2014).

G. Cossu, J. Noaki, S. Hashimoto, T. Kaneko, H. Fukaya, P. A. Boyle and J. Doi,
arXiv:1311.0084 [heplat];
URL :http://suchix.kek.jp/guido_cossu/documents/DoxyGen/html/index.html  S. Borsanyi et al., JHEP 1209, 010 (2012) [arXiv:1203.4469 [heplat]].
 S. Hashimoto, S. Aoki, G. Cossu, H. Fukaya, T. Kaneko, J. Noaki and P. A. Boyle, PoS LATTICE 2013, 431 (2014).
 B. Fahy et al. [JLQCD Collaboration], PoS LATTICE 2015, 074 (2016) [arXiv:1512.08599 [heplat]].
 M. Bruno et al. [ALPHA Collaboration], JHEP 1408, 150 (2014) doi:10.1007/JHEP08(2014)150 [arXiv:1406.5363 [heplat]].
 J. Gasser and H. Leutwyler, Nucl. Phys. B 250, 465 (1985). doi:10.1016/05503213(85)904924
 M. Tomii et al. [JLQCD Collaboration], Phys. Rev. D 94, no. 5, 054504 (2016) doi:10.1103/PhysRevD.94.054504 [arXiv:1604.08702 [heplat]].
 G. Cossu, H. Fukaya, S. Hashimoto, T. Kaneko and J. I. Noaki, PTEP 2016, no. 9, 093B06 (2016) doi:10.1093/ptep/ptw129 [arXiv:1607.01099 [heplat]].
 G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, JHEP 1601, 034 (2016) doi:10.1007/JHEP01(2016)034 [arXiv:1511.02867 [hepph]].