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

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.

We compute the topological susceptibility χ t of lattice QCD with 2 + 1 dynamical quark flavors described by the Möbius domain wall fermion. Violation of chiral symmetry as measured by the residual mass is kept at ∼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 et al. The quark mass dependence of χ t is consistent with the prediction of chiral perturbation theory, from which the chiral condensate is extracted as Σ MS (2GeV) = [274 (13) (29)MeV] 3 , where the first error is statistical and the second one is systematic. Combining the results for the pion mass M π and decay constant F π , we obtain χ t = 0.229(03)(13)M 2 π F 2 π at the physical point.

Introduction
The topological susceptibility χ t is an interesting quantity which characterizes how many topological excitations are created in the QCD vacuum. Witten [1] and Veneziano [2] estimated χ t in the large-N c (number of colors) limit and showed that it is proportional to the square of the η ′ 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 η ′ 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), χ 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 [3][4][5] where Σ denotes the chiral condensate, l = l r 3 − l r 7 + h r 1 − h r 3 is a combination of the lowenergy constants at next-to-leading order (NLO) [6] and M phys (=135 MeV) and F phys (=92 MeV) are the physical values of the pion mass and decay constant, respectively. Here l r i are renormalized at M 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 π and F π ), one can eliminate the chiral logarithm in (1): where l ′ = −l r 4 − l r 7 + h r 1 − h r 3 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 ′ . We can therefore precisely estimate the topological susceptibility at the physical point by measuring χ t , M π , and F π at each simulation point.
It has been a challenging task for lattice QCD to compute χ t , since it is sensitive to the discretization effects and the violation of chiral symmetry [7][8][9] in particular. This is partly because the quark mass dependence of χ t is due to sea quarks, or a small quantum mechanical effect suppressed by O( ), 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 [10]. Due to these difficulties, the study of the quark mass dependence and its comparison with the ChPT formula of χ t has been very limited, and only some pilot works with dynamical chiral fermions on rather small or coarse lattices [11][12][13][14][15][16][17][18][19] are available.
In this work, we improve the computation of χ t in two ways. One is to employ the domain wall fermion [20,21] with an improvement by [22], 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 χ t on the lattice spacing, up to a ∼ 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 [12].
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 π , there is essentially no need to use the global topological charge to compute χ t . The use of sub-volume was tested in our previous simulations with overlap quarks [12,16] (see also [23,24]), where the signal was extracted from finite volume effects, which have some sensitivity to χ t [25,26]. In this work, we utilize a different method, which was originally proposed by Bietenholz et al. [27,28] (similar methods were proposed in [29] and [30]). 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 [31] 1 . 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.
We also employ a modern technique, the Yang-Mills (YM) gradient flow [33][34][35], 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, our data of the topological susceptibility is consistent with the ChPT prediction (1), from which the values of chiral condensate and l ′ are extracted.
The same set of data was also used to calculate the η ′ meson mass [36], 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 1 See also Ref. [32] where a similar method to ours was attempted.
gauge field: it creates the η ′ 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 β 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 β 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.

Lattice set-up
In the numerical simulation 2 of QCD, we use the Symanzik gauge action and the Möbius domain wall fermion action for gauge ensemble generations [37][38][39]. We apply three steps of stout smearing of the gauge links [40] 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 [41] where we use the reference YM gradient flow-time t 0 , defined by t 2 E | t=t 0 = 0.3 [33] with the energy density E of the gluon field.  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 ∼ 1 MeV  Table 1 Parameters of the JLQCD gauge ensembles used in this work. Pion masses are rounded to the nearest 10 MeV.
On generated configurations, we perform 500-1640 steps of the YM gradient flow (using the conventional Wilson gauge action) with a step-size ∆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 [10,43]. 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 [44]. 3 Topological susceptibility in a "slab" We use the conventional gluonic definition of the topological charge density q lat (x), the so-called clover construction [9]. 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 the whole sites gives values close to integers, as shown in Fig. 1.
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, as shown in Fig. 2. At the highest β = 4.47, Q lat drifts very slowly with autocorrelation time of possibly O(1000). It is, therefore, not feasible to estimate χ 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 lat , we attempt to 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 compute the correlation functions in them. After some trials and errors, we find that the so-called "slab" method, proposed by Bietenholz et al. [27] 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 6 in the same sub-volume: 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 Here T ref is an arbitrary reference time. Due to the translational invariance, the slabs sharing the same thickness T cut are physically equivalent, and one can average over This method is statistically more stable than the other sub-volume method we applied in [12,16] because 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 just 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 [26]: where C is an unknown constant, and m 0 is the mass of the first excited state, the η ′ meson 3 .
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 ) = 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.
The original proposal in [27] 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 [12,16]. 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 β At β = 4.17, which corresponds to the lattice spacing a ∼ 0.08 fm, both the global topological charge Q lat and Q 2 slab (T 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 Q 2 lat /V . In Fig. 3, 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 a linear plus constant function given in (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, 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, χ t = m ud Σ/2, with Σ = [270MeV] 3 (solid line) is drawn as an eye guide. In the same plot, we also plotted the estimate for χ t obtained from the global topological charge by circles, which again agree with the results, validating the slab method. The values of χ slab t are listed in Table 2. How we estimate their errorbars is explained in the following two sections. 5 Results at high beta At higher β values, we still find a reasonable slope at the lightest quark mass for each β 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 [10] for example), which reported that the heavier pion mass ensembles show the longer auto-correlation of the topological charge, and the show the reference points t 1 and t 2 taken for determination of the topological susceptibility. Note that the value of t 1 = 20 is the same for the two data.
larger deviation of Q lat from zero. We determine the reference t 1 ∼ 1.6 fm using 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. More details are presented in Appendix B.
Our results are summarized in Fig. 7 (see also Fig. 8 for a comparison with Ref. [8] and Ref. [9]). Although the data at higher β are rather scattered compared to those at β = 4.17, they can be used to estimate the chiral condensate Σ, assuming the linear suppression around the chiral limit. Before going to the details, we discuss the auto-correlation of χ slab t 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 where τ denotes the Monte Carlo time, and the average · · · τ is taken over τ .
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   where the upper end of the summation window W is chosen to where ρ(∆τ ) becomes consistent with zero within the error. The above formula assumes that Γ O (∆τ ) 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 ρ(∆τ ). 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 [10] carefully studied the effect of slow modes, and proposed an improved estimate of the auto-correlation time, where τ ′ int is the same summation as (7) but with a smaller upper bound W ′ where ρ(W ′ ) becomes lower than 3/2 standard deviations. τ exp is the auto-correlation of the slowest mode. The proposal is equivalent to considering a continuation of Γ O (∆τ ) at ∆τ = W ′ to the slowest possible exponential function Γ O (W ′ ) exp(−(∆τ − W ′ )/τ exp ).
In lattice QCD simulations, it is natural to assume that τ exp is equal to the autocorrelation of the global topological charge. In our simulations, τ exp is estimated by τ int (W )    Table 2 and the autocorrelation function ρ(∆τ ) at three different β with a similar pion mass M π ∼ 300 MeV is shown in Fig. 9. At the highest β = 4.47, it is clear that χ slab t has shorter auto-correlation time than that of the global topological charge.
With the measured improved auto-correlation time τ imp , we estimate the statistical errors of χ slab t by multiplying 2(τ imp + ∆τ imp )/τ interval to the naive error estimates, where ∆τ imp is the standard deviation of τ imp and τ 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.  Table 2 Results for the pion mass M π , decay constant √ 2F π , τ exp , τ imp , and χ slab t . τ exp for β = 4.47 is estimated from the first zero-crossing point of Q lat . All the data are shown in the lattice units. For χ slab t , the first error denotes the statistical error, while the second shows the systematic error due to the effect of freezing global topological charge. Auto-correlation function ρ(∆τ ) of χ slab t (pluses) and Q 2 lat (crosses) at three different β with a similar pion mass M π ∼ 300 MeV.

Chiral and continuum limit
where the renormalization factor Z S is nonperturbatively computed in [46]: Z S = 1.037, 0.934, and 0.893 for β = 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 β.
First, we compare our results directly to the ChPT formula (1). We perform a twoparameter (Σ and l) fit to the data at β = 4.17 (solid curve in Fig. 7) and β ≥ 4.35 (dashed curve) separately 4 . The results for Σ 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 β, 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 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 π and decay constant F π together with χ slab t , obtained from each ensemble, we take the ratio given in (2). By a linear one-parameter fit, 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. cWe conclude that all these systematic effects are negligibly small compared to the statistical and systematic errors given above.

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 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  Table 3 Our results for Σ, l, l ′ , and χ slab t /(M π F π ) 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.
Lattice data

Fig. 10
Continuum limit of Σ (we also plot our recent result [47] obtained from the Dirac eigenvalue density), l, l ′ and χ slab t /(M π F π ) 2 estimated by a constant fit at the physical point.
and continuum limits, respectively. The value of Σ is consistent with our recent determination  . The data at M 2 π /F 2 π > 15 are not included in the fit.
through Dirac spectrum [47]. We have also estimated the NLO coefficient where l is renormalized at the physical pion mass, while l ′ is renormalization invariant. It is interesting to note that l and l ′ include a combination of the coefficients h r 1 − h r 3 , which is supposed to be unphysical in ChPT unless θ dependence is considered. These are important for possible couplings of QCD to axions [48].
We thank T. Izubuchi, and other members of JLQCD collaboration for fruitful discussions. We also thank the Yukawa Institute for Theoretical Physics, Kyoto University.

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 [3,4] and the formula for N f = 3 is wherem = m ud m s /(2m s + m ud ), M π ,M K , and M η are the (simulated) pion, kaon and η meson masses, respectively. We have also used notations for M 2 ss = 2m s Σ/F 2 π , andM 2 = 2mΣ/F 2 π . The chiral logarithm is expressed by where µ sub denotes the renormalization scale, and g 1 is finite volume correction (see [4] for the details). In the above formula, we can see three NLO low-energy constants [45]: L r 6 and L r 8 are those renormalized at µ sub , while L 7 is a renormalization scheme independent constant.
One-loop corrections to the pion mass and decay constant were computed in [45]: and where M and F are the tree-level mass and decay constant, respectively. Now let us take the ratio of χ t and M 2 π F 2 π . Notinḡ we obtain where both of strange quark effect, as well as finite volume effects from one-loop diagrams are absorbed in the (re)definition of 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 χ t at the physical point.

B Bias from global topology
In this appendix, we discuss systematics due to freezing of the global topological charge. Combining the formulas in [27] and [26], the slab topological charge squared at fixed topology of Q becomes for 0 ≪ T cut ≪ T . Therefore, if the global topological charge Q were badly sampled and its average of Q 2 in the ensemble deviated from χ t V , we should have a quadratic term in T cut as Q 2 slab (T cut ) biased = (χ t V ) where · · · biased denotes the estimate obtained from a biased sampling of configurations.
Here we have included the term Q 2 biased , which comes from the use of the subtracted operator q lat − Q/V biased in our numerical analysis.
If the correction T cut (χ t V )T Q 2 biased − Q 2 biased − χ t V is small, our original linear + constant formula is still valid. As the correction is proportional to T cut T , if we have a window T cut ≪ T , or the freezing Q 2 biased − Q 2 biased happens to be near the true value of χ t V , we can still extract χ t from the linear slope (this seems to happen on the data at β = 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 , t ′ 2 ) = (t 1 , t 1 +t 2 2 ), and ( t 1 +t 2 2 , t 2 ) 5 , and 2) those obtained without the subtraction of Q /V in the definition of the topological charge density. Then we take the larger deviation as the systematic error. Since a part of Q 2 biased is expected to be canceled by Q 2 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 β = 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 Q 2 = 12(4) in the former half (0-2500 MD time) of the simulation time, is quite different from Q 2 = 40 (17) in the latter half (2500-5000 MD time). But the obtained values of χ slab t show a milder deviation, 1.30(53) × 10 −6 for the former half and 1.89(64) × 10 −6 for the latter, which are consistent within errors. This analysis 6 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 χ 2 /d.o.f. ∼ 1.4 also supports our conclusion.