Disorder-free localization in an interacting two-dimensional lattice gauge theory

Disorder-free localization has been recently introduced as a mechanism for ergodicity breaking in low-dimensional homogeneous lattice gauge theories caused by local constraints imposed by gauge invariance. We show that also genuinely interacting systems in two spatial dimensions can become nonergodic as a consequence of this mechanism. Specifically, we prove nonergodic behavior in the quantum link model by obtaining a rigorous bound on the localization-delocalization transition through a classical correlated percolation problem implying a fragmentation of Hilbert space on the nonergodic side of the transition. We study the quantum dynamics in this system by means of an efficient and perturbatively controlled representation of the wavefunction in terms of a variational network of classical spins akin to artificial neural networks. We identify a distinguishing dynamical signature by studying the propagation of line defects, yielding different light cone structures in the localized and ergodic phases, respectively. The methods we introduce in this work can be applied to any lattice gauge theory with finite-dimensional local Hilbert spaces irrespective of spatial dimensionality.

Introduction. Systems with local constraints play an important role in various physical contexts ranging from strongly correlated electrons [1] and frustrated magnets [2, 3] to fundamental theories of matter such as quantum electro-and chromodynamics [4], where constraints take the form of local gauge symmetries. The equilibrium properties of such systems have been extensively studied over the last decades, but only recently their nonequilibrium dynamics has moved into focus. In particular, local constraints have emerged as a new paradigm for ergodicity breaking, besides the two known archetypical scenarios caused by localization due to strong disorder or integrability. Systems with local constraints can exhibit rare nonergodic eigenstates, termed quantum many-body scars [5, 6], or extremely slow relaxation [7][8][9], whereas dipole conservation can prevent thermalization of large parts of the spectrum in one-dimensional fractonic systems [10][11][12][13]. A particularly generic mechanism for nonergodic behavior is hosted in lattice gauge theories (LGTs) where local constraints emerge naturally due to the local gauge symmetry, leading to an extensive number of local conserved quantities. Specifically, this can lead to the absence of ergodicity in 1D LGTs with discrete [14,15] and continuous [16] gauge symmetries or for higher-dimensional systems in the low-energy limit [17] or when they are free [18]. It has remained, however, a key challenge to rigorously identify non-ergodic behavior in genuinely interacting quantum matter beyond one dimension.
In this work we show that the 2D U (1) quantum link model features both localized and ergodic phases in the absence of disorder. Our proof relies on a mapping onto a classical correlated percolation problem providing a rigorous bound on the localization transition of the quantum model. We identify a distinguishing quantum-dynamical signature of the two phases by studying the propagation of an initial line defect, which leads to two different light cone structures. For the description of the nonequilibrium dynamics of the system we introduce variational classical networks (vCNs), which provide an efficient and perturbatively controlled representation of the quantum many-body wave function in terms of a network of classical spins akin to artificial neural networks (ANNs). The introduced methods can be applied to any LGT with finite local Hilbert spaces irrespective of dimensionality.
Quantum link model. We study the 2D U (1) quantum link model (QLM) [19,20], which has been introduced as a descendant of lattice quantum electrodynamics with spin-1/2 gauge degrees of freedom. In the QLM the spins S r,µ reside on the links of a square lattice connecting vertices r = (x, y) and r + µ (here µ =î,ĵ is one of the two unit vectors of the lattice, Fig. 1(a)), with the Hamiltonian: (1) The sums run over all plaquettes , U = S + r,î S + r+î,ĵ S − r+ĵ,î S − r,ĵ induces a collective flip of all spins on plaquette , and S ± r,μ denote the raising and lowering operators. The first (potential) term counts the number of flippable plaquettes and the second (kinetic) term induces coherent dynamics. For what follows, we will consider periodic boundary conditions and the case of a strong potential term with J/λ = −0.1. The QLM not only appears in the context of high-energy physics, but also shares strong connections to condensed matter systems featuring quantum spin ice phases [21,22] or Spins pointing → or ↑ correspond to S z = +1 and ← or ↓ to S z = −1, respectively. Kinetics is introduced by plaquetteflip operators U , U † (shown for the darkened central plaquette) whenever the spins on a plaquette are oriented clockwise or counterclockwise. Flippable plaquettes are denoted by circular arrows. Background charges with nonzero in-our outflow of electric field at a given vertex are indicated by red a blue dots. b) Spatiotemporal buildup of quantum correlations S z r,y (t)S z r+d,y (t) ≡ S z 0 (t)S z d (t) (d = dî) starting from |ψ(α = 0) = |→ in the localized phase of the QLM for J/λ = −0.1 and a system of size 20 × 20, i.e., 800 spins. quantum dimer models [1,23]. On the experimental side various proposals have explored the potential realization of the QLM in quantum simulators within the last years [24,25].
The local gauge symmetry of the QLM is generated by the operators G r = µ (S z r,µ − S z r−µ,µ ) counting the total inflow of the electric field to the vertex r. Since [G r , H] = 0 for all lattice points and [G r , G r ] = 0, eigenstates of H can be classified by the respective eigenvalues q r ∈ {−2, −1, 0, 1, 2} of G r . The set of q = {q r } defines the so-called superselection sector of states |ψ q with G r |ψ q = q r |ψ q , so that each of the q r can be given a physical meaning in terms of static background charges located at r [16]. The QLM further has global conserved quantities Φ x = y S z r,î , Φ y = x S z r,ĵ , which define the flux sectors.
Disorder-free localization. The existence of these sectors, protected by gauge invariance, can lead to an unconventional scenario for ergodicity breaking. Consider a homogeneous superposition state |ψ = q C q |ψ q involving many superselection sectors. As the Hamiltonian and typical observables are block-diagonal, i.e., H |ψ q = H q |ψ q , the expectation values of an operator O during dynamics become equivalent to O(t) = q |C q | 2 ψ q |e iHqt Oe −iHqt |ψ q resembling an effective disorder average with the disorder strength determined by the random background charges in the typical superselection sectors [16]. This can, in principle, lead to non-ergodic behavior of O(t) , although both the initial state and the Hamiltonian are homogeneous leading to the notion of disorder-free localization [14].
Initial states for time-evolution. We now aim to characterize the nonequilibrium dynamics of the QLM for the following homogeneous initial states |ψ 0 : where |↑ i and |↓ i are the two basis states at link i. This state is distributed over all superselection sectors of the model.
(ii) |ψ 0 = |→ FF which is a projection of |→ to a single "fully-flippable" (FF) sector, defined as the zerocharge zero-flux sector. |→ FF is an equal-weighted superposition of all states from the FF sector (i.e. the Rokhsar-Kivelson state [1] for the FF sector).
The resulting dynamics for α = 0 is displayed in Fig. 1(b), where we monitor the spatiotemporal buildup of quantum correlations. We will identify the limited spatial propagation with nonergodic behavior below.
Variational classical networks. Calculating the dynamics of interacting quantum systems in 2D is an inherently hard problem without a general-purpose computational method available to date. Representing a generic quantum many-body state as |ψ = s ψ( s) | s requires, in principle, the storage of exponentially many amplitudes ψ( s) (here s = (s 1 , ..., s N )). Recently, it has been proposed to use networks of classical spins to solve this problem [26,27] by avoiding to store the ψ( s)'s. The amplitudes ψ( s) ≈ exp[H ( s, W)] are rather generated on the fly when needed via a complex classical spin model with Hamiltonian H ( s, W) determined by a set of couplings W between the involved spins. Here, we construct H ( s, W) using a perturbatively controlled expansion and extend the recently proposed classical networks [27] upon imposing an additional optimization principle. The resulting approach can be interpreted as encoding |ψ in an ANN with a specific simplified network structure.
Within the vCNs we perform an expansion around a classical limit, which in the case of the QLM is the potential term H 0 in (1).
By representing the evolution operator in the interaction picture For the remaining term s| W (t) |ψ 0 we perform a cumulant expansion for time-ordered exponential operators [27][28][29], which, e.g., to the first order yields: Here denotes the sum over all flippable plaquettes in the spin configuration s, and ω = −4, ..., 4 counts the difference between number of flippable plaquettes surrounding the given before and after its flip (i.e. λω gives the potential energy difference before and after the flip). For example, for the configuration in Fig. 1(a) we have ω = 3 − 1 = 2 for the central plaquette. Going beyond previous work [27] we promote λt and its functions such as t 0 dt e iλω t to ten variational parameters W k (t) = ( The local connectivity of the vCN is encoded in the function ω ( s). For the actual shown numerical simulations we use a second-order ansatz and more complex initial states (see [29,30]). The W k (t)'s are determined by a time-dependent variational principle translating quantum dynamics into a system of coupled classical differential equations [29]). We solve these equations using a 4th-order Runge-Kutta integrator with step size ∆t = 0.1λ −1 and sample the observables using Metropolis Monte Carlo (MC) with 10 6 sweeps at each time instance, with single spin-flip updates for |→ and plaquette flips for |→ FF .
While our approach is numerically stable and therefore doesn't face some challenges appearing in ANNs [31], it shares its own limitations due to its perturbative construction, which is guaranteed to work only up to times t |1/J|. We find, however, that the errors remain perturbatively controlled up to much longer times as a consequence of the variational optimization. This can be verified, since the method provides a self-contained way of tracking the error, not referring to any reference solution. We present the details of the error analysis in [29,30] together with benchmarks.
Localized and ergodic dynamics. Using the vCNs we now compute nonequilibrium dynamics in the QLM. We start by studying the spatiotemporal buildup of quantum correlations, measured via S z 0 (t)S z d (t) , upon initializing the system in state |ψ 0 = |→ . The result is shown in Fig. 1(b), where one can see that correlations emerge only over a limited spatial distance suggesting nonergodic behavior. We proceed by further corroborating this observation by other measures.
Namely, we study energy transport in the QLM by creating initial conditions with a spatial energy inho- mogeneity in the form of a line defect with subextensive energy contribution and use the character of energy propagation to distinguish between ergodic and localized dynamics. Concretely, we consider the two initial conditions |ψ 0 = |→ or |→ FF upon applying in addition P = ∈C0 [1 + (U + U † ) 2 ] along all plaquettes in column d = 0; here C d denotes the set of plaquettes in column d. In Fig. 2 . Further, H av = H /L denotes the expected H d (t) in the long-time limit when the system is thermalizing (L is the number of columns).
Comparing Figs. 2(a) and (c) we observe that the dynamics differs qualitatively for the two initial conditions, although the Hamiltonian parameters are identical. While for |ψ 0 = |→ energy transport is highly suppressed and only visible on short distances (Fig. 2(a)), the opposite happens for |→ FF . This becomes even more apparent in Figs. 2(b,d), where ε d (t) relative to the initial value ε d (0) is shown, therefore more directly highlighting energy propagation. While for |→ FF we identify a linearly propagating front, for |→ we observe a strong bending. We argue below that this front for |→ can extend only to a finite region as a consequence of disorder-free localization. Bound on quantum dynamics by unconventional percolation. The qualitative difference in the quantum dynamics for the initial states |→ and |→ FF originates from a dynamical transition, which one can study systematically upon tuning the parameter α for the initial state (2). For this purpose, we employ an unconventional correlated classical percolation problem [32] and establish a bound on the quantum localized-ergodic transition in the QLM providing a numerical proof for an extended nonergodic phase as a consequence of disorder-free localization.
We illustrate the idea for the initial state |→ , distributed over all superselection sectors. Consider a typical (random) sector from this distribution (Fig. 3(a)). Such sector exhibits many background charges q r whenever the "two-in two-out" rule at vertex r is violated. Importantly, these background charges (constants of motion by gauge invariance) impose strong kinetic constraints. For instance, q r = ±2 implies that neighboring spins either all point inwards or outwards, hence the adjacent plaquettes remain unflippable forever. The influence of q r = ±1 charges is more subtle. They make at least 2 adjacent plaquettes unflippable, while their positions might change over time.
The question we address now is whether these constraints are so strong to fragment the square lattice into sets of kinetically disconnected islands or whether one can contain an extensive (percolating) connected cluster. For that purpose we study an unconventional percolation problem using an infinite-temperature classical MC simulation. We start from the initial condition (2), sampling a random basis state (and thus a sector) with a distribution set by the amplitudes in |ψ(α) . Then we determine which parts of the systems are kinetically connected, using MC search with random plaquette flips. The simulation is stopped when every plaquette is flipped either 0 or more than some fixed threshold (= 100) number of times (or after 10 11 MC steps if this condition is still not satisfied). As a result we find the number of performed flips for each plaquette (Fig. 3(b)). Repeating this procedure for different initial configurations at a given α and scanning α, we finally obtain the percolation probability (Fig. 3(c)). Most importantly, one can observe a clear evidence for a percolation threshold α c ≈ 0.25. Although the simulation termination condition is chosen such as to minimize the number of potentially missed "weak connections" between flippable clusters, we cannot exclude the possibility of such misses. While we don't expect a significant impact deep in the respective phases, this caveat might become important in the vicinity of α c , preventing us to obtain a precise value of α c and to study the critical behavior.
Since α c > 0, the initial state |ψ(α = 0) = |→ corresponds to the classically non-percolating side of the transition, while from α c < π/4 it follows that state |ψ(α = π/4) = |FF and all other states from the FFsector (including |→ FF ) lie on the percolating side. This classical threshold is imprinted in the quantum dynamics and ultimately leads to the strong localization observed in propagation of correlations ( Fig. 1(b)) and of the energy (Fig. 2(a),(b)) for |→ . For the FF-sector state |→ FF there is no percolation constraint, which allows propagation of the signal to long distances (Fig. 2(c),(d)). This analysis sets a lower bound onto the critical value of the quantum transition α (q) c , since the quantum system might be still localized due to interference even on the classically percolating side.
Summary and outlook. We have shown that genuinely interacting 2D homogeneous LGTs can become nonergodic as a consequence of disorder-free localization. This is all the more surprising as many-body localization is predicted to be unstable in 2D at elevated energy den-sities [33], implying that gauge invariance represents a more robust mechanism of ergodicity breaking compared to conventional disorder. The key element of our analysis is a bound on the localization-delocalization transition based on a classical correlated percolation problem implying a strong fragmentation of Hilbert space into kinetically disconnected regions. Both the percolation analysis as well as the introduced variational classical networks can be directly applied to other quantum many-body systems with finite-dimensional local Hilbert spaces independent of dimensionality, such as 3D quantum spin ice systems, which might be an interesting scope of the developed techniques in the future. Further, it might be interesting to explore how the quantum and classical percolation thresholds are related to each other as well as to determine their respective critical behaviors, and whether the disorder-free localization scenario holds also in the presence of matter degrees of freedom. Our theoretical analysis appears within reach of future experiments: significant efforts in the last years have explored routes to realize the QLM model experimentally in systems of Rydberg atoms [24,25] as a next step after the recent experimental advances on 1D LGTs [34][35][36][37][38].
Acknowledgments. We are grateful to H. Burau, J. In this Supplementary Section we outline the procedure of construction of a variational classical networks for the quantum link model, which can be straightforwardly generalized to other quantum models with finite local Hilbert spaces [S1].
If the Hamiltonian of the model can be written as H = H 0 +V , where H 0 is diagonal in the computational basis and V H 0 is a small non-diagonal perturbation, then in order to describe the evolution of the wave function we can use a perturbative treatment based on the cumulant expansion [S2]. Moreover, we can build on top of this perturbative description and construct adequate variational Ansätze for the wavefunction.
For the studied quantum link model we have , where the sums run over all plaquettes . In the perturbative limit it is convenient to work in the interaction representation and rewrite the evolution operator as e −iHt = e −iH0t W (t), where W (t) = T exp(−i t 0 dt V (t )). Let further the initial state be |ψ 0 = s ψ 0 ( s) | s (i.e. ψ 0 ( s) = s|ψ 0 ), where | s = |s 1 , s 2 , ... ; and the time-evolved state |ψ(t) = e iHt |ψ 0 = s ψ( s, t)| s .
Below we explicitly describe the construction of the first-order ansatz and sketch an analogous derivation for the second-order ansatz.

A. First-order variational ansatz
Within the above settings the coefficients for the time-evolved state |ψ(t) in the lowest-order in the cumulant expansion [S3] in the perturbation V can be obtained as where E s is the eigenvalue of the unperturbed Hamiltonian: H 0 | s = E s | s and we assumed that s|ψ 0 doesn't contain any special smallness. To the lowest order in V we can reexponentiate the result and obtain We recast our wavefunction coefficient in the form of an effective Hamiltonian of a classical spin model H : ψ( s, t) = e H ( s,t) convenient for Monte Carlo sampling: In order to proceed we need to find V (t), which can be expressed using U (t). Equation of motion for operators U (t) ≡ e iH0t U e −iH0t in the interaction representation is given by Only four plaquettes (indexed by letters a, b, c, d, Fig. S1) adjacent to the plaquette under consideration contribute to the commutator [H 0 , U ] = λΩ U , where Ω is a diagonal operator  plaquettes (a, b, c, d) for a given one (central, darkened). Recall that spin directions →, ↑ correspond to S z = +1, while ← and ↓ correspond to S z = −1 and note that S1 ≡ Sa3, S2 ≡ S b4 and so on. and P + and P − are projectors to S z = +1 and S z = −1 correspondingly: P ± = (1 ± S z )/2. Thus the equation of motion (S4) transforms to −i d dt U (t) = λΩ U (t) and can be readily solved by (note that Ω commutes with U since they do not contain the same spins). Now we're ready to substitute V (t) = J (U (t)+U † (t)) in (S3). In order to do so it's convenient to introduce the diagonal operator F = U U † −U † U ; so its matrix element F ( s) = +1, −1, 0 is the flippability of the plaquette in the state s for counterclockwise (+1) or clockwise-oriented (−1) flippable, or non-flippable plaquettes (0) respectively. By s we denote configuration which differs from s by the flip of single plaquette: | s ≡ (U + U † )| s . It's also convenient to denote ω ( s) = F ( s)Ω ( s). It takes values ω = −4, ..., 4 and counts the difference between number of flippable plaquettes surrounding the given before and after its flip (i.e. λω gives the potential energy difference before and after the flip). Then s|( where prime over the sum indicates summation over only flippable plaquettes. This gives us the lowest-order term of the perturbative expansion of ψ( s, t). Now we want to construct a variational ansatz that should at least reproduce this perturbative result (S8) and hopefully improve it. We use the following strategy: we promote all functions of time to the variational parameters W k in such a way that the effective classical Hamiltonian H ( s, W) was a linear function of them. The non-trivial time dependence sits in the integrals Function f (ω , t) depends on t and on the discrete parameter ω that takes values −4, .., 4. We promote this function to the most general function of the same form introducing nine first-order variational parameters W 4 . In the initial state all variational parameters are zero. Finally along with the 0-th order variational parameter W (0) we obtain the following ansatz Such ansatz works for translationally-invariant initial states. For the inhomogeneous initial state we replicate this set of variational parameters for each plaquette or each column depending on the spatial symmetry of the state. Note that the ansatz (S11) can be explicitly recast to a classical Ising-like spin model with multiple (up to 16) spin interaction terms, W (1) ω ( s) = 4 ω=−4 δ ω,ω ( s) W ω , where for example δ ω=4,ω ( s) = P − 1 P − 2 P + 3 P + 4 P + a1 P + a2 P − a4 P + b1 P + b2 P − b3 P + c2 P − c3 P − c4 P + d1 P − d3 P − d4 + h.c., with P ± i = (1 ± S z i )/2. We emphasize that functions ω ( s) and the corresponding variational parameters W (1) ω ( s) encode this rather involved spin model with a complex network of local spin connections (i.e. interaction terms) in a much simpler way, while containing exactly the same information. We refer to the ansatz (S11) as to a first-order variational classical network (vCN) for the quantum link model. In the main text we show its simplified version for the case ψ 0 ( s) = const.

B. Second-order variational ansatz
Here we sketch the analogous derivation for the second-order variational ansatz used for the calculations in the present work.

III. BENCHMARKS
In this Supplementary Section we provide benchmarks for the variational classical network method. The approximate time-evolved state is constructed by using the 2nd order perturbational ansatz (introduced in Suppl. Sec. I), which is evolved using TDVP (Suppl. Sec. II). During the evolution we sample various observables using the Monte Carlo method. We compare the vCN results with the results of exact diagonalization for the fully flippable sector of 10 × 2 system (the dimension of the Hilbert space is 17906). Below we present a) single-plaquette observables, b) correlation functions, and c) error production rate for translationally-invariant initial state |→ FF and a non-uniform initial state P |→ FF defined in the main text. Figures S2-S4 show that up to time λt |λ/J| ≈ 10, vCN method works very well for both the single-plaquette observables (Figs. S2-S3) and for the long-distance correlators (Fig. S4), which is also manifested in the fact that the error production rate remains small till λt 10 and starts significantly grow after that (Fig. S5). Nevertheless we see that vCN can capture the single plaquette observables up to much longer times, especially the average behavior with smoothed out high-frequency oscillations (note much smaller vertical scale for d = 1 − 5 compared to d = 0 in Fig. S3).