Topological susceptibility of QCD with dynamical M\"obius domain wall fermions

TL;DR

This paper calculates the topological susceptibility in lattice QCD with dynamical M"obius domain wall fermions, confirming chiral perturbation theory predictions and extracting the chiral condensate with high precision.

Contribution

It presents a novel lattice QCD calculation of topological susceptibility using M"obius domain wall fermions with controlled chiral symmetry violation.

Findings

Chiral condensate determined as (274(13)(29) MeV)^3

Topological susceptibility at physical point matches theoretical predictions

Residual chiral symmetry breaking kept below 1 MeV

Abstract

We compute the topological susceptibility $\chi_t$ of lattice QCD with $2+1$ dynamical quark flavors described by the M\"obius domain wall fermion. Violation of chiral symmetry as measured by the residual mass is kept at $\sim$1 MeV or smaller. We measure the fluctuation of the topological charge density in a `slab' sub-volume of the simulated lattice using the method proposed by Bietenholz {\it et al.} The quark mass dependence of $\chi_t$ is consistent with the prediction of chiral perturbation theory, from which the chiral condensate is extracted as $\Sigma^{\overline{\rm MS}} (\mbox{2GeV}) = [274(13)(29)\mbox{MeV}]^3$, where the first error is statistical and the second one is systematic. Combining the results for the pion mass $M_\pi$ and decay constant $F_\pi$, we obtain $\chi_t = 0.229(03)(13)M_\pi^2F_\pi^2$ at the physical point.

Tables

Table 1. Table 1: Parameters of the JLQCD gauge ensembles used in this work. Pion masses are rounded to the nearest 10 10 10 MeV.

Lattice Spacing	$L^3\times T$	$L_5$	$a{m}_{ud}$	$a{m}_s$	$m_{\pi }\text{ [MeV]}$	$m_{\pi }L$
$\beta =4.17,$	${32}^3\times 64$ $\left(L=2.6\text{ fm}\right)$	12	0.0035	0.040	230	3.0
$a^{-1}=2.453\left(4\right)\text{ GeV}$			0.0070	0.030	310	4.0
			0.0070	0.040	310	4.0
			0.0120	0.030	400	5.2
			0.0120	0.040	400	5.2
			0.0190	0.030	500	6.5
			0.0190	0.040	500	6.5
	${48}^3\times 96$ $\left(L=3.9\text{ fm}\right)$	12	0.0035	0.040	230	4.4
$\beta =4.35,$	${48}^3\times 96$ $\left(L=2.6\text{ fm}\right)$	8	0.0042	0.018	300	3.9
$a^{-1}=3.610\left(9\right)\text{ GeV}$			0.0042	0.025	300	3.9
			0.0080	0.018	410	5.4
			0.0080	0.025	410	5.4
			0.0120	0.018	500	6.6
			0.0120	0.025	500	6.6
$\beta =4.47,$	${64}^3\times 128$ $\left(L=2.7\text{ fm}\right)$	8	0.0030	0.015	280	4.0
$a^{-1}=4.496\left(9\right)\text{ GeV}$

Table 2. Table 2: Results for the pion mass M π subscript 𝑀 𝜋 M_{\pi} , decay constant 2 ​ F π 2 subscript 𝐹 𝜋 \sqrt{2}F_{\pi} , τ exp subscript 𝜏 exp \tau_{\rm exp} , τ imp subscript 𝜏 imp \tau_{\rm imp} , and χ t slab superscript subscript 𝜒 𝑡 slab \chi_{t}^{\rm slab} . τ exp subscript 𝜏 exp \tau_{\rm exp} for β = 4.47 𝛽 4.47 \beta=4.47 is estimated from the first zero-crossing point of Q lat subscript 𝑄 lat Q_{\rm lat} . All the data are shown in the lattice units. For χ t slab superscript subscript 𝜒 𝑡 slab \chi_{t}^{\rm slab} , the first error denotes the statistical error, while the second shows the systematic error due to the effect of freezing global topological charge.

$\beta$	$L$	$m_{ud}$	$m_s$	$M_{\pi }$	$\sqrt{2}{F}_{\pi }$	${\tau }_{\exp }$	${\tau }_{\mathrm{imp}}$	${\chi }_t^{\mathrm{slab}}$
4.17	32	0.0035	0.04	0.09369(32)	0.05320(19)	17(04)	25(9)	0.217(64)(14)$\times {10}^{-5}$
		0.007	0.04	0.12604(26)	0.05774(15)	14(03)	30(9)	0.400(78)(21)$\times {10}^{-5}$
		0.012	0.04	0.16267(22)	0.06254(14)	65(40)	62(24)	1.01(32)(46)$\times {10}^{-5}$
		0.019	0.04	0.20329(19)	0.06788(14)	65(40)	56(22)	1.59(38)(09)$\times {10}^{-5}$
		0.007	0.03	0.12629(26)	0.05761(15)	29(07)	54(22)	0.56(15)(53)$\times {10}^{-5}$
		0.012	0.03	0.16179(21)	0.06190(14)	74(50)	71(32)	0.56(19)(22)$\times {10}^{-5}$
		0.019	0.03	0.20302(20)	0.06730(13)	56(35)	42(16)	1.20(24)(23)$\times {10}^{-5}$
	48	0.0035	0.04	0.09203(09)	0.05440(09)	38(30)	21(06)	0.282(34)(42)$\times {10}^{-5}$
4.35	48	0.0042	0.025	0.08299(18)	0.03926(11)	243(153)	208(114)	0.91(40)(12)$\times {10}^{-6}$
		0.0080	0.025	0.11312(14)	0.04291(09)	318(200)	362(234)	2.18(98)(48)$\times {10}^{-6}$
		0.0120	0.025	0.13875(14)	0.04630(08)	173(142)	105(52)	1.21(32)(07)$\times {10}^{-6}$
		0.0042	0.018	0.08219(19)	0.03901(11)	111(49)	158(72)	0.59(37)(12)$\times {10}^{-6}$
		0.0080	0.018	0.11284(15)	0.04275(08)	236(148)	220(126)	0.55(53)(19)$\times {10}^{-6}$
		0.0120	0.018	0.13799(13)	0.04603(09)	97(43)	170(82)	1.70(53)(21)$\times {10}^{-6}$
4.47	64	0.0030	0.015	0.06316(15)	0.03141(09)	[1700]	492(836)	0.20(27)(09)$\times {10}^{-6}$

Table 3. Table 3: Our results for Σ Σ \Sigma , l 𝑙 l , l ′ superscript 𝑙 ′ l^{\prime} , and χ t slab / ( M π ​ F π ) 2 superscript subscript 𝜒 𝑡 slab superscript subscript 𝑀 𝜋 subscript 𝐹 𝜋 2 \chi_{t}^{\rm slab}/(M_{\pi}F_{\pi})^{2} at the physical point. The first error denotes the the statistic fluctuation at each simulation point, including the effect of long autocorrelation of global topology. The second is the systematic error in chiral extrapolation, and the third error denotes that in the continuum limit estimates. See the main text for the details.

$\beta$	${\Sigma }^{1/3}$(MeV)	$l$	$l^{\prime }$	$\frac{{\chi }_t^{\mathrm{slab}}}{M_{\pi }^2{F}_{\pi }^2}$ at physical point
4.17	275(13)(13)	0.003(06)(10)	$-$0.018(03)(03)	0.232(04)(03)
$\ge$ 4.35	261(50)(19)	-0.005(09)(06)	$-$0.025(05)(04)	0.223(05)(04)
continuum limit	274(13)(25)(15)	-0.001(05)(06)(19)	$-$0.019(03)(01)(13)	0.229(03)(01)(13)

Equations

$${\chi }_t=\frac{m_{ud}\Sigma }{2}\left\{1-\frac{3m_{ud}\Sigma }{16{\pi }^2{F}_{\mathrm{phys}}^4}\ln \left(\frac{2m_{ud}\Sigma }{F_{\mathrm{phys}}^2{M}_{\mathrm{phys}}^2}\right)+\frac{4m_{ud}\Sigma }{F_{\mathrm{phys}}^4}l\right\},$$

$$\frac{{\chi }_t}{M_{\pi }^2{F}_{\pi }^2}=\frac{1}{4}\left[1+\frac{2M_{\pi }^2{l}^{\prime }}{F_{\pi }^2}\right],$$

$$⟨Q_{\mathrm{slab}}^2(T_{\mathrm{cut}})⟩$$

$$⟨Q_{\mathrm{slab}}^2(T_{\mathrm{cut}})⟩=({\chi }_{t}V)\times \frac{T_{\mathrm{cut}}}{T}+C(1-e^{-m_0{T}_{\mathrm{cut}}})(1-e^{-m_0(T-T_{\mathrm{cut}})}),$$

$${\chi }_t^{\mathrm{slab}}=\frac{T}{V}\left[\frac{⟨Q_{\mathrm{slab}}^2(t_1)⟩-⟨Q_{\mathrm{slab}}^2(t_2)⟩}{t_1-t_2}\right],$$

$${\Gamma }_O(\Delta \tau )=⟨O(\tau )O(\tau +\Delta \tau ){⟩}_{\tau },$$

$${\tau }_{\mathrm{int}}=\frac{1}{2}+\sum _{\Delta \tau =0}^W\rho (\Delta \tau ),\rho (\Delta \tau )=\frac{{\Gamma }_O(\Delta \tau )}{{\Gamma }_O(0)},$$

$${\tau }_{\mathrm{imp}}={\tau }_{\mathrm{int}}^{\prime }+{\tau }_{\exp }\rho (W^{\prime }),$$

$$m_{ud}^{\overline{\text{mbox}MS}}=(m_{ud}+m_{res})/Z_S,$$

$${\chi }_t$$

$${\Sigma }^{\overline{\mathrm{MS}}}(\text{mbox}2GeV)$$

$$l$$

$$l^{\prime }$$

$${\chi }_t$$

$$\Delta (M^2)=\frac{M^2}{16{\pi }^2}\ln \frac{M^2}{{\mu }_{sub}^2}+g_1(M^2),$$

$$M_{\pi }^2=M^2\left[1-\frac{1}{F_{\pi }^2}\left\{-\frac{1}{2}\Delta (M_{\pi }^2)+\frac{1}{6}\Delta (M_{\eta }^2)+8(L_4^r-2L_6^r)(2M_{\pi }^2+M_{ss}^2)+8(L_5^r-2L_8^r)M_{\pi }^2\right\}\right],$$

$$F_{\pi }^2=F^2\left[1-\frac{1}{F_{\pi }^2}\left\{2\Delta (M_{\pi }^2)+\Delta (M_K^2)-8L_4^r(2M_{\pi }^2+M_{ss}^2)-8L_5^r{M}_{\pi }^2\right\}\right],$$

$$\bar{m}\sim \frac{m}{2}\left(1-\frac{m_{ud}}{2m_s}\right)\sim \frac{m}{2}\left(1-\frac{M_{\pi }^2}{2M_{ss}^2}\right),$$

$$\frac{{\chi }_t}{M_{\pi }^2{F}_{\pi }^2}$$

$$l_{(\mathrm{eff})}^{\prime }=-\frac{1}{4M_{ss}^2}\left(F_{\pi }^2+\Delta (M_K^2)+\frac{1}{2}\Delta (M_{\eta }^2)\right)+36L_7+4L_8^r.$$

$$⟨Q_{\mathrm{slab}}^2(T_{\mathrm{cut}}){⟩}_Q=({\chi }_{t}V)\times \frac{T_{\mathrm{cut}}}{T}+\frac{T_{\mathrm{cut}}^2}{T^2}\left(Q^2-{\chi }_{t}V\right)+C,$$

$$⟨Q_{\mathrm{slab}}^2(T_{\mathrm{cut}}){⟩}_{biased}$$

Full text

1 Introduction

The topological susceptibility $\chi_{t}$ is an interesting quantity which characterizes how many topological excitations are created in the QCD vacuum. Witten Witten:1979vv and Veneziano Veneziano:1980xs estimated $\chi_{t}$ in the large-$N_{c}$ (number of colors) limit and showed that it is proportional to the square of the $\eta^{\prime}$ meson mass. In real QCD with $N_{c}=3$ and light dynamical quarks, however, the argument of Witten and Veneziano is no longer valid. It is not the $\eta^{\prime}$ meson but the (zero-momentum mode of) pion that controls the topological susceptibility.

According to the prediction of $SU(2)$ chiral perturbation theory (ChPT) at leading order (LO), $\chi_{t}$ is expected to be proportional to the quark mass $m_{ud}$, when the up and down quark masses are degenerate. At one-loop, the quark mass dependence is predicted as Mao:2009sy ; Aoki:2009mx ; Guo:2015oxa

[TABLE]

where $\Sigma$ denotes the chiral condensate, $l=l_{3}^{r}-l_{7}^{r}+h_{1}^{r}-h_{3}^{r}$ is a combination of the low-energy constants at next-to-leading order (NLO) Gasser:1983yg and $M_{\rm phys}$(=135 MeV) and $F_{\rm phys}$(=92 MeV) are the physical values of the pion mass and decay constant, respectively. Here $l_{i}^{r}$ are renormalized at $M_{\rm phys}$. In the formula, we have assumed that the strange quark is decoupled from the theory and $SU(2)$ chiral perturbation theory works. In other words, the strange quark mass dependence is assumed to be absorbed in the low energy constants.

By taking a ratio with the ChPT predictions for the pion mass and decay constant (let us denote them by $M_{\pi}$ and $F_{\pi}$), one can eliminate the chiral logarithm in  (1):

[TABLE]

where $l^{\prime}=-l_{4}^{r}-l_{7}^{r}+h_{1}^{r}-h_{3}^{r}$ is again a combination of the NLO low energy constants, which is independent of the renormalization scheme and scale at this order. This ratio also cancels possible finite volume effects at NLO. Moreover, the chiral limit of the ratio, 1/4, is protected from the strange sea quark effects (see Appendix A for the details), as they always enter as a function of the ratio $m_{ud}/m_{s}$, which can be absorbed into the (finite) renormalization of $l^{\prime}$. We can therefore precisely estimate the topological susceptibility at the physical point by measuring $\chi_{t}$, $M_{\pi}$, and $F_{\pi}$ at each simulation point.

It has been a challenging task for lattice QCD to compute $\chi_{t}$, since it is sensitive to the discretization effects and the violation of chiral symmetry Bazavov:2010xr ; Cichy:2013rra ; Bruno:2014ova in particular. This is partly because the quark mass dependence of $\chi_{t}$ is due to sea quarks, or a small quantum mechanical effect suppressed by $O(\hbar)$, to which the discretization error is relatively large. Even if we could simulate QCD on a sufficiently fine lattice, the global topological charge would become frozen along the Monte Carlo history Schaefer:2010hu . Due to these difficulties, the study of the quark mass dependence and its comparison with the ChPT formula of $\chi_{t}$ has been very limited, and only some pilot works with dynamical chiral fermions on rather small or coarse lattices Egri:2005cx ; Aoki:2007pw ; DeGrand:2007tx ; Chiu:2008kt ; Chiu:2008jq ; Hsieh:2009zz ; Aoki:2010dy ; Chiu:2011dz ; Chiu:2011za are available.

In this work, we improve the computation of $\chi_{t}$ in two ways. One is to employ the domain wall fermion Kaplan:1992bt ; Shamir:1993zy with an improvement by Brower:2004xi , or known as the Möbius domain wall fermion, for the dynamical quarks, which enables us to precisely preserve chiral symmetry. Even on our coarsest lattice, the residual mass, related to the violation of the chiral symmetry, is kept at the order of $1$ MeV. As will be shown below, our results show only a mild dependence of $\chi_{t}$ on the lattice spacing, up to $a\sim 0.08$ fm. The use of the domain-wall fermion allows us to sample configurations in different topological sectors, which is also an advantage over the simulation with the overlap fermion where we fixed the global topological charge in Aoki:2007pw .

Another 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_{\pi}$, there is essentially no need to use the global topological charge to compute $\chi_{t}$. The use of sub-volume was tested in our previous simulations with overlap quarks Aoki:2007pw ; Hsieh:2009zz (see also Bautista:2015yza ; Bietenholz:2016ymo ), where the signal was extracted from finite volume effects, which have some sensitivity to $\chi_{t}$ Brower:2003yx ; Aoki:2007ka . In this work, we utilize a different method, which was originally proposed by Bietenholz et al. Bietenholz:2015rsa ; Bietenholz:2016szu (similar methods were proposed in Shuryak:1994rr and deForcrand:1998ng ). The method is based on a correlator, which gives a positive finite value even in the thermodynamical limit, and thus is less noisy than our recent attempt in Fukaya:2014zda 111 See also Ref. Brower:2014bqa where a similar method to ours was attempted. . We confirm that 30%–50% sub-volumes of the whole lattice, whose size is $\sim 2$ fm, are sufficient to extract $\chi_{t}$. Moreover, the new definition shows more frequent fluctuation than that of the global topological charge on our finest lattice.

We also employ a modern technique, the Yang-Mills (YM) gradient flow Luscher:2010iy ; Luscher:2011bx ; Bonati:2014tqa , 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 $\chi_{t}$ on the sea quark mass. In fact, our data of the topological susceptibility is consistent with the ChPT prediction (1), from which the values of chiral condensate and $l^{\prime}$ are extracted.

The same set of data was also used to calculate the $\eta^{\prime}$ meson mass Fukaya:2015ara , which was extracted from the shorter distance region of the correlator of the topological charge density. These two results show a nontrivial double-scale structure of topological fluctuation of the gauge field: it creates the $\eta^{\prime}$ meson at short distances, while describing the vacuum mode of the pion (or the lowest mode, which is constant over space-time) at long distances.

The rest of this paper is organized as follows. First, we describe our lattice set-up in Sec. 2. We then explain the method to extract the topological susceptibility from the slab sub-volume in Sec. 3. Our results at lower $\beta$ are presented in Sec. 4. Comparing the data with those obtained from global topology, we examine the validity of our sub-volume method. The results at higher $\beta$ are shown in Sec. 5, and how we estimate the statistical errors is explained in Sec. 6. Finally, we estimate the chiral and continuum limits in Sec. 7 and give a summary in Sec. 8.

2 Lattice set-up

In the numerical simulation222 Numerical works are done with the QCD software package IroIro++ Cossu:2013ola . of QCD, we use the Symanzik gauge action and the Möbius domain wall fermion action for gauge ensemble generations Kaneko:2013jla ; Noaki:2014ura ; Cossu:2013ola . We apply three steps of stout smearing of the gauge links Morningstar:2003gk 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}\times T=32^{3}\times 64$ and $48^{3}\times 96$, for which we set $\beta$ = 4.17 and 4.35, respectively. The inverse lattice spacing $1/a$ is estimated to be 2.453(4) GeV (for $\beta=4.17$) and 3.610(9) GeV (for $\beta=4.35$), using the input $\sqrt{t_{0}}=0.1465$ fm Borsanyi:2012zs where we use the reference YM gradient flow–time $t_{0}$, defined by $t^{2}\langle E\rangle|_{t=t_{0}}=0.3$ Luscher:2010iy with the energy density $E$ of the gluon field. The two lattices share a similar physical size $L\sim 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}$. Since our data at the lightest pion mass around 230 MeV ($am_{ud}$ = 0.0035 at $\beta$ = 4.17) may contain significant finite size effects, we simulate a larger lattice $48^{3}\times 96$ with the same set of the parameters to check if the finite volume systematics is small enough. We also perform a simulation on a finer lattice $64^{3}\times 128$ (at $\beta=4.47$ [$1/a\sim 4.5$ GeV] and $M_{\pi}\sim$ 285 MeV). For each ensemble, 500–1000 gauge configurations are sampled from 5000 molecular dynamics (MD) time. The ensembles used in this work are listed in Table. 1.

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 $\sim 1$ MeV Hashimoto:2014gta by choosing the lattice size in the fifth direction $L_{5}$ = 12 at $\beta$ = 4.17 and less than 0.2 MeV with $L_{5}$ = 8 at $\beta$ = 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 $\Delta t/a^{2}=$0.01. At every 200–400 steps (depending on the parameters) we store the configuration of the topological charge density. The two-point correlator is measured using the Fast Fourier Transform (FFT) technique to average source and sink points over whole lattice sites.

In the following analysis, we measure the integrated auto-correlation time of every quantity, following the method proposed in Luscher:2005rx ; Schaefer:2010hu . The statistical error is estimated by the jackknife method (without binning) multiplied by the square root of auto-correlation time normalized by the MD time interval of the configuration samples. We will discuss more details about the auto-correlation time of topological fluctuations in Sec. 6.

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 Fahy:2015xka .

3 Topological susceptibility in a “slab”

We use the conventional gluonic definition of the topological charge density $q^{\rm lat}(x)$, the so-called clover construction Bruno:2014ova . Since the YM gradient flow smooths the gauge field in the range of $\sqrt{8t}\sim 0.5$ fm of the lattice, a simple summation $Q_{\rm lat}=\sum_{x}q^{\rm lat}(x)$ over the whole sites gives values close to integers, as shown in Fig. 1.

As is well known, the global topological charge $Q_{\rm lat}$ suffers from long auto-correlation time in lattice simulations, especially when the lattice spacing is small. This is true also in our simulations, as shown in Fig. 2. At the highest $\beta=4.47$, $Q_{\rm lat}$ drifts very slowly with auto-correlation time of possibly $O(1000)$. It is, therefore, not feasible to estimate $\chi_{t}$ without performing much longer runs. The details of the auto-correlation time of the topological charge and its density operator will be discussed in Section 6.

Instead of using the global topological charge $Q_{\rm lat}$, we attempt to extract the topological susceptibility from a sub-volume $V_{\rm sub}$ of the whole lattice $V$. Since the correlation length of QCD is limited by at most $1/M_{\pi}$, the subvolume $V_{\rm sub}$ should contain sufficient information to extract $\chi_{t}$, provided that its size is larger than $1/M_{\pi}$. One can then effectively increase the statistics by $V/V_{\rm sub}$, since each piece of $V/V_{\rm 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 compute the correlation functions in them. After some trials and errors, we find that the so-called “slab” method, proposed by Bietenholz et al. Bietenholz:2015rsa is efficient for the purpose of computing $\chi_{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:

[TABLE]

Here the integration over $x$ and $y$ in the spatial directions runs in the whole spatial volume (since the YM gradient flow is already performed, there is no divergence from the points of $x=y$), while the temporal sum is restricted to the region $[T_{\rm ref},T_{\rm cut}+T_{\rm ref}]$, which is called a “slab”. Here $T_{\rm ref}$ is an arbitrary reference time. Due to the translational invariance, the slabs sharing the same thickness $T_{\rm cut}$ are physically equivalent, and one can average over $T_{\rm ref}$. This method is statistically more stable than the other sub-volume method we applied in Aoki:2007pw ; Hsieh:2009zz because $\langle Q_{\rm slab}^{2}(T_{\rm cut})\rangle$ is guaranteed to be always positive.

If we sample large statistics on a large enough lattice volume, $\langle Q_{\rm slab}^{2}(T_{\rm cut})\rangle$ should be just $T_{\rm cut}/T$ of $\chi_{t}V$. Namely $\langle Q_{\rm slab}^{2}(T_{\rm cut})\rangle$ should be a linear function in $T_{\rm cut}$. Its leading finite volume correction can be estimated using the formula in Aoki:2007ka :

[TABLE]

where $C$ is an unknown constant, and $m_{0}$ is the mass of the first excited state, the $\eta^{\prime}$ meson333 The finite volume effects are due to propagation of the mesons in the flavor singlet channel. As the ground state or the $\eta^{\prime}$ meson is heavy, we neglect the higher order effects. Even if we include them, the structure of linear+constant in Eq. 4 is unchanged since their effect is just an additional constant. . Note that for $1/m_{0}\ll T_{\rm cut}\ll T-1/m_{0}$, the formula gives a simple linear function in $T_{\rm cut}$ plus a constant. Also, note that in the limit of $T_{\rm cut}=T$, $\langle Q_{\rm slab}^{2}(T_{\rm cut}=T)\rangle=\langle Q^{2}\rangle=\chi_{t}V$.

Assuming the linearity in $T_{\rm cut}$, one can extract the topological susceptibility through

[TABLE]

with two reference thicknesses $t_{1}$ and $t_{2}$. In our numerical analysis, $T_{\rm 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^{\rm lat}(x)$ by $q^{\rm lat}(x)-\langle Q_{\rm lat}/V\rangle$ to cancel a possible bias due to the long auto-correlation of the global topology.

The original proposal in Bietenholz:2015rsa mainly used the correlator in a fixed topological sector. The formula corresponding to (4) then contains a subtraction of the contribution from the global topology. We find that the statistical noise is larger with this choice while the results from different topological sectors are consistent. In the following analysis, we use (4) after summing over the topological sectors.

We find that the signal using this slab method is less noisy than the previous attempts in Aoki:2007pw ; Hsieh:2009zz . Moreover, as shown in Fig. 2 and will be discussed in details later, the new definition shows more frequent fluctuation than that of the global topological charge on our finest lattice.

4 Results at low $\beta$

At $\beta=4.17$, which corresponds to the lattice spacing $a\sim 0.08$ fm, both the global topological charge $Q_{\rm lat}$ and $Q_{\rm slab}^{2}(T_{\rm cut})$ fluctuate well, as shown in the top panel of Fig. 2. 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 with $\langle Q_{\rm lat}^{2}\rangle/V$.

In Fig. 3, $\langle Q_{\rm slab}^{2}(t_{\rm cut})\rangle$ observed at the lightest sea quark mass $m_{ud}=0.0035$, $\beta=4.17$ on two different volumes $L=32$ and $L=48$ are plotted as a function of $T_{\rm cut}/T$. The data converge to a linear plus constant function given in  (4) at $T_{\rm cut}=20$, which corresponds to $\sim 1.6$ fm. The slope, or $\chi_{t}^{\rm slab}$, 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, as presented in Fig. 4.

The extracted values of the topological susceptibility from the slope, show a good agreement with the ChPT prediction, as shown in Fig. 5 by open and filled squares. The leading-order ChPT formula, $\chi_{t}=m_{ud}\Sigma/2$, with $\Sigma=[270\mbox{MeV}]^{3}$ (solid line) is drawn as an eye guide. In the same plot, we also plotted the estimate for $\chi_{t}$ obtained from the global topological charge by circles, which again agree with the results, validating the slab method. The values of $\chi_{t}^{\rm slab}$ are listed in Table 2. How we estimate their errorbars is explained in the following two sections.

5 Results at high beta

At higher $\beta$ values, we still find a reasonable slope at the lightest quark mass for each $\beta$ and $m_{s}$, as shown in Fig. 6. 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 Schaefer:2010hu for example), which reported that the heavier pion mass ensembles show the longer auto-correlation of the topological charge, and the larger deviation of $\langle Q_{\rm lat}\rangle$ from zero. We determine the reference $t_{1}\sim 1.6$ fm using data at the lightest quark mass and always choose $t_{2}=T/2\sim 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}^{\prime},t_{2}^{\prime})=(t_{1},\frac{t_{1}+t_{2}}{2})$, and $(\frac{t_{1}+t_{2}}{2},t_{2})$, and 2) those obtained without the subtraction of $\langle Q\rangle/V$ in the definition of the topological charge density. The larger deviation is treated as a systematic error. More details are presented in Appendix B.

Our results are summarized in Fig. 7 (see also Fig. 8 for a comparison with Ref. Cichy:2013rra and Ref. Bruno:2014ova ). Although the data at higher $\beta$ are rather scattered compared to those at $\beta=4.17$, they can be used to estimate the chiral condensate $\Sigma$, assuming the linear suppression around the chiral limit. Before going to the details, we discuss the auto-correlation of $\chi_{t}^{\rm slab}$ and show how we estimate the statistical errors in the next section.

6 Auto-correlation and error estimates

Gauge configurations generated by a Markov chain are generally not independent but have auto-correlations. How much they are correlated depends on observables. We therefore need to carefully measure the auto-correlation of the target observable $O$:

[TABLE]

where $\tau$ denotes the Monte Carlo time, and the average $\langle\cdots\rangle_{\tau}$ is taken over $\tau$.

When the Monte Carlo trajectory is long enough, compared to the auto-correlation time of any observables, one can estimate the so-called integrated auto-correlation time by

[TABLE]

where the upper end of the summation window $W$ is chosen to where $\rho(\Delta\tau)$ becomes consistent with zero within the error. The above formula assumes that $\Gamma_{O}(\Delta\tau)$ converges to a single exponential function well below $W$.

If the observables suffer from long auto-correlation, and the Monte Carlo trajectory is not long enough, on the other hand, the above procedure may underestimate the auto-correlation time, since some very slow decay modes can be hidden in the error of $\rho(\Delta\tau)$. This problem is similar to that of hadron spectroscopy with a short temporal extension, where one does not have long enough fitting range to disentangle the ground state from excited states, which leads to over-estimation of the mass.

The ALPHA collaboration Schaefer:2010hu carefully studied the effect of slow modes, and proposed an improved estimate of the auto-correlation time,

[TABLE]

where $\tau_{\rm int}^{\prime}$ is the same summation as (7) but with a smaller upper bound $W^{\prime}$ where $\rho(W^{\prime})$ becomes lower than 3/2 standard deviations. $\tau_{\rm exp}$ is the auto-correlation of the slowest mode. The proposal is equivalent to considering a continuation of $\Gamma_{O}(\Delta\tau)$ at $\Delta\tau=W^{\prime}$ to the slowest possible exponential function $\Gamma_{O}(W^{\prime})\exp(-(\Delta\tau-W^{\prime})/\tau_{\rm exp})$.

In lattice QCD simulations, it is natural to assume that $\tau_{\exp}$ is equal to the auto-correlation of the global topological charge. In our simulations, $\tau_{\exp}$ is estimated by $\tau_{\rm int}(W)$ of $Q^{2}_{\rm lat}$, except for $\beta=4.47$ where we choose $\tau_{\exp}=1700$ MD time by hand (and assuming 100% error for it), which is a rough order estimate from the first zero-crossing point of $Q_{\rm lat}$. Then we compute the auto-correlation of our target observable $\chi_{t}^{\rm slab}$ by $\tau_{\rm imp}$ to estimate the error. The results for $\tau_{\rm imp}$, $\tau_{\exp}$ and $\chi_{t}^{\rm slab}$ are summarized in Table 2 and the auto-correlation function $\rho(\Delta\tau)$ at three different $\beta$ with a similar pion mass $M_{\pi}\sim 300$ MeV is shown in Fig. 9. At the highest $\beta=4.47$, it is clear that $\chi_{t}^{\rm slab}$ has shorter auto-correlation time than that of the global topological charge.

With the measured improved auto-correlation time $\tau_{\rm imp}$, we estimate the statistical errors of $\chi_{t}^{\rm slab}$ by multiplying $\sqrt{2(\tau_{\rm imp}+\Delta\tau_{\rm imp})/\tau_{interval}}$ to the naive error estimates, where $\Delta\tau_{\rm imp}$ is the standard deviation of $\tau_{\rm imp}$ and $\tau_{interval}$ denotes the interval trajectory between samples. The results, as well as the systematic error from the choice of reference points are listed in the last column of Table 2.

7 Chiral and continuum limit

Figure 7 presents our data for $\chi_{t}^{\rm slab}$ from all ensembles plotted in physical units. The horizontal axis, the quark mass defined in the $\overline{\mbox{MS}}$ scheme at 2 GeV is

[TABLE]

where the renormalization factor $Z_{S}$ is nonperturbatively computed in Tomii:2016xiv : $Z_{S}=1.037,\;0.934$, and $0.893$ for $\beta=4.17,\;4.35$, and $4.47$, respectively. In contrast to the results by other groups with non-chiral fermions, we find no strong dependence on $\beta$.

First, we compare our results directly to the ChPT formula (1). We perform a two-parameter ($\Sigma$ and $l$) fit to the data at $\beta=4.17$ (solid curve in Fig. 7) and $\beta\geq 4.35$ (dashed curve) separately444 Since $\beta=4.47$ is simulated at only one choice of the quark masses, we simply add the data as one of the $\beta=4.35$ ensembles. In fact, $\chi_{t}$ values at $\beta=4.35$ and 4.47 at the pion mass $\sim 300$ MeV are consistent with each other. . The results for $\Sigma$ and $l$ are listed in Table 3. 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 several problems, 1) a strong bias is seen in the global topology, 2) ChPT is less reliable, and 3) mismatch between different $\beta$, we take the result without them as our central values. Note, however, this inclusion/elimination affects $l$ but $\Sigma$ is stable against the change in the fit-range. Namely, the chiral condensate $\Sigma$ is determined by the low quark mass data. We then estimate the continuum limit by a constant fit, as shown in the top two panels in Fig. 10. Comparing our result from the constant fit with the linear extrapolation of the central values, we take the difference as an estimate for the systematic error in the continuum limit. In the plots in Fig. 10, all these errors are added in quadrature.

Next, using our data for the pion mass $M_{\pi}$ and decay constant $F_{\pi}$ together with $\chi_{t}^{\rm slab}$, obtained from each ensemble, we take the ratio given in  (2). By a linear one-parameter fit, we determine $l^{\prime}$ and the ratio $\chi_{t}^{\rm slab}/(M_{\pi}F_{\pi})^{2}$ at the physical point. In the same way as the determination of $\Sigma$ 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_{\pi}L>3.9$ are used and we do not expect any sizable finite volume effects. In particular, our lightest mass point has $M_{\pi}L=4.4$. We have used configurations at the YM gradient flow–time around $\sqrt{8t}\sim 0.5$ fm. We confirm that the flow–time dependence is negligible in the range $0.25$ fm $<\sqrt{8t}<0.5$ fm. cWe conclude that all these systematic effects are negligibly small compared to the statistical and systematic errors given above.

8 Summary

With dynamical Möbius domain wall fermions and the new method using sub-volume 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

[TABLE]

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 $\Sigma$ is consistent with our recent determination through Dirac spectrum Cossu:2016eqs . We have also estimated the NLO coefficient

[TABLE]

where $l$ is renormalized at the physical pion mass, while $l^{\prime}$ is renormalization invariant. It is interesting to note that $l$ and $l^{\prime}$ include a combination of the coefficients $h_{1}^{r}-h_{3}^{r}$, which is supposed to be unphysical in ChPT unless $\theta$ dependence is considered. These are important for possible couplings of QCD to axions diCortona:2015ldu .

We thank T. Izubuchi, and other members of JLQCD collaboration for fruitful discussions. We also thank the Yukawa Institute for Theoretical Physics, Kyoto University. Discussions during the YITP workshop YITP-T-14-03 on “Hadrons and Hadron Interactions in QCD” were useful in completing this work. Numerical simulations are performed on IBM System Blue Gene Solution at KEK under a support of its Large Scale Simulation Program (No. 16/17-14). This work is supported in part by the Japanese Grant-in-Aid for Scientific Research (Nos. JP25800147, JP26247043, JP26400259, JP16H03978), and by MEXT as “Priority Issue on Post-K 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.

Appendix A Effect of strange sea quark

In this work, we have assumed that effect of the strange quark is negligible and used $SU(2)$ ChPT in our main analysis to obtain the chiral extrapolation of the topological susceptibility. In this appendix, we consider $SU(3)$ ChPT and compute possible correction from the strange quark loop. We will show that the chiral limit of the ratio (2) is unchanged even in $SU(3)$ ChPT, which is also protected from finite volume corrections.

The one-loop computation of the topological susceptibility in general $N_{f}$-flavor ChPT was given in Mao:2009sy ; Aoki:2009mx and the formula for $N_{f}=3$ is

[TABLE]

where $\bar{m}=m_{ud}m_{s}/(2m_{s}+m_{ud})$, $M_{\pi}$,$M_{K}$, and $M_{\eta}$ are the (simulated) pion, kaon and $\eta$ meson masses, respectively. We have also used notations for $M^{2}_{ss}=2m_{s}\Sigma/F_{\pi}^{2}$, and $\bar{M}^{2}=2\bar{m}\Sigma/F_{\pi}^{2}$. The chiral logarithm is expressed by

[TABLE]

where $\mu_{sub}$ denotes the renormalization scale, and $g_{1}$ is finite volume correction (see Aoki:2009mx for the details). In the above formula, we can see three NLO low-energy constants Gasser:1984gg : $L_{6}^{r}$ and $L_{8}^{r}$ are those renormalized at $\mu_{sub}$, while $L_{7}$ is a renormalization scheme independent constant.

One-loop corrections to the pion mass and decay constant were computed in Gasser:1984gg :

[TABLE]

and

[TABLE]

where $M$ and $F$ are the tree-level mass and decay constant, respectively.

Now let us take the ratio of $\chi_{t}$ and $M_{\pi}^{2}F_{\pi}^{2}$. Noting

[TABLE]

we obtain

[TABLE]

where both of strange quark effect, as well as finite volume effects from one-loop diagrams are absorbed in the (re)definition of

[TABLE]

We, therefore, conclude that the one-loop formula (2) is valid even when strange quark gives nontrivial effect, and is also stable against possible finite volume corrections. This observation helps us in determining $\chi_{t}$ at the physical point.

Appendix B Bias from global topology

In this appendix, we discuss systematics due to freezing of the global topological charge. Combining the formulas in Bietenholz:2015rsa and Aoki:2007ka , the slab topological charge squared at fixed topology of $Q$ becomes

[TABLE]

for $0\ll T_{\rm cut}\ll T$. Therefore, if the global topological charge $Q$ were badly sampled and its average of $Q^{2}$ in the ensemble deviated from $\chi_{t}V$, we should have a quadratic term in $T_{\rm cut}$ as

[TABLE]

where $\langle\cdots\rangle_{biased}$ denotes the estimate obtained from a biased sampling of configurations. Here we have included the term $\langle Q\rangle^{2}_{biased}$, which comes from the use of the subtracted operator $q^{\rm lat}-\langle Q/V\rangle_{biased}$ in our numerical analysis.

If the correction $\frac{T_{\rm cut}}{(\chi_{t}V)T}\left(\langle Q^{2}\rangle_{biased}-\langle Q\rangle^{2}_{biased}-\chi_{t}V\right)$ is small, our original linear $+$ constant formula is still valid. As the correction is proportional to $\frac{T_{\rm cut}}{T}$, if we have a window $T_{\rm cut}\ll T$, or the freezing $\langle Q^{2}\rangle_{biased}-\langle Q\rangle^{2}_{biased}$ happens to be near the true value of $\chi_{t}V$, we can still extract $\chi_{t}$ from the linear slope (this seems to happen on the data at $\beta=4.47$).

In order to estimate the systematics due to the correction term, we compare the results with 1) those obtained from different reference times $(t_{1}^{\prime},t_{2}^{\prime})=(t_{1},\frac{t_{1}+t_{2}}{2})$, and $(\frac{t_{1}+t_{2}}{2},t_{2})$555 This also tests if the effect from the excited state $m_{0}$ is small or not. , and 2) those obtained without the subtraction of $\langle Q\rangle/V$ in the definition of the topological charge density. Then we take the larger deviation as the systematic error. Since a part of $\langle Q\rangle^{2}_{biased}$ is expected to be canceled by $\langle Q^{2}\rangle_{biased}$, this analysis is rather conservative. As presented in Tab. 2, the deviations are comparable to the statistical errors.

Let us look into our “worst” case, the data at $\beta=4.35$ and $(m_{ud},m_{s})=(0.012,0.018)$ in our ensembles, which shows the strongest curvature. As is expected, the global topological charge sampling is biased: the estimate for $\langle Q^{2}\rangle=12(4)$ in the former half (0-2500 MD time) of the simulation time, is quite different from $\langle Q^{2}\rangle=40(17)$ in the latter half (2500-5000 MD time). But the obtained values of $\chi_{t}^{\rm slab}$ show a milder deviation, $1.30(53)\times 10^{-6}$ for the former half and $1.89(64)\times 10^{-6}$ for the latter, which are consistent within errors. This analysis666 We thank W. Bietenholz and P. de Forcrand for suggesting this analysis of freezing topology effects. shows that the systematics due to freezing topology is under control, at least, at the level of the statistical errors. Our ChPT fit with reasonable $\chi^{2}/d.o.f.\sim 1.4$ also supports our conclusion.