Gauge invariant quantum circuits for $U(1)$ and Yang-Mills lattice gauge theories

Quantum computation represents an emerging framework to solve lattice gauge theories (LGT) with arbitrary gauge groups, a general and long-standing problem in computational physics. While quantum computers may encode LGT using only polynomially increasing resources, a major openissue concerns the violation of gauge-invariance during the dynamics and the search for groundstates. Here, we propose a new class of parametrized quantum circuits that can represent states belonging only to the physical sector of the total Hilbert space. This class of circuits is compact yet flexible enough to be used as a variational ansatz to study ground state properties, as well as representing states originating from a real-time dynamics. Concerning the first application, the structure of the wavefunction ansatz guarantees the preservation of physical constraints such as the Gauss law along the entire optimization process, enabling reliable variational calculations. As for the second application, this class of quantum circuits can be used in combination with timedependent variational quantum algorithms, thus drastically reducing the resource requirements to access dynamical properties.


I. INTRODUCTION
Gauge theories lie at the heart of the Standard Model of particle physics and represent the most successful description of elementary particles and their fundamental interactions [1,2]. For instance, these theories represent a generalization of quantum electrodynamics (QED) which elucidates the behavior of charged particles and photons, and quantum chromodynamics (QCD) which describes the strong interactions between quarks, the elementary constituents of protons and neutrons. While perturbative approaches exist in the study of gauge theories, for example describing weak or electromagnetic interactions [3], other theories can be approached with these methods only in certain limiting regimes. For instance, QCD can be studied perturbatively in the limit of high energy while in the low-energy regime however, the strong interactions grow so large that a perturbative analysis is not possible anymore.
A crucial step towards a non-perturbative analysis of gauge theories such as QCD was taken by Wilson who formulated gauge theories on a finitely discretized spacetime lattice while preserving the exact local symmetry of a gauge theory, thereby introducing the concept of lattice gauge theories (LGT) [4]. Despite the success of path-integral Monte Carlo sampling methods and tensor network approaches [5][6][7][8][9][10][11][12][13][14][15][16][17] for estimating equilibrium properties of LGT [18], out-of-equilibrium and real-time dynamics have remained out of reach, as they are limited by a severe numerical sign problem that exponentially increases the time complexity with increasing system size [19], as well as by the exponential growth of the entanglement with the simulation time in arbitrary dimensions [20].
Quantum computing offers a promising alternative, as it does not suffer from the above limitations. By expressing LGTs in the equivalent Hamiltonian formalism [21], quantum simulators may be built to simulate the LGT dynamics [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40] with resources that only grow polynomially with the system size. However, one main problem that arises in the standard Hamiltonian formulation of a LGT [21] is the emergence of (exponentially many) unphysical states [28,38,41,42]. The origin of such an unphysical subspace is related to the additional freedom associated to a particular choice of the gauge fixing implicit to the canonical quantization of the Hamiltonian. In the case of the temporal (partial) gauge, the physical subspace is defined by additional so-called Gauss law constraints. For instance, in the case of QED the Gauss law takes the familiar form ∇ · E = −e ψ † ψ which relates the electric charge distribution −e ψ † ψ to the electric field E. As a consequence, such constraints resulting from the gauge symmetry of LGT must be incorporated either into the encoding of the relevant degrees of freedom or into the operators which govern the evolution of the system.
In this manuscript, we investigate the task of state preparation while retaining the gauge symmetry constraints imposed by LGT as mentioned above. State preparation is the first step in all quantum algorithms targeting ground state calculations, ranging from the variational quantum eigensolver (VQE) [43] to quantum phase estimation (QPE) [44,45], thus representing a key stage in virtually every kind of quantum simulation. State preparation is also crucial for the initialization of non-equilibrium processes (assuming that the initial state of interest is not trivially represented in the computational basis) originating from a quench [46].
To this end, we introduce a class of parametrized quantum circuits that realize states that are gauge invariant by construction and applicable to generic Yang-Mills LGT theories in any spacetime dimension. Notably, the proposed approach is independent of the choice of a specific qubit encoding. We observe a clear advantage over an unstructured variational form ansatz, both in the estimation of ground state wavefunctions as well as in the arXiv:2105.05870v2 [quant-ph] 20 Jan 2022 implementation of real time dynamics.

II. YANG-MILLS MODEL HAMILTONIAN
We start by introducing a lattice model that reproduces the dynamics of a Yang-Mills gauge theory with dynamical fermionic matter in the continuum limit [21,38,[47][48][49]. In Yang-Mills gauge theories, multiple matter particle species ψ a (representing different types of quarks, for example) can exist such that the theory is invariant under certain local transformations V (x) that might mix the different particle species ψ a (x) → V ab (x)ψ b (x) at every space-time point x. These local transformations are characterized by an Abelian or non-Abelian compact Lie group G, the so-called the gauge group of the theory [2]. More specifically, the transformations V (x) correspond to a finite unitary representation of the gauge group G. In this setting, QED can be regarded as a Yang-Mills theory with the Abelian gauge group U (1) with one single particle species representing the electron, while QCD corresponds to a Yang-Mills theory with the non-Abelian gauge group SU (3) describing three particle species which represent the colors of quarks.
In this work, the Yang-Mills theory is formulated on a d-dimensional discrete spatial lattice with lattice spacing a. We denote lattice sites in d dimensions by their coordinate vector x ∈ R d and edges (also called links) between sites by the tuple (x, k) where x is a lattice site and k ∈ {1, ...d} is a direction. The (single flavor) fermionic particles reside on the lattice sites x, while the gauge fields live on the links. The Yang-Mills Hamilto-nianĤ YM =Ĥ hopp +Ĥ mass +Ĥ plaq +Ĥ elec +Ĥ wilson in the temporal gauge readŝ where γ k are given by the Dirac gamma matrices and where m, g denote the fermionic mass parameter and the gauge field coupling parameter, respectively. This Hamiltonian mainly consists of four contributing terms describing the fermionic matter (Ĥ mass ) and gauge field energies (Ĥ elec andĤ plaq ) and their interactions (Ĥ hopp ). An additional termĤ wilson is introduced via the parameter r which acts as a regulator to avoid the fermion doubling problem [50,51]. The quantitiesψ x and ψ x =ψ † x γ 0 represent the fermionic field operators that carry implicit spinor-and color component indices i, j ∈ 1, . . . , 2 d/2 and b, c ∈ {1, . . . , dim(R)} respectively, where dim(R) denotes the dimension of the representation R of the corresponding gauge group G. For example, in the case of G = SU (N ) in the fundamental representation, we have dim(R) = N . Note that whenever the spinor-and color indices are suppressed, we implicitly sum over all these indices. The fermionic operators satisfy the standard anti-commutation relations The gauge field operatorsÛ (x,k) on a link (x, k) have the same structure as the generators T α with α ∈ {1, . . . , dim(G)} of the gauge group: They are dim(R) × dim(R) dimensional matrices whose entries (Û (x,k) ) mn are operators that act on the link Hilbert space. The plaquette operators are defined byÛ = , and the trace of these operators appearing inĤ YM is only taken in the Lie algebra space of the group generators T α and does not extend to the quantum Hilbert space. The conjugate link flux variablesL α (x,k) andR α (x,k) are associated with the left end and the right end of a link (x, k), respectively, and are defined through the commutation relations [52] (3) Note that in the case of QED, these two variables reduce to the electric flux field operatorsÊ (x,k) used in [38]. In addition, a constant background gauge field is introduced via the constant parameters θ α k [53]. Note that all operators assigned to different lattice sites or links commute with each other. The Gauss law operators are defined bŷ and they specify the gauge-invariant states |φ phys that span the physical sector H phys of the total Hilbert space via |φ phys ∈ H phys ⇔Ĝ α x |φ phys = 0

III. ENFORCING THE GAUSS LAW CONSTRAINTS
In the standard Hamiltonian formulation of Yang-Mills LGT, we have seen that the model Hamiltonian (1) operates on a large Hilbert space H YM containing all possible matter and gauge field states, especially including the unphysical states, which are not gauge-invariant. When studying the spectrum of the Yang-Mills Hamiltonian using for instance a variational approach, one has to remain within the physical subspace H phys with gauge-invariant states |φ phys in order to guarantee the correct physical energy spectrum. On the other hand, in the case of real-time evolution initiated from a gauge invariant state |φ phys (0) , the time-evolved state |φ phys (t) at time t will remain gauge invariant since the Gauss law constraints are constants of motion, [Ĝ α x ,Ĥ YM ] = 0. For this reason, there would in principle be no need to additionally enforce the Gauss law. Nevertheless, errors that kick the state out of the physical Hilbert space H phys can occur due to the Trotter approximation [55,56] of the timeevolution operator or from quantum hardware noise.
Various implementations of the Gauss law constraints have been proposed in the literature, ranging from absorbing these constraints directly into the Hamiltonian [36,[57][58][59][60][61] (thereby effectively removing the explicit gauge or matter degrees of freedom), to their incorporation into the state ansatz (gauge and matter particle) [62,63]. However, these approaches are often limited to specific models and particular dimensions and are not of general applicability for digital quantum simulations.
A more pragmatic approach consists in enforcing the Gauss law constraints by means of an energy penalty added to the system Hamiltonian and proportional to the regularization parameter λ [12,23,38,64]. This term effectively lifts the energy of unphysical states while the physical gauge-invariant states, lying in the kernel of the Gauss law operators, remain unaffected. As a result, for large enough values of λ, the low-lying energy spectrum is solely associated with physical states. It is worth mentioning that modifications of this approach exist in which the penalty term is replaced by simpler terms that become linear in the Gauss law operators [65][66][67].
This approach is general enough to be applied regardless of the lattice dimensionality, the gauge group and the qubit encoding. However, the penalty term implements the Gauss law constraints only approximately, requiring a tuning of the parameter λ, and, in addition, they come at the expense of increasing the complexity of the Hamiltonian. In the context of variational quantum algorithms for groundstate optimization, we further argue in Appendix A and B 2 thatĤ gauge leads to a more corrugated energy landscape (see for instance [68]) and that a large regularization parameter λ induces significant sampling errors in estimating the energy expectation value, thereby aggravating the optimization process in the variational algorithm.
The above considerations thus motivate the preparation of states for which the Gauss law constraints need not be additionally enforced but are automatically fulfilled by construction.

IV. GAUGE-INVARIANT STATE PREPARATION
Heuristic state preparation ansätze [37], which are common for example in quantum chemistry applications [69,70], fail to sample efficiently from H phys (see Appendix A). To circumvent this problem, we construct parametrized trial states tailored to the gauge symmetry of the theory. To do so, we look for unitary operatorŝ O that map the physical Hilbert space H phys onto itself, i.e.Ô |φ phys ∈ H phys for all |φ phys ∈ H phys . We call unitary operators with this property gauge invariant. A parametrized familyÛ(θ), with real parameter vector θ, of gauge invariant operators then allows us to construct trial states from any initial state |φ init ∈ H phys via The definition of H phys implies that a family of unitary operatorsÛ(θ) is gauge invariant if and only if it commutes with the Gauss law operatorsĜ α x in Eq. (4), Û (θ),Ĝ α x = 0 for all θ, α, x . With this choice, it follows that a parametrized trial state |φ(θ) is gauge invariant if and only if the initial state |φ init is chosen to be gauge invariant sinceĜ α x |φ(θ) =Û(θ)Ĝ α x |φ init = 0 . As a result, such a gauge-invariant variational form will sample from the physical Hilbert space only. Furthermore, note that due to the linearity and the product rule of the commutator for arbitrary operators A, B and C, it follows that any linear combination or product of gauge invariant terms will remain gauge invariant. Thus, if a Hermitian operator commutes with the Gauss law operators F ,Ĝ α x = 0, then the exponential ofF is a unitary operator that commutes with the Gauss law operators, exp(iF ),Ĝ α x = 0, too. Following these prescriptions, a parametrized family of gauge-invariant unitary operatorsÛ(θ) can be generated by combining pieces of the Yang-Mills Hamiltonian H YM in Eq. (1) which commute with the Gauss law operators. These pieces can be promoted to the following parametrized terms: where Λ(θ) is a real scalar, and the coupling matrix A(θ) is an arbitrary n spinor × n spinor matrix that mixes the spinor components of the fermionic matter fields. We allow the mixing of the spinor components since these spinors are all affected in the same manner under a gauge transformation and a mixing does not change the gauge invariance properties of these terms. On the other hand, mixing the color components in an arbitrary manner would destroy the gauge symmetry of these terms. Note that we implicitly sum over all color and spinor indices b, c and i, j. For example, Eq. (9) would more ex- In the string-like term, the link variablesÛ n can be chosen such that their product forms an arbitrary path which connects the lattice sites x and x + l along lattice edges. The set C in the loop-like term denotes an ordered sequence of links that form a loop. We then construct a gauge invariant unitary variational formÛ(θ) by defininĝ whereF k (θ) denotes any linear combination or product of the gauge invariant terms (8) - (14). In this respect, the present variational form is conceptually different from the so-called Hamiltonian variational ansatz introduced in [71], inspired by a trotterization of an annealing process.
Conceptually, the set of operators (8) - (14) is to be considered as a toolbox for preparing generic gaugeinvariant states in a particular state space H ⊂ H phys . For instance, considering a one-dimensional QED system in a low-energy regime, one might expect large flux strings to be absent (as the energy scales with the length of the flux string) and thus, one would not include any string-like term but rather only (nearest-neighbor) hopping-like terms in the construction of the variational form, as we will see in the following section. Finally, note that since this set actually contains all terms appearing in the Hamiltonian of the system, it is reasonable to expect [71] that these pieces are sufficient to represent the whole physical Hilbert space of states for a particular choice of the variational parameters.

V. RESULTS AND DISCUSSION
In the following, we specialize to the U (1) case of the Yang-Mills model (1), and study the phenomenon of string breaking in (1+1)-and (2+1)-dimensional lattice QED which is reminiscent of confinement in QCD.

A. Groundstate preparation
First, we aim to represent the groundstates in different (fermionic) mass regimes using the gauge invariant trial states |φ(θ) =Û(θ) |φ init . Here, the initial state |φ init is chosen as the gauge invariant bare vacuum state configuration with zero particles on the lattice sites. The The particle number N of the ground state is shown as a function of varying mass m for d = 1, resulting from an ideal VQE algorithm (red dots) and from a hardware simulation on a superconducting quantum device (blue triangles). The resources reduce to 5 qubits, 16 cnots and 3 variational parameters. We observe the string-breaking phenomenon in a form reminiscent of a phase transition from a particle-antiparticle pair (left inset) to a flux-tube configuration (right inset) at critical mass mc = 3. c) The electric flux Li for each lattice link i and the particle number N (lower-right inset) of the ground state are plotted as a function of varying mass m for d = 2, resulting from an ideal VQE algorithm. Only one exact curve (black) for the electric flux is shown since the ground state flux configurations on each link coincide. The string-breaking phase transition is found to occur at a critical mass mc = 2.5. In the small-mass regime, the ground state corresponds to a symmetric superposition of two particle-antiparticle pair configurations in which the flux-tube can pass along two distinct paths (upper-left inset) of the plaquette.
lattice configurations and open boundary conditions are assumed as sketched in Fig. 2 a). For our specific model we use Wilson fermions with r = a = 1, g = 5 − d and we adopt the Quantum Link Model approach [38] for the gauge fields with a minimal truncation value of S = 0.5 and a background electric field θ = 0.5 . This choice encompasses either a single positive or a vanishing electric flux on each link. Note that due to our specific choice of the boundary conditions, this truncation does not affect the quality of the groundstate predictions (see Appendix B 1). To approximate the ground state, we employ the VQE algorithm [43,69] on the circuit variables θ, for each value of the varying mass m (see Fig. 2).
Starting with the d = 1 dimensional case, we define a specific family of gauge-invariant unitaries bŷ which consists of hopping-like terms in order to mimic nearest-neighbor hopping dynamics, and further contains mass-like terms that essentially act as additional singlequbit post-rotations. The coupling matrices are defined as A(θ) = (0, 0), (θ, 0) and B(λ) = (0, 0), (0, λ) , and θ = {θ k x }, {λ x } denote the variational parameters. Note that the hopping terms assigned to neighboring sites commute and the products in Eq. (16) are thus well-defined.
In the d = 2 case, the variational form (16) is extended by string-like termsŜ(θ) of the form where k, l > 0 and S kl =ψ † x A θ k,l x Û (x,k)Û(x+k,l)ψx+k+l as to provide longer range explicit correlations. The need for such additional correlations is due to the groundstate configurations of the system (see left inset in Fig. 2 c)) where the distance of the particle-antiparticle pair stretches over two links, while hopping-like terms merely produce particle-antiparticle creation on neighboring sites. On the other hand, the squared string-like terms were added to produce the corresponding groundstates that appear in a superposition.
Finally, the variational operatorÛ(θ) is translated into a quantum circuit by means of a Trotter approxima-tion [55,56] with one single Trotter step, and by using the Jordan-Wigner fermion-to-qubit mapping and a logarithmic encoding of the gauge field operators [38].
The results of the variational simulation [72] are shown in Fig. 2 b) and c) for d = 1 and d = 2 dimensions, respectively. The best run (i.e., the one with the lowest groundstate energy E 0 ) out of five independent VQE optimization trials (each with randomly chosen initial parameters θ) is plotted. In fact, due to the small difference between the groundstate and the first excited state in the d = 2 case (as explained below), more optimization trials in the VQE algorithm were needed to obtain the true groundstate. We identify a change of behavior from a particle-antiparticle state to a single flux-string state by plotting the particle number N = N = x ψ xψx + I as a function of the mass parameter m. In d = 2 dimensions, we further plot the expectation value of the electric flux L = Ê (x,k) + θ for each link in Fig. 2 c). In this last case, for small m, the groundstate is a superposition of states where the flux is traversing the plaquette through two possible paths (see inset in Fig. 2 c). The degeneracy between the symmetric and antisymmetric combination of the two paths is broken by the presence of the small (as g 2 1/g 2 ) but sizable plaquette term in the Hamiltonian (see Eq. (1)).
We also perform, as a proof-of-principle, calculations on hardware for the d = 1 model using the ibmq vigo device provided by IBM Quantum (see Fig. 2 b). The resulting circuit (C1) features 5 qubits and 16 cnot gates with 3 variational parameters, see Appendix B 3. Crucially, this strategy allows us to retrieve a qualitatively better description of the problem using much fewer parameters and entangling gates compared to heuristic circuit ansätze (see Appendix A).

B. Real-time propagation
Finally, we investigate the ability of the proposed gauge-invariant quantum circuits to prepare states arising in real-time dynamics, which we will refer to as state representability.
We use the proposed ansatz to study the dynamics of the one-dimensional lattice Hamiltonian with three sites in the small mass regime m m c . The model parameters are the same as for the simulations shown in Fig. 2 b) with explicit mass values m = 1 and m c = 3.5. The trial states are constructed using the unitary (16). At each time step, we numerically maximize the overlap | φ(θ)|ψ(t) | 2 between the variational trial states |φ(θ) and the exact time-evolved state |ψ(t) obtained by using the matrix representation of the time propagator. In Fig. 3, we show how the trial states can capture the correct particle number oscillation characteristic of stringbreaking even though the state space spanned by the trial states does not cover the entire physical Hilbert space H phys .
This result is of particular relevance in view of the possible use of the proposed circuits in variational time evolution algorithms [73]. In fact, the circuit used in this example merely requires 32 cnot gates (a value that in the variation algorithm remains constant for every time t) and 5 variational parameters. On the other hand, the standard Trotterization approach requires a much larger number of cnot gates, as a single Trotter step would feature 216 cnots, while at least ∼ 10 Trotter steps would be needed to reach an accuracy [38] comparable to the one reported in Fig. 3. In such a low mass regime, it was indeed expected from an energetic point of view that particle-antiparticle pairs would mainly form within short distances (i.e., in neighboring lattice sites) and thus, the states formed via the unitary (16) would well approximate the involved states in the evolution. By going to a different (critical or large) mass regime however, one would need to adapt the trial states accordingly by extending the unitary (16) with string-like terms, for instance, in order to account for a better approximation of the trial states, while coming at the expense of a larger amount of gates in order to build up the unitary.
Future work will necessarily include an analysis of various possible combinations of the single components in the variational form toolbox in order to find the optimal trade-off between large circuits implementing the unitary and an effective sampling of the whole physical Hilbert space.

VI. CONCLUSIONS
We introduced a general methodology to construct parametrized families of quantum circuits that preserve the gauge symmetry for U (1) and Yang-Mills theories. When applied to ground-state search algorithms, the method allows for reliable variational calculations as the trial states are always constrained within the physical submanifold of the full Hilbert space, reducing therefore the number of variational parameters. The main advantages of realizing the Gauss law constraint at the circuit level, rather than as an energy penalty term, include: (i) the smoothness of the resulting energy landscape, which allows for faster optimizations, and (ii) a substantial decrease of the energy estimator variance [70,74] due to the absence of the (usually large) energy penalty term. In this work, we demonstrated the accuracy of the method in QED models featuring particles and fields as explicit degrees of freedom.
Going beyond ground-state simulations, we show that the proposed gauge invariant trial states can further efficiently represent quantum states originating from realtime evolution, realizing a significant reduction in circuit depth (gate count) compared to standard approaches based on the Trotterization of the time-evolution operator. We therefore anticipate the usage of the present quantum circuits in combination with variational realtime evolution algorithms [73] to reduce the computational costs associated with the simulation of the dynamics of SU (N ) gauge theories, such as quantum chromodynamics.
In quantum field theories, the VQE [43,69] can be exploited to study various groundstate properties for varying model parameters. For example, one might investigate phase diagrams of Abelian and non-Abelian gauge theories [9,10,63] in the regime of a finite fermion density, i.e. at an unbalance between fermionic matter and antimatter, which remain inaccessible to standard Monte Carlo simulation methods due to the notorious sign-problem.
When applying the VQE algorithm to a general Yang-Mills lattice gauge theory, a challenge arises in the search for an efficient ansatz for the variational form [75]. The reason for this challenge stems from the Gauss law constraint that determines the physical Hilbert subspace, H phys ⊂ H which makes up an exponentially small fraction of the total Hilbert space H. In our implementa-tion (that follows the framework of Ref. [38]), the Gauss law constraint is imposed by adding a gauge penalty term H gauge which raises the energy of the unphysical gaugevariant states. As a consequence, the gauge penalty term leads to a corrugated energy landscape [68] (Fig. 4) whose global minimum is in general hard to find for classical optimization algorithms, as was observed in [75]. In order to quantitatively confirm this hypothesis, we tested the performance of the so-called RY and RYRZ variational form on a small lattice QED example in d = 1 dimensions for a varying mass parameter m, which is mapped to a 5-qubit system. This example is introduced in detail in the next section.
The RY and RY RZ are two standard heuristic variational formsÛ(θ) that are constructed such that a wide range of states in a general Hilbert space can be covered without taking into account the underlying structure of the corresponding physical system. To do so, these variational forms systematically entangle all of the qubits and apply parametrized rotations R y for RY or R y and R z for RY RZ, respectively, on each single qubit.
The performance of these two variational forms is shown in Fig. 5. The relative energy difference ∆E = |E vqe − E 0 | / |E 0 | is plotted as a function of varying mass parameter m and for increasing depth of the RY and RY RZ variational forms. Each data point in the plot corresponds to a single VQE run (in statevector simulation) which achieved the lowest energy estimation E vqe out of a total of 10 independent optimization runs (each starting from different randomly chosen initial parameters) for each mass value. We used the Cobyla optimization algorithm for the classical optimization procedure. The initial parameters θ (0) were randomly chosen from a uniform distribution in [−π, π] and the initial state |φ init was given by the bare vacuum state, i.e. the unique gauge-invariant state with zero particles on each lattice site. For comparison, we further included the performance of the gauge-invariant variational form (16) with A(θ) = (0, 0), (θ, 0) and B(λ) = (0, 0), (0, λ) . Consequently, the matrices A couple only the second spinor component at a site x to the first spinor component at an adjacent site x + k. This implies that hopping terms assigned to neighboring sites commute and thus, the prod- ucts in Eq. (16) are well-defined. Furthermore, for each site that sits in a corner of the finite lattice, there is always one spinor component that is left untouched by this variational form, allowing to effectively remove the qubits that encode these spinor components [75]. The corresponding quantum circuit for a one-dimensional lattice system with two sites is shown in (C1). From the results in Fig. 5 it is visible that the standard heuristic variational forms struggle to find the correct groundstate energy, resulting in a poor approximation of E 0 . Raising the circuit depth merely aggravates the performance in spite of a higher number of degrees of freedom. The origin for such large errors lies in the large regularizing coefficient λ = 20 in the gauge penalty termĤ gauge that leads to high energy contributions in the order of ∼ 10 2 but which is needed to fix the correct lowlying gauge-invariant energy spectrum. Furthermore, it was observed that the performance of the standard variational forms greatly depended on the initial set of parameters θ (0) . Therefore, these standard variational form ansätze are not promising for larger system sizes and nonabelian gauge theories.
In contrast, the gauge-invariant variational form clearly outperforms the standard heuristic ansätze by one or two orders of magnitude. At the same time, it requires much less parameter degrees of freedom which is beneficial for the classical optimization procedure, even for larger sizes of the physical system.
We further tested the performance of the standard heuristic RY variational form on a d = 1 dimensional SU (2) gauge theory model where the Gauss law constraint is incorporated in the Hamiltonian and the gauge degrees of freedom are effectively erased [58]. Since the Hilbert space of this model only contains gauge invariant states, we expect the RY variational form to achieve better convergence to the groundstate energy.
The integrated SU (2) model was implemented on M = 2 lattice sites with lattice spacing a = 1, Wilson parameter r = 1, coupling constant g = 3 and boundary con-ditionsR α 0 = 0, ∀ α. The groundstate energy E 0 was examined as a function of varying (negative) mass parameter m. In the regime where m ≈ −r/a, the model exhibits a specific type of behavior reminiscent of a phase transition from a bare vacuum phase with no particles to a charge-crystal phase where all lattice sites are fully occupied [63]. Such a phase transition occurs when the prefactor m + r/a of the mass term nψ nψn in the system Hamiltonian changes its sign from positive to negative since in the former case, the production of particles costs energy while in the latter case, it becomes energetically favorable to create particles pairs.
For various mass values m, we ran a separate ideal VQE for each m using the RY variational form with linear entanglement, i.e. where the qubits are not all pairwise entangled but only neighboring qubits are entangled, leading to a reduced number of CN OT gates. We used the Cobyla optimization algorithm for the classical optimization procedure. The exact groundstate energy curve was calculated by exact diagonalization. The results are plotted in Fig. 6.
From these results, it is visible that the RY variational form shows a better performance with smaller relative errors ∆E = |E vqe − E 0 | / |E 0 | for the integrated SU (2) model than for the general scalable QED model with a gauge penalty term. However, the error increases in the regime of the phase transition, indicating that this heuristic variational form is not able to capture the essential properties of the groundstate in the region of the phase transition where the quantum fluctuations, arising from the hopping term in the system Hamiltonian, start to become significantly large. It has further been checked that changing the linear entanglement to a full entanglement or raising the depth of the variational form did not lead to significant improvements for the overall performance.
To conclude, this short analysis corroborates the need of physically motivated and structured variational forms which respect the symmetries of the Yang-Mills theory at hand. Furthermore, note that heuristic standard variational forms are usually constructed at the qubit level and thus will depend on the exact encoding scheme for the gauge and matter field degrees of freedom.
Appendix B: Explicit lattice QED example in 1D on real quantum hardware Here we present a detailed description of the (1+1)dimensional lattice QED system exhibiting the phenomenon of string breaking, which was simulated on a few-qubit superconducting quantum device.

The model
In the following, we consider a one dimensional lattice with M = 2 sites and one link with open boundary conditions. It is convenient to relabel the matter fields bŷ ψ x →ψ n , the flux operator byL (x,k) →L n and the link variable byÛ (x,k) →Û n . Recall that the QED Hamiltonian (1) for dimension d = 1 including the gauge penalty term readŝ In addition, we allow for a non-trivial constant background electric field θ to be added to the model [53], which is taken into account by shifting the electric flux field operator by a constant amount, (B2) where θ ≡ θ · I acts as an identity operator times a real constant. Therefore, such a background electric field effectively shifts the allowed electric flux eigenvalues m s by an amount of θ. For example, for a spin value S = 1.5 in the quantum link model formulation of the Abelian gauge fields and a background electric field θ = 0.5, the allowed flux eigenvalues are shifted by m s ∈ {−1.5, −0.5, 0.5, 1.5} → m s ∈ {−1, 0, 1, 2}. As a result, the flux eigenvalues are not centered around the zero flux value anymore, thereby breaking the spatial symmetry of the physical system. Nevertheless, for certain choices of the boundary conditions or for high enough spin values S, the low-lying energy spectrum will not be affected by this shift as will be the case in the model that we are considering. On the other hand, the background electric field will allow us to reduce the number qubits needed to store the relevant states of the model.
In our simulations, we choose the model parameters g = 4, a = 1, r = 1, λ = 20, the spin value S = 0.5 with a background electric field θ = 0.5, and boundary conditions as depicted in Fig. 7. A total of 5 qubits are needed to represent the total Hilbert space of this lattice system with spin value S = 0.5, namely two qubits per site and one qubit for the link in between, thus spanning a Hilbert space H of dimension 2 5 = 32. Furthermore, the background electric field θ = 0.5 shifts the two allowed flux eigenvalues by m s ∈ {−0.5, 0.5} → m s ∈ {0, 1}.
Note that out of these 32 states, only 5 states corresponds to physical states. The possible gauge invariant state configurations that are allowed by the Gauss law constraint with this choice of boundary conditions are shown in Fig. 9 and in Fig. 10.
From these configurations it is clearly visible that our choice of the small spin value S = 0.5 with θ = 0.5 does not affect the groundstate configurations of the model (B1): First, all gauge invariant states contain only fluxes in positive direction which is consistent with having possible flux values m s ∈ {0, 1}. Second, the only physical state that is not representable with this choice of flux eigenvalues is the one with a particle-antiparticle pair connected by two flux strings. Nevertheless, such a state will always correspond to a higher energy state configuration than the two states in Fig. 9 and therefore, will not occur as a groundstate of this model. As a conclusion, choosing the values S = 0.5, θ = 0.5 instead of S = 1, θ = 0 allowed us to effectively remove one qubit in order to represent the gauge fields on the link. Now let us analyze the groundstate for varying mass parameter m. First, it is clear that the only possible groundstates that are to be considered are those depicted in Fig. 9 since all other gauge invariant states shown in Fig. 10 will always have a higher energy independently of the mass value m. As a next step, assume that the Hamiltonian (B1) would not contain any hopping term. In that case, the Hamiltonian is already diagonal and the groundstate would depend on whether the mass term H mass or the electric termĤ elec has a smaller energy. The corresponding critical mass m c is simply calculated by equating the total energy of the two state configurations in Fig. 9, and thus, the groundstate energy E 0 as a function of mass m will have a discontinuity in its first derivative, see Fig. 8. Note that one might expect the groundstate energy for vanishing hopping term in the small mass regime m m c to be equal to zero since both the electric field term and the mass vanish for the gauge invariant particleantiparticle pair state. However, the small negative energy contribution that is visible in the left graph originates from the Wilson term which contains a hopping like term proportional to r and which was not set to zero.
By adding the non-diagonal hopping term, the groundstate energy in the limiting regimes m m c and m m c stays approximately the same while in the critical mass regime m ≈ m c , strong quantum fluctuations are introduced that smooth out the groundstate energy and the particle number curves withN = n ψ nψn + I) as is shown in Fig. 8. As a result, a phase transition is observed which corresponds to the flux string breaking.

Statistical sampling errors due to the gauge penalty term
A problem arises when evaluating the parametrized energy expectation value on a real quantum hardware since it can not be exactly calculated as in an ideal simulation, but the parametrized energy is rather statistically sampled from various measurements of the same circuit. This follows from the postulates of quantum mechanics where a measurement outcome of an observable randomly occurs with a probability determined by the wavefunction and thus, the outcome expectation value corresponds to an average of all the outcomes.
This fact can lead to problematic consequences in a VQE algorithm, especially due to the gauge penalty term H gauge with a large regularizing coefficient λ, since the statistical sampling errors might become too big leading to large non-zero energy contribution and thus, the correct groundstate energy can not be recovered anymore, although the gauge penalty term should in principle vanish for gauge invariant states.
In our lattice gauge theory models, the sampled expectation value φ(θ)|Ĥ gauge |φ(θ) of the gauge penalty term will in practice only approximately equal to zero for physical states. Since non-zero sampled expectation values will be proportional to the regularizing parameter λ which is usually chosen to be much larger than the other model parameters, the energy evaluations in a VQE algorithm will be largely disturbed.  Fig. 9. The state with a particle-antiparticle pair is not contained in the truncated state description for S = 0.5 with θ = 0.5.
To certify this hypothesis, we sampled the expectation value of the gauge penalty term λĤ gauge with λ = 20 with respect to an equal superposition of the two gauge invariant groundstates shown in Fig. 9 with increasing number of shots (i.e. measurements) by using IBM's Qiskit QASM -simulator [72]. Each weighted Pauli string in the gauge penalty term was measured separately. The results are plotted in Fig. 11. It is observed that the statistical fluctuations are of the same order of magnitude as the energy scale of the groundstate energies of the considered model (B1). Note that if one is interested in the groundstate only, the regularization value λ could in principle be chosen smaller (as long as it lifts the energy of the unphysical states to an amount which is bigger than the energy of the groundstate). However, higher excited states become important in a real-time evolution in general and thus, the regularization value must increase correspondingly to ensure that the dynamics remain physical [75].
It is further shown in Fig. 12 how such statistical sampling errors increase for a one-dimensional lattice with increasing number of lattice sites. We observe an enhancement of the standard deviation with increasing system size, hinting that such statistical noise might become intractable for LGT models with large lattices.
The observed issue can be circumvented in different ways. First, it might happen that the weighted Pauli strings P which make up the gauge penalty term H gauge = P λ P P all commute with each other, as is the case in a Jordan-Wigner fermion-to-qubit encoding of this term, for example. Then, the Pauli strings can all be simultaneously measured and therefore, the same state measurement-outcome statistics can be used to evaluate the expectation value of all the Pauli strings in a grouped manner [69,70,74], resulting in a vanishing statistical error. However, the exact form of the Pauli strings will depend on the specific encoding of the gauge and matter fields and the Pauli strings may not pairwise commute in general. Furthermore, the exact choice of grouping the Pauli string bases might affect the overall sampling efficiency of evaluation of the total Hamiltonian.
In our simulations we pursue a different strategy. Since we are working with a variational form that is gauge invariant by construction, there is in principle no need to additionally impose the Gauss law constraint as a penalty term at all, provided that the initial state |φ init used in the VQE algorithm corresponds to a physical gauge invariant state configuration. As a result, the variational  form samples only from the physical Hilbert space leading to the correct low-lying energy spectrum and thus, we will simply set the gauge penalty term to zero. It is observed that with vanishing gauge regularizing parameter λ = 0, the statistical fluctuations remain small compared to the energy scale set by the chosen model parameters. Nevertheless, in practice there might be a small probability that the gauge invariant variational form does sample from a set of unphysical states depending on the form of the coupling matrices A(θ), B(θ) in (16) due to Trotter approximation errors. The Trotter error in the approximation can be made arbitrarily small by increasing the number of Trotter steps at the cost of enlarging the circuit gate complexity. Note that the explicit gauge invariant variational form (C1) that is used in our simulations consists of four Pauli strings as a result of the Jordan-Wigner mapping which all pairwise commute. Thus, there is no Trotter approximation error in our simulations.
There exists a third approach as was already pointed out in [75] that makes use of so-called Gauss law oracles [62]. Such Gauss law oracles can be included in the VQE algorithm as a subroutine to check whether the produced variational form |φ(θ) is actually gauge invariant or not. In that way, one can make sure that the small error produced by the Trotterization does not lead to an unphysical energy spectrum even when the gauge penalty term is set to zero.
This observation will be exploited in our quantum simulation in order to reduce the required qubit resources from 5 to 3 qubits, thereby effectively lowering noise errors that are present in a simulation with real quantum hardware.
With the preliminary work in the previous paragraphs, we consequently ran a separate VQE algorithm in order to search the groundstate of the model for each distinct mass value m by using the specific variational form (16) with three parameters, and whose implementing circuit takes the form (C1). The initial state |φ init was chosen as the gauge invariant bare vacuum state with no particles on the sites and one unit of positive flux on the link.
The complete quantum circuit was run on the ibmq vigo superconducting machine which is characterized in Fig. 13. At each iteration step of the VQE algorithm, the parametrized circuit was repeated 8000 times (8000 shots) in order to achieve a small sampling error for the energy and particle number expectation value. For the classical optimization routine, the SPSA optimization algorithm (100 optimization steps, without calibration) was used which has been proposed for minimizing noisy energy functions [70]. Furthermore, we applied a measurement-error mitigation scheme that is included in IBM's Qiskit software-development framework in order to reduce noise errors due to the measurement readout process [72].
Appendix C: Gauge invariant variational form: Circuit example Here, we display the explicit quantum circuit representing the variational form (16) with 3 variational parameters (λ 0 , λ 1 , θ 0 ) that was used in the VQE simulations for the lattice QED model (B1) on a d = 1 dimensional lattice with M = 2 sites and with a spin value S = 1/2 in the quantum link model truncation of the gauge electric fields.
The first two X gates before the first barrier generate the gauge-invariant bare vacuum state with zero particles on the sites and positive flux on the link in between. The following variational form visibly consists of four Pauli strings I 0 X 1 Y 2 I 3 Y 4 , I 0 X 1 X 2 I 3 X 4 , I 0 Y 1 Y 2 I 3 X 4 , and I 0 Y 1 X 2 I 3 Y 4 which all commute with each other, thus implying a vanishing Trotter approximation error.