Quark mass anomalous dimension and $\Lambda_{\overline{\textrm{MS}}}$ from the twisted mass Dirac operator spectrum

We investigate whether it is possible to extract the quark mass anomalous dimension and its scale dependence from the spectrum of the twisted mass Dirac operator in Lattice QCD. The answer to this question appears to be positive, provided that one goes to large enough eigenvalues, sufficiently above the non-perturbative regime. The obtained results are compared to continuum perturbation theory. By analyzing possible sources of systematic effects, we find the domain of applicability of the approach, extending from an energy scale of around 1.5 to 4 GeV. The lower limit is dictated by physics (non-perturbative effects at low energies), while the upper bound is set by the ultraviolet cut-off of present-day lattice simulations. The information about the scale dependence of the anomalous dimension allows also to extract the value of the $\Lambda_{\overline{\textrm{MS}}}$-parameter of 2-flavour QCD, yielding the value $303(13)(25)$ MeV, where the first error is statistical and the second one systematic. We use gauge field configuration ensembles generated by the European Twisted Mass Collaboration (ETMC) with 2 flavours of dynamical twisted mass quarks, at 4 lattice spacings in the range between around 0.04 and 0.08 fm.


Introduction
The spectrum of the Dirac operator in Lattice QCD provides very interesting information about several important properties of QCD. In particular, the low modes of the Dirac operator are intimately linked to the chiral condensate, the order parameter of spontaneous chiral symmetry breaking, via the Banks-Casher relation [1]. This relation has recently led Giusti and Lüscher [2] to a new method of extracting the chiral condensate -from the slope of the mode number ν(M ), which counts the number of eigenmodes of the Hermitian Dirac operator D † D below some threshold value M 2 . The mode number can be very efficiently evaluated using a stochastic method, the so-called spectral projector method. Using this method, we have recently calculated the chiral condensate in the continuum limit, with 2 and 2+1+1 flavours of dynamical twisted mass quarks [3,4]. A similar study using N f = 2 clover fermions was reported in Ref. [5]. Another way of linking the low Dirac eigenmodes to the chiral condensate is provided by chiral random matrix theory [6][7][8][9][10][11][12][13][14]. Apart from the link to the chiral condensate, the Dirac spectrum is also related to the topological susceptibility [2,[15][16][17][18] and the quark mass anomalous dimension. This has become an important tool for extracting the latter at the infrared (IR) fixed point in conformal field theories [19][20][21][22][23][24][25][26][27][28][29].
The issue of the scale dependence (or running coupling dependence) of quantities like the anomalous dimensions and renormalization constants is most often addressed on the lattice using the Schrödinger functional, combined with step scaling techniques [30,31] for applications to the computation of the running quark mass (governed by the quark mass anomalous dimension) see e.g. Refs. [32][33][34][35][36][37][38][39][40]. However, step scaling methods can also be used to extract the scale dependence of the renormalized quark mass in the RI-MOM scheme -see e.g. Refs. [41][42][43][44].
In this paper, we want to investigate whether it is possible to extract the implicit scale dependence of the quark mass anomalous dimension from the spectrum of the Dirac operator in 2-flavour Lattice QCD. In general, any comparison of lattice-extracted quantities with perturbation theory (PT) needs the existence of an energy scale window such that this scale µ is: 1. high enough for PT to be applicable, i.e. µ should be much larger than a typical low energy QCD scale of the order of a few hundred MeV (µ ≫ O(Λ QCD )), 2. low enough to avoid large cut-off effects, i.e. µ ≪ Λ lat , where Λ lat is the lattice ultraviolet cut-off (inverse lattice spacing).
The present-day lattice simulations are typically performed with lattice spacings between around 0.05 and 0.15 fm, which corresponds to cut-offs of ca. 1.3 to 4 GeV. This means that the lattice window for establishing contact with perturbation is very narrow and even its very existence is limited to the finer lattice spacings reached nowadays. Our aim is to investigate whether such a window exists for the quark mass anomalous dimension and whether the latter can be accessed with methods of Refs. [19][20][21][22][23][24][25][26][27][28][29]. If the answer is positive, the method can allow for a relatively cheap extraction of the running of the quark mass anomalous dimension. As such, it can at least complement more standard approaches to the computation of this quantity on the lattice, e.g. in the framework of the Schrödinger functional [32][33][34][35][36][37][38][39][40] or the RI-MOM method [41][42][43][44]. In addition, the information about the scale dependence of the renormalized quark mass allows to extract the Λ MS -parameter of the theory, by matching of the continuum extrapolated lattice data to 4-loop PT. The paper is organized as follows. In section 2, we describe the theoretical principles of the employed method and discuss its potential limitations. Section 3 presents our analysis strategies and section 4 the lattice setup. In section 5, we show our results. Section 6 concludes. Additional tests are presented in 4 appendices.
2 Theoretical principles 2.1 Spectral density and mode number of the Dirac operator The main aim of this paper is to check if it is possible to extract the running of the quark mass anomalous dimension from the lattice QCD Dirac operator spectrum. The formula that will be used in the numerical part can not be derived from first principles. Therefore, we will present here the arguments that lead to it, coming from different premises: perturbation theory and renormalization group.
Scaling of the spectral density in theories with an infrared fixed point The scaling of the spectral density of the Dirac operator ρ(λ) (where its eigenvalues are denoted by λ) is related to the quark mass anomalous dimension. Let us start with a short résumé of the situation in gauge theories with an IR fixed point. Using the properties of these systems, it is possible to show that the scaling of the spectral density of the Dirac operator is related to the scheme-independent mass anomalous dimension at the fixed point [19,20,25]. This relation can be written in the following form: [26] at leading order, where:ρ 0 -dimensionless constant, µ -renormalization scale, γ * m -quark mass anomalous dimension at the IR fixed point. This can be rewritten using the integrated spectral density, i.e. the (dimensionless) mode number: where ν R (M R ) is the renormalized number of eigenmodes of the Hermitian Dirac operator D † D below some renormalized threshold eigenvalue M 2 R , V is the volume andρ is a constant with the dimension of volume. It was shown in Ref. [2] that the mode number is renormalization group invariant, i.e. ν R (M R ) = ν(M ).
Since the mode number is a quantity easily accessible on the lattice, it is, in principle, possible to use the above formula to extract the quark mass anomalous dimension at the IR fixed point. However, in actual lattice simulations, scale invariance is broken by a non-zero quark mass m 1 , leading to the development of a fermion condensate and a mass gap [20]. The lattice results can then be described within mass-deformed conformal gauge theory. In the presence of a non-zero quark mass, the mode number equation (2.2) is modified to: [26] ν R (M R ) = 2V where m R denotes the renormalized quark mass and λ IR is some infrared scale below which quark mass effects are relevant. The mass deformation introduces an extra term ν 0 (m R ), dependent solely on the quark mass. Its practical consequence is that the anomalous dimension at the fixed point γ * m can not be extracted for too small values of the scale M R , where effects of the mass deformation can be large. However, one can get control over these effects by including the term ν 0 (m R ) in fits and simulating at more than one light quark mass, to explicitly check its influence.

Scaling of the spectral density in chirally broken theories
The situation is somewhat similar in chirally broken systems, like QCD. The infrared behaviour of such systems is very different from theories with an IR fixed point. However, the effects of spontaneous chiral symmetry breaking should be relevant only below some scale Λ χ . Hence, one can conjecture that well above this scale, the scaling of the mode number can resemble the one in conformal systems, making Eq. (2.3) valid. However, since no infrared fixed point occurs in chirally broken systems, the quark mass anomalous dimension at the IR fixed point γ * m is replaced by an anomalous dimension γ m (M R ) dependent on the running coupling (and hence implicitly scale-dependent) and Eq. (2.3) becomes: Such an approach was adopted in Ref. [23], where the scale-dependent mass anomalous dimension 2 was extracted for an SU (3) theory with 4 flavours of quarks. Let us recall here the arguments relating the scaling of the mode number to the quark mass anomalous dimension [23], paying special attention to the differences between IRconformal and chirally broken systems.

Spectral density in perturbation theory
It can be shown [24] in one-loop PT that in systems with asymptotic freedom, the following relation holds: where C is a normalization constant, g R is the renormalized coupling and the schemeindependent one-loop quark mass anomalous dimension γ m (g 2 R ) depends on the gauge group and the fermion representation. Thus, the above relation holds, in principle, both in chirallybroken and IR-conformal systems that are asymptotically free.
Note that, in principle, additional contributions to the spectral density can be present in general, e.g. terms growing with lower powers of M . However, such terms would not be visible in PT, i.e. they would modify the above one-loop expression by adding a term like The absence of such term at one-loop implies that the corresponding term in the mode number can only appear non-perturbatively. However, at high enough M , the importance of such terms will become negligible. Hence, it will not affect the hypothesis for numerical evaluation that at high enough M , the method should allow for the extraction of the mass anomalous dimension. It should also be mentioned that terms of such kind could show up also in the extraction of the chiral condensate from the mode number [2,3,5]. It was shown in the analysis of Ref. [3] that such potential effects are numerically small, if M is small enough. An analogous thing has to hold on the other end of the scale of M -for large enough M , any lower order terms (with lower powers of M ) have to be unimportant. Hence, in the following, we will not consider such terms.

Renormalization group scaling of the quark mass and Dirac operator eigenvalues
Let us now consider the renormalization group scaling of the eigenvalues of the Dirac operator and of the quark mass. We assume that the quark mass is multiplicatively renormalizable: where g 0 is the bare coupling and the mass renormalization constant Z m (g 0 , µ) is massindependent (i.e. defined in a mass-independent renormalization scheme). Let us start with the definition of the quark mass anomalous dimension: where g R (µ) is the renormalized coupling and g 0 on the right-hand side is such that it corresponds to the renormalized coupling g R (µ) on the left-hand side (this will be implied in the following formulae). For two chosen scales µ 1 and µ 2 , this equation can be rewritten as: Close to a renormalization group (RG) fixed point, γ m (g R (µ)) depends very mildly on µ, i.e. g R (µ 1 ) ≈ g R (µ 2 ), and the above equation becomes: . (2.9) One can now consider an RG transformation with a scale factor b, i.e. the scale transforms as: µ → µ/b. The bare quark mass scales as: m → m/b. In the renormalized quark mass, there is an additional effect, coming from the running of the renormalization constant with the scale. Thus: Eq.(8) (integration) Eq.(9) (RG transf.) Figure 1: The ratio of MS renormalized quark masses (N f = 2 QCD) at scales µ 1 and µ 2 such that µ 2 = µ 1 + 500 MeV. The blue solid line is the result from Eq. (2.8) (exact integration using 4-loop γ m (g R (µ))) and the black dashed line from Eq. (2.9) (RG transformation valid close to the (ultraviolet) RG fixed point), using 4-loop γ m (g R ( µ 1 +µ 2 2 )).
be: M → M/b 1+γ M . Since the threshold M and the twisted quark mass m (we anticipate the use of twisted mass fermions in the following) renormalize with the same renormalization constant [2], they have to run with the same anomalous dimension, i.e. γ M (g R ) ≡ γ m (g R ).
In the following, we will use the notation γ M (M ) ≡ γ M (g R (M )), since M plays the role of the renormalization scale (see below). Having shown the above scaling properties in the vicinity of an RG fixed point, we can now relate γ M (M ) to the scaling of the spectrum of D † D. In the free theory ν(M ) ∝ V M 4 , while interactions modify this scaling behaviour to ν(M ) ∝ V M α , where α is close to 4 if the eigenvalues are large, i.e. correspond to the ultraviolet. The renormalization group invariance of the mode number, proved in Ref. [2], implies: V M α = b 4 V M/b 1+γm(M ) α and hence α = 4/ (1 + γ m (M )). Thus, ν(M ) ∝ M 4/(1+γm(M )) and the scale-dependent anomalous dimension γ m (M ) parametrizes deviations of the mode number scaling from the free-field theory value of 4, in accordance with the PT formula (2.5). In this way, we have established the relation of the anomalous dimension γ m (M ) to the scaling of the mode number. At this point, it is important to emphasize that this relation is not universally valid and is subject to several conditions restricting its range of applicability. It is therefore essential to discuss these limitations.

Limitations of applicability of the approach
First of all, let us reconsider Eq. (2.8). In general, i.e. if the implicit scale dependence of γ m (g R ) is non-negligible, Eq. (2.8) can not be written as (2.9) -it becomes: with the running of the coupling given by the β-function. In this general case, Eq. (2.10) is modified such that the relation between the scaling exponent α of the mode number and the mass anomalous dimension γ m receives corrections that make the functional form of this relation more complicated (it then includes the integral of Eq.  Fig. 1, which shows the ratio of renormalized quark masses m R (µ 2 )/m R (µ 1 ), in the MS scheme at 4-loops, in QCD with 2 flavours of dynamical quarks, for scales satisfying µ 2 = µ 1 + 500 MeV. The interval of 500 MeV was chosen to match our intervals for extracting γ m from the lattice. The blue solid line shows the result from the application of Eq. (2.8) (or, equivalently, Eq. (2.11)), which takes into account the g R -dependence of γ m and hence its implicit µ-dependence. In this sense, this result is exact, in contrast to the result from Eq. (2.9) (black dashed line in Fig. 1), which is, in principle, valid only in the vicinity of a RG fixed point (ultraviolet in this case). To compute the ratio m R (µ 2 )/m R (µ 1 ) according to Eq. (2.9), we have inserted γ m evaluated in the middle of the interval from µ 1 to µ 2 , i.e. γ m (g R ( µ 1 +µ 2 2 )). Inspection of Fig. 1 shows that the difference between both considered results is minor if µ 1 is larger than approx. 1 GeV. Numerically, this difference amounts to 0.001% at 4 GeV (i.e. for µ 1 = 4 GeV and µ 2 = 4.5 GeV), 0.016% at 2 GeV, 0.05% at 1.5 GeV and 0.38% at 1 GeV. Only below µ 1 = 1 GeV, the effects start to be larger than our statistical errors. However, at such low energies, non-perturbative effects are beginning to be essential and comparison to PT does not make much sense. The chosen difference µ 2 − µ 1 = 500 MeV is a good compromise between availability of lattice data and the discussed effect, which is below 0.5% for scales above 1 GeV. However, even choosing µ 2 − µ 1 = 2 GeV would lead to effects below 0.5% for scales above approx. 1.7 GeV. This results from the following numerical observation -the contribution to the integral in Eq. (2.8) of the intervals [µ 1 , (µ 1 +µ 2 )/2] and [(µ 1 +µ 2 )/2, µ 2 ] balances out as if g R was constant and equal to the value in the middle of the interval and could hence be taken out of the integral. In this way, the use of Eq. (2.9) to relate the quark mass anomalous dimension to the mode number, instead of the exact Eq. (2.8) does not constitute a problem at relevant scales. Note that this is to some extent surprising, since naively one could expect that Eq. (2.9) is valid only in the immediate vicinity of the fixed point, e.g. at scales above 5 or 10 GeV, where no contact to the lattice would be possible at presently simulated lattice spacings.
The second limitation that we would like to shortly discuss is related to the following issue. Given the equation for the massive Dirac operator D † m D m = D † D + m 2 (where D is the massless operator), it is plausible to expect that, at least in some energy range, the eigenvalues of the Hermitian Dirac operator scale in the same way as the quark mass m or the threshold mass M . However, one has to keep in mind that such scaling is not universally valid for all individual eigenvalues. For instance, the behaviour of eigenvalues near the origin gives rise to the chiral condensate in chirally broken systems according to the Banks-Casher relation. Moreover, some eigenvalues are unphysical, e.g. related to the doubler modes (that decouple only strictly in the continuum limit), and the concept of renormalization is ill-defined for them.

Importance of a numerical check of the approach
Summarizing, the domain of applicability of the argument relating the dependence of the mode number on the threshold mass M and the anomalous dimension γ M has to be checked numerically and is not expected to be universally valid, e.g. the behaviour of ν(M ) for very small M gives information about the chiral condensate [2] and not about the anomalous dimension γ m . A priori, therefore, one expects the results obtained from this analysis to be valid for rather large values of M , i.e. such that the properties of the system are governed by perturbative effects. Hence, it is a priori unclear whether it is possible to address these issues with a lattice calculation with presently available lattice spacings (with inverse cutoffs of the order of a few GeV).
The aim of the present paper is to investigate in practice the above mentioned problems, by looking at the Dirac operator spectrum in lattice QCD with N f = 2 flavours of dynamical twisted mass quarks. In particular, we want to check if it is possible to recover values of the mass anomalous dimension predicted by PT for some range of energies where contact can be established between the latter and the lattice theory. This requires either an appropriate renormalization of the threshold parameter M , which will be the subject of the next subsection, or matching to PT.

Renormalization and quark mass anomalous dimension in perturbation theory
The mass scale (threshold for eigenvalues of D † D) M is a bare scale and has to be renormalized in order to make contact with PT. Let us start with an analogy with the computation of the renormalized chiral condensate from the Dirac operator spectrum [2,3]. The bare chiral condensate Σ is extracted from the slope of the mode number vs. M dependence. Since the product M Σ is renormalization group invariant, to calculate Σ MS,µ , the renormalized condensate in the MS scheme, at some scale µ, the threshold M has to be renormalized with the renormalization constant Z P , in the MS scheme, at the scale µ, i.e. M is renormalized which boils down to an extraction of the slope of the renormalized mode number vs. renormalized mass scale M R .
In the present case, we need to renormalize the threshold scale M in an analogous way: i.e. each value of M has to be renormalized with a separate value of Z −1 P , computed at the scale µ = M R and expressed in the MS scheme, since we want comparisons with mass anomalous dimension defined in this scheme and M Σ = Z −1 This renormalization condition will become a basis for our analysis strategy 1.
We remind here the expression for the quark mass anomalous dimension γ m in perturbation theory [45,46]: where a s (µ) = α s (µ)/π = g R (µ) 2 /(4π 2 ) and γ i are known coefficients for i = 0, . . . , 3. In numerical form, γ m for N f flavours reads in the MS scheme: up to 4 loops 3 . The running of the coupling depends on the Λ-pearameter of the theory. Hence, the continuum extrapolated lattice data for the scale dependence of γ m , or equivalently m R , allow to extract the Λ-parameter, denoted by Λ MS , where the superscript stands for 2 flavours and the subscript for the renormalization scheme.

Spectral projector method of computing the mode number
The mode number ν(M ), i.e. the number of eigenvectors of the massive Hermitian Dirac operator (with quark mass m) D † m D m with eigenvalue magnitude below the threshold value of M 2 , can be computed essentially with two methods. First, one can explicitly compute some given amount of n eigenvectors for each gauge field configuration. This gives the mode number ν(M ) for each value of M below the average eigenvalue corresponding to the n-th eigenvector. However, for our present goals, this method is too expensive in terms of computing time, since we want to reach threshold values of M corresponding to mode numbers of O (10 5 ). This makes the other available method, the method of spectral projectors described in Ref. [2], very attractive. Its main advantage is its good scaling with the volume -the cost of the computation scales with V , instead of V 2 as in the case of an explicit computation. Note that the usage of the spectral projector method (i.e. computation of the mode number only for ≈ 20 selected values of M spanning a range of a few GeV) is allowed for present purposes because the relation between the scaling exponent of the mode number and γ m can be derived using Eq. (2.9) and not explicitly using the integral in Eq. (2.8) -else one would need the full continuous dependence of the mode number on the threshold eigenvalue M (as discussed in Sec. 2.1). In this way, one could reach only mode numbers of at most O(10 3 ) (and even this with a higher computational effort than using the spectral projector method) and hence values of M much smaller than the inverse lattice spacing.
For the spectral projector evaluation of the mode number, we use the implementation in the tmLQCD code [48]. This implementation was extensively tested and used in our chiral condensate computation [3].
Here, we shortly describe the method of spectral projectors. We refer to the original work of Ref. [2] for a more complete account. Let us define the orthogonal projector È M to the subspace of fermion fields spanned by the lowest lying eigenmodes of the massive Hermitian Dirac operator D † m D m , with eigenvalues below some threshold value M 2 . The mode number ν(M ) can be represented stochastically by: where the function: is an approximation to the step function θ(−x) in the range −1 ≤ x ≤ 1 and P (y) is in our case the Chebyshev polynomial of some adjustable degree n Chebyshev that minimizes the deviation: for some ǫ > 0. Computing the approximation to the spectral projector È M requires solving the following equation an appropriate number of times: for a given source field η. The parameter M * is related to the spectral threshold value M and the ratio of M/M * depends on the details of the approximation to the projector. For our choice, M/M * ≈ 0.96334 (as shown in Ref. [2]).

Analysis strategy
3.1 Strategy 1 -using Z P as an input The basic relation (2.4) can be used to extract the scale dependence of the quark mass anomalous dimension. We insert the renormalization condition (2.13) and rewrite this relation as: (3.1) Since we will work far away from the scale of spontaneous chiral symmetry breaking Λ χ , we will ignore the term ν 0 (m R ), having checked that it indeed does not influence the extracted values of the mass anomalous dimension if M R 1.5 GeV (see below). The contribution of this term can also be estimated from the spectral density at the origin ρ(0), which gives the value of the chiral condensate according to the Banks-Casher relation [1]. In the low-energy regime ρ(λ) is approximately constant and equal to its value at the origin (Σ/π from the Banks-Casher relation)ν(M ) grows linearly with M and its slope determines the chiral condensate and hence it gives a rough estimate of the contribution of the term ν 0 (m) to the total mode number (the dependence on the quark mass m is implicit in the value of the chiral condensate): The interval where the effects of spontaneous chiral symmetry breaking are important is assumed to lie between the origin and some scale denoted by Λ χ . In the numerical part of this work, we will consider several values of Λ χ to check the robustness of the results with respect to effects of spontaneous chiral symmetry breaking, by estimating the contribution of the term ν 0 (m) to the total mode number (the value of the chiral condensate for this estimate can be taken from the low part of ν(M ) vs. M dependence). Note that the above equation is just the leading-order chiral perturbation theory expression for the effective chiral condensate -see e.g. Eq. (4.2) in Ref. [2], with the identification Λ = Λ χ . Considering only a range of M such that the term ν 0 (M ) can be safely neglected, one can rewrite Eq. (3.1) as: where we symbolically write that the extracted values of γ m (M R ) depend on the lattice spacing a. We will use the twisted mass Dirac operator (see next section) to extract the mass anomalous dimension. The twisted mass Dirac operator gives automatic O(a)-improvement in physical (R 5 -parity even) quantities [49]. In particular, the mode number of the Dirac operator is automatically O(a)-improved [4]. However, this does not necessarily imply the improvement of γ m . Inspecting Eq. (2.4), it can not be excluded that the mode number is O(a)-improved even if the factors that enter it have O(a) effects. Hence, the O(a)improvement of γ m can not be concluded from the improvement of the mode number. For this reason, the continuum limit extrapolations in the numerical part will be performed under the assumption that O(a) effects can be present. With a rather good precision of the method, it can turn out a posteriori that the coefficient of the O(a)-term in the continuum limit extrapolations of γ m is compatible with zero. Actually, this will not be the case -see Sec. 5.4 and Appendix D.
We emphasize that only the continuum limit extrapolated values γ m (M R ) = lim a→0 γ m (M R , a) can be compared to γ m (a s (M R )), provided that the lattice window exists for the contact of lattice simulations with continuum PT (which can be written schematically as: A further check of the method can be performed by rewriting equation This equation implies that the renormalized mode number scales with the fourth power of the renormalized threshold parameter M R for all values of the latter. In practice, cut-off effects can lead to deviations from the above statement. Therefore, we write Eq. (3.5) as: with an "anomalous dimension" Γ m (M R ) which is purely a lattice artefact, i.e. its continuum limit should be zero. Hence, we will call it the "artefact" anomalous dimension.

Strategy 2 -matching to perturbation theory
An alternative method of analysis was proposed in Ref. [23]. It does not require the knowledge of the renormalization constant Z P and instead matching to PT is performed. It consists in selecting one value of β as the reference value (denoted by β ref ) and rescaling lattice eigenvalues for other bare couplings to express them in terms of a uniform scale a ref , i.e. the lattice spacing corresponding to β ref . The rescaling of eigenvalue M β at a given value of β is as follows: . (3.8) In this way, we set the scale, i.e. we know that M ref,matching corresponds to µ matching , which is known in physical units. A non-trivial test of the approach is provided by comparing the scale dependence of the lattice extracted anomalous dimension with PT prediction. Since the latter is scheme-dependent and our procedure implicitly defines a renormalization scheme, we can only compare to the universal one-loop PT expression.

Lattice setup
Our lattice setup consists of the tree-level Symanzik improved gauge action [50] and the Wilson twisted mass fermion action [49,[51][52][53]. The former reads: g 0 is the bare coupling, P 1×1 , P 1×2 are the plaquette and rectangular Wilson loops, respectively. The Wilson twisted mass fermion action is given in the so-called twisted basis by: where m 0 (m) is the bare untwisted (twisted) quark mass. The renormalized light quark mass is given by m R = Z −1 P m. The matrix τ 3 acts in flavour space and χ = (u, d) T is a two-component vector in flavour space, related to the one in the physical basis by a chiral rotation. The standard massless Wilson-Dirac operator D W is: where ∇ µ and ∇ * µ are the forward and backward covariant derivatives. One of the main advantages of the twisted mass formulation is that it allows for an automatic O(a) improvement of physical observables, provided the hopping parameter κ = (8+2am 0 ) −1 , is tuned to maximal twist by setting it to its critical value, at which the PCAC quark mass vanishes [51,[54][55][56][57]. However, the spectrum of the Dirac operator itself is not improved -hence we expect the extracted values of the quark mass anomalous dimension to be contaminated by O(a) discretization effects.
Gauge field configurations that we have used for this work were generated by the European Twisted Mass Collaboration (ETMC) with N f = 2 dynamical flavours of quarks [58][59][60]. The details of lattice parameters considered for this work are shown in Tab. 1. The linear extents of our lattices are relatively small, with L ≈ 1.3 fm. However, we checked the size of finite volume effects by including a larger physical volume for β = 3.9, am = 0.004, with L/a = 24, i.e. a physical volume of around 1.9 fm.
To make comparisons to PT employing Strategy  We give values of Z P at two scalesµ 1 = 2 GeV (given in Refs. [42,63,64]) and another scale µ 2 for which no perturbative running to 2 GeV has been performed -see text for more details.
The latter deserves a longer comment, as it regards the applicability of the method to theories with spontaneous chiral symmetry breaking. In such theories, the spectral density of the Dirac operator does not go to zero near the origin, but instead tends to a constant ρ(0) and leads to a non-vanishing value of the chiral condensate. Such constant value of the condensate produces a constant contribution to the mode number, denoted by ν 0 (m) and independent of M . We have tried fits including different values of ν 0 (m). In particular, one can follow the discussion in Sec. 3 and use Eq. . However, these apparent differences becomes much smaller after extrapolation to the continuum limit. This will be shortly discussed again in Sec. 5.4 and in Appendix A. For now, we anticipate the conclusion that the results in the continuum limit always agree for different values of ν 0 (m) ≤ 40 if M 600 MeV, even for the highest assumed value of Λ χ = 500 MeV. We also remark that although the term ν 0 (m) depends on the quark mass, in practice this dependence is not relevant from the point of view of this analysis, as it varies the value of the condensate (mass-dependent condensate defined in Ref. [2]) by at most 5% when going from quark mass am = 0.004 towards the chiral limit, thus having negligible influence on the value of ν 0 (m).

Finite volume effects
We have investigated finite volume effects by comparing the anomalous dimensions γ m (M ) calculated from ensembles B40.16 and B40.24, i.e. β = 3.9, am = 0.004 and L/a = 16 or 24, which corresponds to L ≈ 1.3 or 1.9 fm. It was found in Ref. [3] that finite size effects in the mode number density ν/V are small when one reaches a linear lattice extent of ca. 2 fm. This is true if the threshold parameter M 50 MeV, i.e. in the range used for chiral condensate extraction. However, if M is increased, finite size effects tend to decrease. As argued in Ref. [2], the difference between finite and infinite volume results for the chiral condensate and hence also for the mode number density is of O(exp(−M Λ L/2)), with M 2 Λ = 2ΛΣ/F 2 , Λ = M 2 − µ 2 , F is the pion decay constant in the chiral limit. The chiral condensate is extracted at M 50 MeV, while here we work with values of M a factor of 5-50 larger. We can thus expect that we observe significant finite size effects only at the lower end of considered values.
The results of comparison of B40.16 and B40.24 are shown in Fig. 4. We observe that the extracted values of γ m coincide for all values of M , ranging from around 500 to 2500 MeV. However, there is a systematic tendency towards discrepancy between the results from both ensembles for small values of M , in accordance with theoretical expectations. Finite volume effects above M 500 MeV are small and the extracted values of the anomalous dimension above this threshold can be considered to be infinite-volume results, thus justifying the use of ensembles with L ≈ 1.3 fm.

Non-perturbative running of Z P
In order to compare the extracted values of the quark mass anomalous dimension to the predictions of PT, we have to take the continuum limit. However, the values of γ m at non-zero lattice spacing can be used to perform non-perturbative evolution of Z P (since Z P runs with the same anomalous dimension as the quark mass).
The comparison of non-perturbative and perturbative evolution of Z P is shown in Fig. 5. The former is extracted at the scale 1/a in the RI-MOM scheme or at some chosen scale 1/X 0 in the X-space scheme and converted to the MS scheme at this scale. Then, the lattice extracted quark mass anomalous dimension is used for the non-perturbative evolution of Z P . Note that the evolution procedure is self-consistent, i.e. γ m is extracted at some bare scale M and this scale is renormalized using the values of Z P self-consistently evolved with this γ m by numerically integrating the defining equation of γ m : d ln Z P (µ) = 2γ m (µ)d ln µ, (5.1) which leads to the relation: ln Z P (µ 2 ) = ln Z P (µ 1 ) + 2γ m (M = Z P (µ 1 )µ 1 ) ln which is used for small differences in µ 1 and µ 2 such that γ m can be considered equal for µ 1 and µ 2 . The argument in the parentheses of γ m (M = Z P (µ 1 )µ 1 ) ensures that such evolution procedure is self-consistent. In general, the non-perturbative running of Z P is "faster", since lattice-extracted quark mass anomalous dimensions are always larger at finite lattice spacing than their perturbative values (see Fig. 7). Since the difference with respect to the continuum values is an O(a) effect, the non-perturbatively evolved Z P is contaminated with additional O(a) effects from the extraction of the mass anomalous dimension from the mode number. However, we emphasize that the fact that we have the non-perturbative running of Z P will allow to obtain in the end the scale dependence of γ m in the continuum and hence we will be able to compare these values with the prediction of continuum PT.

Continuum limit -Strategy 1
We have repeated the procedure described in Sec. 5.1 to extract the anomalous dimension for four ensembles of ETMC gauge field configurations with N f = 2 dynamical flavours of MS = 315 MeV [62]. The starting point of the evolution of Z P is the scale 1/a (β = 3.9, 4.05, 4.2) or 1/X 0 ≈ 1.78 GeV (β = 4.35), i.e. at this scale Z P is the same for the non-perturbative and 4-loop perturbative evolution. Z P for 1-to 4-loop perturbative evolution is chosen equal at 2 GeV. quarks. All of them correspond to a fixed physical situation of L ≈ 1.3 fm and a pion mass of ca. 330 MeV. As we have argued above, our results for γ m can still be considered to be infinite-volume ones (if M 500 MeV) and pertaining to the chiral limit. Our aim is now to relate them to continuum PT. To achieve this, we have to take the continuum limit of lattice results. Fig. 6 shows the values of γ m (M ) vs. bare threshold parameter M for all 4 lattice spacings. We observe a clear dependence of the results on the lattice spacing. In order to take the continuum limit, we have to renormalize the scale M according to Eq. (2.13). In addition, to obtain the values of γ m (M R ) at arbitrary scales M R (and not only the discrete set related to the values of M chosen for the computation of the mode number), we perform quadratic interpolation. The outcome of renormalization of M and interpolation is shown in Fig. 7. The plot also shows the quark mass anomalous dimension γ m (a s (µ)) (since we β=3.9, L/a=16, am=0.004 β=4.05, L/a=20, am=0.003 β=4.2, L/a=24, am=0.002 β=4.35, L/a=32, am=0.00175 Figure 6: Results for all 4 lattice spacings. The spectral threshold M is unrenormalized.   The continuum limit extrapolations are shown in Fig. 8. In all cases, the scaling is consistent with O(a) cut-off effects. The results in the continuum limit are also shown in Tab. 3, for selected values of M R and in Fig. 9 for the whole interval of M R between 1 and 4.5 GeV. In Tab. 3, we also give values of the quark mass anomalous dimension in perturbation theory. We observe very good agreement between the lattice results (continuum-extrapolated) and perturbative values in the range between 1 and around 3.5-4 GeV. As we discussed above, the lower limit of M that yields a reliable result from the point of view of finite volume effects and the importance of the term ν 0 (m) corresponds to M ≈ 500-600 MeV, i.e. M R ≈ 1.3-1.5 GeV. In addition, the difference between 3-and 4-loop perturbative values below 1.5 GeV suggests that the good agreement in the interval between 1 and 1.5 GeV should not be taken very seriously. The upper limit of M R that can be used to simulate on the lattice is related to the lattice cut-off. In our case, the inverse lattice spacings correspond to cut-offs of 2 to 4 GeV -hence, above these values, one expects enhanced discretization effects. This explains the observation that above ca. 4 GeV the lattice results yield γ m  β=3.9, L/a=16, am=0.004 β=4.05, L/a=20, am=0.003 β=4.2, L/a=24, am=0.002 β=4.35, L/a=32, am=0.00175 1-loop 2-loop 3-loop 4-loop cont.limit (stat.+syst.error) cont.limit (stat.error) Figure 9: Comparison of the quark mass anomalous dimension extracted from the lattice to perturbation theory (for Strategy 1). For the continuum limit results, we show the statistical error and the total one, i.e. the combined statistical error, the error originating from lattice spacing value in physical units and the error coming from Z P (see Tab. 1 and Sec. 2.2 for more details).
above the one of PT. The agreement can be regained by an inclusion of higher-order O(a 2 ) cut-off effects in the fits, however at a price of a significantly increased error of the fits.
As a check of robustness of our results, we have performed three further checks, described in appendices.
First, we checked the influence of including the term ν 0 (m) in the fits on our continuum limit extrapolations -see Appendix A. Here we only state the conclusion -this influence is noticeable only for relatively small values of M R . Depending on the value of ν 0 (m), we observe agreement between lattice and PT starting at M R between ca. 1100 and 1500 MeV for ν 0 (m) in the range between 10 and 40.
Second, we followed the approach explained in Sec. 3 (fits of Eq. (3.6)) to extract the "artefact" anomalous dimension Γ m (M R ) -see Appendix B. We observed that the latter, as expected, is always compatible with zero in the continuum limit.
Third, to confirm that using the non-perturbative running of Z P yields compatible results with the one using perturbative running of Z P , we have repeated the procedure of this section applying the latter, i.e. 4-loop evolution shown in Fig. 5 -see Appendix C. As expected, the results are fully compatible, which confirms that the non-perturbatively evolved Z P differs from the perturbatively evolved one only by cut-off effects (provided, of course, that one stays at scales where PT is applicable).

Continuum limit -Strategy 2
We now present results of the other strategy of analysis, in which the values of Z P for different ensembles are not needed. Instead, one employs a rescaling procedure proposed in Ref. [23] and described in Sec Having performed this matching, the continuum values obtained from the lattice should agree with the ones of PT, i.e. the scale dependence of the anomalous dimension should agree with the one of PT for some range of scales. We find this is really the case for the range of between around 1 GeV (or even somewhat below) and 3 GeV. As discussed above, the upper bound is due to cut-off effects, while the lower one is set by physical effects -nonperturbative effects, in particular spontaneous chiral symmetry breaking effects -therefore the agreement below ca. 1.5 GeV can be coincidental.
In addition, strategy 2 gives independent estimates of the lattice spacings. Working at β ref = 3.9, we know that a β=3.9 M = 0.8 corresponds to 2.25 GeV and hence we obtain an estimate a β=3.9 = 0.071(7)(7) fm, where the first error is the combined statistical and systematic error coming from the estimate of γ m in the continuum limit (i.e. it combines the error of extraction of γ m (M ) for all 4 ensembles, the errors of r 0 /a needed to combine different lattice spacings and the error of the continuum extrapolation) and the second one is an estimate of effects of using only 1-loop PT 4 . Using β ref β=3.9, L/a=16, am=0.004 β=4.05, L/a=20, am=0.003 β=4.2, L/a=24, am=0.002 β=4.35, L/a=32, am=0.00175 cont.limit (stat.+syst.error) cont.limit (stat.error) 1-loop PT β=3.9, L/a=16, am=0.004 β=4.05, L/a=20, am=0.003 β=4.2, L/a=24, am=0.002 β=4.35, L/a=32, am=0.00175 cont.limit (stat.+syst.error) cont.limit (stat.error) 1-loop PT Nevertheless, some feeling for the size of higher order effects can be obtained by comparing γm at one and four loops in the MS scheme, which amounts to a difference of approx. 10 % at the considered energy scale. in Tab. 1, coming from chiral perturbation theory fits of the quark mass dependence of the pseudoscalar meson mass and decay constant. This agreement is reassuring, although the uncertainties are too large to be conclusive.

Determination of Λ (2) MS
The implicit scale dependence of the quark mass anomalous dimension is governed by the running of the strong coupling constant, which, in turn, depends on the value of the Λparameter of the underlying theory. The method analyzed in this paper allows to obtain the anomalous dimension with good precision over a wide range of scales and hence, as will be demonstrated below, makes it possible to determine this Λ-parameter. Specifically, it will be the Λ-parameter of 2-flavour QCD, expressed in the MS renormalization schemehence we will denote it by Λ  (Fig. 11). The two curves correspond to: • the ratio m R (µ)/m R (µ ref ) determined non-perturbatively from the extracted (continuum extrapolated) values of γ m (using (5.2) with µ 2 and µ 1 differing by 10 MeV), • 4-loop PT [45,46], where the strong coupling constant was computed using Λ MS =315(30) MeV [62] (value corresponding to a lattice determination by ETMC, from the match- ing of QQ static potential to PT, here the uncertainty of Λ (2) MS is reflected in the plot).
We observe very good agreement of the two curves within errors. The error of m R (µ)/m R (µ ref ) determined from γ m is comparable to the uncertainty of the PT curve related to the uncertainty of the used value of Λ (2) MS . Although the curves are compatible, it can be noticed that their shape is slightly different, i.e. the latter is a bit steeper, suggesting that the value of Λ (2) MS resulting from the former is actually smaller than 315 MeV (but compatible within errors).
To find the value of Λ MS that gives the best agreement between γ m -extracted Λ MS and the one from PT, we employ the following procedure. Having chosen a reference scale µ ref , we impose the matching condition at some chosen scale µ matching : MS . The total systematic error comes from combining the individual ones in quadrature. All values in MeV. See text for more details.
operator spectrum is valid, the agreement should also ensue for all other scales (down to µ such that 4-loop PT is no longer acceptable).
The above procedure is illustrated in Fig. 12. As before, µ ref = 1.5 GeV, while µ matching is chosen to be 4 GeV. The value of Λ (2) MS that leads to the fulfillment of the matching condition is 303 MeV. It is also striking that the values of m R (µ)/m R (µ ref ) from γ m and from PT agree remarkably well for the whole range of scales µ (for which the continuumextrapolated lattice data for γ m are available).
It now remains to establish the robustness of the result. We have performed a thorough error analysis. The considered sources of uncertainty are: statistical error, uncertainty from the input value of Z MS,µ 2 P (from Tab. 2, i.e. no perturbative running of the non-perturbatively found RI-MOM/X-space values is performed), error from relative and absolute scale setting, uncertainty from the arbitrary choice of µ ref and µ matching and, finally, uncertainty from neglecting higher order terms in PT.
To estimate the statistical error, a bootstrap with blocking (to account for possible autocorrelations) procedure was performed. The bootstrap procedure was also carried out to find the propagation of the error from the input values of Z MS,µ 2 P and from scale setting. For the latter, effects of relative scale setting and absolute scale setting were separated. The relative scale setting uncertainty was estimated by generating artificial bootstrap samples The results of the error analysis are gathered in Tab. 4. The relative statistical error amounts to around 4%, while the total systematic error is around 8%. The dominating sources of the latter are: the uncertainty of the values of Z P (i.e. mostly systematic errors of the non-perturbative extraction of Z P in the RI-MOM or X-space schemes), which are an external input to the present analysis and the choice of the reference scale µ ref and to a lesser extent the uncertainty of the input values of r 0 /a. Note that the absolute scale setting and the choice of µ matching lead to very small systematic uncertainties of Λ (2) MS , while the effects from neglecting higher orders of PT are essentially negligible.
Finally, we quote: Λ (2) MS = 303(13)(25) MeV, (5.4) where the first error is statistical and the second one systematic. This can be compared to other recent determinations in the literature, e.g. 315(30) MeV (ETMC, Ref. [62]), 331 (21) MeV (ETMC, Ref. [70]) and 310 (20) MeV (ALPHA Collaboration, Ref. [71]). We find very good agreement with these values and a similar total error. For a more thorough comparison with earlier determinations, see the conclusions of Ref. [62]. Note, however, that there is a tension between the two ETMC values and the ALPHA value, when the result is expressed as a dimensionless product r 0 Λ  Fig. 11. This curve gives the same information as the curve m(µ)/M in Fig. 4 of Ref. [33], where m is the renormalized quark mass and M is the RG-invariant quark mass. Although the curves are expressed in different renormalization schemes (MS and Schrödinger functional, respectively), their respective agreement with PT suggests also the mutual agreement between the results of this paper and Ref. [33].

Conclusions
The main conclusion of this work can be formulated in the following way -there exists a window in which the quark mass anomalous dimension extracted from the lattice and extrapolated to the continuum limit agrees well with continuum perturbation theory. Although the relation between the quark mass anomalous dimension and the scaling of the mode number can be shown in perturbation theory [24], it was a priori unclear whether the above mentioned window exists with presently simulated lattices.
One of the main aims of this work was to investigate whether the scaling of the mode number for intermediate eigenvalues of the Hermitian Dirac operator can be described in QCD in a similar way as in conformal field theories with an infrared fixed point. We found that the answer is indeed positive -provided one stays well above the non-perturbative regime. We have employed two analysis strategies and overall the results obtained using both of them are consistent. In both cases, we observe very good agreement with continuum perturbation theory in some range of energy scales. What is worth emphasizing is that both strategies correctly predict the scale dependence of the quark mass anomalous dimension, although they obtain it with different means -either by converting the lattice extracted value to the MS scheme using Z P in this scheme, or by matching to perturbation theory.
Taking into account all the systematic effects, we estimate that the above mentioned contact window between the lattice and perturbation theory extends from around 1.5-2 GeV to 3.5-4 GeV at presently used lattice spacings. The lower limit -ca. 1.5 GeV -originates from two physical reasons.
• In the low-energy regime, QCD with N f = 2 flavours of quarks exhibits spontaneous chiral symmetry breaking, signaled by a non-zero value of the chiral condensate. In this regime, the spectral density does not tend to zero near the origin -hence, its scaling does not follow Eq. (2.1) and the mode number obtains an additive correction from the non-zero value of spectral density at the origin. Taking realistic values of this contribution, we estimated that the influence of spontaneous chiral symmetry breaking regime extends to scales of M R ≈ 1.5 GeV.
• The lower end of the lattice window is also related to the applicability of perturbation theory at low energies. Taking the difference of 3-and 4-loop perturbative results as a proxy of this applicability, it can be estimated that 4-loop perturbation theory formula for the quark mass anomalous dimension can be trusted above around 1.5-2 GeV (with strong coupling constant α s ≈ 0.3 at this scale).
The upper limit of the window is related to the values of lattice spacings in present-day simulations and can be, in principle, improved by using finer lattice spacings. On the other hand, the lower limit is dictated by physics, i.e. non-perturbative effects, in particular originating from spontaneous breaking of chiral symmetry, that dominate at energies below the lower limit. The computation of the scale dependence of the quark mass anomalous dimension allowed also for an extraction of the Λ-parameter of 2-flavour QCD. The obtained value: where the first error is statistical and the second one systematic, agrees well with earlier determinations, in particular the recent ones by ALPHA [71] and ETM [62,70]  A Continuum limit extrapolation with the term ν 0 (m) included in the fits The influence of the term ν 0 (m) on extracted values of the anomalous dimension was discussed in Sec. 5.1. In Fig. 13, we show the outcome of an alternative continuum limit analysis -assuming that ν 0 (m) = 20 is included in Eq. (2.4). This plot should be compared to Fig. 9, which shows the results in the case ν 0 (m) = 0. Below M R ≈ 1 GeV, the effect is huge (with a non-physical maximum of γ m around this scale), signaling a total breakdown of the scaling formula for the mode number. The difference in the obtained values of γ m (M R ) at a non-vanishing lattice spacing reaches to a renormalized threshold parameter M R range of 1.5-2 GeV. However, the difference in continuum limit extrapolated anomalous  [2,4] fit [3,5] fit [4,6] fit [5,7] fit [6,8] fit [7,9] fit [8,10] fit [9,11] fit [10,12] fit [11,13] fit [12,14] fit [13,15] fit [14,16] fit [15,17]  dimension can not be seen above M R = 1 GeV (note, however, that the continuum limit extrapolation breaks down below about 1.2 GeV).
This analysis can also be repeated for different values of ν 0 (m). Finally, we can compare results with ν 0 (m) = 0, 10, 20, 30 and 40. Our conclusions are qualitatively the same in all cases -extrapolating to the continuum limit, the results are always compatible above M R ≈ 1.5 GeV, even using the conservative estimate ν 0 (m) = 40. Therefore, we conclude that 1.5 GeV is the lower limit for the applicability of the scaling formula for the mode number.
B "Artefact" anomalous dimension Γ m In Sec. 3, we discussed an alternative approach for the analysis of the scaling of the mode number, which consists in rewriting equation for the dependence of ν R (M R ) on M R at the renormalization scale µ = M R . The obtained Eq. (3.5) implies that ν R (M R ) ∝ M 4 R , i.e. the scaling exponent is independent of the quark mass anomalous dimension γ m (M R ) and equals 4, up to cut-off effects.
Here, we show the outcome of applying this alternative approach. We fitted Eq. (3.6) to the data for ν R (M R ) vs. M R and extracted the "artefact" anomalous dimension Γ m (M R ). An example for the ensemble C30.20 is given in Fig. 14. As expected, the "artefact" anomalous dimension is approximately constant in the whole range and close to zero. However, and gives further indication of the self-consistency of the presented approach.
C Continuum limit extrapolation using perturbatively evolved Z P As we have explained above, in order to compare the continuum limit of lattice extracted quark mass anomalous dimension with PT prediction, it is essential that the threshold parameter M is renormalized with Z P that is evolved without using perturbative expressions for γ m . However, as a check of self-consistency of the approach, we have performed the renormalization of M using also perturbatively evolved Z P . The results of this check are shown in Fig. 16. Minor differences between the cases of non-perturbatively and perturbatively evolved Z P are observed only around 1 GeV, although they are still statistically insignificant. We emphasize that this confirms that the curves showing perturbative and non-perturbative running of Z P differ only by cut-off effects (in the regime where PT is applicable). β=3.9, L/a=16, am=0.004 β=4.05, L/a=20, am=0.003 β=4.2, L/a=24, am=0.002 β=4.35, L/a=32, am=0.00175 1-loop 2-loop 3-loop 4-loop cont.limit (stat.+syst.error) cont.limit (stat.error) Figure 16: Comparison of the quark mass anomalous dimension extracted from the lattice to perturbation theory -Strategy 1, using perturbative running of Z P . This plot should be compared with Fig. 9, where non-perturbative running of Z P is applied. For the continuum limit results, we show the statistical error and the total one, i.e. the combined statistical error, the error originating from lattice spacing value in physical units and the error coming from Z P . β=3.9, L/a=16, am=0.004 β=4.05, L/a=20, am=0.003 β=4.2, L/a=24, am=0.002 β=4.35, L/a=32, am=0.00175 cont.limit (stat.+syst.error) cont.limit (stat.error) 1-loop PT Figure 17: Comparison of the quark mass anomalous dimension extracted from the lattice to perturbation theory, assuming O(a 2 ) scaling towards the continuum. The left plot shows the results from Strategy 1 and the right plot from Strategy 2. In the latter, the matching point is marked with a circle.
D O(a) vs. O(a 2 ) scaling R 5 -parity even quantities computed with twisted mass fermions at maximal twist are O(a)improved [49]. The situation is more complex with off-shell quantities in which contact terms can spoil the automatic O(a)-improvement -an example of such quantity is the mode number. However, it was shown in Ref. [4] that this does not happen for the case of the mode number. Still, as we discussed in Sec. 3.1, this does not imply that the quark mass anomalous dimension extracted with the analyzed method is also O(a)-improved.
In previous sections, we have performed the continuum limit extrapolations under the assumption that O(a) effects can be present and indeed we found numerically that the The results of performing continuum extrapolations using Strategy 1 and Strategy 2 and assuming O(a 2 ) scaling are shown in Fig. 17. Needless to say, both strategies lead to results in contradiction with PT in terms of values predicted at a given physical scale (Strategy 1) or the scale dependence of the anomalous dimension (Strategy 2; matching performed at around 1.6 GeV). This further confirms that indeed O(a) leading cut-off effects are present in our data.