*Review* **Localization of Dirac Fermions in Finite-Temperature Gauge Theory**

**Matteo Giordano 1,\* and Tamás G. Kovács 1,2**


**\*** Correspondence: giordano@bodri.elte.hu

**Abstract:** It is by now well established that Dirac fermions coupled to non-Abelian gauge theories can undergo an Anderson-type localization transition. This transition affects eigenmodes in the lowest part of the Dirac spectrum, the ones most relevant to the low-energy physics of these models. Here we review several aspects of this phenomenon, mostly using the tools of lattice gauge theory. In particular, we discuss how the transition is related to the finite-temperature transitions leading to the deconfinement of fermions, as well as to the restoration of chiral symmetry that is spontaneously broken at low temperature. Other topics we touch upon are the universality of the transition, and its connection to topological excitations (instantons) of the gauge field and the associated fermionic zero modes. While the main focus is on Quantum Chromodynamics, we also discuss how the localization transition appears in other related models with different fermionic contents (including the quenched approximation), gauge groups, and in different space-time dimensions. Finally, we offer some speculations about the physical relevance of the localization transition in these models.

**Keywords:** localization; QCD; lattice gauge theory; finite temperature

#### **1. Introduction**

Quantum Chromodynamics (QCD) is currently our best microscopic description of strong interactions. As is well known, QCD is a gauge theory with gauge group SU(3), coupling six "flavors" of *quarks*, which are spin- <sup>1</sup> <sup>2</sup> Dirac fermions transforming in the fundamental representation of the group, to the eight spin-1 gauge bosons (known as *gluons*) associated with the local SU(3) symmetry. Despite their apparently simple form, the interactions of quarks and gluons, as dictated by the gauge principle and encoded in the Dirac operator, give rise to a wide variety of phenomena. Most notably, the low-energy properties of strongly interacting matter are largely determined by the phenomena of confinement and chiral symmetry breaking (see, e.g., Refs. [1–4]). At zero temperature, quarks and gluons are in fact confined within hadrons by a linearly rising potential, up to distances where a quark-antiquark pair can be created out of the vacuum. Furthermore, there is an approximate chiral symmetry associated with the lightest quarks, which in the limit of exactly massless quarks is broken spontaneously. This determines most of the properties of light hadrons, once the effects of the explicit breaking by the light quark masses is taken into account.

Confinement and chiral symmetry breaking persist also at nonzero but low temperature and densities. It is well established that, at vanishing chemical potential, QCD undergoes a finite-temperature transition to a deconfined, chirally restored phase (*quark-gluon plasma*), around *Tc* ≈ 155 MeV [5,6]. This transition is a rapid but analytic crossover [7], with both confining and chiral properties of the theory changing dramatically in a relatively narrow interval of temperatures.

In particular, the confining properties are determined by the fate of an approximate Z<sup>3</sup> center symmetry, i.e., a symmetry under gauge transformations which are periodic in

**Citation:** Giordano, M.; Kovács, T.G. Localization of Dirac Fermions in Finite-Temperature Gauge Theory. *Universe* **2021**, *7*, 194. https:// doi.org/10.3390/universe7060194

Academic Editor: Dmitri Antonov

Received: 28 April 2021 Accepted: 4 June 2021 Published: 8 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

time up to an element of the group center. At low temperature, center symmetry is only explicitly and mildly broken by the presence of quarks; at higher temperatures, instead, center symmetry is strongly broken spontaneously by the ordering of the Polyakov loop, i.e., the holonomy of the gauge field along a straight path winding around the temporal direction. In the "quenched" limit of infinitely heavy quarks, QCD reduces to SU(3) pure gauge theory, where center symmetry is an exact symmetry of the action, realized at low temperature *à la* Wigner-Weyl in the Hilbert space of the system, while spontaneously broken at high temperatures. Around the same temperature at which the Polyakov loop starts becoming ordered in QCD, corresponding roughly to the spontaneous breaking of center symmetry, also the chiral properties of the system change radically. The approximate order parameter of chiral symmetry, i.e., the chiral condensate *ψψ*¯ , decreases rapidly around *Tc*, corresponding to the disappearance of the effects of spontaneous breaking in the chiral limit. The net effect is an effective restoration of chiral symmetry, up to the explicit breaking due to the quark masses.

While the existence and the nature of the finite-temperature transition are by now well established, the mechanisms of confinement and chiral symmetry breaking, and similarly of deconfinement and chiral symmetry restoration, are not fully understood yet; nor is the apparently close relation between these two phenomena that are in principle completely unrelated. An important role in the breaking and restoration of chiral symmetry is played by the topological properties of the gauge field configurations.

It might be possible to understand the formation of a chiral condensate in the low temperature phase of QCD in terms of instantons and the associated zero modes [8–15]. It is well known that there is an exact zero mode of the Dirac operator associated with an isolated instanton or anti-instanton, and with their finite-temperature versions known as calorons [16–23]. Typical gauge field configurations can be interpreted as a more or less dense medium of instantons and anti-instantons. For a sufficiently high density of topological objects, the associated zero modes will strongly mix and form a finite band of near-zero Dirac modes, which in turn gives rise to a finite condensate via the Banks-Casher relation [24]. At higher temperatures the density of topological objects decreases, and so do the density of near-zero Dirac modes and the chiral condensate, until the symmetry is effectively restored. This is the "disordered medium scenario" for chiral symmetry breaking. The mechanism of confinement is understood less clearly, and various proposals have been put forth: we invite the interested reader to consult Refs. [1,2].

Relatively recently, a third phenomenon has been found to take place in correspondence with deconfinement and chiral symmetry restoration in QCD, namely the localization of the low-lying eigenmodes of the Dirac operator [25–40]. Localization is a widely studied subject in condensed matter physics, since Anderson's work on the absence of diffusion in random lattice system [41]. In his seminal paper, Anderson showed how the presence of disorder causes the spatial localization of energy eigenmodes. For electrons in a disordered medium, such as a conductor with impurities, localized modes appear at the band edge, beyond a "mobility edge" separating extended and localized modes. As the amount of impurities/disorder increases, the mobility edge moves towards the band center, eventually leading to all modes becoming localized, and to the conducting sample turning into an insulator. It is outside the scope of this paper (and frankly quite a Herculean task) to provide an exhaustive account of the developments in the theory of Anderson localization, and we refer the interested reader to the reviews [42–46].

It has been shown that a similar localization phenomenon takes place for the low-lying modes of the Dirac operatorin QCD above the pseudocritical temperature [27,28,30,33,34,39,40]: up to a critical point in the spectrum, i.e., the analogue of the "mobility edge", low modes are spatially localized on the scale of the inverse temperature [33,39]. The mobility edge depends on the temperature *T*, and its extrapolation towards the confined phase vanishes at a temperature compatible with *Tc* [33,40]. At the mobility edge, a second-order phase transition takes place in the spectrum, which has been shown to be a genuine Anderson transition [34–37].

A disorder-driven transition, such as the Anderson transition in the Dirac spectrum, obviously needs a source of disorder. This was first identified in the local fluctuations of the topological charge density, treated as a dilute ensemble of pseudoparticles (calorons) [25–28]. As already mentioned above, these topological objects individually support localized zero modes of the Dirac operator; since they overlap, the corresponding modes mix and shift away from zero, but for a dilute ensemble this effect is small and modes remain localized and near zero. While evidence was produced supporting the connection between localization and topological objects, it turned out that not all localized modes could be explained this way [31], and at least another source of disorder was needed. This was identified in the fluctuations of the Polyakov loop [31], a hypothesis supported by numerical results [31,39,40] and by the critical properties found at the mobility edge [34–37]. This led to the so-called "sea/islands picture" of localization, proposed in Ref. [31] and further elaborated in Refs. [47–49]: Dirac eigenmodes tend to localize on "islands" of Polyakov loop fluctuations away from its ordered value, which form an extended "sea" in the deconfined phase. The sea/islands picture requires only the existence of a phase with ordered Polyakov loop in order for localization to appear, and so leads one to expect the localization of the low Dirac modes in a generic gauge theory with a deconfinement transition. This has been verified in a variety of models [32,49–56], including ones without topology, thus providing further support to the sea/islands picture. This also clearly suggests a strong connection between localization of the low Dirac modes and deconfinement.

The relation between localization and chiral symmetry restoration has received less attention, mostly because of the intrinsic difficulty of studying gauge theories with massless fermions, where chiral symmetry is exact. It is, however, clearly established that localization of the low modes is accompanied by evident changes in the spectral density at the low end of the spectrum. In Refs. [51,54,57,58] it was observed that in the quenched theory a peak of near-zero modes of topological origin forms, followed by a spectral range with low mode density. A similar peak of localized modes was observed in Ref. [38] in the presence of dynamical fermions. (The presence of this peak was discussed before in Ref. [59], and has been studied recently in Refs. [60,61], although the localization properties of the eigenmodes are not studied there.) It is shown in Ref. [58] that the near-zero peak can be explained in terms of a dilute instanton/anti-instanton gas and the associated zero modes. Recently, it has been proposed that the presence of a finite density of near-zero localized modes in the chiral limit can lead to the disappearance of the finite-temperature massless excitations predicted by the finite-temperature version of Goldstone's theorem [62].

While the presence of localization at high temperature, and its connection with the finite-temperature transition in QCD and in other gauge theories are by now fairly well established, the physical meaning of localization of Dirac modes, and a detailed understanding of the aforementioned connection with deconfinement and chiral symmetry restoration, have proved to be quite elusive. The hope is that a better understanding of localization can shed light on the mechanisms of confinement and chiral symmetry breaking, and on the finite-temperature deconfining and chirally restoring transition.

In this review we will discuss developments in the study of the localization of Dirac modes in finite-temperature gauge theories. While the main focus is on QCD, we will discuss also other models, showing in particular that the connection between localization and deconfinement is a general phenomenon. As this review is aimed both at the particle physics and condensed matter communities (and expected to disappoint them both), we provide brief introductions to the subjects of finite-temperature QCD and Anderson localization in an attempt to bridge the gaps. Older results were already reviewed in Refs. [36,63].

Localization in gauge theories also occurs in a few other contexts, albeit those correspond to physical situations very different from the one discussed in the present review. For this reason here we do not discuss them in detail, and only mention them for completeness. A non-exhaustive list of references includes the following papers. For results concerning the relation between the localization properties of the low Dirac modes and

the topological structure of the vacuum in gauge theories, we refer to Refs. [64–68] and references therein. The role played by localization in the Aoki phase of quenched QCD with Wilson fermions, especially concerning the fate of Goldstone modes, is studied in Refs. [69,70]. Localization properties of the eigenmodes of the covariant Laplacian in Yang-Mills theories are studied in Refs. [71,72].

The plan of this paper is as follows. In Section 2 we review finite-temperature QCD and related issues. In Section 3 we review the topic of localization and Anderson transitions in some generality. In Section 4 we discuss localization in QCD at finite temperature. The disordered medium scenario and the sea/islands picture, providing mechanisms for localization, are discussed in Section 5. Localization in gauge theories other than QCD is discussed in Section 6. Finally, in Section 7 we draw our conclusions and show some prospects for the future.

#### **2. QCD at Finite Temperature**

In strongly interacting systems such as QCD localization takes place as the systems cross from the hadronic to the high-temperature quark-gluon plasma state. To put localization in QCD in the proper context, in the present section we summarize some basic facts about this finite temperature transition. For an introduction to gauge theories at finite temperature we refer the reader to the literature (see, e.g., Refs. [73,74]).

In an extended sense, QCD is a gauge theory with *Nf* flavors of quarks transforming in the fundamental representation of the gauge group SU(*Nc*) and interacting via the corresponding gauge field, the excitations of which are the gluons. In a stricter sense, in QCD *Nf* is fixed to six, the number of known quark flavors, and *Nc* = 3. At finite temperature, the theory is formally defined by the Euclidean partition function

$$Z\_{\rm QCD} = \int [dA] \, e^{-S\_{\rm YM}[A]} \prod\_f \det \left( \mathcal{D}[A] + m\_f \right), \tag{1}$$

where the product runs over the quark flavors with masses *mf* , *A* is the gauge field, *S*YM is the Euclidean Yang-Mills action and

$$\mathcal{B}[A] = \sum\_{\mu=1}^{4} \gamma\_{\mu} (\partial\_{\mu} + i \lg A\_{\mu}) \tag{2}$$

is the Euclidean Dirac operator, with *γμ* the Euclidean, Hermitean gamma matrices and *g* the coupling constant. The temporal direction is compactified to a circle of size equal to the inverse temperature, and periodic boundary conditions in the temporal direction are imposed on the gauge fields. In this form of the partition function the quark fields, appearing in the action quadratically, have been explicitly integrated out, resulting in the quark determinant. The Dirac operator is an anti-Hermitean operator with purely imaginary spectrum, which is furthermore symmetric about zero thanks to the property {*γ*5, *D*/} = 0.

It is instructive to consider the theory as a function of its parameters, that in reality are fixed by the observed properties of hadrons. The only such parameters of QCD are the quark masses.1 In particular, the low-energy properties of light hadrons are completely determined by the masses of the lightest quarks, the *u* and *d* quark, and to some extent the heavier *s* quark. The other three known quark flavors are so heavy that they have little influence on the low energy physics.

Quark masses are also important parameters from a theoretical point of view, because they crucially influence some symmetries of the system. Even though in nature these are only approximate symmetries, considering them helps to better understand the finite temperature transition of QCD. In an imaginary world with two massless quark flavors, i.e., when *mu* = *md* = 0, QCD would have an exact SU(2)*<sup>V</sup>* × SU(2)*<sup>A</sup>* × U(1)*<sup>A</sup>* chiral symmetry. Here SU(2)*<sup>V</sup>* is a rotation in the two-dimensional (*u*, *d*) flavor-space that acts identically on all the Dirac components. In contrast, the flavor non-singlet axial symmetry SU(2)*<sup>A</sup>* not

only mixes the two flavors, but also transforms the left and right Dirac components with opposite phases. Finally, the U(1)*<sup>A</sup>* flavor-singlet axial symmetry acts trivially in flavor space and rotates the left and right Dirac components with opposite phases. Even though these are all symmetries of the classical Lagrangian, after quantization the U(1)*<sup>A</sup>* part of the symmetry is anomalously broken. Furthermore, at zero temperature the SU(2)*<sup>A</sup>* axial symmetry is spontaneously broken. The emerging three Goldstone bosons are the analogues of the pions, and the order parameter of the symmetry breaking is the light quark condensate *ψψ*¯ . Finally, the vector part of the symmetry SU(2)*<sup>V</sup>* remains intact even for finite, but equal quark masses.

In reality, the nonzero and non-equal masses of the *u* and *d* quarks explicitly break these symmetries; however, the spontaneous and anomalous breaking inherited from the massless theory both turn out to be much stronger than this explicit breaking. In fact, in an imaginary world with zero *u* and *d* quark masses, the low-energy properties of the light hadrons would be much the same as they are in the real world. The only important exceptions would be the pions, which in that case would be exact Goldstone bosons with zero mass.

If one imagines changing the quark masses, the other interesting limit is the one in which quarks are much heavier than in reality. In particular, in the limit of infinitely heavy quarks the quark determinant in the path integral completely decouples and the back-reaction of the quarks on the gauge field disappears. This is the so-called quenched theory. In this limit QCD has a different exact symmetry, the symmetry group being the center of the gauge group, in the case at hand Z3. The symmetry transformation in question is a gauge transformation that is singular along a spacelike hypersurface, and its singularity is characterized by an element of the center Z3. Recalling that the system is finite in the temporal direction with periodic boundary conditions for the gauge field, this symmetry transformation is a gauge transformation that is not periodic in the temporal direction (hence singular), and it multiplies by the same Z<sup>3</sup> center element all the holonomies (gauge parallel transporters) going around the system in the temporal direction. The holonomies wrapping around the system in the temporal direction along a straight path are also called Polyakov loops.

Gauge invariant local gluonic quantities are defined in terms of holonomies around small loops, and those never wrap around the system. These types of loops cross the hypersurface where the gauge transformation is singular the same number of times in both directions. As a result, the Z<sup>3</sup> factors along such a loop always cancel, and gauge invariant local quantities are invariant with respect to the Z<sup>3</sup> center transformations. In contrast, fermionic quantities are affected, since such a singular gauge transformation essentially introduces an extra Z<sup>3</sup> twist for the temporal boundary condition of the fermions through the covariant derivative in the Dirac operator. The boundary condition affects the spectrum of the Dirac operator and also its determinant that appears in the path integral. At low temperatures where the correlation length is much smaller than the temporal extent of the system, and Polyakov loops fluctuate locally with little correlation, the temporal boundary condition has only a small impact on the Dirac determinant. Consequently, the quarks only mildly break the Z<sup>3</sup> symmetry. However, at high temperature, where the correlation length becomes comparable to or larger than the temporal size of the system, and the Polyakov loops tend to align with each other, this picture changes drastically.

To understand exactly how that happens, let us first recall that in finite temperature quantum field theory the temporal boundary condition for fermions is antiperiodic. For free massless fermions this implies a gap in the spectrum of the Dirac operator equal to the first Matsubara frequency. If by a singular gauge transformation (as defined above) we introduce an additional Z<sup>3</sup> twist in the boundary condition then the gap will decrease, because the twist *π* corresponding to the antiperiodic boundary condition will decrease to *π* ± 2*π*/3 = ±*π*/3 (mod 2*π*). In the interacting theory there is no gap in the spectrum, but through this mechanism the low end of the spectrum is much denser when the spatially averaged Polyakov loop is in the complex center sectors (i.e., close to one of the complex

center elements) than when it is in the real one (i.e., close to the identity), where the effective twist comes only from the antiperiodic boundary condition. As a result, the fermion determinant, that disfavors larger low-mode density, strongly favors the real Polyakov loop sector, and for finite quark mass this is the only sector that contributes to the path integral. This is how fermions explicitly break the Z<sup>3</sup> center symmetry and select the real Polyakov loop sector out of the three sectors that would be equivalent in their absence. This mechanism is at work also at low temperature, but much less effective there since the average Polyakov loop fluctuates around zero.<sup>2</sup>

Now going back to the quenched theory, at zero and low temperature, the exact Z<sup>3</sup> symmetry of its Lagrangian remains intact, while above a critical temperature this symmetry is spontaneously broken and its order parameter, the trace of the Polyakov loop, develops a nonzero expectation value. In fact, in the quenched theory, the logarithm of the expectation value of the Polyakov loop is proportional to the gauge field energy it costs to insert an infinitely heavy static quark in the system. In this way, the vanishing of the Polyakov loop in the low-temperature phase shows that no free quarks can exist there, so quarks are confined into hadrons. In contrast, the nonzero expectation value of the Polyakov loop in the high temperature phase implies that quarks are not confined there.

The nature of the finite-temperature transition in extended QCD is governed by the chiral and the Z<sup>3</sup> symmetries, which—as we have already seen—depend on the quark masses. In the quenched limit (infinite quark mass) lattice simulations have shown the transition to be weakly first order [78,79], and this behavior persists for large enough, but finite quark masses. If the quarks become lighter, the transition weakens and for intermediate quark masses there is a wide region where it is only a crossover. In particular, the light-quark masses in nature fall in this range [7]. For even smaller quark masses, the transition is again expected to become a true phase transition, but its order depends on the number of light quark flavors. For two light flavors (and physical strange quark mass) it is expected to be second order, whereas for three light flavors a first order phase transition is anticipated. However, the presence of these phase transitions, previously predicted based on an epsilon expansion [80] (see also [81] for the role played by the U(1)*<sup>A</sup>* anomaly), have not yet been confirmed by lattice simulations, because simulations close to the chiral limit are technically challenging.

Most of the results discussed in this review are based on numerical calculations on the lattice. Lattice field theory is a nonperturbative approach to the quantization of quantum field theories, based on the discretization of the relevant path integrals that define the theory in the path-integral approach. We provide here only a very brief introduction to this subject, referring the interested reader to the extensive literature (e.g., the books [73,74,82–84]). In the lattice approach to gauge theories devised by Wilson [85], the SU(3) gauge fields of QCD are replaced by unitary SU(3) matrices (*link variables*) associated with the links of a finite hypercubic lattice. In continuum language, these correspond to the parallel transporters of the gauge fields along the paths connecting neighbouring lattice points. After a suitable discretization of the gauge action, the relevant path integrals are obtained by integrating over the gauge fields, which in practice means integrating the link variables over the group manifold with the invariant (Haar) group measure. The desired, continuum field theory is obtained (if this is possible) by properly tuning the parameters in the action, so that the correlation length of the system in lattice units diverges, and the system "forgets" about the underlying lattice. For pure gauge theories, the only available parameter is the lattice inverse gauge coupling (usually denoted by *β*), which ceases to be a freely adjustable parameter and turns instead into a measure of the lattice spacing.<sup>3</sup>

The approach outlined above is easily generalized to other gauge theories based on different gauge groups, by simply replacing the SU(3) link variables and the corresponding Haar measure with elements of the relevant gauge group and the corresponding Haar measure. The inclusion of fermions instead is not straightforward, especially for what concerns the implementation of chiral symmetry. Nonetheless, there are several viable discretization of the Dirac operator, which are expected to all lead to the same results in the continuum limit. Since they appear below in Section 4, we mention Wilson fermions, staggered fermions (possibly rooted), domain wall fermions, overlap fermions, and twisted mass fermions (see Ref. [83,84] and references therein for details). We finally mention that several improvement schemes exist that bring the system closer to the continuum limit, i.e., that reduce the effects due to the finiteness of the lattice spacing. Such schemes exist both for the gauge action and for the fermionic determinant (see Ref. [83,84] and references therein for details).

#### *2.1. Finite-Temperature Transition, Dirac Spectrum, and Localization—An Overview*

We have seen that in the two extreme cases, the quenched limit and the chiral limit, two different symmetries, the Z<sup>3</sup> center symmetry and the SU(2)*<sup>A</sup>* axial symmetry govern the transition. The respective order parameters, the Polyakov loop and the quark condensate, signal spontaneous breaking of the symmetry in the high temperature phase for the Z<sup>3</sup> symmetry and in the low temperature phase for the SU(2)*<sup>A</sup>* axial symmetry. In nature, both symmetries are only approximate, the transition is a crossover and the order parameters have only inflection points in the crossover region. It also follows that in real QCD there is no sharply defined transition temperature. In contrast, regardless of the quark mass, the localization transition, i.e., the appearance of the first localized modes at the low edge of the Dirac spectrum, occurs at a sharply defined critical temperature. Moreover—as anticipated in the Introduction, and as we will see below in Sections 4 and 6—the localization transition occurs in the temperature range of the deconfining and chiral crossover in the case of real QCD, and exactly at the deconfining temperature in the quenched limit. This suggests that there might be a connection between the thermodynamic (chiral and deconfining) transitions on the one hand, and the localization transition on the other hand.

Besides the coincidence of their respective critical or pseudocritical temperatures, these phenomena are also connected through the degrees of freedom playing the most important role in their respective dynamics. When the system crosses into the high temperature phase, the spectral density of the Dirac operator around zero drops considerably, exactly vanishing in the chiral limit. This is how the chiral symmetry, spontaneously broken at low temperature, is restored above the transition. Indeed, through the Banks-Casher relation [24], the order parameter of chiral symmetry breaking, the quark condensate, is proportional in the chiral limit to the spectral density at zero, and generally strongly sensitive to the low end of the spectrum. Lower spectral density also means that eigenmodes close to each other in the spectrum are less likely to be mixed by fluctuations of the gauge field, which might lead to localization at the low end of the Dirac spectrum.

The spectral density, however, is not the only important parameter that influences localization. In the quark mass regions numerically explored so far, where the transition is either a crossover (near and below the physical values of the light quark masses) or a true phase transition governed by (approximate) center symmetry (heavy quark limit), the spectral density does not immediately drop to zero at the (pseudo)critical temperature. In particular, in the quenched limit just above the transition a narrow but tall spike at zero appears in the spectral density (see Figure 1). This is due to near-zero modes associated with a dilute gas of calorons and anticalorons, local fluctuations of the topological charge [58]. Even though the spectrum is dense in the spike, eigenmodes there are localized [86]. A similar peak of near-zero modes is also found for physical, near-physical, and belowphysical light-quark masses [38,59–61]. For near-physical masses these modes are found to be localized [38], and most likely this persists as the mass is decreased.

**Figure 1.** The spectral density, Equation (9) (here normalized by the volume), of the overlap Dirac operator in quenched QCD (i.e., SU(3) pure gauge theory), just above the phase transition at *T* = 1.045*Tc*. The grey band indicates the point separating the lowest, topological modes from the rest (its width equals the corresponding uncertainty). From Ref. [58]. (Figure adapted from R.Á. Vig and T.G. Kovács, arXiv:2101.01498 (2021), and used under a CC-BY 4.0 license (https://creativecommons.org/licenses/by/4.0)).

Recently, the spike in the spectral density received another interpretation. It was argued that it signals the appearance of a new, previously undiscovered "phase" of QCD, intermediate between the low-temperature confined and the high-temperature deconfined phase [87]. In a more recent paper the same authors studied a newly defined infrared dimension *d*IR of the eigenmodes in the low end of the spectrum of the chirally symmetric overlap Dirac operator. They concluded that the exact zero modes have *d*IR = 3, and in the spectral peak *d*IR changes rapidly but smoothly from 2 to 1 as one moves up in the spectrum [88]. This behavior persists up to the bulk of the spectrum, where the spectral density, together with the infrared dimension *d*IR of the modes starts to increase again. This nontrivial change in the infrared dimension all happens in the region where based on the spectral statistics and the scaling of the participation ratio with the volume, the eigenmodes are thought to be localized. It would be interesting to further investigate how *d*IR relates to the usual fractal dimension *D*<sup>2</sup> (see Equation (6)), and what kind of spatial structure in the eigenmodes gives rise to this nontrivial behavior. This could also depend on the chiral and locality properties of the particular discretization of the Dirac operator.

Topological fluctuations and the localization of the eigenmodes are both intimately related to fluctuations of the Polyakov loop, the order parameter of the quenched transition. The spatial localization of low Dirac eigenmodes is found to strongly correlate with local fluctuations of the Polyakov loop away from its symmetry-breaking equilibrium value [31,39,40]. This gives rise to the sea/islands picture of localization that we will discuss in Section 5.2 of the present paper in more details. Localization on calorons and localization on Polyakov loop fluctuations are, however, not mutually exclusive, as calorons always contain large fluctuations of the Polyakov loop. In fact, within a caloron, the Polyakov loop wraps around the gauge group in a topologically nontrivial way. The connection between calorons and Polyakov loop fluctuations is also shown by the strong correlation between the Polyakov loop and the topological susceptibility that can be observed close to the transition in the high temperature phase (see Figure 2).

**Figure 2.** The dependence of the topological susceptibility on the value of the spatially averaged Polyakov loop in quenched lattice simulations. The susceptibility was computed by dividing the configurations in sets according to the spatially averaged Polyakov loop. The averages for the whole ensembles (in two volumes) are shown by the horizontal bands, the widths indicating the uncertainties. From Ref. [58]. (Figure adapted from R.Á. Vig and T.G. Kovács, arXiv:2101.01498 (2021), and used under a CC-BY 4.0 license (https://creativecommons.org/licenses/by/4.0)).

The finite temperature transition of QCD is a result of the interplay of all those mechanisms that we just discussed, involving the Polyakov loop, the topological fluctuations and the spectral density of the Dirac operator around zero. Since localization is intimately related to all these aspects, it might hold the key to a better intuitive understanding of the dynamics of the transition.

#### **3. Localization and Anderson Transitions**

In a classic paper [41], Anderson showed that a sufficiently large amount of disorder in a lattice system prevents quantum-mechanical diffusion. Working in the one-particle tight-binding approximation, and mimicking the effect of disorder by supplementing the tight-binding Hamiltonian with a random potential on the lattice sites, Anderson showed that all the eigenfunctions of the system are localized for sufficiently strong disorder (i.e., for a sufficiently broad distribution for the random potential).

A practical example of this situation is a "dirty" crystal where some of the lattice atoms are replaced by impurities. Anderson's results imply that all the electron eigenstates become localized for a sufficiently large concentration of impurities. This prevents electron diffusion and the associated transport phenomena; in particular, the d.c. conductivity at zero temperature vanishes [89,90]. Localization then provides a possible mechanism for a disorder-induced metal-insulator transition (MIT).

Anderson's original arguments were later scrutinized and clarified by several authors [91–97]. Since then, the topic of disorder-induced localization, or *Anderson localization*, has been extensively studied in the condensed matter community, and it is impossible for us to provide here a comprehensive survey, or even do justice to the related literature. In this section we limit ourselves to a short review of the main aspects of Anderson localization, especially those relevant to gauge theories, discussed in the next Section. We invite the interested reader to consult the reviews [42–46].

#### *3.1. The Anderson Model*

In its simplest form, the (orthogonal) Anderson model Hamiltonian reads

$$H\_{\vec{x},\vec{y}}^{\text{AM}} = \varepsilon\_{\vec{x}} \delta\_{\vec{x},\vec{y}} + \sum\_{\mu=1}^{3} \left( \delta\_{\vec{x} + \mu,\vec{y}} + \delta\_{\vec{x} - \mu,\vec{y}} \right), \tag{3}$$

where *x*,*y* label the sites of a simple cubic lattice with lattice vectors *μ*ˆ, *μ* = 1, 2, 3, and *ε<sup>x</sup>* is a random on-site potential, with uniform probability distribution in the interval [<sup>−</sup> *<sup>W</sup>* <sup>2</sup> , *<sup>W</sup>* <sup>2</sup> ]. The lattice spacing and the hopping energy are set to 1 for simplicity. The width *W* of the distribution is a measure of the amount of disorder in the system, with *W* = 0 corresponding to a perfectly pure crystal. In this case, for a lattice of side *L* with periodic boundary conditions the eigenstates of *H*AM are plane waves with wave vectors *<sup>p</sup>* = <sup>2</sup>*πk L* , with *k<sup>μ</sup>* = 0, 1, ... , *L* − 1. However, as soon as even a small amount of disorder is put into the system, i.e., *W* = 0, the eigenmodes *ψ*(*x*) at the band edge become exponentially localized, i.e., |*ψ*(*x*)| <sup>2</sup> <sup>∼</sup> *<sup>e</sup>*−|*x*−*<sup>x</sup>*0|/*<sup>ξ</sup>* for *<sup>E</sup>* beyond critical energies <sup>±</sup>*Ec*(*W*) called "mobility edges" [90] (see Figure 3). As the amount of disorder *W* in the system increases, the mobility edge moves towards the band center. Eventually, for *W* larger than a critical disorder, *Wc*, all the modes become localized. If Equation (3) describes the conduction band of an electron in some "dirty" crystalline system, for large enough *W* the Fermi energy will lie in the localized part of the band; d.c. transport then takes place through hopping of electrons from one localized state to another, which has an exponentially small probability of happening, and in the limit of infinite size leads to the absence of charge transport. As the amount of impurities increases past the critical value, the system then undergoes a metal-to-insulator transition.

**Figure 3.** Sketch of the density of states *ρ* (see Equation (9)) as a function of energy *E* in the Anderson model, Equation (3). Localized modes are present in the shaded region beyond the mobility edges ±*Ec*.

#### *3.2. Anderson Transitions*

When the energy of the modes crosses the mobility edge *Ec*(*W*) at fixed *W*, or equivalently when the disorder in the system crosses the energy-dependent critical disorder *Wc*(*E*) at fixed mode energy *E*, the nature of the eigenmodes of *H*AM changes from delocalized to localized. As argued in Ref. [98], in three dimensions the associated transition is a second-order phase transition (*Anderson transition*), with divergent correlation length *ξ*(*E*) ∼ |*E* − *Ec*| <sup>−</sup>*<sup>ν</sup>* or *<sup>ξ</sup>*(*W*) ∼ |*<sup>W</sup>* <sup>−</sup> *Wc*<sup>|</sup> <sup>−</sup>*ν*, where the same exponent *ν* is expected.

This prediction is based on the so-called *scaling theory of localization* (see Ref. [43] for an introduction): first proposed in Ref. [98] based on previous ideas exposed in Refs. [42,99–101], it was later put on a firmer basis through a field-theoretical description of disordered systems and Anderson transitions [102–104] (see Ref. [45] for a full list of references). The basic idea is that the change in the conductance *G*(*L*) of the system4 as its size *L* is increased is controlled only by the localized or delocalized nature of the energy eigenmodes, which in turn is measured by the conductance itself, as a proxy for the disorder in the system. This implies a scaling behavior of the conductance, *<sup>d</sup>* ln *<sup>G</sup>*(*L*) *<sup>d</sup>* ln *<sup>L</sup>* = *β*(*G*(*L*)). Using the asymptotics of the *β* function obtained from localized or delocalized modes is then enough to show that in three dimensions there is an unstable fixed point (in the renormalization-group sense), and so a mobility edge in the energy spectrum and a phase

transition at some critical amount of disorder (see Figure 4). In one dimension no Anderson transition is expected as all modes are localized in the presence of disorder [105,106], while the situation in two dimensions is more complicated (see below).

**Figure 4.** The scaling function *β*(*g*) against the dimensionless conductance *g* = 2¯*hG <sup>e</sup>*<sup>2</sup> in various dimensions. From Ref. [98]. (Reprinted figure with permission from E. Abrahams, P.W. Anderson, D.C. Licciardello, and T.V. Ramakrishnan, Phys. Rev. Lett. 42, 673 (1979). Copyright (1979) by the American Physical Society).

The Anderson model that we just described, Equation (3), is but the simplest disordered Hamiltonian in three dimensions, and can be generalized in various ways. One can consider different probability distributions for the on-site disorder (e.g., the Lloyd model [107]), add off-diagonal disorder by making also the hopping terms random [108–110], increase the range of interaction by adding more hopping terms (e.g., Ref. [111]), and so on. However, according to the general theory of the renormalization group (see, e.g., Ref. [112]), critical properties at a second-order phase transition are shared by systems in the same *universality class*, determined only by general properties such as the dimensionality and the symmetries of the system.<sup>5</sup>

From a technical point of view, the Anderson model is a model of (sparse) random matrices. The properties of these models are the subject of Random Matrix Theory (RMT) [113–115]. For random systems, the relevant symmetry classification has been provided by Dyson [116], later extended by Verbaarschot [117], and completed by Altland and Zirnbauer [118–121] (see also Refs. [45,122]). The main *symmetry classes*, relevant to the models discussed in this review, are determined by the existence (or not) of an antiunitary symmetry operator *T* ("time reversal") commuting with the Hamiltonian, [*T*, *H*] = 0, and further specified by whether *<sup>T</sup>*<sup>2</sup> <sup>=</sup> 1 or *<sup>T</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>1. If *<sup>T</sup>* exists and *<sup>T</sup>*<sup>2</sup> <sup>=</sup> 1, the system is in the *orthogonal class* (*O*): this is the case for the model in Equation (3). If *T* does not exist, the system is in the *unitary class* (U). Perhaps the simplest example of a system in this class is the so-called unitary Anderson model (UAM),

$$H\_{\vec{x},\vec{y}}^{\text{UAM}} = \varepsilon\_{\vec{x}}\delta\_{\vec{x},\vec{y}} + \sum\_{\mu=1}^{3} (\delta\_{\vec{x}+\hat{\mu},\vec{y}} + \delta\_{\vec{x}-\hat{\mu},\vec{y}})e^{i\Phi\_{\vec{x},\vec{y}}}, \qquad \Phi\_{\vec{y},\vec{x}} = -\Phi\_{\vec{x},\vec{y}}.\tag{4}$$

which includes also off-diagonal disorder in the form of random phases *φx*,*<sup>y</sup>* in the hopping terms, mimicking the presence of a random magnetic field. Finally, if *<sup>T</sup>* exists and *<sup>T</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>1, the system is in the *symplectic class* (S). This classification is complete as far as the statistical properties in the bulk of the spectrum are concerned.

A refined classification is needed if one wants to discuss statistical spectral properties near the origin. In this case, one has to consider whether also a "particle-hole" symmetry exists, realized in terms of an antiunitary operator *C* obeying {*C*, *H*} = 0, and if so whether

*<sup>C</sup>*<sup>2</sup> <sup>=</sup> <sup>±</sup>1. This gives rise to nine different combinations. The eight combinations obtained when at least *T* or *C* exists correspond to eight different symmetry classes. If both *T* and *C* exist, it automatically follows that a unitary operator Γ = *TC* exists, anticommuting with the Hamiltonian, {Γ, *<sup>H</sup>*} <sup>=</sup> 0, and satisfying <sup>Γ</sup><sup>2</sup> <sup>=</sup> 1. However, a <sup>Γ</sup> satisfying this property can exist also if *T* and *C* are both absent. In this case there are two further symmetry classes, corresponding to whether such a Γ exists or not, for a total of ten. The classification is summarized in Table 1. In particular, if Γ exists and commutes with *T* (if this also exists), the system belongs to one of the chiral classes (chO, chU, and chS). Examples of systems of this type are provided by certain lattice models with random hopping terms, and no on-site potential, on bipartite lattices.

**Table 1.** Symmetry classes of random matrix ensembles. Entries corresponding to time-reversal (*T*) and particle-hole (*C*) symmetry indicate whether the symmetry is absent (0) or, if present, what is its square (±); entries corresponding to chiral symmetry (Γ) indicate whether it is absent or present (0 or 1).


#### *3.3. Detecting Localization: Eigenmode Observables*

A convenient way to study the localization properties of the eigenmodes of a random lattice Hamiltonian and how they change along the spectrum is by means of the *inverse participation ratios* (IPRs),6

$$\text{IPR}\_{\emptyset} \equiv \sum\_{\mathbf{x}} |\psi(\mathbf{x})|^{2q} \,\, \, \, \, \tag{5}$$

where it is assumed that eigenmodes obey the usual normalization condition, i.e., ∑*<sup>x</sup>* |*ψ*(*x*)| <sup>2</sup> = 1, and *x* now labels the sites of the relevant lattice, assumed finite and of volume *V* = *Ld*, with *L* the linear size and *d* the dimensionality. Unless specified otherwise, in the following both the term and the notation IPR, without subscript, will be used to refer specifically to the case *q* = 2. For modes extended throughout the whole system, one has qualitatively |*ψ*ext(*x*)| <sup>2</sup> <sup>∼</sup> 1/*V*, and so IPR*<sup>q</sup>* <sup>∼</sup> *<sup>V</sup>*1−*q*: after averaging over the possible realizations of disorder, which will be denoted with ... , and taking the large-volume limit, one has then IPR*<sup>q</sup>* → 0 as *V* → ∞ (for *q* > 1). For modes localized in a region of volume *V*<sup>0</sup> one has |*ψ*loc(*x*)| <sup>2</sup> <sup>∼</sup> 1/*V*<sup>0</sup> inside the localization region and negligible outside, and so IPR*<sup>q</sup>* <sup>∼</sup> *<sup>V</sup>*1−*<sup>q</sup>* <sup>0</sup> : one has then IPR*<sup>q</sup>* → const. as *V* → ∞. At the mobility edge, instead, the scaling of IPR*<sup>q</sup>* with the volume depends on *q* in a highly nontrivial way. One has in general

$$\text{IPR}\_q \sim L^{-(q-1)D\_q} \text{ },\tag{6}$$

where *Dq* = *d* for extended modes and *Dq* = 0 for localized modes, while at criticality *Dq* is not a constant. This reflects the multifractal nature of eigenmodes at *Ec* [124,125], and leads to define a set of multifractal exponents characterizing the critical behavior at the Anderson transition (see Ref. [45]).

Closely related to the IPR is the *participation ratio*,

$$\text{PR} \equiv \frac{1}{V} \text{IPR}^{-1} \qquad \left(= \frac{1}{V} \text{IPR}\_2^{-1}\right), \tag{7}$$

which measures the fraction of the system effectively occupied by the mode. For localized modes one has in the infinite volume limit PR → 0, while for delocalized modes extended throughout the system one finds PR → a nonzero constant. Another equivalent way to measure the localization properties is to use the mode "size", i.e., *<sup>V</sup>* · PR <sup>=</sup> IPR<sup>−</sup>1, which as *V* → ∞ (after averaging over the disorder) remains constant for localized modes and diverges for delocalized modes. For systems with nontrivial spin and/or internal degrees of freedom, the eigenvectors *ψα*,*c*(*x*) possess extra spin (*α*) and/or internal indices (*c*). In these cases it is convenient to employ a definition of the IPR which is invariant under spacetime and internal (unitary) rotations, i.e.,

$$\text{IPR} = \sum\_{\mathbf{x}} \left( \sum\_{\mathbf{a}, \mathbf{c}} \left| \psi\_{\mathbf{a}, \mathbf{c}}(\mathbf{x}) \right|^2 \right)^2 = \sum\_{\mathbf{x}} \left( \psi(\mathbf{x})^\dagger \psi(\mathbf{x}) \right)^2,\tag{8}$$

where the normalization condition ∑*<sup>x</sup> ψ*(*x*)†*ψ*(*x*) = 1 is understood. For example, for eigenmodes of the continuum, Wilson, or overlap Dirac operators, *α* = 1, ... , 4 is the Dirac index, and *c* = 1, ... , *Nc* is the gauge group ("color") index; for eigenmodes of the staggered operator *α* is absent but *c* is present.

#### *3.4. Detecting Localization: Eigenvalue Observables*

Another useful tool to detect localization are the statistical properties of the eigenvalues *λ<sup>i</sup>* of a random lattice Hamiltonian, which are closely related to the localization properties of its eigenvectors [126]. Localized modes are in fact expected to be sensitive only to local fluctuations in the disorder, and so the corresponding eigenvalues are expected to fluctuate independently. More precisely, after removal of non-universal, model-dependent features by means of the so-called *unfolding* procedure [113] (see below), the unfolded eigenvalues corresponding to localized modes should obey Poisson statistics. Delocalized modes, on the other hand, are expected to be mixed easily by fluctuations in the disorder, and so the corresponding unfolded spectrum should behave like that of a dense random matrix, and display the statistics of the Gaussian ensemble of Random Matrix Theory [113–115] in the appropriate symmetry class.

Unfolding is a monotonic mapping of the eigenvalues *λ<sup>i</sup>* that makes the spectral density equal to 1 throughout the spectrum. The spectral density is defined as

$$\rho(\lambda) \equiv \langle \sum\_{i} \delta(\lambda - \lambda\_{i}) \rangle \,. \tag{9}$$

The unfolded eigenvalues *xi* are given by

*λ<sup>i</sup>* → *xi* = *λ<sup>i</sup> dλ ρ*(*λ*), (10)

and it is easy to see that they have unit density, *ρ*¯(*x*) = *<sup>d</sup><sup>λ</sup> dx ρ*(*λ*) = 1. For random matrix models with dense matrices, the bulk statistical properties of the unfolded spectrum are expected to be universal (i.e., to not depend on the details of the model) and uniform throughout the spectrum. This has been proved rigorously for a large class of matrix ensembles (see Refs. [127,128] and references therein). One can then determine these properties in the exactly solvable Gaussian ensembles in the various symmetry classes (orthogonal, unitary, symplectic) [113]. In particular, the probability distribution of the unfolded spacings *si* = *xi*+<sup>1</sup> − *xi*, or *unfolded level spacing distribution* (ULSD), *p*ULSD(*s*), can be obtained exactly, although not in closed form. A good approximation for the ULSD is provided by the so-called *Wigner surmise* in the appropriate symmetry class [114] (see Figure 5),

$$p\_{\rm WS}^{(\beta)}(s) = a\_{\beta}s^{\beta}e^{-b\_{\beta}s^{2}},\tag{11}$$

where *β* is the *Dyson index* of the Gaussian ensemble, and one has for the various symmetry classes7

$$\begin{array}{llll}\text{orthogonal}: & \beta = 1, & a\_1 = \frac{\pi}{2}, & b\_1 = \frac{\pi}{4}, \\ \text{unitary}: & \beta = 2, & a\_2 = \frac{32}{\pi^2}, & b\_2 = \frac{4}{\pi}, \\ \text{ symplectic}: & \beta = 4, & a\_4 = \frac{262144}{729\pi^3}, & b\_4 = \frac{64}{9\pi}. \end{array} \tag{12}$$

For chiral classes, the same ULSD is found as for the corresponding non-chiral classes.

**Figure 5.** Wigner surmise for the various symmetry classes. The exponential distribution for Poisson statistics is also shown.

For independent eigenvalues, obeying Poisson statistics, the ULSD is the exponential distribution,

$$
\mathcal{p}\_{\text{Poisson}}(s) = e^{-s}.\tag{13}
$$

Both for localized and delocalized modes, exact analytical results are then available for the statistical properties of the unfolded spectrum, and the transition from one type of modes to the other can be easily monitored across the spectrum. This allows in particular to identify the mobility edge, where a different, critical statistics is expected instead of Poisson or RMT statistics [129,130]. Various families of random matrix models have been developed to describe the critical statistics [131–138].

We mention in passing an alternative approach to the study of universal statistical properties of the spectrum, based on the use of the ratio of consecutive level spacings [139]. Since this ratio is independent of the local spectral density, this approach has the advantage of not requiring any unfolding, and has been shown to provide more precise results than those obtained from the unfolded spectrum in a variety of many-body systems (see Ref. [140] and references therein).

#### *3.5. Finite-Size Scaling at the Anderson Transition*

The mobility edge and the correlation-length critical exponent *ν* can be obtained by means of a finite-size scaling study [130] of the average O*L*(*λ*) of suitable observables built out of the unfolded spectrum, and computed locally in the spectrum, for lattices of linear size *L*. Local statistics are defined formally as

$$\mathcal{O}\_L(\lambda) \equiv \frac{1}{\rho(\lambda)} \langle \sum\_i \delta(\lambda - \lambda\_i) \mathcal{O}(\lambda\_i) \rangle \equiv \langle \mathcal{O} \rangle\_{\lambda} \,\, \, \, \tag{14}$$

where O(*λi*) is some function of the eigenvalues (e.g., the level spacing Δ*λ<sup>i</sup>* = *λi*+<sup>1</sup> − *λi*, or the unfolded level spacing *si* = *xi*+1(*λi*+1) − *xi*(*λi*) and its powers), and the average is over the ensemble. The last equality sets an alternative notation for local averages.

In practice, for a finite sample of disorder realizations obtained numerically, local statistics are computed by dividing the spectrum in small bins and averaging over modes within those, as well as over the sample. Unfolding is done by first fitting some smooth *ρ*average(*λ*) to the numerical data and then applying Equation (10). Alternatively, one can sort the eigenvalues of the sample by magnitude and replace them by their rank divided by the number of disorder realizations. If one is interested only in unfolded level spacings, one can divide Δ*λ<sup>i</sup>* by the average level spacing Δ*λ <sup>λ</sup>* in the relevant spectral region, *si* = <sup>Δ</sup>*λ<sup>i</sup>* Δ*λ λ* (notice that in the infinite-volume, infinite-statistics limit one has Δ*λ <sup>λ</sup>* = 1/*ρ*(*λ*)).

For a finite-size scaling study, convenient observables are obtained from the unfolded level spacing distribution, *p*ULSD(*s*), defined above. Commonly used are the second moment, *<sup>s</sup>*<sup>2</sup> = <sup>∞</sup> <sup>0</sup> *ds p*ULSD(*s*)*s*2, and the integrated probability distribution *Is*<sup>0</sup> <sup>=</sup> *<sup>s</sup>*<sup>0</sup> <sup>0</sup> *ds p*ULSD(*s*). As the system size grows, O*L*(*λ*) tends to its value for Poisson statistics in spectral regions where modes are localized, and to its value for (the appropriate) RMT statistics in spectral regions where modes are delocalized. Near the mobility edge *λc*, renormalization-group arguments and the one-parameter scaling hypothesis [98] imply that O*L*(*λ*) depends on *λ* and *L* only through the combination *ξ*(*λ*)/*L*, where *ξ* is the correlation length,8 that diverges at *<sup>λ</sup><sup>c</sup>* like *<sup>ξ</sup>*(*λ*) ∼ |*<sup>λ</sup>* <sup>−</sup> *<sup>λ</sup>c*<sup>|</sup> <sup>−</sup>*ν*. Since <sup>O</sup>*L*(*λ*) is analytic in *<sup>λ</sup>* for finite *<sup>L</sup>*, it must then take the form <sup>O</sup>*L*(*λ*) = *<sup>f</sup>*((*<sup>λ</sup>* <sup>−</sup> *<sup>λ</sup>c*)*L*1/*ν*). Corrections to oneparameter scaling due to irrelevant operators can also be included, and the corresponding critical exponents be measured [141] (see Ref. [142] for an introduction). The goodness of one-parameter scaling can be visualized by means of the so-called "shape analysis" [143], obtained by plotting one spectral observable against another. If the scaling hypothesis is correct, only *ξ*/*L* should determine the statistical properties of the spectrum, and so points corresponding to different *λ* and system sizes should all lie on a single curve, corresponding to a path in the space of probability distributions connecting RMT and Poisson going through the critical statistics.<sup>9</sup> Thanks to the persistence of a remnant of multifractality near the mobility edge [144], one can apply similar finite-size scaling techniques also to the study of eigenmodes near criticality, in order to obtain the multifractal exponents [145], as well as the correlation-length exponent *ν* [146].

#### *3.6. Anderson Transitions in Specific Models: Analytic Predictions and Numerical Results*

Critical properties at the Anderson transition have been extensively studied by means of numerical simulations in the case of the conventional symmetry classes (O, U, and S), see Refs. [45,142,147,148] and references therein. According to the scaling theory of localization [98], a second-order Anderson transition is expected in all the conventional classes in three dimensions. The existence of these Anderson transitions has been confirmed numerically, and measurements of the correlation length critical exponent *ν* have shown that the three classes belong to different universality classes [141,148–152] (see Table 2). The expected nontrivial multifractal structure has also been found [145,146,148,151,152]. Universality has been explicitly demonstrated for the orthogonal class using different disorder distributions [141], and for the unitary class using different Hamiltonians [148,152].



In two dimensions, the predictions of the scaling theory of localization depend strongly on the details of the model. Absence of an Anderson transition is predicted in the orthogonal Anderson model, where all modes are expected to be localized for nonzero disorder, while an Anderson transition is predicted in the symplectic case (Ando model) [153,154]. The inclusion of topological effects in the field-theoretical description of disordered systems led one to expect an Anderson transition also in the theory of the integer Quantum Hall Effect [155], which belongs to the unitary class (see Ref. [156] for a review). While numerical evidence qualitatively supported this idea [157], significant quantitative discrepancies between different microscopic models were observed (see Ref. [158] for a summary), in contrast with the expected universality of the transition. A better understanding of the field theory describing the critical point was obtained only recently, in terms of a conformal field theory deformed only by marginal perturbations, that emerge from the spontaneous breaking of the replica (super)symmetry of the relevant nonlinear sigma model [159]. This proposal is quantitatively supported by numerical results, and can explain the apparent numerical discrepancies [158]. A transition between localized and delocalized modes was observed in the two-dimensional unitary Anderson model [160]. This transition is a disorder-induced transition of topological (Berezinski˘ı-Kosterlitz-Thouless [161–163]) type, with exponentially divergent correlation length, log *ξ* ∼ |*λ* − *λc*| <sup>−</sup>1/2, in contrast to the usual second-order transition. For the unitary Anderson model there are conflicting theoretical predictions: while perturbative contributions lead to all states being localized (see, e.g., Refs. [45,154]), the inclusion of nonperturbative terms can possibly lead to the presence of an Anderson transition (see references cited in Ref. [160]). A similar transition of topological type was also observed in a model for disordered graphene with strong long-range impurities [164], belonging to the orthogonal class in two dimensions.

For our purposes, it is important to discuss the effect of off-diagonal disorder on localization. In the orthogonal class, theoretical arguments [109,110] suggest that offdiagonal disorder alone cannot localize modes at the band center; only increasing the on-site disorder leads eventually to localization of all the modes. This is confirmed by numerical results in three dimensions for the orthogonal Anderson model with random hopping [165–168]. The mobility edges ±*Ec* separating extended and localized modes move towards the band center as the on-site disorder *W* is increased, and all modes are localized for *W* > *Wc*. In the absence of diagonal disorder (*W* = 0) for bipartite lattices, this model belongs to the chiral orthogonal class. The critical exponent *ν* characterizing the Anderson transition at *Ec* = 0 when *W* = 0 is found to be in agreement with that of the (non-chiral) orthogonal class, as well as with the one characterizing the transition at *E* = 0 as the critical on-site disorder *Wc* is reached [166,167].

The origin *E* = 0 is singled out when the system has chiral symmetry (which is always the case when the lattice is bipartite and only off-diagonal disorder is present). In two

dimensions, theoretical arguments predict critical behavior of modes at *E* = 0 (i.e., modes are extended but not fully delocalized) [111,169–175], while all other modes are localized. Numerical results indicate that modes are indeed critical at *E* = 0 [174,176–181], but an Anderson transition to localized modes can also appear [174,182]. It has been argued that such an Anderson transition can be present due to non-perturbative, topological effects [183]. In three-dimensional models with chiral symmetry, Anderson transitions at the origin (*Ec* = 0) are expected to show critical properties differing from those of the correponding conventional class (and from those at *Ec* = 0). Ref. [184] studied the Anderson transition at *E* ∼ 0 in a chiral unitary model with purely off-diagonal disorder, finding multifractal exponents differing from those of the corresponding non-chiral class. In Refs. [185,186] the Anderson transition at the origin was studied in two-band models with on-site disorder in the chiral orthogonal and chiral unitary classes (as well as in other non-conventional symmetry classes), finding correlation length critical exponents differing from those of the corresponding non-chiral classes (and not entirely universal).

The critical properties of Anderson transitions at *Ec* = 0 in systems with chiral symmetry are instead not expected to differ from those in the corresponding conventional classes. Ref. [187] provides evidence of localization near the band center in a chiral unitary model mimicking fermions in a background of correlated spins with antiferromagnetic coupling in three dimensions, with the same critical properties as the non-chiral class; in two dimensions states near the band center seem instead to remain extended. Anticipating the results discussed in the following Sections, Refs. [27,34,35,37,51,53] provide examples of models in the chiral unitary class, both in three and two dimensions, displaying Anderson transitions at finite energy, and showing the same critical properties at the mobility edge as the corresponding non-chiral classes.

#### **4. Localization and Deconfinement in QCD at Finite Temperature**

The study of localization in QCD was initially motivated by the idea that the spontaneous breaking of chiral symmetry could have a similar origin as conductivity in a disordered medium [8–15]. The basic idea of the disordered medium scenario is that the near-zero modes responsible for the breaking of chiral symmetry originate from the mixing of the zero modes associated with overlapping instantons (or, more precisely, calorons at finite temperature). At finite temperature these zero modes are exponentially localized on the scale of the inverse temperature. If instantons/calorons overlap sufficiently, mixing of the corresponding zero modes will transform the zero eigenvalues into a near-zero band of levels, and lead to delocalized eigenmodes [10].<sup>10</sup> This is analogous, for example, to the Anderson-Mott insulator-metal transition<sup>11</sup> driven by the impurity concentration in doped semiconductors (see, e.g., Ref. [189]). Starting from the chirally symmetric phase and decreasing *T*, the overlap of calorons increases and eventually leads to a finite density of near-zero, delocalized modes. This leads to expect a localization-delocalization transition in the near-zero region.

As we will see below, this scenario is most likely only a part of the story, and overlooks the important role played by deconfinement in localizing the low Dirac modes. Instead of sticking to the disordered medium scenario, we prefer to adopt a more general point of view, looking at the Dirac operator as a random matrix, ignoring initially any relation with topological objects and deconfinement. After a few introductory remarks on this approach, we review the available results regarding localization and Anderson transitions in various lattice approximations for QCD, following a chronological order. Some of the results discussed here deal with the pure gauge theory, sometimes for *Nc* = 2, and are included in this section mostly for historical reasons. A summary of the results for QCD proper and organized by topic is provided at the end of the section.

#### *4.1. The Dirac Operator as a Random Matrix*

The Dirac operator in the background of fluctuating gauge fields can be intepreted as a sparse random matrix, and the properties of its eigenvalues and eigenvectors can be studied with the machinery discussed in the previous Section. For the continuum anti-Hermitian Dirac operator, −*iD*/ can be formally treated as the Hamiltonian of a disordered system, with disorder provided by the gauge fields. If an Anderson transition is present in its spectrum, its critical properties are expected to be determined by the symmetry class of the Dirac operator and by the dimensionality of the space-time over which it is defined.

Concerning the symmetry class, the four-dimensional Dirac operator for fundamental fermions in SU(*Nc*) theories belongs to the chiral unitary class for *Nc* > 2, and to the chiral orthogonal class for *Nc* = 2.12 The spectral correlations of Dirac eigenvalues in QCD (*Nc* = 3) are then expected to display GUE-type bulk statistics,<sup>13</sup> as long as the corresponding eigenvectors are delocalized.<sup>14</sup> If localized modes are present, they are expected to obey Poisson statistics, regardless of the symmetry class.

The discretization of the Dirac operator on a lattice is known to be tricky due to the doubling problem (see Refs. [73,74,82–84]), and in some cases its chiral properties are changed (see Ref. [188] and Ref. [115], Section 5.2.1). For staggered fermions [196–198] a remnant of the continuum chiral symmetry preserves the chiral nature of the symmetry class, and the staggered Dirac operator belongs to the chiral unitary class, as the continuum operator, for *Nc* ≥ 3. For *Nc* = 2 the symmetry class is instead changed to the chiral symplectic one.<sup>15</sup> Overlap fermions [199–202] possess an exact lattice chiral symmetry, and belong to the same symmetry class as their continuum counterpart. More precisely, since the overlap operator is not anti-Hermitean, this is true for its anti-Hermitean part. It is then understood, unless specified otherwise, that the imaginary part of the overlap eigenvalues is considered in the following. For the low modes this is an adequate approximation, that becomes exact in the continuum limit. Moreover, since the unfolded spectrum is unaffected by any monotonic mapping, the statistical properties of the low modes are unchanged if one uses other types of projection on the imaginary axis (i.e., eigenvalue magnitude, stereographic projection).16

Concerning the dimensionality of the problem, in finite-temperature field theory the temporal size of the system is fixed (in physical units) in the thermodynamic limit, and only the size of the *d* spatial directions is sent to infinity in the thermodynamic limit. For *d* + 1 spacetime dimensions, the dimensionality of the disordered system described by the Dirac operator in a gauge-field background is then equal to *d*, while the temporal direction can be technically seen as an internal degree of freedom.

#### *4.2. Numerical Results on the Lattice*

The disordered medium scenario was investigated in Ref. [25] by means of numerical simulations of quenched QCD on the lattice on both sides of the finite-temperature transition. They used a single spatial volume, employed the staggered discretization of the Dirac operator, and studied the rotation- and gauge-invariant version of the IPR of an eigenmode *<sup>ψ</sup>*, Equation (8). They observed that in the physical Z<sup>3</sup> sector (real Polyakov loop sector) the IPR of the lowest modes was considerably larger above the transition than below the transition (see Figure 6). Moreover, above *Tc* it was larger in the real sector than in the complex Polyakov loop sectors, where it does not change much across the transition. Sensitivity to the Polyakov loop sector is equivalent to sensitivity to the temporal boundary conditions, and shows that the low modes cannot be localized in the temporal direction on a scale much shorter than the temporal size.17 This suggests the presence of a localization transition in the physical sector, with spatially localized low modes at high temperature. Evidence for some of the localized modes being related to calorons was also provided. More evidence for localization of the low modes in the real sector appeared in Ref. [26], where more volumes and a chirally improved discretization of the Dirac operator were used. The volume scaling of the IPR of the low modes in the real sector was found to be in qualitative agreement with that expected for localized modes.

**Figure 6.** IPR of the low staggered modes on quenched configurations below (**left**), slightly above (**center**), and well above (**right**) the deconfinement transition for real (*θ<sup>P</sup>* <sup>=</sup> 0) and complex (*θ<sup>P</sup>* <sup>=</sup> <sup>±</sup><sup>2</sup> <sup>3</sup> ) Polyakov loop sectors. From Ref. [25]. (Reprinted figure with permission from M. Göckeler, P.E.L. Rakow, A. Schäfer, W. Söldner, and T. Wettig, Phys. Rev. Lett. 87, 042001 (2001). Copyright (2001) by the American Physical Society).

The disordered medium scenario was investigated further by García-García and Osborn in Refs. [27,28]. In Ref. [27] they considered an Instanton Liquid Model (ILM) for the QCD vacuum (see Ref. [204]), and studied the behavior of the instantonic zero modes. Changing the temperature, and so the spatial extension of the zero modes, they observed the appearance of a mobility edge near the origin, both in the quenched approximation and in the presence of fermions. In the quenched case, the multifractal properties of the near-zero modes at the transition were found to be consistent with those of the 3d unitary Anderson transition (see Refs. [45,148,152]). With two massless flavors, the mobility edge appears at the same temperature where the chiral condensate shows a drop (see Figure 7, left). Although the thermodynamic limit was not studied, this was taken as an indication that localization of the low modes coincides with the chiral transition.

**Figure 7. Left**: IPR of the lowest Dirac mode and chiral condensate in an ILM model for QCD with two massless fermions for various system sizes (in fm3). From Ref. [27]. (Reprinted from Nucl. Phys. A, 770, A.M. García-García and J.C. Osborn, "Chiral phase transition and Anderson localization in the Instanton Liquid Model for QCD", Pages 141–161, Copyright (2006), with permission from Elsevier). **Right**: IPR (here times the volume *V* = *L*3) of the low Dirac modes and chiral condensate in 2+1 flavor QCD with staggered fermions. From Ref. [28]. (Reprinted figure with permission from A.M. García-García and J.C. Osborn, Phys. Rev. D 75, 034503 (2007). Copyright (2007) by the American Physical Society).

A test of the disordered medium scenario in a more realistic context was presented in Ref. [28]. There the localization properties of the near-zero modes were studied on the lattice in quenched QCD, i.e., pure gauge SU(3) theory, and "unquenched" QCD, i.e., with 2+1 flavors of dynamical quarks of relatively large masses, leading to heavier-than-physical pions [205,206]. The one-loop Symanzik improved gauge action was used, and the Asqtadimproved [207–210] staggered discretization was employed for the lattice Dirac operator. In both cases, they found indications of critical (i.e., volume-independent) spectral statistics

(from the second moment of the ULSD) at a temperature *T*loc *<sup>c</sup>* , where also IPR · *V* starts increasing with the volume (for the unquenched case see Figure 7, right). In the quenched case, indications of a vanishing spectral density and of an increase of the Polyakov loop are found at a similar temperature *<sup>T</sup>*dec/*<sup>χ</sup> <sup>c</sup>* , identified with the deconfinement temperature (in the physical Z<sup>3</sup> sector). In the unquenched case, *<sup>T</sup>*loc *<sup>c</sup>* is close to the crossover temperature *Tχ <sup>c</sup>* obtained from the chiral susceptibility. Although the use of small lattices does not allow a full quantitative assessment, these indications suggest that an Anderson transition takes place near the origin of the spectrum as the system crosses over from the low-temperature to the high-temperature phase, with the low-lying Dirac modes turning from delocalized to localized, and the formation of a mobility edge that separates them from delocalized modes in the bulk of the spectrum.

Studies of localization in QCD-like settings includes also the case of two flavors of dynamical staggered fermions [29], and that of quenched two-color QCD (i.e., pure gauge SU(2) theory) analyzing overlap [30,31] and staggered [32] spectra. In the two-flavor three-color case Ref. [29] found that IPR · *V* of the low modes was volume-independent below the transition temperature *T*<sup>2</sup> *<sup>f</sup> <sup>c</sup>* , but above that it scaled with the volume in a manner compatible with localization (see Figure 8, left). Ref. [30] shows evidence of absence of correlations in the low-lying overlap spectrum, which is typical of localized modes, at *<sup>T</sup>* <sup>=</sup> 2.6*T*SU(2) *<sup>c</sup>* , where *<sup>T</sup>*SU(2) *<sup>c</sup>* is the deconfinement temperature of the pure gauge SU(2) theory. Ref. [32] shows clear evidence of localization of the low-lying staggered modes, and of the presence of a mobility edge separating them from delocalized bulk modes, again at *<sup>T</sup>* <sup>=</sup> 2.6*T*SU(2) *<sup>c</sup>* . This is obtained by studying how (i) the scaling with the lattice spatial volume of the spatial "size" (IPR<sup>−</sup>1/*Nt*) 1 <sup>3</sup> of the eigenmodes, and (ii) the ULSD of the corresponding eigenvalues change along the spectrum. The spatial extension of the low modes is volume-independent, while higher up in the spectrum it is seen to increase with the lattice size (see Figure 8, right). Looking at the ULSD in different spectral regions, it is observed that it matches the exponential distribution of Poisson statistics for the lowest modes, changing towards the symplectic Wigner surmise<sup>18</sup> as one moves towards the bulk. Finally, in Ref. [31] the transition in the overlap spectrum from localized to delocalized modes is studied via the ULSD at *<sup>T</sup>* <sup>=</sup> 2.6*T*SU(2) *<sup>c</sup>* . A clear change from the exponential to the orthogonal Wigner surmise is observed.19 Moreover, assuming that there are no strong interactions among instantons and anti-instantons, it is argued that the instanton density is too low to match the density of localized modes at this temperature. Indications of correlations between localized modes and local fluctuations of the Polyakov loop away from its ordered value (i.e., 1, in the physical sector) are also reported (see Section 5.2 for a detailed discussion).

**Figure 8. Left**: IPR (here times the volume *V* = *L*3) of low staggered modes for various volumes at *T* = 2*T*<sup>2</sup> *<sup>f</sup> <sup>c</sup>* in QCD with two flavors of dynamical staggered fermions. From Ref. [29]. (Reprinted figure with permission from R.V. Gavai, S. Gupta, and R. Lacaze, Phys. Rev. D 77, 114506 (2008). Copyright (2008) by the American Physical Society). **Right**: Mode size for low staggered modes in the background of quenched SU(2) configurations at *<sup>T</sup>* <sup>=</sup> 2.6*T*SU(2) *<sup>c</sup>* . From Ref. [32]. (Reprinted figure with permission from T.G. Kovács and F. Pittler, Phys. Rev. Lett. 105, 192001 (2010). Copyright (2010) by the American Physical Society).

A comprehensive study of localization in the high-temperature phase of real-world QCD appeared in Ref. [33], using a tree-level Symanzik improved gauge action and a two-level stout smeared [211] staggered fermion action for 2+1 quark flavors with physical mass [212]. Several volumes, aspect ratios and lattice spacings were used, covering the temperature range 1.7*Tc* < *T* < 5*Tc* (here *Tc* = 155 MeV [5]) with lattices of linear size 2 fm ≤ *L* ≤ 6 fm. Localized modes were observed at the low end of the spectrum (see Figure 9, left). A temperature-dependent mobility edge *λc*(*T*) separating low-lying, localized modes from delocalized bulk modes was found in the whole temperature range, studying how the spectral statistics change along the spectrum from Poisson to unitary RMT type. More precisely, *λ<sup>c</sup>* was estimated as the inflection point of the variance of the ULSD, *<sup>s</sup>*<sup>2</sup> *<sup>λ</sup>* − *s* 2 *<sup>λ</sup>* <sup>=</sup> *<sup>s</sup>*<sup>2</sup> *<sup>λ</sup>* − 1, computed locally in the spectrum. As *T* increases, *λc*(*T*) increases as well. Extrapolation to the continuum is studied at *T* = 400 MeV. The mobility edge is expected to renormalize like a quark mass, and the ratio *λc*/*m*ud is indeed shown to be independent of the lattice spacing within numerical uncertainties. The localization length *<sup>l</sup>* <sup>≡</sup> *<sup>a</sup>* IPR<sup>−</sup> <sup>1</sup> 4 of the low modes is also shown to extrapolate to a finite continuum limit, and *lT* is found to be between 0.7 and 0.9 for all the lattice ensembles. A second-order polynomial fit to the RG-invariant quantity *λc*(*T*)/*m*ud shows that it extrapolates to zero at *T*loc *<sup>c</sup>* = 170 MeV (see Figure 9, right), which is within the temperature range where the system undergoes a crossover from the low-temperature phase to the high-temperature phase [5,6]. This is consistent with localization of the low modes appearing as the system changes from confined and chirally broken to deconfined and chirally restored. The density of localized modes (number of modes per unit spatial volume) is seen to increase with *T* (see Figure 13 below).

**Figure 9.** Average PR of the staggered Dirac modes as a function of the eigenvalue *λ* for various volumes at *T* = 394 MeV, *a* = 0.125 fm (**left**) and renormalized mobility edge *λc*/*m*ud as a function of *T* (**right**) in 2+1 flavor QCD with staggered fermions and physical quark masses. From Ref. [33]. (Reprinted figures with permission from T.G. Kovács and F. Pittler, Phys. Rev. D 86, 114515 (2012). Copyright (2012) by the American Physical Society).

The critical behavior of the eigenmodes at the mobility edge was studied in Refs. [34,35,37]. All these references use the same setup as Ref. [33] with *Nt* = 4 and *a* = 0.125 fm, corresponding to *T* = 2.6*Tc*. In Ref. [34] it was established, by means of a finite size scaling analysis of the integrated ULSD *Is*<sup>0</sup> , that the transition from localized to delocalized modes at the mobility edge is indeed an Anderson transition (see Figure 10, left). The critical exponent was found to be *ν* = 1.43(6), in agreement with the one obtained for the 3d unitary Anderson model [149] (see Table 2). In Ref. [35] the critical eigenvalue statistics at the mobility edge was studied in terms of the one-parameter family of deformed random matrix ensembles of Refs. [136,137]. The critical statistics was shown to be indeed volume-independent, and well described by a deformed random matrix ensemble, with deformation parameter consistent with the one found for the 3d unitary Anderson model [148,149,152]. Finally, in Ref. [37] the critical exponent *ν* and the multifractal exponents where studied using the finite size scaling techniques for the eigenmode density developed in Refs. [145,146]. All exponents were found to be in agreement with those of the 3d unitary Anderson model [148,152] (see Figure 10, right).

While mostly focussed on the properties of the spectrum, Ref. [38] briefly discussed the localization properties of the low Dirac modes in QCD near the crossover temperature. The spectrum of the overlap operator was studied in the background of gauge configurations generated with tree-level improved Symanzik gauge action and 2+1 flavors of highly improved staggered quarks (HISQ) [213] with near-physical quark masses (*ml*/*ms* = 1/20, *m<sup>π</sup>* = 160 MeV). Evidence was found of a small peak of localized near-zero modes at *T* = 1.2*Tc* and *T* = 1.5*Tc* (here *Tc* = 154 MeV). Localization was inferred from the smallness of the PR; the volume scaling was not discussed. It was suggested that near-zero modes in the peak correspond to an approximate superposition of the exact zero modes associated with instanton–anti-instanton pairs. Comparison of the PR of zero and near-zero modes shows however that only a fraction of near-zero modes is compatible with this interpretation. The large fluctuations of the PR of the near-zero modes suggests instead a large variability of the number of topological lumps participating in the superposition.

**Figure 10. Left**: scaling of *Is*<sup>0</sup> near the mobility edge in 2+1 QCD with staggered quarks, at *T* = 394 MeV and *a* = 0.125 fm. The linear size *L* of the lattice is expressed in lattice units. Here *λ<sup>c</sup>* 0.3364, *ν* 1.43, and the critical value is *Is*0|*λ<sup>c</sup>* 0.194. From Ref. [34]. (Reprinted figure with permission from M. Giordano, T.G. Kovács and F. Pittler, Phys. Rev. Lett. 112, 102002 (2014). Copyright (2014) by the American Physical Society). **Right**: multifractal exponents *Dq*, governing the scaling of IPR*<sup>q</sup>* <sup>∼</sup> *<sup>L</sup>*−*Dq* (*q*−1) at large linear size *<sup>L</sup>* at the mobility edge, in 2+1 QCD with staggered quarks (*<sup>T</sup>* <sup>=</sup> <sup>394</sup> MeV, *<sup>a</sup>* <sup>=</sup> 0.125 fm) and in the 3d unitary Anderson model. From Ref. [37]. (Reprinted figure with permission from L. Ujfalusi, M. Giordano, F. Pittler, T.G. Kovács and I. Varga, Phys. Rev. D 92, 094513 (2015). Copyright (2015) by the American Physical Society).

Localization in two-flavor QCD was studied in Ref. [39] using tree-level improved Symanzik gauge action and dynamical Möbius domain-wall fermions [214–216] with stout smearing, and looking at the spectrum of the Hermitian operator *γ*5*D*, with *D* the four-dimensional effective Dirac operator of the five-dimensional domain-wall fermion. This operator is in the chiral unitary class. The temperature range was 0.9*Tc* ≤ *T* ≤ 1.9*Tc* with *Tc* 175 MeV the deconfinement temperature estimated from the average Polyakov loop. A range of bare quark masses, two spatial volumes and two temporal extensions (in lattice units) were used. Above *Tc*, the scaling of the PR of the eigenmodes shows that the lowest modes are localized, while moving up in the spectrum modes become delocalized (see Figure 11, left). The size *v* ≡ *V* · PR of the low modes increases along the spectrum and decreases with *T* (see Figure 11, right). The localization length *l* ≡ *v* 1 <sup>4</sup> of the lowest (nonzero) mode was shown to be of the order of the inverse temperature, *lT* ∼ 1.3. The ULSD computed locally in the spectrum was seen to change from Poisson-type to RMT-type in the unitary class as one moves from the lowest modes towards the bulk. By contrast, below *Tc* RMT statistics was observed everywhere in the spectrum. Changing boundary conditions to periodic in time, low modes were found to be delocalized with

RMT statistics also above *Tc*. A clear correlation between the spatial density *ψ*†*ψ*(*x*) of the low modes and the local fluctuations of the Polyakov loop *P*(*x*) away from its ordered value was observed, favoring sites with Re tr*P*(*x*) close to −1, and becoming stronger as *T* increases (see Figure 12 below). Correlation with action (*s*(*x*)) and topological charge (*q*(*x*)) densities was also observed, with localized modes favoring sites with large *s* and *q*, in particular "(anti)self-dual" sites where |*q*|/*s* ∼ 1. The overlap of the left- and right-chirality components of the modes was seen to be the smallest for the lowest modes, and to increase as one moves towards the bulk; it was also seen to increase with temperature, and showed little dependence on the bare quark mass.

In Ref. [40] (see also Ref. [217]), localization was studied in 2+1+1 flavor QCD with physical strange and charm masses but heavy pions (*m<sup>π</sup>* 370 MeV), looking at the (stereographically projected) spectrum of the overlap operator in the background of configurations generated with Iwasaki gauge action and dynamical twisted-mass Wilson fermions [218–220]. Above *Tc* 188 MeV, a very small PR is found for the lowest modes, which increases to around 0.8 in the bulk. The position of the mobility edge was estimated as the inflection point of the PR in the spectrum. As a function of *T*, *λc*(*T*) appears to be linear, and its extrapolation vanishes at a temperature compatible with *Tc*. A strong anticorrelation of the localized modes with Re tr*P*(*x*) was observed.

**Figure 11.** Scaling of the PR (and of *V* · PR in the inset) for *β* = 4.18, *Nt* = 8 (*T* = 257 MeV 1.5*Tc*) and bare mass *m* = 0.01 (**left**), and dependence on *T* and *m* of the PR of the 10 lowest modes (**right**) in two-flavor QCD with Möbius domain-wall fermions. From Ref. [39]. (Figures adapted from G. Cossu and S. Hashimoto, J. High Energy Phys. 06, 56 (2016), and used under a CC-BY 4.0 license (https://creativecommons.org/licenses/by/4.0)).

#### *4.3. Summary*

Let us summarize the results discussed in this section, restricting to QCD "proper", i.e., gauge group SU(3) and dynamical quarks. QCD-like models and other gauge theories are discussed in detail below in Section 6.


localization appears somewhere in the range of temperatures where the crossover takes place.


We now list a few remarks.


by the correlation with Polyakov-loop fluctuations, action and topological density, and chirality. A direct identification of monopoles or a quantitative estimate of their density is, however, unavailable.

**Figure 12.** Density plot of the average local norm *ψ*†*ψ* of low Dirac modes in the Polyakov loop plane (here P = <sup>1</sup> <sup>3</sup> tr *P*) in the high-temperature phase. From Ref. [39]. (Figure adapted from G. Cossu and S. Hashimoto, J. High Energy Phys. 06, 56 (2016), and used under a CC-BY 4.0 license (https://creativecommons.org/licenses/by/4.0)).

#### **5. Mechanisms for Localization**

In this section, we discuss in some detail the two mechanisms, or more precisely the two sources of disorder, proposed so far to explain localization of the low Dirac modes in QCD: the disordered medium scenario, based on topology fluctuations; and the sea/islands picture, based on fluctuations of the Polyakov loop. We have already briefly discussed the disordered medium scenario in the previous section; here we discuss it again, both to keep this section self-contained, and to give more details.

#### *5.1. The Disordered Medium Scenario*

As is well known, the continuum Dirac operator in the background of a gauge configuration of topological charge *Q* has *n*<sup>±</sup> exact zero modes of definite chirality ±1, with *Q* = *n*<sup>+</sup> − *n*<sup>−</sup> (index theorem). In particular, for instantons (resp. anti-instantons) of topological charge 1 (resp. −1) one finds an exact zero mode of positive (resp. negative) chirality. The same holds for the finite-temperature generalization of instantons known as calorons [16–23].<sup>22</sup> The zero modes supported by instantons (i.e., at *T* = 0) are algebraically localized, decaying like 1/*R*<sup>3</sup> with the distance *R* from the instanton center. The zero modes supported by calorons are instead fully delocalized in the temporal direction, and exponentially localized in the spatial directions, decaying like *e*−*rT*, with *r* the spatial distance from the caloron center and *T* the temperature of the system. For a dilute ensemble of these objects, their associated zero modes are not exact Dirac eigenmodes any longer, due to the fact that instantons/calorons overlap. The low-lying Dirac eigenmodes are instead linear combinations of these "unperturbed" zero modes,23 obtained by diagonalizing the "perturbed" Dirac operator, which in the zero-mode basis reads (see, e.g., Ref. [204])

$$i\!i\!\!\!\!\!\/\!\/= \begin{pmatrix} \mathbf{0} & T\_{IA} \\ T\_{AI} & \mathbf{0} \end{pmatrix} \;\prime\,\tag{15}$$

with *TIA* and *TAI* = *T*† *IA* the matrices of the overlap integrals of *iD*/ between the unperturbed zero modes associated with an instanton–anti-instanton pair.24 For a dilute ensemble the total topological charge *Q* is expected to be simply equal to the sum of the individual charges. Out of all the unperturbed zero modes, *Q* are preserved by topology despite mixing,<sup>25</sup> while the remaining ones are not protected by topology and spread around *λ* = 0 forming a band. The extent of this spreading and the resulting density of near-zero Dirac modes for typical gauge configurations are dynamical issues, determined by the typical density and size of topological objects. In the quenched case, a finite spectral density of nearzero modes is expected to survive as long as a non-negligible density of topological objects supports them. In the presence of dynamical fermions, the fermionic determinant tends to suppress configurations with a higher density of near-zero modes, and so suppresses topological excitations, with respect to the quenched case, but a finite spectral density is still possible. In any case, the details of the dynamics, including especially the temperature and the fermion masses, determine whether a nonzero density of near-zero modes is formed, i.e., loosely speaking, whether chiral symmetry is spontaneously broken.

For sufficiently low temperature, and not too many quark flavors, the density of topological objects and the effect of mixing become strong enough to overcome the repulsive effect of the fermionic determinant, and chiral symmetry breaks spontaneously through the formation of a nonzero density of near-zero modes. This is the disordered medium scenario for chiral symmetry breaking [8–15]. The mixing of the unperturbed modes is also expected to spread them out in space, over topological objects that overlap non-negligibly with their original location. If the typical spatial distance *n*<sup>−</sup> <sup>1</sup> <sup>3</sup> between topological objects is large compared to the typical spatial range 1/*T* of the corresponding unperturbed zero modes, one expects the resulting perturbed near-zero modes to remain localized on a few objects only. Here *n* = <sup>N</sup>top *<sup>V</sup>* is the spatial density of Ntop calorons and anti-calorons in a finite spatial volume *V*. At high temperatures both density and range are small, and near-zero modes are expected to be localized. As the temperature decreases, both density and range increase, with more and more topological objects overlapping, and near-zero modes are expected to eventually delocalize over the whole system [10]. It is reasonable to expect that delocalization will take place around the same temperature as chiral symmetry breaking (in the loose sense explained above). It should be clear, however, that finite spectral density and delocalization of modes near the origin are not automatically linked.26

According to the scenario above, near-zero localized modes should be associated with local lumps of topological charge. It is worth noting that (anti)calorons in SU(*Nc*) gauge theory are made up of *Nc* (anti)monopole constituents [20], and that when these are well separated the associated zero mode is localized on a single constituent; which one depends on the holonomy (Polyakov loop) of the gauge field at asymptotic distance from the core [222,223]. For typical high-temperature ordered configurations with Polyakov loop in the trivial sector, the relevant constituent is the type-*L* monopole, which also has the largest action and topological charge densities, as well as the smallest size (see Ref. [221]). This further characterizes the favorable locations for modes according to the disordered medium scenario.

Refs. [25,26,38,39,86] provide evidence that *some* of the modes are indeed localized on topological objects. In particular, Ref. [39] shows that the locations favored by some of the localized low modes have all the features of *L*-type monopoles and antimonopoles: large action and topological charge densities, near (anti)self-duality, and near degeneracy of two eigenvalues of the untraced Polyakov loop (see Figure 12).<sup>27</sup> In Refs. [58,86] it is shown that for pure gauge SU(3) theory the distribution of the number of near-zero modes in the peak of the spectral density near zero (see Figure 1) is consistent with the distribution of a dilute gas of topological objects. This suggests that the peak of near-zero modes indeed originates from the zero modes associated with topological objects. This provides further evidence supporting the disordered medium scenario as a viable mechanism for localization, and its close relation with spontaneous chiral symmetry breaking.

If localization were entirely due to mixing topological would-be zero modes, then the density of localized modes would be equal to that of the topological objects (calorons).28 Since above the transition the density of calorons decreases sharply with increasing temperature, we would expect the same behavior of the density of localized modes. However, as shown in Ref. [33], the density of localized modes actually increases with temperature (see Figure 13). In the pure gauge case, no more than half of the localized modes seem to be of topological origin for temperatures as low as 1.03*Tc* [86] and as the temperature increases, this fraction rapidly decreases. In Figure 14 we show the temperature dependence of the fraction of localized modes that can be associated with near-zero modes of topological origin [51,86]. We show results obtained with the overlap Dirac operator, and with the staggered Dirac operator for three different values of the lattice spacing. All the results are consistent and show that with increasing temperature a rapidly decreasing fraction of the localized modes are of topological origin. We can conclude that while topology-related localized modes may suffice to explain the near-zero peak, they cannot explain all the remaining localized modes found in a typical high-temperature gauge configuration, which therefore require a different supplementary mechanism.

**Figure 13.** Density of localized modes (modes per cubic fermi) in high-temperature QCD. From Ref. [33]. (Reprinted figure with permission from T.G. Kovács and F. Pittler, Phys. Rev. D 86, 114515 (2012). Copyright (2012) by the American Physical Society).

**Figure 14.** The fraction of localized modes in the quenched theory that are topology-related near-zero modes. The figure shows data calculated with the overlap Dirac operator on *Nt* = 6 lattices, as well as with the staggered Dirac operator using three different lattice spacings corresponding to *Nt* = 6, 8 and 10. From Ref. [86]. (Figure adapted from T.G. Kovács and R.Á. Vig, PoS LATTICE2018, 258 (2019), and used under a CC-BY-NC-ND 4.0 license (https://creativecommons.org/licenses/by-nc-nd/4.0)).

#### *5.2. The Sea/Islands Picture*

An alternative mechanism has been proposed in Ref. [31], and further elaborated in Refs. [47–49]. The basic observation is that the eigenvalues of the untraced Polyakov loop at a spatial site *x* effectively change the temporal boundary condition felt locally by the quark eigenfunctions *ψ*(*t*,*x*). Indeed, working for simplicity in the continuum, if one fixes the gauge to the temporal gauge *A*4(*t*,*x*) = **0**, the eigenvalue problem reduces to

$$\mathcal{D}\psi = (\partial\_4 \gamma\_4 + \mathcal{D}\_{(3)})\psi = i\lambda\psi, \qquad \mathcal{D}\_{(3)} = \sum\_{j=1}^{3} \gamma\_j (\partial\_j + i\mathfrak{g}\mathcal{A}\_j) \, , \tag{16}$$

while the effect of a nontrivial (untraced) Polyakov loop *P*(*x*) is to change the temporal boundary condition from antiperiodic to

$$
\psi(1/T, \vec{x}) = -P(\vec{x})\psi(0, \vec{x})\,. \tag{17}
$$

Equations (16) and (17) define the eigenvalue problem in temporal gauge. Clearly, the effective local boundary condition Equation (17) affects the contribution of site *x* to the Dirac eigenvalue, *iλ* = (*ψ*, *D*/*ψ*).

To gain some insight on the effects of the Polyakov loop, it is useful to study a family of configurations for which the eigenvalue problem can be solved exactly, namely those with *Aμ*(*t*,*x*) = **0** everywhere, and with a constant but nontrivial Polyakov loop. This can always be diagonalized by means of a global gauge transformation, so without loss of generality we can take *P*(*x*) = diag(*eiφ*<sup>1</sup> ,*eiφ*<sup>2</sup> ,*eiφ*<sup>3</sup> ), with *ei*(*φ*1+*φ*2+*φ*3) = 1. On these configurations the eigenfunctions of <sup>−</sup>*D*/<sup>2</sup> are plane waves,

$$
\psi\_c^{(a,k,\vec{p})}(t,\vec{x}) = \delta\_{ca} e^{i\left(\omega\_{ak}t + \vec{p}\cdot\vec{x}\right)}\,,\tag{18}
$$

where *c* is the color index,<sup>29</sup> with temporal frequency (*effective Matsubara frequency*) given by

$$
\omega\_{ak} = T[(2k+1)\pi + \phi\_a], \quad k \in \mathbb{Z}, \tag{19}
$$

and corresponding eigenvalues

$$
\lambda\_{ak} \left( \vec{p} \right)^2 = \omega\_{ak}^2 + \vec{p}^2 \,. \tag{20}
$$

Restricting without loss of generality to *φ*1,2 ∈ (−*π*, *π*], *φ*<sup>1</sup> + *φ*<sup>2</sup> + *φ*<sup>3</sup> = 0, the lowest positive Dirac eigenvalue is seen to be

$$
\lambda\_{\min} = T(\pi - \max\_{a} |\phi\_{\mathbb{I}}|) \, \prime \tag{21}
$$

i.e., it decreseas monotonically and symmetrically as one moves away from *φ<sup>a</sup>* = 0 ∀*a*, and vanishes when at least one of the Polyakov loop eigenvalues equals −1.

While the configuration discussed above is obviously unrealistic, the result Equation (21) allows understanding qualitatively which sites will be favored by a low Dirac eigenmode when the Polyakov loop configuration is mostly ordered near *P*(*x*) ≈ **1**, with "islands" of fluctuations in the "sea" of ordered Polyakov loops, as it happens at high temperature. One can in fact interpret Equation (21), now with *x*-dependent phases *φ<sup>a</sup>* = *φa*(*x*), as a sort of three-dimensional local potential for the quarks, to which one should add the appropriate "hopping terms" originating from the spatial dependence of the Polyakov loops, as well as from the spatial components of the gauge potential. From this point of view, fluctuations of the Polyakov loop away from order provide regions of lower potential that can "trap" the quarks.

More precisely, at high temperatures *φa*(*x*) ≈ 0 in an extended region, and neglecting in a first approximation the effect of the islands and of hopping, one finds fully delocalized modes. In the same approximation one finds a spectral gap, with the corresponding (positive) eigenvalues starting at *Tπ*, i.e., the usual lowest (fermionic) Matsubara frequency. The presence of islands and the effect of the interactions are expected to reduce this gap,

but the lowest eigenvalue that can be reached by a delocalized mode is expected to remain separated from the origin. On the other hand, localizing on an island of fluctuations can be "energetically" more favorable, and bring the corresponding eigenvalues inside the gap, as long as the gain in potential energy achieved by avoiding the sea of ordered Polyakov loops is sufficiently larger than the price paid for localization in terms of spatial momenta. This leads to expect the following scenario, at least when the islands are sufficiently distant from each other: a region of low spectral density, or *pseudogap*, opens between the origin and some point *λ<sup>c</sup>* in the spectrum, or *mobility edge*, above which modes are extended throughout the whole space; modes in the pseudogap can exist only if they are localized on energetically convenient islands of fluctuations.

The scenario described above, which has been dubbed the sea/islands picture of localization, shows a clear similarity with the Anderson-type models of condensed matter physics, and the terminology has been chosen precisely to reflect this similarity. From the point of view of random Hamiltonians, the local fluctuations of the Polyakov loop provide a three-dimensional source of on-site disorder. This would naturally explain the fact that the critical behavior found at the mobility edge in QCD is the same as that of the threedimensional unitary Anderson model.<sup>30</sup> There are, however, important differences with the simple unitary Anderson model of Equation (4). In that case, localization starts from the band edges and moves towards the band center as the amount of disorder, as measured by the width of its probability distribution, is increased. In QCD, while localization may as well be present at the band edges (cf. the results of Ref. [27] in the ILM model and Ref. [56] in Z<sup>2</sup> gauge theory), it is its appearance directly at the band center that characterizes the deconfined phase. Moreover, the actual source of disorder are the eigenvalues of the Polyakov loops, which are complex numbers lying on the unit circle, and so the magnitude of the disorder is actually bounded.

Another important difference, which is relevant also to the problem of spontaneous chiral symmetry breaking, is the different structure of the "free" Hamiltonian associated with the two cases in the absence of fluctuations. For the Anderson model this is simply *H*(AM) free <sup>=</sup> <sup>∑</sup><sup>3</sup> *<sup>j</sup>*=<sup>1</sup> *Tj* + *<sup>T</sup>*† *<sup>j</sup>* with *Tj* the translation operator in direction *j*, while for the Dirac operator one has (for a naive lattice discretization) *D*/free = ∑<sup>4</sup> *<sup>μ</sup>*=<sup>1</sup> *γμ*(*T<sup>μ</sup>* <sup>−</sup> *<sup>T</sup>*† *<sup>μ</sup>*). Here periodic boundary conditions are understood in the spatial directions; the effective boundary conditions Equation (17) are assumed for the temporal direction. While the spectrum of *H*(AM) free , *E*(*<sup>p</sup>*) = <sup>∑</sup>*<sup>j</sup>* cos *pj*, with *pj* <sup>=</sup> <sup>2</sup>*πkj <sup>L</sup>* , is dense near the origin, the presence of the gamma matrices in *D*/free leads to *λak*(*<sup>p</sup>*) = *ω*2 *ak* +*p* 2. Even in the case *ωak* = 0, in which there is no sharp spectral gap, the spectral density near the origin is low and vanishes at *λ* = 0.

An important aspect of this scenario is that deconfinement is naturally associated with the two effects that lead to localization of the low modes in high-temperature QCD. The first such effect is of course the formation of a sea of ordered Polyakov loops close to the identity, which can cause the opening of a spectral pseudogap and so make modes that localize on the islands of fluctuations stable against delocalization, as explained above. However, the appearance of the pseudogap requires also a second effect due to the ordering of Polyakov loops at deconfinement, namely the increased correlation between gauge fields on different time slices. The discussion of this effect requires a more detailed description of the Dirac operator in the language of Anderson models, in what can be called the "Dirac-Anderson approach". Here we sketch the discussion in the continuum in the temporal gauge; a more detailed and mathematically more precise analysis is presented in Ref. [48] for staggered fermions on the lattice.

Due to its compactness, the temporal direction can be treated as an internal degree of freedom, in particular by expanding the quark eigenfunctions on a complete basis of plane waves *eiωak <sup>t</sup>* obeying the appropriate effective boundary conditions, Equation (17). Here *ωak* are the effective Matsubara frequencies of Equation (19), now *x*-dependent, which provide a random on-site potential of the form *ωak*(*x*)*γ*4. For every color *a* with associated

Polyakov-loop phase *φa*(*x*), in correspondence to each wave number *k* there is a different branch of the on-site potential, and so a different associated three-dimensional Andersontype model, built by adding the on-site disorder to the spatial part of the Dirac operator (projected on the *a*, *k* subspace). We will refer to each of these models as a Dirac-Anderson model. The full Dirac operator is obtained by putting the various Dirac-Anderson models together, and by including their coupling induced by the hopping terms (i.e., the spatial part of the operator). The strength of the coupling among the Dirac-Anderson models turns out to be inversely related to the correlation of the gauge fields on different time slices. At low temperatures this correlation is small, the Dirac-Anderson models are strongly coupled, and the internal degree of freedom is effectively one more direction in which the modes can extend, thus facilitating their delocalization. This is in agreement with the effectively four-dimensional nature of QCD in the low temperature phase. As a matter of fact, the pseudogap does not open at low temperatures, where the spectral density is finite near zero, and this can only happen if the various Dirac-Anderson models do mix with each other (see the discussion about the free Dirac operator). In the absence of a pseudogap, localization of a mode is generally unstable against mixing with modes of similar energies. At high temperature, instead, the Polyakov loops become ordered inducing stronger correlations among different time slices, and the Dirac-Anderson models decouple making the problem effectively three-dimensional. In particular, the pseudogap is now expected to appear: it would be present for exactly decoupled Dirac-Anderson models (see again the discussion about the free Dirac operator), and their limited mixing is not sufficient to close it. Localized modes near the band center can then be supported by Polyakov loop fluctuations, as discussed above. As shown in Ref. [48] in a toy model where ordering of the Polyakov loop and correlation of the time slices can be varied independently, both effects are required for localization to appear at the band center.

An important aspect of the sea/islands picture is that it is not incompatible with the growing density of localized modes observed in QCD. Differently from the case of topological charge, Polyakov-loop fluctuations are not quantized. As *T* grows in the deconfined phase, the volume *V*fluct occupied by Polyakov-loop fluctuations is expected to decrease, *<sup>V</sup>*fluct <sup>∼</sup> *<sup>T</sup>*−*c*<sup>1</sup> , as the Polyakov loop becomes more and more ordered. On the other hand, the typical size *V*<sup>0</sup> of the islands of fluctuations is also expected to decrease, *<sup>V</sup>*<sup>0</sup> <sup>∼</sup> *<sup>T</sup>*−*c*<sup>2</sup> . The number of localized modes is expected to be directly related to the number of islands, *<sup>V</sup>*fluct/*V*<sup>0</sup> <sup>∼</sup> *<sup>T</sup>c*2−*c*<sup>1</sup> , and whether this number increases or decreases with temperature depends on the details of the dynamics. For example, while increasing in QCD [33] up to *T* ∼ 5*Tc*, it is seen to decrease in 2+1-dimensional SU(3) gauge theory above *T* ∼ 1.1 ÷ 1.2*Tc* [53].

Perhaps the most appealing feature of the sea/islands picture is its simplicity: all that it needs to work is the ordering of the Polyakov loop. This leads immediately to expect that localized modes will appear at the low end of the spectrum whenever an ordering transition takes place, independently of details such as the gauge group and its representation, fermionic content, nontrivial topological features, dimensionality,<sup>31</sup> and so on. This is discussed in the next section.

#### **6. Localization in Other Gauge Theories**

In this section, we discuss localization of the low Dirac modes in gauge theories other than QCD. Some of the references have been already discussed in Section 4 in connection with QCD, where they where treated as approximations. Here they are briefly discussed again, focussing more on the differences than on the similarities with QCD.

The main motivation in studying more general gauge theories is to investigate further the extent of the connection between localization on one side, and deconfinement and chiral restoration on the other. In particular, studying localization in models with genuine deconfining and/or chirally restoring phase transitions allows one to investigate this connection in a more clear-cut setting than in QCD, where it is somewhat blurred by the crossover nature of the transition. Studying more general gauge theories also allows one to test the sea/islands picture discussed in the previous section, and its generic prediction of localization of low modes in the high-temperature, "ordered" phase.

Genuine deconfining phase transitions are found in pure gauge SU(*Nc*) theory. In 3+1 dimensions the transition is second order for *Nc* = 2 and first order for *Nc* ≥ 3, while in 2+1 dimensions it is second order for *Nc* = 2, 3 and first order for *Nc* ≥ 4 (see, e.g., Refs. [78,79,224,225]). As already mentioned in Section 4, localized low Dirac modes have been found in pure gauge SU(2) theory in 3+1 dimensions, both with the overlap [30,31] and with the staggered [31,32] Dirac operator (see Figure 8, right), above the deconfinement temperature *Tc*. (Further details can be found in Section 4 and will not be repeated here.) No sign of localized modes was found instead below *Tc*. From the random-matrix point of view, the SU(2) case differs from SU(*Nc* ≥ 3) as the symmetry class is the symplectic instead of the unitary one. This is reflected in the different behavior of the unfolded spectrum in the bulk, which agrees with the symplectic Wigner surmise, Equations (11) and (12). A detailed study of the Anderson transition was not pursued.

Results for pure gauge SU(3) have been presented in Refs. [51,54,86] for the 3+1-dimensional case, and in Ref. [53] for the 2+1-dimensional case. Localized low Dirac modes are found in both cases in the deconfined phase. In 3+1 dimensions the temperature dependence of the mobility edge *λ<sup>c</sup>* was studied using the Wilson gauge action both with staggered [51] and overlap [54] fermions (using in this case the magnitude of the eigenvalues), smearing the gauge fields with two steps of stout smearing [211] in the staggered case, and two steps of hex smearing [226] in the overlap case. The integrated ULSD computed locally in the spectrum, *Is*<sup>0</sup> (*λ*), was used to determine *λ<sup>c</sup>* as the point where *Is*<sup>0</sup> takes its critical value *I* (*c*) *<sup>s</sup>*<sup>0</sup> [34], i.e., *Is*<sup>0</sup> (*λc*) = *I* (*c*) *<sup>s</sup>*<sup>0</sup> . Here use was made of the universality of the critical properties of the Anderson transition, which should be shared by QCD and pure gauge SU(3) theory, as they are both in the 3d unitary class. This was confirmed by the volume-independence of the resulting *λc*. For both discretizations, *λ<sup>c</sup>* is seen to extrapolate to zero at a temperature which agrees with the deconfinement temperature (see Refs. [78,79] and references therein) within numerical errors (see Figure 15).

**Figure 15.** Mobility edge in pure gauge SU(3) theory in the staggered (**left**) and overlap (**right**) Dirac spectrum. From Refs. [51,54]. (Figures adapted from T.G. Kovács and R.Á. Vig, Phys. Rev. D 97, 014502 (2018) (left), and from R.Á. Vig and T.G. Kovács, Phys. Rev. D 101, 094511 (2020) (right), and used under a CC-BY 4.0 license (https://creativecommons.org/ licenses/by/4.0)).

The 2+1-dimensional case was studied in Ref. [53] using the Wilson gauge action and the staggered discretization (without smearing). Universality arguments lead to expect that the Anderson transition is of BKT type with exponentially divergent correlation length, as found in Ref. [160] for the 2d unitary Anderson model. The results of Ref. [53] support this scenario. In particular, spectral statistics are critical, i.e., volume independent for all *λ* above *λc*, as expected for a BKT-type Anderson transition [227], see Figure 16, left. The mobility edge was determined by means of a finite size scaling study, and found to

extrapolate to zero at a temperature compatible with the deconfinement temperature [225] (although with much larger numerical uncertainty), see Figure 16, right. In the confined phase no localization was found, but low modes were seen to display a nontrivial fractal dimension *D*<sup>2</sup> < 2 (see Equation (6)).

**Figure 16.** Spectral statistics *Is*<sup>0</sup> along the spectrum for various volumes at coupling *β*¯ = 6.25 (**left**), and mobility edge as a function of *β*¯ (**right**) for the staggered operator in 2+1-dimensional pure gauge SU(3) theory. (Here *β*¯ = *β*/3 with *β* the Wilson action coupling.) In the left panel, Poisson and RMT predictions and the position of the mobility edge are also shown. In the right panel, a power-law fit (black solid line) to *λc*, the position (blue solid line) and error band (blue dashed lines) of *β*¯ loc at which *λ<sup>c</sup>* extrapolates to zero, and the error band of the critical *β*¯ *<sup>c</sup>* (magenta dashed lines) are also shown. From Ref. [53]. (Figures adapted from M. Giordano, J. High Energy Phys. 05, 204 (2019), and used under a CC-BY 4.0 license (https://creativecommons.org/licenses/by/4.0)).

Localization of Dirac modes was studied in Z<sup>2</sup> pure gauge theory in 2+1 dimensions in Ref. [56], probed with unimproved staggered fermions. This model has the simplest gauge group, and the lowest dimensionality in which a deconfining transition is found. Studying the fractal dimension *D*2, it was shown that low modes are localized (*D*<sup>2</sup> = 0) in the high-temperature, deconfined phase of the theory in the positive center sector (i.e., positive spatially averaged Polyakov loop), while they are delocalized (with *D*<sup>2</sup> < 2) in the low-temperature, confined phase, and in the high temperature phase in the negative center sector (i.e., negative spatially averaged Polyakov loop). Localized modes are also found at the high end of the spectrum, independently of the phase and of the center sector. Significant correlation between localized modes and both Polyakov loops and clusters of negative plaquettes was observed.

While a genuine phase transition is expected for SU(3) gauge group in the presence of *Nf* = 3, light enough dynamical fermions [80], so far a critical point has been observed only on coarse lattices, and disappears in the continuum limit [228–230]. Although only a toy model for QCD, the SU(3) theory with *Nf* = 3 flavors of unimproved staggered fermions on *Nt* = 4 lattices is nonetheless a well-defined statistical model with a genuine first order transition, affecting both its chiral and confining properties, despite the absence of exact chiral and center symmetries. More precisely, as the coupling *β* crosses the critical value *βc*, the chiral condensate jumps downwards to a much smaller but still finite value; and the average Polyakov loop jumps upwards from its small but nonzero value to a considerably larger value. Evidence of localization of the low staggered Dirac modes was reported in Ref. [50] for bare fermion mass *m* = 0.01, below the critical value *mc* = 0.0259 [230], where genuine first-order phase transitions are present. A mobility edge was shown to be present for *β* > *βc*: it increases with *β*, and extrapolates to zero close to *βc*. The lowest mode was also seen to turn from delocalized to localized at a coupling *β*loc compatible with *βc*, i.e., in correspondence with the finite-temperature transition (see Figure 17, left).

**Figure 17. Left**: average Polyakov loop (red squares), chiral condensate (upward blue triangles) and average PR of the lowest staggered mode (downward magenta triangle) for *Nf* = 3 unimproved staggered fermions of bare mass *m* = 0.01 on *Nt* = 4 lattices. The critical coupling *β<sup>c</sup>* = 5.0985 is also shown. From Ref. [50]. (Figure adapted from M. Giordano, S. D. Katz, T. G. Kovács, and F. Pittler, J. High Energy Phys. 02, 055 (2017), and used under a CC-BY 4.0 license (https://creativecommons.org/licenses/by/4.0)). **Right**: mobility edge as a function of the deformation parameter *h* in trace-deformed SU(3) gauge theory in the high-temperature deconfined phase (*β* = 6.0). Here the critical deformation parameter for reconfinement is *hc* = 0.1. The black diamond is the *h* = 0 result of Ref. [51]. From Ref. [55]. (Figure adapted from C. Bonati, M. Cardinali, M. D'Elia, M. Giordano, and F. Mazziotti, Phys. Rev. D 103, 034506 (2021), and used under a CC-BY 4.0 license (https://creativecommons.org/licenses/by/4.0)).

The relation between localization and deconfinement was tested at a different deconfinement phase transition in trace-deformed [231,232] pure gauge SU(3) theory at finite temperature in Ref. [55]. In this model a deformation term Δ*S* = *h* ∑*<sup>x</sup>* |tr*P*(*x*)| <sup>2</sup> is added to the action, which (for *h* > 0) tends to locally suppress a nonzero trace for the Polyakov loop *P*(*x*). For temperatures above the deconfinement temperature, Δ*S* pushes the theory towards a "reconfined" phase where tr*P*(*x*) ∼ 0. This happens when the deformation parameter *h* crosses a (temperature dependent) critical value *hc*. Ref. [55] studied the spectrum of the two-stout smeared staggered spectrum at *β* = 6.0 on *Nt* = 6 lattices for various volumes and deformation parameters. Results showed that localized modes are present for *h* < *hc*, but disappear as the system crosses over into the reconfined phase. The mobility edge was determined by comparing the fractal dimension of the modes with its value at criticality [148]. While monotonically decreasing with *h*, it is not clear whether it vanishes continuously at *hc* or jumps to zero discontinuously (see Figure 17, right).

Finally, the connection between localization and ordering of the background configuration was studied in spin models in Refs. [47–49,52]. Ref. [47] used a simple 3d Hamiltonian in the orthogonal class with on-site disorder provided by the spins of a continuous-spin Ising-type model. In the ordered phase of the spin model, localization was observed for the low modes, with a mobility edge separating them from higher modes, and critical behavior compatible with that of the 3d orthogonal Anderson model. Refs. [48,49] dealt with the Dirac-Anderson form of the staggered operator (see Section 5.2), so in the 3d chiral unitary class, in the case *Nt* = 2 in the background of Polyakov loops constructed from a spin model. Localized low modes are observed in the ordered phase of the model [48], appearing at the critical temperature [49]. Ref. [52] reports on the *CP*<sup>3</sup> model in 1+1 and 2+1 dimensions. While in 1+1 dimensions localized modes are found in both phases of the model, as expected in one spatial dimension, in the 2+1 case localized modes are found only in the ordered phase. This model belongs to the 2d chiral unitary class.

#### **7. Conclusions and Outlook**

The presence of localized modes in the spectrum of the Dirac operator in the hightemperature phase of gauge theories is by now well established. Numerical studies on the lattice have shown that above the transition temperature the low-lying Dirac modes

are spatially localized on the scale of the inverse temperature, in QCD and QCD-like theories, as well as in several pure gauge theories and related models in 3+1 and 2+1 dimensions [25–40,47–56,86].

The physical significance of localization has so far remained quite elusive. First of all we should emphasize an important difference between localization in electron systems and localization in QCD. In the former case the mobility edge in the spectrum can be "accessed" by tuning a suitable control parameter, such as an electric field or the density of electrons. As the Fermi energy crosses the mobility edge, the system undergoes a genuine phase transition, with the zero-temperature conductivity changing non-analytically. In contrast, the mobility edge in the QCD Dirac spectrum cannot be directly connected to a thermodynamic transition. This is because in that case, in general, there is no control parameter that can be adjusted to make the system sensitive to just the eigenmodes at the mobility edge in the spectrum. The only exception is when the mobility edge is at zero, which happens only at the critical temperature of localization. If at the same time the quark masses are set to zero, the system becomes most sensitive to the lowest Dirac eigenmodes, the ones closest to zero. Thus, only in this double limit when the temperature tends to the critical temperature of localization and the quark mass to zero can one possibly directly connect the localization transition to a genuine thermodynamic phase transition. Unfortunately, this limit is out of the reach of present day lattice simulations and we have no numerical evidence of what happens there.

On the other hand, some progress has been made to understand the physical significance of localization in QCD. A clear connection with deconfinement has emerged: in all the models investigated so far, localization of the low modes shows up when the system transitions from the confined, low-temperature phase to the deconfined, high-temperature phase [28,33,49–51,53–56]. Convincing evidence has been presented for the crucial role played by the ordering of the Polyakov loop and by its fluctuations in the formation of a mobility edge in the Dirac spectrum, separating low-lying, localized modes from the delocalized bulk modes [31,39,40,48,56]. As the Polyakov loop is the (approximate, in the case of QCD) order parameter for confinement, the observed connection between localization and deconfinement has a dynamical explanation, further backed by a viable mechanism (the sea/islands picture [31,47–49], see Section 5.2) relating the two phenomena. This raises the hope that further studies can lead to a better understanding of confinement, and possibly to the uncovering of the mechanism behind this remarkable property of gauge theories. In this context, it would be interesting to further elucidate the relation between localization and center symmetry, since so far only models with nontrivial gauge group center have been investigated.

Localization could also help in explaining the close relation observed between deconfinement and restoration of chiral symmetry. These two phenomena in fact take place at the same temperature, or in a relatively narrow interval of temperatures, where also localized low Dirac modes appear. Localization could then provide the key to understanding this relation between in principle unrelated phenomena. Unfortunately, while the connection between deconfinement and localization can be easily studied in a clear-cut situation by investigating pure gauge theories with a genuine deconfining phase transition, studying the connection between chiral symmetry restoration and localization by means of numerical lattice simulations faces the considerable difficulties involved in taking the chiral limit. Studies of this type would be of great interest, especially in the light of the possible role played by localized modes in suppressing the finite-temperature Goldstone excitations, suggested in Ref. [62]. A particularly interesting case would be that of adjoint massless fermions, for which both chiral and center symmetries are exact, and an intermediate, deconfined but chirally broken phase was observed on the lattice for two flavors in Ref. [233]. This suggests that a nonzero density of near-zero, localized modes is present in this phase, and that no Goldstone excitation is present.

Nonetheless, even in theories such as QCD where chiral symmetry is only approximate, the study of the relation between localized modes and topological fluctuations of the gauge

fields sheds indirectly some light on the interplay of chiral symmetry and localization. Indeed, a peak of near-zero [59–61] localized [38] modes of topological origin appears around the QCD pseudocritical temperature. These modes originate most likely from the mixing of the localized zero modes associated with isolated instantons and anti-instantons, which in the high-temperature phase form a dilute gas of topological excitations (the disordered medium scenario [8–15], see Section 5.1). In contrast, in the low temperature phase these excitations form a dense medium, and the mixing of the associated zero modes leads to a band of near-zero delocalized modes giving rise to a nonzero spectral density near the origin, and so a large increase of the chiral condensate.

An interesting observation is that in the quenched limit of QCD this peak of nearzero modes can be accurately described in terms of a non-interacting gas of topological objects [58,86]. The absence (or near absence) of interactions could be related to why the associated zero modes do not mix efficiently, thus remaining localized and failing to spread in a near-zero band of eigenvalues. On the other hand, as this occurs only in the high-temperature phase, it is natural to expect that deconfinement is responsible for the radical change in the behavior of topological excitations. This aspect surely deserves more attention.

In this review, we summarized what is known (to us, at least) about localization of Dirac modes in the deconfined phase of gauge theories, and highlighted the connections between localized modes, ordering of the Polyakov loop, density of low modes, and topological objects. We hope that this will motivate further investigations of the interplay of confinement, chiral symmetry, topology, and localization in finite-temperature gauge theories.

**Author Contributions:** Writing—original draft preparation, M.G. and T.G.K; writing—review and editing, M.G. and T.G.K. Both authors have read and agreed to the published version of the manuscript.

**Funding:** M.G. was partially supported by the NKFIH grant KKP-126769.

**Acknowledgments:** We thank R. Shindou and M.R. Zirnbauer for useful correspondence.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Notes**


#### **References**


## *Article* **Influence of Fermions on Vortices in SU(2)-QCD**

**Zeinab Dehghan 1, Sedigheh Deldar 1, Manfried Faber 2,\*, Rudolf Golubich <sup>2</sup> and Roman Höllwieser <sup>3</sup>**


**Abstract:** Gauge fields control the dynamics of fermions, and, in addition, a back reaction of fermions on the gauge field is expected. This back reaction is investigated within the vortex picture of the QCD vacuum. We show that the center vortex model reproduces the string tension of the full theory also in the presence of fermionic fields.

**Keywords:** quantum chromodynamics; confinement; center vortex model; vacuum structure

**PACS:** 11.15.Ha; 12.38.Gc

#### **1. Introduction**

The QCD vacuum is highly nontrivial and has magnetic properties, as we have known since Savvidy's article [1]. The QCD vacuum should explain the non-perturbative properties of QCD, including confinement [2] and chiral symmetry breaking [3]. Lattice QCD puts the means at our disposal to answer the question about the important degrees of freedom of this non-perturbative vacuum. In the center vortex picture [4–6], the QCD vacuum is seen as a condensate of closed quantized magnetic flux tubes. These flux tubes have random shapes and evolve in time and therefore form closed surfaces in the dual space. They may expand and shrink, fuse and split and percolate in the confinement phase in all space-time directions and pierce Wilson loops randomly. Thus, Wilson loops asymptotically follow an exponential decay with the area. This is the area law of Wilson loops, which allows attributing the string tension to center vortices. The finite temperature phase transition is characterized by a loss of center symmetry and correspondingly by a loss of percolation in time direction. Therefore, vortices get static and only spatial Wilson loops keep showing the area law behavior.

Color electric charges are sources of electric flux according to Gauss's law. The electric flux between opposite color charges does not like to penetrate this magnetic "medium" of center vortices and shrinks to the well-known electric flux tube. On the other hand, the magnetic flux does not like to enter the electric string. Since fermions carry color charges, their dynamics is controlled by the gauge field. The presence of a fermion condensate is expected to suppress the quantized magnetic flux lines, and as a result the gluon condensate and therefore the string tension are reduced. Since, as usual, the lattice spacing is determined via the string tension, taking into account dynamical fermions leads to a decrease of the lattice spacing. In this article, we show a careful investigation of the string tension within the vortex picture of the QCD vacuum.

SU(2) and SU(3) QCD have equivalent non-perturbative properties. In a first study, we restrict our analysis to the simpler case of SU(2)-QCD. The most important difference between SU(2) and SU(3) QCD is the order of the finite temperature phase transition for a pure gluonic Lagrangian. There is a natural explanation for this difference from the structure of SU(2) and SU(3) vortices. There is only one non-trivial center element in the

**Citation:** Dehghan, Z.; Deldar, S.; Faber, M.; Golubich, R.; Höllwieser, R. Influence of Fermions on Vortices in SU(2)-QCD. *Universe* **2021**, *7*, 130. https://doi.org/10.3390/ universe7050130

Academic Editor: Dmitri Antonov

Received: 2 April 2021 Accepted: 1 May 2021 Published: 4 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

group SU(2) and therefore one type of center vortices, whereas there are two non-trivial center elements for SU(3) and two types of vortices, allowing two vortices of the same type to fuse to the other type. This leads to a more stable structure of the net of vortices for SU(3) and to a first order phase transition, whereas in SU(2) the transition is of second order.

We investigate the fermionic back reaction on the gluonic degrees of freedom in SU(2) QCD. Visualizing the distribution of center vortices, this back reaction can be easily observed (see Figure 1). One can clearly see that dynamical fermions decrease the percolation of vortices. It is difficult to draw a closed surface in four dimensions. Therefore, we restrict ourselves to the three dimensional diagram of a time-slice and indicate the continuation to other slices by line stubs.

**Figure 1.** The closed vortex surface is visualized by showing the dual P-plaquettes of threedimensional lattice slices. Stubs of red lines indicate plaquettes that are not fully part of the lattice slice shown. We clearly see that with fermions (**left**) an overall smaller amount of P-plaquettes is observed compared with the pure gluonic case (**right**). In both cases, one big vortex cluster dominates.

We want to quantify the effect in more detail. We are especially interested in the center vortex model [4–6] and its sensitivity to the fermionic back reaction. We also analyze the influence of fermionic fields on the geometric structure of the center vortex surface. This work compares four different estimates of the string tension, with and without fermions, in the full theory and in the vortex picture:


With this comparison, we study the sensitivity of the center vortex model to the fermionic back reaction.

Our work is based on the QCD path integral which defines the vacuum to vacuum transition amplitude. In lattice QCD, we usually evaluate this amplitude on a lattice periodic in Euclidean time. Inserting a complete set of eigenstates of QCD with the quantum numbers of the vacuum in this amplitude results in an exponential decay of the eigenstates with the physical time extent *aNt* of the lattice, where *a* is the lattice spacing and *Nt* is the number of lattice sites in time direction. The inverse of this time extent therefore acts as a temperature of the ensemble. In a Monte-Carlo simulation, the states are occupied with the corresponding Boltzmann factor. The higher is the excitation, the smaller is the Boltzmann factor and the more difficult is the measurement of its properties. Finally, the excited states are vanishing in the noise.

The potential, as well as the string tension, can be calculated using *Wilson loops*

$$\mathcal{W}(\mathbb{R}, T) = \text{Tr}\mathcal{P}\mathbf{e}^{\mathrm{i}\,\,\mathcal{S}} \, ^{\mathrm{i}}\!\!\!\!\!\!\!^{\mathrm{R}\times T} \, A^{a}\_{\mu}(\mathbf{x}) t\_{a} \mathrm{d}\mathbf{x}^{\mu} \, ^{\mu}\mathrm{.}\tag{1}$$

A loop of size *R* × *T* in space-time represents the world-line of a quark-antiquark system at distance *R* propagating in the QCD vacuum for a time *T*. On a Euclidean lattice in SU(2)-QCD, a path ordered loop is determined by the product of link variables *Uμ*(*x*) ∈ SU(2) along the loop. Inserting a complete set of eigenstates of the quark-antiquark system into the expectation value *W*(*R*, *T*) , the contributions of the eigenstates decay exponentially with Euclidean time *T*. The expectation values of Wilson loops can therefore be expanded in a series of eigenstates of the quark-antiquark system

$$
\langle \mathcal{W}(R, T) \rangle = \sum\_{i=0}^{\infty} c\_i \mathbf{e}^{-\varepsilon\_i(R)T} \, \, \, \, \, \tag{2}
$$

For large times, *W*(*R*, *T*) is dominated by the ground state energy *ε*0(*R*). The more precise we determine lim*T*→<sup>∞</sup> *W*(*R*, *T*) , the better is the precision of the quark-antiquark potential *V*(*R*) := *ε*0(*R*). Since the energy of the quark-antiquark system increases with the distance *R*, it follows from the above discussion that for increasing *R* the signal for *V*(*R*) is vanishing soon in the noise. How we handle this noise and how center vortices are detected is explained in Section 2. We assume that the potential is dominated by a Coulombic part at small *R* but rises linearly for large *R*,

$$
\varepsilon\_0(R) = V(R) = V\_0 + \sigma R - \frac{\alpha}{R}.\tag{3}
$$

We use *W*(*R*, *T*) to approximate *ε*0(*R*), denoted as *1-exp fit*. *V*<sup>0</sup> parameterizes the scale dependent self-energy of the quark-antiquark sources. Wilson loops extracted from the center degrees of freedom are dominated by the long-range fluctuations of the QCD vacuum, hence we describe the potential within these degrees of freedom by

$$V\_{\rm CP}(R) = \upsilon\_0 + \sigma\_{\rm CP} R.\tag{4}$$

The aim of this article is to investigate whether we can understand the string tension and its modification in the presence of fermions in the vortex model of confinement. Further, we present and discuss conceptual improvements to the gauge fixing procedure, required for the center vortex detection.

For systems with dynamical fermions one would expect string breaking when the energy of the system rises above twice the pion mass, but string breaking has been detected only using mesonic channels (see [7]). The center vortex model explains the asymptotic behavior of Wilson loops. There are indications that center vortices are sensitive to string breaking [8,9], but a direct measurement is not possible. From the vortex structure, we do not find any indication for string breaking which could show up as disintegration of the percolating vortex.

#### **2. Materials and Methods**

This section starts with a description of the parameters of the lattice configurations, used for our analysis. Then, our method of detecting center vortices with some novel improvements is discussed. We explain how the information about the geometric structure of the vortex surface can be acquired by smoothing procedures and we end with a detailed explanation of our method to extract the potential from Wilson loops. In each subsection, we list the intermediate results.

#### *2.1. Simulation Specifications*

We study the configurations described in [10] for chemical potential *μ* = 0 with *SG* defined by a tree level improved Symanzik gauge action [11,12]

$$S\_G = \beta \left( c\_0 \sum\_{\Box} (1 - \frac{1}{2} \text{Tr} \, \Box) + c\_1 \sum\_{\Box \Box} (1 - \frac{1}{2} \text{Tr} \, \Box \Box) \right), \tag{5}$$

with coefficients *c*<sup>0</sup> = 5/3 and *c*<sup>1</sup> = −1/12. The first sum corresponds to the Wilson action with indicating single unoriented plaquettes, while the second sum uses rectangular Wilson loops built of 6 links, symbolized by . The inverse coupling is defined as *β* = <sup>4</sup> *g*2 for SU(2).

For the fermionic degrees of freedom, staggered fermions are used with an action of the form

$$S\_F = \sum\_{x,y} \vec{\psi}\_X M(m)\_{x,y} \psi\_Y + \frac{\lambda}{2} \sum\_x \left( \psi\_x^T \tau\_2 \psi\_x + \vec{\psi}\_X \tau\_2 \vec{\psi}\_x^T \right), \tag{6}$$

with *τ<sup>i</sup>* being the Pauli matrices and

$$M(m)\_{xy} = m\delta\_{xy} + \frac{1}{2} \sum\_{\nu=1}^{4} \eta\_{\nu}(x) \left[ \mathcal{U}\_{\mathbf{x},\nu} \delta\_{\mathbf{x}+\mathbf{h}\_{\nu},y} - \mathcal{U}\_{\mathbf{x}-\mathbf{h}\_{\nu},y}^{\dagger} \delta\_{\mathbf{x}-\mathbf{h}\_{\nu},y} \right],\tag{7}$$

where *ψ*¯, *ψ* are staggered fermion fields, *a* is the lattice spacing, *m* is the bare quark mass, *Ux*,*<sup>ν</sup>* is a SU(2) element corresponding to a link at position *x* is in direction and *μ* and *ην*(*x*) are the standard staggered phase factors: *<sup>η</sup>*1(*x*) = 1, *ην*(*x*)=(−1)*x*1+...+*xν*−<sup>1</sup> , *<sup>ν</sup>* <sup>=</sup> 2, 3, 4. The total action is given by *S* = *SG* + *SF*. Integrating out the fermionic degrees of freedom, the partition function with *Nf* = 2 is given by

$$Z = \int \, D\!\!\!\!\!\!\!\!\!\!\!\!\!\!\/\, e^{-S\_G} \, (\det(\!M^\dagger M\!\!\/) + \lambda^2)^{\frac{1}{4}}.\tag{8}$$

The properties of 1000 configurations of size 324 with *β* = 1.8, quark mass parameter *m* = 0.0075 (corresponding to *m<sup>π</sup>* = 740(40) MeV with lattice spacing *a* = 0.044 fm), and *λ* = 0.00075 are compared to 1000 pure gluonic configurations at the same inverse coupling *β*. For both sets of 1000 configurations, we extract the potentials from all available Wilson loop data and compare them with the string tensions resulting from the two sets of 40 × 100 center projected configurations. In this way, we try to answer the question, if in the presence of dynamical fermions the center degrees of freedom determine the string tension of the gluonic flux tube in quark–antiquark systems.

#### *2.2. Center Vortex Detection*

Assuming that center excitations are the relevant degrees of freedom for confinement, we detect these center vortices within the lattice configurations. We first identify gauge matrices <sup>Ω</sup>(*x*) <sup>∈</sup> SU(2) at each site *<sup>x</sup><sup>μ</sup>* maximizing the functional

$$\mathcal{R}\_{\mathbb{F}} = \sum\_{\mathbf{x}} \sum\_{\mu} |\, \text{Tr}[\acute{l}\acute{l}\_{\mu}(\mathbf{x})] \mid \text{ $^2$ } \quad \text{with} \quad \acute{l}\acute{l}\_{\mu}(\mathbf{x}) = \Omega(\mathbf{x} + \boldsymbol{\varepsilon}\_{\mu}) \wrlearrowleft \Omega\_{\mu}(\mathbf{x}) \Omega^{\dagger}(\mathbf{x}). \tag{9}$$

After fixing the gauge, the link variables *U*´ *<sup>μ</sup>*(*x*) are projected on the center degrees of freedom, that is ±1 for SU(2), to neglect short range properties and keep only longrange effects

$$\mathcal{U}\_{\mu}(\mathbf{x}) \to Z\_{\mu}(\mathbf{x}) \equiv \text{signTr}[\mathcal{U}\_{\mu}(\mathbf{x})].\tag{10}$$

After performing the center projection, the center projected plaquettes resulting from the vortex detection are the products of four center elements. The projected plaquettes are non-trivial, known as P-plaquettes, *U* = −1, if one or three links are non-trivial. In the four-dimensional lattice, a given link belongs to six plaquettes. On the dual lattice, the corresponding six plaquettes build the surface of a cube. Therefore, the duals of Pplaquettes form closed surfaces, dual P-vortices which correspond to the closed flux line evolving in time.

This procedure is the original DMCG [13] in which a gradient climb with *over relaxation* was used to maximize the gauge functional. From a few gauge copies only, produced in this way, the one with the highest value of the functional is usually chosen for further analysis. This method leads to promising results, but improvements at maximizing the gauge functional using *simulated annealing* have brought a flaw to light—the many local maxima of *RF* do not necessarily correspond to the same physics. Bornyakov et al. [14] showed that there exist local maxima of the gauge functional that underestimate the string tension. We have been able to resolve these problems for smaller lattices using improved version of the gauge fixing routines based on non-trivial center regions [15–17], but our implementation was not able to handle the big lattices used in this work. Taking a closer look at the problem at hand, we can look for a different approach. We now consider Creutz ratios to estimate the string tension

$$\sigma \approx \chi(R) = -\ln \frac{\langle \mathcal{W}(R+1, R+1) \rangle \, \langle \mathcal{W}(R, R) \rangle}{\langle \mathcal{W}(R, R+1) \rangle \, \langle \mathcal{W}(R+1, R) \rangle^{\prime}} \tag{11}$$

with Wilson loops *W*(*R*, *R*) of size *R* = *T*. Some probability densities for the relation between the values of the gauge functional *RF* and the Creutz ratio *χ*(*R*) for individual configurations are shown in Figure 2. This determination is based on 40 configurations with 100 gauge copies for configurations with (left) and without (right) dynamical fermions. For Creutz ratios of small Wilson loops, we observe a nearly linear relation between the two quantities reflecting the finding of Bornyakov et al. [14]: there exist gauge copies of the configurations with maximal *RF* and very low *σ*. With increasing size of Wilson loops, this correlation weakens. Nevertheless, the request to maximize the gauge functional (9) fails.

Another observation is of high interest: extremely small and large values of the gauge functional are strongly suppressed in the probability densities. Instead of looking for higher local maxima of the gauge functional, we propose a different approach: "ensemble averaged maximal center gauge" (EaMCG). We produce many random gauge copies, approach the next local maximum by the gradient method and take the average of the ensemble. The idea is that not the best local maximum alone carries the physical meaning, but the average over all local maxima does: maxima with a higher value of the gauge functional result in a reduced string tension, but they are not dominating the ensemble. The same holds for lower valued maxima, possibly overestimating the string tension.

Taking again a look at Figure 2, it can be seen that the average values and the most probable values are in good agreement for small loops. This is shown in more detail in Figure 3 for Creutz ratios of different loop-sizes. The fact that differences increase with loop sizes can probably be explained by the lack of statistics for the Creutz ratios of single configurations. Until the values start to deviate from one another, there is a variation of 10% over the whole *R*-region. Despite the low statistics of a single configuration, the intermediate loop sizes already reproduce the asymptotic behavior and let us expect the possibility for a more precise determination. First averaging over Wilson loops and then calculating Creutz ratios gives much more stable results (see *χW*(*R*) in Figure 3). The final estimate of the string tensions in Sections 2.4 and 3 is based on the determination of the potential.

**Figure 2.** These probability densities specify the relation between the values of gauge functional and Creutz ratio for individual configurations. This determination is based on 40 configurations with 100 gauge copies for the configurations with dynamical fermions (**left**) and without (**right**). For Creutz ratios of small Wilson loops, we observe a nearly linear relation. With increasing size of the Wilson loops this correlation weakens. We marked the average values (star) and the most probable values (circle) of the distributions.

**Figure 3.** The average and most probable values of *χ*(*R*) are compared for simulations with and without fermions. This complements the probability densities of Figure 2. The increasing discrepancy between the two quantities with larger *R* can probably be explained by the low precision of Creutz ratios of single configurations. Until the values start to deviate from one another, the variations of *χ*(*R*) are of the order of 10%. For comparison, we show also the more precise Creutz ratios *χW*(*R*) extracted from averages of Wilson loops.

Thus far, we have calculated the Creutz ratios for single configurations of the ensemble and have taken the average afterwards. The EaMCG itself does not average over Creutz ratios, but combines first the Wilson loops of all gauge copies and configurations. From this fact, it is possible to extract the quark anti-quark potential, which allows a more precise determination of the string tension from the center vortex model.

In the respective single configurations, we observe one percolating large cluster that is surrounded and traversed by small fluctuations. These result in an increased number of P-plaquettes that do not contribute to the string tension. Analyzing these distortions, we gain insight on the influence of fermions on the geometric structure of the vortex surface.

#### *2.3. Smoothing the Vortex Surface*

There exist several procedures for smoothing the vortex surface by removing distortions. These procedures are discussed in detail in [18]. They do not modify the long range effects of the configuration. To get information about the smoothness of the vortex surface with and without fermions, we use the smoothing steps depicted in Figure 4. The smoothing 0 is not depicted, which removes unit-cubes.

**Figure 4.** The effect of the smoothing procedures on the vortex surface is depicted, taken from ([19] Figure 5.8). We distinguish *warts* (**left**), *bottlenecks* (**middle**) and *stumbling blocks* (**right**). The *unit cubes* are not depicted, which are simply deleted.

The smoothing steps 1–3 cut out parts of the vortex surface and closes the emerging holes with a flat surface. In this way, short-range fluctuations of the vortex surface are suppressed. We first count the P-plaquettes without any smoothing performed, and then the loss of P-plaquettes for the respective smoothing steps is determined. The results are given in Table 1.


**Table 1.** Reduction of the total count of P-plaquettes for different smoothing procedures.

This quantifies the percentage of the respective structures depicted in Figure 4. When fermions are present, we clearly have a higher proportion of unit cubes and a lower proportion of bottlenecks than without fermions.

By restricting this analysis to the single percolating vortex cluster, we gain information about the long range excitations. The results are given in Table 2. The reduction in the proportion of bottlenecks is also seen here. The presence of fermions leads to a smoother surface of the percolating cluster.

**Table 2.** Reduction of P-plaquettes for the percolating vortex cluster for different smoothing procedures.


*2.4. Potential Fits and Noise Handling*

When extracting the potential from Wilson loops, two effects have to be taken care of:


An example for a 1-exponential fit to Wilson loops *W*(*R*, *T*) for given *R* and *T* ≥ *Ti* (see Equation (2)) is shown in the left diagram of Figure 5. The dependence of this example on the initial *T* = *Ti* is depicted in the right diagram.

**Figure 5.** (**Left**) Example of the optimal 1-exponential fit of Wilson loops for given *R*. (**Right**) Dependence of *ε*0(*R*, *Ti*) on the fit region *T* ≥ *Ti*. The line marks the fit for the optimal value for *Ti*.

At lower *Ti*, an increase of *Ti* causes large changes of the fit parameters, but with growing *Ti* these changes become smaller until a most stationary point is reached which may be hidden behind a strong increase of error bars. With the naked eye, one sees data that result in quite good fits, but finding analytic or numeric criteria for the choice of *Ti* proves difficult. The smaller is the change of the values of the fit parameters, the smaller are the error bars of Wilson loops, and a rapid increase of the *p*-value of the fits often coincide, but this is not a general rule. Our criteria to choose *Ti* is based on identifying the first local minimum of an error quantifier

$$\text{Err} := \frac{2}{3} \langle \Lambda\_{\delta i} \rangle + \frac{1}{3} \langle \Lambda\_{\text{err}} \rangle. \tag{12}$$

Here, Δ*δ<sup>i</sup>* denotes the average change of the fit parameter *<sup>ε</sup>*(*R*, *Ti*±1) when decreasing or increasing *Ti*; and Δerr denotes the average over the error bars of *<sup>ε</sup>*(*R*, *Ti*−1), *<sup>ε</sup>*(*R*, *Ti*), and *ε*(*R*, *Ti*+1). The weight factors are chosen to avoid the choice of occasionally nearly stationary regions with large error bars. For *R* > 3, we prevent any further increase of *Ti*, because with increasing *R* the error bars start to grow earlier. The example in Figure 5 tries to convince that the selection of *Ti* based on the error quantifier results in optimal fits under the boundary conditions of systematic deviations for low *Ti* and increasing error bars for high *Ti*. Using this procedure, we determine the potential for the whole range of *R*-values, which allows extracting the slope of the potential at large values of *R*.

#### **3. Summarized Results and Discussion**

The fermionic back reaction on the string tension is clearly observed in the full theory as well as for EaMCG (Ensemble averaged Maximal Center Gauge), where the link variables of the gauge field are projected to *Z*2. The potentials for the gluonic and fermionic configurations are depicted in Figure 6 and compare the full SU(2) theory with the *Z*<sup>2</sup> theory. The string tension was extracted by fitting the respective Equation (3) or (4) to the data describing the potential. The resulting parameters of these fits are given in Table 3. The relevant parameters to compare are *σ* and *σ*CP.

**Figure 6.** (**Left**) Potential *V*(*R*) in lattice units between two sources in the fundamental representation. There is a large difference between the string tensions for pure gluonic configurations ("gluonic") and in the presence of one species of dynamical fermions. (**Right**) Potentials extracted from Wilson loops after ensemble averaged maximal center projection are depicted for pure gluonic configurations and for configurations with dynamical fermions. Due to the removal of short range fluctuations the potentials are in both cases almost linearly increasing with the lattice distance *R*. Data are fitted by linear functions. For gluonic (fermionic) configurations, only data with *R* ≥ 6(2) are fitted.

Without fermions, both estimates for *σ* are compatible within errors to one another: In the full *SU*(2), we observe *σ* = 0.0756(12) compared to *σ*CP = 0.07691(13) in the Z2 description. With fermions the full SU(2) theory results with *σ* = 0.0199(9) a lower value than the Z2 theory with *σ*CP = 0.02291(5). In all cases, we clearly observe that the presence of fermions reduces the string tension in lattice units: The back reaction is observed in the full SU(2) theory and also reproduced by the center vortices. The determination of the lattice spacing via the usual formula (*a* = √*χ*/2.23 fm, corresponding

to *χ* = (440 MeV)2) results in 0.123(1) fm for the gluonic configurations and 0.0633(15) fm for the fermionic configurations.

**Table 3.** The parameters for the fits according Equations (3) and (4) in Figure 6 allow a direct comparison of the respective string tensions. A strong suppression of the Coulomb part can be seen in the *Z*<sup>2</sup> theory.


Concerning the geometric structure of the vortex surface, we observe that the presence of fermions increases the number of isolated short range fluctuations (see Table 1): without fermions, about 6.98% of the P-plaquettes are part of isolated unit cubes, whereas, with fermions, this proportion increases to 12.45%. The proportion of P-plaquettes belonging to bottlenecks is in total decreased from 27.81% to 24%. Fermions increase the amount of unit cubes, but decrease the amount of bottlenecks.

Restricting the analysis to the long-ranged cluster we observe a decrease of fluctuations, especially bottlenecks, when fermions are present (see Table 2): the proportion of P-plaquettes belonging to bottlenecks is reduced from 28.12% to 24.45%. All other fluctuations are only reduced by about 1%.

From this, we can conclude that the presence of fermions causes short range fluctuations to detach from the vortex surface, resulting in a more smooth vortex surface that is surrounded by an increased number of isolated short range fluctuations.

**Author Contributions:** Conceptualization, M.F. and R.H.; methodology, Z.D., M.F., R.H. and R.G.; validation, S.D.; investigation, Z.D. and R.G.; data curation, Z.D. and M.F.; writing—original draft preparation, Z.D. and R.G.; writing—review and editing, Z.D. and S.D.; visualization, R.G. and R.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** We thank Aleksandr Nikolaev, Nikita Astrakhantsev and Andrey Kotov for their cooperation in the early stage of this investigation and Vitaly Bornyakov for important advice.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

