Topological susceptibility in 2+1-flavor QCD with chiral fermions

We compute the topological susceptibility $\chi_t$ of 2+1-flavor lattice QCD with dynamical M\"obius domain-wall 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"sub-volume 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 $M_\pi$ and decay constant $F_\pi$, we obtain $\chi_t$ = 0.227(02)(11)$M_\pi^2 F_\pi^2$ at the physical point, where the first error is statistical and the second is systematic.


Introduction
Computing the topological susceptibility, χ t , 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 χ t comes entirely from sea quarks, or a small quantum effect suppressed by O( ), 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 χ t and its comparison with the ChPT formula [2][3][4], 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, l = l r 3 − l r 7 + h r 1 − h r 3 is a combination of the low-energy constants (renormalized at M phys ) at next-to-leading order (NLO) [8], and M phys and F phys are the physical values of the pion mass and decay constant, respectively. In this work we improve the computation of χ t in three ways. The first is to employ the Möbius domain-wall fermion [9] for the dynamical quarks. It realizes good chiral symmetry, i.e. the residual mass is kept at the order of 1 MeV or lower. As will be shown below, our results have only a mild Speaker, e-mail: hfukaya[at]het.phys.sci.osaka-u.ac.jp dependence on the lattice spacing, up to a ∼ 0.08 fm. Unlike the overlap fermion that we employed in the previous studies [5], the use of the domain-wall fermion allows us to sample configurations in different topological sectors.
The second improvement comes from the use of sub-volumes of the simulated lattices. Since the correlation length of QCD is limited, at most by 1/M π , there is no fundamental reason to use the global topological charge to compute χ t . The use of sub-volumes was tested in our previous simulations with overlap quarks [5,7], where the signal of χ t 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] (similar methods were proposed in [14] and [15]). 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% sub-volumes of the whole lattice, whose size is ∼ 2 fm, are sufficient to extract χ t . 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 M 2 π F 2 π , a product of the pion mass M π and decay constant F π squared calculated on the same ensemble. This ratio, where l = −l r 4 − l r 7 + h r 1 − h r 3 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 χ t , M π , and F π at each simulation point.
We also employ the Yang-Mills (YM) gradient flow [16][17][18] 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 χ t 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 l , taking the chiral and continuum limits from formulas (1) and (2). The same set of data was also used to calculate the η meson mass [19], 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 [20].

Numerical set-up
In the numerical simulation of QCD, we use the Symanzik gauge action and the Möbius domain-wall fermion action for gauge ensemble generations [21][22][23]. We apply three steps of stout smearing of the gauge links for the computation of the Dirac operator. Our main runs of 2 + 1-flavor lattice QCD simulations are performed with two different lattice sizes L 3 ×T = 32 3 ×64 and 48 3 ×96, for which we set β = 4.17 and 4.35, respectively. The inverse lattice spacing 1/a is estimated to be 2.453(4) GeV (for β = 4.17) and 3.610(9) GeV (for β = 4.35), using the input √ t 0 = 0.1465 fm [24]. The two lattices share a similar physical size L ∼ 2.6 fm. For the quark masses, we choose two values of the strange quark mass m s around its physical point, and 3-4 values of the up and down quark masses m ud for each m s . For our lightest pion mass around 230 MeV (am ud = 0.0035 at β = 4.17) we simulate a larger lattice 48 3 × 96 with the same set of the parameters to control the finite volume effects. We also perform a simulation on a finer lattice 64 3 × 128 (at β = 4.47 [1/a ∼ 4.5 GeV] and M π ∼ 285 MeV). For each ensemble, 5000 molecular dynamics (MD) time is simulated. Numerical works are done with the QCD software package IroIro++ [23].
In this setup, we confirm that the violation of the chiral symmetry in the Möbius domain-wall fermion formalism is well under control. The residual mass is ∼ 1 MeV [25] by choosing the lattice size in the fifth direction L 5 = 12 at β = 4.17 and less than 0.2 MeV with L 5 = 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 step-size a 2 ∆t =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 auto-correlation 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 auto-correlation 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 [26].

Topology fluctuation in a "slab"
We use the conventional gluonic definition of the topological charge density q lat (x), the so-called clover construction [27]. Since the YM gradient flow smooths the gauge field in the range of √ 8t ∼ 0.5 fm of the lattice, a simple summation Q lat = x q lat (x) over all sites gives values close to integers.
As is well known, the global topological charge Q lat suffers from long auto-correlation time in lattice simulations, especially when the lattice spacing is small. This is true also in our simulations. At the highest β = 4.47, Q lat drifts very slowly with auto-correlation time of O(1000).
Instead of using the global topological charge Q lat , we extract the topological susceptibility from a sub-volume V sub of the whole lattice V. Since the correlation length of QCD is limited by at most 1/M π , the subvolume V sub should contain sufficient information to extract χ t , provided that its size is larger than 1/M π . One can then effectively increase the statistics by V/V sub , since each piece of V/V sub sub-lattices may be considered as an uncorrelated sample. Moreover, there is no potential barrier among topological sectors: the instantons and anti-instantons freely come in and go out of the sub-volume, which should make the auto-correlation time of the observable shorter than that of the global topological charge.
There are various ways of cutting the whole lattice into sub-volumes and computing the correlation functions in them. After some trial and error, we find that the so-called "slab" method, proposed by Bietenholz et al. [12] is efficient for the purpose of computing χ t . The idea is to sum up the two-point correlators of the topological charge density, over x and y in the same sub-volume: Here the integration over x and y in the spatial directions runs in the whole spatial volume while the temporal sum is restricted to the region [T ref , T cut + T ref ], which is called a "slab". Since the YM gradient flow is already performed, there is no divergence from the points of x = y. T ref denotes an arbitrary reference time. Due to the translational invariance, the slabs of the same thickness T cut are equivalent, and one can average over T ref . This method is statistically more stable than the other sub-volume method we applied in [5,7] since Q 2 slab (T cut ) is guaranteed to be always positive. If we sample large statistics on a large enough lattice volume, Q 2 slab (T cut ) should be a fraction T cut /T of χ t V. Namely, Q 2 slab (T cut ) should be a linear function in T cut . Its leading finite volume correction can be estimated using the formula in [11]: where C is an unknown constant, and m 0 is the mass of the first excited state, the η meson. Note that for 1/m 0 T cut T − 1/m 0 , the formula gives a simple linear function in T cut plus a constant. Also, note that in the limit of T cut = T , Q 2 slab (T cut = T ) reduces to the global topology Q 2 = χ t V. Assuming the linearity in T cut , one can extract the topological susceptibility through with two reference thicknesses t 1 and t 2 . In our numerical analysis, T ref is averaged over the temporal direction. Since the data at t i and T − t i are not independent, we choose t 1 and t 2 in a range 1.6 fm < t 1 , t 2 < T/2. In the numerical analysis, we replace q lat (x) by q lat (x) − Q lat /V to cancel a possible bias due to the long auto-correlation of the global topology. We find that the signal using this slab method is less noisy than the previous attempts in [5,7]. Moreover, the new method shows more frequent fluctuation than that of the global topological charge on our finest lattice.

Results
At β = 4.17, which corresponds to the lattice spacing a ∼ 0.08 fm, both Q lat and Q 2 slab (T cut ) fluctuate well. The data on this lattice, therefore, provide a good testing ground to examine the validity of the slab sub-volume method, comparing with the naive definition of the topological susceptibility, Q 2 lat /V. In Fig. 1, Q 2 slab (t cut ) observed at the lightest sea quark mass m ud = 0.0035, β = 4.17 on two different volumes L = 32 and L = 48 are plotted as a function of T cut /T . The data converge to the linear plus constant function (4) at T cut = 20, which corresponds to ∼ 1.6 fm. The slope, or χ slab t , is consistent with that from global topology shown by solid and dotted lines for L = 32 and L = 48 lattices, respectively. We also observe the consistency between the L = 32 and L = 48 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 m s , 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 auto-correlation of the topological charge, and the larger deviation of Q lat from zero. We determine the shorter reference t 1 ∼ 1.6 fm using the data at the lightest quark mass and always choose t 2 = T/2 ∼ 2.6 fm. In order to estimate the systematic errors due to non-linear behavior, we compare the results with 1) those obtained from different reference times (t 1 , t 2 ) = (t 1 , t 1 +t 2 2 ), and ( t 1 +t 2 2 , t 2 ), and 2) those obtained without the subtraction of Q /V 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 χ slab t from all ensembles plotted in physical units. The error-bar in the plot represents the statistical and systematic errors due to the choice of the references t 1 and t 2 , added in quadrature. On the horizontal axis, the quark mass defined in the MS scheme at 2 GeV is   where the renormalization factor Z S is nonperturbatively computed in [29]. We find a linear suppression near the chiral limit. We also find that there is no clear dependence on β and m s . First, we compare our results directly to the ChPT formula (1). We perform a two-parameter (Σ and l) fit to the data at β = 4.17 (solid curve in Fig. 3) and β ≥ 4.35 (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 l but Σ is stable against the change in the fit-range. 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 M π and decay constant F π together with χ slab t , obtained from each ensemble, we take the ratio (2). By a linear one-parameter fit as shown in the bottom panel of Fig. 3, we determine l and the ratio χ slab t /(M π F π ) 2 at the physical point. In the same way as the determination of Σ and l, we take the chiral and continuum limits of both quantities. Note that the fixed chiral limit at 1/4 of the ratio helps us to determine these quantities. Finally let us discuss other possible systematic effects. In our analysis, the ensembles satisfying M π L > 3.9 are used and we do not expect any sizable finite volume effects. In particular, our lightest mass point has M π L = 4.4. We have used configurations at the YM gradient flow-time around √ 8t ∼ 0.5 fm. We confirm that the flow-time dependence is negligible in the range 0.25 fm < √ 8t < 0.5 fm.