An Ising-Anderson model of localisation in high-temperature QCD

We discuss a possible mechanism leading to localisation of the low-lying Dirac eigenmodes in high-temperature lattice QCD, based on the spatial fluctuations of the local Polyakov lines in the partially ordered configurations above $T_c$. This mechanism provides a qualitative explanation of the dependence of localisation on the temperature and on the lattice spacing, and also of the phase diagram of QCD with an imaginary chemical potential. To test the viability of this mechanism we propose a three-dimensional effective, Anderson-like model, mimicking the effect of the Polyakov lines on the quarks. The diagonal, on-site disorder is governed by a three-dimensional Ising-like spin model with continuous spins. Our numerical results show that localised modes are indeed present in the ordered phase of the Ising model, thus supporting the proposed mechanism for localisation in QCD.


Introduction
The spectrum of the lattice Dirac operator plays a prominent role in current attempts to improve our understanding of the spontaneous breaking of chiral symmetry in QCD. The key relation in this context is the celebrated Banks-Casher formula [1], which clarifies the relevance of the low-lying Dirac eigenmodes for the development of a non-vanishing chiral condensate in the chiral limit. The Banks-Casher relation also suggests that the (pseudo)critical behaviour of the theory at the chiral transition/cross-over, even away from the chiral limit, will be mostly determined by the behaviour of these modes.
It has long been known that in QCD, below the pseudocritical temperature, T c , the low-lying eigenmodes of the Dirac operator are delocalised, i.e., they extend throughout the whole four-dimensional volume of the system [2] (see also Ref. [3]). In recent years it has been observed that above T c this is no longer true: the lowest-lying modes are localised on the scale of the inverse temperature [4,5]. Delocalised modes are still present, but only above a (temperature-dependent) critical point λ c = λ c (T ) in the spectrum. A similar change in the properties of the low-lying modes is found also in other, QCD-like theories, e.g., with SU (2) gauge group, and/or in the quenched limit [6][7][8][9].
As the temperature is decreased, λ c moves towards the origin, and vanishes at a temperature compatible with T c , the cross-over temperature determined from thermodynamic observables [10,11]. This is certainly not a coincidence. First of all, in QCD [4,7] and in all the QCD-like theories where localised modes have been observed (pure-gauge SU (2) [8,9] and SU (3) [7] theories), they are present only in the "chirally-symmetric phase", independently of the specific lattice discretisation of the Dirac operator (staggered [7,9], overlap [8], and recently also Möbius domain wall [12]). Moreover, further evidence of a very close relationship between the appearence of localised modes and the chiral transition has been found recently [13] in a toy model for QCD, namely N T = 4 unimproved staggered fermions [14], which displays a genuine first-order chiral transition accompanied by the appearence of localised modes at the low end of the spectrum. However, the precise nature of this relationship has not been fully understood yet.
At fixed temperature, the transition in the spectrum from localised to delocalised modes has been shown to be a true, second-order Anderson-type transition [15]. Such transitions have been most extensively studied before using the Anderson model [16][17][18], which is a model for electrons in a crystal with disorder. The Hamiltonian of the Anderson model consists of the usual tight-binding Hamiltonian plus an on-site (diagonal) random potential mimicking the presence of impurities in the crystal. In three dimensions, when the diagonal disorder is switched on, localised modes appear at the band edges, beyond critical energies called "mobility edges". As the width of the distribution of the diagonal disorder is increased, i.e., as the system is made more disordered, the mobility edges move towards the band center, and beyond a critical value for the disorder all the modes become localised.
The model described above is the so-called orthogonal Anderson model. A variant of this model is the unitary Anderson model, in which the hopping terms are modified with the introduction of random phases to account for the presence of a magnetic field. The specifiers "orthogonal" and "unitary" refer to the symmetry class of the model in the language of random matrix theory (RMT) [19]. In this classification, QCD belongs to the unitary class. RMT makes universal predictions for certain statistical properties of the spectrum of a random Hamiltonian with i.i.d. entries, depending on its symmetry class. These predictions correctly describe also the spectral region corresponding to the delocalised modes of the Anderson Hamiltonian, despite this being a sparse matrix. However, this can be understood by noticing that delocalised eigenmodes are easily mixed by fluctuations in the disorder. On the other hand, localised modes are sensitive only to fluctuations taking place within their support, and thus they hardly mix. As a consequence, the corresponding eigenvalues are expected to fluctuate independently, thus obeying Poisson statistics. This behaviour is actually observed in the localised part of the spectrum of the Anderson model. Therefore, a change in the spectral statistics takes place in correspondence to the localisation/delocalisation transition.
The staggered Dirac operator in lattice QCD, being anti-Hermitian, admits a straightforward interpretation as a random Hamiltonian. As in the Anderson model, the statistical properties of the spectrum change from Poisson-type in the spectral region of localised modes, to the appropriate RMT-type in the region of delocalised modes. Besides this, however, at first sight one can find very little in common with the three-dimensional Anderson model. For one thing, the dimensionalities apparently do not match. Moreover, the Dirac operator has only off-diagonal non vanishing terms, i.e., the fluctuating hopping terms representing the gauge fields. In this case the location of the mobility edge was found to be controlled by the physical temperature. Finally, unlike in the Anderson model, the disorder in the gauge links is correlated. It is therefore surprising that the critical exponent of the correlation length has been found to be compatible with that of the unitary Anderson model in three dimensions [20], suggesting that the two models belong to the same universality class [15].
The unitary class is indeed the symmetry class appropriate for QCD, and correlations in the disorder should not affect the critical behaviour since they are short-ranged, as the disorder consists of gauge-field fluctuations. On the other hand, the fact that a four-dimensional model with only off-diagonal disorder behaves like a three-dimensional one with diagonal disorder calls for an explanation. 1 An important step towards clarifying these issues would be the determination of the mechanism responsible for the localisation of the low-lying modes. It was initially suggested [6] that instantons could play the role of localising "defects", but numerical investigations by one of the authors and collaborators [23] found a large mismatch between the density of localised modes and the density of instantons, thus making this explanation unlikely. In the same paper, an alternative proposal was made to explain localisation at high temperature, which involves the behaviour of local fluctuations of Polyakov loops above T c . The basic idea is that, in the ordered phase, local fluctuations of the Polyakov loop provide a localising "trap" for eigenmodes. This is because the "islands" where the Polyakov loop fluctuates away from the ordered value provide "energetically" favourable spatial regions for the quark eigenmodes to live in. Numerical support to this proposal was also obtained, in the case of the SU (2) gauge group.
The purpose of this paper is to refine this proposal, and extend it to the SU (3) case. This simple generalisation has nevertheless the merit of clarifying some important aspects that went previously unnoticed. To further test the viability of our explanation, we construct a simple three-dimensional model which should display localisation precisely through the mechanism proposed for QCD. Our numerical results show that this is indeed the case, and that the qualitative features of the eigenmodes and of the spectrum correspond to those of the eigenmodes and of the spectrum of the Dirac operator. In particular, our toy model displays a transition in the spectrum from localised to delocalised modes, which is likely to become a true phase transition in the thermodynamic limit.
The plan of the paper is the following. In Section 2 we discuss the mechanism for localisation in high-temperature QCD based on fluctuations of Polyakov lines. In Section 3 we construct an effective model ("Ising-Anderson" model) which should produce localised states precisely through the mechanism proposed in Section 2. In Section 4 we report on the results of numerical simulations of the Ising-Anderson model. Finally, in Section 5 we state our conclusions and discuss possible extensions of the present study.

Polyakov lines and localisation
In this Section we argue that high-temperature QCD contains an effectively diagonal and three-dimensional disorder, which explains why the critical behaviour at the localisation/delocalisation transition in the spectrum corresponds to that of the three-dimensional unitary Anderson model. An intuitive argument for the effectively reduced dimensionality of high-temperature QCD is the following: as the size of the temporal direction is smaller than the correlation length, the time slices are strongly correlated and the system is effectively three-dimensional. This implies that the Dirac eigenmodes should look qualitatively the same on all time slices, in particular for what concerns their localisation properties.
This intuition is strengthened when we consider the role of the antiperiodic boundary conditions in the temporal direction. To this end, it is convenient to work in the temporal gauge, U 4 (t, x) = 1 for t = 0, . . . , N T − 2 and ∀ x, in which U 4 (N T − 1, x) equals the local Polyakov line P ( x) ≡ N T −1 t=0 U 4 (t, x). A further time-independent gauge transformation allows one to diagonalise each local Polyakov line: we will refer to this as the diagonal temporal gauge. In any temporal gauge, covariant time differences are replaced by ordinary differences for all t = N T − 1. One can do the same also at t = N T − 1, if at the same time one trades the antiperiodic boundary conditions for effective, x-dependent boundary conditions, which involve the local Polyakov line, (2.1) Since the time slices are strongly correlated, these effective, x-dependent boundary conditions will affect the behaviour at the spatial point x for all times t. Furthermore, P ( x) fluctuates in space (and obviously from one configuration to another). From the point of view of a disordered system, QCD above T c therefore contains effectively a diagonal (on-site), threedimensional source of disorder. To see how these effective boundary conditions affect the quark wave functions, let us first discuss a simplified setting in which the Dirac equation can be explicitly solved, generalising the argument of Ref. [23] to SU (3). We discuss here the naïve lattice Dirac operator for simplicity, but the considerations of the present Section clearly apply to the staggered Dirac operator as well, since the two are essentially connected by a unitary transformation. Consider configurations with constant temporal links U 4 (t, x) = U and trivial spatial links U j (t, x) = 1, j = 1, 2, 3. In the diagonal temporal gauge, one has P ( x) = P = U N T = diag(e iφ 1 , e iφ 2 , e iφ 3 ) with 3 c=1 e iφc = 1, so the Dirac operator is diagonal in colour, and the colour components of the quark eigenfunctions decouple. The eigenfunctions of the naïve lattice Dirac operator have the following simple factorised form, Here and throughout the paper the eigenvalues are expressed in lattice units. The spacetime dependence is fully contained in the plane waves f (ck p) (t, x), while the spin and colour dependence are encoded in the bispinors χ where a is the lattice spacing and N T the temporal size of the lattice, to fulfill the effective temporal boundary conditions. We will refer to ω (ck) as the effective Matsubara frequencies.
For a given value of the phase φ c , there are N T different branches for ω (ck) , corresponding, however, to only N T /2 different eigenvalues, since the effective Matsubara frequencies obey the relation Notice that the lowest positive branch of eigenvalues is described by the function which decreases as φ c moves away from zero, where it is maximal, and it is minimal when e iφc = −1. The bispinors χ (ck ps±) α , s = 1, 2, satisfy the equation where γ µ , µ = 1, . . . , 4 are Euclidean Dirac matrices. The spin index s labels eigenmodes corresponding to the same eigenvalue ±iλ (ck p) . Taking into account the degeneracy with respect to k mentioned above, the eigenvalues are therefore fourfold degenerate. Finally, the colour vectors are trivially ϕ (c) j = δ cj , c, j = 1, 2, 3, and satisfy P ϕ (c) = e iφc ϕ (c) . Let us now consider the qualitative features of the eigenmodes in the general case. Above T c , a typical gauge configuration consists of a "sea" where the Polyakov line P ( x) gets ordered around 1 = diag(1, 1, 1), which percolates through the lattice, with "islands" of "wrong" P ( x) = 1. If P ( x) were completely ordered, and spatial links were trivial, there would be a sharp gap in the spectrum at λ g = M (0) = sin(πaT ) (and a symmetric one at −λ g ), corresponding to the lowest positive branch at φ c = 0 in Eq. (2.2) (see also Eq. (2.5)). This is a "zeroth-order" picture, which is modified by the fluctuations of P (x) and of the spatial links. The "first-order" picture is obtained by allowing P (x) to fluctuate while still keeping the spatial links trivial. It is clear from Eq. (2.5) that the quark eigenfunctions can exploit the "islands" of "wrong" P (x), which are "energetically" favourable, to lower their eigenvalues. In particular, localising the wavefunction on the "islands" can achieve a large eigenvalue reduction, if the momentum required to localise the state is not too large, which is the case, for example, if the "islands" are sufficiently big. 2 Therefore, in the "first-order" picture both 2 This can be seen by considering a configuration with uniform "islands" of "wrong" Polyakov lines, and a test wavefunction of the form given in Eq. (2.2) but localised and constant on one of such "islands". More precisely, take ψtest(t, x) = ψ (ck 0s+) (t, x)C isl ( x), with k chosen so that λ (ck 0) belongs to the lowest positive branch. Here C isl ( x) is 1 on the given "island" and 0 elsewhere. Computing the expectation value of − / D 2 , one finds a value smaller than λ 2 g if the "island" is big enough, so that surface effects at the interface with the "sea" do not overbalance the negative difference λ (ck 0)2 − λ 2 g .
localised and delocalised eigenmodes appear between −λ g and λ g , with a few localised low modes well separated from the bulk of delocalised modes. The full, "second-order" picture is finally obtained by switching on the fluctuations of the spatial links. These fluctuations have most likely a delocalising effect: for example, they allow different colour components to mix (recall that we are working in the diagonal temporal gauge). Nevertheless, the lowest modes can still remain localised, due to the large energy difference with the bulk of delocalised states. The sharp gap of the "zeroth-order" picture has turned into an "effective gap" λ c < λ g , which we identify with the mobility edge separating localised and delocalised modes.
To better understand this picture, it is useful to set up an analogy between the Dirac operator and the Hamiltonian of an electron in a crystal in the tight-binding approximation (TBA). If one discards the spatial links, setting them to zero, the solutions of the Dirac equation are localised at a given spatial site, and are given by Eq. (2.2) with p = 0, using the local Polyakov line as the gauge background. These "free" solutions correspond to the atomic orbitals in the TBA, which are localised on a crystal site. In turn, the effect of the spatial links corresponds to that of the hopping terms in the TBA, which allow the electron to hop between crystal sites. In this analogy, the "islands" of "wrong" Polyakov lines correspond to "defects" in the crystal (i.e., the diagonal disorder in the Anderson model), and they act in the same way as localising "traps" for the eigenmodes. Given the three-dimensional nature of these "islands", the qualitative picture described above allows one to understand why, in the high-temperature regime, the critical properties of the Dirac spectrum at the localisation/delocalisation transition correspond to those of the three-dimensional unitary Anderson model. In order to check the viability of our explanation, in the following Section we study an effective model which should display localisation precisely through the mechanism that we have proposed. In the remaining part of this Section we add a few comments.
The dependence of λ c on the temperature and on the lattice spacing is expected to correspond qualitatively to that of λ g , which should somehow "drag" the mobility edge along. In the "zeroth-order" picture, λ g = sin π N T = sin(πaT ) increases when increasing the physical temperature at fixed a, and it decreases as the lattice spacing is decreased at fixed physical temperature. This behaviour is qualitatively the same as that of the Polyakov-loop expectation value, which is a measure of the ordering of the Polyakov lines in the gauge configurations. Therefore, even though λ g is determined assuming a perfectly ordered system, we can loosely say that it responds in a direct way to the ordering of the Polyakov lines. In the "first-order" picture, the separation point λ ′ g between localised and delocalised modes lies below λ g . For a delocalised mode with low spatial momentum, which sees all the fluctuations of the Polyakov lines equally, the gain in "energy" will average out, so bringing λ ′ g near the average M (φ c ) of the lowest positive branch. In the "second-order" picture, λ ′ g will be displaced further down to λ c by the mode mixing, but the size of this displacement depends on the effectiveness of the mixing, which is difficult to assess. It is likely, however, that an increase in the ordering of the Polyakov lines will also make more correlated the spatial links at a given spatial site but on different time slices. The fluctuations of the spatial hopping should therefore become less pronounced, which in turn should reduce the effectiveness of the mixing, and so reduce the distance between λ c and λ ′ g . In this case, both the "first-order" and the "second-order" effects tend to reduce the distance between λ c and λ g as the system is made more ordered, and to increase it when the ordering is reduced, thus acting in the same direction as the "zeroth-order" effect. The bottom line is that the response of λ c to a change in the ordering of the system is expected to be qualitatively the same as that of λ g , i.e., it goes up with T at fixed a, and it goes down with a at fixed T . This expectation is matched by the results of numerical simulations [4].
The qualitative correctness of the "zeroth-order" picture allows also a qualitative understanding of the dependence of the (pseudo)critical temperature T c on an imaginary chemical potential. As we have said above, in the high-temperature phase where the Polyakov line gets ordered, the position of the mobility edge, λ c , is affected by the effective Matsubara frequency corresponding to the trivial Polyakov line. At φ c = 0, the "twist" of the quark wavefunction required by the effective boundary conditions is maximal, and so the quantity M (φ c ), defined in Eq. (2.5), which "drags" λ c , is maximal. Introducing an imaginary chemical potential µ (small enough in order to avoid the Roberge-Weiss critical line [24]) is equivalent to adding an extra phase e iµN T to the effective boundary conditions, which reduces the "twist". As a consequence, M (φ c ) → M (φ c + µN T ), and in particular M (0) → M (µN T ), i.e., it diminishes (independently of the sign of µ), and so one expects λ c to become smaller. 3 If the relation between the chiral transition and the appearence of localised modes holds true, then one expects T c to increase with |µ|. This is actually what has been observed in lattice simulations [25,26]. The presence of the Roberge-Weiss critical lines can also be understood in this framework, noticing that the system prefers to have an effective gap as large as possible due to the fermionic determinant. Comparing M ( 2π 3 k + µN T ) for k = −1, 0, 1, corresponding to the Polyakov line sectors P = e i 2π 3 k , one finds that it becomes "energetically" favourable to switch to the sector k = −1 when µN T exceeds π 3 , and similarly to the k = 1 sector, and again to the trivial sector, when µN T exceeds π and 5π 3 , respectively. As one would expect, our qualitative picture cannot explain all the features of the localised modes, especially when there are competing effects at play. As it has been shown in Ref. [4], the lowest-lying Dirac eigenmodes should remain localised also in the continuum limit. While λ c goes to zero as a → 0, in agreement with our qualitative explanation, the "renormalised mobility edge", λ c /m ud , with m ud the bare light-quark mass, has a finite continuum limit, 3 It is possible to argue that the expected dependence of λc on µ is the same also in the "first-order" picture.
The "first-order" correction changes the sharp gap λg into an effective gap λ ′ g < λg, corresponding to the mobility edge in the "first-order" picture. As we discussed above, λ ′ g should be close to M (φ) , or, at nonzero µ, to O(µ) = M (φ + µNT ) . The probability distribution function fµ(φ) of the Polyakov-line phases is a periodic function, and f−µ(−φ) = fµ(φ) due to charge-conjugation invariance. This implies O(µ) = O(−µ). In the ordered phase, even at nonzero µ, fµ(φ) is peaked around zero, approximately symmetric, i.e., fµ(−φ) = fµ(φ), and most likely monotonically decreasing in the interval [0, π]. Exploiting these properties, one can show that O(µ) is a decreasing function of |µ|. and the number of localised modes per unit volume, n loc , also remains finite in the continuum. Here ρ(λ) is the spectral density in lattice units. Within our picture, the fact that localisation survives the continuum limit is the nontrivial outcome of two competing effects. In fact, while the ordering of the Polyakov lines tends to disappear as a is reduced, thus lowering λ c , at the same time ρ(λ) at the low end of the spectrum is increased, due to the larger amount of fluctuations, and so of "islands" that can support localised modes. In Ref. [4] it has also been observed that n loc increases rapidly with T . At fixed lattice spacing, increasing the temperature brings λ c up, while at the same time decreasing the spectral density. Both phenomena are in agreement with our explanation, since increasing the ordering of the Polyakov lines, besides pushing up the mobility edge, also reduces the density of fluctuations and thus the spectral density of localised modes. 4 The value of n loc at fixed a results from the balance of these two effects. On top of this, to determine the physical value of n loc one has to take the continuum limit, which has an effect on λ c and ρ(λ) opposite to that of increasing the temperature, as we have said above. The actual behaviour of n loc as a function of the temperature ultimately depends on the relative magnitude of all these effects, which cannot be determined in our simple picture.

Ising-Anderson model
The considerations of the previous Section suggest that it should be possible to understand the qualitative features of the Dirac spectrum and eigenfunctions in QCD, for what concerns the localisation properties, by using a genuinely three-dimensional model. To construct such a model one has to strip off all the features that are irrelevant to localisation. The first step is to get rid of the time direction, reducing the lattice to three dimensions, and replacing the time covariant derivative in the Dirac operator with a diagonal noise term, intended to mimic the effective boundary conditions. This is motivated by the strong correlation among time-slices, as already mentioned in the previous Section. Furthermore, it is known that, in general, off-diagonal disorder is less effective than diagonal disorder in producing localisation [21,22], so it should be safe to replace spatial covariant derivatives with ordinary derivatives. More precisely, in our case the structure of the hopping terms is of the form considered in Ref. [27]. It is shown there that the width of the disorder distribution has to be rather large to produce localised modes near λ = 0 through off-diagonal disorder only. This is definitely not the case in QCD, where the disorder in the hopping terms involves unitary matrices. As we have already remarked, fluctuations of the spatial links have most likely a delocalising effect on the low modes. While this is certainly important for the detailed features of the spectrum and of the eigenmodes observed in QCD, in our picture it can nevertheless be regarded as a "second-order" effect, acting on the localised modes which should be produced by the "first-order" effect, i.e., the fluctuations of the Polyakov lines. To be precise, it is thus the "first-order" picture, discussed in the previous Section, that we are going to test here. Notice that replacing the spatial covariant derivatives with ordinary derivatives decouples the colour components of the quark wavefunction: this changes the symmetry class (in the sense of random matrix models), as we discuss in detail below.
If our "sea/islands" mechanism is viable, then the main features of localisation should still be captured by a genuine three-dimensional model of the following general form, if the diagonal noise term N x has the same features as the "diagonal noise", i.e., the effective boundary conditions, appearing in the true "Hamiltonian" −i / D. The relevant features are the following: • the diagonal noise involves continuous variables, namely the phases of the local Polyakov lines; • the diagonal noise terms are correlated, and are governed by the dynamics of the local Polyakov lines; • as the system is made more ordered, the diagonal noise tends to introduce a gap in the spectrum.
The simplest way to incorporate these features in N x is to base the diagonal noise on a spin model with continuous spins: this obviously satisfies the first requirement, and also the second one as Polyakov lines in the high-temperature phase display indeed a spin-model type dynamics [28,29]. Finally, also the third requirement is easily implemented if we take where s x is the spin variable at point x, and Λ is a constant determining the strength of the coupling of the fermions to the spins. The distribution of the spins s x is determined by the dynamics of an Ising-like model with continuous spins, which we take of the simplest possible form with nearest-neighbour interactions only, Since 1+s x 2 is 1 for "aligned" spins (i.e., for s x = 1), this choice provides indeed an effective spectral gap when the spins are ordered. 5 In the ordered phase of the Ising model there is a "sea" of s x ≃ 1 spins with "islands" of s x = 1 spins, so the underlying configurations have the same features as the Polyakov line configurations in QCD. In QCD we have a single parameter governing the ordering of the configuration and the size of the effective gap. In the effective model the ordering of the spin configuration is governed by β Ising , while the size of the gap is mainly determined by the spin-fermion coupling Λ, although one expects it to be affected also by the magnetisation of the system, and thus by β Ising .
We notice in passing a curious feature of this model: it belongs to different symmetry classes for lattices of even or odd size L, namely to the orthogonal class for L even, and to the symplectic class for L odd (see Appendix A). For our purposes it is convenient to work with even-sized lattices, which allow to get rid of the Dirac degree of freedom through a spin diagonalisation.
Let us now describe the spectrum of our effective model in more detail. The spectrum is symmetric with respect to the origin on a configuration-by-configuration basis, due to the fact that {γ 5 , H eff } = 0. In order to find the eigenvalues, it is convenient to first perform a spin diagonalisation through the following unitary transformation, where 1 is here the identity in Dirac space, and The Hamiltonian H obviously commutes with γ 4 , so one can diagonalise them simultaneously. Denoting with ψ ±n the common eigenvectors, Hψ ±n = λ ±n ψ ±n , γ 4 ψ ±n = ±ψ ±n , one has for the bispinors ψ ±n , (3.7) and the identical single-component eigenvalue equations where H (±) It is immediate to notice that there is an exact twofold degeneracy of the eigenvalues. It is also straightforward to prove that H (−) η 4 = −η 4 H (+) , so setting ξ i.e., ξ It is less straightforward to show that for a given noise configuration N x and a given eigenvector ξ + with eigenvalue λ + , there is another configuration with noiseN x and with the same Boltzmann weight, and an eigenvectorξ + corresponding to the eigenvalue −λ + . Indeed, takingN x andξ + as follows, one can explicitly verify that this is the case. This implies that the average spectrum of H (+) (and so that of H (−) ), obtained by integrating over the disorder (i.e., over the spin configurations), has an exact symmetry with respect to the origin.

Numerical results
We have performed numerical simulations of the effective Ising-Anderson model discussed in the previous Section, Eqs. In Fig. 1 we show the spectral density per unit volume of the system. Low modes have small spectral density, which rapidly increases as one goes up in the spectrum, as expected. Decreasing the temperature, thus making the system more ordered, decreases the spectral density of low modes, as one expects if these modes are localised on fluctuations of the spin variables. Increasing the spin-fermion coupling enlargens the region where the spectral density is small, again as expected, since it should push the effective gap up in the spectrum. The same effect is obtained by increasing β Ising , i.e., by making the system more ordered. Notice the symmetry under λ → −λ of the average spectrum, which has been discussed in the previous Section. There is apparently also a sharp gap in the spectrum, which is however not very sensitive to Λ, β Ising , or to the size of the system. This gap is probably a feature of the effective model, due to the drastic simplifications of the spatial hopping terms; in any case, it is irrelevant for our purposes. In Figs. 2 and 3 we show the average participation ratio PR = ( x |ψ(x)| 4 ) −1 /V , which gives a measure of the fraction of space occupied by a given mode, as a function of its location in the spectrum. 6 The PR changes by two orders of magnitude as one moves up in the spectrum starting from the low-density region; more importantly, when increasing the size of the system it remains almost constant in the bulk of the spectrum, while it visibly decreases near the origin. This signals that modes in the bulk are delocalised, while modes near the origin are localised. Notice in Fig. 3 the upturn of the PR at the low end of the spectrum: a similar phenomenon has been observed in QCD with staggered fermions [4].
Finally, in Fig. 4 we show a suitably defined local spectral statistics I λ across the spectrum. Spectral statistics can be used to detect a localisation/delocalisation transition, since the eigenvalues corresponding to localised or delocalised eigenmodes are expected to obey different statistics, namely Poisson or Wigner-Dyson statistics, respectively. More precisely, I λ is defined as where the unfolded level spacing distribution (ULSD) p λ (s) is the probability distribution, computed locally in the spectrum, of the so-called unfolded level spacing s i , i.e., the level spacing divided by the local average level spacing λ i+1 − λ i λ . The quantity I λ is simply the integrated ULSD up to the crossing points of the exponential distribution, corresponding to Poisson statistics, and the orthogonal Wigner surmise, which accurately describes the ULSD for (orthogonal) Wigner-Dyson statistics. The results confirm that the eigenmodes change from localised to delocalised when one moves up in the spectrum. Moreover, the mobility edge separating localised and delocalised modes goes up in the spectrum when the Ising system is made more ordered by decreasing the temperature, or when the spin-fermion coupling is increased. Our results give also a first indication that the slope of the curve increases as the volume is increased, thus hinting at the existence of a true phase transition in the spectrum. Notice that the transition takes place near the point Λ 1+ sx 2 in the spectrum, thus indicating that the position of the mobility edge is connected to the magnetisation of the system, as expected (see Table 1). This is also in agreement with our expectation for the position of the mobility edge in the "first-order" picture, discussed previously in Section 2.
To further investigate the onset of critical behaviour, we have determined eigenvalues and eigenvectors on larger volumes (L = 24, 32, 40) in the vicinity of the point in the spectrum where the localisation properties show a dramatic change, for β Ising = 1.0 and Λ = 1.5. We used 20k configurations for L = 24, scaling down the statistics for larger volumes by keeping approximately constant the total number of eigenvalues computed in the relevant spectral window. In Fig. 5 we show I λ as obtained from these simulations: one can easily appreciate the tendency of the cross-over from Poisson to Wigner-Dyson statistics to become steeper as the volume is increased.
Although the data are not of sufficiently good quality to allow a proper finite-size-scaling analysis, it is still possible to make a rough qualitative attempt at showing compatibility of the critical properties with those of the three-dimensional orthogonal Anderson model. According to the one-parameter scaling hypothesis, data from different volumes should collapse on a single curve if plotted against the scaling variable (λ − λ c )L  and ν the critical exponent of the correlation length. The pretty large distance between the crossing point of the L = 24, 32 data and of the L = 32, 40 data indicates the presence of sizable finite-size corrections to the position of the critical point, so we limit ourselves to the two largest volumes. In Fig. 5 we show the data collapse for the L = 32, 40 data. We determined λ c as the crossing point of second-order polynomial fits to the data, and used the value ν = 1.57 corresponding to the three-dimensional orthogonal Anderson model [30]. Although far from conclusive, this plot shows at least that agreement between the critical properties of the two models is possible.

Conclusions
In this paper we have discussed a possible mechanism leading to the localisation of lowlying modes of the Dirac operator in the high-temperature phase of QCD. Our proposal is that the localising "traps" are provided by spatial fluctuations ("islands") of the local Polyakov lines away from the ordered value (i.e., the identity in colour space). In this paper we have generalised the original proposal of Ref. [23] to SU (3), which has allowed us to identify the relevant variables for localisation, namely the phases of the eigenvalues of the local Polyakov lines, which govern the "depth" of the localising "traps". These phases enter the Dirac equation in the form of effective boundary conditions for the quark eigenmodes, and provide a three-dimensional source of disorder. The three-dimensional nature of the "islands", which constitute disconnected regions in a "sea" of (almost) ordered Polyakov lines, sheds light on the compatibility of the critical properties of the Dirac spectrum at the localisation/delocalisation transition with those of the three-dimensional unitary Anderson model [15].
To further substantiate our proposal, and check the viability of the "sea/islands" explanation, we have reproduced the qualitative features of the localised modes in a three-dimensional effective model, designed to produce localisation precisely through the proposed mechanism. This "Ising-Anderson" model is built replacing the four-dimensional (anisotropic) lattice, appropriate for finite-temperature QCD, with a three-dimensional one, and replacing the temporal covariant derivative in the Dirac operator with a diagonal (on-site) noise. The dynamics of the diagonal noise is governed by an Ising model with continuous spins in the ordered phase, in order to reproduce the qualitative features of the Polyakov line configurations. Localised modes are indeed present, and respond to changes in the parameters (i.e., the temperature of the Ising system and the strength of the spin-fermion coupling) as expected theoretically, thus reinforcing our confidence in the validity of the "sea/islands" explanation in QCD.
There are several possible extensions of the present work. Besides a thorough check of the critical properties of the model at the localisation/delocalisation transition, to be compared with those of the Anderson model in the appropriate symmetry class, the most interesting issue to study is how to properly model the effect of the spatial hopping. As a matter of fact, the diagonal disorder which we believe is responsible for localisation in the high-temperature phase, i.e., the effective boundary conditions, is present also in the low-temperature phase of QCD, but in that case it is ineffective in producing localisation. Indeed, in the lowtemperature phase the system looks "fully four-dimensional", and is basically insensitive to the temporal boundary conditions. In our opinion, in order to understand why this happens, and how the transition between the two phases is realised, adopting the point of view of "QCD as a random-matrix model" held in the present paper, it is necessary to understand how the spatial hopping acts in the two phases, especially concerning its capability of producing delocalised states around the origin. This could lead to useful insights in the subject of chiral symmetry breaking in QCD.

A Symmetry class of the effective model
The Hamiltonian of our toy model is of the form with A x y = a x δ x y , a x ∈ R, and (∂ j ) x y = 1 2 (δ x+, y − δ x−, y ). In this Appendix Latin indices run from 1 to 3, Greek indices from 1 to 4. We want to find an antiunitary "time-reversal" symmetry of this system. We denote with T the corresponding operator, which we parameterise as T = KΓU , with K the usual complex conjugation and U unitary to be determined, and Γ = γ 5 γ 4 γ 2 such that Γ † γ * µ Γ = γ µ , with Γ * = Γ, Γ † Γ = 1 and Γ 2 = −1. We have to impose T † HT = U † ΓH * ΓU = U † (γ 4 A − i γ · ∂)U = γ 4 A + i γ · ∂ . This imposes the following equations, which are solved by so for odd lattices the symmetry class is the symplectic one. For even lattices, it is possible to choose U x = e iφ (−1) k x k −x j iγ 5 γ j for any j, which yields T 2 = −(γ 5 γ j ) 2 = 1 , (A.13) so that the symmetry class is the orthogonal one (the existence of at least one antiunitary symmetry operator with square one is sufficient to see this). Even more simply, after spin diagonalisation the Hamiltonian reads H x y = η 4 ( x)A x y + i η( x) · ∂ x y , (A.14) and one can easily identify the time-reversal operator T x y = η 4 ( x)δ x y K, which has T 2 = 1.