Eigenstate Thermalization in 2+1 dimensional SU(2) Lattice Gauge Theory

,


I. INTRODUCTION
The question of how thermal behavior emerges in systems governed by the strong interaction has a long history, dating back to Fermi's statistical [1] and Landau's hydrodynamical [2] models of multiparticle production and ultimately culminating in Hagedorn's statistical bootstrap model [3].While these models posited microcanonical or canonical thermalization, they did not explain its origins.More recent attempts to show that systems whose dynamics is governed by nonabelian gauge theory thermalize rapidly are based on three approaches: (1) kinetic theory applicable to weak coupling [4][5][6]; (2) semiclassical dynamics in holographic duals of strongly coupled supersymmetric gauge theories [7][8][9]; (3) classical simulations of Hamiltonian SU (2) and SU(3) lattice gauge theories [10,11].All these approaches have their limitations and do not fully answer the question.
The numerically established fact that SU(2) lattice gauge theory is extensively chaotic at the classical level [12] suggests that the theory also exhibits quantum chaos when quantized.It is generally believed that chaotic quantum systems satisfy the Eigenstate Thermalization Hypothesis (ETH) [13,14], which posits that the matrix elements of generic operators in the energy eigenstate basis are given by: where E = (E α + E β )/2 and ω = E α − E β , the R αβ vary radically with zero mean and unit variance, forming a Gaussian distribution in large systems, ⟨A⟩ mc (E) is the microcanonical expectation value of A, f A (E, ω) is the spectral response function of the operator, and the exponential prefactor accounts for the energy level density in terms of the microcanonical entropy S(E) (see [15,16] for reviews).
It is often stated (see e.g., [15]) that ETH behavior is a generalization of random matrix theory (RMT).It is worth noting that the matrix R αβ may not be a random matrix due to the constraints of the underlying theory.Furthermore, it is important to recognize that RMT expresses statistical properties that may apply to any type of mathematical object, like the zeros of the Riemann zeta function, but the ETH relation Eq. ( 1) expresses a dynamical property that is encoded in the spectral function f (E, ω).
While Eq. ( 1) is at the center of a huge body of work, it raises many important questions that are still under debate.Some of these are: • For which operators A does the defining ETH relation Eq. ( 1) apply?Most likely, this is only the case for certain classes of "physical" operators, such as local operators or sufficiently averaged ones.For nonabelian gauge theories a restriction to gauge invariant and multiplicatively renormalizable operators also seems a natural requirement.As such constrained operators first come to mind and thus are typically studied in numerical investigations, it could be that such numerical studies provide only a biased perspective and should not be blindly assumed to generalize.
• What is the thermalization time for a given operator A? It was shown that if the ETH is satisfied in the weak sense (i.e., most eigenstates satisfy the ETH), a two-point correlation function factorizes at late time [17].Thermalization happens rapidly for most observables [18].However, the thermalization time in general depends on the initial state.
For some pathological initial states with narrow energy spread, it can be arbitrarily long [19,20].
• To which precision must Eq. ( 1) be fulfilled to reach a certain degree of thermalization?For example, if equilibrium is reached at late time, only states in a small ω window are probed.
At present we are unable to answer such questions rigorously, making numerical studies of specific systems the most promising avenue to pursue.Comparing the results with phenomenology allows to check the plausibility of the assumptions made and helps to develop a better understanding of the systems studied, without reaching mathematical rigor.Our contribution is of this kind.
Recently there have been many numerical studies of the ETH using classical digital computers [21][22][23][24][25][26], which are generally limited to rather small systems, e.g., O (20) Ising spins, because the size of the Hilbert space basis generally grows exponentially with the number of degrees of freedom.Even so, these studies provided strong indications that multi-particle states with ergodic dynamics (on the classical level) obey Eq. ( 1) for many operators.The generalization to SU(N) gauge theories [a particular type of quantum field theories (QFTs) that is invariant under local gauge transformation] with their infinite number of degrees of freedom thus is possible if they are discretized, i.e., if one studies lattice gauge theories (LGTs).At the same time, the level spectrum must be dense for RMT and ETH to apply, implying the need for a large number of states and consequently significant computer power.Hence, we have to content ourselves to provide some numerical circumstantial evidence for selected thermalization properties of nonabelian gauge theories.Still, we believe that such limited evidence can be phenomenologically quite valuable.
An example of thermalization in QCD (SU(3) gauge theory with dynamical quarks) is the production of a quark-gluon plasma as investigated at LHC and RHIC with great effort.For spin systems it has been demonstrated that the macroscopic equilibration time can be much smaller than the Thouless time controlling RMT behavior [27].This suggests that, e.g., a system formed by two colliding heavy ions is not characterized by a single thermalization time but possibly by a range of such times, depending on the observable that is being studied.Different equilibration time scales were indeed observed phenomenologically, leading researchers to postulate the existence of two distinct mechanisms, namely fast hydrodynamization and slow thermalization or phase space equilibration.As the relevant operators for hydrodynamics, i.e., local energy and momentum density, are of rather short range compared to the size of the thermalizing fireball, such a distinction is plausible, but an oversimplification in view of the discussion sketched above.An explicit demonstration of the domain size dependence of the thermalization time was presented in the context of a holographic (AdS/CFT) calculation [7,8], which showed for selected operators that the thermalization time scales linearly with the spatial support of the operator.
To study the properties of a QFT we follow [28] and define a continuum SU(N) gauge theory as the limit of SU(N) LGT for vanishing lattice spacing and an infinite number of lattice points (the continuum and infinite volume limits).Current computational limitations force us to attempt to glean some properties of the full QFT from the properties of very small lattices.In spite of their small size, these lattices still involve Hilbert spaces with dimensions 10 5 − 10 6 .
The long success story of LGT has established beyond any reasonable doubt that the Wilson formulation of LGT gives correct results for a vast range of observables, reproduces the renormalization group running of the continuum theory in the ultraviolet (UV), and has a convergent infrared (IR) behavior for large volumes.The Hamiltonian formulation of Kogut and Susskind (KS) [29] is equally accepted as a valid formulation and forms the basis for the development of quantum computing for nonabelian gauge theories [30,31].Thus, it is sufficient to establish ETH behavior for these well-studied discretized formulations of SU(N) LGT, which can be done numerically, although the continuum limit is probably out of reach for currently available computing resources.In the past, this insight made it possible to establish RMT behavior for Euclidean gauge theories [32][33][34].These studies demonstrated that RMT behavior could already be observed for quite small lattices, which makes us optimistic that a demonstration of ETH behavior for LGT is achievable.
This brings us to the question whether quantum or classical computers are optimal for such numerical studies.On one hand, the Wilson (Lagrangian) formulation can be implemented efficiently on classical computers using imaginary time propagation.Algorithms have been developed which allow for a limited description of real time evolution (see e.g., [35]).However, an attempt to establish the validity of Eq. ( 1) within the Wilson framework would have to go far beyond the demonstrated realm of applicability of such methods.An alternative approach to study quantum chaotic behavior in matrix models and lattice gauge theory is the use of Gaussian density matrices with adjustable parameters [36,37].
On the other hand, very small systems could be insufficient to obtain the results we aim at and quantum computers or quantum simulators might be needed to fully profit from a Hamiltonian formulation.Our present calculations are performed on classical digital computers and correspondingly small systems.Our decision was primarily motivated by previous experience in establishing RMT behavior for gauge theories, which was also possible using very small lattices.Furthermore, classical ergodicity of SU(2) gauge theory could be demonstrated with clearly apparent finite-size scaling on rather small lattices of size L 3 for L = 2, 4, 6 [12].
Presently we also limit ourselves to SU(2) gauge theory for which an efficient formulation in terms of angular momentum algebra was developed in [38][39][40] and extended to SU(3) in [38,41].In addition, the studies reported here will be limited to lattices extending in one and two dimensions.The one-dimensional plaquette chain has been studied as a quantum system for very short chains [40,42].The honeycomb lattice was previously considered by [43], and point-splitting methods for square lattices were discussed in [44] as convenient approaches to implement the KS Hamiltonian in higher dimensions.While the present study is still an exploratory one we hope that improved algorithms and larger computer resources will allow us to better control the continuum and infinite volume limits in the future.
For SU(2) there are three lines of investigation which we pursue in our present study within the constraints of the available computer resources: the dimensionality of the lattice, the number of plaquettes, and the maximum angular momentum representation for electric gauge fields j max .To guarantee the exact reproduction of the KS Hamiltonian convergence in j max must be demonstrated.Recently one of us (X.Y.) proposed a mapping of the 2+1 dimensional SU(2) gauge theory Hamiltonian onto an Ising system valid when the Hilbert space for each gauge link is constrained to the electric field representations j = 0, 1  2 [45].This mapping allows one to find a complete set of energy eigenstates by direct diagonalization for chains of N ≲ 20 plaquettes.While the simulated theory is thus not strictly a discretized version of SU(2) gauge theory, but rather a truncated model of it, the differences could be small, especially at strong coupling.This would be indicative of some universality of ETH behavior.
Overall, as we will show, artifacts caused by the smallness of the studied systems are substantial, which limits the parameter range in which we can trust our results.Within these parameter ranges, however, our results clearly provide evidence for ETH behavior.However, we caution the reader at the outset that the question for which operators Eq. (1) in a QFT should be valid is a highly non-trivial one.Physical operators typically renormalize multiplicatively such that we do not expect deviations from the generic ETH distribution properties.For unphysical operators, however, complications can occur such as operator mixing.The operators we study here are expected to have physical meaning in the continuum limit, which gives us hope that the ETH behavior observed on our small lattices survives in that limit.
We present three groups of results.First, in Section II, we focus on the system of a SU(2) plaquette chain with j max = 1 2 , which can be mapped onto an Ising chain.We study the statistical properties of matrix elements between different nearby energy eigenstates and show that their magnitudes form a Gaussian distribution.We also calculate the spectral function f A (E, ω) and study its dependence on energy E. Second, in Section III, we use a similar mapping for a two-dimensional honeycomb lattice to calculate the complete set of eigenstates (within the truncated Hilbert space) for lattices of size 5 × 4 (for periodic boundary conditions) and triangular lattices of 15 hexagonal plaquettes (for confining boundary conditions).For these we perform similar studies of their ETH properties.Third, in Section IV, we explore the convergence of the energy spectrum with respect to the cutoff j max in the electric field representation on small chains of three plaquettes.Also in this case, we search for the presence of ETH behavior in the energy range for which convergence of the basis is reached and investigate the corresponding spectral function f A (E, ω) in more detail.We finally summarize our results and discuss promising avenues of future research in Section V.

II. LINEAR PLAQUETTE CHAIN
The discretized KS Hamiltonian of the 2+1 dimensional SU(2) gauge theory can be written as [29] where g is the coupling constant, a in the denominator the lattice spacing, E a i the electric field operator along the direction i = x or ŷ with the SU(2) index a (both of which are implicitly summed) and □ ≡ Tr[U † (n, ŷ)U † (n + ŷ, x)U (n + x, ŷ)U (n, x)] the plaquette operator at n = (n x a, n y a) that is the trace of the product of four link variables (Wilson lines) along a square lattice.(In Section III, the honeycomb plaquette operator is defined as the trace of the product of six link variables and the prefactor of the magnetic term in the Hamiltonian will differ.)The coupling constant has the mass dimension of [g] = 0.5; E 2 and □ are dimensionless operators.Throughout the article, we express all dimensionful quantities in units of the lattice spacing a and imply the appropriate factor of powers of a without denoting it explicitly.For example, when we state the coupling constant as g 2 = 0.9 and the energy as E = 2, we imply g 2 = 0.9a −1 (g 2 has dimension of energy in the 2+1 dimensional gauge theory) and E = 2a −1 .
The discretized KS Hamiltonian can be represented in the electric energy basis, which labels the state on each link by the angular momentum quantum numbers |jm L m R ⟩.They specify how the link state transforms under SU(2) gauge transformations generated by the electric fields at both the tail (L) and head (R).When the electric field and the link variable operators act on the corresponding link state, one has [38,39] i,a where ⟨j ′ m ′ |j m; J M ⟩ denotes Clebsch-Gordan coefficients and U n L n R denotes one entry of the SU(2) matrix.Physical states satisfy Gauss's law which in this lattice setup means all link states joining the same vertex transform as a SU(2) singlet.One can write out explicitly the matrix elements of the plaquette operator □ between physical states by using Eq. ( 3) and Gauss's law, which has been done for a plaquette chain in [40,42,46] and a honeycomb lattice in [47].This allows one to explicitly write out the Hamiltonian matrix in the physical Hilbert space and exactly diagonalize it numerically.We begin by further analyzing the statistical properties of the linear color-singlet plaquette chain with periodic boundary conditions on a truncated Hilbert space with j max = 1 2 in the electric field representation.j ≤ j max denotes the SU(2) representation of the electric field operator on the lattice links.It was shown in [45] that the KS Hamiltonian for this system can be mapped onto the Hamiltonian for an Ising spin chain with nearest neighbor interaction ) where N denotes the number of plaquettes in the chain, and the coupling constants are related to the gauge coupling g and lattice spacing a as follows: ( The periodic boundary condition is imposed by setting σ α N = σ α 0 .In the following we present results for g 2 = 1.2 and N = 19, except where stated otherwise.It is expected that in the large volume limit all nonzero values of g 2 are equivalent for the purpose of demonstrating the ETH, because the Hamiltonian is non-integrable for all nonvanishing g 2 .However, for the small lattice sizes that are numerically accessible on classical computers, some choices of g 2 serve better to demonstrate the ETH behavior.Other choices will eventually exhibit the ETH behavior when the lattice becomes sufficiently large.Numerical evidence for this statement can be found in Appendix A, where we show the energy spectra in the case with g 2 = 1 for different chain length N . The space of energy eigenstates for this system can be decomposed into 19 sectors with good linear momentum pa = 2πk/N with integer k ∈ {0, 1, • • • , N − 1}.k = i and k = N − i correspond to two sectors that are related by the reflection.Each sector contains (2 N − 2)/N = 27, 594 states with the exception of the k = 0 sector, which contains two additional states, the states with all the spins pointing downward or upward.
The energy level spacing statistics and the statistical properties of the matrix elements of the local 1-plaquette and 2-plaquette operators were studied in [45] with results that supported the assumption that the ETH is manifested in this system.Because these operators are not translation-invariant, they have nonvanishing matrix elements between states with different k.Here we consider matrix elements of the total electric energy operator given by the first two terms in the Hamiltonian Eq. ( 4) which is translation-invariant and only couples states in the same momentum sector.In the following we focus on the k = 1 sector; nearly identical results are obtained for other momentum sectors.(For the k = 0 sector, the demonstration of Gaussian Orthogonal Ensemble (GOE) behavior in the level spacing statistics requires separation of the sector further into two subsectors with definite parity [45].This will be done in Section IV where we study the convergence of the spectrum under changes of the j max cutoff for a short plaquette chain.The distribution of level spacings was already studied in [45] and shown to be well described by the GOE distribution.As an additional sensitive statistical measure of the energy gap distribution, we have calculated the normalized energy gap ratio [48] where δ α = E α+1 −E α is the energy gap between adjacent levels.The advantage of this measure is that it is not sensitive to the change of the level density with energy.
In order to avoid distortions from the tails of the spectrum, we only include states in the range ( , where E min and E max are the lowest and highest energy eigenvalues, respectively, and E av is the mean energy of the spectrum [24].In order to investigate the statistics of off-diagonal matrix elements and calculate the spectral function f el (E, ω) for the electric energy operator, we select energy windows of width ∆E = 1 shown in Fig. 1.In the selection of the energy eigenstates for the study of matrix element statistics one has two options.One choice is to select energy pairs with a constraint on the mean, where Ē is a chosen constant, the other is to require that both energies fall into a common window, |E α − Ē|, |E β − Ē| < ∆E/2.The first option corresponds to selecting a diagonal window in energy pair space, the second chooses a square window.
In the limit |ω| = |E α − E β | ≪ ∆E/2 the two options effectively coincide.Here we choose the second option, because it is more closely related to the long-time behavior of a state defined by a wave packet with a narrow energy spread.
We begin with the window 0 < E < 1 (shown in red in Fig. 1 and corresponding to Ē = 0.5, henceforth called Window 1), where the level density has its maximum.The window contains 2173 eigenstates.Figure 3 Similarly good fits are obtained for other ω windows of the same width in the interval |ω| < 0.1.The width of the Gaussian is related to the spectral function f el (E, ω).Using the ETH relation Eq. ( 1), Tr[R 2 ] ≡ R αβ R βα = 1, and e S(E) = ρ(E) one finds representing a diffusive transport peak superimposed on a flat pedestal.The fits are shown as solid, dashed, and dotted lines, respectively.We shall discuss the shape of the spectral function in more detail and over a wider ω range in Section IV.
The final statistical test of off-diagonal matrix elements we perform is to check whether they form a GOE.It is also worth noting that, although the matrix elements are complex numbers for energy eigenstates with definite momentum, the associated Gaussian ensemble for comparison is the GOE (and not the Gaussian Unitary Ensemble, GUE), because the Hamiltonian is time reversal invariant.We calculate the second and fourth moments of the band matrix obtained by keeping only matrix elements between eigenstates that differ in energy by less than ∆E(T ) = 2π/T [27].The idea is that contributions to matrix elements at time T between wave packets from eigenstate pairs α, β with |E α − E β | > ∆E(T ) cancel out because of their random energy phases, but the randomness of matrix elements between states differing by less than ∆E(T ) may be required for some thermal properties to manifest themselves.As time progresses, only the contributions from an increasingly narrower band of energy eigenstate pairs are relevant.
We define the projected submatrix elements of an operator O as where the diagonal part is always included.Our definition differs slightly from the one adopted in [27] in that we consider all pairs of states within each of the fixed (T -independent) energy windows shown in Fig. 1.As T increases, the matrix O T contains a growing number of vanishing elements.The GOE measure Λ T is defined as where d is the number of eigenstates in the chosen energy window and is the traceless part of the matrix.In our case O T c is complex but Hermitian, which ensures that Tr[(O T c ) n ] is real.For matrices approaching a perfect GOE at late times, one expects Λ T → 1 2 .The Λ T measure stands out by being a sensitive test not only of the Gaussian distribution of the absolute values of the matrix elements, but also for correlations among their signs, which encode the orthogonality of the GOE matrix ensemble.In practice, as a smaller and smaller number of elements of the projected matrix O T is nonvanishing for T → ∞, the statistical significance of the result for Λ T deteriorates.This can be seen if one arbitrarily replaces many of the elements of a large matrix from the GOE by zero, converting it into a band matrix with sparse off-diagonal elements.Our results for the three energy windows for N = 17 (dashed lines) and N = 19 (solid lines) are shown in Fig. 5.As can be seen, the value of Λ T increases with the energy of the window.The values Λ T ≈ 0.45 − 0.47 closest to the GOE limit 0.5 are reached for Window 3, for which the deviation from the GOE value 0.5 is compatible with those expected for a sparse off-diagonal GOE matrix.However, the deviations for the the other two energy windows are too large to be explained in this way.This means that the off-diagonal matrix elements of the operator H el in these energy windows do not form GOEs although the distribution of the absolute values is Gaussian.We attribute this to the fact that the electric energy operator in the j max = 1 2 truncated basis is constrained to have values that are multiples of j max (j max + 1) = 3  4 .We shall see in Section IV that the off-diagonal matrix elements of H el are closer to a GOE when j max is increased, and we expect this property to be realized in the con-tinuum limit where the spectrum of the electric energy operator becomes continuous.It will be interesting to explore whether the matrix elements of other operators exhibit a faster approach to the GOE.

III. HONEYCOMB LATTICE
In this section we will generalize the previous onedimensional (1D) plaquette chain study to two dimensions (2D).The new dissipative mode that contributes to thermalization in 2D, which does not exist in 1D, is the shear mode.Scale invariance, which is exact in gauge theories at the classical level and only broken at the quantum level by the trace anomaly, constrains dissipation of certain modes.This manifests itself by the strong suppression of bulk viscosity in gauge theories at high temperature.On the other hand, dissipation in the shear channel is not suppressed by symmetry and remains finite at strong coupling.Therefore, the new physical insight one can gain from the ETH study in 2D is the thermalization dynamics related to shear viscosity.
We will present the results obtained for the honeycomb lattice with j max = 1 2 .At each vertex, physical states can be fully specified by the j values on the links connecting to the vertex if the vertex has fewer than four links.If the vertex has four links connected or more, as on a square lattice in 2D (spatial dimensions) or a cubic lattice in 3D, physical states depend on the order of how the color representations on different links are coupled and thus cannot be fully specified by these j values alone.This motivates the use of the honeycomb lattice where each vertex is touched by at most three links.Two examples of honeycomb lattice setups are shown in Fig. 6.
The Hamiltonian of the 2D SU(2) gauge theory on a honeycomb lattice with j max = 1 2 has been mapped onto a 2D spin model in [47] for both periodic and closed boundary conditions.We will first consider the periodic boundary condition before discussing the closed boundary condition.

A. Periodic Boundary Condition
The spin Hamiltonian for the original SU(2) lattice gauge theory with j max = 1 2 on a parallelogram under the periodic boundary condition can be written as [47] H = (i,j) where The (i, j) symbol labels the plaquette position at j along the x-direction defined by (1, 0) in the Cartesian 2D plane and i along the y-direction defined by 2 ), as shown in Fig. 6a.We assume there are N x plaquettes along the x-direction and N y along the y-direction.The first term of Eq. ( 14) is the electric part of the Hamiltonian while the second term is the magnetic part.The magnetic part involves an exponential of Pauli matrices but one can rewrite it purely in terms of tensor products of Pauli matrices [47].In our current studies using classical computers, both work equally well.For a quantum simulation, the latter format is more convenient to use since it can be easily mapped onto qubits.
Under the periodic boundary condition, the Hamiltonian in Eq. ( 14) is translationally invariant along the x and y directions by one lattice unit defined above: [H, Tx ] = [H, Ty ] = 0.As a result, the Hamiltonian and the translational operators can be simultaneously diagonalized.Each eigenstate belongs to a particular momentum sector given by 2πk x /N x , 2πk y /N y where k . The Hamiltonian in each momentum sector can be constructed as discussed in [47] and diagonalized exactly for small lattice sizes.
We first study the statistics of the eigenspectrum.We consider the k x = 1, k y = 1 sector in the N x = 5, N y = 4 system, which contains 26163 states.(We choose N x ̸ = N y to avoid discrete hexagonal symmetry.)We plot the density of eigenstates in Fig. 7 for g 2 = 0.75, where we also plot the distributions of the energy gaps δ between nearest eigenstates and the rescaled gaps s ≡ δ/δ (δ is the average of the gap), and the distribution of the restricted gap ratios defined in Eq. ( 7), for the eigenstates in the middle of the energy spectrum (we removed the first 5000 and last 5000 eigenstates when generating these three plots).The red curve in the gap distribution Fig. 7b is   which is a feature of a quantum system whose classical counterpart is chaotic.
We emphasize that level repulsion is very sensitive to any discrete symmetry of the system.If we plot the gap distribution for momentum sectors involving a zero momentum (i.e., k x = 0 and/or k y = 0), level repulsion cannot be seen and the distribution does not exhibit a Wigner-Dyson shape.This is due to the remaining parity symmetry in the zero momentum sectors.If we further separate each zero momentum sector into a parity-even and parity-odd sector and plot the gap distribution in each parity sector, we will see level repulsion clearly again and the distribution can be well fitted by a Wigner-Dyson shape, as shown for the chain case in [45].In Fig. 7c, the red curve is a prediction from the GOE of 2 × 2 random matrices which states In Fig. 7d, the red curve is a GOE prediction as written in Eq. ( 8).The expectation value of the restricted gap ratio can be calculated as ⟨r⟩ ≈ 0.5284, which is very close to the GOE prediction ⟨r⟩ GOE ≈ 0.5307.As can be seen, the statistics of eigenenergies in the middle of the spectrum can be well described by the GOE, up to statistical uncertainties due to the finite size.Before moving on to operator matrix elements, we want to comment on the effect of g 2 values on the energy spectrum.Here we focus on the case with j max = 1 2 .In Section IV, we will discuss how to choose j max with varying values of g 2 in order to obtain physical results.In the large volume limit, it is expected that any value of g 2 will lead to qualitatively similar level statistics that are well described by the GOE.However, with a finite lattice size that is numerically accessible by current classical computers, some choices of g 2 values will give better results than others.For example, the reasonably good results shown above are obtained from the choice g 2 = 0.75.If we choose g 2 = 1, we observe a spiky structure in the level density, which distorts the distributions of gaps and gap ratios away from the GOE predictions.We believe these are finite volume effects and expect them to go away as the lattice size increases.Numerical evidence for this expectation is shown in Appendix A for the plaquette chain (we choose the chain to demonstrate this since we can go to a relatively large lattice size in the quasi-1D case such that momentum level separation is smaller and so is relevant kinetic energy).
Next we study the matrix elements of the electric energy in the energy eigenbasis.The electric energy is given by which is translationally invariant along the x and y directions.As a result, the off-diagonal matrix elements of H el vanish between states in different momentum sectors.In the following, we will only compute the matrix elements of H el in one momentum sector.We consider the k x = k y = 1 momentum sector on the N x = 5, N y = 4 lattice and plot the magnitudes of the off-diagonal matrix elements, i.e., |⟨α|H el |β⟩| in Fig. 8 for states satisfying 2 < E α + E β < 4 for two choices of couplings: g 2 = 1 and g 2 = 0.75.Since the distributions are symmetric around ω ≡ E α −E β = 0, we only plot the distributions for ω > 0. If the ETH applies, the distribution of the off-diagonal matrix elements in a small ω window should be Gaussian when the number of elements is large.We test this in the small ω region by comparing the width obtained in two ways.In the first way, we fit the random variable distribution in an ω window of size 0.02 to a Gaussian to obtain the width, shown in magenta.Alternatively, we calculate the second moment of the matrix constrained to pairs of states with ω min < ω < ω max , shown in black (ω min and ω max are the lower and upper bounds of an ω window).If the off-diagonal matrix elements magnitudes satisfy a Gaussian distribution, the two methods will give the same result.In the case of g 2 = 1 (Fig. 8c) we see a noticeable disagreement between the two methods for small ω, which indicates the magnitudes of the off-diagonal matrix elements do not exactly follow Gaussian distributions.On the other hand, we see a very good agreement between the methods in the case of g 2 = 0.75 (Fig. 8d), which shows the magnitudes of the off-diagonal matrix elements in the small ω region follow Gaussian distributions.The contrast between the two cases of different couplings is consistent with the above studies of the level statistics, where we see the case with g 2 = 0.75 can be well described by the GOE.
We also look at the off-diagonal matrix elements in a wider ω region.As can be seen, the distributions shown in Figs.8a and 8b exhibit bumpy plateaus at small and intermediate ω and decay exponentially at large ω.The case with g 2 = 1 shown in Fig. 8a contains more bumpy structures, in particular at small ω as shown in Fig. 8c, where we fit the ω dependence of the second moment extracted width by the function a/(ω 2 + b 2 ) + c.The peak near ω ≈ 0 is likely related to diffusive transport processes.However, as we decrease the coupling to g 2 = 0.75, we find the peak disappears, as shown in Fig. 8d.The reason for the disappearance of the transport peak at weaker coupling is unclear and may be a small lattice artifact.
We next study the GOE measure Λ T defined in Eq. ( 13), which tracks the approach of the off-diagonal part of the matrix to GOE behavior.We choose four energy windows in the case of g 2 = 0.75 for our analysis: −2 < E < −1, −1 < E < 0, 0 < E < 1 and 1 < E < 2. The results are shown in Fig. 9.If the off-diagonal matrix elements are described by the GOE, the value of Λ T will reach 0.5.We see some deviations from the GOE  prediction, which indicates the off-diagonal matrix elements are still correlated in their signs.(We have shown in Fig. 8 that the magnitudes of the off-diagonal matrix elements satisfy Gaussian distributions in the small ω region.So the deviations seen here reflect the sign correlations and thus non-orthogonality rather than the non-Gaussianity.)It is worth further investigating the Λ T observable on bigger lattices in future work.
Finally, we study matrix elements of operators that break the translational invariance.Following [45], we study Wilson loops corresponding to 1-plaquette O 1 and 2-plaquette O 2 operators.The Pauli matrix representation of these Wilson loop operators and their matrix elements in the momentum basis have been worked out for the honeycomb lattice in [47].We use all the momentum sectors up to (including) k x = ⌊N x /2⌋ and k y = ⌊N y /2⌋ in the g 2 = 1 case.We first investigate lattice size dependence of the diagonal part of their matrix elements by considering lattices of size 3 × 3, 4 × 3, 4 × 4, and 5 × 4. We estimate the deviation of the diagonal part from the microcanonical ensemble average by using a proxy for the microcanonical ensemble made up of 10 eigenstates below and 10 above the eigenstate under consideration [45].The deviation is estimated as We plot the magnitude of ∆ i averaged over all eigenstates except for the 20 lowest and 20 highest states in Fig. 10, as a function of lattice size.We find that |∆ i | decreases approximately exponentially with lattice size, confirming the exponentially decaying factor e −S/2 in the fluctuating part of Eq. ( 1) for most eigenstates, since the entropy is an extensive quantity, S ∝ N x N y .In other words, we demonstrated the diagonal matrix elements of the two Wilson loop operators are exponentially close to the microcanonical ensemble average value, as the system size increases.
We also plot the off-diagonal matrix elements of the two Wilson loop operators in Fig. 11.In this plot, we consider eigenstates within a thin energy shell 1.99 < E α + E β < 2.01 in all the momentum sectors up to (including) k x = 2 and k y = 2 on the N x = N y = 4 lattice with g 2 = 1.The upper envelops of these distributions also exhibit bumpy plateaus at small and intermediate ω and show hints of an exponential fall-off at large ω, although the statistics becomes marginal above ω = 12.5.

B. Closed Boundary Condition
In this subsection, we consider a triangular honeycomb lattice with a closed boundary condition, as shown in where √ 3 8 g 2 , and h x and c ij are the same as in the case of periodic boundary conditions.In the following, we will consider g 2 = 1.As explained earlier, any value of g 2 is expected to give qualitatively similar results for asymptotically large lattice sizes.
As for the other systems investigated here, we exactly diagonalize the Hamiltonian and calculate matrix elements of certain operators in the energy eigenbasis.Here we will omit the studies of level statistics since the system has a discrete symmetry given by the dihedral group D 3 and focus on studying operators corresponding to the 1-plaquette and 2-plaquette Wilson loops, which were constructed in [47].Unlike the periodic case, here the results of operator matrix elements will depend the location of the operator.We consider two scenarios: Wilson loop operators defined at the center of the lattice, and those defined at a corner of the triangular lattice shown in Fig. 6b.
As in the periodic case we first plot the deviation of the diagonal matrix element from the microcanonical ensemble average in Fig. 12.For the operators defined at the lattice center, we clearly see an exponential decay of the averaged deviation magnitude as the system size increases (Fig. 12a), which indicates the diagonal matrix elements are rapidly approaching the microcanonical expectation values.For the operators defined at the corner, we see a cusp in the 1-plaquette operator case (upper panel of Fig. 12b), which could be a boundary effect.However, for the 2-plaquette operator the exponential decrease is visible also for the corner location (lower panel of Fig. 12b).
Then we study the off-diagonal matrix elements of the Wilson loop operators on the N = 5 lattice.When calculating the matrix elements, we use only eigenstates within a thin energy shell 43.999 < E α +E β < 44.001.(One cannot simplify the Hamiltonian in the closed boundary case, i.e., Eq. ( 20) by simply removing the constant terms, which shifts the zero energy reference such that the middle of the spectrum is around zero energy, as done in the case with a periodic boundary condition.This is because after removing the constant terms, one has to consider how the plaquettes (spins) just outside the boundary enters Eq. (20), which is different for the three edges of the triangular lattice.So we choose not to shift the zero energy reference here for convenience in the numerical construction of the Hamiltonian.This is why the values of the energy windows in the two cases with different boundary conditions look so different.) The results for the absolute values of the matrix elements are shown in Fig. 13 as a function of ω = E α − E β .Again, we notice a bumpy plateau in the small ω region (ω ≲ 10) and an exponential fall-off at large ω for both operators, no matter whether they are defined at the center or corner.In the case of the 1-plaquette operator, some off-diagonal matrix elements have much smaller values (≲ 10 −13 ) than others, which are essentially zero within the numerical precision.This is a result of discrete symmetries in the system.The 1-plaquette operator defined at the center is symmetric under the dihedral group D 3 on the N = 5 triangular lattice.As a result, its matrix element between two eigenstates with different symmetry under D 3 transformations vanishes by virtues of the Wigner-Eckart theorem.Similarly, the 1-plaquette operator defined at a corner of the triangular lattice respects the reflection around the perpendicular bisector line going through this corner.Thus its matrix element between two eigenstates that have different reflection symmetry vanishes.We do not see vanishing off-diagonal matrix elements in the case of 2-plaquette operators, since they FIG.14. Gaussian fits for distributions of ⟨n|O2|m⟩ (the operator is defined at the center) in two ω slots: 0 < ω < 4 (left) and 4 < ω < 8 (right).The fitted Gaussian widths are roughly 0.00200 and 0.00174 for the left and right panel, respectively.
break the D 3 symmetry of the lattice.
Finally, we study whether the off-diagonal matrix elements of the Wilson loop operators can be well described by Gaussian distributions.To this end, we consider the 2-plaquette operator defined at the center of the N = 5 lattice studied above and use the eigenstates contained in the same energy shell.The distributions of the off-diagonal matrix elements ⟨n|O 2 |m⟩ in different ω regions are plotted in Fig. 14 together with Gaussian fits.The two different ω windows are 0 < ω < 4 and 4 < ω < 8.We see in the smaller ω region, the distribution of ⟨n|O 2 |m⟩ can be better described by a Gaussian, which is consistent with the results obtained for the plaquette chain and the honeycomb lattice with periodic boundary conditions.

IV. CUTOFF CONVERGENCE
As explained in Section I we are interested in the double limit of large number of plaquettes N and large enough j max to reach convergence.This is not feasible for us with our present computer resources.Therefore we investigate separately the cases N large, j max = 1 2 (Section II and III) and the case N = 3, j max ≤ 7 2 (this section).More precisely we investigate the validity of the ETH relation Eq. ( 1) for a linear chain of three plaquettes with periodic boundary conditions in an energy region where we observe good convergence with the cutoff j max , i.e., we obtain results for the KS Hamiltonian in a converged Hilbert space region.
We start by examining this convergence behavior.A few typical examples of eigenenergies for g 2 = 0.8 are shown in Fig. 15.Eigenstates with different momentum k (which is discretized) and/or parity do not mix.Therefore, to study the characteristic RMT behavior in the eigenenergy spectrum one has to choose a specific sector.We concentrate on the sector with momentum k = 0 and positive parity.An additional symmetry, namely topbottom symmetry, arises for j max > 1 2 .In this case we chose the sector corresponding to the eigenvalue +1.Ob-FIG.15.Some energy eigenvalues for the k = 0, positive parity and positive top-down symmetry sector plotted against the cutoff for g 2 = 0.8 with a corresponding exponential fit of the form of Eq. ( 21).
viously, very good convergence can be achieved for states that are not too highly excited, i.e., E − E 1 ≲ 20, where E 1 denotes the ground state energy of the system.This convergence can be tracked quantitatively by fitting the eigenvalues as function of j by simple exponential curves of the form and extrapolating to infinite j (see dashed lines in Fig. 15).Here λ α is the α-th energy eigenvalue in the limit j → ∞ and A α is a fit parameter which does not depend on j max but on α.We observe a power law behavior of the parameter A α with respect to α which we use to further improve the extrapolation: We fit A α by a threeparameter function A α = a • α b + c.This finally leaves a function E α (j) with only one parameter, namely λ α , allowing us to define a precise criterion for a converged eigenvalue, where in our case j max is 7  2 .Using this criterion we find convergence for eigenvalues up to E max = 24.3.It is clear from Fig. 15 that j max = 1 2 is typically far from convergence at this value of the coupling constant (g 2 = 0.8), which is not surprising since even the ground state energy becomes exact for j max = 1 2 only in the strong coupling limit.
In the selected converged energy regions we proceed with the analysis of ETH properties.In line with our discussion in Section I we choose a "physical" observable, namely the electric energy H el , in a specific symmetry sector to check the validity of Eq. ( 1).Furthermore, we study the structure of the spectral response function f (E, ω) in more detail.Since the possibilities of creating gauge invariant states on our lattice increase with increasing j max , new energy eigenvalues emerge such that the available statistics improve rapidly with increasing j max .
Again, we proceed as in Sections II and III and follow the template of similar investigations of possible ETH behavior (see e.g.[24]): We first check for RMT behavior in the eigenenergy spectrum, before we analyze the diagonal, and then the off-diagonal matrix elements of Eq. ( 1).However, the KS Hamiltonian has a special property one needs to be aware of.As one can see in Eq. ( 2) the KS Hamiltonian has an electric and a magnetic part.A sufficient contribution of the magnetic part to the overall Hamiltonian is essential for observing chaotic behavior.Thus, for fixed lattice spacing a and size N the degree of ergodicity depends on g.For sufficiently large g the theory becomes non-ergodic and the resulting energy level statistics approach a Poisson distribution while for moderate g its RMT behavior should be clearly visible.For example the mean gap ratio should interpolate smoothly between the Poisson and RMT (GOE) values, which is exactly what we observe in Fig. 16.This plot also shows in which range g 2 has to be chosen to observe ergodic properties for N = 3, namely g 2 ≤ 1.2.We choose typically g 2 = 0.8 or g 2 = 1.Smaller values of g 2 are not practical, because the slow convergence with respect to j max would yield an insufficient converged region size.
Overall, requiring convergence with respect to j max , a value of g 2 for which the system behaves ergodically, and the magnitude of discretization and finite size artifacts strongly limit the usable statistics and parameter range.This is illustrated in Fig. 17 for the histograms of the complete spectra for electric field cutoffs j max ∈ { 5  2 , 3, 7 2 }.While there is still a large difference in the level density in a wide energy region between j max = 3 and j max = 7  2 , our convergence studies described above allow us to assess the convergence beyond j max = 7 2 , i.e., confirm that FIG. 17. Histograms of the whole energy spectra for jmax ∈ { 5 2 , 3, 7 2 } and g 2 = 0.8.The blue rectangle depicts the converged region, leaving out states affected by finite size effects.
ρ 7/2 (E) ≈ ρ ∞ (E) for E ≲ E max , where the index of ρ j (E) denotes the j-cutoff.In order to suppress finite size effects, we discard 14% of the states at both edges of the full spectrum.As a result, the low-energy part of the converged spectrum is not taken into account.The truncated converged energy window, used in our analysis, is thus determined by the upper bound derived from our convergence criterion (22) and the lower bound from the finite size truncation.The converged window is highlighted by a blue shaded rectangle in Fig. 17.The remaining number of states and matrix elements of H el is large enough to obtain rather precise results.
Another standard test for GOE behavior in the eigenenergy spectrum is the nearest-neighbor level spacing statistics.Figure 18 shows a histogram of the normalized level spacing distribution P (s) in the converged energy window for g 2 = 0.8 in comparison with the Wigner-Dyson distribution (16).The agreement is good, but not perfect, suggesting that there are still deviations from exact RMT behavior in this energy window.
Next, we analyze the diagonal part of Eq. ( 1) for the matrix elements of the electrical energy operator in the specified symmetry sector.Equation (1) predicts the presence of exponentially suppressed rapidly varying fluctuations superimposed on a smooth function of E. This ensures the existence of a well defined microcanonical ensemble average that is obtained by averaging ⟨E α |H el |E α ⟩ over a small energy window.Figure 19 shows the distribution of the expectation value of the electric energy as a function of the energy for all energy eigenstates of the system in the selected symmetry sector.The values for states inside the chosen energy window are highlighted in red.The spread of the distribution around a nearly linear dependence on energy is seen to be small, i.e., the electric energy operator has a well defined microcanonical average ⟨H el ⟩ mc (E) over a wide range of the spectrum.The width of the distribution is expected to shrink exponentially with the lattice size as observed for the plaquette operators in the truncated (j max = As an aside let us note that Fig. 19 shows that the electrical energy contributes the dominant part to the total energy, with E el ranging between 0 and 56.7, whereas the magnetic energy is limited to the range −15 < E mag < 15 in lattice units.This shows that the system is still far from the thermodynamic limit of the continuum gauge theory.
We now proceed with the analysis of the off-diagonal matrix elements of the electric energy, i.e., ⟨E α |H el |E β ⟩, in the converged energy window.Figure 20 shows these matrix elements in the energy window with Ē = 21.5 and ∆E = 1 for g 2 = 0.  an exponential decay at large ω, a bumpy region at intermediate ω reflecting quasiparticle contributions, and a diffusive plateau at the smallest values of ω.For nonabelian gauge theories an effective description in terms of quasiparticles is typically based on analogs of glueballs; at high excitation energy such quasiparticles may or may not exist depending on the strength of the coupling.Such quasiparticles show up as broad peaks in the spectral function.Typically these structures are quite broad but details are specific to the theory under investigation.These properties are visible in Figs.20a and  20b.
For many systems one also observes a transport peak and a diffusive plateau at very low ω.Whether a prominent transport peak exists for nonabelian gauge theories and, if so, for which operators and in which parameter range, is still under investigation [50][51][52].Interestingly, our results shown in Figs.20 and 21 support the existence of such a diffusive transport peak, at least for the pure SU(2) gauge theory.Therefore we now focus our analysis specifically on the small ω regime, i.e., the transport peak region.We start by discussing the g dependence of the transport peak, shown in Fig. 21.In order to be able to compare the behavior for different coupling constants we define the equivalent energy regions with respect to each ground state.In the following we restrict ourselves to the energy window with Ē − E 1 (g 2 ) = 31 and ∆E = 1.We know from our discussion of Fig. 16 that the chaotic nature of our system is lost if the coupling constant exceeds g 2 ≈ 1.2.As the diffusive transport peak is part of this dynamics, it is natural to expect it to become less well defined when g 2 approaches or exceeds 1.2.Indeed, Fig. 21 suggests that the diffusive plateau at small |ω| disappears right around this value of the coupling constant although there is still a pronounced peak in the spectral function.As explained in Section III the spectral function can be obtained in different ways.On one hand it can be defined via the second moment of matrix elements as already shown in Figs.20 and 21.On the other hand, the distribution of matrix elements magnitudes in small ω windows appears Gaussian, suggested by the central limit theorem, which is clearly supported by Fig. 22.Similar plots are obtained in the whole ω spectrum.Thus we are in the position to define the spectral function from Eq. ( 10), using the Gaussian width of the matrix element distributions.The comparison of these two methods for g 2 = 0.8 is depicted in Fig. 23.Both the second moment and the Gaussian fit methods show very good agreement, even for the diffusive plateau |ω| < 0.02, supporting that the offdiagonal matrix element magnitudes follow a Gaussian distribution.The plot also supports a functional form of a diffusive transport peak, i.e., Eq. ( 11).
Finally we want to discuss certain limitations of the GOE indicator used by considering the Λ T measure defined in Eq. ( 13).Unlike the above analyses, the calculation of the Λ T measure does include the original matrix elements without taking the absolute value.In Fig. 24 we show Λ T for three different sizes of the energy windows, namely ∆E = 1, ∆E = 1.3 and ∆E = 1.6, around the same mean energy Ē = 23.5 (all of which lie in the blue band denoting the converged energy region indicated in Fig. 17).All three plots exhibit the expected monotonic growth with increasing T , i.e., decreasing number of allowed nonzero matrix elements, and convergence to some saturation value Λ ∞ ∆E .In the case of ∆E = 1 at T = 14, 000, we find this value to be Λ T 1 ≈ 0.429, where we are left with 905 nonzero matrix elements.For the two larger energy windows and the same T we obtain Λ T 1.3 ≈ 0.482 and Λ T 1.6 ≈ 0.490 with 1182 and 1474 nonzero matrix elements, respectively.We observe a sustained approach to the GOE prediction Λ GOE = 0.5 for increasing energy window width ∆E.We interpret this as a result of the larger number of nonzero matrix elements at each T as the energy window width ∆E increases.This finding highlights our numerical limitations for small energy windows, e.g., the case ∆E = 1.These results also show similar saturation time scales of Λ T for the electric energy H el of order T ∼ 10 4 as found in Section III.

V. SUMMARY AND OUTLOOK
We investigated several truncated versions of latticediscretized 2+1 dimensional SU(2) gauge theory with respect to the properties underpinning the Eigenstate Thermalization Hypothesis.The systems we studied include linear plaquette chains up to length N = 19, honeycomb lattices with up to 20 plaquettes with the electric field Hilbert space truncated at j max = 1 2 , and the N = 3 plaquette chain on the fully converged electric field Hilbert space.Our results are encouraging: All these truncated versions of lattice gauge theory exhibit clear signs of behavior that is consistent with the ETH: The energy level spectrum obeys GOE statistics within the limits of statistical and systematic uncertainties.The fluctuations of diagonal matrix elements of the operators we studied (the total electric energy as well as 1-and 2plaquette Wilson loops) were found to exponentially decay with the lattice area.We studied the off-diagonal matrix elements of several operators between nearby energy eigenstates and found their magnitudes follow Gaussian distributions.Their signs are correlated and only become more random in smaller ω windows.We also calculated the spectral function f (E, ω) for all three truncated LGT models and found that they conform to the expectations for a system that exhibits quantum chaos.
While encouraging, these pioneering studies were limited to small systems due to lack of computer resources and thus rather constitute a proof of principle than definitive evidence that 2+1 dimensional SU(2) lattice gauge theory exhibits ETH behavior in the continuum limit.However, they motivate hope that along these lines the numerical evidence can be significantly strengthened if more computer resources are invested.In addition to substantial hardware resources this will also require the development and implementation of efficient algorithms for the calculation of matrix elements between energy eigenstates in very large Hilbert spaces, such as the kernel polynomial method [53].
The following bullet points illustrate the range of possible future work extending the results reported here: • Explore the convergence of Hilbert space truncations for lattices with a larger number of plaquettes.
• Extend our investigations to more operators A to better understand for which operators the defining ETH relation (1) holds.
• Extend our studies to three spatial dimensions using point-splitting methods as in [44].
• Determine the Thouless energy for 2+1 dimensional SU(2) gauge theory if it exists, i.e., the value of ∆E below which the RMT behavior in the offdiagonal matrix elements is observed, as a function of lattice size.
• Extend our studies to SU(3) gauge theory using the techniques of [40].
Obviously, realizing this agenda will require years of dedicated work.Our present contribution only marks its beginning.
Figure 2 depicts our result for P (r), shown as a histogram with red bars.The analytical distribution Eq.(8), shown as a blue line, is seen to provide an excellent description, and the mean value ⟨r⟩ = 0.5340 agrees well with the GOE prediction.

FIG. 2 .
FIG. 2. Histogram of the distribution of restricted energy gap ratios rα (red horizontal bars), shown together with the GOE limit distribution PGOE(r) (blue line).
shows a histogram of the values |M αβ | ≡ |⟨α|H el |β⟩| in the ω window 0.02 < |ω| < 0.04 together with a fit to a Gaussian distribution of the form

FIG. 3 .
FIG. 3. Distribution of matrix elements |M αβ | of the total electric energy operator H el in the range 0.02 < |ω| < 0.04, shown together with a Gaussian fit.

FIG. 5 .
FIG. 5.The GOE measure Λ T for the projected matrix elements Eq. (12) of the electric energy operator as function of time T for eigenstates in Window 2 (bottom), Window 1 (middle) and Window 3 (top).Solid lines: N = 19; dashed lines: N = 17.

FIG. 6 .
FIG. 6. Examples of the honeycomb lattice considered in this work: (a) a parallelogrammatic shape with Nx = Ny = 3; (b) a triangular shape with N = 3.
a fit from the Wigner surmise P ws (δ) = aδ exp(−bδ 2 ).The fitted parameter values are a ≈ 10164541, b ≈ 5337217.Level repulsion is clearly seen in the gap distribution,

P
FIG. 7. Density of eigenstates, distributions of gaps, rescaled gaps and gap ratios in the momentum kx = ky = 1 sector on the Nx = 5, Ny = 4 lattice for g 2 = 0.75.The red curve in (b) is a Wigner surmise fit and the other two red curves in (c) and (d) are GOE predictions.

FIG. 8 .
FIG. 8. Off-diagonal matrix elements of the electric energy in the eigenbasis in the kx = ky = 1 momentum sector on the Nx = 5, Ny = 4 lattice for g 2 = 1 (left) and g 2 = 0.75 (right).The plots in the top row are in the full ω region while the bottom row is an enlargement of the small ω region.The black curves are the widths calculated by the second moment of the matrix element distribution.The magenta curves are the widths fitted by Gaussian distributions.The yellow curve is a fit of the second moment results in the g 2 = 1 case with a ≈ 3.50 × 10 −5 , b 2 ≈ 1.34 × 10 −3 and c ≈ 5.77 × 10 −3 .

FIG. 10 .
FIG.10.Averaged magnitude of the difference between the diagonal matrix element and the microcanonical ensemble average decays roughly exponentially as a function of system size for both Wilson loop operators.

FIG. 12 .
FIG. 12. Averaged magnitude of the difference between the diagonal matrix element and the microcanonical ensemble average as a function of the lattice size, for Wilson loop operators defined at the center (left) and the left bottom corner (right).

Fig. 6b .
Fig.6b.The total number of plaquettes of a triangular lattice is fully determined by the number of plaquettes on each side N : N tot = N (N + 1)/2.Closed boundary conditions assume that all links outside the boundary have vacuum quantum numbers j = 0.With j max = 1 2 , the 2D SU(2) lattice gauge theory can be mapped onto a 2D spin model[47]

FIG. 13 .
FIG. 13.Magnitudes of off-diagonal matrix elements of the Wilson loop operators defined at the center (left) or left bottom corner (right).The two eigenstates are picked up from a narrow energy window 43.999 < Eα + E β < 44.001, which leads to the discontinuity in ω here.

FIG. 16 .
FIG.16.Average restricted gap ratio ⟨r⟩ plotted against the coupling constant g 2 for N = 3 plaquettes and jmax = 3.5.The value for a GOE distribution is 0.53 and that for a Poisson distribution is 0.39.Obviously ⟨r⟩ extrapolates between the ergodic and non-ergodic regime as a function of g 2 .

1 2 )FIG. 18 .
FIG. 18.The normalized level spacing s = δα/ δ (as defined in Section III A) distribution in the converged energy window compared to the Wigner-Dyson distribution P (s).

FIG. 20 .
FIG. 20.Absolute values of the off-diagonal matrix elements of the electric energy in the energy basis |⟨Eα|H el |E β ⟩| against ω = Eα −E β in the energy window |(Eα +E β )/2− Ē| < ∆E/2 with Ē = 21.5 and ∆E = 1.Part (a) presents a logarithmic plot of all matrix elements (blue) and their Gaussian widths (orange), obtained by the second moment method, whereas (b) shows a linear plot of the Gaussian widths only.

FIG. 24
FIG. 24.Λ T measure for the electric energy within the energy windows defined by the mean energy Ē = 23.5 and the window sizes ∆E = 1, ∆E = 1.3 and ∆E = 1.6 for g 2 = 0.8.

FIG. 26 .
FIG. 26.Internal energy U (T) as a function of temperature T for N = 19 (solid, blue) and N = 17 (dashed, black) plaquette chains with jmax = 1 2 .U (T) saturates at large T because the system contains only a finite number of degrees of freedom.All quantities are shown in lattice units.