{'fraction_non_alphanumeric': 0.03363683030286656, 'fraction_numerical': 0.009510172101694, 'mean_word_length': 4.661364331599878, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Polyakov loops and spectral properties of the staggered Dirac operator

We study the spectrum of the staggered Dirac operator in SU(2) gauge fields close to the free limit, for both the fundamental and the adjoint representation. Numerically we find a characteristic cluster structure with spacings of adjacent levels separating into three scales. We derive an analytical formula which explains the emergence of these different spectral scales. The behavior on the two coarser scales is determined by the lattice geometry and the Polyakov loops, respectively. Furthermore, we analyze the spectral statistics on all three scales, comparing to predictions from random matrix theory.


I. INTRODUCTION
One of the goals of quantum chromodynamics (QCD), the fundamental theory of quarks and gluons, is the calculation of the hadronic mass spectrum from first principles. For this purpose, a nonperturbative regularization of QCD can be formulated on a space-time lattice, which makes the theory mathematically well defined and accessible for numerical simulations. In order to describe fermionic fields, one needs a discretized Dirac operator, which can be constructed in different ways. In particular, the staggered or Kogut-Susskind Dirac operator [1] is widely used as it is computationally cheaper than other options. (We shall not enter the debate of the rooting issue [2,3] here.) During the past 15 years it has been shown that some universal properties of the QCD Dirac spectrum can be described by a version of random matrix theory (RMT) [4,5] which incorporates the chiral symmetry of the massless Dirac operator and is accordingly called chiral random matrix theory (chRMT) [6,7]. In chRMT, one models the Dirac operator by a block off-diagonal matrix with random entries, respecting the global symmetries of the massless Dirac operator. The path integral over the gauge fields is then replaced by averaging over an ensemble of such matrices. Chiral RMT reproduces low-energy sum rules [8] and yields accurate predictions for, e.g., the microscopic spectral density [9], the distribution of the smallest eigenvalues [10,11], level spacing distributions, and other short-range spectral correlations. Predictions from chRMT have been successfully compared with numerical results from lattice gauge theory in many different settings [12,13,14,15,16,17,18,19,20,21].
When using RMT to predict spectral properties of complex quantum systems, it is essential that the random matrices have the same antiunitary symmetries as the system to be modeled, see, e.g., Ref. [22]. One distinguishes different symmetry classes and corresponding ensembles of random matrices, labeled by their Dyson index β D = 1, 2, or 4 [4].
The staggered Dirac operator on the lattice exhibits the peculiar feature that its symmetry properties can be different from that of the continuum theory. In particular, for gauge group SU (2) with fermions in the fundamental representation and for gauge group SU(N ) with fermions in the adjoint representation, the β D = 1 and β D = 4 cases are interchanged compared to the continuum Dirac operator (see below). This implies that the spectral properties of the staggered operator are different from those of the continuum Dirac operator. Accordingly, a transition is expected to take place in the continuum limit. The first indication of such a transition has been reported in Ref. [23].
In this work, we study a related but different case, namely, the staggered Dirac operator for fermions in the fundamental and adjoint representation of SU(2) in the free limit. This limit is approached by increasing the Wilson gauge action parameter β = 4/g 2 at fixed (or mildly varying) lattice size, i.e., the lattice spacing, and thus the physical volume, shrinks to zero.
In the free limit the dynamics becomes integrable and therefore one no longer expects RMT statistics. Consequently, in this limit we do not expect the same transition between RMT symmetry classes as in the continuum limit described above. Instead, generically, the eigenvalues should be uncorrelated, like numbers drawn from a Poisson process [24,25]. It turns out that in the case we are studying, the situation is somewhat more complicated: Close to the free limit it is possible to disentangle three different scales that appear in the separations between the eigenvalues.
First, the eigenvalues form well-separated plateaux centered at the eigenvalues of the free staggered operator.
Second, we observe an internal structure of the plateaux: The eigenvalues arrange themselves in clusters of eight eigenvalues each. We will show that, for each configuration, the location of these clusters can be predicted from the knowledge of only four real variables, i.e., the averaged traced Polyakov loops in the four lattice directions.
Third, on the finest scale, the eigenvalue fluctuations within clusters can be described in terms of chRMT. We will present numerical results for the level spacing densities which agree with the chRMT prediction for all values of the Wilson gauge action parameter β.
This article is organized as follows. In Sec. II we briefly introduce some basic notions of random matrix theory and discuss the different antiunitary symmetries which are relevant for the Dirac operator in the continuum and for the staggered Dirac operator on the lattice. In Sec. III we derive an analytical prediction for the staggered Dirac eigenvalues in certain gauge field configurations close to the free limit. Section IV begins with a comparison of numerical data from lattice simulations to our analytical prediction and continues with an analysis of spectral statistics on all three scales; i.e., we study the distribution of level spacings within clusters, between clusters, and between plateaux. Wherever appropriate we compare to the predictions from chRMT and from a Poisson process. Throughout Secs. II to IV we discuss the fundamental and adjoint representation side by side. In Sec. V we conclude with a discussion of the implications of our findings. Preliminary results of this study have been presented in Ref. [26].

II. ANTIUNITARY SYMMETRIES AND RMT ENSEMBLES
The notion of universality, as commonly used in the context of RMT and the analysis of spectra of complex quantum systems, means that spectral statistics can be described by an appropriate ensemble of random matrices, which shares the symmetries of the system under consideration. Depending on the presence of antiunitary symmetries, the entries of matrices of the ensemble have to be either real, complex, or quaternion real. The associated Dyson indices are β D = 1, 2, or 4, respectively.
In particular, RMT yields a prediction for the universal quantity P (s), the probability density for the unfolded nearest-neighbor spacings s (see Sec. IV C 1 for a discussion of unfolding). This prediction is well approximated by the Wigner surmise, The Wigner surmise is the level spacing density for ensembles of 2 × 2 matrices with Dyson index β D . For the QCD Dirac operator D the RMT description is formulated in terms of matrices which reflect the chiral, flavor, and antiunitary symmetries of D [27]. For gauge group SU(N ) with N ≥ 3 colors and fermions in the fundamental representation, D generically has complex elements and does not commute with any antiunitary operator. Accordingly, universal spectral correlations are described by the chiral unitary ensemble (chUE), labeled by the Dyson index β D = 2.
However, the Dirac operator enjoys invariance with respect to an antiunitary transformation if the fermions are either in the fundamental representation of the gauge group SU (2), or in the adjoint representation of SU(N ) with N ≥ 2 arbitrary. These two cases are discussed in more detail in the following.

A. Fundamental representation
In a nutshell, the antiunitary invariance of the Dirac operator with SU(2) gauge fields and fermions in the fundamental representation, i.e., QCD with two colors, is based on the fact that the generators are τ a /2, and the Pauli matrices τ a possess the following complex conjugation property: (2. 3) The antiunitary symmetry operator, however, is realized in different ways in the continuum and in the lattice formulation with staggered fermions.

Continuum
In the continuum, the anti-Hermitian massless Dirac operator is defined as with T a = τ a /2 for the fundamental representation of SU (2). It anticommutes with the chirality operator γ 5 , and therefore its nonzero eigenvalues, which are purely imaginary, come in complex conjugate pairs. Using Eq. where C = γ 2 γ 4 is the charge conjugation matrix, and K denotes complex conjugation (in the position representation). Note that Cγ 5 acts on the spinor indices, whereas τ 2 acts in color space. Since (Cγ 5 τ 2 K) 2 = ½, it follows that D can be made real by a basis transformation that does not depend on the gauge configuration [22,28]. Accordingly, the RMT description for two-color QCD and fundamental fermions is formulated in terms of the chiral orthogonal ensemble (chOE), characterized by β D = 1.

Lattice
The staggered or Kogut-Susskind Dirac operator for a hypercubic lattice of finite spacing a in d dimensions is given by with η µ (x) = (−1) P ν<µ xν and U µ (x) ∈ SU(2). On the lattice the operator S = δ x,y (−1) P d ν=1 xν plays the same role as γ 5 does in the continuum: Since {D KS , S} = 0, the eigenvalues of D KS also come in (purely imaginary) complex conjugate pairs.
The staggered Dirac operator, which is widely used in numerical simulations (because it maintains a remnant of chiral symmetry, partially solves the doubling problem, and is computationally cheaper than other lattice Dirac operators) exhibits the peculiar feature that its antiunitary symmetries are different from those of the continuum Dirac operator [12,29,30]. Because of Eq. (2.3) the links obey U µ (x) = τ 2 U * µ (x)τ 2 . Since the γ matrices have been replaced by real numbers, η µ (x), no charge conjugation is required in order to compensate for the complex conjugation. Therefore, D KS is invariant under the following antiunitary symmetry: As (τ 2 K) 2 = −½, it follows that D KS can always be written as a quaternion real matrix [22]. Hence, the chiral symplectic ensemble (chSE), with β D = 4, is used to describe its universal properties. Another consequence of invariance with respect to an antiunitary transformation with square −½ is Kramers' degeneracy, see, e.g., [22,31]; i.e., all eigenvalues have (at least) multiplicity two. This degeneracy has to be removed by hand before one discusses spectral correlations.

B. Adjoint representation
In this case, the antiunitary symmetries are determined by the purely imaginary nature of the matrix elements of the generators in the adjoint representation (being the structure constants) of the gauge group SU(N ). As in the previous case, the way this symmetry is realized in the continuum and for the staggered Dirac operator on the lattice is different.

Continuum
The purely imaginary adjoint generators induce a real covariant derivative in the continuum Dirac operator of Eq. (2.4), which is therefore invariant under the following antiunitary symmetry: (2.8) Since (Cγ 5 K) 2 = −½, it is possible to recast D into real quaternionic form. This implies that the Dirac spectrum can be described in terms of the chiral symplectic ensemble (chSE), labeled by β D = 4.

Lattice
For fermions in the adjoint representation, the staggered Dirac operator D KS is explicitly real, because such are the U µ (x) link matrices appearing on the right-hand side of Eq. (2.6), which now take values in the adjoint representation of SU (2). This implies that the appropriate RMT ensemble is the chOE, characterized by the Dyson index β D = 1.

III. SPECTRAL PREDICTIONS FROM POLYAKOV LOOPS
The spectral properties of the Dirac operator are relevant to the two major unsolved nonperturbative problems in QCD: confinement and chiral symmetry breaking. On the one hand, the eigenvalue density is related to the chiral condensate, the order parameter for the chiral phase transition, by the Banks-Casher relation [32]. On the other hand, the average value of the Polyakov loop, which, in the quenched case, is an order parameter for the confinement-deconfinement transition, can be expressed through sums of the eigenvalues of the lattice Dirac operator with different boundary conditions [33,34,35,36].
On commensurate lattices, i.e., on lattices for which the numbers L µ /2, µ = 1, . . . , d, are rationally dependent (the staggered Dirac operator only makes sense on lattices with even L µ ), the spectrum of D KS in the trivial vacuum, i.e., for the field configuration with all U µ (x) equal to unity, is highly degenerate. Our analysis of the spectrum close to the free limit shows that the way in which this degeneracy is (partially) lifted can be expressed in terms of the traced and averaged Polyakov loops in all directions, where the x average is over the whole lattice. For each configuration, this effect can be predicted by calculating the spectrum obtained from a vacuum configuration with uniform link variables in each direction, taking values in an Abelian subgroup of the gauge group, and yielding the same averaged Polyakov loops as the original configuration.

A. Fundamental representation
For the trivial gauge vacuum (all links set to ½ or any gauge transform thereof), the eigenvalues of the staggered Dirac operator read with the wave numbers taking integer values 0 ≤ k µ < L µ /2 and c µ = 0 c µ = 1 2 for (anti)periodic boundary conditions for the fermionic wave function in direction µ. However, the system also admits vacua in different center sectors, which can be labeled by the traced Polyakov loops P µ = ±1 [for SU (2)]. The sign of the latter can always be inverted by a multiplication of all the links in the µ direction in a given fixed-x µ slab by −½. In the Dirac operator this can be compensated for by switching from periodic to antiperiodic boundary conditions (or vice versa) in the direction µ. Therefore, when we consider a vacuum where the Polyakov loop in direction µ is −1, we can equivalently set c µ = 1 2 (c µ = 0) for (anti)periodic boundary conditions of the Dirac operator in direction µ.
For a generic lattice in d dimensions, the number of possible free spectra is thus equal to the number of allowed topological sectors for the vacuum, i.e., 2 d . When at least two lattice extensions are equal, the number of possible free spectra is reduced. For lattices with L µ = L ν ∀ µ, ν = 1, . . . , d, the different patterns are labeled by the number of c µ values that are equal to 1 2 , thus yielding d + 1 inequivalent possibilities.
Configurations close to the free limit are expected to approach (modulo gauge transformations) one of the possible free vacua. This is confirmed by our lattice simulations, where we find the distribution of P µ to be peaked at ±1. The corresponding free vacua can then easily be identified by the sign of the averaged traced Polyakov loops P µ in the various directions. Accordingly, Eq. (3.2) provides a first approximation to the observed spectrum of the staggered Dirac operator.
This prediction can be refined as follows. For a given configuration, let us introduce a configuration built from uniform links U µ (x) ≡ U µ in each direction, taking values in an Abelian subgroup of SU(2) (for instance the diagonal one), and yielding the same averaged traced Polyakov loops as the original configuration. Since for uniform and commuting links all plaquettes are equal to unity these configurations may also be called vacuum configurations.
For these vacuum configurations, the gauge transformation can be used in order to trivialize all Polyakov loops to P µ = 1 at the expense of introducing periodicity only up to g(x) (with x µ = L µ ). The latter equals the original Polyakov loop P µ and (in the diagonal subgroup) comprises two opposite phases, which behave like the constants c µ (this is a generalization of our previous argument that P µ = −1 can be absorbed by switching c µ between 0 and 1 2 ). The spectrum of the staggered Dirac operator in such a vacuum configuration is given by Now c µ is again fixed to 0 ( 1 2 ) for (anti)periodic boundary conditions in the direction µ. Our expectation is that, close to the free limit, Eq. (3.4) provides a better approximation than Eq. (3.2) to the observed spectrum of D KS .
The eigenvalues in Eq. (3.4) have a multiplicity of 2 d , which we derive in the following paragraph. For dimension d = 4, this implies an eightfold degeneracy in addition to Kramers' degeneracy. A small perturbation of the vacuum generically lifts this eightfold degeneracy but not Kramers' degeneracy, which remains exact. This mechanism gives rise to what we will call clusters of eight eigenvalues below.
In order to explain how the multiplicity of 2 d comes about, we have to digress to a sketch of the derivation of Eq. (3.4). The spectrum is most easily obtained by looking at the square D 2 KS of the staggered Dirac operator and applying it to plane waves. After some algebra, one finds the eigenvalues of D 2 KS to be given by we can restrict the wave numbers to 0 ≤ k µ ≤ L µ /2 − 1 and assign a multiplicity of 2 d to each eigenvalue. An additional multiplicity factor of 2 arises from the color degeneracy, yielding an overall multiplicity of 2 d+1 . Because of the symmetry of the spectrum of D KS about λ = 0 (which arose due to {D KS , S} = 0, see Sec. II A 2 above), the eigenvalues of D KS are given by the positive and negative square roots of Λ n with each eigenvalue λ ± n having half the multiplicity of the corresponding Λ n , i.e., 2 d .

B. Adjoint representation
The situation is similar for the adjoint representation of SU (2). Since the latter is insensitive to the group center, all trivial vacua are equivalent to the configuration with all links equal to unity.
The construction above can be repeated by considering a vacuum configuration built from link matrices of the form where α µ is related to the traced Polyakov loop P µ by L µ α µ = arccos((P µ − 1)/2). The analog of Eq. (3.4) now where n ∈ {−1, 0, 1} and 0 ≤ k µ < L µ /2. In contrast to the situation for the fundamental representation, Eq. (3.7) predicts that one-third of the eigenvalues remains unchanged, i.e., they are identical to the eigenvalues in the trivial vacuum (the configuration with all links equal to unity). An analysis analogous to that for the fundamental representation shows that now the multiplicity is 2 d−1 (as opposed to 2 d above). However, here we have no Kramers degeneracy, and thus a small perturbation of the vacuum again gives rise to clusters of eight eigenvalues in dimension d = 4.

A. Simulation details
Our numerical results are obtained from sets of quenched SU(2) configurations generated using the Wilson gauge action. The simulation algorithm is based on a combination of Metropolis and over-relaxation; center rotations to explore different topological sectors are implemented as well.
We obtain the full spectrum of the staggered Dirac operator from ensembles of configurations on hyper-cubic, isotropic lattices with volumes ranging from V = 4 4 to 16 4 . The spectrum is evaluated using the Cullum-Willoughby version of Lanczos' algorithm [37]; periodic (antiperiodic) boundary conditions are assumed in the spatial (temporal) directions.
For each lattice volume V and for each value of the Wilson gauge action parameter β, our analysis of the full staggered spectrum is based on a number of thermalized and uncorrelated configurations between a few tens and a few thousands. For each configuration the number of distinct eigenvalues with positive imaginary part is V /2 for the fundamental representation and 3V /2 for the adjoint representation. Hence, our data for each pair (V, β) typically contains on the order of one million distinct eigenvalues.
Furthermore, we also generate gauge configurations on much larger lattices (up to 34 × 38 × 46 × 58) in order to investigate the distribution of eigenvalues from Eqs. (3.4) and (3.7). On these lattice we do not diagonalize the Dirac operator but calculate only the averaged Polyakov loops.

B. Separation of spectral scales
When β is increased to large values, on a lattice with a fixed number of sites in each of the four directions, the spectrum of D KS shows structure on three different scales, see Figs. 1 and 2 for examples with fundamental and adjoint fermions, respectively.

Plateaux
The coarse structure is given by the free limit, i.e., the eigenvalues approach the values of Eq. As discussed in Sec. III A, for the fundamental representation there are five different plateau structures. We find all these classes in our simulated configurations and show representatives for each class in Fig. 1. (The spectra shown are the first representatives of each class which we come across in the Monte Carlo history.) For adjoint fermions the plateau structure is uniquely determined by the lattice sizes (and the boundary conditions) since the adjoint representation is center-blind, as discussed in Sec. III B. Thus, configurations with different signs of P µ lead to the same plateau structure, see Fig. 2.

Clusters
A closer look at the staggered spectra reveals that the distribution of eigenvalues inside a given plateau shows additional structure. The eigenvalues are grouped in clusters of eight. When the free limit is approached, the typical separation between nearest clusters within the same plateau is smaller than the separation between different plateaux, but larger than the separation of eigenvalues within each cluster.
The position of the clusters can be predicted by Eqs. (3.4) and (3.7) for the fundamental and adjoint representation, respectively. In Figs. 1 and 2 these predictions are indicated by solid red lines. Thus, the cluster structure is determined by the traces of the averaged Polyakov loops (and the lattice size and boundary conditions). For the adjoint representation, we see that onethird of the eigenvalues forming a plateau does not split into clusters, but stays close to the plateau levels. This is also in agreement with our discussion of Eq. (3.7).
Recall that, approaching the free limit, the distribution of the traced Polyakov loops becomes peaked at ±1, corresponding to the center elements of SU(2). This a posteriori justifies the approximation of the staggered spectrum by the plateaux.
The observation that each cluster contains eight eigenvalues reflects the fact that the eigenvalues of Eqs.  Since the simulated configurations close to the free limit, which we approximate by these vacuum configurations, do not have exactly uniform and commuting links, the eightfold degeneracy is lifted.
In the lower right panel of Fig. 1 we also show an example obtained from simulations on a 6 4 lattice. We see that for larger lattices the expected patterns get more and more complicated, but the agreement persists with just the four parameters P µ of Eq. (3.4) determining the complete cluster structure.

Unfolding
Before one can discuss spectral statistics and compare to, e.g., predictions from RMT, the spectra have to be unfolded, i.e., the (imaginary parts of the) eigenvalues have to be rescaled such that the mean separation between adjacent levels is unity. This can be achieved by defining a new spectrum {x n } with x n := N (λ n ), where N (λ) denotes the mean integrated spectral density, i.e., N (λ) is a smooth function satisfying N (λ) ≈ #{λ n ≤ λ} in some approximate or asymptotic sense. In low-dimensional quantum chaos, one usually has an analytical formula for N (λ), provided by Weyl's law, see, e.g., [22]. For the QCD Dirac operator first steps towards an equivalent asymptotics have been discussed in [38]. Until now, however, this approach does not provide a formula which can be used for unfolding spectra in lattice QCD.
If one has no analytical prediction for the mean spectral density, it has to be extracted from the data themselves. The latter can be done by averaging over several spectra, which is known as ensemble unfolding, thus yielding one mean density for a whole ensemble of spectra. In contrast, fitting an ansatz to the spectral density of an individual spectrum or extracting its mean density by a moving average in λ is known as configuration unfolding. The resulting mean densities will in general differ from each other. This ambiguity can make it difficult to extract reliable information on long-range spectral correlations, see the detailed discussion of different unfolding methods for lattice Dirac spectra in Ref. [39]. The statistical function we are interested in is always the density P (s) of nearest-neighbor spacings. Since P (s) measures short-range spectral correlations, changing the unfolding method will not impair our results, as long as we make sure that each unfolding procedure is stable and consistent.
As we want to discuss spectral statistics for very different values of the gauge action parameter β and on different scales, we are forced to employ different unfolding methods, each one of them tailored to fit the requirements of the particular situation.
For large β the spectra show different scales, see Sec. IV B. Therefore, we unfold separately on each scale.
The different plateau and cluster structures forbid methods of ensemble unfolding, and thus we employ configuration unfolding. The methods we use are all variants of what is called local unfolding in Ref. [39].
When studying spacings within clusters we divide spacings of adjacent levels by the mean level spacing within the cluster. The latter we calculate as the difference between the largest and the smallest level divided by seven. The unfolded spacings then have unit mean as required.
For spectral statistics between clusters within a plateau, we proceed analogously. We divide spacings between adjacent clusters by the difference between the position of the largest and the smallest cluster divided by the number of clusters within the plateau minus one.
Finally, when we want to study spacings between plateaux, we adapt the previous methods as follows. We divide the difference between the position of plateau j + 1 and plateau j by the difference between the positions of plateaux j + 5 and j − 4 divided by the number of plateau spacings in this range, i.e., by nine.
For small β the spectra do not show different scales. Thus, complete spectra could be unfolded in one go, and both configuration and ensemble unfolding are admissible. We experimented with the ensemble unfolding described in Ref. [18], which works reliably for small β. However, when we increase β the method has to break down, due to the different plateau structures which begin to emerge. In order to have an unfolding method which does not exclude an intermediate β-range but instead allows us to go smoothly from small to large β, we decided to unfold spectra for small β in the same way as we unfold spacings within clusters for large β. Effectively, this means that we neglect every eighth spacing and unfold the remaining spacings on the scale of the neighboring eight eigenvalues.

Level spacings for small β
For small values of β, i.e., far away from the free and continuum limits, the level spacing densities have to be compared to the predictions reported in Secs. II A 2 and II B 2. Figure 3 shows the level spacing densities obtained from spectra of the staggered Dirac operator in the fundamental representation of SU (2). The numerical results on all our lattices are consistent with the chSE, in contrast to the chOE, which would be expected in the continuum.
Likewise, in Fig. 4 we show level spacing densities obtained from staggered Dirac spectra in the adjoint representation of SU (2), which now are consistent with the chOE as expected. Again this is in contrast to the continuum situation in which one should have chSE statistics.
This behavior of the staggered Dirac operator in SU(2) gauge fields, having spectra in different universality classes than the corresponding continuum operator, has been observed earlier [12,14,17].

Level spacings for large β: Approaching the free limit
For large β, plateaux and clusters emerge. Therefore, we discuss spectral statistics separately on the different scales.
(a) Spacings within clusters. Generically, the eightfold degeneracy of the levels predicted by Eqs. (3.4) and (3.7) is lifted by the nonuniformity of the gauge field configuration. For this reason we observe small clusters of eight eigenvalues.
Close to the free limit we can think of treating the nonuniformity of the gauge field configuration as a small perturbation of the corresponding vacuum configuration, i.e., of the configuration with uniform and commuting links leading to the same Polyakov loops. A cluster then arises by diagonalizing an 8 × 8 matrix, the perturbation restricted to the subspace corresponding to a degenerate eigenvalue of Eq. (3.4) or Eq. (3.7). This matrix inherits symmetries and effective randomness of the perturbation, i.e., of the gauge field part of the staggered Dirac operator. Therefore, we expect the distribution of level spacings to follow the same chRMT prediction as for small β.
Our data for fermions in the fundamental representation of SU(2) confirm this expectation. As Fig. 5 shows, the level spacing density within each cluster is consistent with the prediction from the chSE. Note that this is true over a very large range of β values.
Also in the adjoint representation the spacings within clusters follow the same pattern as the spacings for small β, in this case leading to a chOE distribution, see Fig. 6. Note that for the adjoint representation, the clusters at the ends of a plateau show some nongeneric features which would require a more sophisticated unfolding procedure. We avoid this complication by restricting the analysis to cluster spacings in the central part of each plateau.
These results confirm once more that the spectral properties of the staggered Dirac operator are different from those of the continuum operator. Moreover, in the sense described above, our analysis demonstrates that this discrepancy persists when approaching the free limit.
(b) Spacings between clusters. For the fundamental representation the cluster positions are, to a good approximation, described by Eq. (3.4). Both Eq. (3.2) and Eq. (3.4) describe spectra of integrable systems, and thus one generically expects uncorrelated levels, like for a Poisson process [24,25].
However, on commensurate lattices, i.e., on lattices with rationally dependent L µ /2, Eq. (3.2) predicts many accidental degeneracies which would lead to nongeneric spectral statistics. Generic behavior can in principle be restored in two different ways: On the one hand by changing to lattices with rationally independent extensions L µ /2, and on the other hand by introducing additional phase shifts, like the Polyakov loops do in Eq. (3.4). However, as we pointed out earlier, close to the free limit the distribution of averaged traced Polyakov loops is peaked at ±1. In Fig. 3 of Ref. [26] it has been observed that these almost equal phase shifts in each direction are not able to restore generic behavior.
Therefore, we choose to study a large incommensurate lattice with extensions L 1 × L 2 × L 3 × L 4 = 34 × 38 × 46 × 58. We refrain from diagonalizing the Dirac operator on this lattice but instead calculate the averaged traced Polyakov loops and determine the approximate cluster spectrum from Eq. (3.4). Then the density of spacings between different clusters within the same plateau agrees well with the prediction from a Poisson process, P Poisson (s) = e −s , see Fig. 7 (left). The same holds true for fermions in the adjoint representation with the approximate cluster spectrum determined from Eq. (3.7), see Fig. 7 (middle).
Note that we analyze the spacing distribution for the spectra (3.4) and (3.7) but not the spacings between the actual cluster positions which one would obtain by averaging over all eigenvalues belonging to a given cluster. As can be seen in Figs. 1 and 2 these actual positions differ slightly from the predicted values. Whether or not these differences could lead to a deviation from Poisson behavior is an open question.
(c) Spacings between plateaux. The positions of the plateaux can be approximately described by the free spectrum, Eq. (3.2). As in the case of spacings between clusters, which we discussed above, we thus once more study the spectrum of an integrable system. Again accidental degeneracies lead to nongeneric statistics on commensurate lattices. Generic behavior can be restored by switching to an incommensurate lattice, which we demonstrate in Fig. 7 (right). Level spacing distributions like for a Poisson process were predicted and observed earlier for incommensurate lattices [40,41].
We add the same note of caution as for the spacings between clusters. Our analysis concerns the level spacings of the spectrum (3.2) but not the spacings between the actual plateau positions, obtained by averaging over all eigenvalues within a given plateau.
Let us also point out that, in spite of the strong similarities between the data describing the cluster spacings and the plateau spacings in Fig. 7, they describe correlations on spectral scales that typically differ by an order of magnitude.

V. CONCLUSIONS AND OUTLOOK
We have studied the spectral properties of the staggered Dirac operator D KS when approaching the free limit for gauge group SU(2), both in the fundamental and in the adjoint representation. With SU(2) gauge fields the staggered Dirac operator on the lattice belongs to a different symmetry class than the corresponding continuum operator.
Our numerical analysis revealed that, when the free  limit is approached, the spectrum of the staggered Dirac operator shows structure on three well-separated scales. The behavior on the two coarser scales is characterized by the formation of plateaux and clusters. We have shown that for a given gauge field configuration the positions of plateaux and clusters can be predicted analytically. To this end we have constructed a vacuum configuration with uniform and commuting links in each direction, chosen such that they reproduce the averaged traced Polyakov loops of the original configuration. The Dirac spectrum in this vacuum configuration yields a good approximation to the plateau and cluster positions. In turn, plateaux and clusters are determined by lattice geometry, boundary conditions, and Polyakov loops alone.
Our model also predicts a systematic degeneracy of the staggered spectra in nontrivial vacuum configurations: All eigenvalues have a multiplicity of eight (in addition to Kramers' degeneracy in the fundamental representation). For typical gauge field configurations this degeneracy is lifted leading to the formation of clusters of eight eigen-values, which we observe numerically.
We have analyzed spectral correlations on all three scales in terms of the distribution of spacings between adjacent eigenvalues, clusters, and plateaux. Spacings between approximate plateau and cluster positions are uncorrelated as for a Poisson process, whereas level spacings on the finest scale, i.e., within clusters, follow the chRMT predictions. For the latter the symmetry class is always that of the staggered operator and never that of the continuum operator, even for very large β, i.e., close to the free limit.
Finally, we briefly comment on the possibility to recover the symmetry class of the continuum operator in spectra of the staggered Dirac operator in SU(2) gauge fields. We restrict ourselves to the fundamental representation. For every finite lattice spacing, the antiunitary symmetries of the D KS operator are different from those of the continuum Dirac operator. When approaching the continuum limit, i.e., when increasing β at fixed physical volume, the eigenvalues should form near-degenerate quartets (more precisely, near-degenerate pairs of exactly Kramers-degenerate pairs), signaling the suppression of unphysical taste-changing interactions of the staggered operator. This has been seen for the first time in Ref. [23] using highly improved staggered fermions. In [23] it was also shown that the distribution of the smallest eigenvalue makes a transition from chSE to chOE, i.e., from the symmetry class of the lattice operator to that of the continuum operator. Related phenomena have been observed for improved staggered fermions in the fundamental representation of SU(3) [42,43,44]. We expect that the distribution of spacings between eigenvalues will also show a similar transition to the continuum behavior. Our present findings do not contradict such an expectation as our study corresponds to a different physical setting: In-stead of the continuum limit we have investigated the behavior in the free limit, where RMT applies only at the finest scale of eigenvalue fluctuations.