Stochastic calculation of the Dirac spectrum on the lattice and a determination of chiral condensate in 2+1-flavor QCD

We compute the chiral condensate in 2+1-flavor QCD through the spectrum of low-lying eigenmodes of Dirac operator. The number of eigenvalues of the Dirac operator is evaluated using a stochastic method with an eigenvalue filtering technique on the background gauge configurations generated by lattice QCD simulations including the effects of dynamical up, down and strange quarks described by the Mobius domain-wall fermion formulation. The low-lying spectrum is related to the chiral condensate, which is one of the leading order low-energy constants in chiral effective theory, as dictated by the Banks-Casher relation. The spectrum shape and its dependence on the sea quark masses calculated in numerical simulations are consistent with the expectation from one-loop chiral perturbation theory. After taking the chiral limit as well as the continuum limit using the data at three lattice spacings ranging 0.080-0.045 fm, we obtain $\Sigma^{1/3}$(2 GeV) = 270.0(4.9) MeV, with the error combining those from statistical and from various sources of systematic errors. Finite volume effect is confirmed to be under control by a direct comparison of the results from two different volumes at the lightest available sea quarks corresponding to 230 MeV pions.


Introduction
Spectrum ρ(λ) of the eigenvalues λ of the Dirac operator D in Quantum Chromodynamics (QCD) reflects the properties of background gauge field. At zero temperature, pairs of quark and antiquark condense in the vacuum as represented by the Banks-Casher relation ρ(0) = Σ/π [1], which is valid in the thermodynamical limit, i.e. massless quark limit after taking an infinite volume limit. In other words, the density of near-zero eigenvalues of the Dirac operator is related to the chiral condensate Σ, which is an order parameter of spontaneous chiral symmetry breaking in QCD. In the chiral effective theory, for which pions play the role of effective degrees of freedom of QCD at low energy, the chiral condensate Σ and pion decay constant F π are the most fundamental parameters appearing at the lowest order in an expansion in terms of pion mass and momenta. The QCD Dirac spectrum can thus be related to physical observables involving pions at low energy.
The chiral effective theory predicts the functional form of ρ(λ) in the low energy regime.
In the limit of infinite volume, the slope of ρ(λ) at λ = 0 was calculated including the loop effect of pions [2], and the dependence on the number of dynamical quark flavors was predicted. In a finite volume, the lowest end of the spectrum is largely affected and exact zero-modes play a special role. Such system is related to the chiral Random Matrix Theory (RMT), with which the distribution of individual eigenvalue can be calculated [3]. (For more results, see a recent review article [4].) The most elaborate calculation to date includes finite volume and finite quark mass corrections in a systematic expansion [5].
In lattice gauge theory calculations, the spectral density has so far been calculated by direct computation of the low-lying eigenvalues or by stochastic estimates of the mode number below some value [6]. The direct computation of individual eigenvalues has an advantage of allowing a comparison of the microscopic distribution with that predicted by chiral RMT. Even with a few lowest eigenvalues, one can then extract Σ assuming the correspondence between the chiral effective theory and the random matrix theory. In our previous works using the overlap fermion formulation, we studied the quark mass and volume dependence of the eigenvalue distribution and extracted the value of Σ in 2-flavor [7,8] and 2+1-flavor QCD [9,10]. Since the overlap fermion preserves exact chiral symmetry, the smallest eigenvalues satisfy the relations derived from chiral symmetry, and the correspondence between the non-perturbative lattice calculation and the analytic prediction of the effective theory and chiral RMT [11,12] has been precisely established.
In order to achieve precise calculation of the physical value of Σ, on the other hand, the direct eigenvalue calculation with the exactly chiral fermion formulation is computationally 2 too expensive. Finite volume effect and discretization effect are best controlled by calculating on sufficiently large and fine lattices. The number of relevant low-lying eigenvalues to be calculated grows as the (four-dimensional) volume V , and the computation of individual eigenvalues rapidly becomes impractical. The stochastic estimate introduced in [6] offers an alternative method in such situations. The method has been successfully applied to extract Σ in 2-and/or 2+1-flavor QCD with Wilson [13,14] and twisted-mass [15] fermion formulations.
In this work we use a slightly different implementation of the stochastic estimate. It is based on a filtering of eigenvalues in a given interval [16]. The method allows us to estimate the number of eigenvalues in any interval once the necessary coefficients have been calculated.
We use the domain-wall fermion formulation, with which chiral symmetry can be maintained at the level that the effective residual quark mass is of order of 1 MeV. We design the eigenvalue filtering such that the number of eigenvalues in a bin of 5 MeV or larger is counted and the possible effect of the residual chiral symmetry violation is harmless.
We calculate the eigenvalue spectrum on the lattices generated with 2+1 flavors of light sea quarks described by the Möbius domain-wall fermion. Sea quark masses in the simulations correspond to the pion mass in the range of 230-500 MeV. Physical volume is sufficiently large, L ∼ 2.6 fm or larger, in order to safely neglect the effect of finite volume which affects the lowest eigenvalues of order λ ∼ 1/(ΣV ) (∼ 1-2 MeV) most strongly while the number of eigenvalues below 10-20 MeV is little affected. The finite volume effect due to the loop effects of light pions is suppressed as exp(−M π L), and is sufficiently small on our lattices satisfying M π L > 4.
Our lattice ensembles are in a range of lattice spacing a between 0.080-0.044 fm. The corresponding lattice cutoff a −1 ranges between 2.45 GeV and 4.50 GeV. On these fine lattices, the discretization effects for the near-zero eigenvalues of order 10 MeV should be negligible. Indeed, we found that the scaling violation is consistent with zero for the spectral function.
The Möbius domain-wall fermion is an (approximate) implementation of the Ginsparg-Wilson relation [17]. The residual mass with our parameter choices is O(1 MeV) or less strongly depending on the lattice spacing, and its effect on the calculation of the eigenvalue spectral density is minor.
Using these data sets we obtain the spectral density, which we then fit with the formula predicted by the chiral effective theory to obtain the value of chiral condensate Σ in the chiral limit of up and down quarks.
The rest of the paper is organized as follows. In Section 2 we review the method of the eigenvalue filtering and the stochastic eigenvalue counting. Section 3 summarizes the lattice 3 fermion formulation, which is followed by the details of our data sets in Section 4. The spectral function in the entire range of eigenvalues is shown in the plots given in Section 5. We then focus on the low-lying eigenvalue spectrum to extract the low-energy constants including the chiral condensate using the chiral perturbation theory, as described in Section 6. Our conclusion is in Section 7. A preliminary report of this work is found in [18].

Stochastic estimate of eigenvalue count
We review the method to evaluate the eigenvalue count of a hermitian matrix in a given interval. More details are described in [16]. In the lattice gauge theory calculations, the method is introduced recently in [19].
Let A be a hermitian matrix and assume that its eigenvalues are distributed in the range [−1, 1]. If not, we can easily rescale the matrix by a linear transformation. We aim at calculating the number of eigenvalues of this matrix in a given interval [s, t]. By introducing a step function h(A) that has a value 1 only in the interval [s, t] and zero elsewhere, the number of eigenvalues is written as n[s, t] = Tr h(A). Then, introducing N v Gaussian random in the limit of large N v . This evaluation can be promoted to the ensemble average as where · · · represents an average over Monte Carlo samples, or the gauge configurations.
With sufficiently large number of gauge configurations, we may even take N v = 1 to obtain a statistically significant signal.
The discrete function h(A) may be constructed approximately using a polynomial function even when the matrix A is large. The best approximation of h(x) in the sense of min-max (smallest maximum deviation) in the interval x ∈ [−1, 1] achieved within a given computation cost is the Chebyshev approximation using the Chebyshev polynomial T j (x). Explicitly, we may write with coefficients γ j , which can be calculated as a function of s and t. See equation (7) of [16], which is reproduced below for convenience: sin(j arccos(s)) − sin(j arccos(t)) j j > 0.
In order to suppress a strong oscillation emerging with this approximation, the so-called Jackson damping factor g p j is introduced, sacrificing the "best" approximation. An explicit formula, equation (10) of [16], is where α p = π/(p + 2).
The formula (3) can then be modified as Using this form, the stochastic estimate of (2) can be approximated as This approximation is convenient, because one can obtain the eigenvalue count in any range [s, t] once we have the set of measurements for ξ † k T j (A)ξ k . The Chebyshev polynomial is constructed using the recursion relation: T 0 (x) = 1, T 1 (x) = x and There is also an useful formula, 2T m (x)T n (x) = T m+n (x) + T |m−n| (x), which in particular reads One can then apply A on ξ k repeatedly to obtain ξ † k T j (A)ξ k . Note that the 2n-th order is obtained from (T n (A)ξ) † (T n (A)ξ) using the formula above. One therefore needs n multiplication of A to obtain the order of polynomial p = 2n.
The accuracy of the approximation depends on the order of the polynomial p. The size of error is discussed in the next section for the application to the spectral function of the domain-wall fermion Dirac operator.

3 Domain-wall Dirac operator
In this work we utilize the Möbius domain-wall fermion formulation [20] to define the Dirac operator on the lattice. It is a generalization of the domain-wall fermion [21,22] introduced to achieve better chiral symmetry within a given computational cost. In this fermion formulation, the fermion field is defined on a five-dimensional (5D) lattice, and a four-dimensional (4D) fermion emerges on the 4D surfaces of the 5D space. The fermion modes of right-handed and left-handed chiralities localize on the opposite 4D surfaces, and thus chiral fermion is realized with exponentially suppressed violation as a function of the extent in the fifth direction L s .
The effective 4D Dirac operator D (4) is constructed combining the 5D Dirac operator D Here P is a certain permutation operator acting on the fifth coordinate s designed to move the physical surface modes (both left-handed and right-handed) to the slice of s = 1. The suffix "11" then means to extract that 4D slice. The term (D GDW (1)) −1 implies an introduction of a Pauli-Villars field, which cancels unnecessary 5D modes in the ultraviolet limit.
The 4D operator D (4) approximately satisfies the Ginsparg-Wilson relation [17] and the eigenvalues of the hermitian operator D (4) † D (4) are constrained in the range [0, 1]. In order to apply the eigenvalue filtering method described in the previous section, we therefore define such that A has eigenvalues between −1 and 1.
The low upper limit (=1) of the eigenvalue of D (4) † D (4) is one of the advantages of using the domain-wall fermion. With the Wilson fermion formulation, for instance, the highest eigenvalue is 8 2 = 64 (or slightly less for interacting cases) and one has to shrink the whole eigenvalue range by multiplying a factor ∼ 30 to fit in [−1, 1] when we map the Wilson operator on A as in (12). The target eigenvalue interval is then much narrower for A, and one needs larger polynomial order p to obtain the same level of accuracy. Although the numerical cost is higher for the domain-wall fermion due to the inversion of the Pauli-Villars operator for each application of D (4) , the difference of the entire eigenvalue range nearly compensates the cost compared to the Wilson fermion. An eigenvalue a 2 λ D † D of D (4) † D (4) can be related to that of D (4) assuming the Ginsparg-Wilson relation (11) as well as the γ 5 -hermiticity property D (4) † = γ 5 D (4) γ 5 . The relation (11) is slightly violated in the actual implementation, and the associated error is discussed later. The eigenvalues λ D of D (4) lie on a circle on the complex plane to satisfy |aλ D − 1/2| = 1/2. We project them to the imaginary axis to obtain the continuum-like eigenvalue λ as This is a convention, and other definitions such as aλ = |aλ D | are equally valid up to the discretization effect of O(a 2 ). For the low-lying modes below 20 MeV, which are the eigenmodes we use to extract the chiral condensate, the discretization error of O(a 2 ) is expected to be very small. Examples of the filtering function are shown in Figure 1 for the order of polynomial p = 8000. Here, the Dirac eigenvalue aλ as defined in (13) is taken on the horizontal axis. The plot shows the function to extract the count in the lowest bin [0, aδ] of bin size aδ = 0.01, 0.005, 0.002 and 0.001. The approximation of the step function is very precise except for the region close to the threshold aλ = aδ. The width where the function varies is nearly independent of aδ, and as a result, the relative error of the approximation is smaller for larger bin sizes. Note that the lowest eigenvalue is the worst case, because it is mapped onto a narrow bin of size 2(aδ) 2 of A.
In order to quantify the size of the error in filtering, we calculate a fraction of leakage from the lowest bin [0, aδ] to the neighboring bin. It is defined as an integral of the filtering function from aδ to infinity, which should vanish for the exact step function. The leakage is equal to the deficit in the bin of interest [0, aδ]. Figure 2 shows the leakage for various widths aδ. The relative error increases for smaller aδ as an inverse power 1/aδ. If we allow an 1% error for the calculation of the spectral function, we may take aδ to be 0.005 when p = 8000. This bin size corresponds to 12 MeV on our coarsest lattice. On finer lattices we take larger values of p so that the error with a fixed δ, which implies a smaller aδ on a finer lattice, is not larger than 0.005. Since the deficit is largely compensated by the leakage from the neighboring bin when the spectral function is nearly constant as it is the case for zero temperature QCD, the actual error would be much smaller than this naive estimate.
According to the general theory of the Chebyshev approximation, the error as measured by the L 2 norm scales as 1/ √ p [16]. With the Jackson damping factor implemented in this work, this bound does not apply, but an actual calculation as outlined above indicates that the leakage decreases as 1/p. This determines the computational cost when one wants to improve the precision using this method.

Lattice ensembles
We calculate the spectral function at three β values on 15 gauge ensembles in total, generated with 2+1 flavors of sea quarks [23], as listed in Table 1. The formulation for the sea quarks is the Möbius domain-wall fermion, which is the same for the lattice Dirac operator used in the eigenvalue counting. The gauge action is tree-level Symanzik improved, and we apply the stout link smearing [24] three times for the link variables entering the definition of the fermionic operators.  Table 1 Lattice ensembles used in the eigenvalue spectrum calculation. Spatial lattice size L/a and sea quark masses am ud , am s are listed in the lattice unit. p is the order of the Chebyshev polynomial, and N meas is the number of measurements. Empty entries are the same as the ones in the previous line. we chose to count the eigenvalues, such effect would be minor; we eventually eliminate the associated error by taking the continuum limit using the three lattice spacings we prepared.
The same set of ensembles is used for a wide variety of applications including a determination of non-perturbative renormalization constant [27], a determination of the charm quark mass from temporal moments of charmonium correlator [28], a calculation of the η ′ meson mass through a gluonic observable [29], and a calculation of D (s) meson decay constant [30]. Numerical calculations of the projects are performed using the code set IroIro++ Averaging over 50 gauge configurations each with only one noise per configuration, we calculated ξ † k T j (A)ξ k . The mode numbern[s, t] is then evaluated by summing over j from 0 to p as (6). The spectral density is obtained with an appropriate normalization, where aλ = s 1/2 /(1 − s) and a(λ + δ) = t 1/2 /(1 − t). The factor 2 in the denominator of (14) reflects the pairing of the eigenvalues, i.e. ±iλ.  Figure 3 shows the spectral density a 3 ρ(λ) in the whole range of aλ. The bin size is aδ = 0.005. One can clearly observe that the spectrum starts from a tiny constant at λ ≃ 0 and increases towards higher eigenvalues. The near-zero modes show some dependence on the sea quark mass (see below), but the high modes are nearly independent of the sea quark mass.
The increase towards the perturbative regime at high aλ is qualitatively consistent with the free-theory scaling ∼ λ 3 , but saturates at around aλ ∼ 1 due to the discretization effect. plotted. We find that they are consistent within the statistical errors. The statistical error is larger for smaller bins since the number of eigenvalues in each bin is fewer.
The volume scaling of the spectral function is demonstrated in Figure 5. For the lightest pion (m π ≃ 230 MeV), there are data on two volumes 32 3 × 64 and 48 3 × 96 available. We calculate the spectral density on both lattices with exactly the same method. The results are consistent with each other within the statistical error, which is about 5% on the 32 3 × 64 lattice. The statistical error is smaller, about 2%, on the larger volume since the number of eigenvalues in a given bin is proportional to the physical volume V . 6 Analysis with chiral perturbation theory The Banks-Casher relation ρ(0) = Σ/π is valid only in the chiral limit after taking the infinite volume limit. Therefore, the effects of finite sea quark masses and finite volume need to be taken into account in the analysis. We use the functional form predicted by the chiral effective theory to analyze the quark mass dependence. The finite volume effect is also The analytic calculation is available at the one-loop order of chiral perturbation theory (χPT), which is valid at the leading non-trivial order of finite quark mass correction, i.e. of order m 2 π /(4πF π ) 2 . The formula is concisely written in the form [5] (see also [10]) where the chiral condensate Σ and pion decay constant F are those in the chiral limit. One of the low-energy constants at the one-loop order, L 6 , appears for this quantity. The functions ∆(0, M 2 ) and G(0, M 2 , M 2 ) are given as They are evaluated at a "pion mass" as determined by the Gell-Mann-Oakes-Renner (GMOR) relation M 2 ij = (m i + m j )Σ/F 2 , where the indices i and j label the sea quark mass or a fictitious valence quark v. For the sea quark mass, it gives a leading-order estimate of the corresponding pion mass. It slightly deviates from the actual pion mass calculated on the lattice with the same quark mass, but the difference is from higher orders of the chiral expansion and thus can be neglected at the order considered for ρ(λ). The "valence quark" mass m v is taken at an imaginary value iλ to obtain the spectral function ρ(λ) at a finite λ, according to the procedure in [5]. The scale parameter µ sub denotes the renormalization scale, which is conventionally taken at the ρ meson mass.
The function g 1 (M 2 ) in (16) represents the finite volume effect and is written in terms of a sum of the modified Bessel function. In the analysis of chiral extrapolation, we ignore the contribution of g 1 (M 2 ), which is a good approximation for our data. The largest possible finite volume effect may arise for the ensemble of lightest pion with the smaller volume, i.e. the 32 3 × 64 lattice of am ud = 0.0035 at β = 4.17, for which our estimate of g 1 (M 2 )/F 2 is ∼0.05 (0.02) at λ ≃ 5 MeV (10 MeV). The maximum finite volume effect appears for smaller λ. Even for this maximum case, the expected error due to neglecting such effects is about the same size as the statistical error. For the analysis of chiral extrapolation, we mainly use a larger bin of size 15 MeV, for which the estimated finite volume effect is well below the statistical error.
In Figures 6-8 we compare the lattice results with those of N f = 2 χPT at one-loop.
The plots for each β and strange quark mass are shown in separate panels. The lattice data are renormalized with the renormalization factor for the scalar-density operator calculated separately using the short-distance current correlator [27]. The renormalization scheme is  From the data we can see a clear dependence on the up-down quark mass near λ = 0. In χPT, the quark mass dependence is induced at the one-loop order through the functions ∆(0, M 2 ) and G(0, M 2 , M 2 ) as well as through the counter term including L 6 . Another prominent feature of the low-mode spectrum ρ(λ) is the increase below λ ∼ 20 MeV, which is more pronounced for heavier sea quarks, while the rise almost disappears at the lightest up and down sea quarks available at β = 4.17 (upper panel of Figure 6).
The one-loop χPT prediction of N f = 2 is shown by curves in Figures 6-8  introduced assuming a linear dependence of Σ 1/3 (2 GeV) on m s . The value of Σ 1/3 (2 GeV) mentioned is at the physical strange quark mass. Figures 6-8 demonstrate that the χPT curves also show the increase toward λ = 0 especially for heavier sea quarks and nicely reproduce the lattice data, which show the increase below λ ∼ 15-20 MeV. This is not due to a tuning of parameters. In fact, the extra parameter L 6 appearing at the one-loop order controls only the overall shift of ρ(λ) without influencing its λ dependence. The functional form of the pion-loop contribution, Re∆(0, M 2 vi ) and ReG(0, M 2 vv , M 2 vv ), is responsible for the increase toward λ = 0. On the other hand, the oneloop χPT formula does not explain the slight growth toward larger λ above λ ∼ 20 MeV.
The higher order calculations would be needed to describe this regime.  With the N f = 3 χPT in which kaons and η are also taken as the dynamical degrees of freedom of chiral effective theory, the number of parameters is reduced as we do not need to separately model the strange quark mass dependence. It turned out that a formula including Σ 1/3 (2 GeV), L 6 and a parameter to describe the discretization effect as fit parameters does not fit the data well. (χ 2 /dof is larger than 3.5.) It is probably due to too large strange quark mass to be treated within the χPT framework. In fact, our data for ρ(λ) deviates from the one-loop χPT results above λ ≃ 20 MeV. The physical strange quark mass 90-100 MeV is far beyond this threshold.
We determine the parameters Σ and L 6 through a fit of the lattice data while fixing F = 90 MeV. The fit is done for the value of with δ = 0.015 GeV. Both the lattice data and the χPT formula are integrated in the region [0, δ]. This value of δ corresponds to 2δΣ/F 2 ≃ 250 MeV, which is well below the kaon mass. It corresponds to the lowest three bins in the plots shown in Figures 6-8. In this region, the χPT formula describes the data quite well.
The strange quark mass dependence of ρ(λ) is introduced assuming a linear dependence of ρ(λ) on m s . In the narrow range of the strange quark quark mass adopted in our simulation and with the mild dependence of ρ(λ) on m s , this approximation should describe the data well. Namely, we multiply as an overall factor to ρ(λ) in (15) to interpolate the data to the physical strange quark mass. The parameter c s is to be determined by a fit. Here, M (phys) ηss = 687 MeV is a mass of fictitious ss pseudo-scalar meson estimated using the GMOR relation. Our lattice ensembles contain those of different strange quark masses while other parameters are fixed. The strange quark masses in the simulations are chosen in such a way that they sandwich its physical value. We in effect interpolate between them by (19).
Similarly, the discretization effect is parameterized by a linear function in a 2 , multiplying 1 + c a a 2 as an overall factor with c a a fit parameter.
The fit for the all available data points yield Σ 1/3 (2 GeV) = 270.0(1.3) MeV, L 6 = 0.00016 (6), as well as c s = 0.50(30) GeV −2 , c a = 0.00(15) GeV 2 , with χ 2 /dof = 1.29. As advertised, the discretization effect is invisible within the statistical error. Chiral extrapolation ofρ[0 : δ] is shown in Figure 9 as a function of sea up and down quark mass m ud . Data points do not lie on a single universal curve because the data at different strange quark masses are put in the same plot. In other words, there is a significant strange quark mass dependence, which seems to be well described by an overall shift of the curve. Dependence on the lattice spacing is not very significant from the plot, as the fit also suggests. The curvature due to the one-loop correction is not strong but still visible, and makes the chiral limit slightly lower than a naive linear extrapolation in m ud .
We list the possible sources of systematic errors in the following. First of all, the renormalization constant Z S (2 GeV) determined in [27] contains some errors. (The numbers are given above.) The size is 1.4%, 1.0% and 0.8% for coarse, medium and fine lattices, respectively. We take the largest error, 1.4%, to be conservative, for the estimate of the error for Σ(2 GeV). When we quote the number for Σ 1/3 (2 GeV), we therefore assign 0.5% as an estimated systematic error from this source.
The discretization effect is well under control in our calculation. In fact, our fit implies that the lattice-spacing dependence is consistent with zero. Although it is insignificant, by keeping the term describing this effect in the fit function, we can take account of possible systematic effect. We therefore do not add extra errors from the discretization effects.
Finite volume effect is explicitly checked on the ensembles with the lightest pion (∼ 230 MeV) as shown in Figure 5. We do not observe any statistically significant difference between the two volumes (32 3 and 48 3 ), which is consistent with an expectation from χPT, i.e. the predicted size of the finite volume effect is about 5% for the smaller lattice and is about the same size as the statistical error. For heavier pions the χPT predicts exponentially suppressed finite volume effects. Therefore, for all the data used in the fit to extract the chiral condensate, this source of error is within our statistical error. (Note that the smaller volume data at M π ∼ 230 MeV are not included in the fit.) Higher-order corrections from χPT may be significant especially for larger λ and heavier quarks. Since we can explicitly confirm the consistency of the lattice data with the one-loop χPT for its λ-dependence in the range of our analysis, we expect that two-loop correction is insignificant below λ = 15 MeV. We checked that the result with a slightly smaller bin size, 10 MeV, is consistent within the statistical error. Also for the quark mass, the one-loop χPT fits the lattice data well up to the data points of heaviest pion masses (∼ 500 MeV).
In order to examine the significance of the higher order effects, we tried to fit the data with a function including the analytic terms of O(M 4 /F 4 ). The coefficient obtained from such an analysis is of order of 3 × 10 −6 and statistically consistent with zero. The best fit value of Σ 1/3 (2 GeV) is shifted by only 0.1 MeV, which is much smaller than the statistical error. We can conclude that such effects are well below the statistical error in our analysis.
There is a potential effect of slightly inaccurate implementation of the Ginsparg-Wilson relation with the Möbius domain-wall fermion. As we already discussed, the 4D effective operator of the Möbius domain-wall fermion violates the Ginsparg-Wilson relation by the amount characterized by the residual mass, which is about 1 MeV on our coarsest lattice and an order of magnitude smaller on finer lattices. It means that the eigenvalue of the Dirac operator is distorted by the amount of O(1 MeV) on the lattices at β = 4.17. Since the bin size in the analysis is much larger (= 15 MeV), the error due to this effect is minor.
Moreover, the effect should be negligible on finer lattices, and it is also taken into account by the continuum extrapolation. We therefore do not introduce additional error budget for this effect.
Finally, our input value for lattice spacing has an error of 1.7%, which affects dimensionful quantities, including the chiral condensate. We therefore add this size of error for Σ 1/3 .
Having these various systematic errors considered, we quote where the errors are those from statistical, renormalization, and lattice scale, respectively.
Adding in quadrature, the total error is 4.9 MeV, which is 1.8%. The Flavour Lattice Averaging Group (FLAG) quotes the chiral condensate for N f = 2+1, Σ 1/3 (2 GeV) = 274 (3) MeV [32], as an average of [33][34][35][36]. They are obtained by fitting meson masses and decay constants with the χPT formulae, where the chiral condensate appears as a coefficient in the Gell-Mann-Oakes-Renner (GMOR) relation. Our result (20) is consistent with the world average and the precision is comparable.

Conclusion
The eigenvalue spectrum of the Dirac operator reflects the quantum effects of QCD. The near-zero eigenvalue regime is special, as it can be connected to the order parameter of spontaneous chiral symmetry breaking in QCD, i.e. the chiral condensate. This relation known as the Banks-Casher relation can be extended to the case of finite λ as well as finite quark masses using χPT. This work provides a direct test of these relations by calculating the spectral function in lattice QCD simulations.
The Möbius domain-wall fermion formulation used in this work to define the Dirac operator possesses an approximate chiral symmetry with an error of order 1 MeV at most, and the accumulation of the eigenvalues above this value is not much affected by this artifact.
We extract the chiral condensate from the spectrum below 15 MeV by fitting the lattice data with the χPT formula. The discretization error is well under control and even extrapolated away to the continuum limit using relatively fine-grained lattices of a = 0.080-0.044 fm.
The remaining uncertainty is at the level of 2% for Σ 1/3 (2 GeV). This provides a precise test of the GMOR relation, since there is no free parameter left for the leading-order equation m 2 π /m = 2Σ/F 2 once m π and F are calculated. The agreement of our result 270.0(4.9) MeV with that of an average of previous results obtained through GMOR gives further evidence supporting χPT as an effective theory of QCD at low energies.
The eigenvalue filtering technique utilized in this work is proven to be effective to obtain the spectral function of the Dirac operator. In this analysis we used only the near-zero regime of the eigenvalues, while the entire spectrum is calculated as a by-product. Such information may be useful to extract the mass anomalous dimension of QCD with a non-perturbative method as discussed in [19].