Spectral density of the Hermitean Wilson Dirac operator: a NLO computation in chiral perturbation theory

We compute the lattice spacing corrections to the spectral density of the Hermitean Wilson Dirac operator using Wilson Chiral Perturbation Theory at NLO. We consider a regime where the quark mass $m$ and the lattice spacing $a$ obey the relative power counting $m\sim a \Lambda_{\rm QCD}^2$: in this situation discretisation effects can be treated as perturbation of the continuum behaviour. While this framework fails to describe lattice spectral density close to the threshold, it allows nevertheless to investigate important properties of the spectrum of the Wilson Dirac operator. We discuss the range of validity of our results and the possible implications in understanding the phase diagram of Wilson fermions.


I. INTRODUCTION
Simulations with Wilson fermions at light quark masses are nowadays feasible thanks to considerable algorithmic improvements [1][2][3][4][5] developed in the last years. Since chiral symmetry is explicitly broken, the spectrum of the Wilson Dirac operator is not protected from arbitrarily small eigenvalues, which might induce instabilities in numerical simulations [6]. It is therefore very important to have a theoretical understanding of the properties of the lowend spectrum of the Wilson Dirac operator. Moreover, spectral observables can be efficiently used to extract relevant quantities such as the quark condensate, like recently implemented in [7], yielding a further strong motivation to investigate the impact of lattice artefacts on the eigenvalues spectrum.
In the continuum and in an infinite volume, the (renormalised) spectral density of the massive Dirac operator has a threshold given by the (renormalised) quark mass. Lattice artefacts are expected to change both the location of the threshold and the shape of the spectral density close to the threshold. Moreover, in a finite box of volume V one expects that when γΣV 1, where γ is an eigenvalue of the massless Dirac operator and Σ the chiral condensate, finite-size effects become important and can also induce relevant deformations of the spectral density with respect to the infinite-volume case.
From the numerical point of view, some information about the low-end spectrum is provided by lattice simulations with Wilson fermions carried out in the past years. For N f = 2 unimproved Wilson fermions, empirical observations [6] indicate that the median of the spectral gap distribution is linearly proportional to the quark mass m, while the width is basically independent on m and scales like ∼ a/ √ V , where a is the lattice spacing. On the other hand, for the O(a)-improved theory, the situation is less clear [8] and those properties have not been confirmed. In [7] the mode number of the O(a)-improved Wilson Dirac operator is computed, finding a nearly linear behavior up to 100 MeV above the threshold. At the low end of the spectrum, a significant deviation from the continuum expectation is observed.
From the theoretical side, Wilson chiral perturbation theory [9,10] (WχPT) is the tool which provides a systematic description of low-energy properties of lattice Wilson QCD including the leading discretisation effects. When approaching the chiral limit at finite lattice spacing, one enters in the regime where m ∼ a 2 Λ 3 QCD , which is where discretisation effects compete with the quark mass to the explicit breaking of chiral symmetry. Lattice artefacts induce a non-trivial phase diagram, and two different scenario have been foreseen: in the so-called Aoki scenario there is a range of quark masses (in the Aoki phase [11]) where there are two massless pions. On the other hand, in the so-called Sharpe-Singleton scenario [9] there is a first-order phase transition, and the three pions remain massive in the chiral limit. From the point of view of WχPT, in this regime the discretisation effects appear already at leading order in the chiral expansion. Lattice artefacts in the infinite volume spectrum of the Hermitean Wilson Dirac operator have been computed in this regime in [12], although a working framework has been only found by assuming additional conditions on the couplings of WχPT associated to discretization effects.
In a recent study [13,14] the same power counting for a has been adopted, but in a finite-box in the -regime 1 . In this case the LO predictions of WχPT can be obtained also by means of a "modified" Random Matrix Theory which includes O(a 2 ) effects. The main difficulty of this computation is that it involves exact integrals over the zero modes which must be defined at fixed topological charge introduced via the number of real modes of the Dirac operator, which nevertheless has an intrinsic ambiguity for Wilson fermions at finite lattice spacing.
In this work we will study the discretisation effects in the spectrum of the Dirac operator in a different regime. We consider WχPT with the power counting m ∼ aΛ 2 QCD (GSM regime), where lattice artefacts appear only at next-to-leading order in the chiral effective theory, and can hence be treated as a perturbation. This simplifies considerably the computation; it allows nevertheless to extract important information about the impact of lattice artefacts in the spectral density.
This paper is structured as follows: in sec. II we recall definitions and properties of the spectral density of the Wilson Dirac operator; in sec. III we present the setup for WχPT needed for the computation; in sec. IV we give the details about the calculation of the spectral density, including the finite-volume corrections (in the p-regime); in sec. V we discuss our results and show a comparison with numerical data. All technical details about the computation are deferred in the Appendices.

II. SPECTRAL DENSITY FOR THE WILSON DIRAC OPERATOR
We start considering QCD for N f = 2 degenerate quarks in infinite volume with current quark mass m. The average spectral density of the massless Hermitian Dirac operator −iD can be defined as where γ k are the eigenvalues of the massless Hermitian Dirac operator and . . . indicates the usual path-integral average. The Banks-Casher relation [18] tells us that the spectral density is related to the chiral condensate Σ in the following way One useful way to determine this relation in continuum QCD [19,20] consists in adding a valence quark ψ v of mass m v to the theory. Using the spectral decomposition for the valence chiral condensate one obtains This relation can be inverted because the spectral density is independent on the valence quark mass where Disc indicates the discontinuity across the imaginary valence quark mass. Strictly speaking eq. (3) is ultraviolet divergent also after renormalisation of the quark masses and the gauge coupling constant and one must introduce a cutoff in the integration range. The ultraviolet divergences turn out to cancel when computing the discontinuity of the valence scalar condensate in eq. (4). Thus it is enough to compute the valence condensate for real masses and analytically continuing the resulting expression for complex masses. Eq. (4) can be used naturally in chiral perturbation theory and one obtains [20] the Banks-Casher relation (2) and the NLO corrections [7,21]. We now discretise our N f = 2 continuum QCD action on a lattice of spacing a and we consider Wilson fermions. All our considerations can be generalised to the case of Wilson twisted mass. To study the spectral density of the Wilson operator one has to take into account that the Wilson operator D W is neither Hermitian nor anti-Hermitian and it has complex eigenvalues. It is thus advantageous to define the Hermitian Wilson-Dirac operator We indicate now the spectral density of Q m as ρ Q (λ, m) and the one of Q 2 m as ρ(α, m), where λ and α label the eigenvalues of the corresponding operators.
If the continuum theory has been regulated on a lattice the ultraviolet divergences of the chiral condensate appear as power law and logarithmic divergences in the lattice spacing a. Moreover, when using Wilson fermions the lack of chiral symmetry implies that the scalar condensate is not directly related to the spectral density as in the continuum. It is thus not obvious how to work out the renormalisation and O(a) improvement of the spectral density starting from eq. (4) or from the corresponding version for Wilson lattice QCD. To study the renormalisation and the O(a) improvement of the spectral density one needs to relate it to correlation functions of local operators where standard arguments on renormalisability and O(a) improvement can be applied.
It has been shown in ref. [7], using chain correlators of scalar and pseudoscalar densities how the spectral density of Q 2 m renormalises. In the following we will be mostly interested in the spectral density of Q m that renormalises as follows where m R is the renormalised quark mass. Moreover, additionally to the standard improvement of the action and local operators there are O(am) cutoff effects which need to be removed to fully improve the spectral density [7]. This analysis guarantees that with Wilson fermions the spectral density is a well defined renormalisable quantity and with a well defined Symanzik expansion. We can thus compute the spectral density in the chiral effective theory using a generalisation of eq. (4) for Wilson chiral perturbation theory.
A further spectral observable which can be defined and measured in a lattice simulation is the integrated spectral density which represents the density of modes in the interval between Λ 1 and Λ 2 ; it contains the same physical information as the spectral density and satisfies N R (Λ 1,R , Λ 2,R , m R ) = N (Λ 1 , Λ 2 , m), i.e. it is a renormalisation-group invariant [7].
To relate a partially quenched condensate with the spectral density of Q m it is convenient to introduce a doublet of degenerate twisted mass fermions χ v with a mass term iµ v γ 5 τ 3 , where τ 3 is the third Pauli matrix. The relation between the condensate and the spectral density is now different but can be worked out and it reads [12] The need to symmetrise in λ the spectral density in the numerator of (8) can be understood from the lack of chiral symmetry of the Wilson operator which renders the spectral density not symmetric when λ → −λ. Inverting eq. (8) one obtains where the untwisted mass in the valence sector coincides with the sea quark mass m. Eq. (9) is our starting point for the computation of the spectral density in Wilson chiral perturbation theory. We stress that in chiral perturbation theory the problem of power law divergences is absent because those are mapped to the presence of the so called high-energy constants usually indicated with H. As we will see in sect. IV this is indeed the case for the expectation value of the pseudoscalar condensate in the chiral effective theory. This implies that in the effective theory the computation of the discontinuity is a well defined procedure. It remains to make sure that the spectral density computed in the effective theory through the discontinuity of the partial quenched chiral condensate does not overlook additional cutoff effects which go beyond the standard O(a) corrections stemming from the action or the local operators. Our previous discussion and the results of ref. [7] guarantees that the only additional cutoff effects are of O(am) and as we will see in sect. III these cutoff effects are of subleading order in our power counting of Wilson chiral perturbation theory.

III. SPECTRAL DENSITY IN WILSON CHIRAL PERTURBATION THEORY
In order to investigate the discretisation effects in the spectral density of the Hermitian Wilson Dirac operator we compute the valence pseudoscalar density in the framework of Partially Quenched Wilson chiral perturbation theory (PQWχPT). According to eq. (8), we consider N f = 2 sea quarks with bare mass m and a doublet of valence twisted mass quarks with (m + iµ v τ 3 ). We formulate the effective theory on a finite volume V = L 3 T and we keep the quark masses in the p-regime, corresponding to the power counting m, µ v ∼ O(p 2 ), 1/L, 1/T ∼ O(p) (10) in terms of the momenta p. The extension of the effective theory to the case of non-zero lattice spacing is done in two steps: after matching the lattice QCD action with the appropriate Symanzik continuum effective action [22,23], one writes down a chiral Lagrangian which contains the standard continuum terms plus additional operators that transform under chiral symmetry as the operators of the Symanzik effective theory. For the Wilson action, this has been studied in [9,10]. When introducing lattice artefacts in the chiral effective theory, we have to define as usual the relative power counting between the quark mass and the lattice spacing a. In this work we adopt a counting corresponding to the so-called GSM (Generically Small quark Mass) regime [24,25], where m, µ v ∼ aΛ 2 QCD : in this case the explicit breaking of chiral symmetry is dominated by the quark mass, and lattice artefacts can be treated as perturbations.
Partially quenching can be implemented in the Chiral Effective theory by means of two techniques, namely the graded-symmetry method [26,27] and the replica method [28]. In the first one, one introduces "ghost" quarks, whose determinant cancels the one of the valence quarks. In the second one, one enlarges the valence sector to N r flavors and one eliminates the corresponding determinant by sending N r → 0. The equivalence of the two methods has been shown at the perturbative level in [28]. In this section we will focus on the gradedsymmetry method; we checked nevertheless that the same results can be obtained with the replica method, which we summarise in App. B. For our specific case we have to consider a chiral Effective Theory with a graded symmetry group SU(4|2) L ×SU(4|2) R spontaneously broken to SU(4|2) R+L . The formulation we adopt for our computation using the SU(4|2) effective theory is based on an extension of the framework proposed in [29] for the SU(3|1) case. We consider a mass matrix of the form The quark mass in the ghost sector has been already set equal to the mass in the valence sector. Moreover, the untwisted part of the valence quark mass has to be equal to the sea quark mass. We parametrise the pseudo Nambu-Goldstone bosons by the field U ∈ SU(4|2) where F is as usual the pseudoscalar decay constant, and T a , a = 1, . . . , 35 represent the generators of the corresponding Lie algebra, which satisfy We refer to App. A for the explicit form of g ab (eq. (A15)), and for a summary of conventions and properties of the SU(m|n) group. The constant field u V represents the ground state of the theory, which can be obtained by minimising the potential in the LO Chiral Lagrangian. The solution of the potential minimisation yields 2 where In the language of twisted mass Wilson theory, m P is the so-called polar mass.
Taking into account our power counting, the full NLO Chiral Lagrangian in PQWχPT can be written as In the GSM regime, the LO Lagrangian can be written in the same form as the continuum one provided we substitute the quark mass with the so-called shifted mass [9], which incorporates the leading O(a) corrections. In the following we assume that m in eq. (11) includes already this shift. The coupling B is related to the chiral condensate Σ, B = Σ/F 2 . From L 2 we can extract the propagator where is the infinite-volume contribution, while g r (M 2 ) represents the finite-volume correction and is UV-finite. A specific representation of g r (M 2 ) will be considered in App. C. The matrix h ab is defined in App. A, eq. (A17). The mass term M 2 a appearing in the propagator is given by 2 In the graded symmetry formulation of PQχPT there are subtleties concerning the minimisation of the potential. We will discuss them in more detail in sec. IV.

The continuum NLO PQ Lagrangian reads [31]
where we have written down only terms which can enter in our specific computation. The PQ Chiral Lagrangian encoding discretisation effects up to O(a 2 ) can be written as [10,32,33] whereâ = 2W 0 a has dimension [energy] 2 and W 0 is the LO low-energy coupling absorbed in the shifted mass m. Also in this case we have disregarded terms which are not relevant for our computation.

IV. CALCULATION OF THE SPECTRAL DENSITY WITH SU(4|2) GRADED GROUP METHOD
In this section we compute the spectral density at NLO in WχPT, making use of the definitions and assumptions introduced in the previous section. The observable from which we start is the partially quenched twisted pseudoscalar condensate made of valence quarks defined in eq. (8). In order to extract it we have to introduce source terms in the Chiral Lagrangian by means of where M is the mass matrix introduced in eq. (11), andτ 3 has non-zero components only in the valence sectorτ The pseudoscalar density is obtained by deriving the action associated to the Chiral Lagrangian in eq. (16) with respect to p 3 At LO (O(p 0 )) we obtain the (continuum) expectation value The NLO result is where δ is an O(p 2 ) quantity that represents the correction to the ground state angle ω 0 due to lattice artefacts. In the following subsection we discuss the computation of this effect.

A. NLO correction to the vacuum
The O(p 2 ) terms in the chiral Lagrangian give rise to a shift in the vacuum angle ω 0 [25], which must be computed by minimising the NLO potential. An important fact is that at this order the continuum L i do not contribute to this realignment, and only the discretisation effects must be taken into account. Since this is an important point of the calculation, we start by recalling how the continuum ground state given in eq. (14) is obtained. Using the equations of motions one can show that u V has to commute with the matrix M † M, which implies that u V is diagonal. Therefore we can parametrise it as where π < φ 1 , . . . , φ 4 ≤ π are standard phase factors. The condition Sdet(u V ) = 1 requires As suggested in [34], we perform an analytic continuation for the ghost components where nowφ 5,6 are real variables taking values along the real axis. The LO potential can be separated in quark and ghost components with The potential is complex, and a minimisation procedure in the usual sense is not applicable. The minimisation of the potential in Euclidean field theory is equivalent to perform a saddlepoint expansion of the functional integral, and the saddle-point expansion can be performed even if the potential is complex. Our prescription, following [35], is to find the saddles for complexφ 5,6 and then deform the contour of integration for each field variable in order to pass through the saddle points. The saddle points are chosen following two criteria: a) in order to maximise the value of the real part of the potential at the saddle; b) in order to have the direction where the real part of the potential rises more steeply consistent with the possibility of deforming the contour integration without encountering the other saddle points. This procedure has been discussed in detail in [35] in the context of the phase diagram of quenched Wilson ChPT. We find that both criteria for the choice of the appropriate saddle points are fulfilled by φ 5 = −φ 6 = ω 0 . In this way one obtains the expected result given in eq. (14). We now repeat this procedure including the NLO lattice artefacts. We consider a ground state u V, N LO of the form given in eq. (30) and we perturb the LO solution by choosing where δ s,v,g are corrections of O(p 2 ). We then perform the same analytic continuation on the ghost components, By expanding the NLO potential at O(p 4 ) we obtain with a 1 = 0 a 2 = mΣ, a 3 = ia 5 = 16â sin w 0 (âW 8 cos w 0 + 2âW 6 + Bm P W 8 + 2BmW 6 ), a 4 = a 6 = m P Σ.
Also in this case the potential can be split in a real (sea + valence) part V q,N LO and a complex (ghost) part V g,N LO . The saddle point expansion yields The ground state of the theory at NLO is then given by with Since the discretisation effects represent a perturbation in our power counting, this solution does not impose constrictions on the absolute value or the sign of the LECs W 6 , W 8 , W 6 , W 8 .

B. Spectral density in infinite volume
The spectral density can be computed from P 3 by taking the discontinuity along imaginary µ v , as given in eq. (9). At LO we obtain the known result which does not depend on the volume and on the lattice spacing in our specific regime. The NLO result in infinite volume is given by where∆ arises from the correction to the ground state and is given bỹ If we define a shifted sea quark mass as followŝ under the condition that λ m, the correction proportional to∆ can be resummed. In this case one could then rewrite eq. (46) by omitting the term proportional to∆ and substituting m →m everywhere. This resummed formula would be equivalent to the one in eq. (46) up to higher order corrections. Notice that the termâ/M 2 ss in eq. (47) is not singular in the GSM regime. The UV divergences arising from G 1 (0, M 2 ) and G 2 (0, M 2 ) have been cancelled by defining renormalised couplings [36]. Using dimensional regularisation in 4 − 2 dimensions and adopting the convention of [7] we define where µ is a renormalisation scale. The integrated spectral density N (Λ 1 , Λ 2 , m) defined in eq. (7) can then be computed straightforwardly. For better readability we report it in App. C, eq.(C1). We stress that in the continuum these expressions coincides with the ones computed in [7]. It is useful to express our results for the spectral density in terms of the PCAC quark mass.
The NLO lattice corrections to the PCAC mass can be computed in the conventional WχPT for N f = 2. The result is [25] m PCAC = m 1 + 16â where W 10 is an extra LEC. Eq. (46) can be then rewritten by substituting everywhere m → m PCAC and∆ → ∆, with As before, the correction given by ∆ can be resummed; the shift in the quark mass given in eq. (48) can be rewritten as a shift in the PCAC mass: Notice that at NLO we can consistently substitute m → m PCAC in M 2 ss . In fig. 1 we plot [ρ Q (λ, m) + ρ Q (−λ, m)] N LO in the continuum theory and in the discretised case with W 6 = W 8 = W 10 = 0, corresponding to a non-perturbatively O(a)-improved theory. We consider in particular a quark mass m PCAC = 26.5 MeV; the values of other parameters used in this plot are specified in the caption and justified by our numerical analysis discussed in sect. V. The central black curve corresponds to the continuum case; W 8 > 0 (∆ < 0) gives rise to the red curve, while W 8 < 0 (∆ > 0) corresponds to the blue curve. The dashed lines for λ 40 MeV are a reminder that our result is not expected to be valid close to the threshold. This issue will be discussed in more details in sec. V.

C. Finite volume corrections
Finite-volume effects arise at NLO and can be computed by taking into account the corrections to the propagators given in eq. (20). We can write down the spectral density as sum of the infinite volume and finite-volume corrections: The explicit expression for the finite-volume correction ∆ρ V Q (λ, m) + ∆ρ V Q (−λ, m) is given in App. C, eq. (C4). By integrating over λ we obtain the corresponding finite-volume corrections for the integrated density In eq. (C5) we give the formula for ∆N V (Λ 1 , Λ 2 , m) for the particular case Λ 1 = m. We stress that in our power counting finite-volume effects are insensitive to lattice artefacts at NLO, and these results are the one of the continuum theory already obtained in [7]. Notice that the finite-volume correction N V (Λ 1 , Λ 2 , m) diverges at the threshold (Λ 1 = Λ 2 = m): this could be interpreted as a signal that, even if the quark mass is in the p-regime, another scale comes into play in this problem, namely √ λ 2 − m 2 , and it needs to be treated with a different power counting when approaching the threshold at finite volume. Since in the -expansion of the chiral theory divergences naturally appear when taking the chiral limit at finite volume, this might be the correct power counting to adopt for this additional scale when √ λ 2 − m 2 ΣV ∼ 1. Having said that, for L 2.0 fm and sufficiently far from the threshold, finite-volume corrections amount to few percents.

V. DISCUSSION OF THE RESULTS AND CONCLUSIONS
The important results of this paper are the expression of the spectral density eq. (46) and its integrated form (C1). They are both obtained in the p-regime of WχPT at NLO. Eq. (46) shows that the spectral density, when computed on a lattice with Wilson fermions, is modified in two aspects with respect to the NLO continuum formula (see for example refs. [7,37] or simply set all the W s in eq. (46) to zero). The lattice artefacts modify the behaviour of the spectral density near the threshold, i.e. the term proportional to∆ in (46), and and its absolute normalisation, i.e. the term proportional to W 6 in (46). The absolute normalisation correction proportional to W 6 vanishes if the theory is non-perturbatively O(a)-improved. If the theory is not improved, beside the term proportional to W 6 , the correction to the normalisation contains in principle an additional O(a) term stemming from the renormalisation constant Z P . We prefer to omit in our final formula (46) this term because it strictly depends on the way Z P is computed and only in very few cases we can have a suitable representation in WχPT [38]. We stress nevertheless that when analysing the spectral density computed in the unimproved theory using eq. (46) the O(a) cutoff effects of Z P have to be taken into account in some way.
More interestingly the cutoff effects modify, the λ and m dependence of the continuum formula. These corrections, expressed in terms of the PCAC mass (see eqs. (50-52)), depends on the following LECs: W 8 , W 10 and W 8 . The modifications induced by cutoff effects to the continuum formula are plotted in fig. 1. To simplify the discussion we have chosen a lattice setup where the theory has been non-perturbatively improved, thus we set all the O(a) LECs W 6 = W 8 = W 10 = 0. With the black curve we plot the λ dependence of the spectral density in the continuum as given by NLO χPT [7,37] at given values for m, Σ and F (see the caption of the plot). In the same plot we also show the modifications to the continuum formula induced by WχPT at NLO. The first observation we make is that the way the spectral density depends on λ, close to the threshold, is rather different from the continuum formula. The second observation we make is that the WχPT formula has a rather different behaviour depending on the sign of W 8 , especially close to the threshold. If ∆ > 0 when λ → m PCAC the WχPT formula has a non-integrable singularity, while the continuum formula has an integrable singularity. If ∆ < 0 the spectral density goes to zero for a value of λ > m PCAC . These behaviours are shown in fig. 1. This λ dependence of the spectral density near the threshold is totally dominated by the lattice NLO corrections induced by the vacuum realignment, indicating that our expansion cannot be trusted for too small λ. Our power counting breaks down near the threshold: more specifically, it is valid only when √ λ 2 − m 2 ∼ aΛ 2 QCD , i.e. when the scale √ λ 2 − m 2 obeys a GSM power counting. By moving towards the threshold at finite lattice spacing one enters in the region where √ λ 2 − m 2 ∼ a 2 Λ 3 QCD (Aoki regime), and lattice spacing corrections can not be treated in a perturbative way as it is done in this work. In fig. 1 we draw dashed line where we assume our formula for the spectral density to break down. The minimal λ 40 MeV is just an estimate obtained fitting the available numerical data, as we will discuss futher on. It seems a reasonable estimate given the fact that for λ 40 MeV the relative lattice spacing NLO corrections are still sufficiently small. Of course this is just an estimate based on the value of ∆ we input, i.e. from the size of the lattice artefacts. Only a thorough analysis of the numerical data can determine the value of ∆, thus the range of validity of our results.
Thus we can try to use our formula (C1) to fit the available numerical data and see if we obtain reasonable estimates of Σ and ∆. To compare our formula with the available numerical data from ref. [7] we decided to integrate the spectral density from λ = Λ 1 and λ = Λ 2 , taking Λ 1 = 40 MeV, in order to avoid the small λ region where the lattice corrections to the χPT formula dominate. In order to compare with the results of [7] we use the mode number defined as ν(Λ 1 , Λ 2 , m) = V N (Λ 1 , Λ 2 , m) and we subtract the value ν(0, Λ 1 , m) from the numerical data. In fig. 3 we show a typical fit of the numerical data of ref. [7] using the infinite volume formulae of continuum χPT and WχPT in the range Λ 2 ∈ [40 : 60] MeV. The numerical data are obtained at a value of the renormalised PCAC mass m R 26 MeV. We have performed fits using the both the NLO formulae of χPT and WχPT. 4 With the continuum χPT formula we have fit Σ and fixed the following parameters: µ = 139.6 MeV, F = 80, 90 MeV,L 6 = 3, 4, . . . , 7 and m PCAC = 26.5 MeV. With the WχPT formula we have set all the O(a) LECs W 6 = W 8 = W 10 = 0 because the numerical data are nonperturbatively O(a)-improved. Additionally we have fixed the value of Σ 1/3 within the range 250 − 300 MeV and we have fit ∆. We have performed fits with several forms for the WχPT formula, all differing by NNLO terms. We have observed that the best fits are obtained with eq. (C1). In order to decide on the 'best fits' we have selected the fit results that gives small values of the squared differences between the input data points and the theoretical formula evaluated at the same points.
A summary of the fit results, satisfying our criterion is given in fig. 4. The plot shows the parameter space of the fit results in the ∆, Σ 1/3 plane. Each coloured band represents a fixed value ofL 6 fig. 3 gives ∆ = −0.55 with a value of Σ 1/3 = 273 MeV, which is in the right ballpark. We will come back to the consequences of this result at the end of this section. We remark that while the fit results obtained using the continuum NLO formula are not disastrous always a better agreement is found using our WχPT result. The WχPT formula is able to describe well the numerical data also in a larger Λ 2 range, i.e.
[40 : 100] MeV, with reasonable values of the chiral condensate. Even if our analysis is rather qualitative and more numerical data are needed, we can try to draw some conclusions. WχPT describes the numerical data better than continuum χPT in the range of applicability of our results, i.e. sufficiently away from the threshold, but we cannot exclude that adding NNLO terms in the continuum formula would improve the quality of the fits. Our fit results seem to prefer a negative value for ∆ implying a positive value for W 8 . In the literature, the LECs W 6 and W 8 are often found combined as [39] which is the combination appearing in quantities computed in full N f = 2 WχPT. The chiral phase diagram for Wilson-type fermions depends on the LEC c 2 (55) and in particular two scenarios can take place depending on its sign. In our case only one of the two LECs appears in the discretisation correction to the spectral density. For large number of colours N c , one finds W 6 /W 8 ∼ 1/N c [12], indicating that W 8 is the dominant contribution in the coefficient c 2 . This implies that independently on the sign of W 6 a positive value for W 8 would suggest a negative value for c 2 . 5 Assuming for instance W 6 = W 8 /3, then the results we obtain from the fit ∆ = −0.55 imply thatâ 2 W 8 = 5.4 · 10 6 MeV 4 which for a = 0.08 fm would correspond to c 2 −(473 MeV) 4 . It would be interesting to extend the numerical data in order to establish more accurately the sign of and the absolute value of W 8 . This will give us valuable informations on the behaviour of the spectrum of the Hermitean Wilson operator and a hint on the scenario for the chiral phase diagram that takes place with clover fermions and Wilson gauge action. 6 We close this section by stressing again an important remark. Our results seem to indicate that apart from the quark mass, the lattice spacing and the lattice size, there is an additional scale entering in this problem, namely the eigenvalues of the massless Dirac operator. In particular, after the analytic continuation of eq. (9) the mass m 2 + µ 2 v becomes √ λ 2 − m 2 , which parametrizes the "distance" from the threshold of the spectral density. Even if we assume, as it is done in this work, that all masses are in the p−regime the value of this additional scale can be arbitrarily small. In order to get meaningful results we have to implicitly assume that this parameter obeys the same power counting as the quark masses with respect to the lattice spacing a and the lattice size L. While this fact seems trivial at this stage, it is not totally obvious to infer it from the partially quenched initial setup, where valence quarks are introduced as probes to obtain the spectral density. This implies that, at fixed a and L, our results will not be valid in the vicinity of the threshold, since from one side the p-regime expansion will fail, and from the other side the GSM counting will not be valid anymore, because close to the threshold, as we have discussed before, the NLO corrections induced by the finite lattice spacing dominate the LO result.
We stress that these two effects are decoupled: close to the threshold, there is a scale which in continuum χPT at finite volume needs to be treated with the -expansion, and which in WχPT in infinite volume needs to be treated with the Aoki power counting. The way the two effects combine close to threshold is a not trivial issue and it deserves further investigations.
We conclude that our formula is a useful tool, combined with numerical data, to understand the behaviour of the spectral density of the Hermitean Wilson operator. We also consider this work a necessary and important step towards a complete theoretical understanding of the behaviour close to the threshold of the spectral density of the Hermitean Wilson operator.

VI. ACKNOWLEDGMENTS
We are indebted to Leonardo Giusti for the access to some crucial private notes on the graded-group method for partially quenched Chiral Effective Theory and for useful discussions. We thank Steve Sharpe and Leonardo Giusti for a critical reading of the manuscript and for precious comments, as well as Martin Lüscher for advices and suggestions and for providing numerical data for the spectral density of the Wilson Dirac operator.
In this Appendix we collect conventions and properties related to SU(m|n) [29,43] which are relevant for the partially quenched Chiral Effective Theory in the graded-group formulation.

Supermatrices
For n, m positive integers, a square even supermatrix is defined as a (m + n) × (m + n) matrix with the structure where A and D have elements in the even subspace of the Grassmann algebra, while B and C have elements belonging to the odd subspace of the Grassmann algebra (see [43] for definitions and properties of the Grassmann algebras). We define the supertrace

Lie supergroups and superalgebras
The Lie Supergroup SU(m|n) is the set of (m + n) × (m + n) even supermatrices U satisfying the conditions By parametrizing U as these conditions convert into The Lie superalgebra associated to SU(m|n) can be constructed by supplementing the associative superalgebra of (m+n)×(m+n) complex matrices defined in the previous subsection with the generalised commutators and zero supertrace Str(T ) a = 0.
The normalisation is given by where with g ab = (−1) deg(a)deg(b) g ba . The elements from m 2 to m 2 + 2mn − 1 are odd matrices, while the remaining ones are even.
In this work we consider the particular case of SU(4|2). With our conventions, T 1 , . . . , T 15 are the (even) generators of the SU(4) subgroup that includes sea and valence quarks, T 32 , . . . , T 34 (even) are associated with the SU(2) ghost sector, T 16 , . . . , T 31 (odd) mix the ghosts and the quarks, and finally T 35 (even) is a diagonal matrix with components in both quark and ghost sectors.
The pseudo Nambu-Goldstone fields are parametrised by where now T a , a = 1, . . . , (2 + N r ) 2 − 1 are the generators of SU(2 + N r ), with the normalisation convention while u V represents the ground state of the theory. By minimising the LO potential we obtain with the definitions already given in eq. (15) sin In this framework it is convenient to write down the pseudo Nambu-Goldstone propagator in explicit components The observable we have to consider in order to extract the spectral density is the expectation value of the pseudoscalar density defined in eq. (8). The pseudoscalar density can be obtained by taking the functional derivative of the action with respect to appropriate sources, which are introduced by the following procedure: whereτ 3 has non-zero elements only in one of the replica of the valence sector This result is independent on N r , and coincide with the one obtained with the graded group method, eq. (28). At O(p 2 ) we obtain − 8â (W 6 (2 + N r cos ω 0 ) + W 8 cos ω 0 ) The shift δ is the O(p 2 ) correction to the ground state angle ω 0 that must be computed by minimising the NLO potential. Unlike the graded group case discussed in sec. IV A, with the replica method the minimisation is trivial, since the fields belong to the conventional SU(N s + N r ) group. By taking only the linear term in N r we get which coincides with the shift calculated with the graded groups in sec. IV A, eq. (44). The final result in the limit N r → 0 for the partially quenched pseudoscalar density at NLO is then This result for the pseudoscalar density is fully equivalent to the one obtained with the graded groups, eq. (29), as expected.

Integrated spectral density at NLO
We report the integrated spectral density computed at NLO in WχPT. Starting from the result presented in eq. (46) for the spectral density, we obtain: where∆ is given in eq. (47).