Variational quantum simulation of U(1) lattice gauge theories with qudit systems

Lattice gauge theories are fundamental to various fields, including particle physics, condensed matter, and quantum information theory. Recent progress in the control of quantum systems allows for studying Abelian lattice gauge theories in table-top experiments. However, several challenges remain, such as implementing dynamical fermions in higher spatial dimensions and magnetic field terms. Here, we map D-dimensional U(1) Abelian lattice gauge theories onto qudit systems with local interactions for arbitrary D. We propose a variational quantum simulation scheme for the qudit system with a local Hamiltonian, that can be implemented on a universal qudit quantum device as the one developed in [Nat. Phys. 18, 1053-1057 (2022)]. We describe how to implement the variational imaginary-time evolution protocol for ground state preparation as well as the variational real-time evolution protocol to simulate non-equilibrium physics on universal qudit quantum computers, supplemented with numerical simulations. Our proposal can serve as a way of simulating lattice gauge theories, particularly in higher spatial dimensions, with minimal resources, regarding both system sizes and gate count.


I. INTRODUCTION
Most quantum information processing platforms are based on qubits, the quantum generalization of classical bits.However, the underlying physical systems representing qubits frequently involve higher-dimensional Hilbert spaces that must be artificially restricted to twolevel systems.Instead of limiting it, however, one can use the Hilbert space that the physical system provides for information processing.This leads to the multi-level analog of the qubit -the qudit, which can be a powerful resource for quantum information processing [1].The additional levels can enable alternative implementations of quantum algorithms [2] the implementation of optimal quantum measurements [3], as well as the native simulation of higher spin models or problems in quantum chemistry [4].Moreover, the fundamentally different coherence [5], dissipation, and entanglement structure [6] of qudit systems can be advantageous in terms of noise resilience [7] or quantum error correction [8].These prospects of qudit systems and recent experimental progress make multi-level systems ideal for advanced quantum information processing.
So far, qudit experiments have been proposed [9] and extensively used in quantum cryptography [10] for increased information capacity and improved resilience to perturbations [7].Beyond photons, almost all quan- * pavel.popov@icfo.eutum technology platforms have demonstrated some degree of qudit control.More recently, superconducting systems [11], single photons [12] and trapped-ion experiments [13] have demonstrated a universal set of gates for qudit quantum computing.This rapid development of qudit hardware allows for the study of state-of-theart quantum algorithms such as quantum simulation on these novel devices.Originally driven by the goal of developing a large-scale quantum computer, quantum simulation has been identified also as an attractive target for devices of the so-called noisy intermediate-scale quantum era [14][15][16].In particular, the quantum simulation of lattice gauge theories (LGTs) has made spectacular progress over the last decade [17][18][19][20][21][22][23].
LGTs are many-body systems with important applications in high-energy physics, condensed matter systems, and quantum information.In high energy physics, LGTs appear as the space-discretized description of the Standard Model of particle physics [24].LGTs can also be found as an effective description in condensed matter physics and are important for quantum error correction, e.g., the toric code [25].This versatility of LGTs makes them a central object of research for different communities.Quantum simulation protocols for LGTs have been proposed for numerous quantum platforms, from cold atoms [26][27][28][29][30][31][32][33][34][35][36][37], through trapped ions [38][39][40][41][42] to superconducting qubits [43][44][45][46] and others [47][48][49].Indeed, first experimental implementations of LGT simulations are now a reality [50][51][52][53][54][55][56][57].Even though it is possible to quantum simulate LGTs in the laboratory, the experimental demonstrations remain constrained to specific scenarios.In particular, systems in one spatial dimension and some specific systems in two spatial dimensions with Abelian symmetry have been experimentally realized.However, the implementation of Abelian LGTs beyond one spatial dimension remains challenging due to the presence of dynamical fermions and magnetic field terms (fourbody interactions on a lattice) in the Hamiltonian of the theory [22].
Here, we propose a quantum simulation protocol for an Abelian U (1) LGT in (1+1) and in (2+1) spacetime dimensions for qudit quantum processors based on trapped ions [13].Following [58,59], we integrate out the fermionic fields and construct an LGT that can directly be mapped to a qudit system.The construction keeps the Hamiltonian of the theory local even in the case of higher spatial dimensions and evades the need of introducing Jordan-Wigner strings in order to encode the fermionic degrees of freedom in the qudits.In contrast to previous works that make use of the procedure of integrating out the fermions [60,61], we circumvent the need for the precise implementation of the resulting interactions in the Hamiltonian on the quantum device by proposing a hybrid quantum-classical variational simulation scheme for both ground-state preparation and quench dynamics.We also elaborate on which quantities have to be measured on the quantum device in order to execute the protocol and how to perform the measurement.We benchmark this approach by performing numerical simulations of the variational algorithm with the specific quantum circuits we propose for preparing the ground state of the LGT as well as simulating quench dynamics.As our results show, employing variational algorithms for static and dynamic properties of Abelian gauge theories is well within reach of current devices also for dimensions larger than (1+1).
This paper is organized as follows.In Section II, we briefly introduce the mathematical description of the target LGT and construct the corresponding qudit Hamiltonian.Then, in Section III, we give an outline of the variational algorithm we use.Here, we also explain how the relevant quantities for the execution of the algorithm can be measured in the experiment.Later, in Section IV, we present specific quantum circuits that need to be implemented on the quantum device and discuss how the simulation protocol would be implemented.In Section V, we present the results of the numerical study and discuss the findings.

II. U(1) LGT WITH FERMIONIC MATTER
We consider a D-dimensional spatial lattice with sites x ∈ Z D .The spinless fermionic matter is localized on the sites x where we denote the creation and annihilation operators by ψ † x and ψ x , respectively.The gauge degrees of freedom are situated on the links between neighboring sites x and y and act on a d-dimensional Hilbert space with the basis |0⟩ , . . ., |d − 1⟩.The electric field operator E x,i and the parallel transport (link) operator U x,i , where the latter keeps track of the phase generated by the gauge field, are defined as where the subscript x, i denotes the link starting from the lattice site x and pointing in direction i, i.e., the link between sites x and x + e i .For a D-dimensional spatial lattice, there are also D distinct directions.We used the notation ).Further, the fermionic degrees of freedom fulfill the algebra and the electric field and the link operator fulfill The dynamics of the quantum many-body system are determined by the Hamiltonian where the pure gauge part H G is given by with the plaquette term and the coupling constant g is given by the charge of the electron.The first term in Eq. ( 5) has the meaning of the electric field energy, whereas the second term, which involves the plaquette terms, has the meaning of discretized magnetic field energy.The fermionic part H M of the Hamiltonian is given by where we introduced the occupation number operator and the staggered charge M denotes the fermionic rest mass.The gauge-matter interaction is The system has local symmetries generated by the Gauss's law operator FIG. 1. Particle content on the lattice in the case of a two-dimensional spatial lattice.Left: The hardcore bosons ηx and the fermions χx reside on the lattice sites (orange circles), whereas the gauge fields Ex,i and the fermions ζx,i reside on the links (blue ovals).Right: The Majorana modes cx, αx,i and βx,i defined by Eq. ( 14) and the corresponding sites/links on which they reside.Thanks to Gauss's law, this configuration permits one to replace the fermionic matter ψx residing in the original gauge theory on lattice sites, such that the final theory contains only the gauge fields Ex,i.
Because of the fact that the Hamiltonian H and the Gauss's law operator G x commute, they have a common eigenvector system and we restrict the dynamics to the zero eigenvalues of each G x , i.e., Precisely these local constraints on the physical states between the matter and the gauge field degrees of freedom will allow us later to formulate a unitarily equivalent system without the matter.
A. U (1) LGT formulated as a qudit system In this subsection, following two previous works [58,59], we construct a gauge theory from a system of hardcore bosons, fermions, and gauge fields residing on a lattice.This construction introduces one fermionic and one hardcore bosonic degree of freedom per site, both constituting the matter of the LGT.The fermionic degree of freedom keeps track of the fermionic statistics and the hardcore bosonic degree of freedom allows to decouple the attachment of the gauge field to the matter.On top of that, on each link of the lattice, gauge field and fermionic degrees of freedom reside.Moreover, a gauge-invariant Hilbert space and a local Hamiltonian for these degrees of freedom are formulated, so that this Hamiltonian has the same matrix elements as the one of Eq. ( 4).The degrees of freedom of the constructed LGT are mapped onto a qudit system, while the locality of the interactions is preserved.This mapping allows for the efficient implementation of the LGT simulation in a qudit-based quantum device.In what follows, we present the details of this mapping.
We consider a Hilbert space consisting of three kinds of particles: χ x , η x , and ζ x,i where χ x and ζ x,i are fermionic operators and η x is a hardcore boson operator (see the defining algebra below).While χ x and η x live on the sites of the lattice, the fermion ζ x,i lives on the links, see FIG. 1. Specifically, these operators fulfill the following commutation and anti-commutation relations: [η x , η † y ] = [η x , η y ] = 0 for x ̸ = y and (13c) For notational convenience, we introduce the following Majorana operators for χ x and ζ x,i By having two modes per site (η x and c x ), and one fermionic mode on each link (ζ x,i ), we have more degrees of freedom than in the original formulation of the LGT.Therefore, we restrict the operators to act on wave functions that fulfill We construct the operator which is fermionic as can be confirmed by direct calculation.By gauging these fermions, i.e., by coupling them to a gauge field residing on the links of the lattice, as explained above, we can define a Hamiltonian of a LGT like the one in Eq. ( 4), but acting on the Hilbert space defined by Eq. ( 15).This Hamiltonian reads and, using the notation introduced above, Gauss's law is In the following, we perform two unitary transformations to map the dynamics to a completely electricfield-dependent operator being a pure qudit system.As shown in Ref. [58], for specific gauge groups, one can define a unitary transformation that replaces the fermionic FIG. 2. Two-dimensional spatial lattice with matter fields on the sites (orange circles) and gauge fields on the links (blue ovals): After performing the transformation from Eq. ( 19), the hopping terms between two neighboring sites (indicated by an arrow) acquire prefactors that depend on the electric fields on the adjacent links.The involved electric fields depend on the hopping direction (horizontal or vertical) and result in direction-dependent phase factors as given in Eq. (23).
operators in the Hamiltonian with operators that act locally on the gauge fields.This unitary transformation is given by the expression Since the different parts of the product do not commute with each other, we have to explicitly specify the ordering.The choice we made by writing Eq. ( 19) coincides with the one in Ref. [58].Transforming the operators of the gauge theory degrees of freedom yields We have used the following abbreviation for the prefactors of the transformed link operators: Moreover, the property [V x , V y ] = 0 allows us to define a unitary transformation V = x V x .This global unitary transforms H GM as where f x,i (E •,• ) is a prefactor that depends on the electric fields around the site x and is also dependent on the spatial direction i.Specifically, in (2+1) dimensions, i = 1, 2 for horizontal and for vertical hopping, respectively, and we have The electric fields involved in Eq. ( 23) are highlighted in FIG. 2. In (1+1) dimensions, since there is only one direction, we omit the index i and we have The outcome of this transformation is that the hopping term H GM does not include the Majorana modes c x anymore.
The plaquette term present in the Hamiltonian in higher spatial dimensions also transforms non-trivially: with the transformed plaquette operator defined as After the transformation, the occupation numbers of the fermionic mode χ x in the physical space are the same as the occupation numbers of the hardcore bosons η x .Furthermore, as easy to see from Eq. (20e), the modes ζ x,i transform identically, therefore they stay in vacuum.The Hamiltonian of the system projected onto the subspace of absent fermions ζ x,i reads The hardcore bosons on each site are furthermore related to the divergence of the gauge fields on the adjacent links due to the Gauss's law, Eq. ( 18).Precisely these local constraints allow for complete elimination of the matter from the Hamiltonian, as shown in Ref. [59].We can rewrite the Gauss's law also as where we introduced We introduce the operator and the unitary transformation for the entire lattice We denote physical states transformed by U by a tilde | ψ⟩ and define them as The matter transforms as Since for a vector of the physical Hilbert space, the operator g x has eigenvalues 0 or 1, we obtain This last equation signifies that all matter will be transformed to its vacuum state.Next, we apply the unitary transformation on the Hamiltonian of the gauge theory.In order to write down the transformed Hamiltonian in a convenient form, we define on each vertex x the operators P g,x , which projects on the g = 0 or g = 1 eigenspace of the operator g x .Applying Eq. ( 30) to the raising and the lowering operators in the interaction part of the Hamiltonian, we get precisely the projector P 1,x .Note, that the ordering of the link operators and the projectors is now important since they do not commute.We end up with the following Hamiltonian in (2+1) dimensions In (1+1) dimensions, the Hamiltonian simplifies to the form Since we have decoupled the matter by a unitary transformation, we are left with a system of qudits that reside on the links of the lattice, with dynamics governed by a local Hamiltonian.This formulation allows for an effi-cient implementation on a qudit quantum device.In summary, we achieved the following: We applied two transformations on the LGT.The first transformation (Eq.( 19)) removed the Majorana modes c x from the Hamiltonian and replaced them effectively with electricfield-dependent prefactors.On the Hilbert space level, it coupled the modes c x and η x so that physical states only contain configurations with equal occupation numbers for both c x and η x .This allowed us to keep the modes η x only.
The second transformation (Eq.( 30)), removed the hardcore bosons η x by using the constraints between mat-ter and gauge field degrees of freedom (Gauss's law).As indicated in Eq. ( 35), physical states after the transformation are the ones with all the matter fields in the vacuum state.Furthermore, the transformed Hamiltonian projects onto the subspace of physical states, making only the gauge fields have non-trivial dynamics.This allowed us to consider a system of gauge fields only.Combining both transformations, we arrived at local qudit Hamiltonian, given by Eq. (36) or Eq. ( 37) in (2+1) or (1+1) dimensions, respectively.
Later in this work, for the numerical studies and the explicit circuits, we will consider a system of qutritsqudits with three internal levels (d = 3).However, our approach is general and can be used (after modification of the quantum circuits) for general d-level qudits.

III. VARIATIONAL QUANTUM SIMULATION
The simulation of quantum many-body systems stands as a top application of quantum computers.However, the task of engineering intricate Hamiltonians poses a significant challenge.In the realm of analog quantum simulation, the construction of the Hamiltonian itself proves to be the challenge, whereas, in digital quantum simulation, the primary difficulties lie in efficient digitization and fault tolerance with minimal resources.To tackle these obstacles in both analog and digital quantum simulation, variational algorithms have emerged as particularly powerful and represent a potent tool for both classical and quantum systems.In this approach, a parametrized family of quantum states is efficiently prepared through a variational Ansatz circuit, which is typically based on the physical properties and symmetries of the system under study.Here, the classical computer is tasked with optimizing a cost function, such as the energy of the candidate states for a given Hamiltonian, while the quantum computer is responsible for measuring the cost function on the system.Through an iterative process, this variational approach can provide an efficient representation of quantum many-body states, such as ground states of a given Hamiltonian.
Recent technological advancements now permit the implementation and study of quantum variational algorithms.In the quantum domain, variational simulation finds frequent application in many-body physics, encompassing disciplines such as solid-state physics, highenergy physics, and chemistry [62].Notably, in the field of LGTs, there have been proposals and experimental implementations of simulation protocols based on the variational principle [41,45,52,[63][64][65].The principle has been used to create approximations to the ground state [45,52,66,67], thermal states [68], excited states [45,69], or the dynamics [70,71] of a quantum many-body system.

A. Parametrized quantum circuits
We consider trial quantum states generated by a parametrized quantum circuit (PQC) of the form where U k (θ k ) are unitary operations acting on the initial state |ψ 0 ⟩.These unitaries are parameterized by a real number θ k and which can be written as The generators g k can be chosen in a hardware-efficient fashion according to the native operations that can be performed on the specific quantum device, where the variational algorithm should be executed.Here, we consider the gate set used in the trapped-ion qudit quantum computer of Ref. [13], where the native gates are two-level single-qudit rotations and entangling Mølmer-Sørensen (MS) gates.In addition, we will consider qudit controlled-rotation (CROT) gates, which can be constructed from the MS gate [13].For more details on the implementation, see Section IV.
In the following, we discuss the variational principle for imaginary-time and real-time evolution, where the exact imaginary time evolution converges to the ground state of the quantum many-body systems at infinite imaginary time and the real-time evolution describes the change of the quantum many-body dynamics in time.In addition, we explain a strategy to determine the imaginary and real-time evolution from a measurement of the quantum device.

B. Variational time evolution
The equation describing the dynamics of a quantum state ρ of a system with Hamiltonian H is where the super-operator L[•] is given by for Hamiltonian dynamics and for imaginary time evolution.The expectation value in Eq. ( 42) is taken in the state ρ(t), leading to non-linear dynamics.
In the following, we parameterize the quantum state ρ(t) by a set of variational parameters θ(t), which are time-dependent on their own, and we can obtain the (approximate) time evolution by minimizing the so-called McLachlan distance, with respect to θ(t), i.e., trajectories of the variational parameters.The minimization leads to the equations of motion for the variational parameters The real symmetric matrix M µν represents the so-called Fubini-Study metric tensor [72] and is defined by Note that M µν does not contain information about the time evolution operator.In turn, the vector contains information about the time-evolution operator L[•], which is different for real-time Eq. ( 41) and for imaginary-time evolution Eq. ( 42).
In the following, we will assume that the variational state ρ(θ) is a pure state defined by a PQC: We will suppress the explicit dependence on θ for notational convenience.Inserting the pure state from Eq. ( 47) in the definition of the metric tensor Eq. ( 45) leads to For imaginary-time evolution, the vector V I µ , once we insert the pure state into Eq.( 46), is given by For real-time evolution, instead, we have The execution of the variational time evolution algorithm requires that the quantities M µν in Eq. ( 48), V I µ in Eq. ( 49) and V R µ in Eq. ( 50) are being measured on the quantum device.In the rest of this Section, we will give possible measurement schemes for these quantities.

C. Direct measurement protocol for Mµν and Vµ in the experiment
We proceed with the discussion on how to measure the quantities of Eq. ( 45) and Eq. ( 46) on the quantum device.First, we define the overlap functions The components of the metric tensor are related to the second derivatives of these functions as For the diagonal elements, since we consider pure states with the property we obtain Therefore, either measuring the derivatives of the functions in Eq. ( 51) directly on the quantum device or measuring the functions and then calculating the derivatives, will allow us to obtain the metric tensor.In Appendix A, we explain how to extract the derivatives by measuring the functions, using the parameter shift rules [72].The components of the vector V µ for imaginary-time evolution are also derivatives of observables that can be measured on the quantum device.From Eq. ( 49), it is evident that the first derivatives of the functions are directly connected to the components of V I µ : In Appendix B, we discuss how to extract V I from the values of the function v I µ (a).For real-time evolution, in contrast, V R µ cannot be represented as a derivative of an observable.However, we can make use of the representation of the unitary operations involved in the PQC in Eq. (39).Performing the derivatives with respect to the variational parameters in Eq. ( 50), we obtain where we defined using the notation The subscript c in Eq. ( 57) means that we consider the connected (anti-)commutator.One possible measuring procedure for such a quantity involves averaging over global random unitaries [73,74] |0⟩ 3. Hadamard test for measuring unequal-time (anti-)commutators: In order to measure the (anti-) commutator, the ancillary qubit needs to be initialized in an equal superposition of |0⟩ and |1⟩ with a relative phase of α = 0 ( π 2 ).We obtain the desired correlation function by measuring the probability of the ancilla to be in the state |+⟩.
Here, • denotes the averaging over random unitaries, drawn from the so-called circular unitary ensemble (see Ref. [73] for details), and N H is the Hilbert space dimension.Hence, instead of having to deal with the unequaltime anti-commutator of Eq. ( 57), we just have to measure expectation values of observables (⟨•⟩ u ) in states that result from the action of global random unitaries on the state to be measured and averaged over these unitaries.
Here, we used a generic operator H representing the whole Hamiltonian.However, the local structure of our particular Hamiltonian from Eq. (36) or Eq. ( 37) is of great use, since we can represent the anti-commutator from Eq. ( 57) as a sum of anti-commutators with each summand of the Hamiltonian.
One challenge related to the execution of this direct measurement protocol is the high number of repetitions that is needed for the reliable estimation of the matrix and vector elements M µν and V µ .This is due to the fact that these quantities contain derivatives of observables, the estimation of which requires the evaluation of the function at many nearby points.Remarkably, this difficulty can be circumvented by performing an alternative measurement procedure using an additional ancilla, as described below.

D. Indirect measurement of Mµν and Vµ through an ancilla
In this subsection, we give an alternative protocol for measuring the components of the matrix M µν and the vector V µ .This protocol is based on performing Hadamard tests on the quantum device in order to extract unequal-time (anti-)commutators.We have already shown in Eq. ( 57) that for real-time evolution, the components of V R µ can be written as unequal-time anticommutators of observables.It is also possible for the components of M µν and V I µ to be written as unequal-time (anti-)commutators.A simple calculation shows that for pure states, we have and We now present a protocol for indirect measurement of such qudit quantities through an ancilla.In the following, we briefly outline the main idea of the protocol and show how to use it for measuring M µν and V µ .In a separate work [75], we will present the algorithm in a more general framework, show how one can use it to extract unequal-time correlation functions of observables in the experiment and compare it to other such algorithms.
We observe that we can write a general Hermitian matrix H as a sum of two unitary matrices, i.e.U H + U † H . Once we insert this expression in the unequal-time (anti-)commutator, we obtain a sum of four quantities.We measure the unequal-time (anti-)commutator of two generators g k and g l with the circuits shown in FIG. 3. In these circuits we have to couple an ancillary qubit to the qudit register that represents our system, using a controlled operation that corresponds to each of the unitaries U H .These circuits are also known as the Hadamard test.In total, we can measure the (anti-)commutator by performing four Hadamard tests, one for each combination of unitaries.Further, we can measure anti-commutator or commutator, if we choose the initial phase α = 0 or π 2 , correspondingly.

E. Measuring the Hamiltonian
In this subsection, we rewrite the Hamiltonian in a form convenient for performing the measurement protocols from the previous subsection.For simplicity, we discuss the case of a (1+1)-dimensional theory and comment on the case of (2+1) dimensions.The first two terms in Eq. ( 37) are diagonal in the computational qudit basis, and thus can directly be measured by projectively reading out the qudit state.In addition, these terms are also local, so that the resulting controlled unitaries for the execution of the Hadamard test represent entangling operations between only the ancilla and the corresponding qudit.
The third term in Eq. ( 37) has to be rewritten as a tensor product of local observables in order to be measured projectively.In the case of qutrits, we have Here, we introduced the projector onto the i-th level of the qudit on site x as P i x .The matrices σ i,j X,x are two-level Pauli-X matrices embedded in the qudit Hilbert space; these are defined as Measuring this part of the Hamiltonian projectively would therefore involve a local change of bases.
In order to implement the measurement procedure involving Hadamard tests for this part of the Hamiltonian in a quantum device, we need to rewrite the different terms as a sum of unitaries.For a Pauli matrix acting on a two-level subspace of a qutrit, the corresponding unitary decomposition is given by where we defined We construct the unitary decomposition of the product of all projectors with The unitary decomposition of the first term in the Hamiltonian H GM is given by The decomposition of the whole Hamiltonian as a sum of local unitary operators has to be inserted in Eq. ( 57) and Eq. ( 62).Furthermore, each term of the resulting expression can be measured by the circuits given in FIG. 3.
In the case of the Hamiltonian in (2+1) spatial dimensions, the number of summands of H GM is larger and we need to perform entangling gates that correspond to the few-body product of projectors.This means some few-body-controlled operations are involved.Recent proposals for the experimental implementation of such gates have been made [76].Alternatively, we could write each projector as a sum of two unitaries; in this case, we would = FIG.4. Individual layers of the variational quantum circuit for seven qudits in (1+1) dimensions: Each building block of the parametrized quantum circuit used for generating the trial states is built from a set of single qudit operations, given by two-level rotations, and entangling operations, given by the two-level Mølmer-Sørensen [77,78] gate.
only have to perform two-body entangling operations at the expense of more Hadamard tests.
In summary, we have proposed two measurement protocols for the extraction of the geometric tensor M µν and of the vector V µ , which make the implementation of the variational time evolution on a qudit quantum device possible.

IV. TOWARDS IMPLEMENTATION ON QUDIT HARDWARE
In this section, we explain how an implementation of the time evolution protocol would look on qudit quantum hardware.We briefly discuss what solution strategy for the equations of motion of the variational parameters we = FIG. 5. Individual layers of the variational quantum circuit for one plaquette in (2+1) dimensions: Each building block of the parametrized quantum circuit used for generating the trial states is built from a set of single-qudit operations, given by two-level rotations, and entangling operations, given by the CROT gate (see Eq. ( 77) for definition).FIG. 6. Variational imaginary-time evolution of a system of seven qudits in (1+1) dimensions.Left: Fidelity F θ (τ ) of the variational state with respect to the exact ground state as a function of the imaginary time for different circuit depths (N = 1, 2, and 3 layers).All of the results converge after a finite time.The highest reached fidelities are ∼ 74%, ∼ 99%, and > 99% for N = 1, 2, and 3 layers, respectively.Right: Expectation value of the Hamiltonian of the system in the trial state as a function of the imaginary time.The error of the ground state energy in the final time of the simulation is below 1% for N = 3 layers.
could employ and we present details about parametrized quantum circuits that can be implemented on the qudit device from Ref. [13].

A. Solving the equation of motion for the variational parameters
For both variational real-and imaginary-time evolution, solving Eq. ( 44) would involve a feedback loop between the quantum device and a classical computer.In practice, one needs to solve a discretized version of Eq. (44).For example, a possible discretization is with t n = n∆t and where ∆t is the time step.Here, we have used Euler discretization, but other discretizations (Runge-Kutta, etc.) for better numerical stability are also possible.The quantities M µν (t n ) and V µ (t n ) have to be measured on the quantum device; then Eq. ( 72) for the values of the parameters for the next time step θ ν (t n+1 ) has to be solved classically, and then the obtained values have to be given to the quantum device for the next iteration.The whole procedure has to be repeated until the end of the simulation.

B. Ansatz for the variational quantum circuit
One important aspect of our work is the choice of the parametrized quantum circuit.On the one hand, we need an expressive ansatz that allows for the parametrization of the gauge invariant Hilbert space.That is, we would ideally like to have a parametrized quantum circuit that can reproduce every physical state of the system we want to simulate, for some set of values for the variational parameters.On the other hand, we aim at near-term realizability of the simulation protocol on qudit hardware.This naturally restricts the types of quantum gates available for implementation and the number of entangling operations we can perform since each entangling operation introduces errors.Furthermore, high expressivity can restrict the efficiency of the classical optimization of the variational parameters.For these reasons, we have to look for the middle ground between expressivity and realizability.
In the case of quench dynamics, the entanglement in the system may grow fast.This fact suggests that the expressivity of the quantum circuit for real-time evolution should be chosen high.Reversely, if the expressivity and in particular the entanglement that the ansatz allows, are not chosen correspondingly, the real-time evolution simulation may fail to reproduce the correct time evolved state already at early times.For ground states, in contrast, it is known for various systems that the entanglement structure follows an area law [79], and a small number of entangling operations as well as a small number of single qudit operations should be sufficient to deliver the desired approximation with high precision.
The depth of the circuit, as well as the number of variational parameters, can be reduced by exploiting the symmetries of the system.For example, the (1+1)dimensional Schwinger model has a CP symmetry that can be used to reduce the number of variational parameters by applying the same gates with the same variational parameters on various qudits in the system, as shown for qubits in Refs.[52,80].
For both variational imaginary-and real-time evolution, we employ a variational quantum circuit structured in layers.Each layer consists of a set of single qudit FIG. 7. Variational imaginary-time evolution of a system of one plaquette in (2+1) dimensions.Left: Fidelity F θ (τ ) of the variational state with respect to the exact ground state as a function of the imaginary time for different circuit depth (N = 1, 2, and 3 layers).All of the results reach high fidelity after a finite time, around ∼ 91%, ∼ 99%, and > 99% for N = 1, 2, and 3 layers, respectively.Right: Expectation value of the Hamiltonian of the system in the trial state as a function of the imaginary time.The error of the ground state energy in the final time of the simulation is below 1% for N = 3 layers.
operations and entangling two-qudit gates, see FIG. 4 and FIG. 5.Each gate in each layer of the circuit has a variational parameter given by the respective angle of rotation.We choose the set of single-qudit operations differently for imaginary-and for real-time evolution.In the case of imaginary-time evolution, we employ single qudit unitaries given by where the subscripts e, o, and r stand for even, odd, and edge links, respectively and mean that for some set of qudits, we apply gates with the same variational parameter (in accordance with the symmetry of the system).In the above unitary, the single qudit rotations are defined as with σ i,j φ = cos(φ)σ i,j X + sin(φ)σ i,j Y .For real-time evolution, we choose the single-qudit unitary operations more expressively as In each layer, after a set of single-qudit unitaries, we apply a set of two-qudit entangling operations.For most of cases, these entangling gates will be Mølmer-Sørensen gates Note that the set of single-qudit unitaries used for the simulation of real-time evolution, together with the entangling Mølmer-Sørensen gate, constitute a universal set of operations for the qudit device.
In the case of real-time evolution of a plaquette in (2+1) dimensions, we use CROT entangling gates, which for qutrits is defined as Due to the connectivity in the trapped-ion quantum device of Ref. [13], entangling MS and CROT gates can be implemented between any pair of qudits, thus not restricting the setup to nearest-neighbour entangling gates.

V. RESULTS
To benchmark the experimental realizability, we perform numerical simulations for the real-and imaginarytime evolution of an Abelian LGT in (1+1) and (2+1) dimensions.As proof of principle, we set the coupling constant g = 1 and the fermion mass to M = 0.1 and work with a number of qudits and gate sets feasible with current hardware [13].In terms of entangling gate count, for simulating a system with L links in (1+1) dimensions, we construct circuits that use L − 1 entangling gates per layer of the variational circuit, or L gates for the (2+1)dimensional case.Hence, for N layers, we need at most LN qudit entangling gates; see Appendix E.

A. Variational imaginary-time evolution
For simulating variational imaginary-time evolution, we employ a quantum circuit with variable number of layers N , where each individual layer looks like the one in FIG. 4 for a system in (1+1) dimensions and FIG. 5 for a system in (2+1) dimensions.The variational circuit is FIG.8. Variational real-time evolution of a system of five qudits in (1+1) dimensions.Left: Fidelity of the variational state with respect to the exactly time-evolved state (see Eq. ( 79)).The fidelity decreases with time as the correct time evolution of the initial state becomes harder and harder to approximate with fixed number of layers.Right: Real-time dynamics after a quench of the gauge invariant fermionic numbers in the middle of the one-dimensional system.The dashed line corresponds to the exact time evolution of the fermion numbers, and the colorful lines correspond to the approximate time evolution with different number of layers N .

A B
FIG. 9. Evolution of the entanglement entropy of subsystem A after a quench in a system of five qudits in (1+1) dimensions: After initial linear growth, the entropy S(ρA) saturates to a finite constant up to fluctuations due to finite size effects.As we increase the number of layers N in the variational circuits, we observe better and better approximations to this behavior of the entropy in the trial state.
applied on an initial state that is easy to implement, for example, the product state of all qudits initialized in the state |1⟩, |ψ 0 ⟩ = i |1⟩ i .For performing the imaginarytime evolution, we are interested in initializing random states captured by the parametrized circuit.To prepare such random states, we have to choose a random set of variational parameters and apply the corresponding circuit on the state |ψ 0 ⟩.By numerically solving the discretized equation of motion for the variational parameters during imaginary-time evolution, we obtain a set of values for the parameters for each time instance.For each such set of parameters, we calculate the fidelity of the variational state with the ex-act ground state and the energy of the variational state ⟨H⟩ θ .
In FIG. 6 and FIG. 7, these two quantities are plotted as a function of the imaginary time for a system of 7 qudits in (1+1) dimensions and for a plaquette in (2+1) dimensions, respectively.The different lines correspond to variational circuits with different number of layers N = 1, 2, and 3. We see that for all circuits, the fidelity and the energy converge.In addition, the quantities converge to the exact ones with the infidelity of less than 1% for three layers.Correspondingly, we can obtain the correct energy with a relative error below 1%.

B. Variational real-time evolution
With the variational time-evolution approach, we can also approximate quench dynamics.Initially, a system of 5 qudits is prepared in the product state |ψ 0 ⟩ = i |1⟩ i , which corresponds to the ground state of the Hamiltonian in Eq. ( 4) for infinite fermion rest mass M → ∞.Then, we apply the time evolution operator with respect to the Hamiltonian in Eq. ( 4).In FIG. 8, the fidelity of the trial state with respect to the ideally time-evolved state, is shown as a function of time.The fidelity can serve as a strict quantifier of how reliable the variational approach reproduces the time evolution of the state of the system.In FIG. 8 also the evolution in time of a gauge invariant observable-the fermion number ⟨n⟩ θ -is given for different circuit depths (number of layers N = 2, 3, and FIG.10.Variational real-time evolution of a system of one plaquette in (2+1) dimensions.Left: Fidelity of the variational state with respect to the exactly time-evolved state (see Eq. ( 79)).The fidelity for N = 2 and 3 layers decreases with time, while the decrease for N = 4 is much slower.This suggests that for a system of 1 plaquette, there might not be a need to gradually increase the circuit depth with time to approximate the time evolution reliably.Right: Real-time dynamics after a quench of the gauge invariant fermionic numbers in the middle of the one-dimensional system.The dashed line corresponds to the exact time evolution of the fermion numbers, and the colored lines correspond to the approximate time evolution with a different number of layers N .
FIG. 11.Evolution of the entanglement entropy of subsystem A after a quench in a system of one plaquette in (2+1) dimensions: We observe that the growth of entropy at very early times is not correct; however, there is a qualitative agreement with the exact diagonalization.This growth of entropy leads to small deviations from the exact result for the fermionic occupation numbers; however, these deviations do not grow with time.
4).The dashed black line represents the exact time evolution of the fermion numbers under the Hamiltonian of the gauge theory.We observe good qualitative agreement between the variationally evolved and the exactly evolved fermion numbers in the case of N = 3, 4 (blue lines).Furthermore, for 4 layers, we also obtain quantitatively high fidelity (> 80%) for times up to two oscillation periods of the fermion numbers.
A quantity that can give us insight about why the variational time evolution starts to deviate after some time is the entanglement (von Neumann) entropy.We introduce a cut in the system of qudits (for example, in the middle of the 1 dimensional chain) and denote the two resulting subsystems as A and B.Then, the bipartite entanglement entropy of subsystem A reads where the reduced density matrix of subsystem A is In FIG. 9, we plot this entropy as a function of time.
As we expect, initially, the entropy grows linearly with time.After that, it saturates and becomes constant, up to some fluctuations due to finite-size effects.We observe a correlation between the ability of the quantum circuit to capture this entropy growth up to the saturated value and its ability to reproduce the correct behavior of the time-evolved system.Therefore we need a circuit capable of generating sufficient entanglement for the correct simulation of the time evolution.This capability grows with the number of entangling operations in the circuit.However, we need to restrict this number to keep the experimental error low.We also simulate the variational real-time evolution for a system of 1 plaquette in (2+1) dimensions.In FIG. 10, the gauge-invariant fermion occupation numbers ⟨n⟩ θ in the lower left corner of the plaquette is shown as a function of time.Also, the fidelity with respect to the exactly time-evolved state is shown as a function of time.Quantitatively, the fidelity drops below 80% for early times, and the time evolution of the fermion occupation numbers slightly deviates from the exact one.However, qualitatively, the time evolution of ⟨n⟩ θ is captured very well.As we see from the fidelity plot, even for later times, the fidelity remains constant.Here again, we investigate the bipartite entanglement entropy (FIG.11) and establish the correlation between quantitative results and entropy growth at early times.In Appendix D, we show how to modify the variational circuit to give very good qualitative results throughout the simulation.The new element is an entangling gate motivated by the plaquette term present in the Hamiltonian of the gauge theory in (2+1) dimensions.Though currently not state of the art, there are considerable efforts towards constructing native multipartite gates in trapped-ion quantum computers [76,81].

VI. CONCLUSION
The core of the presented approach to the quantum simulation of Abelian LGTs is based on three ingredients: the hardware-efficient implementation of the LGT by integrating of redundant degrees of freedom, while preserving the locality of the Hamiltonian; the mapping of the resulting LGT on a qudit model that can be directly implemented on qudit quantum hardware; and the application of a variational quantum simulation protocol for ground state preparation as well as quench dynamics.
In particular, we showed how to encode a U (1) LGT with dynamical fermions in arbitrary dimensions in a qudit system that can be implemented on already existing trapped-ion quantum computers [13].In a proof-ofprinciple numerical simulation, we showed the power of the variational time evolution algorithm in finding the ground state of this LGT in (1+1) and in (2+1) dimensions.Using only gates available in present-day experimental platforms and a small number of entangling operations, we successfully found a good approximation of the ground state of the LGT for various system sizes, including a system in (2+1) dimensions.
Furthermore, we also showed that we can variationally simulate the real-time evolution of an LGT and provided a circuit that is implementable on present-day devices.Gauge invariant quantities such as fermion occupation numbers showed several periods of oscillating behavior, a clear improvement with respect to previous studies [50,56].We also simulated the real-time dynamics of an LGT in (2+1) dimension and showed an improvement over various oscillation periods by extending the circuit with a four-body entangling gate.
So far, we were primarily concerned with the case of Abelian LGTs.In the future, our methods can be extended to gauge theories with a non-Abelian symmetry, such as SU(2) and SU(3), relevant for the Standard model of particle physics, but also some discrete non-Abelian gauge groups such as the dihedral groups D n .The gauge field Hilbert space of the latter is naturally finite-dimensional and thus can be implemented without truncation in a qudit quantum device.By combining the resource-friendly variational approach with a qudit-based quantum computation strategy, we anticipate significant improvement in the quantum simulation of non-Abelian gauge theories in the near future.the author(s) only and do not necessarily reflect those of the European Union, European Commission, European Climate, Infrastructure and Environment Executive Agency (CINEA), nor any other granting authority.Neither the European Union nor any granting authority can be held responsible for them.
FIG. 12. Variational real-time evolution of a system of 1 plaquette in (2+1) dimensions with a circuit that includes a four-body gate.Left: Fidelity of the variational state with respect to the exactly time-evolved state (see Eq. ( 79)).For N = 4, the fidelity does not drop below 80% for the whole simulation period.Right: Real-time dynamics after a quench of the gauge invariant fermionic numbers in the bottom left corner site of the plaquette.We see virtually perfect agreement between the exactly diagonalized (dashed line) result and the result from the variational simulation for N = 4 layers.FIG. 13.Evolution of the entanglement entropy of subsystem A after a quench in a system of 1 plaquette in (2+1) dimensions: The initial growth of the entropy is captured very accurately for N = 2, 3, 4 layers.Furthermore, for N = 4 layers, the time evolution of the entropy is very well approximated over the whole simulation time.A simple calculation shows that every hermitian matrix S can be represented as a sum of two unitary matrices.We construct those unitaries as follows: • Define U S := V U D V † .
Then, the relation holds where we have an explicit construction of the matrix U S .
Appendix D: Variational real-time evolution with a four-body entangling gate Here, we propose how to minimally enhance the variational quantum circuit from FIG. 5, so that we obtain significantly better results for the real-time evolution in (2+1) dimensions.The Hamiltonian of the system includes a four-body interaction, given by the plaquette term.If we include such a gate and assign a variational parameter to it U plaq (θ) = exp −i( ŨP + Ũ † P )θ , (D1) we can extend each layer of the parametrized quantum circuit by one such unitary.Then, performing the variational real-time evolution, we obtain the results shown in FIG. 12 and FIG. 13.

A B
FIG. 15.Evolution of the entanglement entropy of subsystem A after a quench in a system of 7 qudits in (1+1) dimensions: After initial linear growth, the entropy S(ρA) saturates to a finite constant up to fluctuations due to finite size effects.As we increase the number of layers N in the variational circuits, we observe better and better approximation to this behaviour of the entropy in the trial state.
Appendix E: Variational real-time evolution for 7 qudits in (1+1) dimensions In Section V, we studied the variational real-time evolution of 5 qudits in (1+1) dimensions and observed qualitative agreement for the period of two oscillations (see FIG. 8).Here, we complement the study of the variational real-time for a system of 7 qudits.We calculate the fidelity of the variational state at each time step concerning the exact state, see FIG. 14.In addition, we calculate the dynamics of the fermion numbers in the center of the one-dimensional system.We observe that the variational circuit for N = 3 and N = 4 layers can capture the correct behavior of the fermion number for up to one oscillation period.However, the fidelity drops significantly at later times, and the time evolution is no longer accurate.To understand this limitation, we investigate the entanglement entropy in the system of 7 qudits, shown in FIG. 15 and observe a restriction in the growth of the entanglement entropy.The limitation of entangle-ment entropy suggests that the circuits cannot generate the correct form of entanglement to approximate the exact real-time evolution for times later than one oscillation period.In future studies, we will improve this circuit by going to deeper circuits but also by choosing different entangling gates.
Appendix B: How to measure Vµ for imaginary time evolutionSince the components of the vector V I µ from the main text are connected to the derivatives of the function v I µ (a) (see Eq. (55)), we show here how to calculate these derivatives by measuring the function itself.Similar to f µν (a), the function v I µ (a) is of the formv I µ (a) = ⟨ψ µ |Q † µ (a)R † HRQ µ (a)|ψ µ ⟩ ,(B1)with unitary operators R and Q µ (a).The latter is of the formQ µ (a) = e ia qµ (B2)with hermitian q µ .Using the eigenvalue equation forQ µ (a) |q µ,m ⟩ = e iqµ,ma |q µ,m ⟩ , (B3)we can rewrite the function v I µ (a) asv I µ (a) = m,n v I µ,mn e ia(qµ,m−qµ,n) .(B4)There are R µ,v different values of q µ,m − q µ,n , which we call χ µ,rv , allowing us to rewrite the sum in Eq. (B4) asv I µ (a) = rv v I µ,rv e iaχµ,r v .(B5)In order to determine the Fourier coefficients v µ,rv we evaluate the function v µ (a) at R µ,v different values of a and determine the Fourier coefficients v µ,rv by solving the linear systemv I µ (a µ,r ) = rv v I µ,rv e iaµ,rχµ,r v .(B6)Given the analytical form of v I µ (a) the derivative is given by χ µ,rv e iaχµ,r v (B7)Having this derivative, the component of the vector V I µ are obtained by Eq. (56).

Appendix C :
Representing a Hermitian matrix as a sum of two unitaries

FIG. 14 .
FIG.14.Variational real-time evolution of a system of 7 qudits in (1+1) dimensions : (a) Fidelity of the variational state |ψ(θ(t))⟩ with respect to the exactly time-evolved state |ψ(t)⟩.(b) Real-time dynamics after a quench of the gauge invariant fermionic numbers in the middle of the one dimensional system.The dashed line corresponds to the exact time evolution of the fermion numbers and the colorful lines correspond to the approximate time evolution with different number of layers N .