Towards simulating 2D effects in lattice gauge theories on a quantum computer

Gauge theories are the most successful theories for describing nature at its fundamental level, but obtaining analytical or numerical solutions often remains a challenge. We propose an experimental quantum simulation scheme to study ground state properties in two-dimensional quantum electrodynamics (2D QED) using existing quantum technology. The proposal builds on a formulation of lattice gauge theories as effective spin models in arXiv:2006.14160, which reduces the number of qubits needed by eliminating redundant degrees of freedom and by using an efficient truncation scheme for the gauge fields. The latter endows our proposal with the perspective to take a well-controlled continuum limit. Our protocols allow in principle scaling up to large lattices and offer the perspective to connect the lattice simulation to low energy observable quantities, e.g. the hadron spectrum, in the continuum theory. By including both dynamical matter and a non-minimal gauge field truncation, we provide the novel opportunity to observe 2D effects on present-day quantum hardware. More specifically, we present two Variational Quantum Eigensolver (VQE) based protocols for the study of magnetic field effects, and for taking an important first step towards computing the running coupling of QED. For both instances, we include variational quantum circuits for qubit-based hardware, which we explicitly apply to trapped ion quantum computers. We simulate the proposed VQE experiments classically to calculate the required measurement budget under realistic conditions. While this feasibility analysis is done for trapped ions, our approach can be easily adapted to other platforms. The techniques presented here, combined with advancements in quantum hardware pave the way for reaching beyond the capabilities of classical simulations by extending our framework to include fermionic potentials or topological terms.


I. INTRODUCTION
Quantum computers are at the edge of pushing our computational capabilities beyond the boundaries set by classical machines [1][2][3]. One of the fields where quantum computers are particularly promising is the simulation of gauge theories, which describe the interactions of elementary particles. While current numerical simulations have led to several breakthroughs [4][5][6], they are ultimately restricted in their predictive capabilities. On one side, limitations originate from the inherent difficulty faced by classical computers in simulating quantum properties. On the other, the sign problem [7,8] affects simulations of both equilibrium (e.g., phase diagrams) and nonequilibrium (e.g., real-time dynamics) physics. Therefore, classical simulations based on tensor networks [9][10][11][12] or Markov chain Monte Carlo (MCMC) methods face hard numerical problems [13,14]. Quantum simulations are a promising solution, with proof-of-concept demonstrations in one-dimensional (1D) [15] gauge theories [16][17][18][19][20][21][22][23][24][25][26] already being achieved. Extending quantum simulations of fundamental-particle interactions to higher spatial dimensions represents an enormous scientific opportunity to address open questions that lie beyond the capabilities of classical computations. However, the step from one to two (or higher) spatial dimensions is extremely challenging due to an inherent increase of complexity of the models. Moreover, the current limitations of so-called noisy intermediate-scale quantum (NISQ) devices [27] make this leap even more difficult. In this work, we develop a protocol that ease the requirements on the quantum hardware and, as a consequence, allows for quantum simulations of gauge theories in higher spatial dimensions.
Current quantum simulation protocols include analog, digital, and variational schemes [21,31,32]. Analog protocols [33][34][35][36][37][38][39][40][41][42][43][44][45] aim at implementing the Hamiltonian of the simulated theory directly on quantum hardware. While this approach is challenged by the current practical difficulty of implementing gauge invariance and the required many-body terms in the laboratory, the first promising proof-of-concept demonstrations have been realized [16,[23][24][25]. Noteworthy are approaches based on ultracold atoms, which allow for large system sizes and have the advantage that fermions can be used to simulate matter fields. Digital protocols [17,18,20,[46][47][48] face similar practical challenges but have the advantage of being universally programmable and allow for the simulation of both real-time dynamics and equilibrium physics. Lastly, hybrid quantum-classical variational approaches [19,22,49,50] are in an early development stage and can be used to address equilibrium phenomena. These schemes do not require the simulated Hamiltonian to be realized on the quantum device. Along with their inherent robustness to imperfections [51], this feature makes variational schemes suitable for NISQ technology.
As shown in Fig. 1, we use a variational approach to simulate ground-state properties of LGTs. In contrast to previous schemes, our protocols provide the novel opportunity to use existing quantum hardware [22,52] to simulate 2D effects in LGTs, including dynamical matter and nonminimal gauge-field truncations, with a perspective to go to the continuum limit. We consider QED, the gauge theory describing charged particles interacting through electromagnetic fields. In contrast to 1D QED, where the gauge field can be eliminated [17,18,53,54], in higher dimensions both appear to be nontrivial magnetic field effects and the Fermi statistics of the matter fields become important. As a consequence, many-body terms appear in the Hamiltonian and implementations on currently available quantum hardware become challenging. In this work, we outline novel approaches to overcome these difficulties and to render near-term demonstrations possible.
Specifically, we provide effective simulation techniques for simulating quasi-2D and 2D lattice-QED systems with open and periodic boundary conditions. To address the problem of efficiently finding the ground state of these models on NISQ hardware, we develop a variational quantum eigensolver (VQE) algorithm [51] for current qubitbased quantum computers. In the quest of simulating LGT employing this VQE algorithm, we address the crucial points of: (1) developing a formulation of the problem within the resources of NISQ devices, (2) implementing the model efficiently on the quantum hardware, (3) having a clear procedure to scale up to larger more complex systems, and (4) verifying the results in known parameter regimes.
For the first step, we cast the model into an effective Hamiltonian description, as done in Ref. [28]. The total Hilbert space is then reduced to a smaller gauge-invariant subspace by eliminating redundant degrees of freedom (at the cost of introducing nonlocal interactions). In order to measure this effective Hamiltonian on the quantum hardware (second step), we find an encoding for translating fermionic and gauge operators into qubit operators. The quantum circuits for the VQE are subsequently determined, respecting the symmetries of both the encoding and the Hamiltonian. From one side, this allows for an optimal exploration of the subsector of the Hilbert space in which gauge-invariant states lie. From the other, our circuits are Hamiltonian inspired and can be scaled up to bigger systems and to less harsh truncations (third step). The fourth and last step in the above list is taken care of by comparing the outcomes of the VQE algorithm with analytical results, in parameter regimes that are accessible to both. In the following, we resort to perturbation theory and exact diagonalization but, in principle, more advanced techniques [55][56][57] can be used.
Our work allows us to run quantum simulations of 2D LGTs on lattices with arbitrary values of the lattice spacing a. With the results in Ref. [28], it is now possible to reach a well-controlled continuum limit a → 0 while avoiding the problem of autocorrelations inherent in MCMC methods. This offers the exciting long-term perspective of computing physical (i.e., continuum) quantities, such as bound-state spectra, nonperturbative matrix elements, and form factors that can be related to collider and low-energy experiments. The study of these physical observables requires large lattices, the simulation of which is inaccessible to present-day quantum devices. However, various local quantities describe fundamental properties of a theory and can be simulated on small lattices that are accessible today. An example is the plaquette expectation value, as used in the pioneering work of Creutz [58].
Despite its simplicity, this observable can be related to the renormalized running coupling that we consider below. We emphasize that the methods presented here can be used as a launch pad to further developments in the realm of quantum simulations of LGTs that are aiming at higher dimensions, larger lattice sizes, and/or non-Abelian theories. Furthermore, the LGT models can be extended to include topological terms, or finite fermionic chemical potentials, that are currently hard or even inaccessible to MCMC methods due to the sign problem.
We treat the case of open boundary conditions within a ladder system [see Fig. 2(a)]. Although this system does not encompass the physics of the full 2D plane, it allows us to explore magnetic field properties on currently available quantum hardware. For that reason, we provide a simulation protocol for the basic building block of 2D LGTs, a single plaquette including matter, that demonstrates dynamically generated gauge fields. This is shown by observing the effect of particle-antiparticle pair creations on the magnetic field energy. In particular, both positive and negative fermion masses are considered. In . Using Gauss's law, the number of gauge degrees of freedom can be reduced to one per plaquette (blue ellipses). (b) A single plaquette with matter sites at the vertices and gauge fields on the links. The circles (squares) represent odd (even) sites and black (gray) corresponds to unoccupied (occupied) fermionic fields. The positive field direction is to the right and up. (c) A table showing the mapping between fermionic sites, particles (e) and antiparticles (p), and spins. (d) The conventions for Gauss's law. (e) The two gauge-field configurations that minimize the electric field energy for a single plaquette with two particle-antiparticle pairs. As explained in the main text, these configurations are relevant for the 2D effects discussed in Sec. V A (for more details, see Appendix D). 030334-3 the latter case, MCMC methods cannot be applied due to the zero-mode problem [7,8].
Additionally, we consider a single plaquette with periodic boundary conditions (see Fig. 3) to demonstrate that our approach provides an important first step toward calculating the so-called "running coupling" in gauge theories [4,5]. The running of the coupling, i.e., the dependence of the charge on the energy scale on which it is probed, is fundamental to gauge theories and is absent in 1D QED. For example, its precise determination in quantum chromodynamics is crucial for analyzing particle-collider experiments. Here, we propose a first proof-of-concept quantum simulation of the running coupling for pure gauge QED.
The paper is structured as follows (see Fig. 1). In Sec. II, we introduce lattice QED and present the effective Hamiltonian description on which the proposal is based. The VQE algorithm is then outlined in Sec. III, along with a short description of the available quantum hardware on which it may be implemented, along with possible classical optimization routines. In Sec. IV, we demonstrate our VQE algorithm for an ion-based quantum computer. The all-to-all connectivity available on this platform is an excellent match for our approach in the NISQ era, since the elimination of redundant degrees of freedom results in nonlocal interactions. We present a detailed experimental proposal, along with a classical simulation of the proposed experiments, to demonstrate that an implementation with present-day quantum computers is feasible. Finally, we describe the physical 2D phenomena that we aim to study in Sec. V. Our conclusions and an outlook are presented in Sec. VI.

II. SIMULATED MODELS
In this section, we present the models to be simulated in our proposal. In Sec. II A, we review the Hamiltonian formulation of lattice QED in two dimensions along with the truncation applied to gauge degrees of freedom. The specific systems considered in the rest of this paper are then described in Sec. II B [open boundary conditions (OBCs)] and Sec. II C [periodic boundary conditions (PBCs)].

A. Lattice QED in two dimensions
In this work, we consider 2D lattices, with matter and gauge fields defined on the vertices and on the links, respectively. Using staggered fermions [29], electrons and positrons are represented by single-component fermionic field operatorsφ i for each site i. As shown in Figs. 2(b) and 2(c), odd-(even-) numbered lattice sites hold particles (antiparticles) that carry a +1 (−1) charge q i . The gauge fields on the links between sites i and j are described by the operatorsÊ ij (electric fields) and U ij . The electric field operators take integer eigenvalues e ij = 0, ±1, ±2, . . . withÊ ij |e ij = e ij |e ij , whileÛ ij acts as a lowering operator on electric field eigenstates, i.e., Using the Kogut-Susskind formulation [29], the Hamiltonian consists of an electric, a magnetic, a mass, and a kinetic term;Ĥ =Ĥ E +Ĥ B +Ĥ m +Ĥ kin , wherê In the summations, we use i + −→ j to denote the link between lattice sites i and j with positive orientation [see Figs. 2(a) and 2(b)]. We denote the bare coupling by g, m is the fermion mass, a is the lattice spacing, and is the kinetic strength. We use natural units = c = 1 and all operators in Eqs. (1) are dimensionless [28].
In the above Hamiltonian, we introduce the operator P n =Û † ijÛ † jkÛ ilÛlk , where sites (i, j , k, l) form a closed loop clockwise around the plaquette n as in Fig. 2(b) [59]. This allows us to define the plaquette operator where N is the number of plaquettes. At each vertex i, gauge invariance is imposed by the symmetry generators [29,33] i is the charge operator. Relevant-i.e., gauge-invariant-quantum states are defined by Gauss's lawĜ i | phys = i | phys for each vertex of the lattice, where the eigenvalue i corresponds to the static charge at site i. We consider the case i = 0∀i, but background charges can easily be included in the derivations below. In the continuum limit a → 0 and three spatial dimensions, Gauss's law takes the familiar form ∇Ê =ρ, whereρ is the charge density.
Gauss's law can be used to lower the requirements for a quantum simulation. More specifically, within the Hamiltonian formulation of a gauge theory, only an exponentially small part of the whole Hilbert space consists of gauge-invariant states. In QED, this subspace is selected by applying Gauss's law with a specific choice of the static charges. As a result, it is possible to eliminate some of the gauge fields and obtain an effective Hamiltonian that is constrained to the chosen subspace (see also Ref. [60]). Practically, this reduces the number of qubits that are required to simulate the gauge fields-which is particularly important in the NISQ era. As can be seen in the following sections, the trade-off of this resource-efficient encoding is the appearance of long-range many-body interactions in the effective Hamiltonian. Nevertheless, universal quantum computers offer the possibility to perform gates between arbitrary pairs of qubits (long-range interactions). Moreover, the VQE approach is promising when dealing with interactions that cannot be directly implemented in the quantum hardware. Therefore, resorting to the effective Hamiltonians of the models considered for our quantum simulation in Sec. IV is advantageous in terms of experimental feasibility of our protocols.

B. Effective Hamiltonian for open boundary conditions
In this section, we consider a 1D ladder of plaquettes with OBCs [see Fig. 2(a)]. Although this system does not encapsulate the full 2D physics of QED, it allows us to study important aspects of gauge theories that are not present in one spatial dimension, such as magnetic phenomena.
The OBCs for a ladder of N plaquettes appear as dashed lines in Fig. 2(a), indicating a null background field. While there are 3N + 1 gauge fields in the ladder, the application of Gauss's law to all vertices reduces the independent gauge degrees of freedom to N . Given the freedom in choosing the independent gauge fields, we select the (2n, 2n + 1) links for each plaquette n, as shown by the blue ellipses in Fig. 2(a). On each link that does not hold an independent gauge field, the corresponding unitary oper-atorÛ is set to the identity. As a result, the plaquette operator becomesP n =Û 2n,2n+1 and both the kinetic and magnetic terms are consequently simplified. The effective Hamiltonian for the ladder is derived in Appendix A.
We now focus on the basic building block of the 2D ladder system, the plaquette [see Fig. 2 The notation for describing the state of the plaquette with OBCs is a tensor product of two kets, with the first ket representing the matter sites 1 − 4 (v, e, and p are the vacuum, particle, and antiparticle, respectively) and the second ket the state of the gauge field on the We remark that in 1D QED, all gauge fields can be eliminated for systems with OBCs. As a result, only two-body terms remain in the Hamiltonian [17,18,53,54]. For the plaquette, many-body terms are unavoidable since gauge degrees of freedom survive. This represents one of the main challenges in simulating LGTs in more than one spatial dimension.

C. Effective Hamiltonian for periodic boundary conditions
For the second model, we consider a square lattice with PBCs and without fermionic matter. This system has been considered in Ref. [28] but it is beneficial to summarize here the main results that are required for the VQE simulation. As explained in Sec. II A, Gauss's law can be used to eliminate redundant gauge fields, resulting in an unconstrained effective Hamiltonian. The basic building block, a single plaquette with PBCs, includes eight gauge fields [see Fig. 3(a)] and is equivalent to an infinite 2D lattice of four distinct plaquettes [see Fig. 3(b)]. Due to the absence of matter, we consider the pure gauge Hamiltonian H gauge =Ĥ E +Ĥ B , given by Eqs. (1a) and (1b).
Since we are interested in studying ground-state properties, it is sufficient to consider three independent gauge degrees of freedom (see Ref. [28]), called "rotators"R i , i ∈ {1, 2, 3}, shown in Fig. 3(b). Rotators are a convenient representation for the gauge fields, as they correspond to the circulating electric field in a specific plaquette. As such, the operatorsP n (P † n ) in Eqs. (1b) and (2) are directly their descending (ascending) operators. Indeed, the commutation relations between a rotatorR and its associated operatorP are the same as those of an electric fieldÊ and its lowering operatorÛ, presented in Sec. II A. The notation for describing the state of the plaquette with PBCs is thus a tensor product of three kets, each corresponding to one of the three rotators represented by full lines in Fig. 3(b).
Usually, LGT Hamiltonians are formulated in the basis of the electric field, whereĤ E inĤ gauge is diagonal. For large values of the bare coupling g, the termĤ E is dominant and this representation is efficient. For small g, The representation in terms of independent gaugeinvariant (i.e., physical) operators called "rotators." As explained in Ref. [28], the plaquette can be seen as an infinite lattice of plaquettes. Moreover, to explore the ground-state properties, only three independent rotators,R 1 ,R 2 , andR 3 (shown in solid blue), are sufficient for describing the system. however, the magnetic termĤ B is dominant and the eigenstates ofĤ gauge are superpositions of all electric field basis states. Since the Hilbert space of these operators is of infinite dimension, this necessarily leads to truncation errors when the Hamiltonian is mapped to any quantum device (as in Sec. V A). To mitigate the truncation error, we apply the novel technique introduced in Ref. [28] and resort to two different formulations of the Hamiltonian. The electric formulation in Eqs. (1) is applied in the region g −2 1, while for g −2 1, a so-called magnetic basis is used, for whichĤ B is diagonal (see also Ref. [61]). This magnetic representation is obtained using the cyclic group Z 2L+1 (L ∈ N) to approximate the continuous group U(1) associated to each rotator. We first take an approximation in the form of a series expansion that is exact for L → ∞. Then, we apply the Fourier transform and truncate the group Z 2L+1 by only considering 2l + 1 < 2L + 1 states in order to achieve equivalence to a truncated U(1) theory. As such, the basis of the rotator becomes |−l , . . . , |0 , . . . , |l . This leaves us with L as a free parameter, the optimal choice of which, along with the technical details of approximating the U(1) theory with a truncated Z 2L+1 model, is studied in Ref. [28].
We thus arrive at the HamiltonianĤ and we denote the magnetic representation with γ = b, The coefficients f c k and f s k come from the Fourier transform and their analytical form is where ψ k (·) is the kth polygamma function [28]. We note that making use of electric and magnetic representations in the case of the plaquette with OBCs (see Sec. II B) is straightforward. For clarity, however, we use this formulation only when considering PBCs.

III. FOUNDATIONS FOR VARIATIONAL QUANTUM SIMULATIONS
As detailed in Sec. II, the Hamiltonians associated with the considered models contain exotic long-range interactions and many-body terms that are beyond the capabilities of current analog quantum simulations. This is particularly true once the infinite-dimensional operatorsÊ and U, corresponding to the gauge field and its connection [29], respectively, are encoded in terms of qubit operators (see Sec. IV A). Resorting to digital approaches, Trotterized adiabatic state preparation is possible in theory [62]. However, the high complexity of LGTs beyond one dimension renders this approach infeasible for the technology available today. Due to the absence of large-scale and universal fault-tolerant quantum computers, we thus employ a hybrid quantum-classical strategy, using a VQE protocol that renders the observation of the effects described in Sec. V within reach of present-day quantum devices. In fact, those protocols have already proven their validity in quantum chemistry [51,[63][64][65], nuclear physics [49,50], and even classical applications [66][67][68].
In the following, we first introduce the VQE protocol (Sec. III A) and discuss feasibility on currently available quantum platforms (Sec. III B). Then, we describe our minimization algorithm, which is run on the classical device to find the desired state (Sec. III C).

A. Hybrid quantum-classical simulations
The VQE is a closed feedback loop between a quantum device and an optimization algorithm performed by a classical computer [51,65]. The quantum device is queried for a cost function C(θ), which is expressed in terms of a parametrized variational state | (θ) = U(θ )| in . Here, | in is an initial state that can easily be prepared, while U(θ ) entails the application of a quantum circuit involving the variational parameters θ . The classical optimization algorithm minimizes C(θ) and provides an updated set of variational parameters θ after each iteration of the feedback loop. Since our quantum simulations aim at preparing the ground state (see Sec. V), our VQE cost function is defined as the expectation value of the HamiltonianĤ of the system with respect to the variational state, i.e., Importantly, the Hamiltonian is never physically realized on the quantum hardware. Instead, its expectation value is estimated by measuringĤ with respect to the variational state | (θ) [51,65] (see Appendix B). Thus, the VQE is advantageous whenever the studied model contains complicated long-range many-body interactions (as in this case) that cannot be realized in current devices. Moreover, this approach is insensitive to several systematic errors, such as offsets in the variational parameters θ, since these are compensated by the classical optimization routine. Both these properties relax the quantum resource requirements and facilitate future experimental implementations of the VQE routines presented in the remainder of this paper.

B. Quantum hardware considerations
In this work, we consider qubit-based platforms for the realization of a proof-of-principle experiment with present-day technology. Promising quantum computing platforms include configurable Rydberg arrays [69,70], superconducting architectures [1,71], and trapped ions [72,73]. Comparing these approaches, Rydberg-based systems offer the advantage that 2D and three-dimensional (3D) arrays can be realized using optical tweezers, which translates to a large number of available qubits. The ability to implement large qubit registers will make future generations of Rydberg arrays a very promising candidate for our schemes, once higher levels of controllability and gate fidelity become available. Superconducting and ion-based architectures already offer both high controllability and gate fidelities today. For superconducting qubits, entangling gates are inherently of nearest-nearest-neighbor type, which implies that entangling operations between nonneighboring qubits have to be realized through a number of swap gates. While this is entirely possible, especially for next-generation devices, our models entail long-range interactions that result from the elimination of redundant gauge fields (see Sec. II B and C) and thus suggest the necessity of creating highly entangled states, which would require a significant gate overhead on superconducting platforms. In contrast, ion-based quantum computers [74,75] have all-to-all connectivity, allowing for addressed entangling gates between arbitrary qubits. This aligns well with our target models, motivating their use for this proposal. Despite their comparatively slow readout, which limits the measurement budget, we solve this issue with both an efficient measurement strategy (see Appendix B) and an optimized classical algorithm (see Sec. III C). A realistic budget for currently available ion-based quantum computers is given by about 10 7 measurement shots. This includes the time for the quantum computer to initialize, modify, and measure the states and the time for the classical computer to minimize the cost function C(θ).
Here, we describe the properties of ion-based computers in more detail. We consider a string of ions confined in a macroscopic linear Paul trap [72,73]. The qubit states |0 and |1 are encoded in the electronic states of a single ion and can be manipulated using laser light in the visible spectrum. A universal set of quantum operations is available, which can be combined to implement arbitrary unitary operations. More specifically, the available gate set consists of local rotationsÛ j (φ) = exp[−i(φ/2)σ α j ] and addressed Mølmer-Sørensen (MS) gates [76] between arbitrary pairs of qubits. These gates are implemented with fidelities exceeding 98% for both single-and multiqubit gates [74,75,77]. While currently available ion-based hardware provides a sufficient number of qubits to carry our proposed protocols, future LGT quantum simulations addressing larger lattices will require large-scale quantum devices. Ion-based quantum computers can be scaled up using segmented 2D traps [78] and networking approaches that connect several traps together [79,80] and therefore offer a pathway for developing quantum simulations of increasing size and complexity.

C. Classical optimization routine
The classical optimization routine employed for the VQE needs to be chosen depending on the requirements from both the hardware and the cost function C(θ). Indeed, the stochastic nature of the latter has to be taken into account and the experimental repetition rate poses limitations on the number of data points that the classical machine can use for minimization. In the following, we consider different optimization routines, discuss their strengths and weaknesses, and motivate our specific choice.
Algorithms based on the scheme of gradient descent are promising candidates, particularly when the number of 030334-7 variational parameters grows large. Recently, a variety of these algorithms have been modified to take the stochastic nature of the cost function into account [81,82] and there are versions that are tailored to the quantum regime [83]. Methods to measure the gradient of an operator [84] have been proposed but they are costly when applied to multiqubit gates and require additional ancilla qubits [85]. In addition, finite-step-size approximations of the gradient require a large number of experimental evaluations to obtain reliable values for the gradient.
A limitation of the available quantum hardware is the small number of experimental shots, particularly when resorting to ion-based platforms. This drawback is especially limiting for gradient-based optimization. Indeed, to prevent the VQE to settle in a local minimum, we have to perform multiple runs from different initial conditions. Furthermore, already acquired data cannot easily be recycled in scenarios where the algorithm is applied to different parameter regimes (e.g., variations in g in Sec. IV). A family of algorithms that can recycle previous data is the mesh-based algorithms. These evaluate the cost function on a grid in parameter space, with refinements where the global minimum is expected to be [86]. The mesh and the values of the cost function can be stored and reused for a different set of parameters, to accelerate the optimization. The feature of building a data repository is shared with Bayesian optimizers, which build a regression model of the energy landscape based on Gaussian processes [87]. These algorithms explore the space of the variational parameters by performing cost-function evaluations at the location of potential minimizing points, suggested by their regression model. Hence, they are able to reduce the required number of measurements substantially. However, their effectiveness is limited to around 20 variational parameters [88].
Here, we employ an optimization algorithm similar to the one used in Ref. [22], which is a modified version of the dividing-rectangles algorithm (DIRECT) [89][90][91][92][93]. This algorithm divides the search space into so-called hypercells. Each hypercell contains a single sample point representative of the cost-function value in that cell. The algorithm selects promising hypercells, to be divided into smaller cells, based on the cost-function value as well as the cell size. Larger cells contain more unexplored territory and are hence statistically more likely to harbor the global minimum. During the optimization, the algorithm maintains a regression model as used in Bayesian optimization. This metamodel is used for an accurate functionvalue estimation that aids in selecting the hypercells to be divided. Furthermore, at regular intervals, one or more direct Bayesian optimization steps are carried out [88]. Our choice is motivated from the small number of measurement shots that can be performed on an ion-based quantum computer (see Sec. III A), and from the limited number of variational parameters θ required in our models (see Sec. IV).

IV. QUANTUM SIMULATION OF 2D LGTs
In this section, we present the VQE protocol for simulating the two models introduced in Sec. II. We give numerical results from a classical simulation, which includes the projection-noise error. Both proposed experiments prepare the ground state of the theory and measure the ground-state expectation value of the plaquette operator ∼ Ĥ B , as defined in Sec. II A. For OBCs, this allows us to study the dynamical generation of gauge fields by pair-creation processes, while in the case of PBCs it is related to the renormalization of the coupling at different energy scales (see Sec. V).

A. Encoding for quantum hardware
To run a VQE, we first require an encoding of the models outlined in Sec. II for qubit-based quantum hardware. While we follow the Jordan-Wigner transformation [94] for mapping the fermionic operators, there is not a unique procedure to represent the gauge-field operators (alternatives are given in Refs. [47,95]). Here, the encoding is chosen to reduce the complexity of the quantum circuits and to respect symmetries that both suit the simulated models and the quantum platform.
Gauge-field operators are defined on infinite-dimensional Hilbert spaces and the gauge field takes the values 0, ±1, ±2, . . . (see Sec. II A). To simulate them using finite-dimensional quantum systems, a truncation scheme is required. Let us take 2l + 1 basis states into account, i.e., the gauge field can take at most the values ±l. Consequently, we substitute the electric field operatorÊ with This opens up two different paths to implementation of the truncation forÛ. The first employs the spin lowerinĝ S − =Ŝ x − iŜ y and raisingŜ + =Ŝ x + iŜ y operators, such thatÛ →Ŝ − /|l|. The second prescribeŝ For |l| → ∞, both mappings forÛ ensureÛ †Û = 1 and the correct commutation relations betweenÊ and U. The errors introduced by finite |l| have been studied in Refs. [28,96], where it is proven that they are negligible in most scenarios [97]. The specific choice of the truncation scheme depends on the quantum hardware that is employed. As an example, for quantum systems that allow for the implementation of interacting spin chains with large spins, such as l = 1 [98], the definition based on the spin loweringŜ − and raisingŜ + is ideal. For qubit-based quantum hardware, it is convenient to employ Eq. (9) for the representation ofÛ, as done in the following. For a given l, each gauge degree of freedom is then described by the 2l + 1 states |e , e = −l, −l + 1, . . . , 0, . . . , l − 1, l. We map this vector space onto 2l + 1 qubits using where 0 ≤ j ≤ 2l. With this encoding, the gauge-field operators might be replaced by the simple formŝ i are the Pauli operators associated with the ith qubit. From these mappings, we directly recover the relationsÊ|e = e|e and U|e = (1 − δ e,−l )|e − 1 for all −l ≤ e ≤ l. As an example, for l = 1 the states in the gauge-field basis become For our protocol, the required resources scale linearly in terms of both the parameter l and the number of gauge and matter fields. The encoding presented in Eq. (10) requires 2l + 1 qubits for storing 2l + 1 states. In principle, the same information can be stored in log(2l + 1) qubits [47,95]. The reasons for which we choose the qubit encoding in Eq. (10) are the following. First, as demonstrated in Ref. [28], small values of l allow for a good description of the untruncated model. Second, the terms assembling the operatorÛ in Eq. (11b) are easily implementable on different hardware platforms and are even native in trapped-ion systems in the form of MS gates. Finally, the states in Eq. (10) form a subspace of fixed magnetization, which is decoherence free under the action of correlated noise (e.g., a globally fluctuating magnetic field) and allows for detection of single-qubit bit-flip errors. Importantly, the latter is separately true for each individual gauge field and for the fermionic state. Since all the utilized quantum states lie in the single-excitation subspace, a measurement resulting in states outside of this subspace can only be due to errors occurring during the quantum evolution, as long as the applied gates conserve the excitations. Hence, erroneous outcomes can be discriminated from faithful ones.
In the following, we use the term "physical states" to refer to the computationally relevant qubit states. For the gauge fields, this means that the qubit states lie in the computational space spanned by the states in Eqs. (12). For the matter fields the computational states of which are obtained via the Jordan-Wigner transformation, we operate in the zero-charge subsector (see Sec. II A), which translates into matter states of zero magnetization. As a final remark, we highlight that our encoding for the electric field operatorsÊ andÛ applies equally well to the rotatorR and plaquetteP operator used for the plaquette with PBCs (see Sec. II C).

B. Open boundary conditions: dynamically generated magnetic fields
In the case of a plaquette with OBCs, the simulation involves four qubits for the matter fields and 2l + 1 qubits for the gauge field. Here, we consider l = 1 and thus the system consists of seven qubits. According to Fig. 4(c), we number the matter qubits as 1, 2, 3, and 4, and the gauge qubits as 5, 6, and 7. Plugging the encoding presented in Sec. IV A into the Hamiltonian in Eq. (3), we obtain The matter qubit states are also given in Fig. 2(c). Recall that we choose the encoding such that physical states have total magnetization Ŝ z tot = 1. The VQE quantum circuit shown in Fig. 4(a) preserves not only the total magnetization of the system but also the magnetization of each of the gauge and matter subsystems. Hence, as mentioned in Sec. IV A, our magnetizationpreserving quantum circuit used in combination with physical input states confines the VQE to the space of physical states.
The VQE circuit in its unreduced form [i.e., including the gray-shaded gates in Fig. 4(a)] is motivated by the 030334-9 form of the Hamiltonian in Eqs. (3). All qubits are initialized in the input state |0 and NOT gates prepare the bare vacuum |vvvv |0 [see Fig. 2(c)] as the initial state for the VQE. For the gauge-field subsystem, the application of the parametrized iSWAP gates allows for accessing all three gauge-field states |1 , |0 , and |−1 , ensuring that the free ground state for = 0 can be produced. Similarly, the parametrized iSWAP gates on the qubits 1-4 are used to allow for all physical basis states within the matter subsystem and resemble the hopping terms of the kinetic Hamiltonian in Eq. (13d). These gates correspond to particle-antiparticle pair creation or annihilation in the model and, as a consequence, all matter basis states in the zero-charge subsector are made available by this part of the circuit. The kinetic Hamiltonian is likewise responsible for the entanglement between the subsystems, as the pair creation and annihilation processes are combined with a correction of the gauge field. The layer of parametrized controlled-iSWAP gates hence takes the role of the annihilation operatorÛ and entangles the matter and gauge subsystems. In the effective Hamiltonian of Eqs. (3), the gauge degree of freedom lies on the (2, 3) link and is directly coupled to matter sites 2 and 3. Accordingly, the circuit couples the gauge field with only these two fermions, which act as controls in the layer of controlled-iSWAP gates. Finally, parametrized iSWAP gates are applied on the matter qubits again to adjust the state after entangling the two subsystems and a layer of single-qubit z rotations is utilized to correct for relative phases. Other single-qubit operations are avoided, as they are generally not magnetization preserving. We highlight that this circuit (as well as those in Fig. 6), being Hamiltonian inspired, from one side ensures the capability of exploring the subsector of physical states in the Hilbert space. From the other, it avoids redundant gates and limits the circuit depth, ensuring that barren plateaus [99,100] in the energy landscape are prevented [101] when increasing the truncation l and/or the number of plaquettes.
The quantum circuit described above involves a total of 21 variational parameters. While not strictly necessary, it is beneficial for currently available quantum hardware to reduce the number of variational parameters. This can be done by identifying the intrinsic symmetries of the ground state. By classically simulating the circuit in Fig. 4, we find that the solution space is still accessible by fixing θ 16 = θ 19 = θ 21 = π/2 and removing the gates shaded in gray in Fig. 4(a), leaving a total of 11 parameters. Furthermore, θ 11 , θ 12 and θ 13 can be set to zero outside the transition region 2 g −2 5, within which an almost vanishing energy gap between the ground and first excited state requires more precision for correctly estimating the ground state.
Since the circuit design is based on the structure of the Hamiltonian, the same design principles can be applied to larger-scale systems. When adding more plaquettes, additional parametric iSWAP gates are used to populate all basis states within the matter and gauge subsystems. Then, the gauge degrees of freedom are coupled to their respective neighboring matter sites using additional controlled-iSWAP gates [in correspondence with Eq. (A1e)]. Finally, z rotations adjust the relative phases of all qubits. When increasing the truncation l, additional iSWAP gates are inserted for populating the newly introduced gauge-field states and controlled-iSWAP gates are added for entangling them with the respective matter sites. In both cases-adding more plaquettes and increasing the truncation cut-off-a linear increase in the number of qubits and iSWAP gates and a quadratic increase in the number of controlled-iSWAP operations is expected.
A classical simulation of the proposed experiment, including statistical noise on the cost function C(θ) (representative of the probabilistic nature of quantum state measurements-see Appendix B), is shown in Fig. 5. The data points (blue and black dots) correspond to the lowest energies found by the VQE. The energies and the plaquette expectation values are calculated using the exact state corresponding to the optimal variational parameters found for each value of g −2 . We verify that the VQE resorts to statistical errors affecting C(θ) that are always lower than the energy gaps of the ground and first excited states. The black solid lines are obtained via exact diagonalization of the Hamiltonian in Eqs. (13). Using the measurement procedure described in Appendix B and taking statistical error 030334-10 All dots are calculated using the exact state corresponding to the optimal variational parameters found by the VQE. into account, the entire plot corresponds to approximately 10 7 measurements to be performed on the quantum device. Half of this budget is used for the point at g −2 2.18, where the energy difference between the ground and first excited states is much smaller if compared to other values of g −2 (see above).
The ground-state energy found by the VQE approximates well the energy of the exact ground state, as shown in Fig. 5(a). The plaquette expectation value is sensitive to small changes in the variational state, which leads to relatively large deviations with respect to the results obtained via exact diagonalization in Fig. 5(b), even for states the energy of which is very close to the exact energy. Yet, the variational optimization is able to accurately resolve transitions in the order parameter. The fidelity of the variational ground state with respect to the exact ground state is particularly high in the extremal regions, exceeding 98%. All points achieve a fidelity greater than 90%.

C. Periodic boundary conditions: running coupling
In this section, we provide a VQE protocol for simulating the running coupling in LGTs. As explained in Sec. V B, the running coupling is a genuine 2D effect that can be studied experimentally in a proof-of-concept demonstration by first preparing the ground state of a plaquette with PBCs and subsequently measuring the expectation value . Differently from the plaquette with OPC, the electric and magnetic representations of the Hamiltonian [see Eqs. (14) and (15)] are used for different regions of the bare coupling g −2 . This is done to obtain better convergence to the untruncated result, as the two representations are well suited for the strong-and weak-coupling regimes, respectively (see Sec. II C and Ref. [28]). Noting that the definitions in Sec. IV A presented for the electric gauge field and their lowering operators are trivially extended to rotators and plaquette operators, we encode the plaquette with PBCs into nine qubits. Rotator 1 is represented by qubits 1 through 3, rotator 2 by qubits 4 through 6, and rotator 3 by qubits 7 through 9 (see Fig. 6). Thus, the Hamiltonian

030334-11
becomeŝ in case g −2 1, and whenever g −2 1. As for the case of OBCs, the circuit design for a plaquette with PBCs is motivated by the structure of the Hamiltonian and employs the same gate set. Due to the differences between the electric and magnetic representations, we use two different VQE circuits, which are shown in Fig. 6. Contrary to the plaquette with OBCs, the controlled iSWAP gates are not used to entangle the matter with the gauge subsystems. Instead, they are motivated by the coupling of the rotators and plaquette operators in the last terms of Eqs. (4a) and (4b), respectively, and the corresponding ones in the magnetic representation.
The circuit for the electric representation of Eqs. (14) is shown in Fig. 6(a). All qubits are initialized in the input state |0 and NOT gates prepare the vacuum state |0 = |010 for each of the three rotators as the initial state for the VQE. The layer of controlled-iSWAP gates reflects the coupling between the rotators in the electric HamiltonianĤ (e) . This term results from the elimination of rotator 4 as a redundant degree of freedom and introduces an asymmetry between rotator 2 and rotators 1 and 3. When increasing g −2 , the ground state spreads from |0 |0 |0 (in the rotator basis) to all other electric levels and states in which all three rotators have the same sign (|1 |1 |1 and |−1 |−1 |−1 ) receive the strongest negative contribution. To encourage the VQE to prepare the correct superposition of states, the parametrized controlled-iSWAP gates are connected to control the spread of population within rotators 1 and 3 based on the population of rotator 2.
The circuit for the magnetic representation of Eqs. (15) is shown in Fig. 6(b). Its construction is similar to the circuit for the electric representation described above but, rather, encourages the flip-flop interactions between rotators described byĤ (b) E [see Eq. (5a)]. Both the electric and magnetic circuits maintain constant magnetization of each gauge field, which prevents access to unphysical states.
Designing the VQE circuit based on the form of the Hamiltonian allows for a scalable architecture. For systems with additional plaquettes, the coupling between rotators remains pairwise, which translates into the addition of controlled-iSWAP gates between all pairs of coupled gauge fields, as described above for OBCs. When considering larger truncations l, additional iSWAP gates are introduced to allow for all gauge-field basis states. Additional controlled-iSWAP gates are then added to share entanglement in a similar fashion as for the case l = 1 considered here. In both cases, the scaling is the same as for the OBC circuit, i.e., the number of qubits and the number of iSWAP gates scale linearly, while the number of controlled-iSWAP gates scales quadratically in the worst-case scenario.
We simulate the proposed experiment classically, including statistical noise on the cost function C(θ). Our  Fig. 6. The red and blue data points are obtained from variational minimization with a total finite measurement budget of 6 × 10 5 measurements for the entire plot. The electric (magnetic) representation is shown in blue (red) and the black solid lines are determined via exact diagonalization of the Hamiltonians in Eqs. (4) and (5). (a) The energy of the variational ground state using the electric representation in the region g −2 < 1 and using the magnetic representation in the region g −2 > 1. (b) The plaquette expectation value as a function of g −2 . All dots are calculated using the exact state corresponding to the optimal variational parameters found by the VQE. results are shown in Fig. 7. Points obtained with the electric and the magnetic representation of the Hamiltonian are shown in blue and red, respectively, while the black solid lines come from exact diagonalization of the Hamiltonians. As for the OBCs, the energy and the plaquette expectation value are calculated using the exact state obtained with the optimal variational parameters found by the VQE. Using the measurement procedure described in Appendix B and taking statistical errors into account, the entire plot corresponds to 6 × 10 5 measurements to be performed on the quantum device.
The VQE protocol reaches the correct ground-state energy [see Fig. 7(a)] and the expectation value of the plaquette operator is accurate if compared to the exact truncated results [ Fig. 7(b)]. The fidelity of the variational ground state with respect to the exact ground state exceeds 96% for all points-and for the majority of points, it exceeds 99%.

D. Errors in ion platforms
In this section, we discuss the effects of experimental imperfections, as well as statistical noise, on our protocols. Since the results in Secs. IV B and IV C are derived assuming an ion-based quantum computer [77,[102][103][104][105][106], we consider here the main sources of errors in such platforms. Our considerations are based on the experimental apparatus in Refs. [22,107], which has been previously used for LGT simulations in 1D. For further information about using a superconducting device, see Refs. [1,71,108].
In the following, we distinguish between extrinsic and intrinsic errors. Extrinsic errors are determined by the characteristics of the experimental apparatus and include imperfect gate operations, dephasing, and systematic errors such as offsets in the variational parameters θ. By contrast, intrinsic errors are inherent to the VQE protocol and as such unavoidable. These are statistical errors and they follow from the stochastic nature of quantum measurements (see Appendix B).
As discussed in Refs. [22,107], for low circuit depths, intrinsic errors are the dominant source of noise. Since the cost function C(θ) is estimated from independent contributions (see Appendix B), we add up the uncertainties due to each term for separately. This substantially increases the uncertainty of C(θ), even for eigenstates of the considered Hamiltonian. As an example, for each value of the cost function, obtaining a statistical error smaller than the energy gap between the ground and first excited state requires averaging several thousands of experimental shots. This is expected to be the main, and largest, overhead from noise sources in our proposal.
Extrinsic errors have limited impact for the small system sizes considered here. Indeed, VQE schemes are robust against systematic over-or under-rotation in one-and two-qubit gates [51]. This follows from the fact that the classical optimizer determines the variational parameters θ by minimizing the cost function. As such, since the shifts in θ are slowly varying if compared to the time required for the cost-function minimization, these errors are compensated by the classical subroutine. Furthermore, imperfect gates and the finite coherence time limit the number of gates and the circuit depth that can be used, respectively. This results in a limitation in the number of plaquettes and truncation l that can be studied. Such a limitation can, however, be relaxed by technological improvements.

V. ANALYTICAL INTERPRETATION OF THE RESULTS
In this section, we use the effective Hamiltonian description provided in Sec. II to study lattice QED with OBCs (Sec. V A) and PBCs (Sec. V B). While next-generation quantum computers will widen the scope of our approach, we focus on phenomena that can be studied with small system sizes and for which quantum simulations can be carried out on current quantum hardware [1,17,52,[69][70][71]. Here, we motivate the relevance of the results of the quantum simulations in Sec. IV.

A. Dynamical matter and magnetic fields
We describe a minimal example that allows for the study of 2D effects in lattice QED with OBCs. More specifically, we examine the appearance of dynamically generated magnetic fields due to particle-antiparticle creation processes within a single plaquette. This effect manifests itself in an abrupt change of the ground state as a result of the competition between the kinetic and magnetic terms in the QED Hamiltonian of Eqs. (3). There are two parameter regimes in which this phenomenon occurs. In one, the parameter [see Eq. (3d)] is dominant if compared to the mass and the bare coupling g, allowing the kinetic and the magnetic terms to be the leading contributions to the energy. Alternatively, the shift of the ground state appears when the mass is smaller than zero (or with positive mass and nonzero background field). This last scenario is hard to simulate using MCMC methods, in which the use of the inverse of the lattice Dirac operator in combination with a negative mass induces zero modes [109], leading to unstable simulations. We note that in 1D QED, considering that a negative fermion mass is equivalent to a theory with positive mass in the presence of a topological term θ = π (see, e.g., Ref. [110]) [111]. It will be interesting to investigate the relation of the negative fermion mass to a theory in the presence of a topological term in two dimensions.
The proposed experiment consists of preparing the ground states of Eqs. (3) for different values of the coupling g and subsequently measuring the magnetic field energy, which is proportional to [see Eq. (2)]. In the strongcoupling regime g −2 1, the magnetic field energy vanishes = 0, since the ground state approaches the bare vacuum |vvvv |0 . For weak couplings, g −2 1, the magnetic term dominates and the ground state is a superposition of all electric field basis elements. As such, converges to 1 when g −2 → ∞. However, the truncation described in Sec. II A bounds to a smaller value [28]. To explain the underlying physics in the intermediate region between the strong-and weak-coupling regimes, we first examine the case in which the kinetic term of the Hamiltonian can be treated as a perturbation (Fig. 8). Perturbative calculations are valid for 1 |4m| and we show results for = 5 and m = −50. We then proceed to the nonperturbative case ( Fig. 9) with parameters = 5 and m = 0.1. Our analytical results within the perturbative regime are determined using our iterative algorithm described in Appendix C and are unaffected by truncation effects in the parameter l. Therefore, they also serve as a test of validity for exact results.
For = 5 and m = −50, the kinetic term can be treated as a perturbation and we take the sum of the electric and mass terms in Eqs. (3a) and (3c) as the bare Hamiltonian. We can thus well describe the strong-coupling regime g −2 1, in which a 2D effect occurs that is characterized by a jump of the expectation value of the plaquette operator. A plot of as a function of g −2 is shown in Fig. 8(a), where the black dashed line is obtained using exact diagonalization of the Hamiltonian in Eqs. calculations. We begin by analyzing the system for g −2 1 (marked as region 1 in Fig. 8) and examine increasing values of g −2 as we move through regions 2 and 3.
In region 1, the ground state of the Hamiltonian is essentially the bare vacuum |vvvv |0 , which is characterized by null energy and = 0. The dominating electric term prevents the creation of an electric field, despite the mass term incentivizing particle-antiparticle pair creation (m < 0). The magnetic term is negligible in the strongcoupling regime. With increasing g −2 , the cost of creating electric fields is reduced, until it becomes favorable to create particle-antiparticle pairs. In particular, for g −2 = −(4m) −1 , the energy relief 2m of creating a pair compensates the cost g 2 /2 of creating an electric field on the link connecting the pair.
The kinetic termĤ kin is responsible for the energy anticrossing between regions 1 and 2, as it allows the creation of particle-antiparticle pairs and the corresponding electric fields. It couples the vacuum to the plaquette states, which contain the maximum number of particle-antiparticle pairs (we refer to this as a fully filled plaquette). However, Gauss's law allows two different gauge-field configurations for the fully filled plaquette, as shown in Fig. 2(e). These states are given by |f (0) ± = 1 √ 2 |epep (|0 ± |1 ) and while they are degenerate with respect to the bare Hamil-tonianĤ E +Ĥ m , Ĥ kin is minimized for |f − and Ĥ B for |f + . The existence of these two configurations in the presence of the perturbationĤ kin and the magnetic termĤ B creates competition between two quasidegenerate vacua in regions 2 and 3 of Fig. 8(b). The ground states in these regions are described by the corresponding corrected states in perturbation theory |f ± = n |f (n) ± (see Appendix C), as shown in Fig. 8(b). Here, we see that the kinetic term facilitates the creation of particle-antiparticle pairs and drives the ground state from the vacuum in region 1 toward |f − in region 2.
The remarkable 2D feature of the theory is the jump of between regions 2 and 3, which is shown by the black dotted line in Fig. 8(a). This jump corresponds to the sharp energy anticrossing in Fig. 8(b) and follows from the competition between the quasidegenerate vacua in the presence of the kinetic and magnetic terms. As the relative weights ofĤ kin andĤ B change with g −2 , an anticrossing occurs and the ground state changes from |f − to |f + . More specifically, there is a value g c such that, for any Despite the fact that the magnetic Hamiltonian is not included in our perturbative analysis, we can analytically calculate this value g c (see Appendix C). By requiring the magnetic and kinetic contributions from the energy to be equal, we find g c = 0.012 + o( 8 ), which is in excellent agreement with the results obtained from exact diagonalization [see Fig. 8(a)]. Following the jump, in region 3, the ground state is thus |f + , until we reach the weak-coupling regime and the magnetic term becomes dominant compared to all other contributions.
The above analysis shows how pair-creation processes can lead to dynamically generated magnetic fluxes that result in negative values of the magnetic field energy. For negative mass, the considered effect of competing vacua occurs over a wide range of parameters. The strength of the kinetic term broadens the dip and shifts the position of the jump g c .
In the following, we consider the regime in which the kinetic term cannot be treated as a perturbation. For = 5 and m = 0.1, exact diagonalization results are shown in Fig. 9(a). In the intermediate region where g 1, the kinetic term has significant weight, which leads to entanglement between matter and gauge degrees of freedom in the ground state. In particular, in the limit |m|, g 2 , g −2 and for l → ∞, it can be shown for N plaquettes [112] that the energy is minimized for a magnetic flux of π , which is generated by pair-creation processes and corresponds to = −1 [which is not reached in Fig. 9(a)  and |GS (b) is that of the magnetic term (see Appendix D). For g −2 = 1, however, the ground state approximates the ground state of the truncated kinetic term for l = 1 with 92% fidelity (which can be increased by incrementing ). To quantify the ground-state entanglement between matter and gauge degrees of freedom, we calculate the entanglement entropy S(ρ g ) = −Tr[ρ g log ρ g ], where ρ g is the density matrix of the gauge degrees of freedom that remain after eliminating those that are redundant. The entanglement entropy is plotted as the red dashed line in Fig. 9(a).

B. Running coupling
In this section, we consider 2D effects in QED on a lattice subject to PBCs. The method described in Sec. II C and in Ref. [28] is based on the Hamiltonian formalism of LGTs, which allows for simulations that are unaffected by the problem of autocorrelations inherent to MCMC methods. Thus, we can compute physical observables at arbitrary values of the lattice spacing, ultimately allowing for reaching a well-controlled continuum limit. As such, future quantum hardware could be used to efficiently calculate, for instance, the bound-state mass spectrum of the theory, properties of the inner structure of such bound states, or form factors that are important for experiments [4,109].
We are therefore bound to consider local observables that only need a small number of lattice points. We consider the plaquette operator as a simple example, as was done in the pioneering work by Creutz [58] at the beginning of MCMC simulations of LGTs. Despite its local nature, the operator can be related to a fundamental parameter of the theory, namely, the renormalized coupling g ren [113]. Importantly, the 1D Schwinger model is super-renormalizable, meaning that the coupling does not get renormalized, while, as we show, renormalization is necessary for our 2D model.
Renormalization appears in quantum field theories through quantum fluctuations, i.e., the spontaneous generation of particle-antiparticle pairs from the vacuum. This phenomenon leads to a charge shielding (or antishielding for non-Abelian gauge theories) [114], which in turn changes the strength of the charge depending on the distance (or the energy scale) at which it is probed. As such, the charge becomes scale dependent and in this work we choose the inverse lattice spacing 1/a as the scale. The scale dependence, which is a renormalization effect, is referred to as the running of the coupling and its knowledge is fundamental in understanding the interactions between elementary particles. In particular, the running coupling serves as an input to interpret results from collider experiments, such as the Large Hadron Collider. Hence, the ability to compute the running coupling from the plaquette operator can have a direct impact on such experiments and our knowledge of elementary-particle interactions.
As explained in Sec. II C, we study the case of 2D QED without matter. Despite the absence of matter, and in contrast to the previously studied 1D Schwinger model, renormalization is needed; hence the calculation of the running coupling. A definition of the renormalized coupling g ren can be given through the ground-state expectation value of the plaquette [113] g 2 ren = where g is the bare coupling from the Hamiltonian [see Eqs. (4) and (5)]. By looking at Eq. (16), it follows that to cover the scale dependence of the coupling, we need to evaluate over a broad range of the lattice spacing a. Equivalently (we set a = 1 above), this means that we need to perform simulations at many values of g −2 , covering the whole spectrum between the extremal regions of strong and weak couplings, g −2 1 and g −2 1, respectively. In particular, the interesting region where perturbation theory is no longer applicable and bound states can be computed on not too large lattices is where g −2 1. Covering such a broad range of scales is the major problem in standard lattice gauge simulations. In fact, while approaching the weak-coupling regime, autocorrelation effects prevent the calculations from reaching very small values of the lattice spacing, making continuum extrapolations and hence convergence to meaningful results difficult. We remark that in large-scale lattice simulations, much more sophisticated definitions of the renormalized coupling g ren are generally used [4,115]. These alternative definitions allow us, besides other things, to disentangle cut-off effects, inherent to the plaquette coupling, from the true physical running. However, the alternative forms of g ren described in [4,115] are meaningful for large lattices only and hence are beyond the capabilities of current quantum simulators. In this work, we allow for a proof-of-concept demonstration that paves the way for future improvements.
With our method (see Sec. II C and Ref. [28]), we can simulate the system both in the strong-and in the weakcoupling regimes, with very modest truncations. The most difficult region to simulate is characterized by g −2 1, where convergence is studied in Ref. [28]. We demonstrate this in Fig. 10, where the ground-state expectation value of the plaquette operator is plotted against g −2 . For PBCs, the plaquette operator is given by where γ = e (γ = b) refers to the electric (magnetic) representation [see Eqs. (4b) and (5b)]. Accordingly, the blue and red lines are determined with these two representations, while the gray dotted lines are obtained using perturbation theory and describe the exact results in the extremal regions well. The left graph is for a truncation l = 1 and the right for l = 2. In Sec. III, we provide a protocol for quantum simulation of the system with a truncation l = 1 of the gauge field using nine qubits. As shown in Fig. 10, the use of a truncation l = 2 around the region of the plot where g −2 1 brings us much closer to the converged result. Using the same protocol as for l = 1 (see Sec. IV C), the simulation for l = 2 requires 15 qubits.

VI. CONCLUSIONS AND OUTLOOK
In this work, we propose a protocol to observe 2D effects in LGTs on currently available quantum computers. By using the methods in Ref. [28], we provide a practical VQE-based framework to simulate two toy models using NISQ devices. Importantly, we include the numerics for observing 2D phenomena in a basic building block of 2D QED with present-day quantum resources.
The effective models studied here include both dynamical matter and a nonminimal gauge-field truncation, providing the novel opportunity to study several 2D effects in LGTs. More specifically, we show how to observe dynamical generation of magnetic fields as a result of particleantiparticle pair creation and pave the way for an important first step toward simulating short-distance quantities such as the running coupling of QED. While the protocols presented in Sec. III are designed for trapped-ion systems, our approach can easily be adapted to suit different types of quantum hardware.
One of the most appealing characteristics of our approach is that it can be generalized to more complex systems. Immediate extensions of our results include simulations of a 2D plane of multiple plaquettes and of QED in three spatial dimensions. Moreover, the inclusion of fermionic matter in the case of PBCs can be achieved by following the same procedure as for OBCs [28]. It will also be interesting to explore implementations on nonqubit-based hardware that is capable of representing gauge degrees of freedom with spin-l systems (l > 1/2) [98]. Another extension is the development of Trotter-type protocols to simulate real-time evolution using the encoding given in Sec. IV A. As part of the quest to move toward quantum simulations of QCD, another possible extension is to progress from a U(1) gauge theory (QED) to a non-Abelian gauge theory, such as SU(2) and eventually to SU (3). Importantly, simulations of LGTs with larger lattice sizes will become feasible with future advancements in quantum computing, allowing for exciting possibilities, such as the ability to relate the running coupling of a gauge theory to a physical parameter and to make connections between quantum simulations and experiments in high-energy physics. Ultimately, future quantum computers may also offer the potential to simulate models with a topological term, a nonzero chemical potential, or real-time phenomena, effects that are very hard or even impossible to access with MCMC techniques.
The field of quantum simulations of LGTs is in its early exploratory stages and is rapidly developing. The type of problems considered here will act as a test bed for quantum computer implementations. By providing practical and experimentally feasible solutions for simulating gauge theories beyond 1D, our work opens up new pathways toward accessing regimes that are classically out of reach.

APPENDIX A: EFFECTIVE HAMILTONIAN FOR A LADDER OF PLAQUETTES WITH OBCs
In this appendix, we derive an expression for the electric term in the Hamiltonian in Eq. (1) for a ladder of N plaquettes [see Fig. 2(a)]. Following the procedure outlined in Sec. II B, we implement Gauss's law to remove redundant gauge degrees of freedom and obtain an expression with only the independent gauge fields remaining. We then use the Jordan-Wigner transformation and the mapping in Eqs. (8) to arrive at the final Hamiltonian for the ladder. The magnetic, mass, and kinetic terms remain unchanged, except for the identity operators where the degrees of freedom are removed.
To write the electric termĤ E in Eqs. (1), we express the electric field on each link in terms of the chosen independent gauge and matter degrees of freedom (see Sec. II A). Categorizing links according to their orientation and location [defined in Fig. 2(a)], we examine the following cases: horizontal independent links, horizontal dependent links, outermost vertical links, and internal vertical links. Listed in Table I are the electric field contributions of each of these terms for even-and odd-numbered plaquettes, written in terms of the horizontal independent links.
Squaring each contribution in Table I and summing over the plaquettes allows us to arrive at an expression for the electric term in terms of only the independent gauge fields. Using the Jordan-Wigner transformation and Eqs. (8)-(9), the Hamiltonian for the ladder can then be written as follows:   Fig. 2(a). For the vertical outermost links, the cases of "even" and "odd" refer to the total number of plaquettes N .

APPENDIX B: OPTIMAL PARTITIONING FOR HAMILTONIAN AVERAGING
Our VQE algorithm requires the measurement of the energy of the variational state, i.e., the cost function C(θ). Precise estimation of the expectation value of a complicated Hamiltonian demands, in general, a large number of measurements, which translates to extensive computational run time. In this appendix, we show how the number of measurements required for the proposed quantum simulations can be substantially reduced by combining the standard approach [51,65] with the measurement strategy put forward in Ref. [116]. We adapt the latter to our models and provide circuits for their implementation on different types of quantum hardware. The described scheme assumes local measurements only.
By using Pauli operators as a basis, we can express any target HamiltonianĤ T aŝ where K is the number of terms in the decomposition, n is the dimension of the system,σ k i i is a Pauli operator acting on the ith qubit, and k i can either be 0 (identity), x, y, or z for any i. We remark that, for our qubit encoding (see Sec. IV A), decomposition of the Hamiltonian as in Eq. (B1a) can be done efficiently using the relations given in Eq. (11).
The cost function C(θ) is evaluated by measuring the expectation values Q k with single-qubit Pauli measurements and calculating the weighted sum Ĥ T = K k=1 c k Q k . OperatorsQ k that commute with each other can be measured simultaneously and will be grouped together in the following. Thus, we rewrite Eq. (B1a) aŝ Here, the set V m contains indices of commuting operatorŝ Q k . Importantly, the number M of commuting sets can be drastically smaller than the total number of Pauli operators K in Eq. (B1a), leading to a reduction in the experimental run time.
To successfully run a VQE algorithm, we must be able to calculate the average values of the energy for any variational state preparation. Due to the statistical nature of quantum mechanics, the exact average value of any observable is unknown. Thus, we need to approximate the energy with an estimator. For a generic operatorÔ, we denote the estimator of its expectation value Ô withŌ. Depending on the intrinsic variance ( Ô ) 2 and the number N of repeated measurements performed, the estimator O is affected by a statistical error, which is characterized by its variance Var [Ō].
To approximate the energy of the system, we choose beforehand a certain precision, which is generally problem and/or application dependent. In practice, this means that we are required to fix a threshold 2 Besides the chosen threshold 2 , the total number of measurements N is proportional to the intrinsic variances ( R m ) 2 and the number of partitioned sets M . In general, there are many different ways to partition theQ k into theR m and we seek to identify a strategy that minimizes both M and ( R m ) 2 . Importantly, ( R m ) 2 depends on the intrinsic variances and the covariances of theQ k from which eachR m is comprised. While a particular partitioning strategy might achieve M < K, it might also group togetherQ k with positive covariance, thereby adding to ( R m ) 2 and increasing the number of measurements needed to achieve 030334-19 2 precision ofH T . For instance, as pointed out in Ref. [51], the Hamiltonian could be partitioned in the following ways: When measuring the state | = |01 , partition 1 has variance 5 m=1 ( R m ) 2 = 2 and M = 5 commuting sets for a total of N (1) = 10/ 2 measurements [see Eq. (B5)]. Partition 2 has the same variance because the operators within the sets have zero covariance, but here M = 3, so N (2) = 6/ 2 . Although partition 3 has only M = 2 commuting sets, the covariance betweenσ x 1σ x 2 andσ y 1σ y 2 adds to ( R 1 ) 2 , giving N (3) = 8/ 2 measurements. In this example, partition 2 gives the fewest number of measurements even though partition 3 has the fewest number of commuting sets.
The variances and covariances ofQ k are dependent on the state | being measured, which is generally unknown. Thus, to reduce the total number of measurements required, we seek a partitioning strategy that yields the minimum number of commuting sets, i.e., we minimize M without regard for the variance of the partitions. In principle, it is possible to design an algorithm that, while acquiring data, reorganizes the partitions in order to minimize their variances, but this goes beyond the scope of this work [117].
The algorithm presented in Ref. [116] finds the partition of theQ k intoR m that gives the near-minimum number of commuting sets M [118]. Although commuting operators can, in principle, be measured simultaneously, from an experimental point of view, only those operators that bitwise commute can easily be measured simultaneously with the local measurements provided by the quantum device. To exemplify this, considerσ x 1σ x 2 andσ z 1σ z 2 . While these operators commute and thus can be simultaneously measured, there are no single-qubit operations that rotate them into a diagonal form. To achieve this, two-qubit gates are required. Reference [116] provides the procedure for building a quantum circuit that diagonalizes the Pauli oper-atorsQ k within a commuting set (for an example, see Fig. 11). This circuit is to be placed after the VQE circuit and before the local measurements, effectively rotating the measurement basis and allowing us to measure in a complicated entangled basis. As such, all commuting operators within a commuting set, bitwise or not, can be measured simultaneously. The gate set from which the quantum circuits are built can be adapted to suit the quantum computing platform, using, e.g., CNOT operations as entangling gates for superconducting platforms and iSWAP gates for ion-based platforms. While it is possible to express a CNOT gate in terms of iSWAP gates, and vice versa, this leads to linear overhead in circuit depth, which impacts the feasibility of experiments with NISQ devices. Therefore, we adapt the algorithm in Ref. [116] according to the preferred type of two-qubit gate to avoid this linear overhead in the number of gates. As an example, Figs. 11(a) and 11(b) show two possible circuits that diagonalize one of the partitionsR m of the plaquette with OBCs introduced in Sec. II B, with the qubit labeling convention shown in Fig. 11(c). The ele-mentsQ k of this partition are shown in Fig. 11(d), along with the diagonal forms obtained after applying the CNOT and iSWAP circuits (for details, see Refs. [116,117]).
In the worst-case scenario, the depth of the additional circuits scales quadratically with K; however, in practice, we often see linear scaling. The computing time of the algorithm in Ref. [116] also scales quadratically in K, while the number of qubits in the system does not play a significant role in the complexity [119]. By applying this method to the Hamiltonian of Eqs. (13), we reduce the number of measurement sets (i.e., commuting sets) from 49 down to M = 5 in the case of OBCs. With PBCs, the number of measurement sets is reduced from 295 to M = 11 for the electric representation [Hamiltonian in Eqs. (14)] and from 138 down to M = 15 for the magnetic basis [Hamiltonian in Eqs. (15)]. The depth of the circuits used to implement the measurements is around ten in the worst-case scenario.

APPENDIX C: PERTURBATION THEORY
In this appendix, we outline the perturbative algorithm used to determine the analytical results presented in Sec. V A. We consider the system Hamiltonian to beĤ =Ĥ 0 + H kin , whereĤ 0 =Ĥ E +Ĥ m [see Eqs (3a) and (3c)] is the bare Hamiltonian andĤ kin [see Eq. (3d)] is the perturbation. Following standard perturbation theory, the corrected state and energy for eigenvalue i up to order n can be written as The unperturbed state | (0) i can be any of the eigenstates ofĤ 0 . Our goal is to find the corrections for j ≥ 1 in Eqs. (C1) up to an arbitrary order and for any initial state. We require higher orders in perturbation theory because the leading and the second leading contributions are found at fourth and sixth orders, respectively. Furthermore, degeneracies in our system are often lifted at fourth order. Importantly, here we do not truncate the gauge-field operators. Since the action ofĤ kin is of creating and/or annihilating, whenever possible, a fermion-antifermion pair and fixing the gauge field between, we can select all contributing states at each perturbation order.
The eigenstates and corresponding eigenvalues of the bare HamiltonianĤ 0 are denoted as | (0) i,μ and E (0) i , where μ is an index used within the degenerate subspace that contains the unperturbed state | (0) i . As an example, the six lowest-energy eigenstates ofĤ 0 are as follows (see Fig. 12): the nondegenerate vacuum state |vvvv |0 with 0th order 1st order 2nd order 3rd order FIG. 12. A schematic description of our perturbative algorithm when starting from the vacuum. In the first column, we show the bare energies of the corresponding states. Each application of the perturbationĤ kin allows for creating and annihilating particle-antiparticle pairs, as depicted by the tree structure. At each iteration, our algorithm expands the Hilbert space to include new terms contributing to the correction. This allows for efficient use of the resources, ultimately reaching several hundred orders in perturbation theory. null energy; the fourfold degenerate two-particle states |epvv |0 , |vpev |1 , |vvep |0 , and |evvp |0 , each with energy g 2 /2 + 2m; and the twofold degenerate fully filled plaquette states |epep |0 and |epep |1 , each with energy g 2 + 4m [see Fig. 2(e)]. Acting on an eigenstate ofĤ 0 with the kinetic term creates and annihilates particle-antiparticle pairs and fixes the field on the links to satisfy Gauss's law. Successive applications ofĤ kin on the vacuum are shown in Fig. 12. This illustration shows that actingĤ kin on a state an odd number of times results in a state that is orthogonal to the initial state, implying that odd-order perturbative energy corrections are vanishing [see Eq. (C2b)].
At the beginning of our algorithm, the initial state | (0) i is set to an unknown superposition of states within the selected eigensubspace i, i.e., | (0) If the subspace i is degenerate, the coefficients d (0) ν are determined by lifting-when possible-the degeneracy at higher orders (see below). The nth-order perturbative corrections are calculated as where | (0) k with energy E (0) k are the eigenstates ofĤ 0 outside the degenerate eigensubspace that contains the initial state. The overlap d (n) ν of the state correction with the eigenstates | (0) i,ν within the same degenerate subspace can be determined by imposing normalization when the degeneracy is lifted. The specific form of Eqs. (C2) follows from the symmetries that our perturbation H kin satisfies. In particular, E (n) i = d (n) ν = 0 for odd n. To lift a degeneracy, the coefficients d (0) ν need to be found by diagonalizing the matrix equation (C3) Equation (C3) can only be solved at the order m at which the perturbation lifts the degeneracy. For instance, when starting from the fully filled plaquette states |epep |0 and |epep |1 , the twofold degeneracy is lifted at m = 4. However, the fourfold degeneracy between the two-particle states is only partially lifted at m = 2 into two twofold degenerate systems {|evvp |0 , |vpev |1 } and {|epvv |0 , |vvep |0 }. The former is fully lifted at m = 4. The latter remains degenerate for all orders because these states satisfy the same symmetry as the perturbation H kin . As mentioned in Sec. V A, the jump coordinate g c can be estimated using perturbation theory. The jump occurs when the ground-state transitions from |f − to |f + (see Sec. V A), i.e., when the energies of the two states are equal. Given that the bare energies E (0) i are the same, this means that for g = g c the sum of the kinetic E kin and magnetic E B = − /(2g 2 ) energies for the two states are equal. While E kin is the energy correction given from perturbation theory, one can determine E B using the corrected state | i = |f ± . Thus, we can write E (+) where the superscript (±) indicates the perturbative state |f ± used. By solving for g, we arrive at the estimated jump coordinate up to order n in perturbation theory.

APPENDIX D: GROUND STATE OF OPEN-BOUNDARY PLAQUETTE IN THE WEAK-COUPLING LIMIT
In this appendix, we derive the ground state of the plaquette with OBCs in the weak-coupling regime of g −2 1. The following derivations are valid for the nonperturbative regime described in Sec. V A [see Fig. 9(d)], where |m| and g −2 . The Hamiltonian is given in Eqs. (1).
For large g −2 , the Hamiltonian is dominated by the magnetic termĤ B , which contains only gauge-field operators. As such, the ground state is required to be a product state between some matter configuration and ground state of the magnetic term. This is confirmed by the plot of the entanglement entropy in Fig. 9(a), which shows that the ground state is separable between matter and gauge components for large g −2 . Therefore, we set the ground state to be a tensor product | = |ψ |GS (b) between some matter configuration |ψ = 6 j =1 c j |m j and the ground state of the magnetic term |GS (b) . Here, |m j are the allowed matter configurations {|vvvv , |epvv , |vpev , |vvep , |evvp , |epep } for j = 1, 2, . . . , 6, respectively (see Table II). While the magnetic term is dominant, it does not influence the matter configuration. For simplicity, we set the mass m = 0 (the following procedure can be generalized for nonzero m), leaving the remaining candidates for determination of the matter state to beĤ E andĤ kin . We are in the limit 1, g −2 1 and so the kinetic term is dominant over the electric term. By minimizing |Ĥ kin | over the c j , we can estimate the matter component of the ground state at large g −2 .
The expression for |Ĥ kin | can be written as and the optimal c j are shown in Table II for truncation l = 1 and l → ∞. The residual probabilities of the states |vpev and |evvp in Fig. 9(d) are an effect of the l = 1 truncation. Constructing the estimated ground state | = |ψ |GS (b) using the optimal c j , we find that | has over 99.9% overlap with the ground state obtained via exact diagonalization of the Hamiltonian in Eqs. (1) for large g −2 .