{'fraction_non_alphanumeric': 0.029333666552122706, 'fraction_numerical': 0.010121520727250009, 'mean_word_length': 4.515506547208822, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Topological susceptibility with the asqtad action

Chiral perturbation theory predicts that in quantum chromodynamics (QCD), light dynamical quarks suppress the gauge-field topological susceptibility of the vacuum. The degree of suppression depends on quark multiplicity and masses. It provides a strong consistency test for fermion formulations in lattice QCD. Such tests are especially important for staggered fermion formulations that lack a full chiral symmetry and use the"fourth-root"procedure to achieve the desired number of sea quarks. Over the past few years we have measured the topological susceptibility on a large database of 18 gauge field ensembles, generated in the presence of 2+1 flavors of dynamical asqtad quarks with up and down quark masses ranging from 0.05 to 1 in units of the strange quark mass and lattice spacings ranging from 0.045 fm to 0.12 fm. Our study also includes three quenched ensembles with lattice spacings ranging from 0.06 to 0.12 fm. We construct the topological susceptibility from the integrated point-to-point correlator of the discretized topological charge density F-Fdual. To reduce its variance, we model the asymptotic tail of the correlator. The continuum extrapolation of our results for the topological susceptibility agrees nicely at small quark mass with the predictions of lowest-order SU(3) chiral perturbation theory, thus lending support to the validity of the fourth-root procedure.


I. INTRODUCTION
The rich topological structure of the QCD vacuum is known to be responsible for many interesting nonperturbative effects, such as the chiral anomaly and chiral symmetry breaking, instantons, and the large mass of the η ′ meson. Among the wide variety of ways of looking at these phenomena, one may consider the effect that topological charge has on the kernel of the Dirac operator. It has broad implications. For example, it is intimately connected with the value of the chiral condensate [1].
The topological susceptibility χ t characterizes the tunneling rate between topologically distinct vacua by instantons and shows up in low energy phenomenology through the Witten-Veneziano formula [2,3] and in chiral perturbation theory. A gauge configuration with topological charge ν requires at least ν fermionic zero-modes of the Dirac operator. The effect of quark mass on the topological susceptibility can be seen by separating the fermion determinant for a particular gauge field configuration into zero and non-zero modes. For N f flavors we have [1,4] where λ is the imaginary part of the eigenvalue of D / . Thus gauge configurations of nontrivial topology tend to be suppressed as any one of the quark masses approaches zero.
However, this effect is compensated at increasing volume by a growing tendency of gauge field fluctuations to produce nontrivial topology. Chiral perturbation theory tells us [1] that the outcome of the competition is controlled by the parameter x = V Σm ′ , where Σ is the chiral condensate, V is the Euclidean space-time volume, and m ′ is the reduced mass 1/m ′ = 1/m 1 + 1/m 2 + . . . .
When at least one quark mass gets small at fixed volume (the "epsilon" regime, x ≪ 1), gauge configurations with nontrivial topological charge are strongly suppressed, as implied by Eq. (1). In the physical regime, in which x ≫ 1, which is the case for our study, topologically nontrivial configurations are not suppressed. Instead, leading order chiral perturbation theory requires that the mean squared topological charge be equal to the parameter x: where the angle brackets represent an average over gauge fields. Thus the topological susceptibility, remains finite in the large-volume limit. Even so, it is still suppressed as m ′ → 0.
While lattice simulations of QCD have enjoyed considerable success in recent years, with errors on hadronic spectroscopy computations at the 1-2% level, simulations have struggled to reproduce this dependence of χ t on both m f and N f , until recently. This progress has come with improvements in lattice fermion technology, which has given much more control over chiral symmetry and lattice artifacts.
In this article we present results for the dependence of χ t on the quark mass (through the taste-singlet pion mass) using improved staggered fermions (asqtad formulation). Descriptions of the asqtad formulation have been given elsewhere [5]. To eliminate contributions from unwanted fermion doublers, the staggered formulation takes the fourth root of the fermion determinant 4 det[D / + m f ] for each quark ("fourth-root procedure"), which may raise questions about flavor counting. For a discussion of the issues, please see [5] and references therein. The primary purpose of our study, then, is to test the ability of the fourth root procedure to produce the correct number of sea quarks. Since the topological susceptibility is measured directly on the gauge field configuration without the involvement of valence quarks, it is directly sensitive to sea quark effects. We will show that the continuum extrapolation of our results agrees well with lowest-order SU(3) chiral perturbation theory.
This article summarizes results of calculations carried out over the past few years on ensembles with (2 + 1) flavors of asqtad quarks as they were being generated (see the Appendix). We continue to use the methodology of our previous work at larger lattice spacing and quark mass [6][7][8] with some refinements which appear here. The key features of our approach are these: 1. obtaining the square of the topological charge from the integral of the two-point correlator of the topological charge density.
2. reducing the variance of the integral by modeling the asymptotic form of the correlator in terms of known hadronic contributions, and 3. analyzing the quark-mass and lattice-spacing dependence of the resulting susceptibility in terms of predictions from rooted staggered chiral perturbation theory.
In the following section, we discuss the details of our method for calculating the topological susceptibility on the lattice. We present our results and analysis in Sec. III. Finally, we comment on our results in the Conclusions, Sec. IV. The Appendix lists the parameters of the gauge field ensembles used in this study.

II. METHODOLOGY
A. Definition of the topological susceptibility We introduced the topological susceptibility in Eq. (4) as the mean squared charge per unit volume: χ t = ν 2 /V . The net topological charge ν is the integral over Euclidean space-time of the topological charge density, The susceptibility is then given by the integral of the correlator of the charge density, provided the integral is well defined.
where r = |x|. Because the exponential decay of the correlator at large r is set by nonzero hadron masses, we see that the susceptibility is properly regarded as a local observable, i.e., it can be defined in terms of a correlator that has finite physical range. We use this definition of the susceptibility, coupled with a smeared lattice discretization of F F .
In the continuum limit the integral definition above is problematic. The unregulated correlator C(r) is nonintegrable: it has a positive, divergent contact term (at the origin) and, close to the origin, a compensating negative ultraviolet singularity of order (up to possible logarithms) r −8 [9][10][11][12]. Cancellation is required in order to produce the expected finite integral Eq. (4). To circumvent this mathematical difficulty Lüscher formulated a definition of the topological susceptibility in terms of a product of pseudoscalar and scalar densities of Ginsparg-Wilson quarks [13]. Since the definition requires computing all-to-all quark-line disconnected correlators, it is more difficult to implement, and, to our knowledge, it has not yet been put into practice.
For present purposes we resort to the naive definition in Eq. (6) and trust that the lattice cutoff and a smoothed definition of ρ(x) regulate the compensating singularities enough over a range of reasonably small lattice spacings that we can test the expected suppression of the topological susceptibility. In our present scheme we fix the smoothing scale in lattice units as we take the lattice spacing to zero. Our numerical simulation provides a practical test of the limitations of such a scheme. If it fails, as the lattice spacing is decreased, we would expect to encounter a growing variance from contributions to the integral near the origin. This would not invalidate the method: The central value in the continuum limit is finite even if the variance is unbounded. It could, however, require an impractically large computational effort to achieve a desired accuracy as we make the lattice spacing smaller.
We return to this question after presenting our results.
There are a variety of lattice methods for obtaining the topological charge. The traditional "algebraic" method uses a lattice discretization of the density F F , constructed at each lattice site from appropriate closed loops of gauge links. To suppress ultraviolet noise at the cutoff scale, smoothing is required [14]. The Boulder discretization [15,16], which we use in the present study, is a refinement of the traditional definition. It is fattened (smoothed) by first performing some number (we use three) of HYP smoothing sweeps [17] on the gauge field and then constructing the operator from the smoothed links.
A more elegant method defines the topological charge density in terms of a chiral (e.g. overlap) Dirac operator D, as ρ(x) ∝ Tr[γ 5 D] x,x [18,19] (the trace is over color and spin), but using it directly in Eq. (6) is computationally expensive [20,21]. For the overlap operator a more tractable method uses the Atiyah-Singer index theorem to relate the topological charge ν to the net number of zero crossings of the low-lying eigenvalues of a Hermitian Dirac kernel from which the chiral operator is built [22]. This method was implemented in [23]. For the overlap operator, smoothing is inherent in the choice of the Dirac kernel from which the overlap action is built.
Another promising method works with gauge configurations of fixed topology [24,25].
In this case, at large distance the correlator of the topological charge density approaches a constant χ t /V plus other known constants that depend on the fixed topological charge.
One can also use a hadronic flavor-singlet interpolating operator with J P C = 0 −+ as a proxy for F F . This method has been tested at one lattice spacing in the two-flavor case on configurations generated with the overlap action [26].
The Lüscher definition [13], based on a chiral Dirac operator, replaces the integral of F F with the integral of a quark pseudoscalar density. The quark field from which that density is constructed can have arbitrary mass, which sets the localization scale of the operator. The expectation value of that density is regulated with a suitable number of zero-momentum scalar-density insertions on the quark line. At large mass in the hopping parameter expansion, the operator can be expressed as a sum of gauge-link loops analogous to those in the Boulder discretization, which regulates the construction of F F through an extended discretization and HYP-smeared gauge-links. In the Boulder case the localization of the gauge paths is controlled by the number of smearing steps, whereas localization of the Lüscher operator is controlled by the quark mass. Of course, the chiral properties of the underlying action in that case allows an arbitrary choice of scale.
Whatever the definition, the resulting susceptibility is subject in general to multiplicative and additive corrections at nonzero lattice spacing [27]: An additive renormalization is not required for chiral actions that use the same operator in the fermion determinant and the measurement of the topological charge [28]. In our case an additive renormalization is expected. We assume that in the continuum limit M approaches one and A approaches zero. Since with our actions lattice artifacts appear at O(a 2 ) (up to logarithms), we expect that the approach to these limits is as a 2 [6]. With the overlap method one can use the same Dirac operator for the Monte Carlo evolution and the measurement of topological charge. In this case the small instantons and dislocations that are not seen by the overlap operator, so not suppressed by a small quark mass, are then also not seen by the topological charge operator. In our case the Monte Carlo Dirac operator and topological charge operators are unrelated, so we might expect larger lattice artifacts.

B. Predictions from chiral perturbation theory
Our computed topological susceptibility is a function of the quark masses and the lattice spacing. As we have already recalled in Sec. I, in chiral perturbation theory the susceptibility χ t depends on the number of light quarks and their masses in leading order through where Σ is the chiral condensate to this order, m u , m d , and m s are the masses of the up, down, and strange quarks, and the ellipsis represents contributions beyond the cutoff from higher quark masses and the axial anomaly [1]. We see that as quark masses vanish, the susceptibility must vanish. The rate at which it vanishes depends on the number of light flavors.
For equal up and down quark masses we may use the Gell-Mann-Oakes-Renner relation, also from leading order chiral perturbation theory, to rewrite this expression as where m 2 ss = 2m 2 K − m 2 π is the squared mass of the fictitious pseudoscalar meson containing two nonannihilating quarks with masses equal to the strange quark, and in our normalization the pion decay constant f π is approximately 130 MeV.
With nonchiral lattice fermions, at nonzero lattice spacing one should instead use a version of chiral perturbation theory appropriate to the lattice fermion formulation. In this way some of the lattice discretization errors can be modeled. For staggered fermions using the fourthroot procedure, we use rooted staggered chiral perturbation theory (rSχPT) [29]. This theory has a taste multiplet of sixteen pions. Among them, only the taste singlet pion is sensitive to the anomaly and so enters the expression for the topological susceptibility at leading order. At tree level the continuum expression is modified by replacing the pseudoscalar meson masses by their taste-singlet counterparts [30]: where the subscript I identifies the taste singlet, and through the term in m 0 , which is proportional to the η ′ mass at lowest order, we have introduced an explicit anomaly contribution. The standard chiral perturbation theory expression corresponds to m 0 → ∞ (and a → 0); introducing m 0 in Eq. (10) is phenomenological because m 0 is beyond the physical cutoff scale of chiral perturbation theory. At infinite quark mass we get the quenched topological susceptibility χ tq , which suggests an alternative phenomenological form [31], C. Topological charge density operator As before [6], we use the topological charge operator of DeGrand, Hasenfratz, and Kovacs [15] optimized for SU(3) by Hasenfratz and Nieter [16]. The operator is constructed from closed ten-link paths of gauge matrices as follows: Specifically, the operator U 1 is constructed from a product along a path from site x in the se- operator also reproduces accurately the charge of an instanton, provided the instanton radius is larger than the lattice spacing. The finer details of the construction of this operator are unimportant for our purposes, since in the end we take the continuum limit.
We applied this operator to gauge configurations smoothed by three HYP steps [17].
From the point of view of the unsmoothed gauge field, this operation, in effect, enlarges the footprint of the topological charge density operator by a small amount. We have shown in [6] that the topological susceptibility on a coarse lattice (a ≈ 0.12 fm) is constant within statistical errors of 8% for one to four HYP sweeps.

D. Variance reduction method
We calculate the topological susceptibility by integrating the topological charge density correlator in Eq. (6) over the lattice four-volume. In the left panel of Fig. 1 we show a typical correlator C(r). It is expressed in units of the Sommer parameter r 0 ≈ 0.454 fm [32]. As expected, it has a positive peak at the origin next to a negative minimum, and it rises to its asymptotic limit of zero from below as required by CP symmetry. To give a better visual impression of contributions to the susceptibility, in the right panel of Fig. 1 we multiply C(r) by the statistical weight factor w(r) that counts the number of lattice points that, by symmetry, have the same four-radius r, or, where the plotted value is binned, have the same range of four-radii. This is essentially a discretized version of r 3 C(r)dr. The irregular binning inherent in the discretized distance r produces the ragged appearance of the weighted values at small r. On the other hand, statistical fluctuations produce the ragged appearance at large r. The topological susceptibility in r 0 units is simply proportional to the sum of the weighted values.
In Fig. 1, right, the substantial cancellation of positive and negative contributions at small r is more evident. We also see that the large distance contribution to the susceptibility is mostly noise. We have found that it is responsible for the bulk of the variance in the integral. This is to be expected. In a suitably large subvolume V 0 of spacetime, we should be able to determine the topological susceptibility reasonably well by measuring fluctuations of the local topological charge ν 0 . Consider putting together N such volumes to create the total volume V . The overall topological charge ν is then obtained as a random walk of local charges, so its variance grows with N. We can measure the susceptibility in two ways: (1) average the locally determined ν 2 0 /V 0 over the N subvolumes or (2) calculate ν 2 /V over the full volume. With the former method the error in the measured susceptibility decreases as 1/ √ N with increasing N and fixed V 0 , whereas with the latter method the error never improves.
In our case the integral of the correlator C(r) replaces the sum over subvolumes. But we still need to eliminate noise from contributions at large r. To do so, several years ago we introduced a variance-reduction method that fits the large r part of the correlator to its asymptotic form in Eq. (15) and then, for r > r c for a suitable cutoff r c , replaces the numerical sum of the correlator with an integral over the fitted function as follows [7]: where r c is chosen inside the fit range. In the present study we chose r c ≈ 1.2 fm for all ensembles. In Fig. 1, right, we illustrate the fit to the large r part of the correlator and indicate r c . We continue to use this method in the present work.

E. Asymptotic fit model
The topological charge density is a flavor-singlet operator with quantum numbers J P C = 0 −+ , so the asymptotic behavior of the correlator is governed by the η and η ′ mesons and, for sufficiently light sea quarks, by multipion states. That is where the A's are overlap constants and S(m, r) is a scalar propagator with asymptotic form S(m, r) ≈ [1 + 3/(8mr)] exp(−mr)/r 3/2 .
The three-pion continuum is the lightest multimeson state in this correlator. For our ensembles the η meson is always lighter. Furthermore, the coupling of the topological charge density operator to multimeson states is Zweig-rule suppressed. Therefore, we ignore them in the present analysis. Since the topological charge density operator is an SU (3) so Our statistics are insufficient for determining all the parameters of the fit function reliably. Instead, we model the masses of the η and η ′ and the ratio A η /A η ′ , leaving only one fitting parameter A η ′ , which simply sets the normalization of the asymptotic form. We set m 2 η = 2m 2 ss /3 + m 2 π /3 for our measured lattice values of ms s and m π , and we fix m η ′ = 958 MeV (its physical value) since we have not calculated it for these ensembles. Finally, we use a simple chiral model to fix the ratio of couplings A η /A η ′ or equivalently, the singlet-octet mixing angle as a function of quark masses.
where M U I and M SI are masses of unmixedūu (dd) andss meson states, and m 2 0 parameterizes the anomaly. The isosinglet eigenvectors are where and N and N ′ normalize the eigenvectors to 1. Since the flavor singlet state in this basis is just (1, 1, 1)/ √ 3, we obtain the ratio which we apply to the fit model of Eq. (14). To complete the model, we need the value of the anomaly parameter m 2 0 . We set it so that for physical values of M 2 U I and M 2 SI (i.e., values that give the physical masses m π and ms s = 2m 2 K − m 2 π ), we get the standard phenomenological mixing angle θ ≈ −20 degrees [33]. At this "physical" point the mixing model also gives us m η = 493 MeV and m η ′ = 953 MeV, reasonably close to their physical values. Then for unphysical masses we use the lattice values of M 2 U I and M 2 SI on each ensemble, always keeping m 0 fixed. This procedure assures that the η decouples as required in the SU(3) flavor limit m u = m d = m s , and it provides a smooth interpolation between that limit and the physical limit. The taste-singlet masses M 2 U I and M 2 SI are obtained by adding measured or estimated taste splittings to the masses of the lightest members of the taste multiplet. Splittings are listed in Table V below.
The model is applied to all the dynamical ensembles in this study, listed in the Appendix A. The resulting fit parameters are listed in Table I. The mixing parameter A η /A η ′ is shown to three digits. Apart from systematic errors in the model itself, in principle it inherits a statistical error from our determination of the taste-singlet masses, which, in turn depends on the error in the taste splitting. The last error, however, is less than 5%, small enough to have no effect on the mixing parameter to the number of digits reported. The remaining fit parameters do not depend on the taste-singlet masses. Consequently, statistical errors in the determination of the taste-singlet masses have negligible effect on results for the topological susceptibility.

F. Asymptotic fit model for the quenched ensembles
For the three quenched ensembles we use the same methodology, except that the fit model has only one mass. We fix it to the mass of the J P C = 0 −+ ground state lattice glueball from Chen et al. [34], namely 2560 MeV. The parameters are listed in Table II. We chose r c for the quenched ensembles to match our choice for the dynamical ensembles at the same lattice spacing. Since the quenched correlators die so quickly at large r, the contribution to the susceptibility for r > r c is negligible, and the asymptotic model has no effect on the result.

III. RESULTS
We smooth the lattices with three HYP smoothing steps [17] and measure the topological charge density with the Boulder operator at each space-time point. We then construct the point-to-point correlator C(r) for every pair of points in the space-time volume. For r/a < 5 we keep values for every displacement, and for larger r we bin data over small intervals in r. The resulting data is then fit to Eq. (14) over a range [r min , r max ]. We replace the raw data with the fit model for r > r c . The fit range is chosen to give an acceptable χ 2 /df (corrected for autocorrelations) and to vary smoothly as a function of sea quark mass and lattice spacing. As we have noted, however, the topological susceptibility is a local observable. We can get a graphical sense of the autocorrelation affecting the susceptibility by considering the time history of the integral of the correlator In Fig. 3 we show the time history of this variable for the case r = 2r 0 for the same set of ensembles. Clearly the fluctuations in this quantity decorrelate much more rapidly than those of the total topological charge.
We estimate the autocorrelation correction, i.e., the amount by which the naive (uncorrelated) variance should be increased to compensate for autocorrelations. For this purpose we consider the integral of the correlator over the proposed fit range rmax r min We block the data in Monte Carlo time and calculate the variance of the mean as a function of block size, extrapolating to infinite block size. The ratio of the extrapolated variance to the naive variance is the correction factor. We also sum the autocorrelation coefficients to obtain another estimate of the correction factor. These determinations fluctuate as a function of sea quark mass. We averaged them at fixed lattice spacing to obtain the correction factors shown in Table III. We should emphasize that the determination of autocorrelation corrections is notoriously difficult. To develop more confidence in these estimates, we should have considerably longer time histories.
Our fits to the data take into account correlations in r as well. For all ensembles, measurements are taken every six or sometimes every five molecular dynamics time units. We do not bin data in Monte Carlo time before constructing the covariance matrix in r and minimizing the correlated χ 2 [35]. Uncorrected errors are derived from a jackknife analysis.
Thus the resulting χ 2 , based on the naive covariance, must be reduced by the factor in Table III before estimating the confidence level. Furthermore, the naive single-elimination jackknife errors in the fit parameters must be increased by the square root of this factor. We use the same factor to adjust the error in the contribution from the raw data for r < r c .

B. Topological charge density correlator
We expect the topological susceptibility to decrease with decreasing light sea quark mass.
It is interesting to see how the topological charge density correlator itself varies with the light sea quark mass at fixed lattice spacing. In Fig. 4 we examine this dependence for a series of fine lattice ensembles (a ≈ 0.09 fm) for which we have results for four ratios of the light to strange quark mass, 0.05, 0.1, 0.15, and 0.2, corresponding to the range 0.601 to 1.074 in m 2 πI r 2 0 . In the upper panel any variation with light quark mass is evidently much smaller than the plot symbol size. In fact the short distance part of the correlator shows very little sea quark mass dependence. In the lower panel we enlarge the region around the minimum where a small variation is now apparent. In this region light meson states begin to dominate the correlator of the gluonic operators. As the light quark mass decreases, the minimum drops, thus giving a larger negative contribution to the integral. This effect leads to the suppression of the susceptibility. According to the model, the correlator should also decay more slowly at large r, but this effect is too subtle to be visible with our statistics.
We next examine the lattice spacing dependence of the correlator at fixed light quark mass ratio. Comparing the local correlators C(r) obtained on ensembles at different lattice spacing is complicated because sampling is naturally done on a lattice scale. Rather than rebinning the data to a common physical scale, we compute the partial integral χ t (r) of Eq. (22) and plot it in physical (r 0 ) units in Fig. 5. As r increases from the origin, we see a peak at short distance coming from the regulated contact term followed by a decrease coming from the negative correlator. The onset and width of the peak is determined by the effective radius of the topological charge density operator, which is fixed in lattice units.
Thus as the lattice spacing decreases, the expected negative 1/r 8 singularity in the correlator is exposed, and the peak increases in height and decreases in width.
At large r the data approach the asymptotic value of the full susceptibility.  Table IV. This statistical error is adjusted for autocorrelations, sample size (factor of N/N 0 ), and lattice volume (factor of V /V 0 ) for N 0 = 500 and V 0 = 100 fm 4 . We see that the adjusted error grows approximately as 1/a over this range. This trend suggests that it will be increasingly expensive to push to smaller lattice spacing with our scheme. However, the continuum limit is nonetheless finite, and our results demonstrate that the method gives reasonable errors over the range of lattice spacings considered.

C. Topological susceptibility
Our results are summarized in Table VI and Fig. 6. Since chiral perturbation theory predicts the behavior as a function of the mass of the taste-singlet pion, we also list estimates of that mass. Unlike the Goldstone pion mass, the mass of the taste singlet is not measured directly on all of our ensembles. However, to a good approximation, splittings of the squared masses of the pion taste multiplet are known to be independent of the light quark mass at  [36]. So if the splitting is measured for one light quark mass for a given lattice spacing, the taste-singlet pion mass can be reconstructed from the Goldstone pion mass for other light quark masses at the same spacing. Table V lists the splittings for the categories of lattice spacings in this study. They were used to obtain the values in Table VI.
The largest error in the estimated splittings is less than 5%, which bounds the error in the abscissa of the plot. We have chosen r c to lie within the fit range. We have found that within this range our results vary by less than one standard deviation.

D. Continuum extrapolation
To model a continuum extrapolation, we fit our data to the following form: This model assumes that lattice artifacts scale as a 2 . The fit yields χ 2 /df = 8.8/13. In Fig. 6 the resulting fit curves are shown, and three representative points in the continuum extrapolation are also plotted. Also plotted is the prediction of Eq. (11) using f π r 0 = 130 × 0.454 MeV-fm with and without our continuum-extrapolated asymptotic quenched topological susceptibility χ t r 4 0 = 0.0523 (29). From the fit itself we obtain f π = 132(6) MeV, which is better than expected for tree-level chiral perturbation theory.

IV. CONCLUSIONS
We have presented an extensive study of the topological susceptibility on 18 (2 + 1)flavor asqtad lattice ensembles and three quenched lattice ensembles. The susceptibility is  Table VI.) The fit curves for a = 0.06 and 0.045 fm are, accidentally, nearly coincident.
Statistical errors on the solid lines are shown on the right edge of the right panel. They have also been corrected for autocorrelations. The fit error for the smallest lattice spacing has the largest error bar. The improvement in variance is evident. defined as the integral of the correlator of the topological charge density. The topological charge density is constructed from a discretized version of F F with smearing to help regulate ultraviolet fluctuations. To reduce the variance from large distances, we replace the measured values of the correlator at large r by a fit model that builds in the expected spectral contribution.
Our method for determining the topological susceptibility through an integral of the topological charge density correlator avoids singularities at zero separation by smearing the charge density operator over a fixed local set of lattice sites. A study of the variance in the small-distance contribution suggests that as the lattice spacing is decreased the variance grows. At our level of statistics and for the range of lattice spacings we consider in this study, this growth is manageable.
Over the range of lattice spacings and masses in this study, within statistical errors, we find good agreement with tree-level staggered chiral perturbation theory and in the continuum limit with tree-level continuum chiral perturbation theory, in both cases with the expected number of flavors. This agreement supports the assertion that the fourth-root procedure for staggered fermions results in the correct number of sea quark species in the continuum limit.