Faster Digital Quantum Simulation by Symmetry Protection

Simulating the dynamics of quantum systems is an important application of quantum computers and has seen a variety of implementations on current hardware. We show that by introducing quantum gates implementing unitary transformations generated by the symmetries of the system, one can induce destructive interference between the errors from different steps of the simulation, effectively giving faster quantum simulation by symmetry protection. We derive rigorous bounds on the error of a symmetry-protected simulation algorithm and identify conditions for optimal symmetry protection. In particular, when the symmetry transformations are chosen as powers of a unitary, the error of the algorithm is approximately projected to the so-called quantum Zeno subspaces. We prove a bound on this approximation error, exponentially improving a recent result of Burgarth, Facchi, Gramegna, and Pascazio. We apply our technique to the simulations of the XXZ Heisenberg interactions with local disorder and the Schwinger model in quantum field theory. For both systems, our algorithm can reduce the simulation error by several orders of magnitude over the unprotected simulation. Finally, we provide numerical evidence suggesting that our technique can also protect simulation against other types of coherent, temporally correlated errors, such as the $1/f$ noise commonly found in solid-state experiments.

Simulating the dynamics of quantum systems is an important application of quantum computers and has seen a variety of implementations on current hardware. We show that by introducing quantum gates implementing unitary transformations generated by the symmetries of the system, one can induce destructive interference between the errors from different steps of the simulation, effectively giving faster quantum simulation by symmetry protection. We derive rigorous bounds on the error of a symmetry-protected simulation algorithm and identify conditions for optimal symmetry protection. In particular, when the symmetry transformations are chosen as powers of a unitary, the error of the algorithm is approximately projected to the so-called quantum Zeno subspaces. We prove a bound on this approximation error, exponentially improving a recent result of Burgarth, Facchi, Gramegna, and Pascazio. We apply our technique to the simulations of the XXZ Heisenberg interactions with local disorder and the Schwinger model in quantum field theory. For both systems, our algorithm can reduce the simulation error by several orders of magnitude over the unprotected simulation. Finally, we provide numerical evidence suggesting that our technique can also protect simulation against other types of coherent, temporally correlated errors, such as the 1/f noise commonly found in solid-state experiments.

I. INTRODUCTION
Simulating the dynamics of quantum systems is a key application of quantum computers. However, digitalizing the continuous time evolutions to enable execution on gate-based and other programmable quantum computers comes with simulation errors that cause the dynamics of the systems to deviate from ideal evolutions. In particular, the errors may violate the symmetries in the target Hamiltonian for simulation, resulting in unphysical states at the end of the simulations. This digitalization error particularly affects Trotterization-the most common algorithm for near-term quantum simulations [1][2][3]-and persists even in more sophisticated, advanced quantum simulation algorithms [4][5][6].
In this paper, we propose an approach, using the symmetries of target Hamiltonians, to protect the dynamics of the systems against simulation errors. Given a simulation algorithm that decomposes the dynamics of the system into many small time steps (e.g., Trotterization), we interweave the simulations with unitary transformations generated by the symmetries of the systems (Fig. 1). While these additional unitary transformations increase the gate complexity of the simulation, the error of the simulation can sometimes be reduced by several orders of magnitude, ultimately resulting in a faster quantum simulation. In addition, depending on the symmetries, the unitary transformations may be implemented using only single-qubit gates, which are considered relatively inexpensive for implementations on near-term quantum computers.
Our technique is general and potentially applies to any algorithms that simulate the time evolution of Hamiltonians with symmetries by splitting the evolution into many time segments, including Trotterization and the higher-order product formulas [4] and more advanced algorithms such as those based on linear combinations of unitaries [5,6], Lieb-Robinson bounds [7,8], and randomized compilations [9,10]. We also provide evi- Figure 1. For algorithms that simulate the dynamics of quantum systems by decomposing the evolutions into many time steps, we interweave the corresponding simulation circuits (blue) with unitary transformations generated by the symmetries of the systems (orange). These transformations protect the simulations against errors that violate the symmetries, resulting in faster and more accurate simulations. dence that our technique can also protect the simulation against other types of temporally correlated errors, such as the 1/f noise commonly found in solid-state devices [11].
In addition, we draw a connection between our symmetry protection technique and the quantum Zeno effect [12][13][14][15][16][17][18]. In particular, the symmetry transformations, when chosen as powers of a unitary, approximately project the error of simulation into the so-called quantum Zeno subspaces, defined by the eigensubspaces of the unitary. We prove a bound on the accuracy of this approximation, exponentially improving a recent result of Ref. [18].
The structure of the paper is as follows. In Section II, we introduce the general technique and provides intuition for the source of error reduction. In Section III, we derive a bound on the error of Trotterization under symmetry protection. In Section IV, we then benchmark our technique in simulating the dynamics of systems with the Heisenberg interactions, including the XXZ Heisenberg model with local disorder that displays a transition between thermalized and many-body localized phases, and in simulating the Schwinger model in the context of lattice field theories. In particular, we show that interweaving the simulation with random gauge transformations can significantly reduce the probability of a state leaking to outside the physical subspace due to the simulation error, extending the results of Ref. [19] to digital quantum simulation. We then demonstrate in Section V how our technique may protect the simulation against other types of coherent, temporally correlated errors, such as the low-frequency noise typically found in experiments. Finally, we discuss several open questions in Section VI.

II. GENERAL FRAMEWORK
We consider the task of simulating the time dynamics of a system under a time-independent Hamiltonian H. Let U t ≡ exp(−iHt) denote the evolution unitary generated by H for time t. Our technique applies to algorithms that simulate U t by first dividing the evolution into many time steps (also known as Trotter steps), and approximate the evolution within each time step by a series of quantum gates. In this paper, we focus on the first-order Trotterization algorithm for simplicity. To be more precise, let r denote the number of steps and δt = t/r denote the length of each time step. These algorithms then simulate U δt by a series of elementary quantum gates S δt , i.e.
The approximation of U δt by S δt introduces an error that is small for small δt. However, errors typically accumulate after many Trotter steps, resulting in a total additive error U t − S r δt that, in the worse case, scales linearly with the number of Trotter steps r at fixed δt. Equivalently, for a fixed total time t, to reduce the total error, we would have to decrease the Trotter step size δt, effectively increasing the number of Trotter steps r, and thus require more elementary quantum gates to run the simulation.
We refer to the simulation in Eq. (1) as the raw simulation. By exploiting symmetries of the system, we will see that we can substantially reduce the total error ε of the simulation without significantly increasing the gate count, ultimately resulting in faster quantum simulation for the same total error budget. For that, we assume that the Hamiltonian is invariant under a group of unitary transformations, which we denote by S. Explicitly, we assume that The group S represents a symmetry of the system. Instead of simply approximating U δt by the circuit S δt , we "rotate" each implementation of S δt by a symmetry transformation C k ∈ S (i = k, . . . , r) so that the approximation in Eq. (1) now reads We refer to Eq. (3) as a symmetry-protected (SP) simulation. The right-hand side in Eq. (3) represents a circuit that, at first, looks more expensive than Eq. (1) due to the additional implementation of the transformations C k . However, for the same r, the total error in Eq. (3) could be much smaller than the Eq. (1). Effectively, to meet the same error tolerance, Eq. (3) may require a much smaller number of steps r, and hence fewer implementations of S δt , than the raw approximation in Eq. (1). Moreover, because many symmetries-the gauge symmetries in lattice field theories for exampleare spatially local, each C k only involves a small number of nearest-neighboring qubits and can be implemented easily in most architectures of quantum computers. Other symmetries, such as the one responsible for the conservation of the total magnetization in the Heisenberg model, are global but may be implemented as a product of only single-qubit gates, which are usually much "cheaper" to perform in experiments than their multi-qubit counterparts.
In the remainder of this section, we provide some intuition, using lowest-order arguments, for the error reduction in simulations under symmetry protection. We later derive rigorous error bounds in Section III.

A. Lowest-order arguments
To build an intuition for the symmetry protection, we consider the effective Hamiltonian of the simulation. The aim of digital quantum simulation is to simulate the time evolution e −iHt of a Hamiltonian H. Assuming that the simulation errors are coherent, we may end up with the time evolution of a different Hamiltonian, say H eff , that may be close but not the same as the targeted Hamiltonian H: where quantifies the difference between the effective and the desired Hamiltonians [20]. We note that the effective Hamiltonian H eff typically depends on the time step δt [See Lemma 1]. With S δt = exp(−iH eff δt) in Eq. (3), we can rewrite the simulation as where we have used the unitarity of C k to move the unitaries to the exponents and exploited the commutativity [C k , H] = 0 from our assumption to simplify the expression. Assuming that the error V is small, we can use the Baker-Campbell-Hausdorff (BCH) formula to combine the exponents in Eq. (6) (to the leading order): Compared to the desired evolution e −iHt , we can identify the error of the entire simulation (ignoring the error from the BCH approximation for now) as Roughly speaking, the error of the entire simulation, given by Eq. (8), can be interpreted as the average of the error in each step of the simulation. To illustrate the effect of the symmetry protection, we could imagine V as a vector in the space of operators and C † k V C k is a version of the vector rotated around an axis specified by C k . The total error is then analogous to a walker that, in each step, walks a distance V in the space of operators towards the direction corresponding to C k (Fig. 2).
Without the symmetry protection (i.e. C k = I for all k), the walker keeps walking in the same direction and its total distance after r steps scales as O (r), resulting in the averaged error V of the same order as V . On the other hand, under the symmetry protection, the walker walks in a possibly different direction in each step, resulting in a smaller total distance (and thus a smaller averaged error.) In particular, if the walker in each step walks towards a uniformly random direction in the space of operators Figure 2. The total error of the simulation is analogous to the average distance a walker walks in r steps of the simulation. In each time step, the walker walks a small distance along a vector representing the error operator in the space of operators. a) Without any symmetry protection, the walker keeps walking towards almost the same direction, resulting in a total distance that scales linearly with the number of steps r, corresponding to the total error scaling as O (1). b) The symmetry transformations make the walker walk in a possibly different direction in every time step. When the direction is uniformly random (see Section IV A 1 and Fig. 3 for an example), the total distance only scales as O ( √ r), resulting in the total error scaling as O (1/ √ r). c) Sometimes, it is possible to design an optimal set of symmetry transformations that makes the walker return to the origin [See Eq. (38) for an example], resulting in an O (1/r) error for the entire simulation.
(which is sometimes the result of choosing C k at random), its total distance should only scale as O ( √ r V ) after r steps. The averaged error V would then scale as O ( V / √ r), decreasing with the number of steps of the simulation. Additionally, if we could design a set of optimal symmetry transformations that makes the walker return to the origin after a fixed number of steps, we would end up with a total distance that does not increase with r and an averaged error V that decreases with r as O ( V /r). We derive rigorous bounds to support this intuition in Section III.
The aim of our technique is to choose the symmetry transformations C k that minimize the error in Eq. (8). While each C k may be chosen independently of the others, we will sometimes focus our attention on a special construction that requires C k = C k 0 for some C 0 ∈ S. This choice for the transformations result in a simpler simulation circuit, i.e.
which corresponds to applying the same symmetry transformation C 0 alternatively with the implementations of the simulating circuit S δt , followed by a final application of C †r 0 to negate the effect of C 0 on the correct evolution. We could either draw C 0 randomly from the symmetry group S or infer an optimal choice of C 0 from the struc-ture of the error V [See Eq. (38) for an example]. We analyze the error bounds for the simulation under the protection from this special construction in Section III and present similar analysis for the general scenario in Appendix C.
It is worth noting that the symmetry transformation C 0 introduced above is also analogous to the fast pulses (or "kicks") commonly used in quantum control to confine the dynamics of quantum systems [12][13][14][15][16][17][18]. In fact, we also show in Appendices A and B that a restricted version of our symmetry protection technique is exactly equivalent to frequently applying fast pulses to the systems, resulting in the error being approximately projected onto the so-called quantum Zeno subspaces. We prove a bound on the error of this approximation, exponentially improving a recent result of Ref. [18]. This quantum Zeno framework provides an alternative explanation for how quantum simulation can be improved by symmetry protection.

III. FASTER TROTTERIZATION BY SYMMETRY PROTECTION
In this section, we analyze the effect of the symmetry protection on the total error of the first-order Trotterization algorithm. Suppose the Hamiltonian H = L µ=1 H µ is a sum of L Hamiltonian terms H µ such that each e −iHµδt can be readily simulated on quantum computers. For readability, we define the following quantities that depend only on the commutators between the terms of the Hamiltonian. We will also use the standard Bachmann-Landau big-O and big-Θ notations in analyzing the asymptotic scalings of the errors with respect to n, t, and r. For reference, α = O (n) and β = O (n) in a system of n nearest-neighbor interacting particles [20]. Given a set of symmetry transformations C = {C k : k = 1, . . . , r}, we define as the version of an operator A averaged over the rotations induced by C k . The first-order Trotterization algorithm approximates exp(−iHδt) by where L µ=1 U µ ≡ U L . . . U 2 U 1 is an ordered product. We define H eff as the generator of S δt , i.e. S δt = exp(−iH eff δt). We prove the following lemma, providing the existence and the structure of the generator H eff .

Lemma 1.
For all δt such that βδt ≤ α, 2αδt ≤ H , and 8δt H ≤ 1, there exists a generator H eff for S δt and where V(δt) is an operator satisfying V(δt) ≤ χδt 2 and We provide the proof of Lemma 1 in Appendix D. The essence of Lemma 1 is that the error of the simulation, defined as V ≡ H eff − H, is given by and it follows that V ≤ 1 2 αδt + χδt 2 . We now consider the effect of protecting the simulation with a set of symmetry transformations {C k : k = 1, . . . , r}. Under this symmetry protection, each circuit S δt is replaced by where we have used [C k , H] = 0 to simplify the expression. The full simulation becomes In the following analysis, we further assume that the symmetry transformations C k have the form C k = C k 0 , where C 0 is a symmetry transformation drawn from the symmetry group S (We extend these results to general symmetry transformations in Appendix C.) Let {e −iφµ : 1 ≤ µ ≤ m} denote the distinct eigenvalues of C 0 and where is the inverse spectral gap that depends on the eigenvalues of C 0 .
The proof of Lemma 2 follows from Lemma 5 in Appendix B. We note that the bound in Lemma 2 depends on m, the number of unique eigenvalues of C 0 , which could be a constant, e.g. when C 0 is generated by local symmetries, or depend on the system size, e.g. when C 0 corresponds to generic rotations generated by global symmetries. We also note that the inverse spectral gap ξ could be large if C 0 is nearly degenerate and one should take this effect into account when choosing the unitary C 0 .
Lemma 2 says that, up to the error given in Eq. (21), the simulation under the symmetry protection is effectively described by H eff . In particular, the total error of the Hamiltonian under the symmetry protection is where we have replaced the expression of V from Lemma 1. Note that V ≤ V by the triangle inequality. Using the identity we arrive at the following bound on the total error of the simulation.
Theorem 1 (Quantum simulation by symmetry protection). Assuming that βδt ≤ α, 2αδt ≤ H , and 8δt H ≤ 1, the total error of simulation under the symmetry protection from {C k = C k 0 : C 0 ∈ S, k = 1, . . . , r} can be bounded as where m is the number of distinct eigenvalues of C 0 , and ξ is the inverse spectral gap defined in Eq. (22).
The proof of Theorem 1 follows immediately from Lemma 2 and Eq. (25) [See Appendix E for the detailed calculations]. The key feature of Theorem 1 is that, to the lowest-order in t r , the error scales is generally smaller than v 0 when [C k , v 0 ] = 0, we expect a smaller simulation error under the symmetry protection.
For demonstration, we consider the simulation of a Hamiltonian H that is a sum of nearest-neighbor interactions on n particles. It is straightforward to verify that for this Hamiltonian, We will also assume that the number of distinct eigenvalues of the C 0 is m = O (1) (corresponding to local symmetries or highly degenerate transformations) which results in κ = O n 2 . We will estimate the required number of steps r-a good proxy for the gate count [21]-for simulations with and without the symmetry protection.
The first scenario corresponds to an unprotected simulation, where v 0 = v 0 . The total error then scales as To meet a fixed error tolerance ε, we would have to choose the number of steps r = Θ(nt 2 /ε).
On the other hand, with symmetry protection, we later show that it is sometimes possible to make v 0 vanish completely, making the higher order terms the dominant contribution to the total error [See Eq. (38) for an example]. For nearest-neighbor interactions, the total error is now which decreases quadratically with r. As a result, we only need whereΘ(·) is Θ(·) up to a logarithmic correction. Note that this choice of r also satisfies the conditions in Theorem 1 when t/ε > 1. Compared to the unprotected simulation, the symmetry protection results in a factor of t/ε improvement in the required number of steps. At ε = 0.01, the improvement in the scaling with ε alone would result in about a factor of ten reduction in the gate count of the simulation. Finally, we consider a scenario where v 0 ∝ v 0 /r γ for some γ ∈ (0, 1). We provide an example of such a scaling in Section IV A 1, where drawing the unitary transformations C k randomly from the symmetry group results in a scaling with γ = 0.5. This scaling of v 0 results in the total error Hence, we require which is again better than the unprotected simulation by a factor of min{(nt 2 /ε) γ/(1+γ) , t/ε}.
We recall that in deriving Theorem 1, we have assumed that the symmetry transformations have the form C k = C k 0 for some C 0 . We derive in Appendix C a different bound for the general case where each C k may be chosen independently. This general bound, while appearing more complicated, holds the same key feature to the bound in Theorem 1: the total error, to the lowestorder, scales with an averaged version of v 0 (under the symmetry transformations) instead of scaling with v 0 .

IV. APPLICATIONS
In this section, we apply our technique to the simulation of the Heisenberg model (Section IV A) and lattice field theories (Section IV B). In both cases, we show that the symmetry protection results in a significant error reduction and thereby gives faster quantum simulation.
In particular, we use the simulation of the homogeneous Heisenberg model in Section IV A 1 to demonstrate the improvement on the total error scaling as a function of the number of steps r when the simulation is protected by a random set of unitary transformations and by an optimally chosen set. In Section IV A 2, we estimate the required number of Trotter steps as a proxy for the gate count in simulating an instance of the Heisenberg model, commonly found in the studies of the many-body localization phenomenon. Finally, in Section IV B, we consider the probability of the state leaking to unphysical subspaces in the digital simulation of the Schwinger model and show that the symmetry protection from the local gauge symmetries can suppress this leakage by a few orders of magnitude.

A. Heisenberg interactions
In this section, we use the symmetries in the Heisenberg model to protect its simulation using the first-order Trotterization. A Heisenberg model of n spins can be described by the Hamiltonian represent the interaction strengths between the spins, and h i correspond to the strengths at site i of an external magnetic field pointing in the z direction. The Heisenberg model provides a good description for the behavior of magnetic materials in the presence of external magnetic fields. Depending on several factors, including the signs of the interactions and the dimensions of the system, the Heisenberg model may undergo a quantum phase transition as we increase the strength of the external magnetic field. Several important instances of the Heisenberg model includes the homogeneous Heisenberg model (J (x) = J (y) = J (z) ), the XXZ model (J (x) = J (y) ) with local disorder, and the Ising model (J (y) = J (z) = 0). In the following subsections, we will consider two pedagogical instances of Eq. (33) with SU(2) and U(1) symmetries respectively and demonstrate how our technique helps in simulating the dynamics of these systems even as they move across critical points.

Homogeneous, random Heisenberg interactions
We first consider a pedagogical toy model where interactions in Eq. (33) are homogeneous, i.e. J 1]. In addition, we assume that h i = 0 ∀i, i.e. there is no external magnetic field. In this case, Eq. (33) simplifies to The combination of homogeneous interactions and no external magnetic field make Eq. (34) invariant under S = {W ⊗n : W ∈ SU(2)}, which contains unitaries that-in the Bloch sphere-simultaneously rotate each spin by the same angle.
To simulate the evolution U t under Eq. (34), we could use the first-order Trotterization to approximate by a product of evolutions of individual terms of the Hamiltonian. The number of Trotter steps r and the time step δt = t/r determine the error of the simulation. We refer to this approach as the raw Trotterization. To protect this simulation, we insert unitaries drawn from the symmetry group S in between the Trotter steps, resulting in the simulation where {C 1 , . . . , C r } ≡ C is a subset of the symmetry group S. Recall that the total error of this symmetryprotected simulation is given by Theorem 1, with the lowest-order error being Figure 3. The total error in simulating the Hamiltonian Eq. (34) at n = 4 for a fixed evolution time t = 1 as a function of the Trotter number r using four different schemes: the raw first-order Trotterization ("Raw"), the first-order Trotterization protected by a random set symmetry transformation ("SP-Rand"), the first-order Trotterization protected by the optimal set in Eq. (38) ("SP-Det"), and the random-ordering scheme in Ref. [10] ("Random Ordering"). We indicate the scalings obtained from power-law fits to the right of the plot. We repeat the simulation 100 times, each with a different set of randomly generated interactions Jij. The dots correspond to the median of the errors at each value of r and the bars represent the corresponding 25%-75% percentiles regions. where comes from the leading contribution to the error in one Trotter step. Different choices of the set C lead to different total error of the simulation. For minimal calculational overhead, we could choose each C k independently and uniformly at random from S (i.e. C k = W ⊗n k where W k is a Haar random unitary on the single-qubit Bloch sphere.) The sum in Eq. (37) is then the sum of v 0 , each rotated under a random unitary. This is analogous to the total error being a random walker that, in each time step, "walks" a distance v 0 in a random direction (See Fig. 2). From this analogy, we then expect v 0 ∝ v 0 / √ r (to the lowest-order). Therefore, we expect the total error of this scheme to decrease with r as O r −3/2 (at fixed total time t).
While randomly choosing the unitary transformation set C requires little to no knowledge about the error operator v 0 , one can expect that this choice of C is not optimal. Indeed, by further exploiting the structure of v 0 , we can construct a set of transformations C that makes Eq. (37) vanishes entirely. One such choice is and U H is the single-qubit Hadamard matrix. Alternatively, we could also write for k = 1, . . . , r. Since the Hadamard matrix switches X ↔ Z and Y ↔ −Y , it is straightforward to verify that Eq. (37) vanishes for all even values of r. Therefore, the total error of the simulation is given by the next lowest order in Theorem 1, which scales with r as O 1/r 2 . In Fig. 3, we plot the total error of the simulation at n = 4, t = 1 as a function of the Trotter number r for the three aforementioned scenarios: the first-order Trotterization without symmetry protection ("Raw"), with symmetry protection from a randomly chosen C ("SP-Rand"), and with symmetry protection from the optimal set C ("SP-Det"). The scalings of the errors as functions of r agree remarkably well with our above prediction. In addition, we also compute the total error using the randomized simulation scheme in Ref. [10], which decreases the Trotter error by randomizing the ordering of the Hamiltonian terms in between Trotter steps. Our numerics shows that this scheme performs similarly to the simulation protected by random symmetry transformations, which are both outperformed by the optimal symmetry protection scheme.

Many-body localization
The homogeneous Heisenberg interactions without external fields considered in the previous section provides a good testbed for benchmarking our technique. In this section, we consider a more physically relevant instance of the Heisenberg model: where we again assume homogeneity for the coupling strengths, but J ij = 1 only when i, j are nearest neighbors and J ij = 0 otherwise. We also adopt the periodic boundary condition and identify the (n + 1)th qubit as the first qubit. In addition, we add an external magnetic field with the field strength h i , each chosen randomly between [−h, h]. This model describes homogeneous Heisenberg interactions with a tunable local disorder strength h. At low disorder h, the system evolved under Eq. (40) thermalizes in the long-time limit, in agreement with the Eigenstate Thermalization Hypothesis (ETH). However, as h increases, the system transitions to a many-body localized (MBL) phase where it no longer thermalizes (See [22] for a review of the manybody localization phenomenon.) To simulate the dynamics of H, we again divide the terms of H into groups of mutually commuting terms: and use the first-order Trotterization similarly to Eq. (35). To symmetry-protect this simulation, we note that the field term breaks the SU(2) symmetry of the Heisenberg interactions, leaving the system invariant While selecting the unitary transformations C k from this U(1) symmetry is no longer sufficient to completely eliminate the lowest-order error-as we have done in the previous section-we can still expect significantly reduction of the total error due to the symmetry protection and thus a lower gate count for the simulation. In Fig. 4, we plot the number of Trotter steps r in simulating the dynamics of Eq. (40) for time t = n at different values of the disorder h that correspond to the ETH and the MBL phases. The required numbers of steps are computed at each n by binary searching for the minimum r such that the total error of the simulation does not exceed ε = 0.01. Figure 4 shows that protecting the simulation with the U(1) symmetry results in several times reduction in the number of Trotter steps for all values of n. In addition, the Trotter number under symmetry protection also appears to scale better with the system size than in the raw simulation, suggesting an even greater advantage from the symmetry protection for simulating larger systems.
Out of curiosity, we study how the symmetry protection performs as the Hamiltonian moves across the ETH-MBL phase transition. In Fig. 5, we plot the required number of steps r in simulating the Hamiltonian of n = 8 qubits for time t = n and error tolerance ε = 0.01 as we tune the Hamiltonian from the ETH to the MBL phase [23]. The improvement due to the symmetry protection appears to be unaffected by the phase transition, suggesting that our technique can be useful for future numerical and experimental studies of the transition.

B. Simulation of lattice gauge field theories
Quantum field theories provide another key target for quantum simulation [24]. In particular, the quantum simulation of real-time Hamiltonian dynamics, for example scattering processes [25], has attracted much attention. An important class of field theories are models with local gauge symmetry, including quantum electrodynamics, chromodynamics, and the Standard Model of particle physics in addition to many condensed matter systems. Substantial effort has gone into the study of analog [26][27][28] and digital [29][30][31] quantum simulation of these models.
In a gauge theory, the system is invariant under a symmetry group which acts separately at each point in space and time (see eg. [32] for a review, as well as the lattice Hamiltonian formulation, of these models). This symmetry is fundamentally a redundancy of our description of the physics which we have introduced to give a local description. The Hilbert space H we use to describe the system contains a subspace H phys of the physical states, those annihilated by the gauge constraints. For example, in electrodynamics, we have the charge and gauge field degrees of freedom, and the physical states are those annihilated by the Gauss law constraint G = ∇·E−ρ, where E is the electric field operator and ρ is the charge density operator. There are many states in the full Hilbert space H which do not live in the kernel of G, and these states are not allowed in nature. Although one can in principle work with a description strictly within the physical Hilbert space, it is in general computationally difficult to do the reduction. More importantly, this description would necessarily have a highly spatially non-local set of interactions, a major drawback in practice.
Thus in the simulation of a gauge theory we are faced with a fundamental source of possible errors: what if our dynamics takes us away from the physical Hilbert space? Although the exact Hamiltonian commutes with the gauge constraints, and thus leaves the physical space invariant, an approximate (for example, Trotterized) version of the Hamiltonian will typically induce leakage into the unphysical space. In this section, we apply our symmetry protection technique and use the gauge symmetry itself to protect the simulation against this undesirable leakage [33].
Explicitly, we consider the one-dimensional Schwinger model consisting of n sites and n − 1 nearest-neighbor links between the sites. We use the formalism outlined in Ref. [34]. The Hamiltonian H = H 0 + H 1 consists of two terms: where Here, H 0 describes the on-site and on-link terms, H 1 describes the site-link interaction, and F i is the electromagnetic field operator for the link that connects the ith and (i + 1)th particles. In a simulation, we have to put a cutoff Λ specifying the maximum excitation number for the bosonic degree of freedom on a given link.
The Hamiltonian is subjected to local symmetries generated by the gauge operators: where Q i = 1 2 −Z i + (−1) i counts the electric charge at site i. In particular, only states |ψ that satisfy G i = 0 for all i are considered physical.
The physical states form a subspace H phys which can be constructed from the kernels of the gauge operators: where Ker(G i ) = {|φ : G i |φ = 0} is the kernel of G i . Due to various errors, an initially physical state may leak to unphysical subspace during the simulation. Formally, we define the leakage of a state |ψ(t) at time t as where P H phys is the projector onto the physical subspace H phys .
To simulate e −iHδt for a small time δt, we first decompose it into e −iH0δt e −iH1δt using the first order Trotterization. Since both H 0 , H 1 commute with G i , this decomposition respects the gauge symmetries and does not result in leakage from the physical subspace. However, to simulate the evolution under H 1 , we need to further decompose it into elementary quantum gates. For that, we follow the steps in Ref. [34] and write where This representation allows us to decompose the evolution into a product of three-qubit gates that can be readily implemented on quantum computers [34].
Note that the cost of simulating e − 1 4 ixtÃiXiXi+1 is that of simulating e − 1 4 ixtAiXiXi+1 , plus the cost of simulating The entire raw first-order Trotterization simulation of e −iHt becomes Similarly to the Heisenberg model, we could protect this simulation by interweaving the Trotter steps with symmetry transformations of the system: where C k are generated by the gauge operators in Eq. (46). Specifically, we choose C k = C k 0 , where and the angles φ i are independently and uniformly chosen at random from [0, 2π]. In Fig. 6, we plot the leakage outside the physical subspace due to the Trotter error. We use the Schwinger model with 4 sites and 3 links (all initialized in state |0 ) with and without symmetry protection. We repeat the simulation 100 times, each with a different choice of the transformation angles φ i . The figure shows that the symmetry protection can reduce the leakage by several orders of magnitude. In addition, while the leakage builds up in a raw simulation, it appears to be bounded during the course of the symmetry-protected simulation.

V. ADDITIONAL PROTECTION AGAINST EXPERIMENTAL ERRORS
So far, we have demonstrated that symmetries in quantum systems can be used to suppress the simulation error of the Trotterization algorithm. In this section, we discuss how our technique may also protect the simulation against other types of error, including the experimental errors that may arise in the implementation of Trotterization.
In our earlier derivation, we show that the lowest-order contribution to the total error is where v 0 is the lowest-order error from the simulation algorithm. This derivation applies equally well for the case when the error v 0 comes from sources other than the approximations in the simulation algorithms. However, in our analysis, we require that v 0 remains the same for different steps of the simulation. In other words, the error v 0 for different Trotter steps are correlated in time. In particular, an error with temporal correlation lengths being longer than the time step δt would enable us to choose the symmetry transformations such that the errors from several consecutive steps interfere destructively. Therefore, we expect our technique to help reduce low-frequency noises, such as the 1/f noise typically found in solid-state qubit systems.
We provide numerical evidence for this argument by adding temporally correlated errors to the simulation of the Schwinger model. Specifically, after each step k of the simulation, we apply single-qubit rotations exp(−iη σ ·n k ) on the system, where η = 0.01 is a small angle, around a random axisn k . These rotations mimic the effect of a depolarizing channel and violate the gauge symmetries, resulting in the state leaking to the unphysical subspace. To impart temporal correlations into this noise model, we choose the random unit vectorsn k again only after every λ consecutive Trotter steps. The param-   eter λ therefore plays the role of the correlation length of the noise.
In Fig. 7, we plot the probability that the state leaks to unphysical subspace (due to the simulation error) as a function of time for several values of the correlation length λ. To study the effect of our symmetry protection scheme on the added experimental noise, we use the fourth-order Trotterization in the simulation to suppress the algorithm error, making the added noise the main contributor to the leakage observed in Fig. 7. As expected, at λ = 1, the experimental error varies too fast between Trotter steps and is immune against our technique. However, our symmetry protection scheme begins to suppress the experimental error as soon as the noise becomes temporally correlated (λ > 1) and becomes more effective as the correlation length λ increases. Even at λ = 4, we have managed to reduce error by about an of magnitude.

VI. DISCUSSION & OUTLOOK
In this paper, we propose a general technique to suppress the error of quantum simulation using the symmetries available in quantum systems, ultimately resulting in faster digital quantum simulation. We have analyzed the technique when applied to the Trotterization algorithm and derived bounds on the total error of the simulation under symmetry protection. The bound provides insights for choosing the set of unitary transformations that optimally suppress the simulation error. We then benchmarked our technique in simulating the Heisenberg model and lattice field theories. Both examples showed that the symmetry protection results in significant reduction in the total error, and thus the gate count, of the simulation. Finally, we argue that our technique can also protect digital quantum simulation against temporally correlated noise in experiments.
An immediate open question is how well our technique can protect other, more advanced quantum simulation algorithms, such as the higher-order Suzuki-Trotter formulas [4], the truncated Taylor series [5], or qubitization [6]. The error structures of those algorithms are typically more complicated than the first-order Trotterization analyzed in this paper. Therefore, it is more difficult to infer the set of symmetry transformations that optimally protects the simulation. Nevertheless, extensive analytical and numerical studies of the effectiveness of our technique for protecting these advanced algorithms, especially when applied to the simulations of various physically relevant systems, such as the lattice field theories [29][30][31], or the electronic structures [35][36][37][38], would be useful for the long-term development of digital quantum simulation.
When the error structure of the algorithm is not readily available, an alternative promising approach for optimizing the set of symmetry transformations is to parameterize the transformations, variationally minimize the error of the first few simulation steps, and apply the same set of transformations repeatedly for the rest of the simulation [39]. Understanding when such a variational approach can suppress the error in a long simulation could provide a path towards a scalable symmetry protection with a minimal calculation overhead.
In addition, our analysis in this paper focuses primarily on the error of the simulation algorithm under the symmetry protection in the full Hilbert space. It would be interesting to, for example, build upon the recent result of Ref. [40] and analyze the symmetry-protected simulation error in a low-energy subspace.
Lastly, we would like to note that, although our analysis focuses on digital quantum simulation, we expect the symmetry protection technique to apply equally well for analog quantum simulation and classical simulation of the dynamics of quantum systems.  . The frequent kicks confine the dynamics of the system (solid arrows) to the so-called quantum Zeno subspaces, defined by the projectors Pµ in the spectral decomposition of the kicks U kick = µ e −iφµ Pµ. In particular, the kicks suppress the probability for the system to travel between the subspaces (dashed arrow). By generating the kicks from the symmetries of the system, we can target the simulation error-the sole contributor to possible violations of the symmetries in an ideal simulation-for suppression.
Theorem 2 (Faster convergence of quantum Zeno effect). Let U kick be the unitary defined in Eq. (A2) with m distinct eigenvalues, inverse spectral gap ξ, and a set of orthogonal projectors {P µ }. Let G Zeno = µ P µ GP µ denote the projection of a Hamiltonian G onto the subspaces defined by {P µ }. We have To prove Theorem 2, we rewrite the evolution as where we have defined Letting G [1,r] ≡ G 1 + · · · + G r , the first step of our proof is to establish the error bound This is the spectral-norm error of the first-order Trotter formula [20]. However, a naive error analysis in terms of the commutators between G j (see [20,Proposition 15] for example) gives a bound that does not decrease with r and thus fails to establish the desirable bound. Instead, we seek a better analysis that exploits the spectral information of U kick [18]. The starting point of our analysis is the established von Neumann's ergodic theorem whose proof is included for completeness.
Theorem 3 (Von Neumann's ergodic theorem). Let U be a unitary operator and U = m µ=1 e −iφµ P µ be its spectral decomposition, with φ 1 = 0 and φ µ distinct. Then, where Proof. The bound follows from We note that the condition φ 1 = 0 is not restrictive as we can always make φ 1 = 0 by adding a global phase to U kick [14]. Corollary 1. Let U be a unitary operator and U = m µ=1 e −iφµ P µ be its spectral decomposition. Then, for any operator G, where Proof. The claimed bound follows from where the first inequality follows from the bound As aforementioned, a naive analysis of the Trotter error fails to provide the desirable bound for quantum Zeno effect. Instead, we use a recursive approach to estimate the Trotter error Eq. (A11).
k=k0 G k for k 0 ≤ k 1 . For any s ≥ 1 and δt, we have s k=1 e −iG k δt − e −iG [1,s] Note that at s = r and δt = t/r, Lemma 3 reduces to Eq. (A11). We prove Lemma 3 by induction on s. Suppose Lemma 3 holds for s = s 1 and s = s 2 such that |s 2 − s 1 | ≤ 1, we shall prove that it holds for s = s 1 + s 2 . Using the triangle inequality [1,s 1 +s 2 ] δt − e −iG [s 1 +1,s 1 +s 2 ] δt e −iG [1,s 1 ] where we have used the inductive hypothesis and the Trotter error bound [20,Eq. (143)] in the last inequality. To bound the commutator norm, we use the following lemma.
where we have used Corollary 1 to prove the second inequality. Therefore, the lemma follows.
We now apply the above equation repeatedly to prove Lemma 3. Note that Lemma 3 holds trivially for s = 1. Suppose that it holds for all s ≤ s 0 for some s 0 ≥ 1. We shall prove that it holds for s = s 0 + 1. First, we consider the case where s is even, i.e. there exists an integer l ≥ 1 such that s = 2l. Applying Eq. (A23) with s 1 = s 2 = l, we get s k=1 e −iG k δt − e −iG [1,s] Therefore, Lemma 3 holds if s is even. When s is odd, there exists an integer l ≥ 1 such that s = 2l + 1. Applying Lemma 3 with s 1 = l and s 2 = l + 1, we have s k=1 e −iG k δt − e −iG [1,s]δt ≤ (2l log 2 l + 2(l + 1) log 2 (l + 1) + 2l + 1)ξ √ m G 2 δt 2 .

Appendix B: Symmetry protection by quantum Zeno effect
In this section, we make a formal connection between our symmetry protection technique and the quantum Zeno effect. In particular, we show how the quantum Zeno framework provides an alternative explanation for the suppression of simulation error under symmetry protection.
We first note that the symmetry transformations in our scheme are analogous to the kicks in the quantum Zeno framework. Suppose that the symmetry transformations have the form C k = C k 0 , where C 0 ∈ S is also a symmetry transformation. Let be the spectral decomposition of C 0 , with e −iφµ being the distinct eigenvalues and P µ being the projectors onto the respective eigensubspaces. The condition on e −iφµ being distinct ensures that C 0 satisfies the definition of U kick in Eq. (A2). With e −iHδt being approximated by a circuit S δt in each time step, our symmetry-protected simulation becomes where H eff is the generator of S δt and exists for a small enough δt (see Lemma 1). Comparing Eq. (B2) with Eq. (A3), we identify C 0 = U kick . Therefore, by Theorem 2, the symmetry protected simulation is effectively described by in the large r limit, where H eff,Zeno = µ P µ H eff P µ .
Recall that H eff is the effective Hamiltonian corresponding the Trotterized evolution S δt . For small δt, it is a sum of the true Hamiltonian H that we are simulating and a small error term V (due to the use of Trotterization): Therefore, under the symmetry protection, the effective Hamiltonian is replaced by its projection onto the Zeno subspaces: where V Zeno = µ P µ V P µ is the corresponding projection of V . In particular, if the error V does not respect the symmetry, the projection V Zeno could be much smaller than the error V in an unprotected simulation. The quantum Zeno framework therefore provides alternative intuition for the error suppression from the symmetry protection. We note, however, that choosing the symmetry transformations C k independently, instead of C k = C k 0 considered in this section, could lead to more reduction of the simulation error, and we demonstrate this advantage in Section IV.
We make these arguments rigorous by proving a bound analogous to that in Theorem 2 for symmetry-protected quantum simulation. Specifically, we consider G = H eff = H + V , where [H, U kick ] = 0. Note that under this assumption, the distinctiveness of the eigenvalues of U kick ensures that [P µ , H] = 0 for all µ in the spectral decomposition of U kick . We will also denote by V k = U †k kick V U k kick = G k − H. Theorem 4 (Symmetry protection by quantum Zeno effect). Let U kick be the unitary defined in Eq. (A2) and suppose that G = H + V such that [H, U kick ] = 0. Let G Zeno = µ P µ GP µ = H + µ P µ V P µ denote the projection of G onto the subspaces defined by a set of orthogonal projectors {P µ } in the spectral decomposition of U kick . We have where ξ is the inverse spectral gap defined in Eq. (A7).
Note that this bound is stronger than Eq. (A8) in that the dependence on the norm of the Hamiltonian is improved from G 2 to G V .  [1,s] Again, we prove Lemma 5 by induction on s. Suppose Lemma 5 holds for s = s 1 and s = s 2 such that |s 2 − s 1 | ≤ 1, we shall prove that it holds for s = s 1 + s 2 . Using the triangle inequality s1+s2 k=1 e −iG k δt − e −iG [ Proof. We have we obtain Eq. (B6).

Appendix C: A general bound on the Trotter error
In Section III, we prove a bound on the simulation error under the protection from a special class of symmetry transformations C k = C k 0 . In this section, we prove a similar, but more general, bound without making such an assumption.
Given a fixed total evolution time t, we first estimate the number of Trotter steps r required to simulate exp(−iHt) so that the total additive error of the simulation meets a threshold ε. Suppose the Hamiltonian H = L µ=1 H µ is a sum of L Hamiltonian terms H µ such that each e −iHµδt can be readily simulated on quantum computers. Again, we define the following quantities that are independent of t, r. The first-order Trotterization approximates exp(−iHδt) by where L µ=1 U µ ≡ U L . . . U 2 .U 1 is an ordered product. To get an accurate scaling of the gate count with the error tolerance, time, and the system size, we extend the approach in Ref. [20] to estimate the higher-order contributions to the total error. First, we estimate the higher-order contributions to the additive error in one Trotter step. Lemma 7. Assuming βδt ≤ 2α and α 2 δt ≤ γ + β, the Trotter error in approximating U δt = exp(−iHδt) by S δt in Eq. (C4) is given by where v 0 is defined in Eq. (15) and V(δt) is an operator bounded by with Λ = 5 6 (γ + β).
As a result of Lemma 7, we can bound the additive error in one Trotter step: Therefore, we arrive at a bound for the total error for the simulation where U kδt = exp(−iHkδt) and we have assume r E δt ≤ 1/2 to bound the sum over j. This bound again has the same feature as the bound in Theorem 1: the total error, to the lowest-order, scales with v 0 -an averaged version of v 0 under the symmetry transformations-instead of scaling with v 0 . Note, however, that the definition of v 0 here, with the addition of the transformations under U kδt , is slightly different from Theorem 1.
Appendix D: Proof of Lemma 1 In this section, we prove Lemma 1, which provides a bound on the error in one Trotter step.

(D4)
We note that G(τ 1 ) = O (τ 1 + τ 2 ) [Recall that O () is the standard Bachmann-Landau big-O notation.] Therefore, we can rewrite it (using either [20] or a direct differentiation) as with v 0 and F(τ 1 ) given above. Next, we rewrite the time-ordered exponential into a regular exponential using the Magnus expansion.
Lemma 8 (Magnus expansion [43][44][45]). Let A(τ ) be a continuous operator-valued function defined for 0 ≤ τ ≤ t such that with the sum being taken over all permutations σ of {1, . . . , j} and d b is the number of descents, i.e. pairs of consecutive numbers σ k , σ k+1 for k = 1, . . . , j − 1 such that σ k > σ k+1 , in the permutation σ. Furthermore, Ω j (t) are all anti-Hermitian if A(τ ) is anti-Hermitian. It is worth noting that the first two Ω j (t) are where the first-order Magnus term is To bound the higher-order terms in the Magnus expansion, we first note that (D28)