Dirac spectral density and mass anomalous dimension in 2+1 flavor QCD

We compute the Dirac spectral density of QCD in a wide range of eigenvalues by using a stochastic method. We use 2+1 flavor lattice ensembles generated with Mobius domain-wall fermion at three lattice spacings ($a=0.083, 0.055, 0.044$ fm) to estimate the continuum limit. The discretization effect can be minimized by a generalization of the valence domain-wall fermion. The spectral density at relatively high eigenvalues can be matched with perturbation theory. We compare the lattice results with the perturbative expansion available to $O(\alpha_s^4)$.


Introduction
Eigenvalue spectrum of the Dirac operator carries rich information of QCD dynamics. A well-known example is the Banks-Casher relation [1], which relates the chiral condensate, a consequence of nonperturbative QCD dynamics, to the low-lying eigenvalue density. Through this relation, the chiral condensate has been calculated on the lattice and precise determination is achieved [2].
In this work, we calculate the Dirac spectrum in the high energy region where the eigenvalue is much larger than the QCD scale Λ QCD . The spectral function is mostly perturbative in this region and the perturbative coefficients are known up to O(α 4 s ). In the past, the Dirac spectral density has been used to extract the mass anomalous dimension of many-flavor QCD, which is expected to exhibit conformal scaling [3][4][5]. In this work we utilise the same quantity but for the case of non-vanishing β function, in order to study the convergence of perturbation theory and to extract the strong coupling constant α s .
By generalizing of the Banks-Casher relation, one can write the Dirac spectral density in terms of the chiral condensate at an imaginary mass [6], To calculate the Dirac eigenvalue denisity on the lattice, we use a stochastic method with a filtering function approximated by the Chebyshev polynomial. The number of eigenvalues n[v, w] of hermitian matrix A in a range [v, w] can be estimated with a filtering function h(x), which is 1 in the range [v, w] and zero otherwise, as where ξ i is a normalized gaussian noise vector and N is its number. We approximate the filtering function h(x) by a polynomial of the form is the Chebyshev polynomial and k i is a coefficient uniquely determined depending on [v, w]. The Chebyshev polynomial may be constructed by a recursion relation T 0 (x) = 1, Once we calculate the inner products ⟨ξ k T i (x)ξ k ⟩, the bin size [v, w] can be chosen afterwards, so that the whole spectrum is obtained by one pass. The details are in [2,8]. Then, the Dirac spectral density can be written as with lattice volume V, bin size δ, and lattice spacing a. Using this relation, we can calculate the Dirac spectral density by identifying A = 2a 2 D † D − 1, which is constructed by the Dirac operator D.
We need to carefully address the discretization effects since the high modes are sensitive to them, and they strongly depend on the fermion formulation. We investigate the discretization effect for the domain-wall fermion formulation and find that the effect on the eigenvalue largely depends on the parameters in the formulation. In particular, by adjusting the value of the Pauli-Villars mass, the discretization error appearing in the Dirac eigenvalue is greatly reduced without sacrifysing the locality and chirality. This generalization of the domain-wall fermion formulation is discussed in Section 2.
Perturbative expansion of ρ(λ) is another important element of this study. We write down the expansion up to O(α 4 s ) using the renormalization group equation for the spectral density as described in Section 3. Finally we compare our lattice calculation with the perturbative calculation in Section 4.

Domain-wall fermion with generalised Pauli-Villars mass
In the lattice calculation of the Dirac spectral density we try to reduce the artifact due to discretized space-time. Discretization effects become more significant when we consider the physics in the perturbative scale, where the bulk of the effect originates from the fermion formulation. We use the Mobius domain-wall fermion in this work.
The domain-wall fermion may be considered as an implementation of the Ginsparg-Wilson relation with a particular approximation of the sign function. It is defined in five dimensional Euclidean space, and the four-dimensional theory is obtained from its surface modes. The details of the Mobius domain-wall fermion may be found in [9,10].
The four-dimensional effective operator D ov for the Mobius domain-wall fermion may be written as with a fermion mass m f , the Pauli-Villars mass m p , and "sgn" an approximated sign function. Here, the Mobius kernel operator D M is defined as in the momentum space when the background gauge field is absent. The domain-wall height M 0 corresponds to a large negative mass in D W , and the kernel parameters b and c are chosen to realize good approximation of the sign function in finite fifth dimension. We use the standard choice b − c = 1, M 0 = 1, and r = 1, same as the Shamir type domain-wall fermion. In this section, except for Figure 2, we take L s → ∞ for simplicity. It allows us to ignore the dependence on b + c since the sign function becomes exact.
The free field case for the massless Dirac eigenvalue may be written as  Figure 1 shows the eigenvalue λ(m p ) as a function of the absolute value of the momentum |ap|. The variation due to different momentum orientation for the same |ap| is shown by the error bar. We find a clear difference between the choices of m p = 1 (left panel), which corresponds to the standard domain-wall fermion, and another possible choice m p = 3 (right panel). It turned out that m p = 3 is much closer to the continuum relation λ(m p ) = |ap|. The difference may be written as which is valid for arbitrary background gauge field, since D † ov D ov 's with different m p commute with each other.
We therefore choose m p = 3 for the study of the Dirac eigenvalue spectrum. Since the ensembles are generated with m p = 1, it corresponds to a partially quenched setup. It is nevertheless harmless because the relation (7) is one-to-one and it can be understood as a slightly modified observable. In the continuum limit, the difference vanishes. We note that the choice of m p has no effect on the pole structure of fermion propagators since m p changes only the denominator of (4) other than the overall scale.
Theoretically, D OV (m f , m p ) with m p 1 has the properties required for the overlap fermion such as the Ginsparg-Wilson relation, and the exponential locality. The Ginsparg-Wilson relation is slightly modified depending on the Pauli-Villars mass, The exponential locality is guaranteed along with the discussion in [11]. The key idea is that the sign function is not strictly local, but has an exponential locality, and |sgn(x)| = 1. Then any polynomial of sgn(x) has an exponentially dumped upper bound.  In Figure 2 we show the non-perturbative results on our generated lattice configurations at three lattice spacings. The results with m p = 1 (left panel) strongly depends on the lattice spacing, whereas the results with m p = 3 (right panel) shows much milder dependence.

Perturbative calculation of spectral density
On the perturbative side, we construct the O(α 4 s ) coefficient in the MS renormalization scheme for the exponent of the Dirac spectral density. In [7], an O(α 3 s ) calculation of the Dirac spectral density is given. At the renormalization scale µ set to µ = λ(µ = λ), it reads for n f = 3. From now on, for simplicity, we sometimes suppress the renormalization scale dependence of λ(µ). Now we focus on the exponent of ρ(λ), i.e. Because an integral of the spectral density ∫ M 0 dλ ρ(λ) with an upper limit M is scale invariant [12], the renormalization equation may be written as with the mass anomalous dimension γ m (α s ) and the beta function β(α s ), defined as Here the beta function β is known up to O(α 6 s ), and the mass anomalous dimension γ m is known to O(α 5 s ) [13,14]. Using (11) we can construct the O(α i+1 s ) exponent of the Dirac spectral density from the O(α i s ) Dirac spectral density since the beta funciton β and mass anomalous dimension γ m start from O(α 2 s ) and from O(α s ), respectively, and both of them are known at least at O(α 4 s ). We thus obtain with the coefficients F (k) for n f = 3 given as Here, ζ 5 = 1.03692, c 3 ≃ 15993.5/(64π 3 ) − 11292.4/(256π), and L λ ≡ ln(λ/µ). At µ = λ(µ = λ) it is numerically written as The leading order is 3 as expected, because ρ(λ) scales as λ D−1 in D dimensions. We note that the expansion coincides with the relation ρ(λ) ∝ λ 4/(1+γ m )−1 up to the O(α s ) level, which is suggested for conformally invariant theories [4,5]. At the O(α s ) level, the mass anomalous dimension is γ m = 2α s /π, then λ 4/(1+γ m )−1 = λ 3−8α s /π at this order, while using (14,15), ρ(λ) is proportional to λ 3−F (1) α s /π with F (1) = 8. Beyond this order, such a simple relation between F(λ) and γ m is lost because of non-zero β function.

Lattice result
We calculate the Dirac spectral density stochastically on the lattice using (1). We set the bin size to 0.05 GeV and the scale to µ = 2 GeV on each lattice ensemble. Our lattice emsembles are generated with 2 + 1 flavor Mobius domain-wall fermion with lattice spacings 1/a = 2.45, 3.61, and 4.50 GeV. Lattice size is chosen such that the physical volume is constant at about 2.6−2.8 fm, i.e. 32 3 ×64, 48 3 × 96, and 64 3 × 128, respectively. Finite volume effects are negligible since our pion mass is 230 − 500 MeV in this calculation. The Dirac spectral density has only tiny sea quark mass dependence, which can be safely ignored. Our calculation is done with the Pauli-Villars mass m p = 1, and the results are transformed to m p = 3 by (7). To match our lattice calculation with the continuum MS scheme, we use the renormalization constant determined through the short-distance vacuum polarization function analysis [15].   Figure 3 shows the exponent of the Dirac spectral density calculated at each lattice spacing compared with the four-loop perturbative calculation. As we already discussed for the spectral density, the Pauli-Villars mass m p = 3 reduces the discretization effect, and enables us to extrapolate to the continuum limit with much smaller uncertainty. After the extrapolation, the exponent of the Dirac spectral density is shown in Figure 4. The error bar represents the sum of the statistical error, the uncertainty from lattice spacings, and from renormalization constant in quadrature. Gray dots without error bars have χ 2 /d.o.f > 2, i.e. the continuum limit in this region are not reliable because of the large lattice artifact. In particular, the results with m p = 1 are difficult to extrapolate for the eigenvalues λ ≥ 1 GeV, because of the large discritization effects. Those with m p = 3 are reliable up to λ ≃ 3 GeV.
We also pay attention to the reliability of the perturbative calculation depending on the scale of the eigenvalue. We estimate the error from truncation of the perturbative expansion by using a difference between 3-loop and 4-loop results. It should be minimal for λ ≃ µ since the perturbative calculation is written with L λ = ln λ/µ. Figure 4 clearly shows that the difference is small around λ = 3 GeV, which is our choice for the renormalization scale µ = 3 GeV, while the larger error is observed around λ ≃ 2 GeV.
In order to see the argeement quantitatively, we choose the bin at λ = 2.92 GeV, as an example, close to µ = 3 GeV. The exponent of the Dirac specral density is obtained as dlnρ(λ)/dlnλ = 2.34 (4)  Eigenvalue density in the high energy region can be calculated perturbatively, while it is also accessible by lattice calculation. We test the perturbation theory by the lattice calculation after m p = 3 an extrapolation to the continuum limit. The precision of ∼ 5% is obtained for the exponent of the Dirac spectral density.
The result may also be used for the determination of the strong coupling constant α s , an analysis of which is underway.