Lattice computation of the Dirac eigenvalue density in the perturbative regime of QCD

The eigenvalue spectrum $\rho(\lambda)$ of the Dirac operator is numerically calculated in lattice QCD with 2+1 flavors of dynamical domain-wall fermions. In the high-energy regime, the discretization effects become significant. We subtract them at the leading order and then take the continuum limit with lattice data at three lattice spacings. Lattice results for the exponent $\partial\ln\rho/\partial\ln\lambda$ are matched to continuum perturbation theory, which is known up to $O(\alpha_s^4)$, to extract the strong coupling constant $\alpha_s$.


I. INTRODUCTION
The Dirac operator D is a fundamental building block of gauge theories such as Quantum Chromodynamics (QCD), the underlying theory of strong interaction. It defines the interaction between quarks and the background gauge field in a gauge invariant manner.
Any observable consisting of quarks in QCD can be written in terms of its eigenmodes, i.e. eigenvalues and their associated eigenfunctions, after an average over background gauge configurations with a weight determined by the path-integral formulation. The eigenvalues are gauge invariant and any quantity built from them is invariant under gauge transformations.
The eigenvalue distribution of the Dirac operator reflects the dynamics of QCD. The best-known example is the Banks-Casher relation [1], which relates the near-zero eigenvalue density to the chiral condensate ψψ , the order parameter of spontaneous chiral symmetry breaking in QCD. In the lattice simulation of QCD, the determination of the chiral condensate can therefore be performed by numerically counting the number of near-zero eigenmodes on finite-volume lattices. In order to correctly identify the near-zero eigenmodes, chiral symmetry needs to be precisely realized in the lattice fermion formulation.
A previous work [2,3] utilized a chirally symmetric lattice operator, the so-called overlap-Dirac operator, and explicitly calculated the individual eigenvalues to obtain the near-zero eigenvalue density. The eigenvalue density can also be estimated stochastically. The first such attempt was performed using Wilson fermions [4]. More recently we applied a slightly different stochastic approach to calculate the eigenvalue spectrum of domain-wall fermions, which are also chiral lattice fermions, and achieved a precise determination of the chiral condensate [5].
Apart from the near-zero eigenvalues, the eigenvalue spectrum ρ(λ) carries more information about the dynamics of the system. In particular, its scale dependence dρ(λ)/dλ is precisely estimated by perturbation theory at sufficiently high energy scales where the strong coupling constant α s is small, i.e. λ ≫ Λ QCD , with Λ QCD the QCD scale. It is therefore of interest to test the perturbative expansion of QCD against non-perturbatively calculated lattice results. So far, tests of high-order perturbation theory against non-perturbative lattice calculations have been performed for the vacuum polarization function at the order of α 4 s [6,7] and for the charmonium correlation function at α 3 s [8,9]. Additionally, stochastic perturbation theory has been applied to simple quantities such as the plaquette [10] and the static quark self-energy [11]. In this work, we perform another such test at α 4 s for the scale dependence of ρ(λ). Since the spectral function can be calculated very precisely on the lattice, it also serves as an input for the determination of α s .
There have been a number of lattice studies of the spectral function aiming at finding the so-called conformal regime of many-flavor QCD or related models, see e.g. [13,14] for instance, and [22] applied the same technique to two-flavor QCD. For the systems that have conformal (or scale) invariance, the scale dependence of the spectral density may be parametrized by a constant γ * , which corresponds to the mass anomalous dimension of associated fermions. For QCD, for which the coupling constant runs and eventually diverges in the low-energy region, this relation does not hold beyond the one-loop level. See the discussions in Section II for details. This paper presents a lattice calculation of the Dirac spectral density in the perturbative regime. We calculate the eigenvalue density in the whole energy range from zero up to the lattice cutoff with the domain-wall fermion formulation using a stochastic technique to evaluate the average number of eigenvalues in small intervals. The lattice results for the eigenvalue density are then extrapolated to the continuum limit using data at three lattice spacings in the range a ≃ 0.044-0.080 fm. We then investigate the consistency with continuum perturbation theory, which is available up to O(α 4 s ). We also attempt to extract α s using the spectral function as an input.
Since we are interested in a relatively high-energy region where the Dirac eigenvalue λ is not much smaller than the lattice cutoff 1/a, discretization effects need to be eliminated as far as possible, in order to achieve precise results. At tree level we calculate the eigenvalue of the Dirac operator constructed from domain-wall fermions and identify discretization effects.
It turns out that discretization effects can be significantly reduced by choosing a value of Pauli-Villars mass different from its standard value. Domain-wall fermions are first defined in five-dimensional space and the eigenmodes on their four-dimensional surface are taken as the physical fermion modes. Our construction corresponds to a slightly different prescription to cancel the bulk effects of the five-dimensional fermion modes leaving the physical modes on the four-dimensional surface. Using this scheme, the remaining discretization effects are made small, so that we are able to extrapolate the lattice data to the continuum limit with a linear ansatz in a 2 .
The rest of this paper is organized as follows. In Section II we briefly discuss the per-turbative results for the Dirac spectral density in QCD. It contains a review of perturbative results for the spectral function ρ(λ) and its exponent. A discussion of the discretization effects for domain-wall fermions is given in Section III. We then describe the lattice setup as well as the method to calculate the eigenvalue density in Section IV. The lattice results are presented and compared with the perturbative expansion in Section V. The extraction of the strong coupling constant is discussed in Section VI. Our conclusions are in Section VII. Some details of the perturbative coefficients are in Appendix A, and a discussion of the locality of the modified Dirac operator we used in this work is found in Appendix B.
A preliminary version of this work was presented in [23].

II. DIRAC EIGENVALUE DENSITY
The eigenvalue density ρ(λ) of the Dirac operator D is defined as where · · · stands for an average over gauge field configurations. The eigenvalue λ k of D depends on the background gauge field. In the free quark limit, λ k may be labeled by the four-momentum p µ as |λ k | 2 = p 2 µ , and the eigenvalue density is given by a surface of a three-dimensional sphere in momentum space: ρ free (λ) = N c (3/4π 2 )|λ| 3 . (N c is the number of colors, N c = 3.) In theories that exibit conformal invariance, there exists a relation ρ FP ∼ |λ| 4/(1+γ * )−1 with γ * the mass anomalous dimension of the theory [12]. It also applies to the case where the theory approaches a renormalization group fixed point. This relation has been utilized to extract the mass anomalous dimension at the fixed point in many-flavor QCD-like theories [13,14].
In QCD, the spectral density has been perturbatively calculated in the MS scheme to order α 3 s [15,16] using where the coefficients ρ (i) depend on λ and are calculated as Here, L λ ≡ ln(λ/µ) and µ denotes the renormalization scale. The numerical constants are  [16]. Numerical results for N f = 2 are also available in [16].
The Dirac eigenvalue λ is renormalized in the same way as the quark mass. It is therefore scale dependent, and the scale is set to µ, i.e. λ(µ). Here and in the following sections, we suppress this scale dependence for brevity. λ denotes its absolute value |λ|. (λ is pure imaginary in the Euclidean continuum theory.) One can reconstruct the scale-dependent coefficient of the next order, i.e. O(α 4 s ), through the renormalization group equation. Namely, we start from which follows from the scale invariance of the mode number, an integral of the spectral density M 0 dλ ρ(λ) with an upper limit M [4]. This can also be understood from the scale invariance of the scalar density operator mqq. Since the spectral function is given by the chiral condensate qq with the valence quark mass set to an imaginary value i|λ|, the renormalization group equation is identical to that forqq. Here the beta function β(α s ) and mass anomalous dimension γ m (α s ) are defined as and known up to O(α 6 s ) and O(α 5 s ), respectively. Their explicit forms are summarized in Appendix A for convenience. Since the µ-dependence of ρ(λ)/λ 3 appears only through the form L λ = ln λ/µ, we may rewrite (II.6) as with K(α s , L λ ) defined by ρ(λ) = (3/4π 2 )λ 3 K(α s , L λ ). By solving this equation at each order in α s , one can determine the scale dependence of ρ(λ), i.e. the terms containing L λ .
Since β(α s ) and γ m (α s ) start at O(α 2 s ) and O(α s ) respectively, the L λ dependent term of ρ(λ) may be determined up to order α n+1 s when a constant term, i.e. the coefficient at L λ = 0, is known at order α n s . The scale-dependent terms thus obtained are summarized in Appendix A.
We also introduce the exponent of the spectral function F (λ), defined by Its perturbative expansion is calculated through order α 4 s using the formula for ρ(λ) summarized in Appendix A: where the coefficients F (k) for N f = 3 are Here, ζ 5 ≃ 1.03692, andρ (3) is the value of ρ (3) at λ = µ, which is given as −ρ (3) = −c 43 /(64π 3 )+c 41 /(256π) using (II.5). F (λ) is defined in the MS scheme at a renormalization scale µ, and the coefficients F (k) have a dependence on L λ = λ/µ. We may also obtain the term at O(α 5 s ) except for the constant term, which comes from an as yet unknown constant ρ (4) defined, analogue withρ (3) , as Let us attempt to estimate the size of the uncertainty in F (λ) due to the unknown constantρ (4) . As we will see below, the value of λ we are going to use to determine α s (µ) is in the range 0.8-1.2 GeV. We therefore take λ = 1.2 GeV as a representative value in the following estimate. By setting µ = λ, the series (II.11) is simplified to where the last term contains the unknown constant d 5 = (−18.0000ρ (4) + 3751.43)/π 5 . Since the coefficients grow as fast as a factor of 2 from one order to the next in (II.17), it is natural to assume that the size of d 5 is about 25. To be conservative, we allow ±50 for the uncertainty in d 5 . Using α s (µ = 1.2 GeV) = 0.406, the size of the O(α 5 s ) term is then ±0.55, which is 18% of the total.
If we choose a higher value for the renormalization scale, e.g. µ = 2.5λ or 5λ, the series becomes F (λ) µ=2.5λ = 3 − 2.54648α s − 6.31432α 2 s − 4.43721α 3 s + 24.5484α 4 s − (d 5 − 123.421)α 5 s + · · · , (II. 18) or respectively. We find that some coefficients are larger than those for µ = λ but not so much as to spoil the convergence. In fact, with α s (3 GeV) = 0.244 and α s (6 GeV) = 0.191, the size of the uncertainty from the unknown O(α 5 s ) term is 1.4% and 0.4% for µ = 2.5λ and 5λ, respectively, if we assume the same ±50 uncertainty for d 5 . Such reduction in uncertainty is of course achieved at the price of potentially large higher order effects due to logarithmic enhancements, i.e. ∼ ln λ/µ, of the coefficients. We roughly estimate the growth of the µ-dependent coefficient to be a factor of 4 from one order to the next for the case of µ = 5λ, and give an estimate ∼ ±700α 6 s for the contribution of the next order. With α s (6 GeV) = 0.204, this amounts to an uncertainty of 2.6% for F (λ). We use this estimate when we quote the value of α s (6 GeV) determined using (II. 19). hand, we find the series converges well when µ is taken to be 3 GeV or higher. In particular, the convergence for λ = 3 GeV or larger is reasonably good for both µ = 2.5λ and 5λ.
The range of convergence extends to lower values of λ, down to ≈ 1 GeV, when µ is set to the higher value, i.e. 6 GeV. This is slightly counter-intuitive because it is generally recommended that the renormalization scale µ is set to be close to the scale of interest, which in this case is µ ≈ λ. The optimal scale does, however, depend on the quantity; for F (λ) it turns out that the higher renormalization scale yields better convergence.
Below λ = 0.8 GeV, the convergence gets worse even with µ = 6 GeV, i.e. the difference between third and higher orders becomes more significant. The main lattice results we use in this analysis are at λ = 1.2 GeV or slightly lower, so the choice of µ = 6 GeV has the advantage of better convergence.

III. DISCRETIZATION EFFECTS WITH LATTICE DOMAIN-WALL FERMION
Since we are interested in the relatively high energy region (∼ 1 GeV or higher) of the Dirac eigenvalue spectrum, discretization effects could be sizable on the lattices with inverse lattice spacing 1/a = 2.4-4.6 GeV. In this section we examine the size of lattice artifacts that affects the Dirac spectrum using the free field limit. The lattice artifacts depend on the details of the lattice formulation. We consider here the domain-wall fermion formulation [17,18], which we used in the numerical simulations.
The domain-wall fermion and its Möbius generalization [19] are formulated on a fivedimensional (5D) Euclidean space, and the mapping to four-dimensional (4D) space is given by in the limit of large fifth dimension L s → ∞. The overlap operator [21] is realized with a polar approximation to the sign function "sgn". Note that the conventional form is recovered for am PV = 1. The Möbius kernel operator D M is given by The general form of the 4D overlap operator (III.2) satisfies the Ginsparg-Wilson relation It can also be shown that the exponential locality of the 4D effective operator is satisfied for any non-zero value of am PV . (See Appendix B for details.) Therefore, domain-wall fermions with any am PV = 0 are a theoretically valid choice. In this work we use the standard choice am PV = 1 in the numerical simulation of sea quarks, and reinterpret the results to the cases of am PV = 1 for valence quarks.
We calculate the eigenvalues of the Hermitian operator a 2 D † ov (m f , m PV )D ov (m f , m PV ) in the massless limit, m f = 0. Leaving the dependence on am PV , we denote its eigenvalue as a 2 λ 2 (am PV ). Since a 2 D † ov D ov 's with different values of am PV commute with each other, their eigenvalues are simply related by (III.5) In the limit am PV → ∞ this corresponds to the projection of the eigenvalue of aD ov in the complex plane to the imaginary axis, which was adopted in [2,3,5].
At tree-level, it is straightforward to calculate the eigenvalue a 2 λ 2 for a plane wave state with a given four-momentum ap µ . To obtain the spectral density, we count the number of states that give an eigenvalue in an interval [a 2 λ 2 , a 2 (λ ± δλ) 2 ] from randomly chosen points of momenta ap µ , each component of which is in the range ap µ = [−π, +π]. We generate a large number of points (up to 10 12 ) so that the statistical error becomes negligible for our choice of bin size aδλ = 0.0025 for the computation of ρ(λ).
The results for domain-wall fermions with various am PV as well as that for Wilson fermion are shown in Figure 2. All formulations coincide with the continuum limit ρ(λ) = 3λ 3 /4π 2 below aλ ≃ 0.3, as expected. Above this value, on the other hand, the spectral density for the domain-wall fermion with Pauli-Villars mass am PV = 1 overshoots the continuum curve very rapidly. This is understood to be because the maximum eigenvalue with am PV = 1 is aλ max (am PV ) = 1, so high eigenvalues tend to rapidly increase in frequency aλ increases.
With am PV > 1, the maximum eigenvalue is stretched by a factor of am PV (when b − c = M 0 = 1) as indicated by (III.5). This allows the spectral densities for am PV = 3 and ∞ follow that of the continuum theory up to aλ ∼ 0.5-0.6. The spectral density for Wilson fermions also resembles the continuum spectrum in the same region.
We now consider the exponent F (λ) of the spectral density (II.10). We approximate the differential ∂/∂λ by a symmetric difference: which is valid up to an error of O(δλ 2 ). For the continuum spectrum ρ(λ) = 3λ 3 /4π 2 , the leading relative error for F (λ) is given by 2δλ 2 /λ 2 . In order that this source of error is below 0.01, δλ/λ < 0.07 has to be satisfied. With our choice of aδλ = 0.005, adopted for the study of F (λ), the reliable range of the approximation is aλ > 0.07. Figure 3 shows F 0 (λ), which is F (λ) evaluated at the tree level, with bin size aδλ = 0.005.
It is compared with the corresponding continuum result, which is a constant, i.e. 3. In the plot, the continuum theory is also calculated with the same finite difference (III.6) (black dots). Some deviation from 3 can be seen near aλ = 0, which indicates the size of the systematic error due to the discretized derivative. Our purpose is to extract the exponent in the high energy region, so we ignore the points below aλ < 0.025, thus this source of error is negligible.
The discretization effects we find for the spectral density ( which reproduce the data well for aλ < 0.6. The smaller coefficients for am PV = 3 confirm the suppression of the discretization effect. In the following analysis we choose am PV = 3. Namely, we convert the eigenvalues from λ(am PV = 1) to λ(am PV = 3) using (III.5). We then use the parametrization (III.7) to further correct the lattice data by multiplying by 3/F 0 (λ) to eliminate the leading discretization effects. Finally we take the continuum limit of the lattice data assuming that the remaining effects are linear in a 2 . This assumption is justified by inspecting the real data.
Strictly speaking, the fermion formulation is different between sea and valence when we choose different values of am PV , and the action of the whole system represents that of a partially quenched theory. However, since the correspondence between the eigenvalues of each formulation is one-to-one even in the interacting case, the continuum limit is guaranteed to be the same. It can be seen explicitely in the relations (III.5), as a 2 λ 2 (am PV ) = a 2 λ 2 (1) + O(a 4 λ 4 (1)). We focus on the region small aλ(1) for which the results from aλ(am PV ) are reliable in the continuum limit.
Another concern may be about the validity of the formulation with am PV = 1. As discussed in Appendix B the localization length for the 4D effective operator (III.1) stays finite for finite am PV . Our choice, am PV = 3, keeps the localization length about the same as that of am PV = 1, and the formulation is valid with respect to the locality property.
In Figure 3 we also plot the result with the standard Wilson fermion (crosses), which shows a smooth and mild continuum limit. Its polynomial approximation gives We use 15 gauge ensembles with different parameters as listed in Table I Table I, where the error given for a −1 is from the statistical error in the evaluation of the Wilson flow time. The error from the input value is taken into account separately.
Details of the ensemble generation are available in [32,33]. The same gauge ensembles have so far been used for a calculation of the η ′ meson mass [34], analysis of short-distance current correlators [6,7,35], a determination of the charm quark mass from charmonium correlators [9], as well as calculations of heavy-light meson decay constants [36], and semileptonic decays [37][38][39]. We use the IroIro++ code set for lattice QCD [40]. The spectral function is also renormalized in the MS scheme using the same renormalization constant Z S (2 GeV). This, however, is implicitly done in our calculation when we calculate the eigenvalue density by dividing the number of eigenvalues in a given bin of λ by the corresponding bin size.
When we compare the results with perturbation theory, we evolve the scale at which we renormalize the eigenvalue from 2 GeV to 6 GeV in order to use the more convergent perturbative series as discussed in Section III.

B. Stochastic calculation of the spectral density
On each ensemble, we calculate the eigenvalue density of the Dirac operator by using a stochastic estimator. We outline the method in this section. The details are available in [5].
Suppose we are able to construct a filtering function h(x) that gives a constant (=1) only in a given range [v, w] and is zero elsewhere. The number of eigenvalues of a hermitian matrix A in that range can then be estimated stochastically as where N is the number of normalized random noise vectors, ξ k . An approximation of the filtering function can be constructed using the Chebyshev polynomial T j (x) through with numerical coefficients γ j , which are calculable as functions of v and w. The approximation is valid in the range −1 ≤ x ≤ 1. In practice, one also introduces a stabilization factor g p j in addition to γ j . The Chebyshev polynomial of a given matrix A can be easily calculated using a recursion relation: We constructed the polynomial up to order p = 8,000-16,000 depending on the ensembles.
Details of the method are found in [5,24]. See also [25] for a related work.
The mode number n[v, w] can then be calculated by An important point is that the computationally intensive part ξ † k T j (A)ξ k is independent of the range [v, w]. Once the inner products ξ † k T j (A)ξ k are calculated for all j below some upper limit p, they can be combined with the known factors g p j γ j to construct the estimate for n[v, w] for arbitrary [v, w]. We use this property to obtain the whole eigenvalue spectrum from λ = 0 to the upper limit at the order of the lattice cutoff.
We apply this technique to a hermitian operator A = 2a 2 D ov (m f , 1) † D ov (m f , 1) − 1, whose range is [−1, +1]. The spectral function ρ(λ) can be obtained from the mode number Here we introduce a bin size δ.
We then obtain the estimate of the spectral function at λ by aδ . (IV.5) The error due to the finite bin size is taken into account when we calculate the derivative of ρ(λ). The systematic effect due to the truncation of the Chebyshev polynomial grow for smaller bin sizes. In our calculation, we confirm that this error is less than the statistical error for our choice of bin size.

V. LATTICE RESULTS FOR THE SPECTRAL FUNCTION AND ITS EXPO-NENT
Following the method outlined in the previous section, we calculate the spectral density ρ(λ) on all ensembles in Table I. We introduce one noise vector per configuration. We present the results in the MS scheme. The renormalization scale µ is set to 2 GeV for the spectral density, and is transformed to 6 GeV when we discuss the matching of the exponent F (λ) with perturbation theory. The bin size δ = 0.05 GeV is chosen in physical units (for the renormalization scale of 2 GeV) such that the error due to the truncation of the Chebyshev polynomial is less than the statistical error for any single bin used in the analysis. 3 and am PV → ∞, the scaling violation is milder than anticipated from the tree-level analysis discussed in Section III.
As described in Sections II and III, we extract the exponent F (λ) of the Dirac eigenvalue spectral density ρ(λ). The derivative in terms of λ is numerically performed using (III.6).
For bin size δλ = 0.05 GeV, the systematic error due to the discretized derivative is about 0.008 (0.3%) in the lowest bin used in the final analysis, λ min = 0.8 GeV. This is smaller than the statistical error of the corresponding data point.
In Figure 5 we plot F (λ) obtained at each lattice spacing. The renormalization scale is set to µ = 6 GeV. The results at three lattice spacings are corrected by a factor 3/F 0 (λ), derived from (III.7), to eliminate the leading discretization effects as discussed in Section III. We find strong discretization effects even after including this correction for the leading discretization effect, especially for the standard Pauli-Villars mass, am PV = 1 (top panel). It is much improved with am PV = 3 (middle panel) and with am PV → ∞ (bottom panel). In the plots, the data points shown by gray symbols have aλ > 0.5, for which the discretization effects seem significant even after the correction. In the plots we also show results from the perturbative expansion calculated at each order, from O(α s ) to O(α 5 s ). The strong coupling constant is taken from the world average Λ (3) QCD = 332(17) MeV for three-flavor QCD [48], and the unknown fifth-order constant is taken to be d 5 = 0 as a representative value. The lattice data are transformed from µ = 2 GeV to 6 GeV using the renormalization group equation with the same value of Λ One can see that the lattice data are in good agreement with the perturbative estimate for am PV = 3 (middle panel) except for those points that are grayed out.
Before comparing the lattice results with perturbation theory in more detail, we extrapolate the lattice data for F (λ) to the continuum limit. We choose the data at am PV = 3 and am PV → ∞ in the following and extrapolate to the continuum limit assuming a linear dependence in a 2 . Figures 6 and 7 show the continuum extrapolation for some representative values of λ. With our lattice parameters, the eigenvalues at λ = 1.22 GeV or below satisfy the condition aλ < 0.5 for all three lattice spacings. For these bins of eigenvalues, the remaining a 2 dependence is well under control as demonstrated in Figure 6. The slope in a 2 is small and the extrapolated results of am PV = 3 and am PV → ∞ agree well with each other. This is consistent with a naive expectation of the size of remaining discretization effect of O(α s a 2 λ 2 ), which is about 5%. The remaining error after the continuum extrapolation is therefore much smaller. For our data, it is estimated to be smaller than the statistical error for each bin of λ.
Above λ ≃ 1.22 GeV, the coarsest lattice fails to satisfy the condition aλ < 0.5. For these bins we ignore the data from the coarsest lattice and extrapolate the other two finer lattice data to the continuum assuming no dependence on a 2 , which is equivalent to taking a weighted average. An example is shown in Figure 7. Here, the results from am PV = 3 and am PV → ∞ agree with each other on the two finer lattices, from which the continuum limit is estimated. The data points at the coarsest lattice spacing, at around a 2 = 0.16 GeV −2 , deviate substantially from the straight lines that represent the average of the finer two lattices.  Again in this plot we overlay the perturbative expansion of F (λ) up to order α 5 s with unknown coefficient d 5 in (II.17) set to zero. Taking α s (6 GeV) = 0.191, which is converted from the world average, we find reasonable agreement between the lattice results and perturbation theory between λ ≃ 0.8 GeV and 2.2 GeV, which suggests that α s (6 GeV) extracted from the lattice data are in fair agreement with the world average.
When we compare the perturbative expansion with the lattice data, we need to consider the non-perturbative effect that may arise as a power correction to the spectral function. According to the general form of the operator product expansion (OPE), the spectral function ρ(λ) receives a correction of the form ∼ ψ ψ at the first non-trivial order. It is suppressed by three powers of λ relative to the leading order contribution ∼ λ 3 . This is parametrically consistent with the Banks-Casher relation ρ(0) = − ψ ψ /π, which is valid in the limit of λ → 0, but it may not be a simple consequence of the OPE because the expansion breaks down for small λ. In any case, the suggested form of the power correction, i.e. a constant in ρ(λ), does not contribute to the exponent F (λ). In other words, the power correction to F (λ) is highly suppressed, and the numerical result from lattice QCD, Figure 8, supports this expectation.

VI. EXTRACTION OF α s
Using the lattice results obtained for the spectral function ρ(λ) as an input, we attempt to extract the strong coupling constant α s . We use the perturbative expansion of F (λ) up to order α 4 s as well as an estimate for the O(α 5 s ) term. The explicit formula is given in (II.11)-(II.16), depending on the value of L λ = λ/µ. We solve the equation for α s (6 GeV) with an input F (λ) on the right-hand side. The error due to the unknown parameter in the O(α 5 s ) term is estimated as discussed below. The determination of α s (6 GeV) may be done for each value (or bin) of λ by solving equation (II.11). The results are obtained as a function of λ, which is plotted in Figure 9. The results for α s (6 GeV) are nearly independent of λ for λ ≃ 0.8 GeV or higher. This is expected from the plot of F (λ), Figure 8, with which we demonstrate that the perturbative estimate follows the lattice data down to λ ≃ 0.8 GeV. Below that, the perturbative expansion is expected to rapidly break down. The statistical error is larger in the high energy region, the statistical signal is not as good as for other lattice spacings. We take the central value from the bins between λ = 0.8 and 1.25 GeV, since the data at all three lattice spacings are included in the continuum extrapolation. The statistically averaged value α s (6 GeV) = 0.204(2) thus obtained is also plotted in Figure 9 on the right panel.
In the following, we describe the possible sources of systematic error and our estimates for them.
The leading discretization effects are removed by the continuum extrapolation, and we estimate the remaining error from the difference between the results with Pauli-Villars masses am PV = 3 and ∞. From an explicit calculation, we estimate the error for α s (6 GeV) as

±0.003.
Discretization errors may also be estimated by varying the number of points included in the continuum extrapolation. For the data points below λ ≃ 1.2 GeV, for which three β values are included in the continuum extrapolation, we attempt to fit the data without the coarsest point and take its difference from the three-point extrapolation as an estimate of the systematic error of this source. This leads to an estimate of ±0.006 for α s (6 GeV). We therefore combine these estimates in quadrature for a conservative estimate of discretization effects as ±0.007.
The perturbative error from unknown higher order coefficients is estimated by allowing ±50 for d 5 of the O(α 5 s ) term (see the discussions in Section II). This amounts to an uncertainty of ±0.012 for F (λ) and gives a variation of ±0.002 for α s (6 GeV). We also estimate the uncertainty due to the unknown O(α 6 s ) contribution as discussed in Section II. This adds another ±0.052 for F (λ) and thus ±0.007 for α s (6 GeV). Overall, the perturbative uncertainty for α s (6 GeV) is estimated to be ±0.007.
The error due to the scale setting enters as an overall shift of the physical scale, which affects α s (6 GeV) due to a change of the QCD scale Λ QCD . The leading scale dependence α s (µ) ≃ (−2β 0 ln(µ 2 /Λ 2 QCD )) −1 implies that the relative uncertainty 1.7% for the input value of t 0 amounts to an uncertainty of ±0.001 for α s (6 GeV). Similarly, the uncertainty in the renormalization scale Z S (2 GeV) gives an overall shift of the eigenvalue, which acts as a shift of the scale. The maximal uncertainty of Z S (2 GeV) is 1.5%, which leads to ±0.001 for the estimated error in α s (6 GeV) from this source.
The error estimates given above are summarized in Table II. By summing them in quadrature, we estimate the total error to be ±0.010, which includes the statistical error. Thus we obtain α s (6 GeV) = 0.204 (10). The total uncertainty is dominated by discretization effects and higher order corrections of the perturbative expansion.
The result may be converted to the value at the Z boson mass scale α (5) s (M Z ) using the four-loop β function. The threshold effect is taken into account when switching to the renormalization scheme of four and then five dynamical quark flavors. We obtain α (5) s (M Z ) = 0.1226 (36), which may be compared to the average of the lattice results 0.1182 (12) [42] (the original works contributing to this average are [43][44][45][46][47]) as well as the world average of the Particle Data Group, 0.1181 (11) [48]. Our result is slightly higher but consistent within about one standard deviation. Compared to our own result, 0.1177(26) from the charmonium current correlator [9], which is obtained on the same set of lattice ensembles, the present result is higher by about 1.1 standard deviations.

VII. CONCLUSIONS
In this work, we perform a lattice calculation of the Dirac spectral density in QCD using the technique to stochastically estimate the number of eigenvalues in given intervals. We use lattice ensembles generated with 2+1 flavors of dynamical quarks. Discretization errors are subtracted as far as possible using an estimate in the free-quark limit. The results are extrapolated to the continuum limit from data taken at three lattice spacings.
The lattice results for the exponent of the spectral density are in good agreement with the perturbative QCD calculation available to order α 4 s . This agreement is highly non-trivial since the value of the exponent is about 30-40% lower than its asymptotic value and the scale (or the eigenvalue λ) dependence is well reproduced. This observation adds another piece of evidence of the validity of QCD in both perturbative and non-perturbative regimes.
We extract the strong coupling constant using the lattice calculation of the spectral function as an input. As in other determinations of α s , the control of discretization effects is essential because we are working at a relatively high energy scale. The results are in agreement with other determinations, although the error is not competitive. The main reason for the relatively large uncertainty is the remaining discretization effects. Namely, the region of λ is limited to 1.2 GeV or below to fully control the continuum limit, and this energy region is not optimal for the perturbative expansion to rapidly converge. A lattice calculation with finer lattice spacing and/or with improved discretization errors would be necessary to improve the precision. In this appendix, we summarize useful formulae for the perturbative calculation of the Dirac spectral density.
We can confirm this for the free field case. We calculate the Fourier transform of the domain-wall fermion operator and obtain its tail at long distances |x|. In general, it shows an exponential fall-off ∝ e −θ|x| . Its rate θ is plotted as a function of am PV in Figure 10.
It turns out that the localization length is optimized in the range am PV ≃ [1,4]. Our choice am PV = 3 is within this range, and it's an equally valid choice as a domain-wall fermion formulation.