An Anderson-like model of the QCD chiral transition

We study the problems of chiral symmetry breaking and eigenmode localisation in finite-temperature QCD by looking at the lattice Dirac operator as a random Hamiltonian. We recast the staggered Dirac operator into an unconventional three-dimensional Anderson Hamiltonian ("Dirac-Anderson Hamiltonian") carrying internal degrees of freedom, with disorder provided by the fluctuations of the gauge links. In this framework, we identify the features relevant to chiral symmetry restoration and localisation of the low-lying Dirac eigenmodes in the ordering of the local Polyakov lines, and in the related correlation between spatial links across time slices, thus tying the two phenomena to the deconfinement transition. We then build a toy model based on QCD and on the Dirac-Anderson approach, replacing the Polyakov lines with spin variables and simplifying the dynamics of the spatial gauge links, but preserving the above-mentioned relevant dynamical features. Our toy model successfully reproduces the main features of the QCD spectrum and of the Dirac eigenmodes concerning chiral symmetry breaking and localisation, both in the ordered (deconfined) and disordered (confined) phases. Moreover, it allows us to study separately the roles played in the two phenomena by the diagonal and the off-diagonal terms of the Dirac-Anderson Hamiltonian. Our results support our expectation that chiral symmetry restoration and localisation of the low modes are closely related, and that both are triggered by the deconfinement transition.


Introduction
The low end of the spectrum of the Euclidean Dirac operator plays an important role in determining the properties of hadronic matter in Quantum Chromodynamics (QCD). In particular, the spectral density around the origin is closely tied to the fate of chiral symmetry, and entirely determines it in the chiral limit [1]. It is thus not surprising that the low end of the Dirac spectrum behaves differently in the broken and in the restored phase. The most important difference is obviously that while in the broken phase eigenvalues accumulate around the origin, in the restored phase the spectral density vanishes there. Besides this, or perhaps as a consequence, the low-lying eigenmodes display different localisation properties and statistical behaviour.
Numerical simulations of QCD on the lattice have shown that the low-lying Dirac eigenmodes, while delocalised on the entire lattice volume at low temperature [2,3], become spatially localised at high temperature [4][5][6][7][8][9][10][11][12][13], above the chiral crossover [14,15]. Although most of the studies of this phenomenon have been carried out using the staggered discretisation of the Dirac operator [4,5,[7][8][9][10][11], there is also evidence in simulations with overlap [6,8] and domain wall [12,13] fermions. Let us summarise the current knowledge about it (see Ref. [16] for a review), focussing on the case of the staggered operator. In this case the eigenvalues iλ are purely imaginary and the spectrum is symmetric with respect to zero, so it suffices to A possible explanation of the fact that the two models share the same universality class has been proposed in Ref. [29], elaborating on a previous proposal made in Ref. [8] for the mechanism leading to localisation. The main idea is that high-temperature QCD effectively contains a source of 3D on-site disorder, which consists of the phases of the local Polyakov line at a given spatial point. Above T c the Polyakov lines get ordered along the identity, with local fluctuations away from it. The Polyakov line phases affect the quark eigenmodes through effective boundary conditions in the temporal direction, and make favourable for the quark wave function to "live" on the "islands" of unordered Polyakov lines. As a random matrix model, the Dirac operator in high-temperature QCD is thus effectively 3D with diagonal noise. Support to the viability of this mechanism was obtained in Ref. [8] by studying the correlation of the Dirac eigenfunctions with the fluctuations of the Polyakov loop on SU (2) gauge configurations. In Ref. [29] we looked for a different kind of evidence: we constructed, and studied numerically, a QCD-inspired toy model which should display localisation precisely through the proposed mechanism, but in a much simplified setting. This "Ising-Anderson" model is essentially obtained by removing the temporal direction, thus working in threedimensional space (i.e., a single time slice), and by mimicking the effect of the Polyakov lines on the quark wave functions with a (continuous) spin model of the Ising class in the ordered phase, used to generate the appropriate diagonal noise. This model adequately reproduces the main features of localisation, in particular their qualitative dependence on the amount of "islands" of "wrong" spins in the "sea" of ordered spins, as expected from the proposed explanation.
Although the Ising-Anderson model of Ref. [29] provides a satisfactory qualitative description of localisation in the high-temperature phase of QCD, it fails completely at describing the low-temperature phase. Indeed, simulations in the disordered phase of the underlying Ising model (not reported in Ref. [29], but see Section 4 below for similar results) fail to reproduce the most distinctive feature of the Dirac spectrum at low temperature, namely the presence of a nonzero spectral density around the origin, which leads to the spontaneous breaking of chiral symmetry [1]. Instead, a sharp gap separates the lowest eigenvalues from the origin, both when the underlying spin model is in the ordered phase, and when it is in the disordered phase. While in the former case this is in rough qualitative agreement with the small density of near-zero modes in QCD at high temperature, in the latter case this is at odds with the existing numerical results. On the other hand, the mechanism through which the Polyakov lines affect the quark wave functions, i.e., the effective boundary conditions, is in principle at work at all temperatures, and only the presence of order in the (relevant) Polyakov-line configurations, or lack thereof, distinguish the two phases. As we have already said above, the effectiveness of the "sea/islands" mechanism devised in Refs. [8,29] can explain the presence of localised modes at high temperature. For the explanation to be complete, the ineffectiveness of this same mechanism should also explain the absence of localised modes at low temperatures. This requirement is even more compelling in light of the apparently very close relation between localisation and the QCD crossover: the reason for the ineffectiveness of the mechanism at low temperature is likely to be also the reason for the finiteness of the spectral density at the origin.
Above we said that the effective boundary conditions are in principle at work at all temperatures but, as a matter of fact, they turn out to be irrelevant at low temperatures, where QCD looks "effectively 4D". The ineffectiveness of the boundary conditions at low temperatures clearly entails the ineffectiveness of the "sea/islands" mechanism for localisation, and so it may seem hopeless to try to give a simultaneous description of the QCD crossover and of the appearance of localised modes by means of a 3D model. Still, one can technically treat the quark degrees of freedom along the compactified temporal direction as internal degrees of freedom of a quark living in 3D space. In doing this, one can in principle make the temporal extent of the lattice as small as possible in lattice units, thus reducing the amount of such internal degrees of freedom to a minimum. In this setting, whether the system is effectively 3D or 4D becomes a question about the correlation among the internal degrees of freedom. For this question to make any sense at all, the number of such degrees of freedom has clearly to be at least two. Perhaps a more natural way of phrasing this is in terms of the temporal correlation length, and the minimal number of time slices for which one can ask if they are strongly correlated is obviously two.
With this (in hindsight rather obvious) insight, our purpose in this paper is to build a refined toy model, aimed at describing the (supposedly) simultaneous appearance of localised modes and recovery of chiral symmetry, as signalled by the vanishing of the spectral density at the origin. 1 The motivation is twofold. On the one hand, we want to implement the "sea/islands" mechanism in a model that reproduces QCD more faithfully (at the qualitative level), in order to make the case of Refs. [8,29] stronger. On the other hand, we want to investigate the connection between localisation and the deconfinement/chiral transition in a simple and controllable setting. This could also lead to some insight into the chiral transition, and into its relation to deconfinement, from the point of view of the QCD Dirac operator as a random matrix model, independently of the issue of localisation. Indeed, since localisation of the low modes and restoration of chiral symmetry take place together, and this happens near the deconfinement transition, it is very likely that they are both triggered by the ordering of the gauge configurations, so that some mechanism could be devised which would explain both phenomena.
Here is the plan of the paper. Before building the toy model, in Section 2 we cast the staggered Dirac operator into a three-dimensional Hamiltonian, with internal degrees of freedom corresponding to colour and to the lattice temporal momenta. In this way the connection with Anderson-type models is made transparent. In Section 3 we write down explicitly our toy model, which we study numerically in Section 4. Finally, in Section 5 we discuss our results, state our conclusions and comment about future prospects. Some technical details are discussed in Appendix A.

Staggered Dirac operator as a random matrix model
In this Section we recast the Dirac operator in lattice QCD, / D, as the Hamiltonian of an unconventional three-dimensional Anderson model. We work here with the staggered Dirac operator for simplicity.
The basic idea is to split the "Hamiltonian" H = −i / D into a "free" and an "interaction" part, H = H 0 + H I , and then work in the basis of the "unperturbed" eigenvectors of H 0 . At finite temperature the temporal direction is compactified and therefore singled out, and so the physically most sensible choice is to identify the free Hamiltonian with the temporal hoppings. This leaves the spatial hoppings as the (spatially isotropic) interaction part. We thus define and depend only on the spatial coordinates, 0 ≤ x i ≤ L − 1, i = 1, 2, 3. Both N T and L must be even. We take the gauge links U µ (t, x) ∈ SU (N c ) for generality, and denote the backward links by x −), j = 1, 2, 3. Periodic boundary conditions are understood for the gauge links, while on the quark wavefunctions antiperiodic boundary conditions in the temporal direction, and periodic boundary conditions in the spatial directions, are imposed. In Eq. (2.1) colour indices are suppressed.
The eigenvectors of H 0 are easily determined. To this end, let us define the following gauge transporters in the temporal direction, with P N T ( x) = P ( x) the usual (untraced) Polyakov line starting at t = 0 and winding around the temporal direction. Let furthermore ϕ a ( x) be the (ortho)normalised eigenvectors of P ( x), Here x runs over the whole 3-volume, a = 1, . . . , N c with N c the number of colours, and each ϕ a has N c colour components, (ϕ a ) i . The eigenvalues e iφa( x) have unit absolute value and satisfy a e iφa( x) = 1. The eigenvectors of the Polyakov line can be used to build the eigenvectors ψ x a k 0 of H 0 . These are localised on a single spatial point, x, have a well-defined temporal momentum, k, and carry a colour quantum number, a, and read (in the coordinate basis) with colour-and space-dependent effective Matsubara frequencies, ω ak ( x), The form of ω ak ( x) results from imposing temporal antiperiodic boundary conditions on the fermions, and from the presence of nontrivial Polyakov lines which modify the free-field result.
The ψ x a k 0 are eigenvectors of H 0 with "unperturbed" eigenvalues λ x a k 0 given by In the basis {ψ x a k 0 }, the operator H I has vanishing diagonal elements, and only nearest-neighbour hopping terms (i.e., connecting eigenvectors localised on nearest-neighbour sites), ±j (t, x) is just the usual link variable in the temporal diagonal gauge, i.e., the temporal gauge where the Polyakov lines have been diagonalised. 2 The full Hamiltonian in the basis of "unperturbed" eigenvectors, {ψ x a k 0 }, will be denoted by H, and it carries space, colour, and temporal-momentum indices, [H( x, y)] ak,bl . Suppressing the latter two indices, we will write where we have introduced the following notation, (2.12) Hermiticity implies V † ±j = V ∓j , as one can also verify explicitly. Moreover,Ũ ±j (t, x) are unitary matrices in colour space, and V ±j ( x) are unitary matrices in (joint) colour and temporalmomentum space. The phases φ a ( x) are defined only modulo 2π, so that in order to fully specify H one needs to pick a convention. The simplest and most sensible possibility is to take φ a ( x) ∈ [−π, π) for a = 1, . . . , N c − 1, and impose 2πq( x) ≡ a φ a ( x) = 0 for all x. This is what we do from now on, unless otherwise specified. In this way one avoids introducing spurious non-uniformities in the Hamiltonian that might obscure the important features. In principle, however, any other choice, possibly different on different lattice sites, is legitimate. One can show that the Hamiltonians obtained with different choices are related by unitary transformations, so that the spectrum does not depend on the convention used, as it should be. If one and the same convention is adopted for q( x) ∈ Z on all lattice sites, thenŨ ±j (t, x) and V ±j ( x) are also unimodular. More details on these issues are given in Appendix A.
The Hamiltonian Eq. (2.11) looks like that of a 3D Anderson model, but with an antisymmetric rather than symmetric hopping term, and moreover carrying internal degrees of freedom. We will sometimes refer to it as the Dirac-Anderson Hamiltonian. Due to the fluctuations of the gauge links from configuration to configuration, both diagonal and off-diagonal disorder are present. The amount of disorder is thus controlled by the size of such fluctuations, and thus ultimately by the temperature of the system (and by the lattice spacing).
The diagonal disorder originates entirely from the Polyakov lines. In contrast with the usual Anderson Hamiltonian, changing the amount of disorder does not change the strength of the diagonal term, since the "unperturbed" eigenvalues are always bounded by one in absolute value. On the other hand, on the two sides of the deconfinement transition the shape of their distribution is different, with an enhancement of "unperturbed" eigenvalues corresponding to the trivial phase at high temperature. Moreover, in the high-temperature phase there is long-range order in the diagonal term, as a consequence of the long-range order in the Polyakov-line configuration.
The off-diagonal disorder in the hopping terms is mostly determined by the spatial links. As with the diagonal disorder, while the overall "size" of the hopping term does not change with temperature, as it is a unitary matrix in any case, its typical matrix structure changes considerably across the transition. Indeed, the most interesting property of the hopping term is that whenŨ (2.13) For this to happen we need both ±j (0, x) ∀t. This is the case when the neighbouring Polyakov lines are both close to the identity, 3 which in turn causes a strong (local) correlation among spatial link variables. At low temperatures the Polyakov lines are disordered, and so this does not occur often: we thus expect that typically V ±j ( x) will have non-negligible off-diagonal terms in temporal-momentum space, which leads to strong mixing of the wave-function components corresponding to different temporal momenta. At high temperature, on the other hand, the Polyakov lines get ordered, and there are large spatial regions where this is approximately true: in these regions V ±j ( x) gets ordered along the identity in temporal-momentum space, and so in these regions the components of the wave function corresponding to different temporal momenta are coupled only weakly. 4 In other words, at high temperature we expect the "correlation length" in temporal-momentum space to become shorter than the size of the system (again, in temporal-momentum space): this is what typically happens when a transition to a disordered phase takes place. Keep in mind that here the system under consideration are the quark eigenfunctions, and not the gauge fields, and moreover that this shortening of the "correlation length" is a local effect (although taking place in the whole "sea" of ordered Polyakov lines).
Finally, we want to remark that the spectrum of H is obviously symmetric with respect to λ = 0. In the new basis, this is seen as a consequence of the anticommutation relation {Q, H} = 0 of the Hamiltonian with a certain unitary matrix Q. Moreover, in the case of gauge group SU (2) the Hamiltonian is invariant under an antiunitary transformation T with T 2 = −1, which implies a double (Kramers) degeneracy of the spectrum, and moreover that the Hamiltonian belongs to the symplectic class of random Hamiltonians [23]. These are nothing but the analogues of the properties of the usual staggered Dirac operator, only expressed in a different basis. Details on the form of Q and T in the new basis are summarised in Appendix A.
i.e., a given temporal-momentum component k mixes only with the "opposite" component k

Construction of the toy model
In this Section we want to build explicitly a simple model describing the change both in the localisation properties and in the spectral density of the low Dirac modes, starting from the Dirac-Anderson Hamiltonian, i.e., the reformulation of the Dirac operator described in the previous Section.

Brief review of gauge dynamics
We begin by reviewing the features of the dynamics of gauge theories that we expect to be relevant for the two above-mentioned phenomena. As we have already pointed out in the Introduction, the Polyakov line is expected to play an important role in the localisation of low modes at high temperature. As is well known, in pure gauge theory the center symmetry of the action is spontaneously broken in the deconfined phase at high temperature, where the Polyakov lines get ordered along one of the center elements. Fluctuations away from the ordered value form "islands" within a "sea" of ordered Polyakov lines. In the low temperature, confined phase the center symmetry is restored and Polyakov lines are disordered. Dynamical fermions in the fundamental representation break the center symmetry explicitly but softly, so that they do not change this picture, besides selecting the trivial vacuum in the hightemperature phase.
As we observed at the end of the previous Section, the ordering of the Polyakov lines in the high-temperature phase of QCD will induce strong correlations between spatial links at the same spatial point but on neighbouring time slices. This is a local effect, i.e., it depends on the ordering of the Polyakov lines at the spatial points connected by the given link. To see this, it is convenient to work in the temporal diagonal gauge [see Eq. (2.10)], where the contributions to the Wilson gauge action read (up to a constant factor) away from the temporal boundary, while at the temporal boundary we have Here SP stands for the contribution of the spatial plaquettes. Although only U interact directly with the Polyakov lines, this effect propagates to the spatial links on the other time slices, as they are coupled according to Eq. (3.1). Moreover, Eq. (3.2) shows that the dynamics of the Polyakov lines is affected by the backreaction of the spatial links.

From QCD to the toy model
Let us now turn to the explicit construction of the model. We want to study the spectrum of a random Hamiltonian of the form with d( x) and v ±j ( x) respectively diagonal and unitary in joint colour and temporal-momentum space. This is of the same form as the staggered Dirac operator in the basis of "unperturbed" eigenvectors, Eqs. (2.11) and (2.12), and both the diagonal and the hopping terms will be modelled on Eq. (2.12), replacing the Polyakov line phases φ a ( x) and the spatial links U (td) j (t, x) with similar quantities in the toy model. This implies in particular that the toy model Hamiltonian will have a symmetric spectrum, and will belong to the same class of random Hamiltonians as the one found in (N c -colour) QCD.
Since our purpose is to build a model simpler than QCD, yet displaying the same behaviour concerning the phenomena of localisation of eigenmodes and accumulation of eigenvalues near the origin, we will take QCD as a starting point and eliminate all those features that we deem irrelevant. We will first of all neglect the backreaction of the quark eigenmodes in the partition function, omitting the fermion determinant (i.e., making the quenched approximation), since it is known that these phenomena are present in pure gauge theory as well.
Next, since we aim at reproducing these phenomena only qualitatively, we will simplify the dynamics of the toy model analogues of φ a ( x) and U (td) j (t, x) with respect to QCD. As in Ref. [29], the main simplifying idea is to mimic the effect of the Polyakov line phases on the quark wave functions by spin-like variables. However, in this work we want to achieve a closer resemblance to the actual dynamics of the phases. To this end, we want to design the spin model so that the effective potential for the magnetisation, in the ordered phase, is similar to that for the Polyakov line phases in QCD [30,31], or more generally in an SU (N c ) gauge theory. The potential should therefore develop N c minima in the ordered phase, corresponding to the N c Polyakov-loop vacua. A possibility is to choose N c complex spin variables s x , corresponding to the N c eigenvalues of an SU (N c ) Polyakov line, that satisfy and which obey the dynamics determined by the following Hamiltonian, where p(x) = diag (s 1 x , . . . , s Nc x ). The N c -dependent factors are chosen for convenience. The first term corresponds to a lattice sigma-model possessing a global [U (1)] Nc−1 symmetry. The second term mimics the absolute value squared of the trace of the Polyakov line (i.e., the Polyakov loop), and at h = 0 breaks the symmetry down to Z Nc . 5 This residual symmetry can hold dynamically or be spontaneously broken, with precisely N c vacua s a x = e i 2π Nc k ∀ x, a, with k = 0, . . . , N c − 1. 6 Parameterising the spins as follows, the Hamiltonian can be recast as up to an irrelevant additive constant. The dynamics of the phases φ a x resembles qualitatively that of the Polyakov line phases: while at low β the (complex) magnetisations m a = L −3 x s a x vanish on average, for large enough β the system transitions to an ordered phase, with φ a x aligning to one of the vacuum values discussed above. Small and large β in the spin model thus correspond to small and large temperatures in QCD. It is then natural to take the diagonal term d( x) to simply be D( x) with the Polyakov line phases replaced by Notice that φ a x will obey their own independent dynamics, unaffected by the analogues of U x and with appropriate SU (N c ) matrices u j (t, x), respectively: . The last step is to define the dynamics of u j (t, x). The important feature we want to mimic from QCD are the local correlations between spatial links across time slices induced by the Polyakov lines. On the other hand, we expect the correlations induced by the spatial plaquettes in Eqs. (3.1) and (3.2) to be less important. We will thus drop them from the action, keeping only the contributions of the temporal plaquettes. The 5 Further unbroken symmetries are the one under "charge conjugation", s a x → s a x * , and the one under permutations of the s a x with respect to a. 6 The minimum of Hnoise is achieved for spatially uniform phases (modulo 2π), s a x = e iφ a , satisfying φ a = φ b mod 2π ∀a, b. For a general parameterisation of the phases the constraint reads a φ a = 0 mod 2π, which leads to φ a = 2π Nc k mod 2π ∀a, k = 0, . . . , Nc − 1.
Boltzmann weight for the configurations of the spatial links thus factorises, with each factor involving only the spatial links at a given spatial point and along a given spatial direction. Explicitly, the dynamics of the toy model links u j (t, x) will be governed by the following action, whereβ is a constant playing the role of the gauge coupling. Expectation values are defined as follows: , there is no backreaction of the link variables on the spins, and φ a x acts as a background field for u j (t, x).

Minimal toy model
Let us describe the model in detail in the simplest case, namely taking N T = 2 and N c = 2, i.e., the minimal possible values. This is the model we have employed in the numerical study discussed in the next Section. For N c = 2, the basic variables are the complex spins s x ≡ s 1 x , and the SU (2) link variables u j (t, x), t = 0, 1. The noise Hamiltonian governing the spin dynamics, Eq. (3.7), simplifies to The U (1) symmetry of the XY model, represented by the first term, is broken to Z 2 by a nonzero "external field", h, appearing in the second term. We thus expect our spin model to belong to the universality class of the three-dimensional Ising model. Concerning the dynamics of the link variables, for N T = 2 there is a single term δS j (0, x) = δS j (1, x) ≡ δS j ( x) determining the Boltzmann weight for j-links at x, namely where p( x) = diag (e iφ x , e −iφ x ), and the action for the links reads simply S u = − x,j δS j ( x). Averages are defined according to Eq. (3.12) above. It is worth reminding the reader that while in QCD there is a single coupling that enters the dynamics of both the Polyakov lines and the spatial links, in the toy model β andβ can be varied independently.
The toy model Hamiltonian H toy ( x, y), Eq. (3.3), mimicking the QCD Dirac operator, consists of a diagonal and a hopping term, both containing disorder. For the on-site, diagonal noise terms d( x), since we have simply As for the hopping terms v j ( x), they read explicitly

Numerical results
In this Section we discuss our numerical results for the toy model defined in the previous Section, both in the ordered and in the disordered phases of the underlying spin model. For simplicity, we have studied the "minimal" case N c = N T = 2. The "gauge coupling",β, was fixed toβ = 5.0. Since we are interested mostly in the dependence on β, we have fixed the coefficient h of the symmetry-breaking term to h = 1.0. We have first studied the spin model on its own to determine the corresponding phase structure. In Fig. 1 we show the magnetic susceptibility χ of the spin model, as a function of β. A phase transition is expected to occur in the thermodynamic limit at a critical β c , with β c ≈ 0.3. For β 0.29 one is safely in the disordered phase, while for β 0.31 one is in the ordered phase, and finite-size effects should not affect the qualitative behaviour of our random toy Hamiltonian H toy . In Fig. 2 we show the distribution of phases in a single typical configuration below and above the transition. The tendency of the system to get ordered is evident. Since our toy model is quenched, in the ordered phase of the spin model we have to select the appropriate vacuum by hand. The appropriate vacuum is of course that in which the phase of the magnetisation m = N −1 c L −3 x,a s a x = |m|e iϕm is zero, corresponding to the trivial Polyakov loop sector selected by fermions in QCD. This is done in practice only at the level of the random Hamiltonian, by "rotating" the spins, i.e., aligning their phase to zero, when the spin model is in a different vacuum. In the case N c = 2 considered in this paper this is easy to implement. For typical configurations in the ordered phase, the (complex) magnetisation is usually close to being real, i.e., e iϕm ≃ ±1. When cos ϕ m < 0, we "rotate" all the spins by replacing s x → −s x . In terms of the phases φ x ∈ [−π, π), this is implemented through The spectral density ρ(λ) of the random Hamiltonian in the two phases is shown in Fig. 3: while in the disordered phase (β < β c ) there is an accumulation of eigenvalues near the origin, so that (presumably) L −3 ρ(0) = 0 in the thermodynamic limit, in the ordered phase (β > β c ) this region is depleted, and ρ(0) = 0. The dependence on β is rather mild in the disordered phase, while the spectral density gets rapidly suppressed as β grows in the ordered phase. The expected connection between magnetisation in the spin system and "chiral" transition in the spectrum of the random Hamiltonian indeed shows up. In order to understand the nature of the lowest eigenmodes in the two phases, it is convenient to study the statistical properties of the corresponding eigenvalues. In fact, localised modes are expected to fluctuate independently, so that the corresponding eigenvalues should obey Poisson statistics. Delocalised modes, on the other hand, mix strongly under fluctuations and are expected to obey the appropriate Wigner-Dyson statistics, which in the case at hand is the one corresponding to the Gaussian Symplectic Ensemble (GSE). A convenient observable to distinguish the two cases is the spectral statistics I 0.5 [32,33], where P ULSD (s) is the probability distribution of the unfolded level spacings s j = λ j+1 −λ j λ j+1 −λ j , where λ j+1 − λ j is the average level spacing in the spectral region corresponding to the given level. Unfolding is a mapping of the eigenvalues that makes the spectral density equal to 1 throughout the spectrum. In practice, we have ordered the eigenvalues obtained in all the explored spin configurations and replaced them by their rank divided by the number of configurations. 7 The P ULSD (s) is known both in the case of Poisson statistics, where it is the exponential distribution P P (s) = exp(−s), and in the case of the GSE, where it is very precisely approximated by the symplectic Wigner surmise P GSE [23], The quantity I 0.5 is sensitive to the behaviour of P ULSD (s) near s = 0, and so it is very different for Poisson (where P P ∼ 1 near s = 0) and GSE (P GSE ∼ s 4 ) statistics. Indeed, I GSE 0.5 ≃ 0.0487 for the GSE, while for the Poisson ensemble I P 0.5 ≃ 0.393. The choice of the upper limit of integration in Eq. (4.3) is made in order to maximise the difference between these two values, as P P and P GSE cross near s = 0.5. In Fig. 4 we show the behaviour of I 0.5 as one moves along the spectrum, i.e., computing I 0.5 locally, using only eigenvalues in disjoint bins of fixed width, and assigning the result to the average of the eigenvalues in that bin. While in the disordered phase one finds Wigner-Dyson statistics throughout the whole spectrum, in the ordered phase the lowest modes have near-Poisson statistics, and become more independent as the volume is increased. Above a β-dependent point λ c in the spectrum, the modes have near-Wigner-Dyson behaviour, and more and more so when the volume is increased. This hints at a localisation/delocalisation transition in the spectrum taking place in the thermodynamic limit. As the system is made more ordered, λ c increases, which is in agreement with our expectations: qualitatively, λ c should behave like the spatial average of the lowest effective Matsubara frequency [29], which indeed grows as the ordering of the system is increased. Our toy model is therefore able to reproduce localisation of the low modes in the ordered phase, and the qualitative dependence of the "mobility edge" λ c on the ordering of the system.
The results discussed above show that our toy model successfully reproduces the important features of QCD for what concerns localisation and "chiral symmetry" breaking/restoration, i.e., the accumulation or not of eigenvalues near the origin. This indicates that we have indeed kept all the important aspects of the Dirac operator and of the gauge dynamics in our model, which can thus be used to gain some reliable qualitative insight into the properties of the Dirac eigenvalues and eigenmodes. Therefore, although the following discussion deals explicitly with the toy model, it can be directly translated to the physically relevant case of QCD by replacing "spins" with Polyakov line phases.

Variations on the toy model
We want now to test how the features of the toy model that we borrowed from QCD affect the "chiral" transition, i.e., the drop in the spectral density near the origin, and localisation of the lowest modes. To this end, we study some ad hoc modifications of the toy model.
The "chiral" transition and localisation of the lowest modes are clearly connected in our toy model, as they are evidently tied to the magnetisation of the spin system. The magnetisation affects our toy model Hamiltonian in two ways: it creates a "sea" of ordered spins, where "islands" of fluctuations provide an "energetically" convenient place for an eigenmode to localise [8,29]; and it locally correlates the gauge links on different time slices, thus leading to the approximate local decoupling of different temporal-momentum components of the eigenfunctions. It is an interesting question how important each of these two effects is for "chiral symmetry" breaking and localisation.
As a matter of fact, the appearance of "islands" alone is not enough to produce either "chiral symmetry" restoration or localisation. In Fig. 5 we show the spectral density and the spectral statistic I 0.5 obtained in the ordered phase at β = 0.32 when settingβ = 0, i.e., imposing no correlation between the time slices. In this case, despite the presence of "islands" in the spin configurations, the spectral density at the origin remains finite, and the low-lying eigenmodes do not localise. This means that to restore "chiral symmetry" one needs that the mixing of different temporal-momentum components of the quark wave functions be suppressed to some extent. Only then can the "islands" effectively act as localising centers for the low modes, since, loosely speaking, the lowered spectral density would make it difficult to mix for modes localised in different regions.
On the other hand, in the disordered phase the absence of "islands" is not sufficient to ensure "chiral symmetry" breaking and prevent localisation of the lowest modes. An essential ingredient for the accumulation of eigenvalues around the origin is the fact that the hopping term has typically sufficiently large off-diagonal components in temporal-momentum space. In Fig. 6 we show the spectral density obtained in the disordered phase at β = 0.29 when setting the off-diagonal part of the hopping term to zero. 8 The Hamiltonian in this case is block-diagonal in temporal-momentum space, with blocks H [±] of the form where [c( x)] ab = cos φ x 2 δ ab , and v (+) ±j have been defined in Eq. (3.17). Apart from the presence of colour, and of disorder in the hopping terms, this is precisely the Hamiltonian considered in Ref. [29]. In this case the spectrum displays a sharp gap near the origin. In Fig. 6 we also show how the spectral statistic I 0.5 changes along the spectrum. Our results show that the lowest-lying modes are localised, while higher up in the spectrum they delocalise, as signalled by the decrease of I 0.5 . The reason why I 0.5 does not tend to the GSE value in the bulk, but remains clearly above it, can be ascribed to the fact that the spectra of H [±] are separately symmetric on average with respect to the origin. The proof of this property can be easily obtained extending the one reported in Ref. [29] to the case of nontrivial hopping terms and in the presence of colour. Combining the two spectra together to obtain that of H toy , one expects that the latter will have an approximate double degeneracy (on top of the Kramers degeneracy) of the eigenvalues. This naturally leads to an enhancement of the ULSD near the origin (levels like to lie closer than on average) with respect to the Wigner-Dyson distribution, and therefore to an enhancement of I 0.5 .
The results above show that the mixing of different temporal momentum components is a necessary condition to have chiral symmetry breaking in the disordered phase, and that its (local) attenuation is a necessary condition to have restoration in the ordered phase. This entails that the "minimal" model studied here is indeed minimal, if one wants to reproduce qualitatively the features of the Dirac spectrum and of the corresponding eigenmodes. It is thus not possible to neglect entirely the temporal direction, as in Ref. [29], if one wants to correctly describe the disordered phase. Furthermore, in the present setting one cannot neglect the correlation of spatial links across time slices if one wants restoration of chiral symmetry and localisation of the low modes in the ordered phase. 9 However, one cannot certainly conclude that those described above are also sufficient conditions, especially because of the extreme nature of the two modifications of the toy model that we have studied. In a realistic situation, correlations across time slices are always present to some extent, but never perfect, so that mixing of the temporal components of the quark wave functions is never completely free or strictly forbidden. One thus expects some other effect to compete against mixing, with the fate of "chiral symmetry" and localisation being decided by the strongest of the two effects. The obvious candidate is the nature of the diagonal disorder, i.e., of the distribution of the unperturbed eigenvalues: for a low density of small unperturbed modes it should be more difficult to accumulate eigenvalues around the origin, and conversely for a high density of such unperturbed modes a weak mixing of temporal momentum components could be sufficient.
The latter case is naturally illustrated by changing the boundary conditions for the fermions in the temporal direction from antiperiodic to periodic. The derivation of the Dirac-Anderson Hamiltonian proceeds exactly in the same way, and leads to the same results, Eqs. (2.11) and (2.12), up to replacing the "fermionic" Matsubara frequencies ω ak ( x) with the "bosonic" frequencies ω PBC ak ( x), The construction of the toy model is also unchanged, up to a similar replacement in the diagonal terms, Eq. (3.8). In the minimal setting N T = N c = 2, this amounts to replace cos φ x 2 → sin φ x 2 . In the ordered phase the distribution of unperturbed eigenvalues is now peaked around zero, while the correlations between spatial links and thus the strength of the mixing of the temporal momentum components are unchanged. Numerical results for this model are shown in Fig. 7: "chiral symmetry" is broken and localisation is absent. Notice that this is precisely what one expects, if the toy model is to correctly reproduce the features of QCD [34,35], which further supports the viability of our model to study the qualitative behaviour of the Dirac spectrum.
To illustrate the former case we have employed a more artificial construction, increasing "by hand" the ordering of the unperturbed eigenvalues while keeping unchanged the hopping terms. This has been achieved by mapping the absolute value of the unperturbed eigenvalues, x with 0 < τ < 1. As we show in Fig. 8, in this way one can deplete the spectral region around zero: the strength of the temporal-momentum-component mixing is not sufficient to compensate for the lowered density of small unperturbed modes.

Discussion
Summarising our findings, the fate of "chiral symmetry" and of localisation depends on the competition between ordering of the unperturbed modes and mixing of temporal momentum components of the wave functions. In the disordered phase, when there is no ordering of the unperturbed modes and a sizeable amount of small such modes is present, there is actually no competition and one expects chiral symmetry to be broken. In the ordered phase, the reduction of mixing can be compensated by changing the boundary conditions, so increasing the amount of small unperturbed modes. Although we have checked here only the case of periodic boundary conditions, one can in principle change the effective boundary conditions in a continuous manner by introducing an imaginary chemical potential, µ: periodic boundary conditions thus correspond to |µN T | = π, and for |µN T | sufficiently close to π we still expect to find a finite spectral density at the origin. Sticking to N c = 2 and feeding this back into the partition function by including the fermionic determinant, in the ordered/deconfined phase one is lead to expect that if at µ = 0 the trivial center sector is selected, at |µN T | = π the system is again in the ordered/deconfined phase but in the other center sector. For sufficiently high temperature one thus expects to find a transition from one vacuum to the other when moving along the µ direction in the phase diagram, which then repeats due to periodicity in µ: these are nothing but the Roberge-Weiss transitions [36].
Our results shed some more light on the "sea/islands" mechanism in QCD, and on its ineffectiveness at low temperature. As we have said above, in the disordered phase the hopping term has typically sizeable off-diagonal entries, which can effectively mix different temporalmomentum components of the quark wave function. This is due to the presence of extended spatial regions where the local correlations across time slices are sufficiently weak, because of the lack of order in the underlying spin system. As we saw above, this is necessary for the accumulation of small eigenvalues, and one can think of mixing as somehow "pushing" the eigenvalues towards zero. On top of that, the lack of order also provides a sizeable density of small unperturbed modes, which also favours such an accumulation. The combination of these two effects thus leads naturally to a finite spectral density near the origin, and to delocalisation of the low modes. The ineffectiveness of the "sea/islands" mechanism at low temperatures is thus due to the fact that in that case there simply are no "islands". In contrast to this, in the ordered phase the regions where mixing is effective are localised on the "islands" where the spins fluctuate away from the ordered value. In the "sea" of ordered spins the unperturbed eigenvalues are large (i.e., close to 1) and the mixing of different temporalmomentum components is suppressed, so that the "push" towards the origin is weaker, and the "sea" does not contribute to the accumulation of eigenvalues around the origin. Low modes thus originate from the "islands", and this naturally leads to their localisation and to a small spectral density near zero.
There is another important aspect involved in chiral symmetry breaking or restoration. For any finite volume, the spectral density decreases to zero within a sufficiently small distance from the origin. What actually determines the fate of chiral symmetry is how this distance scales with the volume of the system. In the disordered phase, the small eigenvalues originate from small unperturbed modes that occupy a finite fraction of the volume. The "push" caused by mixing of the temporal-momentum components is thus expected to scale with the system size, eventually leading to a finite spectral density at the origin in the thermodynamic limit. In the ordered phase, the small eigenvalues originate from small unperturbed modes localised on the "islands", and the strength of the "push" coming from mixing is expected to correspond to the size of the "islands", which does not scale with the system size, so that the spectral density at the origin remains zero also in the thermodynamic limit.
Before concluding this Section, we want to list a few issues which remain open. We have not checked whether the "chiral" transition in our model is a genuine phase transition in the thermodynamic limit, and if so whether it takes place exactly together with the transition in the spin system. We also have not checked the dependence on the strength h of the U (1) Nc−1 -symmetry-breaking term in the Hamiltonian of the spin system. As this would change the depth of the symmetry-breaking potential, it can in principle affect the nature of the transition in the spin system, and so in turn that of the "chiral" transition. Moreover, we have not investigated in details how our results change when changing the parameter β, which affects the strength of the correlation between time slices. Since the amount of small unperturbed modes is unaffected by such a change, it is reasonable to expect that a largerβ, increasing the correlations and thus reducing the mixing between wave function components corresponding to different temporal momenta in the "sea" region, will make the change in the spectral density when going over to the ordered phase even more dramatic. Similarly, a smallerβ is expected to make this change less dramatic, and for small enougĥ β we expect that "chiral symmetry" is broken also in the ordered phase of the spin model, by continuity with theβ = 0 result discussed above. However, it is not clear if the value of β can affect the nature of the "chiral" transition, i.e., whether it is a true phase transition or a crossover, and the value of β at which it takes place. Nevertheless, we think that the connection between correlation of spatial links across time slices and fate of "chiral symmetry" is strongly supported by our findings.

Conclusions
In this paper we have studied the problems of chiral symmetry breaking and localisation in finite-temperature QCD by looking at the lattice Dirac operator as a random Hamiltonian. We have explicitly recast the staggered Dirac operator at finite temperature in the form of a nonconventional 3D Anderson Hamiltonian ("Dirac-Anderson Hamiltonian"), which describes fermions carrying colour and an extra internal degree of freedom, corresponding to the lattice temporal momenta. The on-site noise is provided by the phases of the Polyakov lines, and their ordering, or the lack thereof, reflects in the distribution of the diagonal entries of the Hamiltonian. The Polyakov lines affect the hopping terms as well. Indeed, in the deconfined phase they induce strong correlations among spatial links on different time slices, in the region where they are aligned with the identity. This in turn weakens the coupling among the components of the wave function corresponding to different lattice temporal momenta on neighbouring sites. In the confined phase, on the other hand, such strong correlations are absent due to the absence of order, and these wave function components can mix effectively. We think that this difference in the hopping terms is essential to explain the accumulation or not of eigenvalues near the origin, and ultimately for the spontaneous breaking or not of chiral symmetry. The other important difference between the two phases concerns the density and spatial distribution of small unperturbed eigenvalues, i.e., of small diagonal entries, which is also caused by the different ordering properties of the Polyakov lines. We think that these two properties of the small unperturbed eigenvalues are also essential in explaining the fate of chiral symmetry. Furthermore, the fact that the small unperturbed eigenvalues appear in localised spatial regions at high temperature leads to the localisation of the low Dirac eigenmodes. This suggests that the confinement/deconfinement transition triggers both the chiral transition and the localisation of the low modes.
To test this picture we have constructed a toy model, made up of a spin system with dynamics similar to that of the Polyakov line phases in QCD; of unitary matrices obeying dynamics analogous to the dynamics of spatial gauge links in the background of fixed Polyakov lines; and of a random Hamiltonian with the same structure as the Dirac-Anderson Hamiltonian discussed above, with on-site noise provided by the spins. This toy model is designed to keep precisely the features of QCD which we believe to be relevant to the phenomena of chiral symmetry breaking (which means here a nonzero spectral density at the origin) and localisation. A numerical study of the toy model, in the simplest case of two colours and two time slices, shows that it indeed displays both these phenomena, with the same qualitative dependence on the ordering of the noise source (spins/Polyakov lines) as in QCD. When the noise source is disordered there is chiral symmetry breaking but no localisation; when the noise source is ordered chiral symmetry is restored and low modes are localised, up to a point in the spectrum ("mobility edge") which is pushed towards larger values as the ordering is increased. If, on the other hand, one artificially removes the correlation between the time slices in the ordered phase, then chiral symmetry remains broken and there is no localisation. Moreover, if the mixing between components corresponding to different time-momenta is artificially removed in the disordered phase, then chiral symmetry is restored and the lowest modes localise. These findings support, and actually suggested, that the correlation among time slices and the related mixing of temporal-momentum components of the quark wave functions play an essential role in the chiral transition and in the appearance of localised modes. The importance of the role played by the small unperturbed eigenvalues is made evident by the results obtained when one imposes on the fermions periodic rather than antiperiodic boundary conditions in the temporal direction. In this case the properly modified toy model again reproduces qualitatively the QCD results, with accumulation of eigenvalues around the origin and no localisation also in the ordered phase. Since the hopping terms are exactly the same as with antiperiodic boundary conditions, one is led to conclude that even with weak mixing between temporal-momentum components one can achieve chiral symmetry breaking, if the density of small unperturbed modes is large enough. Conversely, by artificially increasing the unperturbed modes without touching the hopping terms, one can restore chiral symmetry in the disordered phase and make the low Dirac modes localised, which indicates that strong mixing may be insufficient to accumulate eigenvalues around zero if the density of small unperturbed modes is too low.
The results obtained in the toy model support our expectation that in QCD the fate of chiral symmetry and of localisation are closely related, and furthermore that both depend on the amount of small unperturbed modes and on the mixing of temporal-momentum components of the quark wave functions, therefore ultimately on the distribution of the phases of the Polyakov lines.
The picture discussed here makes no direct reference to topology. As is well known, in the "topological" explanation of chiral symmetry breaking the finite density of near-zero modes originates from fermionic zero modes supported by topological objects, which broaden into a band due to mixing. The localised nature of these modes would also explain localisation at high temperature. In the Dirac-Anderson picture, the "unperturbed" small modes have a different origin, being the eigenmodes of the temporal part of the Dirac operator, and moreover the way they mix (i.e., the nature of the hopping terms) is also expected to play an important role in the accumulation of near-zero modes and in their localisation properties. It is of course well possible that the two pictures are just complementary point of views on the same phenomenon, corresponding to a different way to separate the full Dirac operator into a "free" and an "interaction" part. In light of the close connection between the chiral transition and localisation, and of the central role played by the Polyakov lines in both phenomena, we expect that if this is the case, then the topological objects relevant to chiral symmetry breaking at low temperature, and to localisation at high temperature, would also be relevant to the deconfinement transition. Indeed, there are numerical results pointing to a close relation between localisation and certain topological objects which are expected to play a role in the deconfinement transition [13]. This issue certainly deserves more work. Attention should also be paid to the possible relation between localisation and chiral symmetry restoration, on one side, and "non-topological" approaches to confinement like, e.g., "fluxons" [37,38].
It would be interesting to further investigate in our toy model the behaviour of the spectrum in the vicinity of the phase transition in the spin model. This would clarify if the "chiral transition" seen there is actually a genuine phase transition, and how close it takes place to the magnetic transition. This could provide useful insight in the critical properties at the chiral transition of the "parent" physical system, namely, QCD. Furthermore, one could check how the transition is affected by the strength of the coupling between spatial links on different time slices, and by the depth of the minimum of the spin potential, which are parameters of the model besides the temperature of the spin system.
An obvious extension of this work would be to check the ideas presented above directly with QCD gauge configurations, using the Dirac-Anderson form of the Dirac operator to tweak the hopping terms independently of the underlying Polyakov-line dynamics. While the toy model studied here is quenched, with no backreaction of the quark eigenvalues in the partition function, the main ideas are expected to apply in the presence of dynamical fermions as well.
It would be interesting to try to apply the ideas discussed in this paper in the case when a constant (Abelian) magnetic field is turned on. This could shed some light on the issue of (inverse) magnetic catalysis of the quark condensate [39,40]. Another interesting application would be to the case of nonzero imaginary chemical potential, already very briefly discussed here.
Another interesting testing ground for the proposed mechanism is the explanation of the separate occurrence of the deconfinement and chiral transitions in SU (3) gauge theory with adjoint fermions [41]. Since the derivation of the Dirac-Anderson Hamiltonian did not use in an essential way that we were considering fundamental fermions, the same form holds for adjoint fermions, replacing gauge links with their adjoint counterpart, and the N c phases of the fundamental Polyakov line with the N 2 c − 1 phases of the adjoint Polyakov line. In conclusion, we believe that the "Dirac-Anderson" approach of the present paper to the study of the quark eigenvalues and eigenfunctions can lead to a better understanding of the phase structure of QCD and related theories.

A Properties of the hopping term
In this Appendix we describe in some detail the properties of the QCD Dirac-Anderson "Hamiltonian", and in particular of its hopping terms, V ±j . First of all, notice that as it should. From now on we will often omit matrix indices, so we remind the reader that U ±j (t, x) has only colour indices, while D( x) and V ±j ( x) have both colour and temporal-momentum indices. The identity in these spaces will be denoted by 1 c and 1 tm , respectively. SinceŨ ±j (t, x) are unitary matrices, and since V ±j ( x) is the Fourier transform with respect to time of the unitary matrixŪ ±j ( x), we have that V ±j ( x) is also unimodular up to a sign. If we choose one and the same convention for the phases of the local Polyakov lines, i.e., we fix q( x) = q ∀ x, then V ±j ( x) ∈ SU (N c × N T ), ∀ x, j. However, one can choose different phase conventions at different spatial points, and still obtain the same physical results. We will return below to this issue. We observe also the following cyclicity property of V ±j ( x), [V ±j ( x)] a(k+n) N T ,b(l+n) N T = [V ±j ( x)] ak,bl , ∀n ∈ Z . (A.7) One can easily verify that this property is preserved under multiplications, so V ±j ( x) belong to the "(N T × N T )-block cyclic" subgroup of SU (N c × N T ).
Let us now return to the issue of the choice of phase conventions. All φ a ( x) are defined modulo 2π, so after a redefinition φ a ( x) → φ a ( x) + 2πq a ( x) with q a ( x) ∈ Z one should obtain equivalent results. We will denote quantities after the redefinition with the superscript {q}. so that the temporal-momentum indices of D( x) and V ±j ( x) undergo a colour-dependent cyclic permutation k → k + q a ( x) mod N T . This can be written formally using the permutation matrices Π (n) , Π where (η 4 Hη 4 )( x, y) = η 4 ( x)H( x, y)η 4 ( y). We conclude that Q ≡ η 4Z satisfies {Q, H} = 0, which implies that the spectrum is symmetric with respect to λ = 0, as it should be for staggered fermions. One final remark is in order concerning the case of gauge group SU (2). In this case one has φ 1 = −φ 2 ≡ φ. It is known that in this case the Dirac operator has an antiunitary symmetry T with T 2 = −1 [2]. In the new basis, taking the complex conjugate of H has the effect of (i) exchanging the indices k, l in temporal-momentum space, (ii) changing the sign of the phases φ appearing inŨ [Eq. (2.12)] and taking the complex conjugate of the SU (2) matrices U (td) ±j , and (iii) changing the overall sign of the hopping terms. Point (ii) can be "undone" by taking the matrix conjugate with respect to σ 2 in colour space (this remains true also in the presence of nontrivial phases φ, as can be directly checked); this also leads to the diagonal element k being switched with the element N T 2 − 1 − k. Indeed, matrix conjugation by σ 2 exchanges the diagonal terms corresponding to φ and −φ, and this is equivalent to switching temporal-momentum components, since (A. 15) This corresponds to a permutation Π of the temporal-momentum components defined so that k → N T 2 − 1 − k mod N T (notice Π 2 = 1). Since the hopping term depends on k, l only through k − l, one has l − k = N T 2 − 1 − k − N T 2 − 1 − l , and so by applying Π we undo both point (i) and the above-mentioned effect on the diagonal term. Finally, taking the matrix conjugate with respect to η 4 in spatial-coordinate space we undo point (iii). All in all, T = η 4 Πσ 2 K, with K the complex conjugation, is an antiunitary symmetry of the Hamiltonian with T 2 = −σ 2 2 = −1.