A resource efficient approach for quantum and classical simulations of gauge theories in particle physics

Gauge theories establish the standard model of particle physics, and lattice gauge theory (LGT) calculations employing Markov Chain Monte Carlo (MCMC) methods have been pivotal in our understanding of fundamental interactions. The present limitations of MCMC techniques may be overcome by Hamiltonian-based simulations on classical or quantum devices, which further provide the potential to address questions that lay beyond the capabilities of the current approaches. However, for continuous gauge groups, Hamiltonian-based formulations involve infinite-dimensional gauge degrees of freedom that can solely be handled by truncation. Current truncation schemes require dramatically increasing computational resources at small values of the bare couplings, where magnetic field effects become important. Such limitation precludes one from `taking the continuous limit' while working with finite resources. To overcome this limitation, we provide a resource-efficient protocol to simulate LGTs with continuous gauge groups in the Hamiltonian formulation. Our new method allows for calculations at arbitrary values of the bare coupling and lattice spacing. The approach consists of the combination of a Hilbert space truncation with a regularization of the gauge group, which permits an efficient description of the magnetically-dominated regime. We use $2+1$ dimensional quantum electrodynamics as a benchmark example to demonstrate this efficient framework to achieve the continuum limit in LGTs. This possibility is a key requirement to make quantitative predictions at the field theory level and offers the long-term perspective to utilise quantum simulations to compute physically meaningful quantities in regimes that are precluded to quantum Monte Carlo.

Gauge theories establish the standard model of particle physics, and lattice gauge theory (LGT) calculations employing Markov Chain Monte Carlo (MCMC) methods have been pivotal in our understanding of fundamental interactions. The present limitations of MCMC techniques may be overcome by Hamiltonian-based simulations on classical or quantum devices, which further provide the potential to address questions that lay beyond the capabilities of the current approaches. However, for continuous gauge groups, Hamiltonian-based formulations involve infinite-dimensional gauge degrees of freedom that can solely be handled by truncation. Current truncation schemes require dramatically increasing computational resources at small values of the bare couplings, where magnetic field effects become important. Such limitation precludes one from 'taking the continuous limit' while working with finite resources.
To overcome this limitation, we provide a resource-efficient protocol to simulate LGTs with continuous gauge groups in the Hamiltonian formulation. Our new method allows for calculations at arbitrary values of the bare coupling and lattice spacing. The approach consists of the combination of a Hilbert space truncation with a regularization of the gauge group, which permits an efficient description of the magneticallydominated regime. We use 2 + 1 dimensional quantum electrodynamics as a benchmark example to demonstrate this efficient framework to achieve the continuum limit in LGTs. This possibility is a key requirement to make quan-Jan F. Haase: jan.frhaase@gmail.com, contributed equally Luca Dellantonio: luca.dellantonio@uwaterloo.ca, contributed equally Alessio Celi: alessio.celi@uab.cat, contributed equally titative predictions at the field theory level and offers the long-term perspective to utilise quantum simulations to compute physically meaningful quantities in regimes that are precluded to quantum Monte Carlo.

Introduction
Gauge theories are the basis of high energy physics and the foundation of the standard model (SM). They describe the elementary interactions between particles, which are mediated by the electroweak and strong forces [1][2][3], making the SM one of the most successful theories with tremendous predictive power [4]. Still, there are numerous phenomena which cannot be explained by the SM. Examples include the nature of dark matter, the hierarchy of forces and quark masses, the matter antimatter asymmetry and the amount of CP violation [5]. Answering these questions and accessing physics beyond the SM, though, often requires the study of non-perturbative effects. A very successful approach to address nonpertubative phenomena is lattice gauge theory (LGT) [6][7][8], as proposed by Kenneth Wilson in 1974 [9]. In LGT, Feynman's path integral formulation of quantum field theories (QFTs) is employed on an Euclidean space-time grid. Such a discretised form of the path integral allows for numerical simulations utilizing Markov Chain Monte Carlo (MCMC) methods. The prime target of LGT is quantum chromodynamics (QCD), i.e. the theory of strong interactions between quarks and gluons. In this field, LGT has been extremely successful, allowing for example the computation of the the low-lying baryon spectrum [10], the structure of hadrons, fundamental parameters of the theory and many more [11][12][13][14].
However, many of the aforementioned open questions in modern physics cannot be addressed within the standard approach, due to the sign-problem [15][16][17] that renders MCMC methods ineffective. A possible solution is to employ a Hamiltonian formulation of the underlying model. Classical Hamiltonian-based simulations using tensor network states (TNS) have been successful [18][19][20][21][22], but are so far restricted to mostly one spatial dimension. Consequently, there is a necessity for new approaches to both access higher dimensions and address problems where standard MCMC methods fail. It is presently not known whether efficient classical methods can be developed to overcome this problem.
Hamiltonian-based simulations on quantum hardware provide an alternative route, since there is no such fundamental obstacle to simulating QFTs in higher dimensions [24][25][26]. Therefore, this approach holds the potential to address questions that cannot be answered with current and even future classical computers. The rapidly evolving experimental capabilities of quantum technologies [27,28] have led to proof-of-concept demonstrations of simulators tackling one-dimensional theories [29][30][31][32][33][34][35]. Yet, extending these results to higher dimensions represents a crucial step for this field, and realisations on 'Noisy Intermediate-Scale Quantum' devices [36,37], i.e. current quantum hardware, requires novel approaches to make this leap.
To meet this challenge, we provide a resourceefficient approach that facilitates the quantum simulation of higher dimensional LGTs that would otherwise be out of reach for current and near-term quantum hardware, which is exemplified Table 1. In addition, purely classical simulations based on the Hamiltonian formalism also benefit from our resource-optimised approach. Hence, we bring both quantum and classical calculations closer to developing computational strategies that do not rely on Monte Carlo methods, and thus circumvent their fundamental limitations.
Our new approach addresses the important problem of reaching the continuum limit (in which the lattice spacing approaches zero) with finite computational resources. Since QFTs are continuous in their time and space variables, the need to take a controlled continuum limit is inherent to any lattice approach and necessary to extract physically relevant results from a lattice simulation.
Taking QCD as a concrete example, we require an accurate description for particles interacting at both short and long distances. Lattice QCD and other LGTs offer the unique tool to investigate both regimes. At long distances, e.g. the bound state spectrum can be computed. At short distances, and after taking the continuum limit, it is possible to connect the perturbative results derived with QFTs with nonperturbative simulations, thus assessing the range in which perturbation theory is valid. However, taking the continuum limit is in general computationally expensive. MCMC methods, for instance, have the intrinsic problem of autocorrelations, that become more and more severe when decreasing the lattice spacing. This drawback in turn leads to a significant increase in the computational cost, and fixes the smallest value of the lattice spacing that can be reached. On one hand, Hamiltonian approaches circumvent this problem. On the other, however, Hamiltonian-based formulations face the challenge that for continuous (Abelian and non-Abelian) gauge groups, local gauge degrees of freedom are defined in an infinite dimensional Hilbert Table 1: Computational cost for different approaches. We estimate the number of states required to reach a 1% accuracy in the expectation value of the two-dimensional plaquette in QED (see Sec. 4.3) when compared to the value we obtain with our method considering a maximum of 9261 basis elements. The three columns refer, from left to right, to the standard approach described in Sec. 2.1, our approach (see Sec. 3) using a fixed group ZN , and finally our optimised strategy, in which the order of the group N is scaled with the bare coupling g (see Sec. 4). The shown savings in computational resources bring quantum simulations with current technology within reach. Note that 125 states correspond to seven qubits. We present a robust implementation strategy for ion-based quantum computers in [23].
space. As a consequence, any simulation -classical or quantum -requires a truncation of the gauge fields, which is inherently at conflict with the required continuum limit.
In this work, we present a practical solution to overcome this crucial bottleneck and to allow for resourceefficient Hamiltonian simulations of LGTs. Although our approach is general and applicable to LGTs of any dimension, we consider two-dimensional quantum electrodynamics (QED) as a benchmark example.
Truncation of the gauge fields is typically performed in the 'electric basis', i.e. the basis in which the electric Hamiltonian and Gauss' law are diagonal. As such, truncation preserves the gauge symmetry, and the resulting models are known as gauge magnets or link models [38][39][40], which are of direct relevance in condensed matter physics [41][42][43][44]. As recently shown in Ref. [45], spin-1/2 truncations are within the reach of current quantum simulators. From the perspective of fundamental particle interactions, electric truncations can result in an accurate description of the system in the strong coupling regime. However, by decreasing the value of the coupling or equivalently the lattice spacing, the magnetic contributions to the energy become increasingly important and the number of states that have to be included in the electric basis grows dramatically (a similar increase can be realised by adding an auxiliary spatial dimension to the lattice [46]). An alternative approach to describe the gauge degrees of freedom is to approximate continuous gauge groups with discrete ones, for instance, to approximate U (1) with Z 2L+1 (L ∈ N) [47,48]. Such approaches also face similar limitations as the ones described above, as L has to be progressively increased with decreasing coupling.
A natural solution to simulate the weak coupling regime consists of exploiting the self-duality [49] of the electric and magnetic terms by Fourier transforming the Hamiltonian and working in the 'magnetic basis', i.e. the basis in which magnetic interactions are diagonal, as suggested in [50]. However, the fact that the magnetic degrees of freedom are continuous variables with a gapless spectrum, poses intricate challenges for a resource-efficient truncation scheme, that have yet (to the best of our knowledge) to be addressed. In this work, we provide a practical solution by combining state truncation with a gauge group discretisation that is dynamically adjusted to the value of the coupling. This approach allows for controlled simulations at all values of the bare coupling, smoothly connecting the weak, strong and intermediate coupling regimes. As a proof-of-principle of this new approach and its ability to faithfully simulate non-perturbative phenomena, we target the renormalised coupling in QED in 2 + 1 dimensions.
To observe non-perturbative phenomena such as confinement, the simulated physical length scale needs to be larger than the scale at which confinement sets in. As a result, large lattice sizes are required and the number of lattice points grows rapidly when approaching the continuum limit of the theory. This results in computational requirements that cannot be satisfied using current classical and quantum computers. Still, as previously done in the pioneering work by Creutz [51], we can study the bare coupling dependence of the local plaquette operator. This quantity allows us to benchmark our formalism and to show that a smooth connection between the weak and the strong coupling regimes can be established. In addition, our method allows for estimating the precision with which a given truncation approximates the untruncated results.
The paper is organised as follows. In Sec. 2, we review lattice QED in 2 + 1 dimensions as an example of Abelian and non-Abelian gauge theories with continuous groups and magnetic interactions. We consider lattices with periodic boundary conditions and reformulate the lattice Hamiltonian in terms of gaugeinvariant degrees of freedom. By eliminating redundant variables, we obtain an effective Hamiltonian description that allows for simulations at a low computational cost. In Sec. 3, we introduce a new magnetic representation of lattice QED that is equipped with a regularisation in terms of a Z 2L+1 group and an effi-cient truncation scheme. In Sec. 4, we study the performance of our method and benchmark its precision by calculating the expectation value of the plaquette operator on a periodic plaquette in the static charge limit. We show that both the truncation cut-off parameter, i.e. the maximum number of gauge basis states included in the simulation, and L, the dimension of the Z 2L+1 group, can be used as adjustable variational parameters. Both are used to optimise the simulation and estimate its accuracy. In the following Sec. 5, we present the generalisation to an arbitrary, two-dimensional periodic lattice with dynamical matter. Finally, we outline the prospects of this method for classical and quantum simulations is in Sec. 6 2 Minimal encoding of LGTs with continuous gauge groups In this chapter, we provide a Hamiltonian formulation for LGTs with continuous gauge groups that allows for resource-efficient classical and quantum simulations. First, we review the standard Kogut-Susskind Hamiltonian subject to Gauss' law (the local constraints ensuring gauge invariance) in Sec. 2.1, considering QED on a square lattice as a paradigmatic example. In Sec. 2.2, we proceed to provide a minimal formulation of the QED lattice Hamiltonian, in which redundant degrees of freedom have been removed.

QED in two dimensions
We review here the bottom-up construction of the lattice Hamiltonian as originally presented in [52]. For the sake of simplicity, we consider QED in 2 + 1 dimensions which displays key features of phenomenologically relevant theories like QCD, including chiral symmetry breaking and a renormalisation of the coupling constant [53], features that are absent in one spatial dimension.
The Hamiltonian of Abelian and non-Abelian gauge theories in two (or more) dimensions is constructed in terms of electric and magnetic fields, and their coupling to charges. In continuous Abelian U (1) gauge theories like QED (and similarly for non-Abelian gauge theories like QCD), electric and magnetic fields are defined through the vector potential A µ , with Here t, x, y are the time and space coordinate in two dimensions, and µ = x, y.
Gauge invariance, i.e. invariance of the Hamiltonian under local phase (symmetry) transformations of the charges, follows directly from the invariance of E µ and B under A µ → A µ +∂ µ θ(x, y), where θ(x, y) is an arbitrary scalar function. The electric field is sourced by the charges through Gauss' law, µ ∂ µ E µ = 4πρ, where ρ is the charge density. In LGTs [9], the charges occupy the sites n = (n x , n y ) of the lattice while the electromagnetic fields are defined on the links. The links are denoted by their starting site n and their direction e µ (µ = x, y), as shown in Fig. 1. The electric interactions are defined in terms of the electric field operatorÊ n,eµ , which is Hermitian, possesses a discrete spectrum and acts on the link connecting the sites with coordinates n and n + e µ . For each link, one further defines a Wilson operatorÛ n,eµ , as the lowering operator for the electric field, [Ê n,eµ ,Û n ,eν ] = −δ n,n δ µ,νÛn,eµ . The Wilson operator measures the phase proportional to the bare coupling g acquired by a unit charge moved along the link (n, e µ ) of length a, i.e.Û n,eµ ∼ exp{iagÂ µ (n)}. The magnetic interactions are given by (oriented) products of Wilson operators on the links around the plaquettes of the lattice. These operators are used to construct the Kogut-Susskind Hamiltonian asĤ =Ĥ gauge +Ĥ matter . Let us discuss first the pure gauge part that describes the limit of static chargeŝ Here, the sums run over both components of the sites n = (n x , n y ) and P n =Û n,exÛn+ex,eyÛ † n+ey,exÛ † n,ey (2) is the plaquette operator. It is easy to check that Eq.
that determines what states are physical for a given distribution of charges. Here,q n is the operator measuring the charge on the site n and |Φ represents the state of the whole lattice, including both links and sites. The eigenstates of the electric field operatorŝ E n,eµ |E n,eµ = E n,eµ |E n,eµ , E n,eµ ∈ Z (4) form a basis for the link degrees of freedom. In particular, the physical states can be easily identified in this basis via Eq. (3).
Let us now consider moving charges. To ensure gauge invariance, their motion is required to respect Gauss' law, i.e. a charge q moving between two sites has to change the electric field along the path by −q. In other words, the lowering operatorÛ has to be applied q times to the links on the path to preserve gauge-invariance. SinceÛ q = exp{iqagÂ}, the socalled minimal coupling condition [7] is recovered in the continuum limit a → 0, which is equivalent to replacing derivatives of matter fields by the covariant derivatives, i.e. shifting the particles' momentum by a gauge field contributionp µ →p µ − qgÂ µ .
In QED, charges are represented by Dirac fermions. In the staggered representation [52], they are represented on a square lattice as ordinary fermions at half filling, with staggered chemical potential that plays the role of the mass term. Their Hamiltonian iŝ H matter =Ĥ M +Ĥ K , whereĤ M andĤ K are the mass and kinetic contributions, respectivelŷ Here, m and q are the particles' effective mass and (integer) charge, κ the kinetic strength andΨ ( †) n the fermionic lowering (raising) operator for site n. Sincê H M identifies the Dirac vacuum with the state with all odd sites occupied, creating (destroying) a particle at even (odd) site is equivalent to creating a (−)qcharged "fermion" ("antifermion") in the Dirac vacuum. Thus the tunneling processes in the kinetic term correspond to the creation or annihilation of particleantiparticle pairs and the corresponding change in the electric field string connecting them. The charge operatorq n is given bŷ where q is an integer number which we set to one in the following. Note that we rescaled the fermion field by a factor √ α as discussed in App. A, which establishes the relations with M being the bare mass of the particles. We conclude this section with a few comments on the structure of the pure gauge part of the Kogut-Susskind Hamiltonian in Eq. (1). There, the electric and magnetic terms show an apparent asymmetry that obscures the electromagnetic duality in QED in the continuum and in the Wilson lattice formulation [9].
The symmetry between electric and magnetic fields in QED and in Wilson's action-formulation is due to the fact that time and space are treated on the same footing. Wilson's lattice action theory is formulated on a space-time grid with lattice spacing a µ , µ = t, x, y, z. In this case, an isotropic continuum limit is taken in which the lattice spacings in both the temporal and the spatial directions approach zero. In the Hamiltonian formulation, time is a continuous parameter. Accordingly, the above procedure is broken into two steps. Firstly, the continuum limit with respect to time a t → 0 is taken, which results in the Hamiltonian lattice formulation used here. In a second step, the continuum limit has to be taken with respect to space a x,y,z → 0 to obtain physical results.
In the Hamiltonian formulation, the electric field operatorsÊ form an algebra and are non-compact, as their integer spectrum takes values from minus infinity to infinity. By contrast, the Wilson operatorsÛ and hence the plaquette operatorsP , form a group. The moduli of their expectation value is one, as is the case for the Wilson action. More specifically, in Wilson's action formulation, it follows from operatorŝ U = exp{iga µÂµ } thatÂ µ is compact as it is defined from −π/(ga µ ) to π/(ga µ ), where a µ is the lattice spacing in the µ-direction, which includes both space and time as µ = t, x, y, z. The Kogut-Susskind Hamiltonian can be obtained from the Wilson action by taking the continuum limit in the time direction a t → 0 [51]. Thus, the asymmetry between the electric and magnetic terms in the Hamiltonian formulation disappears when the continuum limit is taken in the spatial direction.
While a fully non-compact formulation of Hamiltonian LGT is possible [54] (for the different outcomes of the two approaches see e.g. [55,56]), we do not discuss this approach here as it is not advantageous for quantum simulations. As we show in Sec. 3, it is instead convenient to write the electric term in a compact form.

QED Hamiltonian for physical states
As outlined in the previous section, gauge-invariance constrains the dynamics to the physical states only, i.e. those satisfying Gauss' law in Eq. (3). Practically, unphysical states have to be suppressed, e.g. via energy penalties [57]. In any case, quantum states that are not physical represent an exponential overhead for classical and quantum computation (also after a proper truncation, see Sec. 4). Furthermore, in noisy near term quantum devices or simulation protocols where the Hamiltonian has to be split up, e.g. to simulate time-evolutions employing a product formulas such as the Trotter expansion [58], implementing or imposing Gauss' law during the simulation may be complicated, or even impossible.
It is thus convenient to eliminate the redundant degrees of freedom by solving the constraint at each lattice site. In one dimension, such a procedure allows one to completely eliminate the gauge field, yielding an effective Hamiltonian containing only matter terms (but long-range interactions) [59,60]. A similar approach is applicable in higher dimensions, with the difference that the gauge field has also physical degrees of freedom. Here, we show how to formulate an effective Hamiltonian that directly incorporates the constraints of Eq. (3) by employing a convenient Figure 1: Two-dimensional lattice gauge theory with periodic boundary conditions. A single cell of the periodic 2D lattice in (a) is made of four links, oriented towards the positive x and y directions. Each lattice site is indicated by a unique vector n, which marks the lower left corner of each single plaquette. The associated operatorPn accounts for the electric field quanta circulating along the sketched path. The periodic lattice spans the surface of a torus, shown in the middle, whose minimal instance is assembled by four sites and the corresponding electric fields [thick lines, same color coding as in (a)]. Unwrapping this minimal torus yields the geometry shown in (b). We identify the stringsRx andRy and the four rotatorsRj, j = 1, 2, 3, 4. The eigenstates of the strings and three of the rotators (we arbitrarily removeR4, dashed loop) form a basis for the physical states of the pure gauge theory. To describe the physical states for a generic charge configuration we add three charge strings (dotted green arrows) that correspond to a conventional physical state for the given charge configuration.
parametrization of the physical states that yields an intuitive description of the system.
For the sake of clarity, we consider the minimal instance of a two-dimensional gauge theory: a single plaquette with periodic boundary conditions. The generalisation to an arbitrary lattice on a torus is derived in Sec. 5. Due to the periodic boundary conditions, this minimal system can equivalently be represented as a torus with four faces, or as four distinct plaquettes consisting of eight links [see Fig. 1 Describing the system in terms of these five links, however, entails serious drawbacks. The effective Hamiltonian contains many body interactions which are challenging or even impossible to be implemented using available quantum hardware (see App. B.1). To circumvent this problem, we consider a natural basis for the physical states in terms of small loops around each plaquette, and large electric loops around the whole lattice. In such a basis, the electric and magnetic interactions take a simple form. To conveniently describe these interactions, we introduce a set of operators, rotators and strings (see Fig. 1), that are diagonal and label the loop basis. As we show in [23], the Hamiltonian formulated in terms of these operators can be simulated with current quantum hardware.
With the notation and conventions presented in where the chargesq n are required by Gauss' law. An intuitive way to understand the effect of the charge terms in Eq. (9) is to consider them as sources of additional electric strings (whose concrete choice is just a matter of convention but consistent with Gauss' law), as displayed by the green lines in Fig. 1(b). We remark that this becomes evident from the link formulation in App. B.1. As mentioned above, rotators and strings automatically satisfy Gauss' law, which can be readily verified by observing that at any site, the incoming fields are always balanced by the outgoing ones. Moreover, by recalling the plaquette operatorP n in Eq. (2), it becomes clear whyR i andR µ are a convenient choice to represent the electric gauge field components. The operatorP n increases the anticlockwise quanta of the electric field circulating in the n-th plaquette. Consequently, it does not act on strings and takes the form of the lowering operator of the associated rotator. This fact can be formally proven by examining the raising and lowering operators of rotators and strings. From the commutation relations of the links and the relations shown in Eq. (9), it follows that R y ,Û (0,0),eyÛ(0,1),ey =Û (0,0),eyÛ(0,1),ey ≡P y , (10) whereP j , j = 1, 2, 3 is the plaquette operator of plaquette j as denoted in Fig. 1. Moreover, we defined the string lowering operatorsP x ≡Û (0,0),exÛ(1,0),ex andP y ≡Û (0,0),eyÛ(0,1),ey .
The magnetic Hamiltonian for the periodic plaquette in Fig. 1(b), is proportional to the sum of four plaquette operators, while there are only three independent rotators. The fourth rotator can be written as a combination of the others, since the effect of lowering (raising) all other rotators, i.e.R 1 ,R 2 andR 3 , amounts to raising (lowering)R 4 . This relation can be understood by examining Eq. (9): By lowering all of the three rotatorsR 1 ,R 2 andR 3 , we manipulate the electric fields on the links constitutingR 4 in exactly the same way as an increment of the latter would do. As such, the magnetic Hamiltonian becomeŝ while, by inserting Eq. (9) into Eq. (1), the electric term takes the form: Once the effective gauge HamiltonianĤ gauge = H E +Ĥ B has been derived in terms of rotator and string operators, we must further modify the matter HamiltonianĤ matter =Ĥ M +Ĥ K [for the description in terms of field operators, see App. B.1]. While the mass term in Eq. (5) is independent of the gauge fields, the kinetic contribution has to be rephrased. The kinetic contribution in Eq. (6) corresponds to the creation or annihilation of a particle-antiparticle pair on neighbouring lattice sites and the simultaneous adjustment of the electric field on the link in between. The green lines in Fig. 1(b) mark the fieldŝ E (0,0),ex ,Ê (0,0),ey andÊ (1,0),ey which are automatically adjusted when charges are created. This fact follows from our arbitrary choice of enforcing the three Gauss' law constraints on exactly those links. For any other link, we require combinations of raising and lowering operatorsP j andP µ (j = 1, 2, 3 and µ = x, y) such that the specific link is adjusted, while all others remain unchanged. As an example, let us consider the generation of a particle in position (1, 1) and an antiparticle in (1, 0). This choice implies either that the electric fieldÊ (1,0),ey has to decrease [which is automatically adjusted through the creation of a charge string], or that the electric fieldÊ (1,1),ey has to increase and hence the rotatorR 3 has to decrease. However, this action changes the electric fieldsÊ (1,1),ex , E (0,1),ey andÊ (1,0),ex as well. To remedy that, we lower the rotatorR 2 , adjustingÊ (1,1),ex andÊ (1,0),ex , and raise the stringR y to compensate for the change inÊ (0,1),ey . Following the same procedure, the rules for translating the kinetic Hamiltonian of Eq. (6) into the language of rotators and strings read Inserting these into Eq. (6), we obtain the kinetic contribution to the total Hamiltonian aŝ In conclusion, with the gauge partĤ gauge of the Hamiltonian described by Eqs. (12) and (13) and the matter partĤ matter by Eqs. (5) and (15), the single plaquette is fully characterised.
The effective Hamiltonian we derive here for a periodic plaquette can be extended to a torus of arbitrary size [see Sec. 5.2]. We will use this Hamiltonian to compute the expectation value of the plaquette operator where |Ψ 0 is the ground state, and V the number of plaquettes in the lattice, V = 4 in this case. The expectation value of the operator is defined as a dimensionless number, which is bounded by ±1 and proportional to the magnetic energy.

Transformation into the magnetic representation
In the following, we describe a scheme that allows switching from the so-called electric representation, whereĤ E is diagonal, to the magnetic one, whereĤ B is diagonal. Our method is based on the replacement of the U (1) gauge group with the group Z 2L+1 , and an accompanying transition from the compact formulation to a completely compact formulation, where both field degrees of freedom are treated as compact variables. While this procedure is general, we illustrate it for a single, pure-gauge plaquette and consider generalisations in Sec. 5.
Before presenting the scheme, we discuss the following observations about the considered Hamiltonian, that is now reduced to the sum of Eqs. (12) and (13), while all chargesq n in Eq. (13) are set to zero. In particular, the lowering (raising) oper-atorsP acting on the strings are solely contained in the now absent kinetic Hamiltonian in Eq. (15). The total Hamiltonian thus commutes witĥ and as a consequence the strings become constants of motion. The dynamics induced by the pure-gauge Hamiltonian are thus restricted to different subspaces defined byR µ |r µ = r µ |r µ , for µ = x, y. Starting in Sec. 4, we will be interested in a ground state property, therefore we restrict ourselves to the subspace where both strings are confined to the vacuum. The effective Hamiltonian of this subspace can be readily obtained by settingR x =R y = 0 in Eqs. (12) and (13) which yieldŝ where we introduced the superscript (e) to emphasise is the electric representation.
Since the three rotators possess discrete but infinite spectra, any numerical approach for simulating the Hamiltonian in Eq. (17) requires a truncation of the Hilbert space. In the following, l denotes the cut-off value which is identified bŷ R j |r j = r j |r j ∀ r j = −l, −l + 1, . . . , l. (18) Thus, the action of the truncated lowering operators is given asP Note that the total dimension of the Hilbert space is reduced to (2l + 1) 3 , which is still challenging to simulate even for relatively small values of l. In particular, calculations in the weak coupling regime suffer from this severe limitation and until now, no practical methods to solve this issue have been available.
Let us now introduce a formulation that allows for an efficient representation of the Hamiltonian's eigenstates in the weak coupling regime, where g 1. It is based on the exchange of the continuous U (1) group with the discrete group Z 2L+1 , which provides a discrete basis for the vector potential operatorsÂ n,eµ and enables a direct transformation into this dual basis via a Fourier transform. The approach is motivated by the key observation that, in the electric representation, the Hamiltonians of the continuous U (1) group and the discrete Z 2L+1 group after truncation (l < L) are equivalent. The group Z 2L+1 consists of 2L + 1 elements, thus the parameter L indicates the size of the Hilbert space. In particular, the rotatorŝ R j and lowering (raising) operatorsP The only difference between the truncated U (1) group and untruncated Z 2L+1 group is the cyclic property of the lowering (raising) operator, which transforms |L into |−L (and vice versa). However, after a truncation of Z 2L+1 with l < L, this property is lost, meaning that Eqs. (19) and (20) correspond to each other and the two truncated groups become equivalent. For now, consider the Hamiltonian which employs the complete Z 2L+1 group. Importantly, the relations in Eq. (20) resort to a compact description of the electric field since the spectra of the rotators and strings are constrained to the compact interval [−L, L]. We now introduce the following replacement rules for these operators, which reassemble Fourier series expansions. Crucially, this replacement is exact, i.e. there is no truncation of the Fourier series. Employing the fact that the spectrum ofR is discrete and takes integer values, the periodicity of the trigonometric functions can be exploited, which allows one to perform a summation over all coefficients where the sine (cosine) is equivalent. Hence, a finite number of 2L coefficients remain, which take the form Here, ψ k (•) is the k-th polygamma function. Let us further remark that these rules can be extended to higher powers in the variablesR than considered in (21).
This replacement turns out to be convenient for the basis transformation explained below. Introducing the convention |r = |r 1 |r 2 |r 3 and recalling Eq. (20), the electric contribution from Eq. (17) readŝ Note that we use the notation We are now in a position to perform the switch to the dual basis. As shown in App. C, for any γ ∈ N, the discrete Fourier transformF 2L+1 diagonalises the lowering operators aŝ Hence, by applying the discrete Fourier transform to the total Hamiltonian we diagonalise the magnetic contributions, while sacrificing the diagonal structure in the electric part, i.e.  (27) Note that we introduced the superscript (b), which refers to the magnetic representation of the Hamiltonian. Using this representation, computations in the weak coupling regime g 1 can be performed efficiently, as a truncation l now chooses a cut-off for the magnetic field energy.
However, the parameter L now affects the accuracy of the simulation. In fact, while L is completely irrelevant in the electric representation (truncated U (1) and truncated Z 2L+1 are equivalent), it strongly influences the results derived in the magnetic representation. While examining the relationship between L and l in more detail in Sec. 4, we qualitatively discuss our procedure to simulate the U (1) group with the two representations of Z 2L+1 in the following. To be more precise, for any g we might always formulate a sequence of approximating representations for any quantum state of the system in the computational basis defined by |r , i.e., g, r)|r .
Here, p (e) denotes the expansion coefficients in the electric representation, with the subscript indicating the group to which they are referring to (no subscript stands for the truncated Z 2L+1 ). The first approximation in Eq. (28) is due to the transition from U (1) to the Z 2L+1 group, while the second approximation represents the truncation from (2L+1) 3 down to (2l+1) 3 states.
The same scheme exists for the magnetic representation, where the weights p (b) (g, r, L) are used for the state |ψ (b) (g, L) . In this case, however, the choice of L is important. While the truncated electric representation directly corresponds to the truncated and compact U (1) formulation, the completely compact formulation employed in the magnetic representation is crucially affected by the level of discretisation L. This relation is examined further in Sec. 4.2, where we study the convergence of the two representations to U (1) for intermediate values of the coupling g. Hence, in the remainder of this manuscript we consider the completely compact formulation for the magnetic representation only and resort to the compact formulation for the electric representation, i.e. to Eq. (17) for the case of pure gauge.
The interplay of the parameters L and l can be intuitively understood by employing a geometrical illustration. In Fig. 2, the black circles represent the continuous U (1) group, which is approximated by 2L + 1 possible states (blue lines) of the Z 2L+1 group. For l = L, we faithfully describe the untruncated Z 2L+1 group, and use both the solid and the dashed blue lines in the figure. By choosing l < L, we select the states marked with solid blue lines that lie symmetrically around |r = 0 and achieve a binned approximation of any continuous p U (1) lying in the grey area. Furthermore, for any fixed l, the parameter L controls the spread of the available basis states (or bins) around |0 . Since we are interested in the convergence of the truncated Z 2L+1 to U (1) which occurs for L → ∞, we only consider the 2l + 1 states that are important for the dynamics. In particular, we disregard cyclic effects from the lowering operatorP that are a distinctive feature of Z 2L+1 with respect to U (1) [see Eq. (20)].
As an example of this relationship between l and L, consider the two distributions drawn from U (1), represented by the red and green dashed lines in Fig. 2. Clearly, the combination L = 3, l = 1 is insufficient to approximate the broad red distribution. Hence, we increase l to completely cover the target distribution within the grey shaded area. A reduction of L is also a practicable option, especially given the situation where l could not be increased further because of a lack of computational resources. By raising L instead, our binned approximation has a higher resolution around the zero state. For the more localised green state, therefore, it is advantageous to choose the higher value L = 6 when l = 2 is an available option. In fact, the combination L = 6, l = 1 leads to a worse approximation of the green state than the choice L = 3, l = 1, which is therefore preferable when l is limited to unity.
To that end, it is clear that the interplay of L and l represents a crucial point when estimating results and their error due to the performed discretisation. Note that the spread of the distribution determines the value of the free parameter L, while l will, for all practical purposes, be limited by the amount of physical resources, e.g. by the memory of a classical computer or the number of available qubits in a quantum simulator.

Performance and application of the new approach
In the previous section we outlined the transformation from the electric to the magnetic representation, suited to describe the strong and the weak coupling limits, respectively. Here, we develop a protocol which allows for assessing convergence of the truncated representations.
First, we qualitatively describe the system's behaviour at different values of the bare coupling, which will be useful for motivating the protocol. Later, we consider the non-asymptotic cases, where it is not known whether a given truncation is sufficient to describe the considered system with the desired precision. In our example dealing with a pure gauge U (1) theory, this happens for g ≈ 1, where both representations might be inaccurate. Finally, we discuss convergence in the weak and strong coupling regimes, respectively and apply our protocol to estimate the plaquette expectation value. In particular, we consider the interplay between the parameters l and L, which plays an important role when g 1. From now on, we work within a unit lattice spacing, i.e. a = 1, but emphasise that this represents no restriction for the following results.

Phenomenological analysis
As mentioned above, the Hamiltonians in the electric and magnetic representations are related via the Fourier transform, i.e. the magnetic and electric fields are the dual of one another. This relation consequently holds true for the eigenstates, illustrat-ing the difficulty of expressing the ground state in the extremal regimes via the representation in which the dominant term of the total Hamiltonian is nondiagonal. For example, in the regime g 1 the ground state is either determined by the bare vacuum |0 or a superposition of all basis elements, depending whether the electric or magnetic representation is employed (both cases Z 2L+1 , l = L) In other words, the coefficients p(g 1, r) = (2L + 1) −3/2 in Eq. (28) represent a uniform distribution assembling an equally weighted superposition of all basis states. This demonstrates the issue when truncating the Hilbert space by choosing l < L and motivates the choice of switching to the magnetic representation. Note that in the limit g 1 where the electric Hamiltonian is dominant, the roles of the two representations are interchanged.
While it is known that in the limit g → ∞ (g → 0) the ground state in the electric (magnetic) representation is the vacuum |0 , it is not clear what happens when g deviates from these limits. In the following, we employ perturbation theory to estimate the required resources in order to describe the system. For any value of the bare coupling g, we determine the minimal truncation l and resolution L which suggest a high agreement with the untruncated and U (1) theory.
Throughout the whole range of g the system's ground state is center-symmetric, meaning that p (e) (r) = p (e) (−r) in Eq. (28). This follows from the fact that the total Hamiltonian in Eq. (17) is per-Hermitian [61] (Hermitian with respect to the secondary diagonal; higher excited states can also be center-antisymmetric). One can hence infer that the spread of the distribution |p (e) (g, r)| in the electric representation is centred around |0 and decreases with g, again motivating the developed basis transformation of the Hamiltonian. Equivalently, the same holds in the magnetic representation where the centersymmetric ground state becomes less localised by increasing g.
Employing the magnetic representation, we estimate the influence of the electric Hamiltonian with perturbation theory. For g → 0, the unperturbed ground state is |GS (b) (g = 0, L) = |0 , while the first order correction |GS (b) corr takes the form Here, we require L 1 for the inequality, and we introduced the unperturbed ground state energy E 0 = −4/g 2 . Note that the chosen truncation l determines the maximal length of r as |l| = √ 3l, while the population in the states |r is proportional to (gL/|r|) 8 . The upper bound on the population in each |r allows one to determine the part p r>l of the population that is left is out by the truncation at l, which yields p r>l ∝ g 8 L 8 /l 5 . Hence, in order to cover the whole distribution by our truncation, we require l 5 > g 8 L 8 such that large |r| states which are not covered by the truncation are only marginally populated. Respectively, if the truncation l is fixed, we infer that a resolution change L ∝ g −1 is required. In fact, if g 8 L 8 /|r| 8 takes large values for all r, the chosen discretisation is not able to capture the spread of the true distribution, i.e. one would encounter the situation illustrated in the first row of Fig. 2.
This can be intuitively understood by observing that the transition amplitudes between the ground state and states |r induced byĤ (b) E are suppressed by the respective energy gap, i.e. the denominator in Eq. (30). The gap itself is controlled by L, which is a direct consequence of Eq. (27).
The analogous calculation for the electric representation yields l > g −2 / √ 3 which is independent of L. Here, the energy gap is not affected by L and hence deviations from the continuum result have to be associated with a truncation l which is insufficient for the state one aims to approximate.
Concluding, decreasing (increasing) g requires additional computational resources in the electric (magnetic) representation.

Fidelity and convergence of the two representations
This section is devoted to a convergence analysis, which examines the agreement between the two representations.
Although we have developed a scheme that allows one to represent, discretise and truncate the Hamiltonian in the weak coupling regime, the optimal choice of the parameter L is not clear a priori. Clearly, l should usually be chosen according to the availability of the computational resources, which then determines the most suitable L depending on the bare coupling. Furthermore, there is an uncertainty regarding which representation to choose if one is not explicitly considering one of the extremal regimes, g 1 and g 1.
We first develop a criteria to estimate the agreement of the two representations. Therefore, we employ their relation via a unitary transformation and define the Fourier fidelity F f with respect to the same state derived in both representations, e.g. an eigenstates belonging to the same eigenvalue of some observable, such as the ground state derived in both (truncated) representations for a fixed value of g. We write F f as where the Fourier transform is truncated, i.e. the indices of the sums are limited by ±l instead of ±L 1 . Due to the truncation, the features captured in both states are not necessarily equivalent which results in low values of the Fourier fidelity. Vice versa, high values indicate that -for the considered state -the representations are equivalent and yield the same result. Note that this further suggest that the result is close to the hypothetical one derived within the untruncated theory, since the unification of both representations nearly covers the total Hilbert space 2 . Clearly, a low Fourier fidelity is not the only decisive criteria whether a derived result is robust against changes in l or L, especially in the extremal regimes of the bare coupling where the truncation effects of the non-appropriate representation are severe. We thus employ the so called sequence Fidelity (33) Here, µ = e, b indicates considered representation while L is only present in the magnetic case (in the electric representation we can use the truncated U (1) model). Since the truncated models converge to the untruncated U (1) model in the limit l → ∞, high values of F s indicate, under a suitable assumption, that the chosen truncation l is able to capture the whole distribution of the wave vector (as for the case l = 2, L = 3 in Fig. 2). Such a conclusion can not be drawn in the case where the distribution is multimodal with disjoint fractions that lay outside the covered space. Then, the sequence fidelity yields high values for subsequent values of l but would not for larger differences of the considered truncation. Nevertheless, this represents a common issue present 1 Note that this is a consequence of the truncated Hilbert space. However, the operator in (32) not being unitary is no limitation, since the states in Eq. (31) could be embedded in the Hilbert space required by the full Fourier transform. Here, the coefficients for each basis element cn with n > l are set to zero in both representations. 2 Recall that under the Fourier transform, local features are transformed into global ones and vice versa, e.g. a Gaussian is transformed into a Gaussian with inverse width. in all approaches employing truncation techniques that lack the exact true solution.
Let us now return to the ground state of the pure gauge model. Due to their diagonal forms the electric (magnetic) representation yields more accurate results in the strong (weak) coupling regime, however there is no intuition for the intermediate regime g ≈ 1. We hence calculate the Fourier fidelity to obtain an indicator whether results obtained via the different representations at finite values of l are compatible, that is, whether the chosen truncation is enough to capture the local and non-local properties of the state vector. Fig. 3(a) illustrates the Fourier infidelity 1 − F f (l) of the ground state for different values of g −2 . The global maximum of F f arises in from the compromise of having truncation l and resolution L big enough to both contain and resolve the details of the state's distribution. For example, in Fig. 2, it becomes clear that an increase in resolution L reduces the available domain to accommodate a distribution with too high spread if l is not increased accordingly. This relation is the origin of the kink in the Fourier fidelity appearing for larger l in Fig. 3(a), from where the decrease in the Fourier fidelity is solely attributed to an increase of the resolution. Note that for l = 10 we exceed a fidelity of 99.99%.
In the remainder of this section, we will focus on the strong and weak coupling regime where the Fourier fidelity is no meaningful quantity due to the inability to express the state within a truncated basis in both representations. For the electric representation, the sequence Fidelity has a simple interpretation (L is absent here) as it quantifies the overlap between the ground state obtained within different truncations. Since the energy spectrum is fixed and does not depend upon L, a unit value of F (e) s (l) implies that the considered state is unaffected by an increase in l. This suggests that higher truncations do not improve the result and that the model converged to the untruncated U (1) ground state, which can be further motivated by examining the behaviour of 1 − F (e) s (l) in Fig. 3(b). As expected, in the strong coupling regime the sequence Fidelity approaches unity, indicating convergence to the untruncated model, where it is helpful to recall that the ground state in this limit is given by a single basis state, |0 . Approaching the intermediate regime g ≈ 1, F (e) s (l) reduces to a l-dependent constant value, which indicates that the truncation is insufficient to describe all features of the ground state appropriately.
In the magnetic representation, the situation is substantially more complicated, since the approximation of the continuous U (1) group with the discrete Z 2L+1 group introduces the intricate interplay of l and L. As mentioned above, higher values of L allow for a better local approximation of the state which comes at the expense of the tails, which are cut off if l is too small (see Fig. 2). In terms of the sequence fidelity F (b) s (l, L), this implies that for each value of l there exists a unique optimal value L opt of L. Let us stress that this is only true for the ground state of the pure gauge theory considered here. In a more general setting, possibly including matter and higher excited states, F Another complication is given by the fact that L opt does not necessarily corresponds to the global maximum of the sequence fidelity. In particular, a freezing effect can occur for highly localised distributions, where the resolution L is insufficient to capture any of its features. Consequently |ψ (b) (l, L) and |ψ (b) (l + 1, L) are practically the same state and thus yield high values of the sequence fidelity in Eq. (33). In the scenario examined here, the freezing mechanism can be observed in the weak coupling regime, where the ground state is highly localised around |0 . If L is too small, i.e. the bin belonging to the latter state is to wide, all population is accumulated there and the state does not change if g is decreased while L is kept constant. However, it is possible to identify the freezing effect by an educated interpretation of F (b) s . Fig. 3(c) illustrates that the sequence infidelity s (l, L opt ) in both regimes saturates at a l-dependent value. Analogous to the electric representation, it saturates in the strong coupling regime (g 1), however the saturation for weak coupling stems from the limited ability to approximate a continuous approximation with a fixed number of discrete levels. To be more precise, for every l the optimal L opt is chosen as the best compromise of resolution around |0 and a proper representation of the tails of the distribution. In Sec. 4.1, we demonstrated that L opt increases as g is decreased, which we now confirm numerically in Fig. 3(d) (see App. F for more details). Note that as soon as L increases, it does so as L ∼ g −1 , supporting the perturbative results before. Physically speaking, L increases with g −1 since the spread of the population distribution in the magnetic representation decreases, and thus more resolution nearby the state |0 is required.
The black dashed line in Fig. 3(c) corresponds to the global maximum of F (b) s (l = 1, L) for all g −2 . It does not saturate and vanishes in the limit g −2 → ∞. Comparison with the black dashed line in Fig. 3(d), which indicates that L opt ≡ l + 1 reveals this as a characteristic of the mentioned freezing effect.
Concluding, both the Fourier and the sequence fidelities in Eqs. (31) and (33) are two tools to assess the convergence of and agreement between the two representations. While the sequence fidelity must be applied in the extremal regimes, the Fourier fidelity yields a valuable quantification of the combined capabilities of the two representations for intermediate values of the bare coupling.

Estimation of
We now apply the tools developed in Sec. 4.2 to calculate the expectation value as defined in Eq. (16). The value of with respect to the system's ground state is an important quantity in LGTs, as it can be related to the running of the coupling [23].
In the absence of dynamical matter, the total Hamiltonian solely consists of the two gauge field contributions. Therefore, we may determine a value g m separating the regimes where either of the respective representations is advantageous.
Let g m be the value of g for which the Fourier fidelity in Eq. (31) is maximal with respect to the ground state, i.e. (34) Since the electric (magnetic) representation shows exceeding performance in the strong (weak) coupling regime, we can assume that for a given truncation l, the best approximation is achieved by considering the electric representation in the range g ∈ [g m , ∞) and the magnetic one for g ∈ [0, g m ] (compare also Sec. 4.2 and Fig. 3). representation. In the latter, we obtained the L opt values that have been used for plotting via the sequence fidelity as described above. From here, the true curve as it would be obtained from the untruncated U (1) theory can be estimated via the asymptotic values of the different representations when the truncation l is increased, since in the limit l → ∞ both representations converge to the full theory. We exemplify such an estimation with the inset in Fig. 4(a), that contains the results for different l at g −2 = 10. The convergence can be clearly observed, and both representations yield the same result up to the fourth decimal at l = 10 ( = 0.9572 ± 0.0001). Note that this convergence is not necessarily monotonic. However, in the extremal regimes, we observe that the expectation value of increases with the truncation l when employing the electric representation, while it decreases with the magnetic one, for which we will provide analytical arguments in App. D.
To summarize this section, we recall that a naive approximation of U (1) with Z 2L+1 (with L fixed) leads to dramatically increasing computational costs when working on a wide range of g-values. As explained intuitively in Sec. 3, the problem originates from the fact that Z 2L+1 converges not uniformly but pointwise to U (1). For fixed resolution L and fixed computational resources l, there is always a coupling g small enough such that the Z 2L+1 description displays freezing and hence cannot approximate the U (1) continuum physics accurately. This can be understood by noting that the magnetic field Hamiltonian is gapless in both the continuum theory and in the U (1)-lattice description, but gapped in the Z 2L+1formulation. For fixed L and decreasing g, the offdiagonal elements in the HamiltonianĤ E decrease with respect to the energy gap inĤ B (as explained in more detail in Sec. 4.1). If the energy in the system becomes comparable to the gap, the difference between Z 2L+1 and the true gauge group U (1) becomes noticeable, which leads to the freezing effect (see Fig. 3). Crucially, working with a value of L suitable for the regime g 1 will lead to exploding computational costs, i.e. will require very large values of l, in the intermediated coupling regime g ≈ 1 to capture the relevant physics there. Our solution to this problem is the dynamical adjustment of the parameter L with the coupling g, that allows us to approximate U (1) well for a wide range of couplings while including only a minimal number of states in our simulation (see Fig. 2).

Generalisations: Dynamical matter and arbitrary torus
In the following, we extend the results presented in Sec. 4 by including staggered fermions in the numerical simulations. In particular, we show that matter does not introduce any fundamental complication for the completely compact formulation introduced in Sec. 2. Moreover, to pave the way for further developments in the field, we derive the Hamiltonian for an arbitrary number of plaquettes on a torus with matter and periodic boundary conditions, and explain how to include static charges.

A single plaquette with charges
Since the completely compact formulation only affects the gauge fields, the inclusion of matter is straightforward. Recall first the electric Hamiltonian in Eq. (13) and the substitution rules in Eq. (21). By using the relations for the Fourier transform derived in App. C, the magnetic representation of the electric term in Eq. (13) is found to bê For the sake of clarity, we defined the shorthand no-tationsK The magnetic field HamiltonianĤ remains the same as in Eq. (27), since it does not involve fermionic terms. However, the kinetic Hamiltonian in Eq. (15) is modified in the presence of matter, yieldinĝ +H.c. |r r|. (37) In order to simulate fermionic matter, we recall the Jordan-Wigner transformation [62] where the vectorial relation l < n is defined by (0, 0) < (0, 1) < (1, 1) < (1, 0) to satisfy the fermionic commutation relations. In higher dimensions, it might however be useful to consider alternative approaches, such as fermionic projected entangled pair states or the elimination of the fermionic matter [63,64]. While we do not insert these equations into Eq. (37), we remark that the mass Hamiltonian in Eq. (5) is simplified tô which is independent of the chosen representation. Since these simulations are costly, i.e. the dimension of the truncated Hilbert space is given by 2 4 (2l + 1) 5 (four charges, three rotators and two strings), we estimate the plaquette expectation value employing a harsh truncation of l = 2, while fixing L to the optimal values L opt found in Sec. 4.2 for the pure gauge case. This is a reasonable assumption, since strings and fermionic matter only play a role in the intermediate regime. Therefore, we can recover our previous results for the pure gauge theory in the strong and weak coupling limits, and focus our attention to the differences where g ≈ 1. We further introduce the mass and kinetic energy parameters as m = κ = 10.
In Fig. 5, we display the ground state expectation value as a function of g −2 , together with the Fourier infidelity 1 − F f (l). While the asymptotic regimes g 1 and g 1 show no qualitative difference if compared to the pure gauge case, the situation changes in the intermediate regime. There are novel features in both the electric and magnetic representations, such as the appearance of a negative dip. Nevertheless, we stress that conclusions drawn from this plot have to be taken with care, as the employed truncation limits the Fourier fidelity below 90%.
Concluding, we demonstrated that our method is suitable to tackle simulations with matter, and can be scaled up to more complex systems. A detailed analysis of novel effects and an accompanying study of the convergence is beyond the scope of this manuscript and left for future works.

Hamiltonian for an arbitrary torus and charges
Here, we generalise the Hamiltonian of a single periodic plaquette to any two-dimensional lattice with periodic boundary conditions. As shown in Fig. 6, we extend the strategy used above to a torus of size (N x , N y ). By removing redundant DOF, we obtain the effective Hamiltonian in terms of two strings and N x N y − 1 rotators. As before, we indicate each plaquette with its bottom-left site n = (n x , n y ), where n x (n y ) = 0, . . . , N x − 1 (N y − 1). In addition, the rotator associated with the plaquette n is denoted bŷ R n , and the two strings byR x andR y (see Fig. 6). This leads to N x N y pairwise expressions for the electric fields, E n,ex = δ ny,0Rx +R n −R n−ey +q n,x , E n,ey = δ nx,0Ry +R n−ex −R n +q n,y , (40) where δ n,m = 1 for n = m, and zero otherwise. Moreover,q n,x andq n,y are the electric field's corrections due to the presence of dynamical charges, in accordance with Gauss' law. Since there are several ways to implement Gauss' law, a possible choice forq n,x andq n,y is (see the green lines in Fig. 6) whereq n is the charge operator as defined in Eq. (7). Note that also in this general case it is convenient to explicitly fix one of the rotators to zero, for instancê R (0,Ny−1) = 0.
Moving to the kinetic term, we employ the string convention presented in Eq. (41), which yields the replacement rules (for details, see App. B.2) We remark that these equations are also valid for particles following bosonic statistics [65], where the charge operator for site n is defined aŝ In this case, the only required modification concerns the mass Hamiltonian, which becomeŝ Importantly, employing the relations derived in App. C directly allows for the transformation between the electric and magnetic representations. As a final remark, we highlight that including static background charges can be done by explicitly adding Q n into the electric termĤ E in Eq. (43). This can be achieved with a simple redefinition of the charge operatorsq n →q n + Q n .

Conclusions and outlook
We developed a new strategy for studying gauge theories. Our method is suited for simulations of fundamental particle interactions in all coupling regimes on near-term quantum computers (see [23]), as well as on classical devices. As a testbed, we applied our method to the lattice Hamiltonian formulation of QED in 2+1 dimensions.
The key insight is the approximation of the continuous U (1) gauge group with finite truncations of the Z 2L+1 group, where L ∈ N can be arbitrarily large and is scaled with the value of the bare coupling g. This strategy allows us to work with fixed computational resources, i.e. including only a constant number of states in our simulation, for any value of g. At weak couplings we truncate the gauge fields in the magnetic representation of Z 2L+1 , while the truncation is performed in the electric representation for strong coupling. This choice ensures a minimal truncation error Figure 6: Periodic torus with charges. We extend the construction of the periodic plaquette to a generic torus. We fix the rotator at (0, Ny − 1) to zero and choose the links' corrections to the electric field introduced by charges in accordance with the green dotted line. In particular, for any chargeqr, we connect the origin to the site r by moving first horizontally and then vertically.
in any regime. We benchmarked this novel regularisation scheme by computing the expectation value of the plaquette operator on a small periodic lattice, with and without dynamical matter, and estimated the accuracy of the computation.
Since our methods allows us to work at all values of g and therefore at arbitrarily small values of the lattice spacing a, it provides the perspective to access, in principle, non-perturbative physics close to the continuum limit.

Quantum simulations:
With regard to simulations of LGTs, there are two different lines of work. The first research line studies lattice models in their own right. Lattice gauge theories are for example relevant in condensed matter physics or can be interesting per se. The second line of research considers lattice gauge theories with the aim to study the underlying continuum theory which describes for example fundamental particle interactions and the standard model. For simulations of the latter, it is indeed of crucial importance that one is able to take the continuum limit of a lattice theory. In the field of quantum simulations, this challenge is mostly unanswered.
So far, the only practical route to approach the continuum limit were analog quantum simulators using infinite degrees of freedom to represent the gauge fields. Neutral atoms in optical lattices offer a very good solution in this respect, as the gauge field can be represented as a spinor condensate while the charges are identified with moving fermions, and the gaugematter interaction by spin changing collisions [66][67][68][69][70][71]. A basic building block [32] and a one-dimensional proof-of-principle experiment [33] that exploit some of these ingredients have been already performed. Beyond one dimension, such approaches -although very promising -are fundamentally limited since magnetic (plaquette-type) interactions are realised via higherorder interactions. This results in low effective coupling strengths and thus in extremely challenging experimental requirements.
While the analog approach based on bosonic degrees of freedom outlined above has the potential to achieve the continuum limit, the experimental realisation of two-dimensional theories involving magnetic terms is currently out of reach experimentally. This type of interaction is easier to realise digitally [72][73][74] and in qubit-based platforms [45,75], such as trapped ions, Rydberg atoms and superconducting qubits. These simulation strategies however, currently lack the feature to reach the continuum limit.
Our new framework provides a route towards reaching the continuum limit in quantum simulators on different platforms and therefore opens a new perspective for meaningful simulations that address physical (i.e. continuum) phenomena that can be related to experiments in high energy physics. It will be interesting to explore the use of interacting chains of spins larger than spin-1/2 [76] to simulate gauge fields and to investigate the use of ultracold fermions, which will allow one to simulate fermions on a system that naturally displays the right quantum statistics. As a first step towards proof-of-concept demonstration using our new method, we show how to apply our scheme to current ion-based quantum hardware [23].

Tensor network calculations:
Resource minimisation is especially relevant for classical simulations based on tensor network states. Our regularisation scheme can also be reinterpreted as a coupling dependent variational ansatz for the gauge-field states, where the ratio between the truncation parameter and the dimension of the discrete group l/L is optimised. With this perspective in mind, our scheme could be combined with other variational approaches, e.g. with the method put forward in [77], with the aim to extend the continuum limit calculations beyond one dimension or for addressing toy models of high-energy physics, such as CP (N − 1) theories [78].
Extensions to higher-dimensional non-Abelian gauge groups: We note that both the minimal Hamiltonian formu-lation (in which redundant degrees of freedom have been removed) and our new regularisation scheme can be extended to higher dimensions and to non-Abelian gauge theories, also beyond the Hamiltonian approach. The solutions we proposed here are interesting for Lagrangian-based formalisms like the tensor approach [79,80] or for the development of novel Monte Carlo approaches to avoid the sign problem [81,82]. From a geometric perspective, once U (1) is identified with S 1 , and the magnetic vacuum with the north pole [see Fig. 2], our regularisation scheme at weak couplings can be interpreted as a lattice discretisation of a circle around the north pole. One could exploit the map of SU (N ) groups to higher dimensional spheres considered in [81,82] and repeat a similar procedure. An alternative discretisation of non-Abelian groups for classical and quantum simulations has also been considered in [83][84][85][86].
In conclusion, our work opens new perspectives for resource efficient Hamiltonian-based simulations and provides a concrete step in the 'quantum way' to nonperturbative phenomena in high energy physics.
where we might define the new effective mass m = M/α and effective kinetic energy scale κ = K/α valid for the rescaled fermionic fields as employed in the main text. Moreover, this rescaling implies a rescaling of the charge operator, which we define aŝ where q = Q||Ψ † nΨn || max = Q/α ∈ Z and q = 1 in the main text. Note that this operator is dimensionless as well, since it is required to be of the same dimension as the electric fieldÊ.
B Hamiltonian in the link formalism and link-to-rotator translation rules B.1 Effective Hamiltonian for a single plaquette in the link formulation As explained in Sec. 2, there are several ways to express the Hamiltonian of one or more plaquettes. In the main text, we employed electric loops for parametrizing the physical states and defining the corresponding operators, rotators and strings. These represent a natural choice since the rotator's lowering operators are directly identified with the plaquette operatorsP n [see Eq. (2)]. Here, we derive the Hamiltonian of a single plaquette with staggered fermions in terms of the linksÊ n,eµ . In particular, we choose three out of the eight electric fields on the links, and express them in terms of the others to minimise the number of degrees of freedom. In this appendix, we consider the compact U (1) formulation of the QED lattice Hamiltonian reviewed in Sec. 2.1. The completely compact Z 2L+1 QED formulation for both the electric and the magnetic representation can be obtained following the procedure outlined in Sec. 2.2. Generalisations to multiple plaquettes are straightforward.
Employing Eq. (3), we can directly assess Gauss' laws for the single plaquette in Fig. 1 where µ = x, y and f n,eµ is a function of both rotator and string operators. Depending on the particle pair to be created or annihilated, there are four distinct rules for building up the functions f n,eµ . These are shown in the four corresponding panels displayed in Fig. 7, where we use a yellow (green) circle to indicate an (anti)fermion, and light blue arrows for the corresponding link in which the electric field is required to raise [Û † n,eµ in Eq. (63)]. Employing the notation presented in the legend of Fig. 7, we can directly build . The particles (green and yellow circles) imply the creation of electric fields on the links marked with the corresponding dashed arrows. Only in (c) these contributions are enough to create the gauge field while maintaining gauge-invariance in all other links. Otherwise, the gauge field has to be created by annihilating the plaquette operators marked with red circular arrows, which also counteracts the action of the particles' charges. Note, that the field on a link is effectively unchanged if it is modified by an equal number of arrows in both directions, which in (b) and (d) requires the introduction of the strings (grey and pink) when the particles are created on the boundary condition of the lattice. the functions f n,eµ in Eq. (63) and recover the rules presented in Eq. (42).
As an example, consider Fig. 7(a), which describes all cases in which a pair is created at locations n and n + e x , with n x = N x − 1. If n y = 0, the electric field is automatically corrected by the chosen charge strings that ensure Gauss' law (dotted green and yellow lines in Fig. 7; see Sec. B.1 and Sec. 5.2), implying f (nx,0),ex = 1. Otherwise (n y = 0), to increase the electric fieldÊ n,ex on the link between the two charges, we lower the rotatorR (nx,ny−1) below by applyingP (nx,ny−1) . This, however, affects the electric fields on all other links formingR (nx,ny−1) . While the vertical ones are taken care of by the charge strings, the bottom link is not (unless n y = 1). As graphically explained in Fig. 7(a), to increase only the desired link, we can lower all rotatorsR (nx,ry) with r y = 0, ..., n y − 2 below, yielding f (nx,ny),ex = ny−1 ry=0P(nx,ry) . The remaining panels in Fig. 7 further illustrate the cases of pairs that are connected by vertical links and pairs that are created on links closing the periodic boundary conditions, where we require the stringŝ R x,y . By following the same procedure above, it is then possible to determine the functions f n,eµ for all allowed choices of n and µ = x, y, yielding Eq. (42).

C Diagonalisation of the magnetic gauge field Hamiltonian
The magnetic HamiltonianĤ B is composed of the lowering and raising operatorsP j ,P † j , j = 1, 2, . . . , N − 1, where N is the total number of plaquettes. In the Z 2L+1 group, these operators are the so-called cyclant matrices, which can be diagonalised exactly. Before truncation, the lowering operators are defined according to Eq. (20), For the sake of simplicity, we drop the index j and note that the procedure is equivalent for all subsystems. The spectrum of the lowering operators is while the corresponding eigenvectors are given by with k = −L, . . . , L. Hence,Û is diagonalised by the matrix which is the discrete Fourier transform. Hence, it is straightforward to show that Here, we use r = (r 1 , r 2 , . . . , r J ) T and γ = (γ 1 , γ 2 , . . . , γ J ) T , while we waived to denote that the Fourier transform is now understood as the product of the Fourier transforms in the separate N −1 spaces. Note that, in particular (P γ ) † =P −γ and therefore: Let us finally remark that these relations hold interchangeably forP in the rotator formalism andÛ for the electric fields as introduced in Sec. 2.1. employ the following ansatz to obtain the truncated electric field Hamiltonian We hence define the operatorV to study the effects of the truncation.
In the untruncated Z 2L+1 formulation we can decompose anyP j into four terms, whereP j = l rj =1−l |r j − 1 r j | andP j is the rest. Notice thatP j is the truncated operator as defined in Eq. (19). We now have to collect all contributions fromV j in the electric HamiltonianĤ  E by applying the rules P j + P † k →V j + (V k ) † , P j P k →V j (V k ) † +V j (P k ) † +P j (V k ) † . (77) Note thatV has as well a diagonal contribution given by stemming from the ν = 0 terms in Eq. (26). The result of a numerical procedure to find the elements of ground state is shown in Fig. 8(a). We denote the amplitudes of a state obtained with the untruncated Hamiltonian with the red dashed line. Clearly, the truncation shifts population from the high |r| states to lower ones. Note that each |r| can appear multiple times. A crucial fact is shown in Fig. 8(b), where we removed only the cyclic elements from the operators P j , i.e., the first terms in Eq. (76). This cyclic property is the reason why the distributions p (b) (r) of the ground state's coefficients (see Sec. 4.1) are uniform in the untruncated case. Hence, their removal has a large impact on the derived results.

F Numerical determination of L opt
The sequence fidelity calculated in Sec. 4.2 for the ground state of the pure gauge QED Hamiltonian involves and optimization over L. We plot the sequence infidelity varying the parameters l and L for g −2 = 100 in Fig. 9. For any value of l, a kink is clearly visible, corresponding to L = L opt . Notice that this kink does not always correspond to a global minimum, as can be seen from the points characterised by l = 2. As discussed in the main text, this is the signature of the freezing effect. In fact, the global minimum is found for L = l + 1 = 3 (i.e. its minimal value), where the resolution is insufficient for both l = 1 and l = 2 to capture the distribution of the untruncated ground state. By increasing g −2 , the position of the kink is shifted to higher values of L. In particular, L opt takes its minimal allowed value in the strong coupling regime, and starts to increase at g −2 ≈ 5 [see Fig. 3(d)]. This follows from the fact that, approaching the weak coupling regime, the distribution of the untruncated U (1) ground state becomes more and more localised, and the tails less important. Note that the value L opt can be determined by a greed search, starting at the lowest allowed value of L.