Realistic scheme for quantum simulation of $\mathbb{Z}_2$ lattice gauge theories with dynamical matter in $(2+1)$D

Gauge fields coupled to dynamical matter are ubiquitous in many disciplines of physics, ranging from particle to condensed matter physics, but their implementation in large-scale quantum simulators remains challenging. Here we propose a realistic scheme for Rydberg atom array experiments in which a $\mathbb{Z}_2$ gauge structure with dynamical charges emerges on experimentally relevant timescales from only local two-body interactions and one-body terms in two spatial dimensions. The scheme enables the experimental study of a variety of models, including $(2+1)$D $\mathbb{Z}_2$ lattice gauge theories coupled to different types of dynamical matter and quantum dimer models on the honeycomb lattice, for which we derive effective Hamiltonians. We discuss ground-state phase diagrams of the experimentally most relevant effective $\mathbb{Z}_2$ lattice gauge theories with dynamical matter featuring various confined and deconfined, quantum spin liquid phases. Further, we present selected probes with immediate experimental relevance, including signatures of disorder-free localization and a thermal deconfinement transition of two charges.


I. INTRODUCTION
It has been a long sought goal to faithfully study lattice gauge theories (LGTs) with dynamical matter in the realm of strong coupling.Since their discovery, Z 2 LGTs have sparked the interest of physicists from various different fields including high-energy [1], condensed matter [2][3][4] or biophysics [5].The seminal work by Fradkin and Shenker [6] in 1979 predicted the existence of two phases in their model, in which Z 2 charged particles are either confined or deconfined in (2 + 1)D.This insight made it a particularly promising candidate theory that could capture some of the essential physics of quark confinement in QCD [1] while hosting a much simpler gauge group.Likewise, it provides one of the most fundamental instances of the Higgs mechanism.Since then the study of Z 2 LGTs has inspired physicists because of their intimate relation to topological order [7], quantum spin liquids [8,9] and quantum information [10], to name a few.While the physics of these models could give insights into outstanding problems, e.g., how to define confinement in the presence of dynamical matter, the numerical (e.g.Refs.[11][12][13][14][15]) and experimental exploration is at the same time very challenging beyond (1 + 1)D (e.g.Refs.[16][17][18][19]).
The experimental developments over the past years have driven the field of analog quantum simulation towards exploring many-body physics in system sizes out of reach for any numerical simulation and offering a new toolbox to approach complex, physical phenomena such as quantum spin liquids [20].The difficulty to implement gauge constraints and robustness against everpresent gauge-breaking errors in analog quantum simulators, however, has hindered the field to push forward into the aforementioned direction and a scalable, reliable implementation of LGTs with dynamical matter in (2+1)D remains a central goal.
The rich structure of gauge theories emerges from locally constraining the Hilbert space.This constraint can be formulated by Gauss's law, which requires all physical states |ψ to fulfill Ĝj |ψ = g j |ψ .For the Z 2 LGT with dynamical matter (Z 2 mLGT) we consider in this work the symmetry generators Ĝj are given by Ĝj = (−1) nj i: i,j where nj = â † j âj is the number operator for (hard-core) matter on site j and the Pauli matrix τ x i,j defines the electric field on the link between site i and j; hence g j = ±1.Our starting point throughout this work are link and site qubits on a two-dimensional honeycomb lattice, see Fig. 1a.
We propose to realize matter and link variables as qubits, implementable e.g. by the ground |g and Rydberg |r states of atoms in optical tweezers [20][21][22][23][24][25], see Fig. 1a-c.Thus, the product in Eq. ( 1) measures the parity of qubit excitations of matter and links around vertex j.
By encoding the degrees-of-freedom in qubits the enlarged Hilbert space contains physical (g j = +1) and LGT with matter ℤ 2 electric field matter < l a t e x i t s h a 1 _ b a s e 6 4 = " q / L h 9 g i q E U N Q 4 6 P i C 9 e h i p z 1 A vertex contains matter âj qubits (blue) and shares link τ x i,j qubits (red) with neighboring vertices.All qubits connected to a vertex interact pairwise with strength 2V .In a Rydberg atom array experiment the qubits are implemented by individual atoms in optical tweezers, which are assigned the role of matter or link depending on the position in the lattice.Here, the ground-and Rydberg state of the atoms, |g and |r , encode qubit states, which are coupled by an off-resonant drive Ω to induce effective interactions.To realize equal strength nearest neighbor, two-body Rydberg-Rydberg interactions, the matter atoms can be elevated out of plane.In panel b) we introduce the notation for the Z2 mLGT, for which the Hilbert space constraint is given by Gauss's law Ĝj = +1.We illustrate the electric field τ x i,j = +1 (τ x i,j = −1) with flat (wavy) red lines and the matter site occupation n j = 0 (n j = 1) with empty (full) blue dots.Panel c) shows the notation for the QDM subspace with exactly one dimer per vertex.Panel d) illustrates how the distinct subspaces are energetically separated by the LPG term V Ŵj .The two quantum dimer subspaces are disconnected when the matter is static, which can be exactly realized by the absence of matter atoms in panel a) and setting (2â † j âj − 1) = ±1 in V Ŵj .
unphysical (g j = −1) states: The latter do not fulfill Gauss's law.Since any local perturbations present in a realistic quantum simulation experiment mix the two subspaces, quantum simulations can become unreliable, effectively breaking gauge-invariance.Nevertheless, by energetically separating the physical from unphysical states transitions into the latter can be suppressed and the gauge structure emerges from the enlarged Hilbert space.
The simplest way, theoretically, to achieve such gauge protection, is by adding −V j Ĝj to the Hamiltonian with large V > 0 [26][27][28].But since this would require strong four-body interactions, it is experimentally not feasible in current experimental platforms.
Here we demonstrate that simple two-body Ising-type interactions, which are readily available in e.g.Rydberg tweezer arrays [20][21][22][23][24][25], combined with longitudinal and weak transverse fields provide a minimal set of ingredients which allow to robustly implement a variety of LGTs with dynamical matter [9].The scheme we propose not only offers inherent protection against arbitrary gaugebreaking errors; it also provides a surprising degree of flexibility, including cases with global conserved particle number, global number-parity conservation, and quantum dimer models on a bipartite lattice which map to U (1) gauge theories.
In the following, we show that readily available Isingtype two-body interactions, in addition to local fields, are sufficient to protect Gauss's law on experimentally relevant timescales by employing the so-called local pseudogenerator (LPG) method [29].Moreover, we show that the proposed protection scheme provides a generic means to engineer a variety of effective Z 2 mLGT Hamiltonians by weakly driving the qubits.As an example, we demonstrate how this allows to realize the celebrated Fradkin-Shenker model [6], and discuss the phase diagrams of several related effective Hamiltonians.Finally, we elaborate on some realistic experimental probes that we view as most realistic in state-of-the-art quantum simulators.

II. RESULTS
Local pseudogenerator on the honeycomb lattice.-Themain ingredient of the experimental scheme proposed in this Article is the local pseudogenerator (LPG) interaction term V Ŵj .As shown in Fig. 1a, V Ŵj consists of equal-strength 2V interactions among all qubits (matter and gauge) around vertex j, taking the form We assume that V defines the largest energy scale in the problem, which separates the Hilbert space into constrained subspaces.This overcomes the most challenging step, imposing different gauge constraints in the emerging subspaces (Supplementary note 1).We obtain three distinct eigenspaces of the LPG term: 1) Two (distinct) quantum dimer model (QDM) subspaces with static matter at low-energy, 2) physical states of a Z 2 mLGT at intermediate energies, and 3) trivial, polarized states at high energy, see Fig. 1b-d.
The LPG method requires that V Ŵj acts identical to the full protection term on all physical states in the target gauge sector, i.e.Ŵj |ψ = Ĝj |ψ .For unphysical states, instead, the LPG term splits into many manifolds that can be energetically above and below the target sector [29].This construction allows to reduce experimental complexity from four-to two-body interactions.
Experimentally, we propose to implement strong LPG terms in the Hamiltonian such that quantum dynamics are constrained to remain in LPG eigenspaces by large energy barriers enabling the large-scale quantum simulation of Z 2 mLGTs in (2+1)D.To introduce constraint-preserving dynamics within the LPG subspaces, the latter are coupled by weak on-site driving terms of strength Ω V as discussed below.Through the constrained dynamics, a Z 2 mLGT emerges in an intermediate-energy eigenspace of V Ŵj , which is accessible in quantum simulation platforms and which dis-tinguishes our work from previous studies on emergent gauge symmetries, e.g.[30][31][32].
The LPG method is built upon stabilizing a highenergy sector of the spectrum, which comes with the caveat that a few unphysical states are resonantly coupled when considering the entire lattice.In particular, there is a subset of unphysical states that violate Gauss's law on four vertices with energy lowered on three vertices and raised on one vertex; hence these states are on resonance with physical states.However, numerical simulations in small systems suggest that these gauge-breaking terms only play a subdominant role and gauge-invariance remains intact (Supplementary note 2).
Ultimately, the problem of resonances with a few unphysical states can be remedied by promoting V → V j to be site-dependent such that high-energy sectors can be faithfully protected [33,34] against potential gauge noninvariant processes described above (see Methods section).Site-dependent protection terms do not require any additional experimental capabilities in our protocol described below.Even more, experimental imperfections inherently give disorder stabilizing the gauge sectors further.It is also important to note that the presence of only weak disorder (compared to the energy scale V ) is enough, which does not alter the effective couplings in the emergent gauge-invariant effective Hamiltonian.
In the following, we introduce the microscopic model that we propose to implement in an experiment.From the microscopic model, effective Hamiltonians for the Z 2 mLGT and QDM subspaces can be derived by a Schrieffer-Wolff transformation (Supplementary note 2 and 4).On realistic timescales of experiments, the effective models are gauge-invariant by construction and studied further below.
Experimental realization in Rydberg atom arrays.-Here,we propose the microscopic model Ĥmic which can be directly implemented in state-of-the-art Rydberg atom arrays in optical tweezers, see Fig. 1a.
The constituents are qubits, which can be modeled by the ground |g and Rydberg |r states of individual atoms.As shown in Fig. 1a, we label the atoms as matter atom or link atom depending on their position on the lattice.The Z 2 gauge structure then emerges from nearest-neighbor Ising interactions V realized by Rydberg-Rydberg interactions and hence the real space geometric arrangement plays a key role.The dynamics is induced by a weak transverse field Ω m (Ω l ), which corresponds to a homogeneous drive between the ground and Rydberg states of the matter (link) atoms.Moreover, tunability of parameters defining the phase diagram is achieved by a longitudinal field or detuning ∆ m (∆ l ) of the weak drive.
The interesting physics emerges in different energy subsectors of the LPG protection term ∝ V Ŵj in Eq. ( 2); in particular the Z 2 mLGT is a sector in the middle of the spectrum of Ĥmic .The suitability for Rydberg atom arrays comes from the flexibility in geometric arrangement required for the LPG term as well as from the natural en-ergy scales V Ω in the system, which we use to derive the effective models below, see Eqs. ( 4) and (5).
Matter atoms âj form the sites of a honeycomb lattice and we map the empty |n j = 0 (occupied |n j = 1 ) state on the ground state |g j (Rydberg state |r ) j ) of the atoms.Link atoms τ x i,j are located on the links of the honeycomb lattice, i.e. a Kagome lattice, and analogously we map the τ x i,j = +1 (τ x i,j = −1) state on the atomic state |g i,j (|r i,j = â † i,j |g i,j ).Moreover, we want the matter and link atoms to be in different layers and those layers should be vertically slightly apart in real space to ensure equal two-body interactions between matter and link atoms (Supplementary note 5).Using the out-of-plane direction has the advantage that it only requires atoms of the same species and with the same internal states.However, the equal strength interaction can also be achieved in-plane by using e.g. two atomic species or different (suitable) internal Rydberg states for the matter and link atoms.
We first propose a non gauge-invariant microscopic Hamiltonian from which we later derive an effective model with only gauge-invariant terms.To lowest order in perturbation theory and on experimentally relevant timescales, the system evolves under an emergent gaugeinvariant Hamiltonian.The microscopic Hamiltonian is given by Ĥmic where bosonic operators â † j and â( †) i,j annihilate (create) excitations on the matter and link atoms, respectively; Ŵj is the LPG term introduced in the main text Eq. ( 2).The last two terms describe driving of matter (|g j ↔ |r j ) and link atoms (|g i,j ↔ |r i,j ) in the rotating frame.Rewriting (3) in the atomic basis yields Rydberg-Rydberg interactions of strength 2V and renormalized, large detunings ∆m = −3V + ∆ m and ∆l = −3V + ∆ l .In a Rydberg setup the driving terms can be realized by an external laser, which couples |g ↔ |r , while the detunings ∆ m , ∆ l of the laser relative to the resonance frequency controls the electric field ∆ l and chemical potential ∆ m in the rotating frame.
In the limit Ω m , Ω l V , the energy subspaces defined by the LPG term V Ŵj , Eq. ( 2), are weakly coupled by the drive to induce effective interactions and it is convenient but not required to choose Ω m = Ω l = Ω.The Z 2 mLGT emerges as an intermediate-energy eigenspace of the LPG term V Ŵj .The effective interactions in the constrained Z 2 mLGT and QDM subspaces of Ŵj can be derived by a Schrieffer-Wolff transformation (Supplementary note 2 and 4) and yielding the models discussed in the next section.
In the experiment we propose, the Rydberg-Rydberg interactions are not only restricted to nearest neighbours but are long ranged.We emphasize that beyond nearest neighbour interactions are inherently gauge invariant and hence do neither influence the LPG gauge protection scheme nor the Schrieffer-Wolff transformation.However, the long-range interactions can have strong influence on the Z 2 invariant dynamics.While the interaction strength decreases as 1/R 6 , where R is the distance between atoms, the interaction is still comparable to the effective perturbative dynamics (Supplementary note 5).We note that the dynamics might be slowed down but the qualitative features of the Z 2 mLGT remain intact.
Effective Z 2 mLGT model.-Amodel is locally Z 2 invariant if its Hamiltonian Ĥ commutes with all symmetry generators Ĝj , i.e. [ Ĥ, Ĝj ] = 0 for all j.This ensures that all dynamics is constrained to the physical subspace without leaking into unphysical states.In Eq. ( 2), the target sector is g j = +1 for all j but our scheme can be easily adapted for any {g j } j (Supplementary note 1).
In the presence of strong LPG protection, the system is energetically enforced to remain in a target gauge sector and unphysical states are only virtually occupied by the drive Ω.To be precise, resonant couplings to unphysical sectors are suppressed by the (experimentally feasible) disorder protection scheme discussed above and in the Methods section.Otherwise emergent gauge-breaking terms appear in third-order perturbation theory.However, in small systems we have numerically confirmed that even without disorder in the LPG terms Gauss's law is well conserved (Supplementary note 2), which in larger systems we expect to crossover to an approximate gauge invariance.In the following we assume disorder protection or small systems, where leading order gaugebreaking terms are absent or can be neglect, respectively.
For the proposed on-site driving terms discussed above and shown in Fig. 1a, we derive the following effective Hamiltonian from the microscopic model (3) in the intermediate-energy LPG eigenspace (Supplementary note 2): The first terms in Eq. ( 4) describe gauge-invariant hop-ping of matter excitations with amplitude t and (anoma- Quantum-matter + gauge We show two qualitative sketches of phase diagrams for the effective model (4).In panel a), we consider U (1) matter (∆1 = ∆2 = 0) coupled to a dynamical Z2 gauge field as discussed in the main text.Along the vertical direction the filling is tuned, which yields an even (odd) Z2 pure gauge theory in the vacuum (Mott insulator) illustrated by the grey regions.In between the matter and gauge degrees-of-freedom interplay, for which we examined the limiting cases.Above the deconfined region, we expect a superfluid regime (yellow), while above the confined region composite mesons of Z2 charges may condense (red).In panel b), we show the phase diagram for an Ising Z2 LGT as proposed by Fradkin and Shenker [6].The 2D quantum Hamiltonian of the Ising Z2 mLGT has equal hopping t and pairing ∆1 strength and can thus be mapped on a classical 3D Ising theory.Because our model with quantum Z2 matter coupled to dynamical Z2 gauge fields has slight anisotropy between hopping and pairing, t = ∆1, as well as additional anomalous pairing terms ∆2, the classical mapping can only work approximately.We anticipate that the phase diagram should be qualitatively very similar to panel b).lous) pairing ∝ ∆ 1 (∝ ∆ 2 ).The term ∝ J is the magnetic plaquette interaction on the honeycomb lattice.The last two terms are referred to as electric field term h and chemical potential µ, respectively.Note that deriving Hamiltonian (4) from the microscopic model in Eq. ( 3) yields additional higher-order terms ∝ τ x τ x , τ x n, etc.In the effective model Ĥeff Z2 we treat these higher-order terms on a mean-field level of the electric field and matter density (Supplementary note 2).Moreover, we emphasize that the effective model is solely derived from the microscopic Hamiltonian, which only requires a simple set of one-and two-body interactions between the constituents.
For any site j, one can take âj → −â j and τ z i,j → −τ z i,j ; hence the effective Hamiltonian (4) has a local Z 2 symmetry, [ Ĥeff Z2 , Ĝj ] = 0 ∀j, qualifying it as Z 2 mLGT in (2 + 1)D.In particular, in our proposed scheme we do not have to apply involved steps to engineer Z 2 -invariant interactions but rather we exploit the intrinsic gauge protection by dominant LPG terms, which enforces any weak perturbation to yield an effective Z 2 mLGT.This approach also inherently implies robustness against gauge-symmetry breaking terms in experimental realizations.
In the following, we discuss the rich physics of the effective model (4).However, due to the complexity of the system, it is challenging to conduct faithful numerical studies in extended systems.As a first step, we examine well-known limits of the model and conjecture T = 0 phase diagrams of the effective Hamiltonian when the Z 2 gauge field is coupled to U (1) or quantum-Z 2 dynamical matter, respectively.We note that the strength of the plaquette interaction can only be estimated (Supplementary note 2) and competes with the long-range Rydberg interactions.Moreover, the disorder protection scheme underlying the derivation of the effective Hamiltonian ensures gauge-invariance of the leading order contributions but higher-order gauge breaking terms can in principle appear and affect the physics at very long timescales.
Our effective model describes the physics of experimental system sizes and timescales; the efficiency of the LPG gauge protection in the thermodynamic limit is a subtle open question.Hence, in the following we discuss phases of the effective model (4) that may (or may not) emerge from the microscopic model (3).
U (1) matter.-Byfixing the number of matter excitations in the system, i.e. ∆ 1 = ∆ 2 = 0 in Hamiltonian (4), the model has a global U (1) symmetry of the matter (hard-core) bosons, which can be achieved by choosing the detuning at the matter sites ∆ m comparable to V in our proposed experimental scheme Eq. ( 3).Here, we consider the phase diagram when the filling of matter excitations is controlled by the chemical potential µ.To map out different possible phases, we fix the hopping t and study limiting cases.
First, we consider the pure gauge theory with no matter excitations (µ → −∞), see Fig. 2a (bottom).The Hamiltonian then reduces to the pure Ising LGT [2] with matter vacuum -an even Z 2 LGT.The dual of this model exhibits a continuous (2+1)D Ising phase transition, corresponding to a confined (deconfined) phase below (above) a critical (J/h) c , respectively [2,4].At the toric code point (J/h = ∞) the system is exactly solvable [35] and the gapped ground state has topological order.
Because for J/h = ∞ the gauge field has no fluctuations, we can fix the gauge by setting τ z i,j = +1 and map out the pure matter theory in Fig. 2a (right).For finite µ we find a model with free hopping of hard-core bosons, for which the filling can be tuned by changing the chemical potential µ.Hence, for increasing µ and results based on the square lattice [36,37] we expect two continuous phase transitions: vacuum-to-superfluid and superfluid-to-Mott insulator.The Mott insulator phase is an odd Z 2 LGT because the matter is static and acts as background charge and thus can be treated as a pure gauge theory with g j = −1 [9].In the opposite limit J/h = 0, the same Mott state gives rise to a hard-core quantum dimer constraint for the Z 2 electric field lines.On the square lattice, the quantum dimer model and odd Z 2 LGT exhibit a phase transition from a confined to deconfined phase [15].The honeycomb lattice and next-nearest neighbor Rydberg-Rydberg interactions might feature additional symmetry-broken phases.Hence it requires a sophisticated analysis to map out the substructure of the Mott insulating phase in Fig. 2a.
In the limit of low fillings and small but finite J/h 1, the matter excitations form two-body mesonic bound states [15], which are Z 2 -charge neutral and can be considered as point-like particles.We can derive an effective meson model yielding hard-core bosons on the sites of a Kagome lattice (Supplementary note 3).
At T = 0 and sufficiently low densities, the mesons can condense and spontaneously break the emergent global U (1) symmetry associated with meson number conservation.To determine the phase boundary of the meson condensate, we consider a single pair of matter excitations doped into the vacuum.This pair cannot alter the pure gauge phases and thus the two charges can be considered as probes for the (de)confined regime.For the latter, the matter excitations are bound into mesons, in contrast to free excitations above the deconfined regime.Hence, the effective description of bound mesonic pairs breaks down at the phase transition of the pure gauge theory indicating the phase boundary of the meson condensate phase at small filling.
At higher densities, dimer-dimer interactions and fluctuations of the gauge field play a role, requiring a more sophisticated analysis to predict the ground state.We emphasize that the rich physics in this model emerges from the gauge constraint generated by the LPG terms.Moreover, we note that by lifting the hard-core boson constraint, which is beyond our experimental scheme, the model maps onto a classical XY model coupled to a Z 2 gauge field [9].This model has been studied on the square lattice in the context of topological phases of matter [9] and high-Tc superconductivity [38][39][40], to name a few.
Classical mapping.-Fort = ∆ 1 and ∆ 2 = 0 the model is well-studied and maps onto a classical Ising lattice gauge theory coupled to Ising Z 2 matter [6].In our experimental proposal ∆ 1 and ∆ 2 cannot be independently tuned, but due to the relevance of the model and its proximity to our effective model we briefly summarize the most important results for the square lattice here, see Fig. 2b.
In the limit with frozen gauge fields (pure matter axis, J/h = ∞) the resulting pure matter theory corresponds to a transverse field Ising model with a global Z 2 symmetry, which maps to a classical 3D Ising model and exhibits a continuous phase transition.On the pure gauge axis (t/µ = 0) the model exhibits a topological phase transition without local order parameters [2].Instead, the scaling of non-local Wegner-Wilson loops with their area/perimeter distinguishes the confined from the deconfined phase.Remarkably, the pure gauge model is also dual to a classical 3D Ising model, rendering the pure gauge axis dual to the pure matter axis.The same pure gauge phases are realized for µ → −∞ in the case with U (1) matter.
For more general J/h, the model's self-duality yields a symmetry in the phase diagram, which allows to study the pure gauge and matter theory in Fig. 2b but does not reveal the interior away from the axis.Fradkin's and Shenker's accomplishment was to show the existence of two distinct, extended phases: the confined and deconfined "free charge" phase, which have been confirmed numerically [12,13].From today's perspective, the latter would be characterized as topological phase of matter in the toric code universality class.
Quantum-Z 2 matter.-Now,we consider the full effective Hamiltonian (4), where hopping and pairing are anisotropic t = ∆ 1 and the pairing strength can depend on the electric field configuration ∆ 2 = 0, and relate it to Fig. 2b.Here, the pure matter theory can no longer be mapped on the classical 3D Ising model.Hence, we introduce the term quantum-Z 2 matter, which emphasizes the matter's Z 2 symmetry group but points out that a mapping to a known classical model is lacking.
We note that close to the toric code point (J/h = ∞ and t/µ = 0) in Fig. 2b, the expectation value of the electric field vanishes, τ x i,j = 0, and thus in mean-field approximation the anomalous terms should be negligible and renormalize the pairing ∆ 1 → ∆1 .For the pure gauge theory it has been shown [11] that the expectation value τ x i,j continuously changes by tuning the electric field term h.Hence, by performing a mean-field approximation in the electric field, the quantum-Z 2 mLGT maps onto the classical Ising Z 2 mLGT (Supplementary note 2 C).
Due to its proximity to the Ising Z 2 mLGT and its common symmetries generated by the proposed LPG term, we anticipate that the phase diagram of the quantum-Z 2 mLGT shares all essential features of the Ising Z 2 mLGT as shown in Fig. 2b.
Quantum dimer model (QDM).-Rokhsarand Kivelson introduced the QDM in the context of high-T c superconductivity, which has the constraint that exactly one dimer is attached to each vertex [41,42].The QDM is an odd Z 2 LGT, i.e. a pure gauge theory with g j = +1 replaced by g j = −1 ∀j, with h → ∞, and its fundamental monomer excitations are gapped and can only be created in pairs.
Our proposed scheme allows to directly implement the gauge constraint of the QDM experimentally by preparing the system in the ground-state manifold of the LPG term as shown in Fig. 1b and d.Note that the LPG term splits the ground-state manifold into two distinct subspaces, QDM 1 and QDM 2 , which can be seen by entirely removing the matter atoms and setting nj = 0, 1 in Eq. ( 2), such that only the link atom Kagome lattice remains; hence it can be implemented in-plane.A dimer then corresponds to either τ x i,j = −1 (QDM 1 ) or τ x i,j = +1 (QDM 2 ).Due to the LPG protection the QDM subspaces are energetically protected and monomer excitations cost a finite energy 2V .
By weakly driving the system, the motion of virtual, gapped monomer pairs perturbatively induces plaquette terms of strength J QDM , and we can derive an effective model (Supplementary note 4) given by Here, the NNN link-link interaction K can be tuned by the blockade radius of the Rydberg-Rydberg interactions.
Experimental [20] and theoretical [32,[43][44][45] studies of QDMs in Rydberg atom arrays for different geometries and parameters regimes have shown to be an promising playground to probe Z 2 spin liquids.Our proposed setup is a promising candidate to further study QDMs due to its versatility and its inherent protection by the LPG term and the phase diagram of Hamiltonian (5) remains to be explored Here, we examine two limiting cases of Hamiltonian (5).For J QDM /K 1, the system is in the socalled plaquette phase [46], which is characterized by a maximal number of flippable plaquettes and resonating dimers.On the other hand, for J QDM /K 1 we find a classical Ising antiferromagnet on the Kagome lattice with NN and NNN interactions from the hard-core dimer constraint and K-term, respectively.
Experimental probes.-In the following, we discuss potential signatures of the rich physics that can be readily explored with the proposed experimental setup Eq. ( 3).
Disorder-free localization.-Recently,the idea of disorder-free localization (DFL), where averaging over gauge sectors induces disorder, has sparked theoretical interest [47,48].DFL is an example where the entire Z 2 mLGT Hilbert space participates in the dynamics including sectors with g j = +1.It has been demonstrated that the (2 + 1)D U (1) quantum link models can show DFL [49,50]; further it was proposed that in a (1 + 1)D Z 2 LGT, LPG protection leads to enhanced localization [34].However, experimental evidence is still lacking.The scheme we propose is suitable to experimentally study ergodicity breaking without disorder in a strongly interacting (2 + 1)D system with U (1) matter.
In Fig. 3a we show results of a small-scale exact diagonalization (ED) study using realistic parameters for the experimentally relevant microscopic Hamiltonian (Supplementary note 6).The system is prepared in two different initial states: 1) A gauge-invariant state |ψ inv , and 2) a gauge-noninvariant state |ψ ninv , both with (without) localized matter excitations in subsystem A (B).
We find distinctly different behaviours for the timeaveraged matter occupation imbalance between subsystem A and B (Supplementary note 6): While the gaugeinvariant state |ψ inv thermalizes, the gauge-noninvariant state |ψ ninv breaks ergodicity on experimentally relevant timescales.Experimentally much larger systems can be addressed.
Schwinger effect.-TheSchwinger effect describes the creation of pairwise matter excitations from vacuum in strongly-coupled gauge theories [51].Here, we use the Schwinger effect to test the validity of our LPG scheme.Starting from the microscopic model ( 3), we time-evolve the vacuum state with no matter excitations and extract the maximum number of created matter excitations in the initial gauge sector g j = +1 ∀j.As shown in Fig. 3b, by tuning the electric field and chemical potential we find resonance lines, where many matter excitations are produced in the system, and we verify that gauge-invariant processes dominate (Supplementary note 7).
Phase transitions in a ladder geometry.-Ourproposed scheme is suitable for any geometry with coordination number z = 3; hence one can experimentally study square ladders of coupled 1D chains.Here, we have examined the ground state of Hamiltonian (4) with U (1) matter using the density matrix renormalization group (DMRG) technique [52] (Supplementary note 8) on a ladder and we find signatures of a quantum phase transition.As shown in Fig. 3c, both the average density of matter excitations and the plaquette terms, which are experimentally directly accessible by projective measurements, change abruptly by tuning the electric field h indicating a transition into the vacuum phase.We emphasize that the ladder geometry is different from the (2 + 1)D model studied in Fig. 2a, however numerical simulations suggest the presence of a phase transition and hence the ladder geometry offers a numerically and experimentally realistic playground for future studies of our model.
Thermal deconfinement from string percolation.-Weexamine a temperature-induced deconfinement transition in a classical limit of our effective model ( 4), which neglects charge and gauge dynamics t = ∆ 1,2 = J = 0. We use Monte Carlo simulations on a 35×35 honeycomb lattice (Supplementary note 9).
To study thermal deconfinement, we consider exactly two matter excitations which, due to Gauss's law, have to be connected by a string Σ of electric field lines; i.e.Σ is a path of links with electric fields τ x i,j = −1 for i, j ∈ Σ.This setting can be used as a probe of a deconfined (con-

Disorder-free localization
Schwinger effect

d)
Signatures of finite-T phase transition  3) with experimentally realistic parameters in a system with coordination number z = 3 (see inset).In panel a) we observe disorder-free localization by initializing the system in a gauge-invariant (blue curve) and gauge-noninvariant (red curve) initial state with two matter excitations localized in subsystem A and calculating the time-averaged imbalance between subsystem A and B as shown.In panel b), we probe the Schwinger effect by quenching the vacuum state with the microscopic model for different experimentally relevant parameters: matter detuning ∆m (chemical potential) and link detuning ∆ l (electric field).We find lines of resonance, where the production of matter excitations out of the vacuum is large.In panel c) we plot the average U (1) matter density (blue curve) obtained from DMRG calculations on a ladder with J < 0. We can qualitatively understand the sharp decay of matter as a transition into the vacuum phase as discussed in Fig. 2a.Additionally, a kink in the plaquette expectation value (red curve) signals a phase transition.In panel d), we use two fluctuating test charges to probe a temperature-induced deconfinement transition in a classical limit of our effective model using Monte Carlo simulations.Both in the percolation strength (red curve) and the Euclidean distance of two matter excitations (blue curve), we find that above a certain temperature T /h the system undergoes a percolation transition.fined) phase, in which the Z 2 matter is free (bound) [53].
To determine the classical equilibrium state, we note the following: 1) Due to the electric field term h in the Hamiltonian, a string of flipped electric fields τ x i,j = −1 costs an energy 2h• , where is the length of the string.2) Gauss's law enforces that at least one string is connected to each matter excitation.
Hence, in the classical ground state the two matter excitations form a mesonic bound state on nearest neighbor lattice sites.Therefore, the matter excitations are confined by a linear string potential.In the co-moving frame of one matter excitation, this model can approximately be described as a particle in a linear confining potential.
At non-zero temperature T > 0, the entropy contribution to the free energy F = E − T S must also be considered.Even though the electric field term h yields an approximately linear string tension, the two charges can separate infinitely in thermal equilibrium provided that E( ) < T log(N ) for → ∞, where log(N ) = S denotes the entropy S of all the string states N with length (setting k B = 1) and E( ) is their typical energy [53].This happens beyond a critical temperature T > T c , when a percolating net of Z 2 electric strings forms.
At the critical temperature T c we anticipate a thermal deconfinement transition, where matter excitations become free Z 2 charges (bound mesons) for T > T c (T < T c ).To study this transition we use the percolation strength -a measure for the spatial extend of a global string net (see Methods) -as an order parameter for the deconfined phase.For experimentally realistic parameters, we find a sharp transition for both the percolation strength and Euclidean distance between two mat-ter excitations around (T /h) c ≈ 2 as shown in Fig. 3d.Although our classical simulation neglects quantum fluctuations, we expect that the revealed finite-temperature deconfinement transition is qualitatively captured.
For a finite density of matter excitations in the system, the Euclidean distance is not a reasonable measure anymore.However, we speculate that a percolation transition might be related to (de)confinement at finite densities.How this transition is related to the quantum deconfinement transition at T = 0 [54,55], driven by quantum fluctuations, will be subject of our future research.Hence, experimentally exploring this transition not only in the classical case, but also in the presence of quantum fluctuations could give insights in the mechanism of charge (de)confinement.

III. CONCLUSION
We introduced an experimentally feasible protection scheme for Z 2 mLGTs and QDMs in (2 + 1)D based on two-body interactions, where the Z 2 gauge structure emerges from well-defined subspaces at high and low energy, respectively.The scheme not only allows reliable quantum simulation of gauge theories but provides an accessible approach to engineer gauge-invariant Hamiltonians.We derived an effective Z 2 mLGT, Eq. ( 4), and QDM, Eq. ( 5), and discussed some of their rich physics.In particular, we suggested several experimental probes, for which we provide numerical analysis using ED of the experimentally relevant microscopic model (3) as well as DMRG and Monte Carlo simulations of the effective models.Experimentally, we anticipate that significantly larger systems are accessible.
Our proposed scheme is not only suitable and realistic to be implemented in Rydberg atom arrays, see Eq. ( 3), but it is also of high interest for future theoretical and numerical studies.Hard-core bosonic matter coupled to Z 2 gauge fields in (2 + 1)D plays a role in theoretical models, e.g. in the context of high-Tc superconductivity [38].While certain limits such as the fine-tuned, classical limit studied by Fradkin and Shenker [6] or coupling to fermionic matter [14,15] are well-understood, surprisingly little is known about the physics of our proposed model.What are the implications of anisotropic hopping and pairing t = ∆ 1 or anomalous pairing terms ∆ 2 , i.e. when the classical mapping fails?How can (de)confinement in the presence of dynamical matter be captured?Is disorder-free localization a mechanism for ergodicity breaking in ( 1.The physical Hilbert space of gauge theories is highly constrained and given by the gauge constraint Ĝj |ψ physical = g j |ψ physical .In contrast the Hilbert space of the experimental setup is larger and also contains unphysical states |ψ unphysical , which do not satisfy Gauss's law.Therefore, the dynamics of the system is fragile in the presence of experimental errors which couple physical and unphysical states.However, it has been shown that this can be reliably overcome by energetically gapping the physical from unphysical states using stabilizer/protection terms in the Hamiltonian [27,28].These strong stabilizer terms can be understood as strong projectors onto its energy eigenspaces, which are chosen to be the physical subsectors of a Z 2 gauge theory in our case; hence the effective dynamics is constraint to quantum Zeno subspaces [56].Note that here the quantum Zeno effect is fully determined by a unitary timeevolution and not driven by dissipation, in agreement with the original effect [56]. The obvious choice of such a protection term is the symmetry generator, Eq. ( 1).However, this requires strong and hence unfeasible multi-body interactions.In contrast, the LPG term Ŵj , Eq. ( 2), only contains two-and one-body terms and is engineered such that an energy gap between the physical and unphysical states is introduced under the reasonable condition that only one (target) gauge sector is protected.In particular, the LPG term in the 2D honeycomb lattice fulfills the condition where V is the strength of the LPG term.The spectrum of Ŵj for the gauge choice g j = +1 is illustrated in Fig. 1c.
2. To study gauge theories, a Z 2 -invariant Hamiltonian has to be engineered first, e.g. the Hamiltonian (4) discussed in the main text.In our scheme we exploit the LPG term with its large gap between energy sectors to construct an effective Hamiltonian perturbatively as explained in Supplementary note 2.
To faithfully stabilize large systems for -in principleinfinitely long times, we want to discuss the stabilization of high-energy sectors by considering undesired instabilities/resonances in the spectrum V j Ŵj .The eigenvalues of V Ŵj are w j = (0, V, 4V ) and we want to protect a sector with intermediate energies.If the interaction strength V is equally strong at each vertex gaugesymmetry breaking can occur.For example, by exciting vertex j 0 and simultaneously de-exciting three vertices j 1 , j 2 and j 3 .This process has a net energy difference of ∆E = +3V − 3 • V = 0 and the resonance between the two states can lead to an instability towards unphysical states, hence gauge-symmetry breaking (Supplementary note 2 G).
Therefore, the LPG method without disorder cannot energetically protect against some states that break Gauss's law on four vertices.An efficient way to stabilize the gauge theory even against such scenarios is to introduce disorder in the coupling strengths by Ŵ = j V j Ŵj with V j = V + δV j .The couplings δV j are random and form a so-called compliant sequence [27,29].In 1D systems, this has been shown to faithfully protect Z 2 LGTs also for extremely long times, see Ref. [29] for a detailed discussion of (non)compliant sequences.Moreover, we note that for small system sizes and experimentally relevant timescales even noncompliant sequences such as the simple choice V j = V ∀j lead to only small errors (Supplementary note 2 G).
For our (2+1)D model, we illustrate the effect of disordered protection terms in Fig. 4, which shows that only the gauge non-invariant states are shifted out of resonance.Moreover, we propose to use weak disorder such that the overall perturbative couplings remain unchanged in leading order.We emphasize that the disorder scheme does not require any additional experimental capabilities but only arbitrary control over the geometry as well as local detuning patterns.Even more, an experimental realization will always encounter slight disorder, i.e. the gauge non-invariant processes might already be sufficiently suppressed in experiment.
We further note that the example above, where Gauss's law is violated on four vertices, yields gauge-breaking terms in third-order perturbation theory.Ensuring that none of the protection terms V j have gauge-breaking resonances within such a nearest-neighbour cluster, these terms can be suppressed.However, now it remains space for fifth-order breaking terms on next-nearest neighbour with Ω = 0 and plot all eigenstates around energy E = 4V .Green (red) dots are states that fulfil (break) Gauss's law as illustrated with two examples in the inset of panel a).Without disorder, i.e.V j = V for all j, the physical and unphysical states are on resonance.In panel b), we show the effect of disordered protection terms V j = V + δV j , which only shifts the unphysical states out of resonance and hence fully stabilizes the gauge theory.We note that even without disorder, the emergent gauge structure is remarkably robust (Supplementary note 2 G).
vertices.Hence, the non-resonance condition is now desired on a larger cluster and so forth.Therefore, systematically choosing the disorder potentials can suppress gauge-breaking terms to arbitrary finite order and stabilize gauge invariance up to exponential times.Its fate in the thermodynamic limit, however, is an open question beyond the scope of this study.
Percolating strings from classical Monte Carlo.-Thefinite temperature percolation transition in Fig. 3d is obtained from classical Monte Carlo simulations on the honeycomb lattice with matter and link variables.In this section, we discuss the percolation strength order parameter [57] and details of the numerical simulations in more detail.
The classical model we consider is motivated by the microscopic Hamiltonian (3) and its effective model (4) -in particular we used the precise effective model as derived in Eq. (S13) of Supplementary note 2 for Ω/V = 1/8, ∆ m = V /2 and ∆ l /V ≈ 0.044.For elevated temperatures T V , we expect that classical fluctuations dominate in the system while the Gauss's law constraint is still satisfied due to the LPG protection.Therefore, we neglect quantum fluctuations and set t = ∆ 1 = ∆ 2 = J = 0. Hence, the resulting matter-excitation conserving Hamiltonian is purely classical and a configuration is fully determined by the distribution of matter and electric field lines under the Gauss's law constraint, i.e. {(n j , τ x i,j ) | (−1) n j = g j i: i,j τ x i,j ∀j} and we consider the sector with g j = +1 ∀j.
From the numerical Monte Carlo simulation, we want to quantify the features discussed in the main text: 1) string net formation and 2) bound versus free matter excitations.To this end, we define the percolation strength as the number of strings in the largest percolating cluster of Z 2 electric strings, normalized to the system size.Furthermore, we consider the Euclidean distance between two matter excitation and show that an abrupt change of behaviour in this quantity indicates the disappearance of the bound state.
The Monte Carlo simulations are performed on a 35×35 honeycomb lattice (in units of lattice spacing) using classical Metropolis-Hastings sampling (Supplementary note 9).Further analysis of the obtained samples allows to extract the number of strings in the largest percolating cluster to calculate the percolation strength.As shown in Fig. 3d, we find that for low temperatures T the percolation strength vanishes.At a critical temperature (T /h) c ≈ 2, the percolation strength abruptly increases, i.e. the string net percolates.Moreover, at the same critical temperature (T /h) c ≈ 2 the Euclidean distance shows a drastic change of behavior and saturates at about 30 for high temperatures.This saturation can be explained by the finite system size.

IV. DATA AVAILABILITY
The datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request.

V. CODE AVAILABILITY
The data analysed in the current study has been obtained using the open-source tenpy package; this DMRG code is available via GitHub at https://github.com/tenpy/tenpy and the documentation can be found at https://tenpy.github.io/#.The code used in the exact diagonalization and Monte Carlo studies are available from the corresponding author on reasonable request.

VI. AUTHOR CONTRIBUTIONS
LH, JCH and FG devised the initial concept.LH proposed the idea for the two-dimensional model, worked out the main analytical calculations and performed the exact diagonalization studies.LH, AB and FG proposed the experimental scheme.SL performed the Monte Carlo simulations.AB conducted the DMRG calculations.All authors contributed substantially to the analysis of the theoretical results and writing of the manuscript.

VII. COMPETING INTERESTS
Authors declare that they have no competing interests.S1.Spectrum of the local pseudogenerator for g j = −1.We show the LPG term for Z2 mLGTs with g j = −1.

I. LOCAL PSEUDOGENERATOR ON THE 2D HONEYCOMB LATTICE
In the following, we discuss local pseudogenerators (LPG) for arbitrary Z 2 mLGT gauge sectors as well as for QDMs.
A. LPG for Z2 mLGTs and gj = −1 The LPG term in the main text, Eq. ( 2), can be easily generalized to protect any of the two g j = ±1 sectors by choosing: The case g j = +1 is shown and discussed in the main text, Fig. 1c, while the case g j = −1 is illustrated in Fig. S1.

B. Quantum Dimer Models
Rokshar and Kivelson [41] introduced the QDM as a toy model to study short-range resonating valence bond (RVB) states on the square lattice.Their model has two phases: a columnar and a staggered phase.At the phase transition, the so-called Rokshar-Kivelson point, the model becomes exactly solvable and has deconfined monomer excitations.The experimental challenge is to impose the hard-core dimer constraint and to gap out monomers -the fundamental, fractionalized excitations of the system.Here, the LPG term overcomes both challenges.
As shown in Fig. 1c the ground-state manifold of the LPG term allows for six different configurations per vertex j.The subsector with n j = 0 (n j = 1) should be called QDM 1 (QDM 2 ) and we want the two subsectors to be decoupled.This can be exactly fulfilled by entirely eliminating the local matter degrees-of-freedom, i.e. experimentally only the gauge sector Perturbative derivation of the effective Z2 mLGT Hamiltonian.From the perspective of Z2 mLGTs, the LPG protection term energetically splits a target gauge sector (yellow) from other sectors (orange, green) as shown in panel a).
A gauge-noninvariant perturbation Ĥdrive with strength Ωm, Ω l V leads to virtual processes to unphysical sectors of the Hilbert space, which can be treated in perturbation theory and which ultimately yield the effective Hamiltonian Ĥeff Z 2 in the main text.In panel b), we illustrate an example for a second-order (left) and third-order (right) process.By using projection operators on the initial (final) state, Πinit ( Πfinal ), the operator form of Ĥeff Z 2 can be determined.Panel c) introduces the notation for sites and links on the 2D honeycomb lattice with lattice vectors shown in light red.
link atoms on the Kagome lattice are implemented, see Fig. 1a.Hence, the LPG term for the two subsectors read In contrast to the Z 2 mLGT, we note that the QDM 1 (QDM 2 ) subspaces are now the lowest-energy eigenspaces of the LPG term.Therefore, any state violating the hard-core dimer constraint has a larger energy, which qualifies the LPG term as a full-protection scheme [27] for QDMs.

II. DERIVATION OF THE EFFECTIVE Z2 mLGT HAMILTONIAN
In this section, we explain the derivation of the effective Hamiltonian (4) in terms of a Schrieffer-Wolff transformation [58].The derivation of the effective QDM is discussed in SI IV.Starting point is the experimentally motivated microscopic Hamiltonian (4), where Ĥ0 = ĤLPG + Ĥdetuning is the unperturbed Hamiltonian and Ĥdrive is a small perturbation [59], i.e.V Ω m , Ω l , see Fig. S2a.Note that the perturbation is a gauge-symmetry breaking term, [ Ĥdrive , Ĝj ] = 0 ∀j.However a state prepared in the physical subspace, g j = +1 ∀j, will only virtually occupy unphysical states under Ĥmic because of the large energy gap V between the sectors in the limit of weak driving, Ω m , Ω l V .Hamiltonian Ĥ0 is diagonal in the matter density and electric field basis and hence the unperturbed eigenstates are product states |α = j |n j i,j |τ x i,j .Since Ĥdrive only contains off-diagonal elements, there are no first-order contributions, α| Ĥdrive |α = 0.The derivation of the second-and third-order terms are explained in the following together with an explicit example, see Fig. S2b.We note that the second-and third-order contributions require to calculate 16 + 32 + 3 • 2 • 16 = 144 amplitudes.
The second-order terms are given by where |α (|β ) are the initial (final) state and |δ are virtual states.Because Ĥdrive has only off-diagonal elements, it always couples to states outside the physical energy sector and hence in second-order the initial and final state coincide, |α = |β , in order to remain within the same energy subspace.In Fig. S2b (left) we show one example process for the parameters While the amplitude of the process can be calculated using Eq.(S8), the operator form can be expressed in terms of projectors, which for the example in Fig. S2b left is given by Πinit and hence only diagonal terms appear in second-order perturbation theory.Here, we have used the notation introduced in Fig. S2c: j ν = (j x , j y , ν) corresponds to an explicit site on the honeycomb lattice with two-basis unit cell and j 2 , j 1 + x or j 2 − y, j 1 describe links, where x and y are the unit vectors and ν = 1, 2 is an intracell index.The notation i, j should still be used when all links are addressed.
In third-order perturbation theory, coupling between different states, |α = |β , occurs, which yields dynamical hopping and pairing terms.The coupling elements in the effective Hamiltonian can be calculated by evaluating where the sum runs over two virtual states |δ , |δ .As shown in the example in Fig. S2b (right), we can write down the operator corresponding to the coupling (S11) by projecting onto the initial state Πinit , then acting with an operator coupling to the final state followed by a projection on the latter by Πfinal .In our example, the projector reads Πfinal = |β β| Executing the above steps to all states in the target energy subspace yield the effective Hamiltonian (4).Note that the plaquette terms are not appearing directly in third-order perturbation but would require to go to sixthorder perturbation theory.Hence, we discuss them separately in SI II B and II D. First, we want to give an explicit expression of the effective Hamiltonian up to third-order perturbation theory and distinguish the cases with and without global U (1) symmetry in SI II A and II C, respectively.
To enforce conservation of matter excitations, we introduce an additional energy gap between different particle number sectors by choosing ∆ l , Ω m , Ω l .This strong chemical potential term suppresses creation and annihilation of matter excitations induced by Ĥdrive .
The effective model for U (1) matter coupled to a Z 2 gauge field in the sector g j = +1 ∀j is given by  U (1) FIG. S3.Effective couplings.We plot the effective couplings for the Z2 mLGT as derived in perturbation theory up to third order, see also Tab.SII.In panel a) we show the U (1) matter case for two different choices of matter detuning ∆m = ±V /2.We do not plot the effective chemical potential term because it only contributes as a constant term in Hamiltonian (S13).In panel b), we show the couplings for the effective quantum-Z2 matter Hamiltonian, Eq. (S23).Note that small detunings ∆m and ∆ l can fully tune the chemical potential µ and electric field h without affecting the other couplings in perturbation theory.
The operator form and its corresponding coupling amplitudes for the second-and third order processes can be found in the fourth column of Tab.SII and are plotted in Fig. S3a.The plaquette interaction ∝ J is a sixth-order perturbative term, which is discussed separately in SI II B. Note that Gauss's law, Ĝj = +1 has been used to simplify, collect and eliminate higher-order terms.
The terms ∝ M , ∝ χ 1 , ∝ χ 2 and ∝ χ 3 are (nearest neighbor density-density), (next-nearest neighbor density-electric field), (next-nearest neighbor electric field-electric field) and (nearest neighbor density-electric field) interactions, respectively.In the main text, Eq. ( 4), we treat these terms on mean-field level in the electric field τ x i,j and matter density nj , which is well-defined since both quantities are gauge invariant.To be explicit, we perform for example a mean-field decoupling of M i,j ni nj → M ni j nj , which simplifies the effective Hamiltonian.We want to perform an order of magnitude estimation of the plaquette interactions J in Eq. ( 4).The goal is to find the matrix elements corresponding to plaquette interactions e.g.J eff ( +h.c.), which we derive by a Schrieffer-Wolff transformation from Eq. ( 3) with Ω m = Ω l = Ω, see below.In general, the effective coupling strengths J eff n = J eff ({n j }, {τ x i,j }) depend on the configuration of matter and electric fields, yielding n max = 416 independent J eff n after taking the 6-fold symmetry of the plaquette and Gauss's law into account.
Hence, the effective plaquette interaction Hamiltonian takes the form The different coupling elements J eff n are calculated in degenerate perturbation theory (see below) and plotted in Fig. S4a for several driving strength Ω/V .Because the plaquette interaction involves six links, we expect that the effective couplings scale as (Ω/V ) 6 .
Here, we want to estimate and simplify the plaquette interaction Eq. (S14) by averaging over all configurations, i.e. we consider J eff U (1) = 1/n max n J eff n .To this end, we extract J eff U (1) = J eff U (1) (Ω/V, ∆ m , V ) and perform a fit with the expected scaling function f (Ω/V, ∆ m , V ) = α(∆ m , V ) (Ω/V ) 6 as shown in Fig. S4b.By extracting the fit parameter α(∆ m , V ) for ∆ m = V /2, we can estimate the strength of the plaquette terms as FIG. S4.Estimation of plaquette terms (U (1) matter).The effective plaquette interaction derived by a Schrieffer-Wolff transformation depends on the matter and electric field configuration within each plaquette.In panel a) we plot the absolute value of the coupling strength for all 416 different configurations for various driving strength Ω/V and ∆m = V /2 (from dark to bright shade: V /Ω = 15, 18, 20, 25, 30).In panel b), we have averages (with the correct sign taken into account) over the 416 couplings elements for each driving strength Ω/V , which we plot on a log-log scale.The linear behaviour indicates a power-law behavior and we fit the expected sixth-order perturbation scaling.From the fit we can extract the prefactor, which yields the effective coupling J eff U (1) (∆m = V /2), Eq. (S16).
Let us now discuss the detailed derivation of Eq. (S14) in terms of a Schrieffer-Wolff transformation.The microscopic model is given by Hamiltonian (S4), where V, |∆ m | ∆ l , Ω m , Ω l and Ω m = Ω l = Ω.Further we set ∆ l = 0.The drive Ω couples the Z 2 mLGT sector to the gapped, virtual energy sectors defined by the LPG term.Since we expect the effective plaquette terms to arise in sixth-order perturbation theory, we also need to consider couplings of Ω within the highly-degenerate virtual energy sectors.Hence, it is required to apply degenerate perturbation theory and to diagonalize all energy sectors with respect to the perturbation Ω first to lift the degeneracies and afterward perform standard perturbation theory.
To gain an intuitive understanding, we want to consider the following path in the perturbative calculation: for instance we start from a state with no matter excitations and all links in the τ x i,j = +1 configuration.Then, the drive flips one link, which costs an energy of −2V because Gauss's law is violated on two vertices.Now, we can consecutively flip all links in clockwise direction.Since these processes at the same time break and restore Gauss's law at different vertices, they are all degenerate and the denominators of the perturbative expansion vanish.To circumvent this non-physical divergence, we first need to diagonalize the degenerate subspaces, which renormalizes all couplings and energy gaps.
The system is perturbed by a weak drive Ĥdrive , Eq. (S7), and diagonalization of the degenerate subspaces yields the transformed Hamiltonian Ĥmic = Û † Ĥmic Û that is diagonal within the energy blocks but couples states from different blocks in a non-trivial way.The off-diagonal terms in Ĥmic now become the perturbation Ĥdrive in the new basis.Note that the states have also transformed and should be denoted by |α = Û |α in the new basis.
Since we have access to the full one-plaquette spectrum, we can now explicitly construct the unitary operator Ŝ of the Schrieffer-Wolff transformation by calculating the matrix elements β| Ŝ|α = β| Ĥdrive |α where Ĥdrive = Û † Ĥdrive Û and E α, E β are the unperturbed energies in the transformed basis.Because we completely diagonalized the degenerate subspace, divergences of the denominator only appear for uncoupled states, i.e. the nominator vanishes, for which we define the matrix element of Ŝ to be zero.In the Schrieffer-Wolff formalism we can now write down a well-defined expansion in Ω/V : Note that in the transformed basis the energy denominator in Eq. (S17) can depend on V and Ω.Since we require Ω V , we can expand the expressions and find in leading order sixth-order contributions for any 2 ≤ n ≤ 6.
In Fig. S4a, we plot the strength of all non-zero plaquette flip matrix elements in the gauge sector g j = +1 for different driving strength Ω/V .Note that the couplings can be positive and negative while we only plot their absolute value; In the average J eff U (1) their signs are properly included, however.
C. Quantum-Z2 matter: V ∆m, ∆ l , Ωm, Ω l In this section, we discuss the derivation of the effective Hamiltonian ( 4) with quantum-Z 2 matter coupled to a Z 2 gauge field.In contrast to SI II A, we do not enforce a global U (1) symmetry for the matter but otherwise the derivation is completely analogous.This leads to the additional pairing terms ∆ 1 , ∆ 2 in Eq. ( 4).The effective model we find is invariant under the local transformation However, the 2D quantum Hamiltonian cannot be mapped exactly on a classical 3D Ising LGT [6], which is origin of the term "quantum-Z 2 mLGT".
In the gauge sector g j = +1 ∀j, the effective model reads Ĥeff, (3) The operator form and its corresponding second-and third-order coupling amplitudes can be found in the fifth column of Tab.SII and are plotted in Fig. S3b, while the discussion of the sixth-order plaquette terms is dedicated to SI II D.
Compared to (S13), we now find pairing terms ∆ 1 and ∆ 2 , which also appear in Fradkin & Shenker's Ising Z 2 mLGT in a similar fashion.As explained in SI II C, the terms ∝ M , ∝ χ 1 , ∝ χ 2 and ∝ χ 3 can be treated on mean-field level yielding the effective model ( 4) discussed in the main text.
In particular, the electric field term −h i,j τ x i,j can be fine-tuned by changing the detuning ∆ l , which in the limit ∆ l V does not alter the results obtained from perturbation theory.On mean-field level, this allows to tune the expectation value to τ x i,j = −1/2.Then, the effective coupling renormalizes to ∆1 = ∆ 1 − τ x i,j ∆ 2 = ∆ 1 −∆ 2 /2 = t.At this particular point, we retrieve the (2 + 1)D model studied by Fradkin & Shenker [6], where it is known to map on a classical 3D Z 2 mLGT with continuous phase transitions in the Ising universality class.Note that our model is defined on the honeycomb and not square lattice.For a detailed discussion of the duality between a Z 2 mLGT on a honeycomb and triangular lattice, we refer to the Supplementary Information of Ref. [32].Because of this duality, the results obtained in Ref. [6] should be still valid, however the phase diagram might not be symmetric across the diagonal as illustrated for simplicity in Fig. 2b.The effective plaquette interaction derived by a Schrieffer-Wolff transformation depends on the matter and electric field configuration within each plaquette.In panel a) we plot the absolute value of the coupling strength for all 416 different configurations for various driving strength Ω/V (from dark to bright shade: V /Ω = 8, 10,12,14,20,30).We find a plateau with strongest coupling for certain staggered and polarized electric field configurations as shown in the inset.Note the sign of J stag eff differs from J pol eff .In panel b) the absolute value of the effective coupling of the staggered, polarized and averaged configurations versus the driving strength V /Ω on a log-log scale is shown.The linear behaviour indicates a power-law behaviour and we fit the expected sixth-order perturbation scaling.From the fit we can extract the prefactor, which yields the effective couplings Eqs.(S26), (S29) and (S31).D. Plaquette terms for quantum-Z2 matter: V ∆m, ∆ l , Ωm, Ω l Similar to the case with U (1) matter discussed in SI II B, we want to estimate the strength of the plaquette terms J eff

Z2
in the quantum-Z 2 matter model.We perform a Schrieffer-Wolff transformation with V ∆ m , ∆ l , Ω m , Ω l and Ω m = Ω l = Ω and ∆ l = 0.In Fig. S5a, we plot the extracted coupling matrix elements between flippable plaquettes.We find that there is a plateau with 14 distinct coupling elements, which are an order of magnitude larger than the remaining couplings.As indicated in Fig. S5b, these couplings correspond to 1) a staggered matter and electric field configuration with J eff stag and to 2) configurations with a polarized electric field J eff pol , where all links are either τ x i,j = +1 or τ x i,j = −1.Note that these coupling elements might give rise to additional phases and we want to include them in the discussion of the plaquette terms here.However, averaging over all coupling elements as in SI II B should give a useful estimation of the overall strength J eff Z2 of the plaquette terms.As discussed in SI II B, we can extract the strength of the plaquette interaction by fitting the coupling elements for different driving strengths Ω/V .We want to examine the three cases 1) staggered, 2) polarized and 3) averaged as shown in Fig. S5c: (1) For plaquettes with a staggered matter and electric field, we find Ĥeff stag,Z2 = −J eff stag P stag (S24) . Dynamics of microscopic model for different parameters.We initialize the system in a gauge-invariant state Ĝj |ψ init = init ∀j and time-evolve for time T under the microscopic Hamiltonian (3) using exact diagonalization in a small-scale system, see inset in panel b)-d).We find strong dynamics of matter, link and plaquette degrees-of-freedom while the error in Gauss's law is small and remains constant for long, experimentally relevant timescales.The results are discussed in detail in SI II E. Note that the labels in the plots refer to the color scheme but not necessarily to the linestyle and we plot e.g. the expectation values of all four matter site with the same color but all four curves with a different linestyle (some curves overlap).
(2) For plaquettes with polarized electric field, we find Ĥeff pol,Z2 = −J eff pol J eff pol ≈ −5.81 (Ω/V )6 (S29) (3) By averaging over all couplings (as in SI II B), we find Small-scale exact diagonalization study of the microscopic Hamiltonian In this section, we present results from time-evolution studies obtained by exact diagonalization of the full microscopic Hamiltonian (3) in a minimal model with coordination number z = 3, i.e. four matter sites and six links as shown in the inset of Fig. S6b-d.While this model has a tetrahedron structure and triangular plaquettes, it is different from the honeycomb lattice.However, because the model has coordination number z = 3 -similar to the honeycomb lattice -the physics of the LPG protection should be correctly modeled in this numerically feasible 2D system.
We demonstrate that Gauss's law is indeed very well conserved, Ĝj ≈ +1, even for relatively strong drive Ω/V (we set Ω m = Ω l = Ω throughout this section).Moreover, the matter and link degrees-of-freedom show dynamics on the expected timescales.The results are summarized in Fig. S6 and we want to elaborate on the different cases here: • Fig. S6a: We plot the expectation value of Gauss's law after time-evolving different initial states and different parameters.If not specified otherwise, the initial state contains two localized matter excitations and fulfills Gauss's law, Ĝj |ψ init = +|ψ init .While Ĝj has an initial fast drop, the gauge-symmetry violation equilibrates around a constant value determined by the driving strength Ω/V .For Ω/V = 0.125 (Ω/V = 0.2), the violation is about 5% (15%).
• Fig. S6b: We consider U (1) matter, i.e. we have strong detuning/chemical potential ∆ m = ±V /2 and plot the expectation values of the matter density nj , the electric field τ x i,j and plaquette terms i,j ∈P τ z i,j as well as the total number of matter excitations and its variance.Note that the total number of matter excitations only fluctuates marginally as anticipated for U (1) matter.Calculating the effective hopping from Tab. SII, we expect oscillations with a timescale T V = 2π × 2520/13 ≈ 1220, which matches the timescales in Fig. S6b approximately.
• Fig. S6c: Next, we consider quantum-Z 2 matter, where pairs of matter excitations can be created and annihilated.Since the initial state has already two matter excitations (and two holes), the pair creation dynamics is not as heavy as in Fig. S6d, we start from vacuum.Because of the interplay between hopping and pairing, it is not straightforward to read off timescales from Rabi oscillation-like behaviour.From hopping and pairing, we would expect timescales of approximately T V ≈ 2400, respectively.However, we find an emergent timescale in this finite size model of about T V = 1000.Note that in Hamiltonian (S23), we have (anomalous) pairing terms, which also influence the propagation of matter excitations.
• Fig. S6d: Here, we initialize the system in the vacuum state and otherwise time-evolve with the same parameters as in Fig. S6c.We find strong particle number fluctuations due to the creation of matter excitations.The expected timescale T V ≈ 800 (on mean-field level) is in agreement with the overall timescale of oscillations we observe in the system.

F. Microscopic versus effective model
In this section, we want to confirm the effective model by comparing the energy of the microscopic and effective model as a function of V and Ω (Ω m = Ω l = Ω) and show that for both quantum-Z 2 matter, Eq. (S23), and U (1) matter, Eq. (S13), the spectra converge in the limit V /Ω → ∞ as expected from perturbation theory.To this end, we perform exact diagonalization calculations of a minimal system (four matter sites and six links) as in SI II E. We set the LPG protection strength to V = 1 and vary the drive Ω/V in the microscopic model (3) or correspondingly use the derived effective couplings, see Tab.SII.Moreover, we set the link detuning ∆ l = 0 and choose ∆ m = 0 (∆ m /V = 0.5) in the quantum-Z 2 matter (U (1) matter) case.
As a first step, we need to identify the correct target sector of the microscopic model since this has no exact Z 2 gauge symmetry or global U (1) symmetry.Therefore, we diagonalize the microscopic Hamiltonian (3) and calculate the expectation value of the symmetry generators g = j Ĝj for each eigenvector as shown in Fig. S7a and b.Because we choose the LPG term to protect the target sector g j = +1 ∀j, we want g = 4 in our numerical study.Additionally, for the case of U (1) matter, we need to select a matter excitation sector by evaluating n = j nj and we choose n = 2 in the following discussion.
Fig. S7 illustrates that the target gauge sectors for both cases, quantum-Z 2 and U (1) matter, form well-separated subspaces.We want to emphasize the efficiency of our proposed LPG protection scheme: As discussed in SI I, there are instabilities because we work in a high-energy sector of the LPG term.These instabilities are resonant processes, where Gauss's law is violated in a way that on three vertices the LPG term lowers the energy while on one vertex the energy is increased.If the instabilities would play a dominant role in the effective dynamics, we would expect no  3) to the effective Hamiltonian with quantum-Z2 matter, Eq. (S23), and U (1) matter, Eq. (S13).For the first (latter) case we choose ∆ l = ∆m = 0 (∆ l = 0, ∆m = 0.5V ).In panel a) and b), we calculate the expectation value g = j Ĝj for each eigenvector of the microscopic model and our target sector is g = 4.In panel b), we additionally require the number of matter excitations to be conserved.For Ω = 0.1V , we find that both the local and global symmetry emerge in the microscopic model.In panel c) and d), we now consider the eigenenergies in the target sector for different driving strength Ω/V for the microscopic (blue) and effective (red) model.We find agreement of both spectra which supports the validity of our perturbative approach discussed in the main text and SI II.
well-defined gauge sectors but a hybridization of all sectors which would broaden the clusters we find in Fig. S7a and  b.
As a next step, we show that the spectra in the target sectors of the effective and microscopic model converge as V /Ω → ∞.To this end, we diagonalize the microscopic model ( 3) and the effective model for different V /Ω, which yields eigenenergies E n mic (V, Ω) and E n eff (V, Ω).To compare the spectrum at different driving strengths, we normalize the eigenenergies by the ground-state energy E 0 eff (V, Ω) of the corresponding effective model at each point V /Ω.In Fig. S7c and d, we plot the spectrum for the quantum-Z 2 and U (1) matter case as described above.We find that by using the derived effective couplings in Tab.SII the effective models, Eqs.(S23) and (S13), very well describe the microscopic models justifying our perturbative analysis.Note that we did not take the above derived plaquette interactions into account here, because the small-scale system we use in the exact diagonalization study has plaquettes with three edges instead of six edges on a honeycomb lattice.
3rd order process FIG.S8.Third-order gauge-breaking process.The state illustrated on the left (right) fulfills (breaks) Gauss's law at every vertex.The two states are coupled resonantly via a third-order process.

G. Gauge non-invariant processes
So far, we have only considered resonant processes that conserve Gauss's law.However, as discussed in the Methods section and Fig. 4, the LPG method without disorder suffers from unwanted resonances with a few unphysical states.Here, we want to discuss the effect of such resonances with respect to the numerical results from section SI II F.
As shown in Fig. S8, it is possible to raise the energy by +3V on one vertex and at the same time lowering the energy by −V on three neighbouring vertices.This process is resonantly coupled to the physical states in a thirdorder process if and only if the four vertices are arranged in a "star" geometry, see Fig. S8.Otherwise, this type of resonance only appears in fifth-order perturbation theory.In the following, we argue that these processes do not alter the emergent gauge structure such that the effective model ( 4) is valid.
First, we note that the above processes can be entirely suppressed by applying weakly disordered protection terms, V → V j + δV j and δV j V .The disorder only shifts the gauge non-invariant states out of resonance, see Fig. 4, and its efficiency in (1 + 1)D has been demonstrated numerically [29].
We point out that our minimal model simulation in the Mercedes star is able to capture the third-order terms described in Fig. S8.Hence, the system is susceptible to resonant non-gauge invariant terms that potentially could lead to a complete breakdown of gauge invariance.However, reconsideration of the numerical results presented in Fig. S7a) and b) show a well-defined target sector g = +4, which is energetically in resonance with a well-defined g = −4 sector.I.e. the two sectors only very weakly hybridize and the eigenstates of the microscopic Hamiltonian are almost exact Gauss's law eigenstates.In contrast to the time-evolution of as shown in Fig. S6a), which depends on the choice of the initial state, the plotted spectrum in Fig. S7a) and b) is a very sensitive probe to validate the emergent gauge Moreover, this holds true for even stronger drivings Ω/V = 1/5, we use in further numerical simulations below.
Both thermalization dynamics well as hybridization between the physical and unphysical sectors is highly suppressed despite comparable Hilbert space dimensions.This robust gauge structure further suggest an additional mechanism that stabilizes the gauge sectors such as Hilbert space fragmentation, which should be investigated in future studies.
To summarize, we identify potential third-order processes and we an easily implementable disorder-based protocol such that gauge invariance remains fully intact.Furthermore, we observe from our numerical results that resonant physical and unphysical sectors show only very weak mutual coupling giving almost perfect gauge invariance even without disorder.

III. EFFECTIVE MESON MODEL
For U (1) matter coupled to a Z 2 gauge field, we predict the existence of a meson condensate phase, see Fig. 2a.Here, we want to derive an effective meson model, which captures the condensate phase.
In the limit J/h, t/h → 0 and dilute U (1) matter in the ground state, electric field strings are minimized under the constraints imposed by Gauss's law, i.e. number of links with τ x i,j = −1 is minimized.To fulfill Gauss's law g j = +1 ∀j matter excitations are bound into pairs connected by an electric field string, see Fig. S9a.Gaugeinvariant hopping of matter excitations prolongs the string and thus kinetic energy t competes with the string tension h.Since h t it is unfavourable for single matter excitations to be mobile, which justifies to describe the constituents as tightly bound mesons.
Nevertheless, the mesons can gain kinetic energy in two distinct processes: 1) a second-order hopping process t eff = −t 2 /2h, in which the entire pair moves from one link to a neighboring link as shown in Fig. S9b and 2) plaquette interactions ∝ J induce fluctuations between the two different meson configurations on a plaquette as illustrated in Fig. S9d.Additionally, the mesons gain dispersive energy shifts δ eff = −t 2 /2h if the matter excitation hops back and forth on neighboring sites, see Fig. S9c.However, this process is only allowed for an empty neighboring site and therefore the dispersive energy shift leads to repulsive interactions between mesons.To summarize, we find an effective model of Z 2 neutral, hard-core bosonic mesons bn hopping on the sites of a Kagome lattice, with infinitely strong nearest neighbor (NN) repulsion, finite next-nearest neighbor (NNN) repulsion and plaquette interactions, see Fig.For experimentally relevant parameters, see SI II, we can choose J/t t/h and J/h 1 and thus neglect the plaquette interaction term.In the limit of dilute mesons b †b ≈ 0, we can treat PNN on a mean-field level yielding free hard-core bosons on the Kagome lattice.In the ground state the mesons b condense as indicated in Fig. 2a.
Taking the plaquette interactions into account, i.e.J/t ≈ t/h, phase separation by clustering of mesons has been discussed [15] for spinless fermions on the square lattice.However, the NNN repulsive interaction should suppress clustering and hence a more sophisticated analysis is required.Away from the above discussed limit J/h 1, the meson pairs have some finite extend 2 , which alters the effective model (S32).However, for a sufficiently dilute gas of matter excitations, i.e. 2  1/ b †b , we expect the description of free hard-core bosons to be still valid which we indicate by the finite extend of the meson condensate phase in Fig. 2a.We note that the phase boundary is expected to end at the deconfinement-confinement transition of the vacuum since in the deconfined phase the picture of bound mesonic pairs breaks down.
Furthermore, at sufficiently large filling the interplay between kinetic energy and repulsion on the Kagome lattice might lead to additional, exotic phases of matter.A more detailed phase diagram is beyond the scope of this Article and is a topic for future studies.3) on one plaquette for different driving strengths in the flippable plaquette subspace.The blue dashes show the unperturbed system.The two low-energy states are the two flippable plaquette configurations in the QDM sector.The first (second) excited manifold has two (four) monomers and the two high-energy states contain a maximum number of six monomers.To perform perturbation theory, we first need to diagonalize the degenerate subspaces because the drive couples within the two and four monomer excitation blocks, i.e. the drive can move around the monomer excitations without energy cost, thus making them mobile.The orange and green dots show the spectrum for different driving strength, where the blocks are diagonalized.However, there are still off-diagonal couplings between the blocks.These couplings are the starting point for the second step, which is the actual Schrieffer-Wolff perturbation theory.The inset illustrates the smallest energy gap and its renormalized coupling after the first step.This allows to compare the effective driving strength to the energy gap in order to determine a regime of validity for the perturbation theory.

IV. DERIVATION OF THE EFFECTIVE QUANTUM DIMER HAMILTONIAN
In this section, we want to derive the effective Hamiltonian (5) of the quantum dimer model, from the microscopic model (3).Here, the physical subspace is given by the QDM 1 low-energy subspace of the LPG term, see Fig. 1c.for strong protection V Ω m , Ω l states in the non-physical sectors are only virtually occupied and we can derive the effective model perturbatively, which yields the plaquette terms ∝ J QDM .
Let us first consider the simpler terms ∝ K.These terms are introduced to drive potential quantum phase transitions in J/K.In our proposed scheme, the strong LPG terms arise from nearest-neighbor interactions.However, the Rydberg-Rydberg interactions decay as R −6 , where R is the distance between two atoms in optical tweezers.Hence, there are small but finite next-nearest neighbor interactions between links of the honeycomb lattice (next-nearest neighbors on the Kagome lattice), which give rise to the term ∝ K in Eq. (S33).Now, we want to elaborate on the degenerate perturbation theory to derive the plaquette terms in sixth order.In the QDM subspace there are only two "flippable" configurations that can be coupled by plaquette terms and hence we can rewrite the interaction by where we have used the electric field string ↔ dimer mapping shown in Fig. 1b.The two configurations shown in (S34) now span the low-energy manifold and are the starting point for our perturbation theory.Above the lowenergy manifold, we have three high-energy subspaces with energies 2V , 4V and 6V since excitations (=monomers) can only be created in pairs.The unperturbed subspaces are shown as blue dashes in Fig. S10.We perturb the system by a weak drive Ĥdrive , Eq. (S7), coupling not only states between subsectors but also within the highly-degenerate manifolds with energy 2V and 4V .Hence, we want to apply the same Schrieffer-Wolff formalism with degenerate subsectors as explained in SI II B.
In Fig. S10, we show the full spectrum for different driving strengths Ω/V (Ω m = Ω l = Ω).In general, the validity of a perturbation theory is determined by the coupling strength divided by the energy gap in the unperturbed system.In degenerate perturbation theory, this quantity has to be evaluated after the transformation Û .As shown in the Hence, the detuning on the matter and link atoms along the boundary has to be adjusted.
inset of Fig. S10, the gap between the low-energy manifold and the first excited states becomes Ṽ = 2V − 2 √ 3Ω and the matrix element between the two states is Ω = Ω/ √ 2. Hence, we find and e.g.Ω/ Ṽ = 1/4 for Ω/V ≈ 1/3, which allows to have relatively strong driving strength in the lab frame.
From Ĥmic , we can now calculate the Schrieffer-Wolff transformation as explained in SI II B [see Eqs.(S18)-(S21)].To summarize, by evaluating (S18) we can derive the leading order contribution of the plaquette interaction in the quantum dimer subspace and find which e.g.yields J QDM /V ≈ 0.01 for Ω/V = 1/3.The effective coupling is surprisingly strong despite the small prefactor 1/144 in the perturbative expansion, Eq. (S21).Intuitively, we can understand these strong many-body interactions to be induced by the highly mobile and gapped monomer excitations.

V. EXPERIMENTAL REALIZATION
The scheme we propose in this Article is particularly suitable for Rydberg atom arrays because we require 1) control over the real-space configuration of atoms, 2) strong nearest-neighbor density-density interactions and 3) driving a two-level system.The microscopic Hamiltonian has been introduced in the main text in Eq. ( 3) and contains only nearest-neighbor interactions as well as an on-site drive, which is detuned from the ground state to Rydberg transition by ∆ m and ∆ l , see Fig. S11a.The detuning has two contributions ∆ m,l = −3V + ∆m,l : on the one hand side the term −3V is essential for the LPG protection and appears when rewriting the Pauli spins as τ x i,j = 2n i,j − 1.On the other hand, the term ∆m,l V enables to arbitrarily tune the electric field h and chemical potential µ in Hamiltonian (4).
The LPG term requires the Rydberg-Rydberg interaction between matter and link atoms at each vertex j to have the same strength 2V , see Fig. 1a.The interaction strength between two atoms both excited to a Rydberg state scales as C 6 /R 6 for large distances, where C 6 is a constant containing details of the internal atomic structure and R is the distance between the two atoms.Therefore, to have equal interaction strength between link and matter atoms, we require a tetrahedron geometry at each vertex as shown in Fig. 1a and S11b.We define a as the length of the edge of the honeycomb lattice.Link atoms are located on the center of edges and neighboring link atoms have a distance of a and thus we require the plane of matter and link atoms to have a distance of √ 3a/2.Since the honeycomb lattice is bipartite, we suggest to lift matter atoms in sublattice A (B) up (down) to decrease next-nearest neighbor interactions between matter atoms.These undesired next-nearest neighbor interactions can be estimated, see Fig. S11b, and we find that the nearest neighbor matter-matter interaction has strength V m−m = V /64 = 0.02V and the next-nearest neighbor link-link interaction has strength V l−l = 64V /729 ≈ 0.09V .Note that in an experimental setup the system has boundaries and here we want to discuss the LPG term at the boundary, see Fig. S11c.In particular, we want to consider the case where the boundary cuts through links of the honeycomb lattice, i.e. the boundary vertices have missing link atoms.To this end, we introduce dummy atoms in Fig. S11c which are not real atoms but instead substitute the missing link atom in the LPG term Eq. ( 2) by setting the corresponding dummy link to a constant value τ x dummy = +1.This yields additional detuning terms for all atoms on the boundary vertices.For technical purposes, instead of individually addressing the detuning on only boundary sites one could introduce real auxiliary atoms on the boundary, which are excited to a Rydberg state different from the matter and link atoms.This allows to shift the Rydberg state of the matter and link atoms on the boundary, i.e. adding the required detuning, while the auxiliary atom is not affected by the driving field Ω m,l .

VI. DISORDER-FREE LOCALIZATION
Disorder-free localization (DFL) is a phenomenon that has been studied in theories with local symmetries [33,34,[47][48][49][50].Here, we use the microscopic model (3) -as it would be implemented in an experimental setup -and show DFL behaviour in a small-scale exact diagonalization (ED) simulation for the case of U (1) matter.In particular, the observation of DFL would be an accessible experimental probe since it only requires to prepare the system in two different initial product states and time-evolve them under the microscopic Hamiltonian.
The key idea is the following: Consider a system with two subsytems A and B and an initial state, where all matter sites in subsystem A (B) are occupied (empty).We let the state time-evolve under a Z 2 invariant Hamiltonian and ask whether the matter excitations stay localized in subsystem A or delocalize equally across subsystem A and B. Hence, the quantity of interest is the time-averaged imbalance I(t) of matter excitations at time t between subsytems given by where nA(B) (τ ) is the expectation value of total matter excitations in subsystem A (B) at time τ and L denotes the system size.The eigenstate thermalization hypothesis (ETH) claims that I(t) eventually approaches its thermal equilibrium value.In DFL, this hypothesis is believed to be broken for gauge-noninvariant initial states |ψ ninv .In the following sections we want to examine this behaviour on the example of 1) the minimal "Mercedes star" lattice with coordination number z = 3 [(2 + 1)D] as presented in the main text and 2) on a Zig-Zag chain [(1 + 1)D], which would be experimentally easier to reach large system sizes.In both examples, we performed ED studies of the microscopic Hamiltonian (3) in small, numerically accessible systems.

A. Mercedes star
Let us first consider the gauge-invariant initial state |ψ inv with Ĝj |ψ inv = +|ψ inv and matter excitations being distributed as described above.The Mercedes star model with four matter sites and six links is illustrated in the inset of Fig. 3a.For this initial state, we can compute the imbalance from the thermal ensemble as predicted by ETH (see also Ref. [33,34], SI B) and find that indeed the system fully delocalizes I thermal = 0. Comparing this to the ED results in Fig. 3a, we find that the time-average imbalance quickly vanishes as expected.
The situation changes for the gauge-noninvariant initial state |ψ ninv .For this state the matter excitations should be located again only in subsystem A but the links should be in τ z i,j = +1 eigenstates as indicated in Fig. 3a.Therefore, |ψ ninv is an equal superposition of all possible gauge sectors g j = ±1.While still I thermal = 0, as we have verified numerically, we find that the state does not thermalize under time-evolution with Ĥmic as shown in Fig. 3a.
For dynamics under a perfectly gauge-invariant Hamiltonian, [ Ĥ, Ĝj ] = 0 ∀j, the intuitive picture is the following: When the system is initialized in a superposition of all gauge sectors, the system independently time evolves in each of these (uncoupled) gauge sectors.In each gauge sector the system has different background charges, which affects the spreading of the matter excitations, and ultimately the average is taken over all these gauge sectors.This effectively induces disorder in the system, which can lead to localization [47].3).In panel a) we show the same calculation as in Fig. 3a but for much longer, experimentally inaccessible times.We see that the plateau at short times is a pre-thermal behaviour followed by several smaller plateaus.Eventually the imbalance of the gauge-noninvariant states slowly decays.This can be explained by the approximate but not exact local symmetry of the system.In panel b), we show analogous calculations to those in panel a) or Fig. 3a for a Zig-Zag chain.We find again distinctly different thermalization behaviour for the two different initial states.
While this interpretation holds for Z 2 gauge-invariant Hamiltonians, we want to note the differences to our scheme here.Firstly, we generate the local symmetries in our system through local pseudogenerators Ŵj , which yield an even enriched symmetry structure because the system has three emerging local symmetry sectors as shown in Fig. 1c.This has been discussed previously [34] for (1 + 1)D systems.Secondly, the microscopic Hamiltonian explicitly breaks the local symmetry generated by the LPG terms, [ Ĥmic , Ŵj ] = 0 ∀j, to induce dynamics in the system (note that the gauge symmetry is only approximate and not exact).Hence, the weak drive has to be considered as an error term, which eventually leads to thermalization of the system for long times.However, at experimentally relevant timescales we find a clear pre-thermal plateau indicating DFL as shown in Fig. S12a.The parameters used in the ED calculation are shown in the inset of Fig. S12a.

B. Zig-Zag chain
The Mercedes star model is a numerical toy model with coordination number z = 3.A truly (2 + 1)D model with z = 3 can be realized on the honeycomb lattice but requires a large number of qubits/atoms.Therefore, as a first step to probe the proposed model, we suggest to implement a Zig-Zag chain with periodic boundary conditions, where each site of the chain is connected to a dummy atom, see Fig. S12b and SI V.This additional dummy atom ensures coordination number z = 3 such that the LPG protection scheme becomes fully applicable.
In Fig. S12b, we show results of an ED study in a Zig-Zag chain with four matter sites and four links with periodic boundary conditions.Again, we can define a subsystem A (B), where matter excitations are located at time t = 0. Furthermore, the gauge-invariant state |ψ inv and gauge-noninvariant state |ψ ninv differ by the configuration of the link atom.For both initial states, we expect the thermal expectation value of the imbalance, Eq. (S37), to vanish I thermal = 0 Again, we find a clearly different behavior after time-evolving under the microscopic Hamiltonian Ĥmic and evaluating the time-averaged imbalance.However, the observed dynamics is slower than in the Mercedes star model.The parameters used in the ED calculation are shown in the inset of Fig. S12b.

VII. SCHWINGER EFFECT
The Schwinger effect is a non-perturbative effect from quantum electrodynamics that describes the production of matter excitations from vacuum [60].Due to the weak coupling constant of quantum electrodynamics this effect is only expected to appear at very strong electric fields and has not been observed.LGTs were originally introduced to study the effects of strong coupling in gauge theories [1] and are therefore a candidate theory to also study the physics We have checked that the maximal evolution time tV = 6000 presumably captures the peak of each time trace.In the main text, Fig. 3b, we plot the peak value for all parameters (∆m, ∆ l ).
In panel b), we show the same result as in Fig. 3b, but we do not project into the target gauge sector.For large positive detunings ∆ l , we can see additional, unphysical resonance lines.In particular, the region where we claim gauge-invariant pair creation processes to appear remains unchanged.Therefore, the appearance of resonance lines in these regions are not driven by gauge-symmetry breaking processes.
of the Schwinger effect.Recently, digital quantum simulation of the (1 + 1)D Schwinger model on the lattice have examined the Schwinger effect [51].Here, we want to present an experimentally measurable quantity by considering the pairwise production of matter excitations.While the potential between two static charges is of theoretical interest and used as a signature of the Schwinger effect [61], it is numerically challenging to extract in (2 + 1)D and is a topic for future studies.Our effective model (4) allows to explore the Schwinger effect in a Z 2 mLGT in (2 + 1)D, which has not yet been observed.The experimental protocol starts by initializing a gauge-invariant state |vac without any matter excitations in one of two vacua, i.e. either all links in τ x i,j = +1 or τ x i,j = −1, which is a simply product state.Then, the system is quenched with the microscopic Hamiltonian (3) for time t yielding |ψ(t) = e −i Ĥmict |vac .The effective model (4) with quantum-Z 2 matter, which we expect to correctly capture the physics, contains pairing terms ∝ (â † τ z â † + H.c.) yielding pair creation processes from vacuum.As soon as the matter excitations are created, they move apart due to the hopping term and interact with the gauge field, which makes the prediction of the dynamics very challenging.
Therefore, an easily accessible quantity to probe the Schwinger effect is the gauge-invariant production rate of matter excitations.To this end, we let the system time-evolve and calculate which projects onto the gauge sector g j = +1 ∀j and gives the expectation value of the total number of matter excitations in the system.We calculate P(t) for t ∈ [0, 6000/V ] for different electric fields ∆ m and chemical potentials ∆ l , see Eq. ( 3), for Ω m = Ω l = V /8.In Fig. S13a we plot the timetraces for some set of parameters (∆ m , ∆ l ) and we have checked that the maximal time-evolution time of 6000/V captures the peak of each timetrace.
Each timetrace has different amplitude and timescale and therefore only considering P(t) at a fixed time t = t 0 is not sufficient to extract the productivity of creating matter excitations from vacuum.To this end, we take the maximum value of each timetrace max t P(t) for each (∆ m , ∆ l ), which is plotted in Fig. 3b in the main text.
We find that for some (∆ m , ∆ l ) the production of matter excitations is significantly higher than for others.An intuitive picture is that a pair of matter excitations costs an energy 2µ (mass of matter excitations) and due to Gauss's law, the two matter excitations have to be connected by an electric string, which costs 2h.Therefore, if 2µ + 2h = 0 the process is on resonance we expect a high production rate of matter pairs.Note that µ, h have to be determined from Tab. SII and are not simply given by (∆ m , ∆ l ).
Besides the resonances described above, there are several more processes that would have to be taken into account.E.g. two neighboring matter excitations repel each other with strength M ; Eq. (II C); matter excitations are mobile; plaquette interaction compete with the electric field; finite size effects etc. Due to these competing interactions, it is very hard to gain a complete picture and thus large-scale numerical simulations as well as experimental observations are needed.Nevertheless, small-scale numerical calculations of the microscopic model show promising signatures of the Schwinger effect.
An additional feature we monitored in our ED study is the role of gauge-noninvariant dynamics.In Eq. (S38) we project out unphysical states.In contrast we can do the same analysis as before but without projecting on the target gauge sector, which is shown in Fig. S13b.We find regions with additional resonances, which are caused by gauge-symmetry breaking processes.However, the resonances we find with projecting on the gauge sectors are not altered, which is an evidence that the physics is purely determined by Z 2 -invariant dynamics.

VIII. DMRG IN THE LADDER
In Fig. 2 in the main text, we map out limiting cases of the ground-state phase of the effective model (4).A numerically more accessible and experimentally realizable model is the Z 2 mLGT coupled to U (1) matter on a ladder.While the ladder geometry is not 2D but mixed-1D, it has coordination number z = 3 and therefore is applicable for our proposed LPG term (2).Thus, the effective Hamiltonian only has to be modified for the plaquette interaction because the plaquettes on a ladder have four instead six edges; note that we anticipate the plaquette terms to be even stronger in the ladder and we choose J < 0 as found in SI II B in Eq. (S16).
For our numerical simulations, we use Hamiltonian (4) with ∆ 1 = ∆ 2 = 0, i.e. we have a global U (1) symmetry for the matter, and we tune the electric field h and chemical potential µ for fixed tunneling and plaquette interactions t = 1 and J = −1.Using the density matrix renormalization group (DMRG) technique, we calculate the ground state of the above described Hamiltonian on a ladder with L = 19 plaquettes.
From the ground state calculated in DMRG, we obtain the average matter excitation density -an experimentally directly accessible quantity, e.g. by taking snapshots in the atomic ground state and Rydberg basis (see Fig. 1a).As shown in Fig. 3c, the system is in a matter vacuum for large h/J similar to the 2D case.For decreasing h/J, we find a sharp increase of matter excitations indicating a phase transition as shown in Fig. S14a.At the same critical electric field value, the plaquette term shows a sharp feature; note that J < 0 and thus the plaquette expectation value is negative indicated by the reversed sign on the plot label in Fig. S14b.
Moreover, we show a similar scan of parameters but now we fix the electric field h/J = 0.7 while scanning the chemical potential µ/J.We find consistent behaviour with Fig. 2a, i.e. a sharp transition into the vacuum phase, see Fig. S14c and d.
To characterize the different phases and its phase transitions requires more elaborate studies of our effective model on the ladder and is a topic for future studies.We emphasize that due to its numerical accessibility and experimental feasibility, the ladder model is a promising playground to probe Z 2 lattice gauge theories coupled to dynamical matter beyond (1 + 1)D.
In the following, we discuss numerical details of the DMRG calculation.We use the TeNPy package [62,63] to find the ground state of Hamiltonian (4) with ∆ 1 = ∆ 2 = 0. Note that while this Hamiltonian conserves the particle number, we do not run the DMRG simulation in a fixed particle number sector.For a given set of Hamiltonian parameters, we find the global ground state, and can therefore use the average matter density as an observable.In the DMRG simulation, we enforce the system to be in the target gauge sector by adding a large energy penalty term proportional to Gauss's law (1) to the Hamiltonian.The ground state will therefore always fulfill Gauss's law by construction.We carefully checked our numerical results for convergence, see Fig. S14.We have confirmed that in our simulations the system thermalizes after 200×L 2 steps for all T /h.In panel b), we show the autocorrelation of the percolation strength.We plot the average over 15 runs with 10 4 samples, respectively.Between each sample we perform 2×L 2 steps.We find negligible autocorrelation between samples.We show results for Hamiltonian (S39) at M/h = 0.2917, χ1/h = 0.1483, χ2/h = 0.073, χ3/h = 0.4347 and two matter excitations.We plot the number of strings in the largest string-cluster, the percolation strength, the Euclidean distance between the two matter excitations and the string number (from left to right).The size of the honeycomb lattice is 10×10 (top) and 35×35 (bottom).We can clearly identify features for thermal deconfinement at (T /h)c ≈ 2, which become sharper for increasing system sizes.
the Euclidean distance between the two matter excitations.We illustrate two snapshots from the deconfined and confined regime in Fig. S15b.
In Fig. S17 we show Monte Carlo results for system size 10×10 and 35×35.We observe a clear change of behaviour in all above discussed quantities, which signals a thermal deconfinement phase transition.For low temperatures, the percolation strength vanishes.At a critical temperature (T /h) c ≈ 2, the percolation strength abruptly increases, i.e. the string-net percolates and the matter excitations are deconfined.At the same critical temperature (T /h) c ≈ 2, the Euclidean distance between the two matter excitations drastically increases to roughly the system size.We note that finite-size effects strongly influence (T /h) c for small system sizes.However, the transition becomes generally sharper for larger lattice sizes as expected.

FIG. 1 .
FIG.1.Constraint-based implementation of Z2 mLGT with qubits.The Z2 gauge structure emerges from the dominant local-pseudogenerator (LPG) interaction on the honeycomb lattice introduced in panel a).A vertex contains matter âj qubits (blue) and shares link τ x i,j qubits (red) with neighboring vertices.All qubits connected to a vertex interact pairwise with strength 2V .In a Rydberg atom array experiment the qubits are implemented by individual atoms in optical tweezers, which are assigned the role of matter or link depending on the position in the lattice.Here, the ground-and Rydberg state of the atoms, |g and |r , encode qubit states, which are coupled by an off-resonant drive Ω to induce effective interactions.To realize equal strength nearest neighbor, two-body Rydberg-Rydberg interactions, the matter atoms can be elevated out of plane.In panel b) we introduce the notation for the Z2 mLGT, for which the Hilbert space constraint is given by Gauss's law Ĝj = +1.We illustrate the electric field τ x i,j = +1 (τ x i,j = −1) with flat (wavy) red lines and the matter site occupation n j = 0 (n j = 1) with empty (full) blue dots.Panel c) shows the notation for the QDM subspace with exactly one dimer per vertex.Panel d) illustrates how the distinct subspaces are energetically separated by the LPG term V Ŵj .The two quantum dimer subspaces are disconnected when the matter is static, which can be exactly realized by the absence of matter atoms in panel a) and setting (2â † j âj − 1) = ±1 in V Ŵj .

FIG. 3 .
FIG.3.Experimental probes.We analyze several observables that could be probed experimentally.Panel a) and b) show results from ED simulations of the time-evolution of the microscopic model (3) with experimentally realistic parameters in a system with coordination number z = 3 (see inset).In panel a) we observe disorder-free localization by initializing the system in a gauge-invariant (blue curve) and gauge-noninvariant (red curve) initial state with two matter excitations localized in subsystem A and calculating the time-averaged imbalance between subsystem A and B as shown.In panel b), we probe the Schwinger effect by quenching the vacuum state with the microscopic model for different experimentally relevant parameters: matter detuning ∆m (chemical potential) and link detuning ∆ l (electric field).We find lines of resonance, where the production of matter excitations out of the vacuum is large.In panel c) we plot the average U (1) matter density (blue curve) obtained from DMRG calculations on a ladder with J < 0. We can qualitatively understand the sharp decay of matter as a transition into the vacuum phase as discussed in Fig.2a.Additionally, a kink in the plaquette expectation value (red curve) signals a phase transition.In panel d), we use two fluctuating test charges to probe a temperature-induced deconfinement transition in a classical limit of our effective model using Monte Carlo simulations.Both in the percolation strength (red curve) and the Euclidean distance of two matter excitations (blue curve), we find that above a certain temperature T /h the system undergoes a percolation transition.

FIG. 4 .
FIG.4.Disorder-based protection scheme.We calculate the spectrum of the minimal model studied in Fig.3a)-b) with Ω = 0 and plot all eigenstates around energy E = 4V .Green (red) dots are states that fulfil (break) Gauss's law as illustrated with two examples in the inset of panel a).Without disorder, i.e.V j = V for all j, the physical and unphysical states are on resonance.In panel b), we show the effect of disordered protection terms V j = V + δV j , which only shifts the unphysical states out of resonance and hence fully stabilizes the gauge theory.We note that even without disorder, the emergent gauge structure is remarkably robust (Supplementary note 2 G).

B.
Plaquette terms for U (1) matter: V = 2|∆m| ∆ l , Ωm, Ω l FIG. S5.Estimation of plaquette terms (quantum-Z2 matter).The effective plaquette interaction derived by a Schrieffer-Wolff transformation depends on the matter and electric field configuration within each plaquette.In panel a) we plot the absolute value of the coupling strength for all 416 different configurations for various driving strength Ω/V (from dark to bright shade: V /Ω = 8,10,12,14,20,30).We find a plateau with strongest coupling for certain staggered and polarized electric field configurations as shown in the inset.Note the sign of J stag FIG.S7.Microscopic versus effective model.We show results of exact diagonalization studies in a minimal model and compare the microscopic model Eq.(3) to the effective Hamiltonian with quantum-Z2 matter, Eq. (S23), and U (1) matter, Eq. (S13).For the first (latter) case we choose ∆ l = ∆m = 0 (∆ l = 0, ∆m = 0.5V ).In panel a) and b), we calculate the expectation value g = j Ĝj for each eigenvector of the microscopic model and our target sector is g = 4.In panel b), we additionally require the number of matter excitations to be conserved.For Ω = 0.1V , we find that both the local and global symmetry emerge in the microscopic model.In panel c) and d), we now consider the eigenenergies in the target sector for different driving strength Ω/V for the microscopic (blue) and effective (red) model.We find agreement of both spectra which supports the validity of our perturbative approach discussed in the main text and SI II.

|
FIG. S9.Effective meson model.Panel a): For J/h, t/h → 0, the matter excitations are tightly bound into mesonic pairs b on the honeycomb lattice.These mesons are again hard-core bosons and Z2 charge neutral.In panel b) and c), we describe the leading order second-order processes derived from the Z2 mLGT, which gives rise to hopping t eff of mesons as well as repulsive interactions |δ eff | due to the absence of dispersive shifts for next-nearest neighbor mesons.Moreover, the plaquette interaction of the Z2 mLGT yields to fluctuating mesons for plaquettes with exactly three mesons as depicted in panel d).Panel e) summarizes the effective meson model.The model is described by hopping ∝ t eff of hard-core bosons b (black circles) on a Kagome lattice with infinitely strong nearest-neighbor repulsion (from the hard-core constraint on the honeycomb lattice), finite NNN repulsive interactions ∝ |δ eff | and plaquette interactions ∝ J. S9e . The infinite repulsive term comes from the hard-core boson constraint of single matter excitations.Therefore, the effective meson model is given by Ĥmeson = t eff n,m PNN b † n bm + H.c. PNN − J P PNN + H.c. PNN + |δ eff | n,m b † n bn b † m bm (S32) where n, m denote sites of the Kagome lattice as shown in Fig. S9e and the projector PNN ensures the constraint that no nearest neighbor mesons can exist.Here, the notation •, • ( •, • ) describes (next-)nearest neighbors on the Kagome lattice.
FIG.S10.Derivation of plaquette terms for the QDM.We show the spectrum of the microscopic model (3) on one plaquette for different driving strengths in the flippable plaquette subspace.The blue dashes show the unperturbed system.The two low-energy states are the two flippable plaquette configurations in the QDM sector.The first (second) excited manifold has two (four) monomers and the two high-energy states contain a maximum number of six monomers.To perform perturbation theory, we first need to diagonalize the degenerate subspaces because the drive couples within the two and four monomer excitation blocks, i.e. the drive can move around the monomer excitations without energy cost, thus making them mobile.The orange and green dots show the spectrum for different driving strength, where the blocks are diagonalized.However, there are still off-diagonal couplings between the blocks.These couplings are the starting point for the second step, which is the actual Schrieffer-Wolff perturbation theory.The inset illustrates the smallest energy gap and its renormalized coupling after the first step.This allows to compare the effective driving strength to the energy gap in order to determine a regime of validity for the perturbation theory.

|
FIG. S11.Rydberg atoms in tweezer array.Panel a): The ground state |g and Rydberg state |r of the matter (link) atom is mapped on the matter field n (electric field τ x ) of the Z2 mLGT.Panel b): To stabilize a gauge sector of the Z2 mLGT or to enforce the hard-core dimer constraint in the QDM, we have introduced the LPG protection (2) which requires two-and one-body interactions.The interaction strength can be adjusted by the geometry in a Rydberg atom array because the strength of the dipolar Rydberg-Rydberg interaction depends on the interatomic distance.Here, we show the suggested geometry to 1) realize the LPG term and 2) minimize the effects of long-range interactions.Panel c): A boundary in the finite size experimental setup cuts through links of the lattice.We introduce dummy atoms, which are not real and only used to illustrate the modification of the LPG term on the boundary.The LPG term on the boundary has less link atoms than in the bulk.Hence, the detuning on the matter and link atoms along the boundary has to be adjusted.

V 5
FIG.S12.Disorder-free localization.We show results from small-scale ED studies, where we time-evolve gaugeinvariant |ψ inv and gauge-noninvariant initial states |ψ ninv under the microscopic Hamiltonian (3).In panel a) we show the same calculation as in Fig.3abut for much longer, experimentally inaccessible times.We see that the plateau at short times is a pre-thermal behaviour followed by several smaller plateaus.Eventually the imbalance of the gauge-noninvariant states slowly decays.This can be explained by the approximate but not exact local symmetry of the system.In panel b), we show analogous calculations to those in panel a) or Fig.3afor a Zig-Zag chain.We find again distinctly different thermalization behaviour for the two different initial states.
FIG. S13.Schwinger Effect.We provide complementary plots to the results shown in Fig. 3b calculated from exact diagonalization studies of the microscopic model (3) on the Mercedes star.In panel a), we show the number of matter excitations in the target gauge sector, Eq. (S38), for some exemplary parameters (∆m, ∆ l ) versus time tV ; the parameters are chosen along the arrow on the left hand side of panel b).We have checked that the maximal evolution time tV = 6000 presumably captures the peak of each time trace.In the main text, Fig.3b, we plot the peak value for all parameters (∆m, ∆ l ).In panel b), we show the same result as in Fig.3b, but we do not project into the target gauge sector.For large positive detunings ∆ l , we can see additional, unphysical resonance lines.In particular, the region where we claim gauge-invariant pair creation processes to appear remains unchanged.Therefore, the appearance of resonance lines in these regions are not driven by gauge-symmetry breaking processes.
FIG.S14.Convergence of DMRG in the ladder.We plot results of the DMRG calculations for different bond dimensions χ = 100 and χ = 200.Panel a) and b) show the expectation value of the matter density and plaquette terms, respectively, as presented in the main text in Fig.3cfor fixed chemical potential µ/J = 1; note that J < 0 and hence µ < 0. In panel c) and d), we keep the electric field fixed at h/J = 0.7 and vary the chemical potential.
FIG. S16.Monte Carlo thermalization & autocorrelation.We show results for Monte Carlo simulations of Hamiltonian (S39) at M/h = 0.2917, χ1/h = 0.1483, χ2/h = 0.073, χ3/h = 0.4347 and two matter excitations for T /h = 3.In panel a) we show the thermalization of the percolation strength (top) and the Monte Carlo weights averaged over 15 runs (bottom).We have confirmed that in our simulations the system thermalizes after 200×L 2 steps for all T /h.In panel b), we show the autocorrelation of the percolation strength.We plot the average over 15 runs with 10 4 samples, respectively.Between each sample we perform 2×L 2 steps.We find negligible autocorrelation between samples.
FIG. S17.Monte Carlo results.We show results for Hamiltonian (S39) at M/h = 0.2917, χ1/h = 0.1483, χ2/h = 0.073, χ3/h = 0.4347 and two matter excitations.We plot the number of strings in the largest string-cluster, the percolation strength, the Euclidean distance between the two matter excitations and the string number (from left to right).The size of the honeycomb lattice is 10×10 (top) and 35×35 (bottom).We can clearly identify features for thermal deconfinement at (T /h)c ≈ 2, which become sharper for increasing system sizes.
Local pseudogenerators for Z 2 mLGTs.-Theimplementation of LGTs in quantum simulation platforms have two inherent challenges to overcome: METHODS

TABLE SI .
Overview and Summary of the Supplementary Information.