Reliability of lattice gauge theories

Currently, there are intense experimental efforts to realize lattice gauge theories in quantum simulators. Except for specific models, however, practical quantum simulators can never be fine-tuned to perfect local gauge invariance. There is thus a strong need for a rigorous understanding of gauge-invariance violation and how to reliably protect against it. As we show through analytic and numerical evidence, in the presence of a gauge invariance-breaking term the gauge violation accumulates only perturbatively at short times before proliferating only at very long times. This proliferation can be suppressed up to infinite times by energetically penalizing processes that drive the dynamics away from the initial gauge-invariant sector. Our results provide a theoretical basis that highlights a surprising robustness of gauge-theory quantum simulators.

Introduction.-In modern physics, gauge theories assume a central role, ranging from the Standard Model of Particle Physics [1] to emergent exotic solid-state phases [2,3]. The defining feature of any gauge theory is its local conservation laws, such as Gauss's law for a U(1) gauge symmetry. Despite their elegance, the computation of gauge theories on classical computers is a daunting task [4,5], particularly for out-of-equilibrium phenomena. Currently, a complementary tool to probe gauge theories is emerging in the form of low-energy tabletop devices, so-called quantum simulators [6][7][8][9]. Although still in the developmental phase, these experiments are rapidly advancing [10][11][12][13][14][15][16][17][18], making one open issue all the more pressing: how can we ensure the reliability of quantum simulators once they scale beyond problems that can be benchmarked by classical computers [19,20]? For gauge theories, this issue is particularly subtle, since a faithful quantum simulation necessitates not only the correct engineering of the Hamiltonian dynamics, but crucially also of the defining local gauge symmetry. It becomes thus an outstanding challenge to understand how gauge-invariant dynamics may be faithfully simulated -without unrealistically fine-tuned interactions between the constituents of the quantum simulator and in the consequently unavoidable presence of gauge-violating errors.
As we show in this Letter, quantum simulators can reliably reproduce the out-of-equilibrium dynamics of gauge theories even if the prohibitive restriction of perfect gauge invariance is relaxed. Our results are based on numerical studies of Z 2 and U(1) gauge theories as well as analytic proofs. In the presence of gauge invariance-violating terms, the short-time dynamics deviates only perturbatively slowly from the idealized scenario, and observables reproduce the ideal dynamics during a time frame controlled by the strength of the gauge-violating errors. Hamiltonian H0 on L matter sites, undesired processes ∝ λH1 that break gauge symmetry drive the dynamics away from the initial gauge-invariant sector (green bubble) given by the analog of Gauss's law, Gj |ψ = 0 ∀j to other sectors in the total Hilbert space (blue bubble). At late times, the average gauge violation assumes a random value. Left: Through the introduction of a protection term ∝ V HG that energetically penalizes gauge-violating processes, gauge invariance is retained up to infinite times. Two sharply distinct regimes emerge: an uncontrolled-violation regime where V is too small to energetically isolate the initial gauge-invariant sector, and a controlled-violation regime for sufficiently large V where the infinite-time violation scales as ∼ (λ/V ) 2 . The axes are both log-scale. See Fig. 4 for a detailed quantitative presentation on Z2 (p = 1) and U(1) (p = 2) gauge theories.
Upon adding a protection term that energetically penalizes such errors, gauge invariance becomes bounded from above even up to infinite times; cf. Fig. 1. As we explain below and in the Supplemental Material (SM) [21], this protection can be mathematically understood as an arXiv:2001.00024v1 [cond-mat.quant-gas] 31 Dec 2019 emergent, deformed gauge symmetry [22]. In surprising contrast to expectations from perturbation theory, our numerics suggests that the strength of the protection term does not need to scale with system size.
Our results thus form a theoretical basis for some previous works, which have found evidence that quantum simulators may approximately retain gauge invariance [23][24][25][26][27]. Moreover, our work complements existing results on equilibrium theories: Gauge-invariant equilibrium phases can emerge in a low-energy effective theory, even if the microscopic description breaks gauge invariance, e.g., in topological phases of matter [28], when the gauge degree of freedom decouples because of a large mass [29], or when gauge-noninvariant terms at a small scale renormalize away at large distances ("light from chaos") [30][31][32].
Model.-Here, we focus on out-of-equilibrium dynamics, which is highly pertinent for current quantum simulators in ultracold atomic gases [14,[16][17][18] and for which emergent gauge symmetry has received considerably less attention. Although our discussion is general, for concreteness, we consider a Z 2 gauge theory [33][34][35][36] coupled to matter inspired by a recent experiment [16] we mostly consign similar results for a U(1) gauge theory to the SM [37].
The Z 2 gauge theory lives on a lattice of L matter sites in one spatial dimension [38], described by the Hamiltonian with periodic boundary conditions. Here, τ {x,z} j,j+1 are spin-1/2 Pauli matrices, where the z (x) component stands for the Z 2 gauge (electric) field on the link between matter sites j and j + 1. The matter fields are represented by hardcore bosons with the annihilation operator a j on site j. The dynamics of the matter field couples to the Z 2 gauge field with strength J a . The electric field has energy J f . Gauge invariance is encoded in a set of local symmetry generators with eigenvalues g j and where Q j = 1 − 2a † j a j is the charge of the matter field. As a Z 2 lattice equivalent of Gauss's law, the G j commute with the model Hamiltonian [H 0 , G j ] = 0, ∀j. For a perfectly gauge-invariant system, the Hilbert space thus decouples into different symmetry sectors with fixed local charge g j = 0, 2.
In the Standard Model of Particle Physics, such local gauge invariance is postulated at a fundamental level. In a quantum simulator, however, gauge invariance needs to be engineered. Some experiments have used Gauss's law to integrate out either matter or gauge fields, which yields an effective spin theory, suitable for implementation in a quantum computer and with encoded gauge invariance [10,11,15]. However, this approach is viable only in one spatial dimension. If, in contrast, both matter and gauge fields are retained as active degrees of freedom, such as in recent ultracold-atom experiments [16][17][18], exact gauge invariance would require fine-tuning at unrealistic levels of accuracy. Consequently, such quantum simulators will always have some inherent gauge invariance-breaking processes.
To study the severity of such terms, we add for concreteness the error Hamiltonian These terms are inspired by Ref. [16], though our conclusions do not depend on their precise form. Here, λ is a parameter representing an adjustable gauge-noninvariant error strength and the constants c l are modeled after experimental parameters [39]. Generically, the error term of Eq. (3) will drive the dynamics out of the gaugeinvariant subspace, such that after a certain time the expectation value of the Gauss-law operator will be the same as in a random state from the entire Hilbert space (i.e., the space that contains gauge-invariant states as well as all states with G j |ψ = 0). Quench dynamics with gauge violation.-We mimic a typical quantum simulator experiment, where the system is initialized in a simple product state |ψ 0 in the gauge-invariant sector G j |ψ 0 = 0, ∀j, which is then quenched with the Hamiltonian H 0 + λH 1 . We choose |ψ 0 such that a † j a j = [1 + (−1) j ]/2 and τ x j,j+1 = (−1) j+1 , and, following Ref. [16], we set J a = 1 and J f = 0.54 (our conclusions do not depend on these precise choices). Our numerical results are obtained using the QuTiP [41,42] and QuSpin [43,44] exact diagonalization toolkits. For time evolution, we have opted to use our own exact exponentation routine in lieu of these toolkits' time-evolution functions, given the extremely long times we simulate, for which solvers based on integrating ordinary differential equations are not optimal.
To evaluate the effect of H 1 on the gauge-invariant dynamics, we consider the dynamics of the spatiotemporal averages of the violation of Gauss's law, the magnetization in the x direction (the 'electric field') as well as the staggered boson number In the above, |ψ(s) = exp[−i(H 0 + λH 1 )s] |ψ 0 . The respective deviations from the ideal gauge-invariant case are denoted by ∆N stag (t) and ∆m x (t). Further observables are discussed in the SM [40]. Note that the gauge violation in Eq. (4) suffices for the Z 2 gauge theory, but would require an even power in G j or its absolute value for the U(1) gauge theory [37]. Our numerical results are presented in Fig. 2. The gauge violation grows only gradually, ε(t) ∼ (λt) 2 at short times, before it saturates at long times. This subleading increase of the gauge violation directly stems from the fact that j G j is gauge invariant and also commutes with H 0 , which necessarily leads to a vanishing  sites. We see two clear regimes: at small V the gauge-noninvariant behavior dominates after a timescale ∝ 1/λ. When V is sufficiently large (controlled-violation regime), the gauge-invariance violation in (a) is indefinitely suppressed by (λ/V ) 2 , while deviations in the electric field in (c) are suppressed by λ/V up to a timescale ∝ V /λ 2 . The deviation in the staggered boson number becomes uncontrolled only at times beyond t ∝ V /λ 2 . See SM [40] for corresponding results at λ = 0.005 and 0.5.
first-order contribution in perturbation theory; cf. SM [45]. In contrast, the short-time scaling behavior can change for gauge-invariant observables that do not commute with H 0 , such as N stag and m x . Even though ∆m x scales as (λt) 2 at short times, as seen in the inset of Fig. 2(b), ∆N stag shown in the inset of Fig. 2(c) scales instead as λt 2 , with the latter emenating from the firstorder term in perturbation theory; cf. SM [45]. Importantly, N stag (t) remains nevertheless close to the gaugeinvariant dynamics up to t ≈ 1/λ. Thus, there is a clear time frame over which gauge-noninvariant terms do not compromise observable properties, and by decreasing λ this time frame can be improved in a controlled manner.
Quench dynamics with energy protection.- We see two separate regimes where at large V the violation is controlled and scales as (λ/V ) 2 , while at small V the violation is uncontrolled. The dotted lines indicate the largest V at which there is an energy overlap between the initial gaugeinvariant sector, and other symmetry sectors, which we protect against. In the controlled-violation regime, the violation is system size-independent. The vertical dotted lines correspond to the largest value of V (for a given system size L) at which the initial gauge-invariant sector still has spectral overlap with other sectors.
Some promising proposals have suggested to perturbatively generate the desired gauge-theory Hamiltonian by adding a term proportional to Gauss's law, which energetically penalizes violations of gauge invariance [7,8,24]. A similar term has been discussed for protecting stored quantum information [46]. We now systematically address the question of how such a term restores the ability to quantum simulate gauge-invariant dynamics.
To analyze this scenario, we introduce a protection term which energetically penalizes all violations of G j |ψ 0 = 0, with adjustable protection strength V . Note that Eq. (7) is not sufficient for the U(1) gauge theory, and should rather be V H G = V j G 2 j [37]. Our quench Hamiltonian is thus given by H = H 0 + λH 1 + V H G . Consequently, in Eqs. (4), (5), and (6) we now have Fig. 3, we show the effect of V on ε(t), ∆N stag (t), and ∆m x (t). As can be shown in degenerate perturbation theory [47], the violation ε(t) of Gauss's law is suppressed by (λ/V ) 2 at sufficiently large V ; cf. SM [40]. One might suspect that such an energy penalty would only suppress gaugenoninvariant processes at short times before the system eventually accesses all possible gauge-invariant sectors at long times. Interestingly, however, once V is on the order of a few λ, a "controlled-violation" regime is reached, where ε(t) is protected for all simulated times (which are many orders of magnitude larger than what is reachable in current experiments).
In contrast, for sufficiently large V the deviation of m x (t) is suppressed only by λ/V , as shown in Fig. 3(b), with the suppression ceasing at a timescale ∝ V /λ 2 ; cf. SM [40]. At the same timescale, N stag starts deviating from the gauge-invariant case if the values of V lie in the controlled-violation regime; see Fig. 3(c). Furthermore, the terms that drive these deviations can be described as gauge-invariant processes in degenerate perturbation theory, and could thus be absorbed into a renormalized, gauge-invariant H 0 . See SM [47] for more details. The existence of these timescales is a considerable improvement, since outside the controlled-violation regime deviations from the ideal gauge-invariant dynamics proliferate already at an earlier timescale ∝ 1/λ.
Physically, V protects gauge invariance by opening a large energy gap between the sector G j |ψ = 0, ∀j, and all other sectors. Since the bandwidth of each sector increases with system size, one could expect the necessity to scale V with L, which would invalidate this protection mechanism for large-scale quantum simulators. However, our numerics suggest that this is not the case. We illustrate this in Fig. 4 for the infinite-time per-site violation of Gauss's law as a function of λ/V (with fixed λ = 0.05). We see two clear regimes. The first at small values of V shows an uncontrolled violation that heavily depends on how the different gauge-invariant sectors are coupled to one another. The second regime, however, displays a controlled violation that scales as (λ/V ) 2 . This regime is expected once the gauge-invariant sector we start in is energetically well separated from other sectors, as can be shown in degenerate perturbation theory and made mathematically rigorous by adapting the results of Ref. [22] (see SM [47]): For any unitary symmetry that is broken on a scale ∼ λ, the opening of a gap generates an emergent symmetry that is perturbatively in (λ/V ) 2 close to the original one. Surprisingly, however, the scaling as (λ/V ) 2 sets in much before full separation between sectors is achieved. This onset appears to be largely independent of system size, contrary to the analytic arguments based on perturbation theory and the emergent gauge symmetry (see SM [47]). The reason is that the relevant gap to gauge-violating sectors is not to be counted relative to the entire gauge-invariant sector, but only to the energy region that is populated during the quench. A similar effect has been observed for the robustness of the ground-state degeneracy in topological matter [28].
Our results in Figs. 3 and 4 illustrate how once the controlled-violation regime is reached at a sufficiently large V , extrapolations to the ideal case become possible in both the gauge-invariance violation and gaugeinvariant observables. Indeed, Fig. 3 shows a clear timescale before which deviations in a gauge-invariant observable are well-determined as a function of λ and V . Even better, in this regime the control in the gaugeinvariance violation is not limited by any timescale, but rather persists indefinitely ∝ (λ/V ) 2 .
Conclusions and outlook.-We have carried out a thorough analysis, through exact diagonalization and perturbation theory, of the reliability of lattice gauge theories in out-of-equilibrium dynamics. We have found that small gauge-nonivariant processes (of strength λ) do not compromise the desired dynamics of observables up to a clear time frame ∼ λ −1 . Moreover, when introducing a sufficiently expensive energy penalty of strength V for such processes, the gauge-invariance violation enters a controlled-violation regime where it scales as (λ/V ) 2 , up to infinite times, and is also robust with respect to system size. This is a very encouraging result for experimental efforts on quantum simulators as it indicates that introducing an energy penalty leads to an indefinite protection of gauge invariance. This enables a well-defined extrapolation of a perfect gauge theory from gauge-invarianceviolating data.
The suppression in observables' deviations from their gauge-invariant dynamics presents a more varied picture. This is to be expected because when a gauge theory is broken by a small parameter, another gauge theory perturbatively close to it emerges [22] that, even though gauge-invariant, is still different from its initial counterpart. The original theory's exact dynamics can, however, be recovered by an appropriate absorption of the new terms into renormalized parameters [30,31]. tailed derivation in perturbation theory of the short-time dynamics due to a gauge invariance-breaking term of strength λ > 0. In this Supplemental Material, we discuss our analysis of the U(1) gauge theory, which supports our conclusions on the Z 2 gauge theory discussed in the main text. Moreover, we add supplemental results on the Z 2 gauge theory that further substantiate our conclusions. Further, we present details of our perturbation theory derivations, in addition to analytic arguments demonstrating the emergent gauge theory in the controlled-violation regime.

SUPPLEMENTAL RESULTS FOR Z2 GAUGE THEORY IN SUPPORT OF MAIN CONCLUSIONS
In the main text, we present in Fig. 3 results for the quench dynamics in the Z 2 gauge theory of various observables upon adding a protection term in the Hamiltonian in order to mitigate gauge invariance-breaking processes of strength λ = 0.05. In Fig. S1 we repeat these results for two more values of λ = 0.005 and 0.5. These results affirm the conclusions of the main text, and further corroborate the scaling of the violation and the suppression of the deviation of observables from their ideal gauge-invariant dynamics.
In the main text, we have resorted to using the spatiotemporal average ε(t) of the gauge-invariance violation j G j (t) given in Eq. (4). Another, more stringent, way to measure gauge-invariance breaking is through the projector where |j, 0; q is an eigenstate of the local gauge-invariance operator G j with eigenvalue 0, and where q denotes all remaining good quantum numbers. Now the violation in gauge invariance can be alternatively measured as (c) (f) Figure S1. (Color online). Same as Fig. 3 of the main text, but for gauge invariance-breaking strength of λ = 0.005 (left-column panels) and λ = 0.5 (right-column panels). The qualitative picture is unchanged form that of the main text.
When the initial state |ψ 0 is quenched by H 0 + λH 1 , i.e., when |ψ(s) = exp[−i(H 0 + λH 1 )s] |ψ 0 , the dynamics of ε proj (t) is qualitatively identical to that of ε(t). After an initial growth that scales as (λt) 2 at short times, ε proj (t) saturates at a timescale of 1/λ 2 , as shown in Fig. S3. Upon adding an energy protection term V H G , i.e., |ψ(s) = exp[−i(H 0 + λH 1 + V H G )s] |ψ 0 in Eq. (S3), we find that for sufficiently large V there is a controlledviolation regime, where the violation is suppressed as (λ/V ) 2 up to infinite times. Indeed, as shown in Fig. S3, the scaling behavior of the gauge-invariance violation as measured by 1 − P 0 (t → ∞) is qualitatively identical to that measured by j G j (t → ∞) /L in Fig. 4 of the main text.
Let us now look more closely at finite-size effects. In the main text, Fig. 4 indicated that in the controlled-violation regime the infinite-time gauge-invariance violation becomes system size-independent. It is thus interesting to see how the spatiotemporally averaged gauge-invariance violation ε(t) in Eq. (4), electric-field deviation ∆m x (t) in Eq. (5) and staggered boson-number deviation ∆N stag in Eq. (6) behave dynamically with respect to system size. As we can see in Fig. S4(a) for λ = 0.05 and V = 5 in the controlled-violation regime, ε(t) actually goes down with system size at long times, while ∆m x (t) in Fig. S4(b) exhibits volume convergence at larger system sizes. The picture is more subtle when it comes to ∆N stag in Fig. S4(c), but nevertheless it seems that the deviation does not significantly increase with system size when comparing the largest system sizes L = 6 and L = 8 matter sites.
Furthermore, we note that the coefficients c 1 , c 2 , c 3 , c 4 of Eq. (3) in the main text are inspired by the error terms discussed in Ref.
[S1]. The error terms are dominated by pair tunneling -processes where both species used in the two-component ultracold-atom experiment tunnel -and by detuning. We set the dimensionless driving parameter χ to various values, including ones that are similar to those in the experiment of Ref.
[S1], while always ensuring that c 1 + c 2 + c 3 + c 4 = 1. This is done to ensure a systematic control of the error strength through only the parameter λ in Eq. (3).

ANALYSIS OF THE U(1) QUANTUM LINK MODEL
We now show how the conclusions of the main text, which mostly focus on the Z 2 gauge theory, also hold for the U(1) gauge theory. We consider the U(1) QLM given by the Hamiltonian where again a j , a † j denote the ladder operators of the hardcore bosons representing the matter field, but now the gauge field is given by the spin-1/2 matrix s + j,j+1 linking sites j and j + 1, and the electric field by s z j,j+1 . The simplicity of this model makes it a perfect benchmark for numerical calculations. Yet, it hosts a number of relevant physical phenomena, such as Coleman's phase transition [S2], a dynamical quantum phase transition [S3], and extremely slow dynamics generated by gauge invariance [S4]. The Hamiltonian H 0 is invariant under unitary operations generated by the local Gauss's-law operator [H 0 , G j ] = 0. The gauge-invariant subspace is defined by states that satisfy G j |ψ = 0 ∀j. Proposals for quantum simulating this QLM have been put forward, e.g., in trapped ions [S5, S6], where H 0 is realized as an effective theory with gauge invariance only approximately preserved. To account for this situation, we add the following error term, inspired by (but more general than) Ref. [S6]: Here, in order to show that the features are independent of a specific choice of parameters, x j , y j , z j ∈ N * (z j > 1) and κ j , η j , γ j ∈ [0, 1] are randomly selected from a uniform distribution. Generically, as in the case of the Z 2 gauge theory discussed in the main text, the error term λH 1 will drive the dynamics out of the gauge-invariant subspace. In order to enforce the gauge invariance of the entire H, we will later also add H G to the Hamiltonian, which differs from that of the Z 2 gauge theory in that G j is squared. The reason is that G j of Eq. (2) of the Z 2 gauge theory has only two eigenvalues g j = 0, 2. However, G j of Eq. (S5) of the U(1) gauge theory has four eigenvalues g j = −1, 0, 1, 2, one of which is negative. As such j G j may be zero even though locally G j itself may not vanish. More generally, using j G j as a measure of the violation in Gauss's law would underestimate it in the U(1) gauge theory, and V j G j would not energetically penalize all violations of Gauss's law. Consequently, G j needs to be raised to an even power in Eq. (S7), where, for simplicity, we choose this power to be 2. Of course, one can alternatively use V H G = j |G j |.
We now prepare the system in the ground state of Eq. (S4) at J = 0, µ = 1, λ = 0, and V = 1, and then quench it with H 0 + λH 1 + V H G at J = 1, µ = 0, and some λ and V , and average the results over 1000 configurations of the integers x j , y j , z j and coefficients κ j , η j , γ j . In the following, however, we note that our conclusions remain the same even without this disorder averaging. In Fig. S5 we show results for λ = 0.05 (left-column panels) and λ = 0.5 (right- Figure S4. (Color online). Finite-size effects on the dynamics of the spatiotemporal averages of the gauge-invariance violation, electric field, and staggered boson number in the Z2 gauge theory for λ = 0.05 and V = 5, which is in the controlled-violation regime as shown in Fig. 4 of the main text. We see that in the controlled-violation regime larger system sizes indicate better protection of gauge invariance.
column panels) at various values of V . We consider the spatiotemporal averages of the gauge-invariance violation the magnetization in the z direction (the 'electric field') as well as the boson number where for the latter two we look at their deviation from the ideal gauge-invariant case. Here |ψ(t) = exp[−i(H 0 + λH 1 + V H G )] |ψ 0 . As in the case of the Z 2 gauge theory, the violation is suppressed by (λ/V ) 2 up to infinite times in the controlled-violation regime as shown in Fig. S5(a) for λ = 0.05 and Fig. S5(d) for λ = 0.5. Also identically to the Z 2 gauge theory, the gauge-invariant observables m z (t) and N (t) deviate from their ideal dynamics at a timescale ∝ V /λ 2 in the controlled-violation regime, as shown in Fig. S5(b,c) for λ = 0.05 and Fig. S5(e,f) for λ = 0.5. In conclusion, we see two distinct regimes: one of uncontrolled violation when V is too small to counter the effects of Eq. (S6), and a controlled-violation regime when V is sufficiently large where the gauge-invariance violation is suppressed by (λ/V ) 2 up to infinite times, as shown in the bottom panel of Fig. 4 in the main text. As such, our results are independent of whether the model is a Z 2 or U(1) gauge theory, and we expect our conclusions to be valid for more general lattice gauge theories.
(a) (d)  It is clear from this plot that the conclusions of the main text are unchanged whether looking at the spatially averaged gauge-invariance violation itself or its temporal average. Nevertheless, we opt to present the spatiotemporal average as that mitigates oscillations due to finite system size.

SPECIFIC CHOICE OF OBSERVABLES
We opt for temporal averages of our signals to reduce oscillation due to the small system sizes we are limited to in exact diagonalization. This does not alter the conclusions of our work, but merely presents them in a more eye-friendly manner. As an example, let us look at the suppression of the spatially (but not temporally) averaged gauge-invariance violation j G j (t) /L in the Z 2 gauge theory, as shown in Fig. S6. The qualitative picture remains unchanged when comparing j G j (t) /L to its temporal average in Eq. (4), shown as dashed lines in Fig. S6; cf. Fig. 3(a). Figure S7. (Color online). The ideal dynamics (λ = V = 0) of the spatially averaged electric field's raw value (blue), its temporal average (red), and mx(t), which is the temporal average of the absolute value (yellow). From the plot, it is straightforward to see that mx(t) not only offers a cleaner signal than the raw data, but it also captures dynamics -such as fluctuations -that cannot be captured by the mean of the raw data.
Another choice we make is to take the spatiotemporal average of the absolute value of the observable. The reason behind this is twofold. Firstly, this provides a way to also look at the fluctuations in the observable. For example, if we consider Fig. 2(b), we see a peculiar behavior in the ideal dynamics of m x (t) defined in Eq. (5). This picture is clarified when considering Fig. S7. The oscillations in the expectation value of the electric field, τ x j,j+1 (t) increases in amplitude with time, albeit its average remains zero. As such, m x (t), in capturing these fluctuations, offers a more stringent comparison when it comes to nonideal dynamics due to breaking of gauge invariance. Secondly, when gauge invariance is violated, even if slightly, this may lead to phase shifts that will never vanish even after adding a strong protection term that controllably suppresses the violation. By taking the spatiotemporal average of the absolute value of a quantity, such "superficial" differences are neglected.

SPECIFICATIONS WITH REGARDS TO THE NUMERICS
For our exact diagonalization simulations we use, in addition to our in-house code, the QuTiP [S7, S8] and QuSpin [S9, S10] toolkits in order to cross-check our results. For time evolution, we use our own exact exponentiation routine, because this does not require an iterative calculation that may hinder our ability to reach the long times we achieve. The solvers in QuTiP and QuSpin usually rely on solving an ordinary differential equation, which requires a small enough time-step in order to converge. This in turn demands a large computational time in order to achieve the long times we consider here.
Moreover, as can be seen in our results, the maximum size used for the U(1) gauge theory is L = 6 matter sites (plus 6 links), while for the Z 2 gauge theory the largest chain is L = 8 matter sites (plus 8 links). The reason is because in the Z 2 gauge theory, even after breaking gauge invariance according to Eq. (3), the particle number is conserved, thereby leading to a reduction in the effective Hilbert system in which the dynamics lives. In the case of the U(1) gauge theory, no such conservation laws are present, which therefore limits the largest L we can reach.

MATHEMATICALLY RIGOROUS BOUNDS ON EMERGENT GAUGE SYMMETRY AT V > 0
In this and the following section, we corroborate our numerical results by analytic arguments. Assume V is sufficiently strong to open an energy gap ∆ ∝ V from the gauge-invariant sector G j |ψ = 0, ∀j, to gauge-violating sectors. By adapting the results of Ref. [S11], we can then mathematically rigorously show that the gauge-invariant sector is protected by an emergent symmetry that is close to the original one.
Following Ref.
[S11], we define a ground symmetry as an operatorŨ that commutes with the ground space projector Π, [Ũ , Π] = 0, and acts unitarily on the ground space ΠŨ †Ũ Π = ΠŨŨ † Π = Π. In our case, Π is the projector onto the sector G j |ψ = 0, ∀j. Further, a unitary U is called an -approximate symmetry if it approximately commutes with the Hamiltonian with respect to a given unitarily invariant norm ||U, H|| ≤ . In our scenario of a gauge theory, U = exp i j ν j G j with ν j ∈ R, i.e., the relevant unitary is a product of O(L) local symmetries. The relevant non-commutativity is governed by the gauge-violating error strength, so ∝ λ.
In a generic many-body system, H 1 is extensive and consists of few-body operators, so its non-commutativity with U will give = O(λL). Moreover, typically the bandwidth of the ground manifold increases with L, which reduces the gap to the other symmetry sectors. Thus, from the above we may expect ||Π U −Ũ Π|| to deteriorate as a power (larger than 2) of L. As our numerics indicates, however, this requirement -though mathematically rigorous -is in practice too stringent, and a scaling with system size that is much more benign can be achieved.

PERTURBATION THEORY
In this section, we further explain our numerical results through analytic arguments based on perturbation theory, first for V = 0 and then for V being the dominant energy scale.

Time-dependent perturbation theory for V = 0
In the case of V = 0, our Hamiltonian is given by where [H 0 , G j ] = 0 and [H 1 , G j ] = 0, ∀j, with G j being a local gauge generator at position j. Consequently, the gauge invariance encapsulated in H 0 is broken up to a strength λ by the gauge-noninvariant term H 1 . We wish to see how the resulting violation grows in time using time-dependent perturbation theory in λ.
The total Hilbert space contains gauge-invariant sectors, where in each such sector Gauss's law takes on a unique value locally. If our initial state is in a given sector, H 0 propagates its dynamics only within that sector, while H 1 includes processes that drive the system into other sectors within the total Hilbert space of the system. Now since [H 0 , G j ] = 0 and [G j , G x ] = 0, we can find a common eigenbasis {|α, j } for H 0 and all G j , where α = (α 1 , α 2 , . . .) denotes the gauge-invariant sector defined by the unique set of local values α, and q stands for all remaining good quantum numbers. Thus, we have H 0 |α, q = E α,q |α, q and G j |α, q = α j |α, q . Our initial state at t = 0 is |ψ 0 , which resides in the sector 0, i.e., G j |ψ 0 = 0 ∀j. As such, our initial state can be written as where in the second equality we have indicated that |ψ 0 is prepared in the gauge-invariant sector 0.
where in the second equality we have utilized the gauge invariance of O. Another useful quantity to calculate, and which we will encounter later, is e i(Eα,q−E α,l )t α, q| O |α, l |α, q α, l| .
On the other hand, gauge-invariant observables that do not commute with H 0 do not necessarily scale as ∼ (λt) 2 at short times. For example, this can be seen in the Z 2 gauge theory where the deviation from ideal gauge-invariant dynamics of the spatiotemporally averaged staggered magnetization of Eq. (6) scales at short times as ∆N stag ∼ λt 2 ; cf. Fig. 2(c). This can be understood by noting that the zeroth-order contribution in Eq. (S20) from perturbation theory vanishes in ∆N stag , whereas in Eq. (S21) the term 0, m| j (−1) j a † j a j |0, n does not necessarily vanish because [H 0 , j (−1) j a † j a j ] = 0. As such, we can look at what happens to Eq. (S21) at short times for O = N = j (−1) j a † j a j by approximating it as We note that the only nonvanishing contribution from 0, n| H 1 |0, m is due to the terms whose coefficients are c 1 and c 2 in Eq. thereby explaining the short-time scaling ∼ λt 2 of the deviation in the staggered boson number from ideal dynamics in Fig. 2(c). As a summary, in case H 1 contains only terms that bring us completely out of the initial gauge-invariant sector, all gauge-invariant observables scale as λ 2 . For j G p j and other gauge-invariant observables that commute with H 0 , this scaling remains true even if H 1 contains also terms that generate dynamics within the initial gauge-invariant sector, but not necessarily for other gauge-invariant observables.