Continuum limit of the mobility edge and taste-degeneracy effects in high-temperature lattice QCD with staggered quarks

We study the effects of taste degeneracy on the continuum scaling of the localization properties of the staggered Dirac operator in high-temperature QCD using numerical simulations on the lattice, focusing in particular on the position of the mobility edge separating localized and delocalized modes at the low end of the spectrum. We find that, if the continuum limit is approached at fixed spatial volume, the restoration of taste symmetry leads to sizeable systematic effects on estimates for the mobility edge obtained from spectral statistics, which become larger and larger as the lattice spacing is decreased. Such systematics, however, are found to decrease if the volume is increased at fixed lattice spacing. We argue that spectral statistics estimate correctly the position of the mobility edge in the thermodynamic limit at fixed spacing, and support this with an independent numerical analysis based directly on the properties of the Dirac eigenvectors, that are unaffected by taste degeneracy. We then provide a theoretical argument justifying the exchange of the thermodynamic and continuum limits when studying localization. This allows us to use spectral statistics to determine the position of the mobility edge, and to obtain a controlled continuum extrapolation of the mobility edge for the first time.


INTRODUCTION
The microscopic mechanism behind the finitetemperature transition in QCD is still the subject of intense research activity.One of the open questions is how the approximate restoration of chiral symmetry and the deconfinement of quarks and gluons into the quark-gluon plasma, both taking place as a rapid crossover in the same range of temperatures [1,2], are related to each other.
An intriguing aspect of the transition, revealed by firstprinciples nonperturbative numerical studies on the lattice, is that it is accompanied by a radical change in the localization properties of the low-lying eigenmodes of the Dirac operator, from extended over the whole system below the crossover, to localized on the scale of the inverse temperature above the crossover [3][4][5][6][7][8][9][10][11] (see Ref. [12] for a review).At high temperatures, a "mobility edge" separates localized low modes and delocalized bulk modes; this decreases as a function of temperature, eventually vanishing at a temperature in the crossover range.A similar behavior has been observed in other gauge theories, both pure gauge [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27] and in the presence of dynamical fermions [28,29] and scalars [30], with the mobility edge vanishing at the critical point when the low-and hightemperature regimes are separated by a genuine phase transition [19-22, 25, 26, 28, 29].
While the Dirac eigenmodes obviously encode the whole dynamics of quarks and antiquarks interacting * claudio.bonanno@csic.es† giordano@bodri.elte.huwith the non-Abelian gluon fields, their relation with physical observables is often not straightforward, as observables are generally obtained by integrating suitable eigenvector correlators over the whole spectrum.A notable exception is the density of near-zero modes, that in the chiral limit entirely determines whether or not a chiral condensate is developed [31], and even at physical quark masses is mainly responsible for the fate of chiral symmetry.On the other hand, it has been argued that the change in the localization properties of nearzero modes is mostly due to the ordering of Polyakov loops above the transition [12,18,26,[32][33][34], hinting at a close relation between localization and the confining properties of the theory.Understanding the behavior of the low Dirac modes across the transition is then key to unveiling the connection between chiral symmetry restoration and deconfinement.A thorough investigation of the "geometric" transition associated with the appearance of localized low modes and of the corresponding mobility edge in the Dirac spectrum, and of its relation with the thermodynamic crossover, cannot dispense with the extrapolation to the physical, continuum limit.Although so far no dedicated study of the continuum scaling of localization properties has been performed in the literature, there are both theoretical arguments [4,35] and numerical evidence that low-mode localization is not a lattice artifact, and survives the continuum limit.In general, localization in QCD has been found with a variety of discretizations of the Dirac operator (including staggered [3][4][5]7], Möbius domain wall [8], overlap on twisted mass [9,10] and on clover [11]).Moreover, the mobility edge in units of the arXiv:2312.02857v2[hep-lat] 12 Mar 2024 quark mass is a renormalization-group invariant quantity [4,35], that numerically shows little dependence on the lattice spacing [4,29].Such numerical evidence, however, relies on the connection between the localization properties of eigenvectors and the statistical properties of the corresponding eigenvalues [36], and was obtained using the staggered discretization of the Dirac operator.This combination is potentially problematic.
On the one hand, a lattice Dirac operator in the background of the fluctuating gauge fields appearing in the path-integral formulation of gauge theories is formally a (sparse) random matrix.As such, its eigenvalues are expected to display the universal types of correlations appearing in these systems [37,38].This allows one to identify the mobility edge as the point in the spectrum where the spectral statistics switch between the universal type corresponding to localized modes, and the universal type corresponding to delocalized modes (which depends only on the symmetry class of the lattice Dirac operator according to the random matrix theory classification [38,39]).On the other hand, in four space-time dimensions the staggered operator describes four "tastes" of fermions that become exactly degenerate in the continuum limit, which in turn leads the staggered spectrum to develop degenerate quartets of eigenvalues as the lattice spacing decreases toward zero.The obvious additional correlations between the nearly-degenerate eigenvalues appearing on fine lattices then spoil the expected universal behavior.These deviations are difficult to control theoretically, and make it more difficult to reliably determine the localization properties of eigenmodes using spectral statistics on finer lattices.This issue has been known for a long time [4,40,41], 1 but it has never been fully addressed in the lattice literature so far.A careful lattice study would be extremely interesting to assess possible systematics affecting the calculation of the mobility edge.
In this respect, an important point to consider is that the influence of the approximate taste degeneracy on the spectrum is reduced as the volume is increased at fixed lattice spacing.Indeed, while the typical distance between would-be-degenerate modes within the same multiplet is controlled by the lattice spacing, the expected typical distance between neighboring eigenvalues is controlled by the inverse of the spectral density (i.e., the number of Dirac modes per unit interval in the spectrum), and so by the inverse volume.When multiplets are clearly distinguishable, this quantity actually mea- 1 For certain gauge groups and choices of representation for the fermions, the symmetries of the continuum Dirac operator differ from those of the staggered one [38]: In these cases, the deviation from the expected universal behavior is actually needed, in order for the statistical properties of the lattice spectrum to approach those of the continuum spectrum (once the degeneracy is removed).However, deviations from the universal behavior will appear even when the staggered and continuum operator are in the same symmetry class.
sures the distance between them (up to the degeneracy factor).As soon as the lattice size becomes large enough and this distance becomes comparable with the multiplet splitting, the multiplets overlap and the corresponding structure gets washed away from the spectrum.However, taste symmetry should still affect the short-distance correlations between neighboring modes, as these still carry a remnant of the multiplet structure.Finally, when the inverse spectral density becomes much smaller than the splitting within would-be multiplets, taste symmetry effects disappear from short-distance correlations; the inverse spectral density measures now the distance between neighboring eigenvalues that, loosely speaking, belong to different multiplets.For the same reason, taste degeneracy effects are less prominent where the spectral density is larger, such as in the bulk of the spectrum, as well as near the origin at temperatures below the pseudocritical temperature.On the other hand, the low end of the staggered Dirac spectrum at high temperatures shows a low spectral density, and taste-degeneracy effects on spectral correlations are strong. 2 The discussion above shows that if we take the thermodynamic limit before the continuum one, we can expect the distortions of the spectral statistics due to taste degeneracy to disappear, and the expected universal behavior based on the exact symmetries of the staggered operator should emerge.This would in turn allow one to employ standard methods to determine the localization properties of the eigenmodes using spectral statistics, that can then be extrapolated to the continuum.On the other hand, since the residual fermion doubling of sea staggered quarks is usually dealt with using the so-called "rooting trick" [50][51][52], one expects the correct order of limits to be continuum first and thermodynamic afterward.In this way, one recovers an exact taste degeneracy in the sea lattice action, so that taking the fourth root of the determinant in the path integral is justified.However, as explained above, the uncontrolled effects of taste degeneracy on the spectral statistics of the valence staggered operator make the determination of localization properties and of the position of the mobility edge using these quantities unreliable on finer lattices if the volume is not large enough.In this context, a non-trivial question is whether one can actually exchange the order of limits and still obtain the correct results.It is important to note that exchanging limits could be allowed 2 There is evidence [19,42,43] that also the staggered spectrum develops a near-zero peak of eigenvalues on sufficiently fine lattices, similarly to what is observed using the overlap operator in the valence [6,11,23,24,[44][45][46][47][48].In the latter case, near-zero peak modes also show peculiar localization properties [11,23,24,27].
The existence of such a peak in the Dirac spectrum in the presence of dynamical fermions is also supported by the arguments and the model calculations of Ref. [49].However, the mobility edge discussed so far in the literature and in this paper is found in a spectral region far above this peak, beyond the depleted region.
when determining the mobility edge even though it is not when studying spectral statistics.In fact, the mobility edge depends only on the localization properties of the eigenmodes, which are well defined entirely independently of the statistical properties of the spectrum.While localization properties and spectral statistics are related, it is the former that strongly affect the latter, and not the other way around.
Clearly, one could entirely bypass the problem by turning to the direct study of the eigenvectors themselves.This is, however, numerically more demanding, as it requires the use of several lattice volumes, larger statistics, and a more sophisticated finite-size-scaling analysis (see Refs. [53][54][55]) to achieve the same accuracy.Being able to reliably exploit the statistical properties of the spectrum is then desirable, and this requires the exchange of limits to be possible.In order for this to be also numerically efficient, one needs that the dependence of the mobility edge on the lattice spacing be mild, so that relatively coarse lattices, where the distortion due to taste degeneracy is negligible, would suffice for a continuum extrapolation.If this is the case, the results of Refs.[4,29], that were obtained for sufficiently large aspect ratios so as to avoid contamination from eigenvalue multiplets, but at the same time on numerically manageable volumes, could be fully trusted.
In this paper we study in detail the effects of taste degeneracy on the statistical properties of the staggered Dirac spectrum, and how they affect the determination of the mobility edge using spectral statistics, with the ultimate goal of providing a controlled investigation of the behavior of the mobility edge in the continuum limit.After briefly reviewing the properties of the staggered operator, and discussing localization and how to detect it using both eigenvector properties and the statistical properties of the spectrum, in Sec. 2 we discuss the complications due to the formation of nearly degenerate eigenvalue multiplets, and how the lattice spacing and the lattice volume affect them.We then argue that the continuum and thermodynamic limits can be exchanged when studying the mobility edge.In Sec. 3 we present numerical determinations of the mobility edge in the staggered spectrum of high-temperature QCD for various lattice spacings and volumes, obtained by a standard analysis of spectral statistics, and study the dependence on the lattice spacing and on the lattice volume.We then compare the results with a determination based on a direct study of eigenvector properties, unaffected by taste degeneracy, and with a heuristic estimate from the spectral statistics of a suitably "taste-symmetrized" spectrum, where the approximate taste degeneracy is removed by hand.Finally, we discuss the extrapolation to the continuum limit.We draw our conclusions in Sec. 4.

LOCALIZATION PROPERTIES OF STAGGERED EIGENMODES
This section is devoted to briefly summarizing the properties of the staggered discretization of the Dirac operator [56][57][58], as well as the main techniques and quantities used to study the localization properties of its eigenmodes.We then discuss in some detail the issues related to taste degeneracy.

A. The staggered lattice Dirac operator
The massless staggered operator in lattice QCD reads where U µ (x) ∈ SU(3), with x = (x 1 , x 2 , x 3 , x 4 ) and µ = 1, . . ., 4, is a gauge link variable attached to the link connecting the sites x and x+aμ of a N s ×N s ×N s ×N t hypercubic lattice with lattice spacing a; T µ are unit translation operators with periodic boundary conditions in the spatial directions µ = 1, 2, 3 and antiperiodic boundary conditions in the temporal direction µ = 4; and η µ are the staggered phases, η µ (x) = (−1) α<µ xα .The operator D stag is anti-Hermitean, so its eigenmodes obey with λ n ∈ R. For notational simplicity we will generally drop the dependence on the gauge configuration.These eigenmodes carry a spacetime index x, running over the lattice sites, and a colour index c = 1, 2, 3. Thanks to the chiral property {ε, D stag } = 0, where ε(x) = (−1) α xα , the spectrum of D stag is symmetric about zero, with The operator D stag is formally i times a random Hamiltonian, like those used in condensed matter physics to model systems with disorder (see, e.g., Refs.[59,60]).Here the disorder is provided by the fluctuations of the gauge fields, over which one integrates in the lattice formulation of gauge theories.The probability distribution of the entries of our random Hamiltonian is then determined by the specifics of the discretization of the Yang-Mills action for the gauge fields, as well as by those of the improvement techniques employed to speed up the approach to the continuum.In particular, the link variables U µ (x) corresponding to the discretized non-Abelian gauge fields need not be (and are not in current numerical practices) the same ones appearing in the lattice Yang-Mills action.These details are not relevant to the general discussion, and are presented below in Sec. 3 A. For our purposes it suffices to say that the integration over gauge fields defines an expectation value ⟨. ..⟩, corresponding to the disorder average in the language of disordered systems, with the usual properties ⟨1⟩ = 1 and ⟨O[U ]O[U ] * ⟩ ≥ 0 for a generic observable O that depends only on the gauge fields.

B. (Inverse) participation ratio and fractal dimension
The localization properties of the eigenmodes of the staggered Dirac operator are studied in full analogy with those of the eigenmodes of random Hamiltonians [59,60].In qualitative terms, the localization properties of the eigenmodes of such systems are defined by the scaling of their effective size, averaged over disorder realizations, with the size of the system.For localized modes, the average mode size remains constant in the large volume limit, while for delocalized modes it grows proportionally to the volume.
To make these statements quantitative, one introduces the inverse participation ratio (IPR) of the eigenmodes, defined here in a gauge-invariant way as where the latter sum is performed over color indices.It is understood that modes are normalized, x ∥ψ n (x)∥ 2 = 1.The effective mode size is simply IPR −1 : indeed, for a mode with constant amplitude square ∥ψ(x)∥ 2 = 1/V 0 in a region of size V 0 (in lattice units), one finds IPR −1 = V 0 , while for a fully delocalized mode with ∥ψ(x)∥ 2 = 1/(V s N t ), where V s = N 3 s is the spatial volume in lattice units, one finds IPR −1 = V s N t .Instead of the mode size, it is sometimes convenient to use the participation ratio (PR), i.e., the fraction of lattice volume effectively occupied by the mode: Notice that ψ −n and ψ n have the same IPR, so the localization properties are exactly the same for the positive and negative part of the spectrum.The localization properties of eigenmodes in a given spectral region are determined by the scaling of the PR in the thermodynamic limit, after averaging over gauge configurations (i.e., over disorder realizations), Here is the (non-normalized) spectral density, which is expected to scale proportionally to V s N t in the largevolume limit.In Eq. ( 6) we have made the dependence of PR on N s explicit, while leaving out that on N t , as the thermodynamic limit is taken here at fixed N t .In spectral regions where modes are localized, PR tends to zero as 1/V s as the spatial volume is increased, while for delocalized modes it tends to a constant.The large-volume behavior of PR determines the fractal dimension α of the modes (D 2 in the notation of Refs.[55,61]), defined locally in the spectrum as with α = 0 for localized modes and α = 3 for delocalized modes.Any intermediate behavior corresponds to "critical modes" in the language of condensed matter physics [61].
The discussion above deals with the lattice system at fixed N t .To study localization in the continuum limit at fixed temperature, one can either proceed as above, taking N s → ∞ at fixed a and T = 1/(aN t ); or take first a → 0 at fixed L = aN s and T to define a continuum system, and only after that study the scaling of the PR with L.

C. Spectral statistics
The localization properties of the eigenmodes of a random Hamiltonian determine the statistical properties of the corresponding eigenvalues [36].To see this, consider a specific disorder realization and the corresponding basis of eigenvectors, and perturb the disorder configuration by adding some local fluctuation.Once represented in the basis of the unperturbed eigenvectors, the perturbed Hamiltonian will have nonzero off-diagonal matrix elements for any pair of delocalized modes, while off-diagonal matrix elements involving localized modes can be non-negligibly different from zero only if they are localized where the fluctuation is introduced.In other words, delocalized modes are easily mixed by disorder fluctuations, with a dense matrix describing this mixing, while localized modes are sensitive only to disorder fluctuations appearing near their localization region.The statistical properties of eigenvalues associated with delocalized modes are then expected to match those of the appropriate Gaussian ensemble of Random Matrix Theory (RMT), according to the symmetry class of the disordered Hamiltonian [37,38].For uncorrelated disorder, or for disorder with only short-range correlations, localized modes are instead expected to fluctuate independently and thus to obey Poisson statistics.
A convenient observable in this context is the probability distribution of the so-called unfolded level spacings, i.e., the distance between neighboring eigenvalues measured on the scale of the average distance between neighboring eigenvalues found in the relevant spectral region [37,38].Formally, one defines the unfolded spec-trum via the mapping which, by construction, has unit spectral density anywhere, and the corresponding unfolded spacings as s n ≡ x n+1 − x n .When the volume becomes large, the probability distribution p(s; λ, N s ) of the unfolded spacings, computed locally in the spectrum, (10) is expected in general to depend only on the localization properties of the modes in the spectral region of interest, and on the symmetry class of the corresponding random Hamiltonian.For localized modes and Poisson statistics one expects an exponential distribution, while for delocalized modes one expects the distribution to be the same as the one found in the appropriate Gaussian ensemble of RMT.These are known but not available in closed form, and often approximated by the so-called "Wigner surmise", where β = 1, 2, 4 for the orthogonal, unitary, and symplectic symmetry class, respectively.The coefficients a β , b β are determined by the normalization conditions ds p(s) = 1 and ⟨s⟩ = ds p(s)s = 1, the latter following from the fact that the unfolded spectrum has unit spectral density.In the case of lattice QCD with staggered quarks, transforming under the fundamental representation of the gauge group SU(3), the relevant symmetry class is the unitary class. 3n exception to this classification is represented by the mobility edge separating localized and delocalized modes.At the mobility edge the localization length diverges, and the system undergoes a second-order phase transition in the spectrum, known as Anderson transition [61].Correspondingly, eigenmodes become critical, with a characteristic fractal dimension α * that depends on the symmetry class [55], and the corresponding eigenvalues obey a characteristic type of statistics p * (s), intermediate between Poisson and RMT-type [61][62][63].
The discussion above is valid only in the limit of infinite volume, where any overlap between different localized modes becomes entirely negligible, and where delocalized modes can fully spread out on a system of infinite size.In any finite volume there are instead corrections to the Poisson or RMT-type behavior, that tend to zero as the volume is increased.In the presence of a mobility edge one then expects p(s; λ, N s ) to interpolate between a near-exponential behavior and a near-RMT behavior, passing through the critical behavior at the mobility edge, where the properties of the eigenmodes are expected to be scale-invariant.Then, as N s is increased, p → p P in the localized spectral region, p → p RMT in the delocalized spectral region, and p = p * at the mobility edge λ = λ c .This can be used to determine λ c and the critical exponents of the Anderson transition by means of a finite size scaling analysis [63], as was done in Ref. [5] for QCD.Alternatively, if the critical distribution p * is known (as is the case, to some extent, for the unitary class [64]), one can use it to determine λ c as the point where the spectral statistics matches the critical one.This can be done very efficiently using a single lattice volume: by virtue of the scale-invariance of the critical point, such an estimate is expected to suffer only from very little finite-size effects.
As a concrete example, one can monitor how the integrated distribution I s0 , for a suitably chosen s 0 , varies along the spectrum, and where it matches the critical value.For the unitary class one customarily chooses s 0 ≃ 0.508 to maximize the difference between the expected values for Poisson and RMT-type statistics, i.e., I s0,P ≃ 0.398 and I s0,RMT ≃ 0.117; for this choice of s 0 the critical value I s0,c has been determined in Ref. [5].One could similarly use the second moment of the distribution, ⟨s 2 ⟩, or equivalently the variance Var(s) ≡ ⟨s 2 ⟩ − ⟨s⟩ 2 , or any other feature extracted from the unfolded level spacing distribution.For Var(s), the expectations for Poisson and for RMT statistics in the unitary class are Var(s) P = 1 and Var(s) RMT = 3π 8 − 1, respectively.The mobility edge can then be estimated as the point λ stat

D. Effects of taste degeneracy
The discussion above ignored entirely the possible presence of additional correlations between eigenvalues due to peculiarities in the structure of the random Hamiltonian.Such correlations can be induced by the presence of an approximate symmetry, leading to near-degenerate multiplets of eigenvalues.In this case, the corresponding symmetry-breaking parameter controls the typical splitting δ within the would-be-degenerate multiplets, which, in turn, gives the spectral scale at which their effects are felt.When this scale is smaller than or comparable to the expected distance ∆ = 1/ρ 0 between neighboring eigenvalues, it affects the unfolded level spacing distribution, skewing it toward lower values: Indeed, neighboring levels belonging to the same approximate multiplet prefer to stay closer than what one would expect simply based on the inverse of the spectral density, as their distance is controlled by a different mechanism with a corresponding smaller spectral scale.More precisely, one has clearly separated multiplets of n eigenvalues if their size (n − 1)δ is much smaller than the typical distance n/ρ 0 = n∆ between them, i.e., if δ ≪ ∆.In this case, out of n level spacings, n − 1 are of order δ, and one is of order n∆ − (n − 1)δ.This means that the average level spacing results from the averaging of fluctuations around two different typical values; for this reason, one expects also a larger variance for the unfolded level spacing distribution.
On the other hand, the spectral density is proportional to the volume, and so in the thermodynamic limit at fixed symmetry-breaking parameter many approximate multiplets will overlap.Indeed, when δ ∼ ∆, the size of a multiplet becomes comparable with the distance between multiplets, and their clear separation becomes impossible.Nonetheless, at this stage the eigenvectors are likely to still carry a remnant of the multiplet structure, that reflects on the correlations between the corresponding eigenvalues, including neighboring ones.When δ ≫ ∆, many multiplets overlap: the approximate symmetry will still affect the eigenvalue correlations, but the corresponding effects will be visible only on a spectral scale much larger than the typical separation between neighboring eigenvalues.Loosely speaking, assuming that an assignment of eigenvalues to multiplets still makes sense in such a dense environment, one finds that neighboring eigenvalues almost certainly belong to different multiplets.This means that short-range correlations involving eigenvalues at a fixed separation in mode number, such as the ones governing the unfolded level spacing distribution, are entirely unaffected by the approximate symmetry in the thermodynamic limit.Loosely speaking again, one would be measuring the distribution of the level spacing between members of different multiplets, for which there should be no distortion from the universal expectation.
The staggered operator provides precisely an example of the situation discussed above: while in the continuum limit it has an exact symmetry under the exchange of degenerate tastes, on a finite lattice this symmetry is only approximate.This manifests in the formation of nearlydegenerate eigenvalue multiplets as one gets sufficiently close to the continuum.The symmetry-breaking parameter is here the lattice spacing a, and the splittings δ in the nearly-degenerate eigenvalue multiplets are of order δ ∼ (aΛ) 2 Λ (see, e.g., Ref. [65,66]), with Λ some physical mass scale.In the thermodynamic limit at fixed a, the argument above applies: the multiplet structure gets washed out, and its effects on short-range eigenvalue correlations disappear.On the other hand, as a → 0 the splittings decrease and multiplets become more degenerate, and more distinguishable from each other: this effect competes with the wash-out effect of the thermodynamic limit, and so the order in which the continuum and the thermodynamic limits are taken becomes important.
In finite-temperature calculations, it is customary to take the continuum limit a → 0 at fixed temperature T = 1/(aN t ) and fixed aspect ratio LT = N s /N t , corresponding to keeping the lattice spatial volume (as well as its temporal extension) fixed in physical units.In this case, the relative size of the multiplet splittings δ and the typical level spacing ∆ i.e., it is expected to decrease quadratically with a.This means that short-range correlations become more sensitive to the effects of taste degeneracy as a → 0, leading to significant systematic effects in the study of spectral statistics, and sizeable deviations from the expected universal behavior of the unfolded spectrum, as soon as δ ≲ ∆.The distortion of the unfolded level spacing distribution caused by the emergence of nearly-degenerate eigenvalue multiplets leads to values of I s0 and Var(s) larger than one would expect based on the localization properties of the eigenmodes, and an estimate λ stat c of the mobility edge using the intercept with the critical value of some spectral statistic, as explained above, would generally lead to an overestimation of the actual value λ true c .The actual physical value of the mobility edge is, of course, independent of the use of spectral statistics to determine it, as it originates in the properties of the eigenvectors.In this context, its most important aspect is its characterization as the point in the spectrum where the localization length diverges and a second-order Anderson transition takes place.This allows one to determine it, for example, by using the intercept of the local fractal dimension α(λ) with the corresponding critical value α * .
For this reason, the determination of λ c from spectral statistics is expected to converge to the true position of the mobility edge if one takes the thermodynamic limit at fixed lattice spacing.In fact, as explained above, by taking limits in this order one gets rid of the taste-multiplet problem, and the spectral statistics are expected to behave like those of an ordinary system without approx-imate symmetries.Then, the scale-invariant nature of the physics of critical Dirac eigenmodes should reflect in the statistical properties of the spectrum in the usual way.This means that the spectral statistics at the mobility edge should be independent of the volume for sufficiently large lattice sizes (as soon as distortions due to taste symmetries become negligible), and agree with the universal critical statistics expected for the given symmetry class (irrespectively of approximate symmetries of the system).Conversely, a scale-invariant behavior of the spectral statistics is a sign of a second-order transition in the spectrum, and so of the presence of a bona fide mobility edge.In other words, one expects that at fixed N t , shows a plateau in N s , entirely unaffected by taste-degeneracy effects.This, of course, should be verified explicitly.
However, even if this is the case, this approach rapidly gets numerically very intensive both for the generation of gauge configurations and for the diagonalization of the staggered operator as one gets closer to the continuum, until it becomes impossible to reach sufficiently large aspect ratios to see the plateau, and the extrapolation of the results for the mobility edge toward the thermodynamic limit introduces larger systematic effects.As the lattice spacing is reduced the statement that the approach based on spectral statistics is numerically more efficient becomes questionable.
Clearly, the ideal scenario would be that indeed λ stat c (T ; N t , N s ) tends to the true mobility edge as the volume is increased; that the continuum and thermodynamic limits commute; and that the true mobility edge in units of the quark mass, λ true  (T ; N t )/m s , so that relatively coarse lattices, for which distortion effects due to taste degeneracy are negligible, can be used for a continuum extrapolation.

E. Continuum and thermodynamic limits
While the discussion above shows that continuum and thermodynamic limits do not commute for spectral statistics (unless possibly if one takes care of the approximate taste degeneracy), it does not exclude this possibility for the mobility edge.In fact, the mobility edge is directly related to the localization properties of the Dirac eigenmodes, and only indirectly to spectral statistics, through the influence that localization properties have on them; one could identify it by looking only at eigenvector properties, without any reference to spectral statistics.We now argue that the two limits can indeed be exchanged when studying localization properties and the position of the mobility edge.
If one switched off the taste-breaking part of the staggered operator, one would find degenerate eigenspaces of "unperturbed" modes related by taste symmetry trans-formations.Since these act on the scale of an elementary hypercube, mixing the corresponding components of the quark wave function, the modes in a degenerate eigenspace will have similar localization properties; in particular, all the modes in the eigenspace corresponding to a localized mode will be localized in the same region and will have similar sizes.In the symmetric theory, one could then straightforwardly define an "unperturbed" mobility edge.The question is how this changes when taste-breaking effects are taken into account.
To answer this question, we set H = −iD stag /Λ = H 0 + g∆, where H 0 is the taste-symmetric part, ∆ the taste-breaking part, Λ a typical QCD scale, and g = (aΛ) 2 plays the role of dimensionless coupling.In dimensionless units, we have E n = λ n /Λ for the eigenvalues of H, and L Λ = LΛ and V Λ = (LΛ) 3 = V Λ 3 for the system size.To find the exact spectrum {E n } and eigenvectors {ψ n } of H, we work in the basis {ψ In this basis we have to diagonalize the matrix For definiteness, within degenerate subspaces we choose eigenvectors that diagonalize ∆ restricted to the subspace, i.e., ∆ n .We work with fine and large lattices, so that the scale of the taste-breaking interaction, g, and that of the inverse spectral density, 1/V Λ , are both small, but we do not initially specify which one is smaller.As discussed in the previous subsection, these are the two scales relevant to the status of taste-symmetry effects in the exact spectrum, and their relative size distinguishes two regimes: if g < 1/V Λ , quartets are well-formed, with O(g) intraquartet splitting and O(1/V Λ ) distance between quartets; if 1/V Λ < g, quartets are overlapping and level spacings are generally O(1/V Λ ).The level spacing between nondegenerate unperturbed eigenvalues is, in any case, of order O(1/V Λ ).
The unperturbed eigenvectors {ψ n } are in part localized and in part delocalized, but working in a finite volume the distinction is, of course, ambiguous.One can nonetheless provide a finite-volume estimate of the unperturbed mobility edge E (0) c using, e.g., the fractal dimension, or even spectral statistics, that are unproblematic here since one can remove the exact degeneracy by hand.Using standard scaling arguments (see, e.g., Ref. [5]), the difference between such estimates and the infinite-volume value is expected to be not more than of the order O(1/L 1/ν Λ ), where ν is the correlation length critical exponent of the Anderson transition [59][60][61], and in practice much less than that.In fact, the localized or delocalized nature of modes in a finite volume is ambiguous at most within the scaling region 1) around the mobility edge, and truly so only within the much smaller region where corrections to scaling are important, and one cannot decide the nature of the modes based on the leading behavior.
Switching the taste-breaking interactions back on, the unperturbed eigenvectors will mix, and localization properties may change.To see to what extent this is the case, we need to estimate the order of magnitude of the matrix elements ∆ n ′ n .For n and n ′ both localized, ∆ n ′ n receives contributions only from the spacetime region where the two localized modes overlap.Since localization is typically exponential one expects with ξ n , ξ n ′ the localization lengths of the two modes, and d n ′ n the distance between their localization centers. 4The largest matrix elements correspond to pairs of modes whose localization regions are close, but whose eigenvalues are typically far away from each other in the spectrum.Conversely, eigenvalues that are close in the spectrum correspond typically to modes that are far away in space: since on average they are distributed uniformly in space, the chance that they are nearby is , where N loc is the number of localized modes.
For n localized and n ′ delocalized, ∆ n ′ n receives contributions only from the localization region of ψ (0) n , and since ψ (0) c , where localized modes have size ξ 3 n = O(V ), this can be as large as O(1); however, H n ′ n is further suppressed by a factor of g.Finally, for n and n ′ both delocalized, ∆ n ′ n receives contributions from the whole volume, and is of order In the delocalized part of the unperturbed spectrum, the effect of the taste-breaking interaction is expected to wash out the quartet structure, despite the smallness of g, if 1/V Λ < g: in this case the typical unperturbed level spacing is smaller than the effect of the perturbation, and mixing affects members of different quartets.The tastebreaking interaction is otherwise expected not to change much neither the localization properties of the modes, nor to shift much the spectrum.In fact, this interaction amounts in practice to adding some weak off-diagonal disorder to the system, as the perturbation ∆ has only off-diagonal matrix elements of order 1 (as they are built out of the unitary gauge link variables), and their order of magnitude does not fluctuate across the lattice.This type of disorder is known to be largely ineffective in localizing modes except in the tails of the spectrum, see Refs.[67][68][69].Here it is the size of the fluctuations of ∆ that is relevant, and not the relative size of g and 1/V Λ , so taste symmetry plays no role.One then expects that if the mobility edge moves, it does not do so because of interference effects between strongly mixing delocalized modes leading to localization, but because of mixing between localized and delocalized modes, or between localized modes alone leading to delocalization.
Indeed, using the obvious identity and denoting by Π the projectors on the localized or delocalized unperturbed eigenvectors, one has for any ψ n the exact bounds where minimization is restricted to the localized unperturbed spectrum, and where minimization is restricted to the delocalized unperturbed spectrum.
Due to the weak localizing power of the taste-breaking interaction, one expects that if ψ n is a localized exact mode, then it does not come entirely from the mixing of delocalized unperturbed modes, so that 1 − ∥Π (0) deloc ψ n ∥ 2 is a nonzero number of order O(1) (except possibly in the tails of the spectrum).One finds then from Eq. ( 20) that for a localized exact mode, E n < E (0) c + δ + , with δ + a positive number of order O(g).
On the other hand, Eq. (21) shows that if ψ n is a delocalized exact mode, then it can be found below the unperturbed mobility edge at distance of order O(g , if a delocalized ψ n is obtained essentially only by forming a linear combination of localized unperturbed modes, up to a contribution from delocalized unperturbed modes of order O(g ϵ ), i.e., (0) loc,deloc linear combinations of localized or delocalized unperturbed modes only, with norm of order 1. (If the contribution of delocalized unperturbed modes comes with a higher power of g, the bound is ineffective.)This seems unlikely: in fact, matrix elements H n ′ n are of order O(g) • O(1) only for pairs of localized modes with close or overlapping spatial support, which on the other hand are typically well separated in the spectrum.One then expects their mixing to be negligible when computing the exact eigenmodes, and so one expects to find localized exact modes only, if only localized unperturbed modes are involved.
To see that this is the case, we impose that ψ n is an exact eigenvector with eigenvalue E, and projecting on the localized unperturbed mode ψ (0) n , we obtain the equation where loc ) and r n = (ψ deloc ).Since g is small, we can attempt to solve it perturbatively, checking afterward that the resulting mixing coefficients are small.Setting E = E 0 + tE 1 + . . .and c n = c n + . .., and replacing g∆ → tg∆ and r n → tr n , with t a formal expansion parameter that will be eventually set to 1, one finds yielding to lowest orders in t the equations These coefficients are generally if the unperturbed eigenvalues are not close, and only for order O(1) modes (i.e., those with small d n ′ n /ξ n ′ n ) one finds a numerator of order O(1).If the unperturbed eigenvalues are nearby in the spectrum, i.e., they differ by O(1/V Λ ), then there is only a O(1/V Λ ) chance that also the corresponding eigenvectors are localized near each other in space, and so for the typical coefficient one has the estimate 5 independently of the volume, which is small, as it should be for perturbation theory to be valid, and negligible compared to c ñ ≃ c ñ but n ̸ = ñ are not determined at this stage.To do that we write Eq. ( 24), which in this case starts from O(t 2 ), and reads having assumed full lifting of the degeneracy (i.e., and n ̸ = ñ).For ϵ < 1, the 5 For a more accurate estimate when g > 1/V Λ , spatially close modes with exactly or nearly degenerate (i.e., E n − E (0) eigenvalues should be treated together as degenerate modes, putting any small difference into the perturbation.In this case the magnitude of the perturbation (after including in it the small eigenvalue differences) remains of order O(g).The perturbation should then be diagonalized exactly in the corresponding subspace, which has dimension O(1), since the probability of finding k modes nearby both in space and in the spectrum is suppressed as 1/V k Λ .The mixing of nearby localized modes leads, of course, to similarly localized modes.For modes in different subspaces, which are now spectrally well-separated, the perturbative analysis discussed above works without any issue.Further mixing of modes within a degenerate subspace cannot change their localization properties, see also discussion below.
second term dominates, but it still gives a small contribution when g ≪ 1, so perturbation theory can be trusted. 6(The fact that the solution is not analytic in g is irrelevant.)Furthermore, mixing between degenerate modes, no matter how large, cannot lead to delocalization, since these modes have the same localization properties and similar spatial support.To fully determine the eigenvector ψ n one would need to determine Ψ (0) deloc selfconsistently, but this is not relevant for our goal: for us, it is enough having shown that a delocalized mode of the required type cannot be found, no matter what Ψ (0) deloc is.While one generally expects perturbation theory to break down in a large volume, in this case it does not.The reason is that the number of localized modes with non-negligible spatial overlap remains always of order O(1), no matter how large the volume is; and that even though the level spacing becomes of order O(1/V Λ ), neighboring modes in the spectrum have only a chance of order O(1/V Λ ) to be also close in space (and so to have non-negligible overlap).The two factors essentially cancel out, and mixing still involves O(1) modes with a factor of order g, which is small in our setup.In other words, mixing of an unperturbed localized mode with different unperturbed localized modes will make up a tiny fraction of the perturbed mode's weight, and so cannot change its localized nature.
In conclusion, no delocalized mode as in Eq. ( 22) can be built out of localized modes alone, showing that if the exact mode ψ n is delocalized, then E n > E (0) c − δ − , with δ − a positive number of order O(g ϵ ), for any 0 < ϵ ≤ 1. (Notice that the case ϵ = 0 is excluded, so that delocalized exact modes including an O(1) contribution from delocalized unperturbed modes can exist.These, however, cannot be found farther away than O(g) from the unperturbed mobility edge.)Together with the upper bound found above, this leads us to conclude that for the exact mobility edge E c one has All in all, the possible effect of taste-breaking interactions on the mobility edge (in a finite volume) is to shift its unperturbed value by at most O(g).As these are short-ranged, UV effects, there are no further volume dependencies expected.Besides taste-breaking effects, other finite-spacing effects amount only to renormalizing the whole spectrum like a quark mass [70][71][72][73][74].
The argument above shows that when g < 1/V Λ , changes in the eigenvalues of would-be degenerate localized modes induced by introducing a localized fluctuation in the gauge field will be strongly correlated, as they are mainly driven by the O(1/V Λ ) change in the corresponding unperturbed eigenvalue.This gives a microscopic explanation for the change in local spectral statistics in the presence of an approximate taste symmetry, discussed in the previous subsection.On the other hand, while the approximate taste symmetry strongly affects the determination of the mobility edge from spectral statistics, it does not affect much the localization properties of modes and the true position of the mobility edge, even in this regime.
As already mentioned above, on top of finite-spacing effects there are finite-size effects that can shift the finite-volume estimate of the mobility edge by at most O(1/L 1/ν Λ ).This can occur because of mixing of localized and delocalized modes around the mobility edge, where the size of localized modes is of the order of the volume of the system.These are IR effects controlled by the ratio of the localization/correlation length of eigenmodes and the linear size of the system, which is a ratio of physical scales and, as such, a renormalization-group-invariant quantity, that will become independent of the spacing as this goes to zero.That the correlation length of eigenvectors is physical follows from the fact that correlators such as ⟨ n δ(λ − λ n )∥ψ n (x)∥ 2 ∥ψ n (0)∥ 2 ⟩ renormalize multiplicatively, up to renormalizing also the spectrum [35].
In conclusion, returning to physical units and renormalizing the spectrum, for the exact and unperturbed mobility edges in quark mass units, λc = λ c /m s = ΛE c /m s and λ(0 where for the behavior of finite-size corrections we conservatively used its expected upper bound.Notice that taking the thermodynamic limit first one finds lim L→∞ λc − λ(0) c = O (aΛ) 2 .Since λ(0) c is by construction not affected by taste-breaking effects, one expects that continuum and thermodynamic limits commute for it, i.e., where corrections, vanishing in the continuum and thermodynamic limit taken in any order, are expected to be at most O(1/L 1/ν ) in the system size, and likely smaller, and O (aΛ) 2 in the lattice spacing thanks to the lattice chiral symmetry of staggered fermions.Then Eq. (31) shows that the two limits commute also for λ c /m s .

NUMERICAL RESULTS
In this section, after briefly summarizing our lattice setup, we present our results for the mobility edge, discussing in detail the effects of taste degeneracy on the spectral statistics.

A. Lattice setup
In this work we considered 5 ensembles of gauge configurations of N f = 2 + 1 QCD at the physical point, originally generated for the investigation of the QCD topological susceptibility at finite temperature reported in Ref. [75], obtained on N 3 s × N t lattices with lattice spacings a ranging from ∼ 0.107 fm down to ∼ 0.054 fm.In all cases, the temperature was fixed to T = 1/(aN t ) ≃ 230 MeV, and the lattice size to L = aN s ≃ 3.4 fm, corresponding to an aspect ratio N s /N t = LT = 4.For the next-to-finest lattice spacing, corresponding to N t = 14, we also considered 2 additional ensembles with different lattice volumes, with aspect ratios LT ≃ 3.4 and ≃ 4.6, corresponding respectively to L ≃ 2.9 fm and ≃ 3.9 fm.In order to keep our paper self-contained, we briefly summarize below the lattice setup employed in Ref. [75] to generate configurations, and we refer the reader to the original reference for further technical details.
Gauge configurations were generated adopting the tree-level Symanzik improved Wilson gauge action [76][77][78][79] to discretize the pure-Yang-Mills term, and rooted stout-smeared staggered fermions to discretize the Dirac term.Expectation values are then defined schematically as Here S g denotes the discretized pure-Yang-Mills term, whose exact form is not relevant; U µ (x) ∈ SU( 3) is the gauge link variable attached to the link connecting the sites x and x+aμ, and [dU ] = x,µ dU µ (x) is the product of the corresponding Haar integration measures.Moreover, , M [U (2) ; m] = det D stag [U (2) ] + m which is real positive thanks to the symmetry of the spectrum [see after Eq. ( 2)].Here U µ (x) ∈ SU(3) is obtained from U µ (x) after 2 steps of isotropic stout smearing [80] with smearing parameter ρ = 0.15.The bare (inverse) coupling β, the bare light and strange quark masses m l and m s , and the lattice spacing a were fixed in Ref. [75] according to the results of Refs.[81][82][83] so as to stay at the physical point, defined by a Line of Constant Physics with the pion mass m π = 135 MeV and the strange-tolight quark mass ratio R = m s /m l = 28.15fixed to their physical values.All simulation parameters are summarized in Tab.I.
The Monte Carlo simulations of Ref. [75] employed the Rational Hybrid Monte Carlo (RHMC) updating algorithm [84][85][86], used, for some simulation points, in combination with the multicanonical algorithm [75,[87][88][89][90][91].As a matter of fact, it is well known that, above the QCD chiral crossover T c ≃ 155 MeV, the topological susceptibility χ = ⟨Q 2 ⟩/V 4 , where V 4 = L 3 /T , is rapidly suppressed as a function of the temperature [75,[92][93][94][95][96][97], leading to ⟨Q 2 ⟩ = V 4 χ ≪ 1 for the lattices employed in  typical simulations.Multicanonical simulations are employed to enhance the probability of visiting nontrivial topological sectors when this probability is suppressed, without spoiling importance sampling.For what concerns the present investigation, the multicanonical algorithm allows one to avoid possible systematic effects in the low-lying staggered spectrum, which is connected to topological excitations.Concerning the technical details of multicanonical runs, we refer the reader to the original work [75].
For each gauge configuration we then analyzed the lowest 150 positive eigenvalues of D stag [U (2) ] [i.e., the same operator appearing in the fermionic determinant, Eq. ( 34)], obtained in Ref. [75] using the PARPACK library [98], from which we computed the corresponding unfolded eigenvalues.We complemented these data by performing a new calculation to obtain the IPRs of the Dirac eigenvectors for the gauge ensembles with N t = 14.We calculated error bars from a standard binned jackknife performed over the whole analysis, including the unfolding procedure (when studying spectral statistics), to take into account possible correlations introduced by the latter.An example of a staggered Dirac spectrum on a single, typical configuration is shown in Fig. 1.We performed unfolding by ordering all the available eigenvalues λ n of all the configurations of a given lattice ensemble, and replacing them with their rank divided by the number of configurations [17].This automatically yields unit spectral density, while the deviation of ⟨s⟩ from 1 can be used to test the accuracy of the procedure.
To compute I s0 and Var(s) locally in the spectrum we approximated Eq. ( 10) by binning the spectrum in bins of size 0.178m s in physical units, averaging inside each bin, and assigning the result to the center of the bin.Since, loosely speaking, the Dirac spectrum renormalizes like the quark mass [70][71][72][73][74], keeping the bin size fixed in units of m s is the appropriate choice.An unfolded spacing was included in a bin average if the corresponding lowest (not unfolded) eigenvalue fell in the bin.
We similarly computed ⟨s⟩, and found it compatible with 1 within errors in the relevant spectral region λ/m s ≥ 1.5, see Fig. 2.This reassures us on the validity of our unfolding procedure.Between 1 and 1.5 there is a small but significant deviation from 1, systematically increasing as one goes down in λ.This is well understood [20], and it is due simply to the low but rapidly changing spectral density, that requires the use of large bins to have a decent signal, at the price of having a non-constant density inside the bin.This leads to the smaller eigenvalue spacings corresponding to modes at the higher, and denser, end of the bin, lowering the numerical estimate ⟨∆λ⟩ bin of the local average spacing between neighboring eigenvalues in a spectral bin below the expected value 1/ρ 0 bin , leading to ⟨s⟩ bin ≃ ⟨∆λ⟩ bin ρ 0 bin < 1 for the numerical estimate of ⟨s⟩.Between 1 and 1.5 one can also see a systematically increasing deviation from 1 as N t increases (see Fig. 2, top panel).This is due to the formation of distinct tastedegenerate multiplets in this spectral region which, being most likely found at the higher end of a spectral bin, are also more likely to spread across neighboring bins.This leads to lowering ⟨∆λ⟩ bin further, as spacings within a taste multiplet are more likely to contribute than spacings between multiplets.This explanation is supported by the fact that increasing N s at fixed N t , so bringing the multiplets closer, reduces the distance of ⟨s⟩ bin from 1, see Fig. 2, bottom panel.Finally, for λ/m s < 1 the spectrum is very sparse, and it is hard to make any reliable statement.

B. Effects of taste degeneracy on mode countings
As a preliminary piece of evidence of the strong effects of taste degeneracy on the statistical properties of the spectrum, in Fig. 3 we show with red circles the probability distribution of the number ν(M/m s ) of eigenmodes found below a fixed cutoff, i.e., λ n ≤ M , with M cho- sen safely below the mobility edge (see below Sec. 3 F).Notice that ν is a renormalized quantity if the renormalized physical value of M is kept fixed [70][71][72][73][74], e.g., if M/m s is kept fixed.While the expectation is that this counting follows a Poisson distribution with parameter equal to the average number of modes below the cutoff, ⟨ν(M/m s )⟩, the data show otherwise.
Quite striking at first sight is the oscillating behavior that leaves the odd bins almost empty.This, however, is easily understood if we take notice that the low-lying spectrum shows the formation of eigenvalue doublets (see Fig. 1), the first step in the formation of the quartets expected in the continuum limit (see, e.g., Ref. [99]).If we correct for this by counting the number of doublets, ν/2, instead of the number of modes, and compare with the Poisson distribution with parameter equal to the average number of doublets (i.e., half the average number of modes) below the cutoff, then we find excellent agreement between the two curves, shown with blue squares and yellow triangles in Fig. 3.This is a simple but convincing demonstration that spectral statistics are heavily distorted by tastedegeneracy effects on fine lattices.It also shows that the taste doublets of eigenmodes fluctuate independently, as expected for localized modes in the absence of neardegeneracy.We now proceed with a quantitative assessment of these effects on the unfolded level spacing distribution.

C. Determination of the mobility edge from spectral statistics
We estimated the mobility edge using spectral statistics as the point where these match their critical behavior using two features of the unfolded level spacing distribution, namely I s0 and Var(s).For I s0 , the critical value was obtained in Ref. [5] by means of a finite-size-scaling analysis, and reads I s0,c = 0.1966 (25) . ( Here we determined the critical value of Var(s) by means of a similar analysis on the same lattice data (2 + 1 QCD at T ≃ 2.6T c , a ≃ 0.125 fm, physical m l,s ):7 technical details about the procedure can be found in Ref. [5].Our estimate is Var(s) c = 0.3702 (98) .The central value is obtained using only lattices of spatial size N s ≥ 40 and data points in a range of width a∆λ = 0.026 around the mobility edge.The error includes the contribution of the statistical error from the fit, performed with the MINUIT routine [100]; the systematic finite-volume effect estimated as the change of the fit result due to including also N s = 36 data; and the systematic of the fitting range estimated as the change of the fit result due to shifting the fitting range down in the spectrum by 10% of its width. 8 The qualitative behavior of the dependence of λ stat c /m s on the lattice spacing at fixed aspect ratio can be easily deduced from Fig. 4: as a → 0, both I s0 and Var(s) tend to increase throughout the spectrum, as expected, 8 We also find aλ (T =2.6Tc) c = 0.33602 (63) for the mobility edge and ν = 1.406 (98) for the localization-length critical exponent, in agreement with those found in Ref. [5].leading to a larger estimate for the mobility edge.In particular, I s0 overshoots the Poisson expectation at the low end of the spectrum already on our coarsest lattices (even on the 32 3 ×8, where this happens outside of the spectral window displayed in the plot), showing that taste degeneracy is already distorting the spectral statistics there.While this is not a problem as it does not affect the region where the mobility edge actually is (see below Sec. 3 F), as the lattice becomes finer I s0 also overshoots its RMT expectation deeper in the bulk of the spectrum, signaling that the effects of taste multiplets on spectral statistics become important there, too.
The data show that the two coarsest lattices (N t = 8, 10) give compatible results for the spectral statistics in all spectral regions in the bulk of the spectrum, and down to the first bin below the mobility edge.Barring a conspiracy, this indicates two things: that the formation of taste multiplets does not have significantly large effects in that spectral region yet; and that, when this is the case, further lattice artifacts, on top of those in- troduced by taste-degenerate multiplets, are small when considering localization properties at fixed physical volume.Clearly, there could still be finite-volume systematics, but these would involve only effects unrelated to taste multiplets (such as the localization length in the spectral region near the mobility edge being too large compared with the spatial size of the system), which are expected not to affect much the determination of the mobility edge via the matching with the critical statistics (see above in Sec. 2 C).Indeed, previous numerical results show that for lattices of similar size and spacing one is already close to the one-parameter scaling regime near the mobility edge [5].This shows that the determination of the mobility edge from the critical spectral statistics is reliable on these lattices.
The dependence of λ stat c /m s on the aspect ratio at fixed lattice spacing is visible instead in Fig. 5, where TABLE II: Summary of our results for the mobility edge in units of the strange quark mass, obtained by matching the behavior of the spectral statistics with the expected critical behavior.
N t = 14 and the spatial size is varied.Here a larger aspect ratio drives I s0 and Var(s) down, as expected, so leading to an estimate for the mobility edge that decreases with the spatial volume.Notice that the scaleinvariant nature of the mobility edge is still masked here by the distortions of the unfolded level spacing distribution, and no volume-independent crossing point of the various curves is present.This shows that the effects of taste degeneracy are still strong even on our largest lattice at N t = 14.This is in agreement with the argument discussed in Sec. 2 D, as, according to Eq. ( 15), we would need an aspect ratio of LT ∼ 9 at this lattice spacing in order to have the same δ/∆ that we have for N t = 8 with LT = 4.We then expect that the plateau in λ stat c /m s against LT is not visible yet, and would appear for an aspect ratio roughly twice as big as our largest one.
In order to estimate λ stat c /m s quantitatively, we interpolated our numerical results for I s0 and Var(s) with splines, defining an uncertainty band by interpolating the central values augmented or reduced by the error.We then determined λ stat c /m s as the center of the interval where the band crosses the critical value, and the corresponding error as the half-width of this interval.This procedure is visualized in Fig. 6 for our 32 3 × 8 lattice.Our final results are summarized in Tab.II.D. Taste-degeneracy effects on the mobility edge in the thermodynamic limit In order to show that the effects of taste degeneracy become irrelevant for the determination of the mobility edge with our method when taking the thermodynamic limit at fixed lattice spacing, we have compared the N s → ∞ extrapolation of our N t = 14 data with an independent determination of λ c based directly on the properties of the eigenvectors.Distortive effects on the statistics are due to approximate taste degeneracy, which is visible when (aΛ) 2 < 1/(V Λ 3 ).We then expect finite-size corrections to λ stat c /m s of order O(1/N /ms assuming leading linear corrections in the inverse spatial volume (Nt/Ns) 3 = 1/(LT ) 3 .The filled points at zero are the results of our extrapolations.Empty and filled starred points, and the shaded area, are our estimates of the mobility edge from the fractal dimension, see text for more details.Points have been slightly shifted horizontally to improve readability.
appear.As these are expected to be about an order of magnitude larger than our largest lattices, we can ignore the onset of the plateau in our extrapolation in (N t /N s ) 3 .The infinite-volume extrapolation is shown in Fig. 7; the corresponding results for the mobility edge are: To determine λ c directly from the eigenvectors, we estimated the local fractal dimension numerically as Our results are shown in Fig. 8.We then looked for the point in the spectrum where α num takes its critical value α * = 1.173 +32 −26 [55], using a spline interpolation of the numerical data.This is visualized in Fig. 9.To quote a final value, we considered the results obtained with the pairs (N s1 , N s2 ) = (48,56) and (48,64), and took as our final estimate a symmetric confidence band including both error bars.The two values are compatible within errors (see Fig. 7), showing that finite-size effects on this estimate of the mobility edge are reasonably small.As a further check, we verified that the pair (N s1 , N s2 ) = (56,64) gave compatible results, although within much  larger statistical errors.In the end, we obtain: 75), (from fractal dimension).(40) This is in good agreement with the thermodynamic extrapolation of our estimates from spectral statistics, Eq. ( 37), see Fig. 7.This shows that these are indeed converging to the actual position of the mobility edge, as we argued in Sec. 2 D.Moreover, this further confirms that finite-size effects on the estimate based on the critical fractal dimension are reasonably small.We finally summarize the results of the previous two subsections in Fig. 10, showing together all our determinations of λ c /m s .It is clear that our two coarsest lattices give compatible estimates via λ stat c /m s , while on finer lattices with the same aspect ratio these estimates rapidly deviate.On the other hand, the infinite-volume extrapolation at N t = 14 gives again results compatible with our coarsest lattices, as well as with the determination from the local fractal dimension.
Before attempting an extrapolation of λ c /m s to the continuum, it is important to check that our numerical results are in agreement with our theoretical argument developed in Sec. 2 E that the thermodynamic and continuum limits can be exchanged.While a direct check is not possible with the available data, indirect support is provided by our results for the estimate based on α(λ).Being unaffected by taste degeneracy, this gives us a valid estimate of the position of the mobility edge at both finite aspect ratio and finite spacing.We have already pointed out above that this direct estimate of the mobility edge does not depend strongly on the aspect ratios employed, and agrees with the indirect estimate from spectral statistics extrapolated to infinite volume.Moreover, our results obtained at N t = 14 are not very different from the estimates obtained using spectral statistics on our coarsest lattices N t = 8, 10, where taste-degeneracy effects on spectral statistics are under control and estimates based on critical statistics should accurately reflect the posi-tion of the mobility edge (within our numerical accuracy).These observations suggest that the true position of the mobility edge, as intrinsically defined by the properties of the eigenvectors, does not depend strongly neither on the aspect ratios nor on the lattice spacing.This is in agreement with our argument that the continuum and thermodynamic limit commute for what concerns the mobility edge.We can then use indirect estimates of the mobility edge from spectral statistics extrapolated to the thermodynamic limit, followed by the continuum limit, to obtain the correct value of the mobility edge in the continuum theory.
Further numerical support to the correctness of our theoretical argument can be obtained if we manage to remove the effects of taste degeneracy from the spectrum.For a suitably "taste-symmetrized" spectrum, one expects a more accurate correspondence between localization properties of the eigenvectors and spectral statistics already at finite volume.Agreement of the resulting estimates for the mobility edge with the ones obtained by using the spectrum without further processing and taking the thermodynamic limit, would provide further support to the correctness of our procedure.Clearly, the taste-symmetrization of the spectrum is an ambiguous procedure, even if one could correctly identify the eigenvalue multiplets, which is not easy when the spacing is too coarse; even on our finest lattices, only doublets are visible.This ambiguity leads to hardly quantifiable systematic effects, that make this approach less reliable for a controlled continuum extrapolation.
To obtain a heuristic estimate we have then opted for the simplest approach, namely replacing pairs of neighboring eigenvalues with their arithmetic average.This has a number of obvious theoretical issues, first of all the likely misidentification of at least part of the doublets.This problem is particularly severe in the bulk of the spectrum, where doublet misidentification is almost certain, leading to interfering with the relevant spectral correlations and possibly to destroying them completely.Nonetheless, one expects that for finer and finer lattices at fixed aspect ratio this procedure should work better and better (and worse and worse on bigger and bigger lattices at fixed lattice spacing), and eventually provide the correct result in the continuum limit (if one extends it to quartets when these show up); and one can hope that for sufficiently fine but still manageable lattices it works reasonably well up to and slightly above the mobility edge.We have then averaged the spectrum as discussed above, and performed on it the same analysis as with the untreated spectrum, looking for the point where spectral statistics match their critical value.We restricted the spectrum to regions where the unfolded symmetrized level spacings satisfy ⟨s⟩ ≃ 1, as a minimum requirement for the reliability of the procedure.This is illustrated in Fig. 11.
With this long list of caveats in mind, we show our results in Figs.11 and 12.As expected, RMT correlations in the bulk of the spectrum are lost.However, in spite  of our rather crude approach, the agreement between estimates based on the taste-averaged and untreated spectra is surprisingly good already for our N t = 14 lattices.While uncertainties are hard to quantify, we believe that these results provide strong numerical support to our theoretical argument for the exchangeability of continuum and thermodynamic limits.
For future utility, we report also the results obtained with the taste-averaging procedure on our finest, N t = 16 lattice, i.e., λ c /m s = 1.754 (33) from I s0 and λ c /m s = 1.766(54) from Var(s).

F. Continuum extrapolation
In the light of this discussion, we have performed an extrapolation to the continuum limit using our estimates of λ c /m s from spectral statistics for N t = 8, 10 and N t = 14, where for N t = 14 we have used the results extrapolated first toward the thermodynamic limit.For the continuum extrapolation we have assumed the behavior λc ms = λc ms cont + C/N 2 t , based on our theoretical expectations discussed in Sec.12: Estimates for the mobility edge from the spectral statistics of the spectrum after taste averaging (see text for details), obtained using results from our Nt = 14, 16 lattices.
For comparison, we show also our results from the fractal dimension and from standard spectral statistics extrapolated to the thermodynamic limit on Nt = 14 lattices, and our continuum extrapolation based on standard spectral statistics (see Sec. 3 F).
from the two observables I s0 and Var(s) are: λ c m s cont = 1.748(75), (from I s0 ), (41) λ c m s cont = 1.816 (43), (from Var(s)), (42) perfectly compatible within errors.We note that, while affected by uncontrolled theoretical uncertainties, the results obtained with our taste-averaging procedure on our finest, N t = 16 lattice, reported above at the end of Sec. 3 E, are in excellent agreement with these estimates (see Fig. 12).
To quote a final number, we did a weighted average of our two estimates and their statistical errors, and used their difference as an estimate of systematic finite-size effects.We find: To our knowledge, this result is the first fully controlled extrapolation of λ c /m s to the continuum.

CONCLUSIONS
Localization of the low-lying Dirac eigenmodes in the high-temperature phase of QCD and other gauge theories  is closely related to the change in the confining properties of these theories taking place at the transition [12, 18-22, 25, 26, 28, 29, 32-34], and its study can lead to a better understanding of the mechanism behind deconfinement, and of its relation with chiral symmetry restoration.
The strong connection between deconfinement and localization is exemplified by the fact that the mobility edge, i.e., the point in the spectrum separating localized and delocalized modes in the spectrum, decreases as one approaches the pseudocritical temperature from above, and vanishes in the crossover range.In theories with a genuine deconfinement transition, this takes place exactly at the critical point [19-22, 25, 26, 28, 29].A more accurate quantitative determination of the "geometric" transition temperature where the mobility edge vanishes and localization disappears, and so of the tightness of the connection between localization and deconfinement, requires the full control of systematic effects, including finite volume and, especially, finite spacing effects.
Most of the numerical studies of localization combine the relation between the localization properties of eigenmodes and the statistical properties of the corresponding eigenvalues [36] with the use of the staggered discretization of the Dirac operator [3-5, 7, 17-20, 22, 25, 26, 28-30].However, this approach faces serious technical problems due to the restoration of taste symmetry in the continuum limit: as the lattice becomes finer the spectrum tends to organize in nearly-degenerate multiplets of eigenmodes, which in turn distort the spectral statistics away from their expected universal behavior, leading to hard-to-control systematic effects in the determination of the mobility edge.
In this paper we have carried out a dedicated study of the effects of taste degeneracy on the statistical properties of the staggered spectrum in high-temperature lattice QCD.We focused in particular on how these affect the numerical determination of the mobility edge, and how these effects change as the lattice spacing is reduced, or the aspect ratio is increased, with the main goal of providing a controlled continuum limit of the mobility edge.To this end, we studied in detail the possibility of exchanging the order of the continuum and thermodynamic limits, arguing that it applies to the study of localization properties and the determination of the mobility edge.
Our findings are in line with theoretical expectations, with a systematic overestimation of the mobility edge using spectral statistics on finer lattices, where the effects of taste degeneracy become sizeable also in the bulk of the spectrum.For larger aspect ratios at fixed lattice spacing these effects are reduced and the overestimation of the mobility edge is mitigated.In the thermodynamic limit, this estimate tends to the correct value of the mobility edge, obtained independently by studying the fractal dimension of the eigenvectors.Moreover, it agrees with finite-volume estimates obtained from spectral statistics using a suitably taste-symmetrized spectrum, for which the effects of taste degeneracy should be reduced.This supports our expectation that continuum and thermody-namic limits can be exchanged when studying the localization properties of Dirac eigenmodes.
From a practical perspective, the most important result of this analysis is that the infinite-volume extrapolation of the mobility edge on a finer lattice is in good agreement with the mobility edge found on coarser lattices, where taste-degeneracy effects do not reach into the bulk of the spectrum and a reliable estimate can be obtained from spectral statistics already at lower aspect ratios.This shows that accurate values for the mobility edge can be obtained on relatively coarse lattices for reasonable aspect ratios using spectral statistics, which is the ideal combination from the numerical point of view.
Moreover, our findings allowed us to perform the first, fully controlled extrapolation of the mobility edge to the continuum limit.This confirms the theoretical expectation [4,35] that the mobility edge in units of the quark mass is a renormalized quantity with physical meaning, and not only a lattice artifact.It also lends support to the numerical evidence for this fact presented in Refs.[4,29].There it was shown that the mobility edge depends only mildly on the lattice spacing.In the light of our results, this is because the calculations of Refs.[4,29] use relatively coarse lattices, where taste-degeneracy effects on estimates of the mobility edge based on spectral statistics are negligible.One might have wondered if the situation could have changed on finer lattices; our results show that this mild dependence would still show on finer lattices, provided one extrapolated first to the thermodynamic limit.(Incidentally, the extrapolation in T shown in Ref. [4] is in good qualitative agreement with our result.) In conclusion, we have shown that one can use staggered fermions efficiently and reliably to study the mobility edge in the Dirac spectrum of high-temperature lattice gauge theories, provided that the aspect ratio is sufficiently large to avoid sizeable taste-degeneracy effects in the relevant region of the spectrum.The possibility to bypass these effects through this procedure is due to the possibility of exchanging the continuum and thermodynamic limits when studying localization properties.Moreover, we have numerically demonstrated that the mobility edge in units of the quark mass is a renormalized quantity, in agreement with theoretical expec-tations; and that it depends only mildly on the lattice spacing.
It would be interesting to understand if and how one could avoid the well known problems of the staggered discretization (i.e., lack of good chiral properties, difficult relation with topology) by relating the staggered and the overlap spectrum on the same gauge configurations.In this context, the mobility edge could be used in two different ways.On the one hand, renormalizing the overlap spectrum (after matching a suitable observable to find the renormalization constant) would allow one to obtain another estimate for λ c /m s from spectral statistics unaffected by taste degeneracy, to be compared with the one obtained here after extrapolating to infinite volume.On the other hand, thanks to its mild dependence on the lattice spacing the mobility edge itself could be efficiently used to match the staggered and overlap spectra.
c (T ; N t )/m s , depends only mildly on the lattice spacing, λ true c (T ; N t )/m s ≃ λ cont c (T )/m cont s = lim Nt→∞ λ true c ) Notice that the contribution of delocalized unperturbed modes plays no role here.The first equation tells us that to leading order, E = E (0) ñ for some ñ, with nonzero c (0) n only in the corresponding degenerate unperturbed subspace.We then choose one of the four unperturbed basis vectors and set c n ̸ = ñ; the other three choices can be treated analogously.The second equation then tells us that E 1 = g∆ ññ , and for n such that E (0)

= 1 .
This surely does not lead to delocalization, and changes only marginally the size of the modes.While Re c (1) ñ = 0 due to the normalization condition, and Im c (1) ñ = 0 with a suitable choice of phases, c

8 FIG. 1 :
FIG. 1: Eigenvalues λ k /ms of the staggered operator in the background of a typical gauge configuration on a 32 3 ×8 lattice at T = 230 MeV.

3 14 FIG. 2 :
FIG.2: Top panel: comparison of ⟨s⟩ as a function of λ/ms for all the explored lattice spacings at fixed value of the aspect ratio LT = 4. Bottom panel: comparison of ⟨s⟩ as a function of λ/ms for all the explored aspect ratios, corresponding to different values of Ns, at fixed value of the lattice spacing, corresponding to a temporal extent Nt = 14.Our continuum estimate for the mobility edge [see Eq. (44) in Sec. 3 F] is marked in both panels by a red circle on the horizontal axis.

75 T 2 FIG. 3 :
FIG. 3: Mode counting in the localized regime of the spectrum.Red circles show the probability P (ν(M/ms) = n) of finding n modes below the cutoff M = 0.75ms, deep in the localized regime.Blue squares show the probability P (ν(M/ms)/2 = n) of finding n mode doublets below the cutoff M .Yellow stars show the corresponding Poisson distribution with parameter µ = ⟨ν(M/ms)/2⟩.

3 3 FIG. 4 :
FIG. 4: Comparison of Is 0 (top panel) and Var(s) (bottom panel), computed locally in the spectrum on lattices at fixed temperature T = 230 MeV, for different values of the lattice spacing at fixed spatial lattice size L = 3.4 fm in physical units (i.e., fixed LT = 4).Our continuum estimate for the mobility edge [see Eq. (44) in Sec. 3 F] is marked by a red circle on the horizontal axis.

14 FIG. 5 :
FIG. 5: Comparison of Is 0 (top panel) and Var(s) (bottom panel), computed locally in the spectrum on lattices at fixed temperature T = 230 MeV, for different spatial volumes at fixed lattice spacing, a = 0.0615 fm (i.e., Nt = 14).Our continuum estimate for the mobility edge [see Eq. (44) in Sec. 3 F] is marked by a red circle on the horizontal axis.

FIG. 6 :
FIG.6: Determination of the mobility edge for our 32 3 × 8 lattice from the dependence of Is 0 (top panel) and Var(s) (bottom panel) on the spectral region.

T 64 FIG. 8 :
FIG.8: Behavior of the fractal dimension as a function of the spectral bin computed for Nt = 14 for several choices of the lattice sizes Ns,1 and Ns,2.

FromFIG. 10 :
FIG.10: Dependence of our estimates for the renormalized mobility edge as a function of the lattice spacing at fixed temperature T = 230 MeV.The dashed lines represent continuum extrapolations of λc/ms assuming leading linear corrections in (1/Nt) 2 = (aT )2 , see text for details.Points have been slightly shifted horizontally to improve readability.

FIG. 11 :
FIG.11: Determination of the mobility edge from the spectral statistics of the spectrum after taste averaging (see text for details), estimated using results from 48 3 × 14 lattices.
FIG. 12:Estimates for the mobility edge from the spectral statistics of the spectrum after taste averaging (see text for details), obtained using results from our Nt = 14, 16 lattices.For comparison, we show also our results from the fractal dimension and from standard spectral statistics extrapolated to the thermodynamic limit on Nt = 14 lattices, and our continuum extrapolation based on standard spectral statistics (see Sec. 3 F).

TABLE I :
[81][82][83]arameters used in Ref.[75]to generate the physical-point gauge ensembles considered in this work.The bare parameters β, ams, am l = ams/28.15andthelattice spacing a have been fixed according to the results of Refs.[81][82][83], and correspond to physical pion mass and physical strange-to-light quark mass ratio.Simulations marked with * have been performed without the multicanonical algorithm.
3s ), up to lattice sizes where the effects of taste degeneracy dis-Dependence of our estimates for the mobility edge on the spatial size of the lattice, at fixed temporal extension Nt = 14.The dashed lines represent thermodynamic extrapolations of λ stat c