Digital Quantum Simulation of the Schwinger Model and Symmetry Protection with Trapped Ions

Tracking the dynamics of physical systems in real time is a prime application of digital quantum computers. Using a trapped-ion system with up to six qubits, we simulate the real-time dynamics of a lattice gauge theory in 1+1 dimensions, i.e., the lattice Schwinger model, and demonstrate non-perturbative effects such as pair creation for times much longer than previously accessible. We study the gate requirement of two formulations of the model using the Suzuki-Trotter product formula, as well as the trade-off between errors from the ordering of the Hamiltonian terms, the Trotter step size, and experimental imperfections. To mitigate experimental errors, a recent symmetry-protection protocol for suppressing coherent errors and a symmetry-inspired post-selection scheme are applied. This work demonstrates the integrated theoretical, algorithmic, and experimental approach that is essential for efficient simulation of lattice gauge theories and other complex physical systems.

The physical system we consider is a low-dimensional lattice gauge theory.Gauge field theories are the underlying formalism describing interactions among elementary particles in the Standard Model, and are prime candidates for modeling physics beyond the Standard Model.They further provide a powerful theoretical framework for describing low-energy excitations in condensed-matter systems.There has been tremendous progress in applying non-perturbative methods to solve lattice gauge theories in various systems and coupling regimes [46][47][48][49][50][51].However, evaluating the real-time dynamics of gauge field theories remains challenging in the strong-coupling regime, where a notorious sign problem halts Monte-Carlo-based classical simulations [52].Both fermionic and bosonic degrees of freedom in gauge theories, when defined on a discretized spacetime, can be mapped to spins and, in principle, efficiently simulated using quantum simulators.A simple Abelian lowdimensional gauge theory is quantum electrodynamics in 1+1 dimensions, or the Schwinger model [53].It exhibits non-trivial dynamics similar to those seen in quantum chromodynamics in 3+1 dimensions, i.e., the theory of the strong force in nature, including particleantiparticle pair creation, chiral symmetry breaking, confinement, and a non-trivial θ-vacuum [54,55].The Hamiltonian formulation of the lattice Schwinger model, i.e., the Schwinger model defined in a discretized space and with continuous time, has served as a testbed for numerous computational techniques in recent years, including quantum simulation and computation.In particular, there have been several theoretical proposals for analog quantum simulation [56][57][58][59][60][61][62] and gate-based quantum algorithms [22,24,[63][64][65][66][67][68][69][70] of the Schwinger model, along with experimental implementations on various quantum platforms such as trapped ions [71,72], Rydberg atoms [60,73], ultracold atoms [74,75], and superconducting qubits [64,76].
Inspired by the first digital simulation of the Schwinger model in Ref. [71], we revisit the simulation to answer the following questions: i) Are larger and longer simulations of lattice gauge theories viable on present hardware?ii) Given recent progress in quantifying the theoretical and empirical bounds on simulation algorithms, what are the resource requirements of simulating the lattice Schwinger model in the purely fermionic formulation with long-range interactions versus the fermion-boson formulation with only local interactions?iii) What are the theoretical considerations in decomposing the time evolution operator using product formulas, e.g., how does the Trotter error depend on the ordering of terms in each Trotter step of the evolution, what order of the product formula should one use, and how small should the Trotter steps be, considering the anticipated size of the experimental error?iv) Can recent cost-efficient symmetry protection protocols provide insights into the nature of experimental errors?
To address these questions, we realize experimentally the time dynamics of the lattice Schwinger model in its staggered formulation [77], and within its purely fermionic representation, for two-, four-, and six-site theories.Along with the experimental demonstration, different term orderings and product formulas are studied, and the gate complexity of the algorithm used here is compared with another algorithm in the local formulation of the same theory [22].The symmetry protection of Ref. [29] is implemented for the first time in experiment, in addition to a simple symmetry-inspired post-selection.

II. THEORY AND ALGORITHM CONSIDERATIONS
In the staggered formulation of the lattice Schwinger model introduced by Kogut and Susskind [77], the twocomponent matter field at one spatial site is split into two one-component fields, ψ, each occupying one site of the staggered lattice.This staggering corresponds to placing electrons at odd sites and positrons at even sites.The electric field, Ê, and the corresponding gauge-link variable, Û , are defined on the link connecting the two adjacent staggered sites.The Hamiltonian in natural units is: where N is the number of staggered sites, m is the mass of the fermions, a denotes the lattice spacing, and g is the coupling constant.The first term in the Hamiltonian involves gauge-boson-assisted fermionic hopping between nearest-neighbor sites.The second term is the energy stored in the electric field.The last term represents the rest mass of the fermions in the staggered formulation.
The theory exhibits a local gauge symmetry, which leads to a Gauss's law constraint on the allowed physical states |φ phys .This constraint enforces the net electric field at each site to be balanced by the electric charge present at the site, i.e., Ĝn ] for all n.The matter-field operators in Eq. ( 1) can be mapped to spin operators by the Jordan-Wigner transformation: With open boundary conditions (OBCs), the gauge links and the electric fields can be eliminated upon a gauge transformation and the application of Gauss's law [78].Without loss of generality, the electric field coming into the lattice is set to zero, and the resulting spin Hamiltonian reads Here, the original Hamiltonian is rescaled by 2 ag 2 , i.e., Ĥ = 2 ag 2 ĤKS , so that Ĥ in Eq. ( 2) is dimensionless with parameters x = 1 g 2 a 2 and µ = 2m ag 2 .Ĥx denotes the term proportional to the coupling x, ĤZZ consists of (nearly) all-to-all spin-spin interactions proportional to σZ n σZ m , and ĤZ is the sum of the terms proportional to σZ n .The constant terms are ignored in the following as they do not affect the evolution.By convention, σZ |0 = |0 , σZ |1 = − |1 , and the presence of a particle (antiparticle) is encoded as |1 (|0 ) at an odd (even) site.Therefore, for an N -site system, the bare-vacuum state, i.e., the ground state of the theory in the x = 0 limit, is The goal is to study the real-time dynamics of barevacuum fluctuations and particle-antiparticle pair creation.One can therefore perform a quench experiment where the simulation is initiated in the bare-vacuum state, |ψ(0) = |ψ 0 , which is then evolved via the Hamiltonian in Eq. ( 2) with x = 0.The survival probability of the bare-vacuum state is and the particle-number density is ν In the limit of N → ∞, the two quantities are related to each other via ν(t) = − 1 N log(P vac ) [53].The local charge density is defined as where the particle's (antiparticle's) charge is −1 (1).The local charge density is related to the local particlenumber density by Q n (t) = (−1) n ν n (t).As there is no interaction in the Hamiltonian that changes the total net charge of the system, n Q n is conserved.Therefore, the model has a global symmetry operator Ŝz = n σZ n , which is consistent with the symmetry of an XXZ Heisenberg spin Hamiltonian.
To perform the unitary evolution Û(t) = e −it Ĥ with Ĥ ≡ K k=1 ĥk on a quantum computer, with noncommuting Hamiltonian terms ĥk , the evolution can be broken into r smaller time steps of size δt = t/r.For each time step, the unitary Û(δt) is approximated by a product formula, which can then be implemented in terms of the available native gates.For the first-order product formula, while for the second-order product formula, where the second product is realized in reversed order.Higher-order formulas Ŝp (δt) can also be constructed recursively for even integers p [79].
While it is generally difficult to obtain an exact estimate for the Trotter error in the pth-order product formula, progress in recent years has resulted in tighter error bounds for them.For example, a nearly optimal bound is derived in Ref. [3], expressed in terms of the nested commutators between different ĥk terms.For the error of the second-order formula: Û(δt) − Ŝ2 (δt)

≤
(δt) 3  24 The norms of non-vanishing nested commutators can be further upper-bounded using the triangle inequality [A, B] ≤ 2 A B , resulting in closed-form, albeit looser, error bounds.Based on this approach, an upper bound for the pth-order formulas in simulating a twobody Hamiltonian can be obtained, where κ p = (4 × 5 p/2−1 ) p+1 is a constant, λ is the maximum strength of the sum of the interactions that involve any particular site, and γ = K−1 k=1 ĥk .Since γ only involves the first K − 1 terms of the Hamiltonian, labeling the index k such that ĥK is the term with the largest norm typically results in the smallest error bound.
Given a small δt, a higher-order formula results in a smaller simulation error.However, since the gate count of Ŝp (δt) scales as O 5 p/2 , increasing p increases the gate complexity of the simulation exponentially.The optimal choice for p is a balance between the gate depth, the evolution time t, and the error tolerance of the simulation.Since the experimental error can dominate the Trotter error at long evolution times, the first-order product formula, which minimizes the gate count per Trotter step, turns out to be a more suitable choice for the experiment in this work.

A. Term ordering
Generally, the Hamiltonian terms in Eq. ( 2) do not commute.Therefore, one needs to find the optimal ordering in the application of the product formula to minimize the Trotter error.There are two typical orderings of the terms in a nearest-neighbor Heisenberg spin model.The odd-even ordering is defined as having the odd-leading two-spin terms, e.g., στ 2k−1 στ 2k with τ = X, Y, Z and 1 ≤ k ∈ Z, applied first, then applying the even-leading two-spin terms, e.g., στ 2k στ 2k+1 , followed by the singlebody terms (if any).In the XYZ ordering, the terms involving σX are applied first, followed by the terms involving σY , and finally those involving σZ .It is known that the odd-even ordering introduces less Trotter error than the XYZ ordering [11].Here, we investigate whether the odd-even ordering is also a better choice for simulating the Schwinger model, which includes both nearestneighbor and (nearly) all-to-all spin-spin interactions.
Since Ĥx in the Hamiltonian in Eq. ( 2) includes only nearest-neighbor interactions, one way to define an oddeven-ordered product formula is to apply the odd-even ordering to terms in Ĥx only.Defining Ĥx n,m to be the term in Ĥx acting on sites n and m, and, similarly, ĤZZ n,m to be the term in ĤZZ acting on sites n and m, this ordering can be written as: Alternatively, the odd-even ordering can be applied to both Ĥx and the nearest-neighbor terms in ĤZZ , followed by the application of the non-nearest neighbor terms in ĤZZ as well as ĤZ ,  11) (blue dots) and the XYZ ordering defined in Eq. ( 14) (red diamonds).The blue line denotes the exact evolution.
The two odd-even-ordered evolution operators are related by a rotation around the Z-axis: As a result, if the state is initialized and measured in the Z-basis, one arrives at the same measurement outcome using either scheme.Therefore, in the following we will only consider Ŝoe1 1 because it requires implementing fewer single-qubit gates.
The symmetry operator Ŝz commutes with σX ., which is not broken up by the odd-even ordering scheme.Therefore, this ordering does not result in any leakage to the symmetry-forbidden subspace.This is unlike the XYZ ordering that implements where Ĥ(XX) k,k+1 and Ĥ(Y Y ) k,k+1 are the terms in the Hamiltonian in Eq. ( 2) proportional to σX k σX k+1 and σY k σY k+1 , respectively.The size of the leakage to the symmetryforbidden subspace for both schemes are verified numerically in Fig. 1.

B. Gate complexity
There are different approaches to digitally simulating the Schwinger model.Ref. [22] truncates the gauge-boson degrees of freedom in the electric-field basis, |E| ≤ Λ.The electric-field Hilbert space at each link is then encoded using log 2 (2Λ + 1) qubits.With OBCs and a zero incoming electric-field flux, the exact theory is recovered for Λ = N/2.To simulate the model with N sites for time t, this approach requires O(N + N log 2 N ) qubits and, up to logarithmic corrections, O N 5/2 t 3/2 twoqubit gates, using the second-order Suzuki-Trotter formula.In contrast, our approach integrates out the gauge bosons, leaving only the fermionic degree of freedom associated with the matter fields [71,78].Therefore, it requires O(N log 2 N ) fewer qubits to simulate the same model.However, in this approach, one trades the local gauge-matter interactions for long-range matter-matter interactions that correspond to a Coulomb force.Here, we investigate whether this purely fermionic model increases the gate complexity of the simulation compared with the local formulation involving both fermions and bosons.
Given a fixed error tolerance, one can use Eq. ( 10) to estimate the required number of Trotter steps in the purely fermionic formulation for simulation time t, system size N , and at fixed x and µ: Here, λ = O N 2 , which is determined by the interactions in ĤZZ , and γ = O (N ), which results from summing over only the interactions in Ĥx . 1 Since each Trotter step requires O N 2 two-qubit gates, the gate complexity of the pth-order product formula is O N 4+1/p t 1+1/p .For the second-order product formula (p = 2), the gate complexity reduces to N 9/2 t 3/2 , which is a factor of N 2 larger than the scaling of Ref. [22].The bound on the Trotter error derived above is an upper bound and the required gate count may be much smaller in practice.In Fig. 2, we plot the empirical value of the number of two-qubit gates required for simulating the time-evolution operator in the purely fermionic formulation using the second-order formula for time t = N .The empirical gate count is obtained through a binary search for the minimum number of Trotter steps such that the difference between the Trotterized and the exact evolution is at most = 0.01.This value is compared to two bounds: the commutator bound obtained from Eq. (10), following the discussion above, and the exact The number of two-qubit gates required to simulate the time evolution under the Hamiltonian in Eq. ( 2) on an Nsite lattice for time t = N given an error tolerance = 0.01, and using the second-order product formula in Eq. ( 8).The commutator bound on the gate complexity (green dots) is estimated from the error bound in Eq. (10).We also obtain a tighter estimate (orange dots) by exactly computing the nested commutators in Eq. ( 9).Finally, the empirical gate count (blue dots) is obtained through a binary search for the minimum number of time steps t/δt such that the total error is at most .The straight lines are linear fits that result in the polynomial scaling given in the figure .commutator bound derived from computing the norms of the nested commutators in Eq. ( 9) exactly.The gate complexity O N 6.2 from the commutator bound agrees with our earlier estimate of O N 6 for this formulation of the model.However, Fig. 2 shows that this bound, obtained by applying the triangle inequality on the norms of the nested commutators, is rather loose.Instead, by computing Eq. ( 9) exactly, hence invoking cancellations between the commutators, the gate complexity reduces to O N 4.6 , very close to the empirical scaling of O N 4.2 .Notably, the empirical scaling is nearly the same as the scaling of Ref. [22] for the fermion-boson formulation.However, we note that one may also be able to obtain tighter error bounds on the algorithm of Ref. [22] empirically.

C. Symmetry protection
Digitizing a quantum evolution may introduce errors that populate states forbidden by the symmetry of the target system.Ref. [29] proposes a method to reduce this leakage when using product formulas and applies it to the fermion-boson formulation of the Schwinger model.In this section, we discuss whether the proposed method is effective at protecting the symmetry of the purely fermionic formulation, and in reducing the Trotter error.
As mentioned earlier, the Schwinger-model Hamiltonian in Eq. ( 2) is invariant under a global rotation around the Z axis of the qubits, i.e., [ Ĥ, Ŝz ] = 0. Therefore, the expectation value of Ŝz is conserved if the evolution is exact.However, due to Trotterization errors, Ŝz may deviate from its initial value during the simulation.To mitigate this error, one can insert rotations generated by the symmetry, i.e., Ĉ(α) = e −iα Ŝz , in between the Trotter steps [29]: 16) By choosing suitable angles α k for different Trotter steps, the errors from each step that do not commute with Ŝz are rotated by different amounts and can interfere destructively.This results in smaller leakage to the sector with the global charges that differ from the initial state.Therefore, the method mitigates errors in symmetryviolating Trotterization schemes, such as the XYZ ordering.Errors that commute with Ŝz remain intact under these rotations Ĉ(α k ).As the result, this method has no effect on symmetry-preserving Trotterization schemes, such as the odd-even ordering.Additionally, since the symmetry operator Ŝz is diagonal in the measurement basis for the observables considered here, one can simply mitigate symmetry-violating errors by post-selecting on measurement outcomes that preserve Ŝz .We compare the two mitigation methods in the following.
For simplicity, the angles α k in Eq. ( 16) can be chosen such that α k = kα 1 , where for each t, α 1 is determined by numerically minimizing the leakage to the symmetry-forbidden subspace after t/δt steps, as detailed in Appendix A. As a result, the optimal value of α 1 depends on the number of Trotter steps and the simulation time.The top panel in Fig. 3 plots the leakage to the symmetry-forbidden sector for the odd-even ordering, the XYZ ordering before and after post-selection, as well as the XYZ ordering with symmetry protection but without post-selection.As expected, the odd-even ordering results in no population in the symmetry-forbidden states, and the leakage seen in the XYZ ordering trivially goes to zero after post-selection.Importantly, the symmetry-protected XYZ ordering leads to a strong suppression of the leakage.
However, this improvement in preserving the symmetry does not guarantee a smaller Trotter error.Given two operators E sym and E sym , corresponding to the symmetrypreserving and symmetry-violating errors in the simulation, respectively, it is generally not true that E sym < E sym + E sym .In particular, depending on the observable of interest, the effect of the two types of errors may interfere destructively.Therefore, eliminating E sym may actually increase the error in the expectation value of the observable.For example, as shown in the lower panel of Fig. 3, the symmetry protection sometimes increases the error in the population of the bare-vacuum state instead of decreasing it.Overall, symmetry protection does not improve the accuracy of the XYZ-ordering scheme over the odd-even ordering scheme.Notably, post-selection is more successful than the symmetry-protection scheme at The leakage to the symmetry-forbidden sector, Psym, defined as the population in the states with a nonvanishing total charge given a bare-vacuum initial state (upper panel), as well as the error in the bare-vacuum population (lower panel) are shown for the odd-even ordering (blue dots), the XYZ ordering before (red diamonds) and after (green squares) post-selection, as well as the XYZ ordering with symmetry protection but without post-selection (yellow stars).The parameters used for the plot are µ = 0.1, x = 0.6, N = 4 and δt = 1.The optimal angles for the symmetryprotected simulation are provided in Appendix A.
mitigating the Trotter error in P vac .
According to Ref. [29], the symmetry protection scheme mitigates time-correlated experimental errors as well.Since the errors in our experiment are expected to be dominated by uncorrelated noise, it is interesting to investigate how well the scheme performs in the experimental implementation.To isolate the effects of symmetry protection on experimental rather than Trotter errors, we implement the odd-even ordering, which preserves the symmetry.The results are presented at the end of the next section.

III. EXPERIMENT AND RESULTS
The experiment consists of a chain of up to nine 171 Yb + ions held in a linear Paul trap, with up to six ions used as qubits.The two qubit states are realized in the hyperfinesplit ground level, |0 = | 2 S 1/2 F = 0, m F = 0 and |1 = | 2 S 1/2 F = 1, m F = 0 .The qubits are initialized to |0 by optical pumping, and read out on a multi-channel photo-detector using state-dependent fluorescence [80].We use two counter-propagating Raman beams to coherently manipulate the states of the qubits [81].One Raman beam is split into individual beams to address each qubit separately.The native gates in this setup are single-qubit Z rotations, Z i (θ) = e −iθ/2 σZ i , singlequbit rotations around the equatorial plane of the Bloch sphere R i (θ, φ) = e −iθ/2 (σ X i cos φ+σ Y i sin φ) , and two-qubit gates X i X j (χ) = e −iχσ X i σX j , where φ is the angle of the axis of rotation, and θ and χ are the rotation angles.The Z rotations are implemented as classical phase advances, and are therefore practically noise free.The single-qubit gates R i (θ, φ) are realized by driving a qubit on resonance, with the laser phase set to φ and the duration proportional to θ.The two-qubit gates use the shared motional modes of the ion chain to mediate interactions between pairs of ions [82][83][84].Carefully designed amplitude modulation of the laser pulse leaves the qubit states decoupled from the motion at the end of the gate [85].The individual addressing beams allow any pair of ions to be entangled, making the system a fully-connected programmable quantum computer.More details on the setup and gate fidelity are presented in Ref. [81,86,87] To simulate the evolution under the Hamiltonian in Eq. ( 2), the product formula needs to be decomposed into the native gates in the experiment.Specifically, to implement the gate, each qubit is rotated around the Z (Y ) axis before and after applying an X i X j (χ) gate.For the odd-even ordering, i.e., S oe1 1 given in Eq. ( 11), the circuit can be broken down into three parts: e −iδt Ĥx , e −iδt ĤZZ , and e −iδt ĤZ .An example of the circuit for the Trotterized evolution with N = 6 is shown in Fig. 4. The number of gates needed in each Trotter step for the different values of N in the experiment are summarized in Table I.While the number of single-qubit rotations can be further reduced in the circuits, such an optimization is not considered here since they have much lower error rates than the two-qubit gates.
The experiment starts by preparing the qubits in the vacuum state of the Hamiltonian in the limit of x = 0, |vac = |ψ 0 , using R i (π, 0) on alternate qubits.The model parameters are set to µ = 0.1 and x = 0.6 to ensure non-trivial dynamics starting from the bare vacuum.We measure all the qubits in the Z basis to get the state populations.The state preparation and measurement errors are corrected by applying the inverse of an independently measured state-transfer matrix.The observables of interest are then calculated from the state populations using Eqs.(3)- (5).
Algorithmic and experimental errors can, in principle, produce measurement results that break the global  2) for N = 6 lattice sites, with odd-even ordering of the Trotter decomposition introduced in Eq. (11).The interaction term e −iδt Ĥx is implemented with nearestneighbor XiXi+1 and YiYi+1 gates.The e −iδt ĤZZ term is implemented with XiXj and Ri rotations.The e −iδt ĤZ term involves only Zi rotations, with the angles µ1 = −(µ + 3)δt, µ2 = (µ − 2)δt, µ3 = −(µ + 2)δt, µ4 = (µ − 1)δt, µ5 = −(µ + 1)δt, and µ6 = µδt.Qubits are initialized in the bare-vacuum state |010101 , then evolved by repeating the circuit in the parentheses t/δt times, and measured individually in the Z basis in the end.charge conservation of the system.This means that from the bare-vacuum state, the probability amplitude for evolving to states with non-vanishing total charge may not be negligible.Therefore, to improve the results, one can post-select the measurements so that only outcomes in the relevant symmetry sector are kept.Since the odd-even ordering implemented in the experiment does not violate the global charge conservation, these errors result entirely from the experimental imperfections.
Fig. 5(a) plots the bare-vacuum population P vac , Eq. (3), and the particle-number density ν, Eq. ( 4), as a function of time for N = 2 and δt = 0.5.For δt = 0.5, the theoretical Trotterized evolution (blue dots) has no significant deviation from the exact evolution (blue line).The experimental results after post-selection agree well with the theory even after 39 Trotter steps corresponding to t = 19.5.Post-selection is shown to significantly mitigate the experimental errors for P vac , especially at long evolution times.However, it does not appear as effective for the particle-number density or the survival amplitude of other initial states in general, as plotted in Appendix B. Fig. 5(b) plots the local charge density Q n , Eq. ( 5), derived from post-selected measurement results and their deviation from the theoretical expectation for the same set of parameters.The local-charge profile is consistent with the global dynamics shown in Fig. 5(a) as the pair creation coincides with the destruction of the bare vacuum.I, and we only perform ten Trotter steps for this case.Since the experimental error dominates the Trotter error after a few Trotter steps for δt = 0.5, the Trotter step size is increased to δt = 1 in Fig. 7, doubling the simulation time.Even though the Trotter error is larger in Fig. 7 than in Fig. 6, the Trotterized evolution still qualitatively follows the exact solution.In both cases, the experimental data after postselection agree reasonably well with the numerical simulation.
Next, we increase the number of staggered sites to N = 6, with the results displayed in Fig. 8. Since each Trotter step now requires 20 two-qubit gates and 38 single-qubit gates, only four Trotter steps could be run before decoherence dampens the evolution.Nevertheless, the qualitative behavior, including the first revival of the bare-vacuum amplitude, can still be observed from the results.
Finally, we study the effect of active symmetry protection by inserting rotations generated by the symmetry operator Ŝz , as discussed in Sec.II.In the numerical study of Sec.II, a set of optimal angles was found to minimize the leakage to the symmetry-forbidden subspace caused by algorithmic error.Since such an optimal set cannot be found a priori, a straightforward strategy is to use a random angle at each Trotter step to average out the leakage after many Trotter steps [29].Since the effectiveness of the scheme depends on the nature of the experimental error, it is interesting to see if symmetry protection can improve our experimental implementation.
Figure 9 plots the result of an experiment using the odd-even term ordering.As before, the initial state is the bare vacuum.The unitaries e −iα k Ŝz , with random angles α k given in Appendix A, are inserted between Trotter steps k and k + 1.While the population in states forbidden by the symmetry, denoted as P sym in the upper panel, decreases with symmetry protection, this reduction is not significant.Furthermore, while the deviation of the barevacuum population from the Trotterized theory generally decreases, post-selecting symmetry-preserving measurements appears more effective in mitigating the error in this quantity than the symmetry protection as shown in the lower panel of the figure.This indicates that the experiment is dominated by noise that is not correlated in time.Note that by construction, the symmetryprotection scheme only mitigates time-correlated errors.

IV. DISCUSSION
We have digitally simulated the time evolution of the lattice Schwinger model with up to six qubits using the purely fermionic formulation.For a four-qubit simulation, we observe four oscillations of the particle density, and the simulated time is almost four times longer than previously demonstrated using a Trotterized scheme [64,71].Given the long circuit depths required for dynamical gauge-theory simulations, gate fidelity, rather than qubit number, is the main limitation of our implementation.Efforts to overcome such a technical limitation are well underway [88].To mitigate the timecorrelated errors, we have applied a symmetry-protection scheme [29] but found negligible effects in suppressing the errors.This symmetry-protection investigation indicates that the dominant noise in experiment is incoherent and uncorrelated.Incoherent errors can be mitigated by post-selection of the experimental measurements using symmetry considerations.Better-motivated and furthertailored schemes for mitigating incoherent errors are desired in future simulations.Furthermore, it is found that the symmetry-violating and symmetry-preserving errors can destructively contribute to given quantities, and removing only one of these errors can decrease the overall accuracy.It is therefore important to also develop and apply symmetry-preserving error-suppression schemes in future experiments.
An avenue for improving the quality of the simulation is reducing the gate depth, e.g., by performing gates in parallel instead of sequentially.In our model, e −iδt Ĥx , consisting of only nearest-neighbor interactions, can be applied in a fixed circuit depth of 4 instead of 2N by performing all the X 2i X 2i+1 terms, then all the X 2i+1 X 2i+2 terms, in parallel.The all-to-all interactions in e −iδt ĤZZ can be reduced to depth of N instead of N 2 if gates X i X i+n , for all i and fixed n, are performed in paral- lel.With trapped ions, parallel operations can be done either in multi-zone architectures [89,90], or in linear chains with advanced control schemes [91].
Alternatively, the gate depth can be reduced by using M -body Mølmer-Sørensen gates MS(χ, M ) ≡ e −iχ M i=1 M j=i+1 σX i σX j [82][83][84].This approach was implemented in Ref. [71] to reduce the number of MS operations in the simulation of the Schwinger model from O(N 2 ) to O(N ).In general, a non-trivial optimization of both frequency and amplitude modulation of the beams may be required to implement an M-body gate with the desired rotation angles, as demonstrated in Ref. [59].Furthermore, one should note that since an M -body MS gate has roughly comparable fidelity as that of M 2-body MS (X i X j ) gates [92], the overall fidelity of the simula-Figure 9. Leakage to the symmetry-forbidden subspace (upper panel) and deviation of the experimental results from theory for the bare-vacuum population (lower panel) when different schemes to mitigate errors are applied: no mitigation (orange triangles), post-selection (green squares), and symmetry protection (blue stars).Note that by definition, the leakage to the symmetry-forbidden subspace is zero after post-selection.In both plots, N = 4 and δt = 1.0, the system is initiated in the bare-vacuum state, and is evolved via the odd-even ordering scheme.tion would likely be similar for both schemes.
When the fermion-boson formulation of the lattice Schwinger model is concerned, a trapped-ion-specific approach to reduce the gate count is to encode the gauge degrees of freedom into the motion of the ions as explained in Ref. [68].Besides the standard set of gates, the proposed hybrid digital-analog scheme involves both spin-phonon and phonon-phonon gates.This approach leads to a reduction in both the number of qubits and the number of entangling gates compared with a fully-digital algorithm involving the gauge bosons [22].Future experimental implementations will determine the realistic fidelity of the operations involving dynamical phonons.
To make implementations of more complex gauge theories possible, including non-Abelian and higherdimensional models, unifying physics insights, algorithm optimization, hardware implementation, and postprocessing is required, as demonstrated in this work.In this context, it would be interesting to investigate whether more resource-efficient encodings of such theories exist, if optimal Trotter decompositions and term ordering schemes can be found, to what degree these preserve local gauge symmetries, whether information regarding the initial state and the symmetries can be incorporated to further tighten the algorithmic error bounds [70,[93][94][95][96][97][98], how to balance these errors with experimental errors, and whether symmetry-protection schemes are advantageous in suppressing algorithmic and experimental errors.While progress along these lines is already being made [23,24,27,30,32,36,67,[99][100][101][102][103][104][105][106][107], further technological advances in quantum hardware are essential to enable advanced gauge-theory simulations in the upcoming years.

Figure 1 .
Figure1.Numerical simulation of the projection on the bare-vacuum state Pvac (upper panel) and the population of symmetry-forbidden states Psym (lower panel) when the system is initialized in the bare-vacuum state for N = 6, µ = 0.1 and x = 0.6.Different term orderings for the Trotterized evolution are considered: the odd-even ordering defined in Eq. (11) (blue dots) and the XYZ ordering defined in Eq. (14) (red diamonds).The blue line denotes the exact evolution.

Figure 3 .
Figure3.The leakage to the symmetry-forbidden sector, Psym, defined as the population in the states with a nonvanishing total charge given a bare-vacuum initial state (upper panel), as well as the error in the bare-vacuum population (lower panel) are shown for the odd-even ordering (blue dots), the XYZ ordering before (red diamonds) and after (green squares) post-selection, as well as the XYZ ordering with symmetry protection but without post-selection (yellow stars).The parameters used for the plot are µ = 0.1, x = 0.6, N = 4 and δt = 1.The optimal angles for the symmetryprotected simulation are provided in Appendix A.

Fig. 6
plots the same observables as a function of time for N = 4 and δt = 0.5.Compared with the N = 2 case, each Trotter step for N = 4 requires seven extra two-qubit gates, see Table

N = 2 ,Figure 5 .
Figure 5. Experimental results for N = 2 and δt = 0.5.(a) The upper plot shows fluctuation in the bare-vacuum population, Pvac, while the lower plot shows particle-number density, ν, as a function of time, indicating the creation and annihilation of the particle-antiparticle pairs.The dashed lines are a guide to the eye.(b) The upper plot shows the local charge density Qn as measured in the experiment after postselection, while the lower plot shows its deviation from theory as a function of time.

N = 4 ,Figure 6 .
Figure 6.Experimental results for = 4 and δt = 0.5.(a) The upper plot shows fluctuation in the bare-vacuum population, Pvac(t), while the lower plot shows particle-number density, ν(t).(b) The upper plot shows the local charge density Qn(t) as measured in the experiment after post-selection, while the lower plot shows its deviation from theory.

N = 4 ,Figure 7 .
Figure 7. Experimental results for N = 4 and δt = 1.(a) The upper plot shows fluctuation in the bare-vacuum population, Pvac(t), while the lower plot shows particle-number density, ν(t).(b) The upper plot shows the local charge density Qn(t) as measured in the experiment after post-selection, while the lower plot shows its deviation from theory.

N = 6 ,Figure 8 .
Figure 8. Experimental results for N = 6 and δt = 1.(a) The upper plot shows fluctuation in the bare-vacuum population, Pvac(t), while the lower plot shows particle-number density, ν(t).(b) The left plot shows the local charge density Qn(t)as measured in the experiment after post-selection, while the right plot shows its deviation from theory.At t = 4, we reach the gate-depth limit of the hardware.

Figure 10 .
Figure 10.The population of the symmetry-forbidden subspace as a function of the angles α1 for t = 2, 4, 8, for µ = 0.1, x = 0.6, N = 4 and δt = 1.We choose α1 to be the smallest angle that minimizes Psym at each t.

Table I .
Gate counts for simulating each Trotter step of the time evolution in the Schwinger model with the odd-even term ordering in Eq. (11), along with the largest number of Trotter steps t/δt implemented in the experiment for N = 2, 4, 6 staggered sites.

Table II .
Values of α1 used for the symmetry-protected XYZ ordering in Fig.3at different times t, together with the corresponding leakage to the symmetry-forbidden subspace.

Table III .
Randomly chosen values of α k used in the experiment to study the effect of symmetry protection on experimental error.The result is presented in Fig.9.