A SUPPLEMENTAL MATERIAL

# Non-Abelian Su(2) Lattice Gauge Theories in Superconducting Circuits

## Abstract

We propose a digital quantum simulator of non-Abelian pure-gauge models with a superconducting circuit setup. Within the framework of quantum link models, we build a minimal instance of a pure gauge theory, using triangular plaquettes involving geometric frustration. This realization is the least demanding, in terms of quantum simulation resources, of a non-Abelian gauge dynamics. We present two superconducting architectures that can host the quantum simulation, estimating the requirements needed to run possible experiments. The proposal establishes a path to the experimental simulation of non-Abelian physics with solid-state quantum platforms.

###### pacs:
03.67.Ac, 03.67.Lx, 85.25.Cp,11.15.Ha

Gauge invariance is a central concept in modern physics, being at the core of the standard model of elementary particle physics. In particular, invariances with respect to and gauge symmetries characterize the weak interaction and quantum chromodynamics (1); (2). In this sense, gauge theories represent a cornerstone in our understanding of the physical world and lie at the heart of diverse phenomena, such as the quark-gluon plasma or quantum spin liquids. In condensed matter physics, gauge fields can also emerge dynamically in relation to exotic many-body phenomena, like quantum Hall systems, frustrated magnets, or superconductors (3); (4); (5).

Lattice gauge theories (LGT) are non-perturbative discrete formulations that contribute to the analysis of key features of these models, such as color confinement or chiral symmetry breaking. Starting from the seminal work by Wilson in 1974 (6); (7), LGT have attracted a significant attention across several branches of theoretical physics. In the last decades (8), quantum Montecarlo simulations have achieved unprecedented accuracies in determining the whole hadronic spectrum of the standard model. However, understanding its full phase diagram from first principles, or simulating dynamical processes, remain out of reach of current numerical computations.

Quantum simulators (9) provide a new approach to solve complex long-standing problems in quantum physics. In a quantum simulator, a proper encoding of LGT can allow for the retrieval of information about ground state and dynamics, in a wide range of regimes. Previous works have considered Abelian (10); (11); (12); (13); (14); (15); (16); (17); (18); (19); (20); (21); (22) and non-Abelian LGT in optical lattices (25); (23); (24) (see also (26); (27) and references therein), both as analog and digital simulations (28). In these implementations, matter-gauge interactions are modeled as a second-order process from Hubbard-like interactions, while the simulation of pure-gauge dynamics remains more demanding.

In the last years, superconducting circuits have proven to be reliable devices that can host quantum information and simulation processes (29). The possibility to perform quantum gates with high fidelities, together with high coherence times, makes them ideal devices for the realization of digital quantum simulations (30); (33); (31); (32); (34), previously considered in ion-trap systems (35); (36); (37).

In this Letter, we propose a digital quantum simulation of a non-Abelian dynamical gauge theory in a superconducting device. We start by building a minimal setup, based on a triangular lattice, that can encode pure-gauge dynamics. The degrees of freedom of a single triangular plaquette of this lattice are encoded into qubits. We propose two implementations of this quantum simulator, using two different superconducting circuit architectures, as depicted in Fig. 1. We consider a setup in which six tunable-coupling transmon qubits are coupled to a single microwave resonator, and a device where six capacitively-coupled Xmon qubits stand on a triangular geometry, coupled to a central auxiliary one. We compute experimental requirements necessary to perform the simulation on one plaquette, and provide arguments for scaling to large lattices.

Lattice gauge theories.- LGT are discretized versions of a gauge theory. In a conventional approach of lattice gauge theories, space-time is discretized while ensuring a covariant formulation of the theory. In quantum simulations, there is no direct access to the time direction, that is fixed and continuous. Hence, the equivalent Hamiltonian formulation of lattice models is used. In this case, the space is discretized and the action of a local gauge invariant Hamiltonian characterizes a continuous dynamical evolution. This fictitious asymmetry of the time direction forces us to define a physical Hilbert space or Gauss law of the model. We make use of a set of discretized space points , and operators associated to each link, connecting two adjacent sites . A single operator has two implicit color indices , , that connect the -th component of the fermionic field at position with the -th component of the field at position .

In order to construct a non-Abelian theory built out of these operators, one has to build a Hamiltonian that is invariant under a generic gauge transformation

 u(→x,^μ)→eiθa(→x)τau(→x,^μ)e−iθa(→x+^μ)τa, (1)

where are the generators of the algebra in the corresponding representation, and the relative phase angles. It is straightforward to verify that the following Hamiltonian on a square lattice is gauge invariant

 H=J∑→xTr[u(→x,^μ)u(→x+^μ,^ν)u†(→x+^ν,^μ)u†(→x,^ν)]+H.c., (2)

where the trace is performed over the color indices, has energy units, and the directions and span the two-dimensional lattice. This Hamiltonian is a pure gauge operator, as it does not involve any fermionic operator and it is the minimum instance of a Wilson loop around a plaquette (1); (2).

We focus now on the construction of a non-Abelian quantum link model that mimics this gauge-invariant behavior in a finite-dimensional Hilbert space, suitable for quantum simulations. Notice that the framework of the quantum link models is valid for any compact Lie group (38). We use a four-dimensional Hilbert space or, equivalently, two qubits associated with each link. One can define link operators acting on this four-dimensional Hilbert space with the associated gauge generators that satisfy the algebra

 [Ga(→x),Gb(→y)]=iδ→x→y∑cϵabcGc(→x). (3)

To define the set of operators, we use a quantum link formulation (38), where the Hilbert space of a link is given by a set of bosonic modes that implement the Schwinger representation of the algebra, with , acting on two different sites on the link between the two adjacent vertices and . In this space, we build two sets of right and left generators , , on this finite-dimensional Hilbert space, using the bosonic modes  (39). We define the gauge generators as , where the sum over is taken among all the links of the lattice converging to a single vertex .

The representations of are quasi-real, therefore, ordinary and their dual representation (i.e., particle and antiparticle) are equivalent. Therefore, there are two multiplets with well-defined transformations (39). One can finally define the link operators, acting on the finite-dimensional Hilbert space of one link, in terms of Schwinger bosons on a given link (47); (48); (49),

 Uαβ=⎛⎝c↑L−ic†↓Lc↓Lic†↑L⎞⎠⋅(c†↑Rc†↓Ric↓R−ic↑R). (4)

These operators have the following commutation rules with the left and right operators

 [U(→y,^ν),Ra(→x,^μ)]=−U(→y,^μ)σa2δ→x→yδ^ν^μ,[U(→y,^ν),La(→x,^μ)]=σa2U(→y,^μ)δ→x→yδ^ν^μ, (5)

where we have defined the usual Pauli matrices , see also (39), while is defined as the identity operator. The commutation rules for a general gauge transformation follow in a straightforward way,

 ∏→ye−iθa(→y)Ga(→y) U(→x,^μ)∏→zeiθa(→z)Ga(→z) =eiθa(→x)σa2U(→x,^μ)e−iθa(→x+^μ)σa2, (6)

where the products over and are extended over the whole lattice. This equation shows that quantum link models are formulations of gauge-invariant models. In fact, one can notice how the action of a gauge transformation in Eq. (6) mimics Eq. (1).

A generic state in the local Hilbert space of a quantum link is given by the occupation of the operators . Since the total occupation per link, , is a constant of motion [see Eq. (4)], we restrict to the subspace with one occupied mode per link.

Notice that matter-gauge interactions in dimensions can be derived (39) in a second-order perturbation theory, considering extensions of previous works (50); (51). In the following, we will focus instead on pure-gauge two-dimensional interactions.

Triangular lattice.- The continuum limit of quantum link models and their thermodynamical properties are topics under active research. In related condensed matter models, the importance of the lattice geometry has been shown in fundamental aspect of the phase diagram, as the existence of confinement and de-confinement phases (52); (53). Quantum simulations on large lattices may provide insights on the continuum limit of LGT, and its dependence on the geometry of the underlying lattice. In this spirit, we consider a minimal implementation of a pure invariant model in a triangular lattice, by using triangular plaquettes, as depicted in Fig. 1c. In this case, the pure-gauge Hamiltonian on a single plaquette reads

 HT=−J Tr[U(→x,^μ)U(→x+^μ,^ν)U(→x+^μ+^ν,−^μ−^ν)]. (7)

This interaction corresponds to the magnetic term of a gauge invariant dynamics, which acts on closed loops. We focus on the magnetic term, since the representation of the electric counterpart is trivial.

Since this Hamiltonian commutes with gauge generators, , these are constant of motion, and one can define different gauge sectors with the initial values of the generators. Due to the triangular geometry and having one occupied mode per link, there are no global states that have zero eigenvalue for the gauge generators at every vertex in a single plaquette. In general, the absence of the zero-eigenvalue sector depends on the topology of the lattice, and is avoided, for example, in the periodic lattice of Fig. 2d. The different gauge sectors can be characterized by the eigenvalues of the operator , as in Fig. 2a,  2b, and 2c, see the supplemental material for additional details (39).

The local Hilbert space of a link is four dimensional, and it can be faithfully spanned by two qubits, called “position” and “spin” qubit . In this subspace, it is useful to define the operators , , such that the total Hamiltonian is written as  (39)

 HT= −J{Γ012Γ023Γ031+∑abcϵabcΓa12Γb23Γc31 (8) −∑a[Γ012Γa23Γa31+Γa12Γ023Γa31+Γa12Γa23Γ031]}.

Superconducting circuit implementation.- In order to simulate the interaction of Eq. (8), one can decompose its dynamics in terms of many-body monomials, and implement them sequentially with a digitized approximation (54). In a digital approach, one decomposes the dynamics of a Hamiltonian by implementing its components stepwise, (here and in the following ), for a total of  gates, with an approximation error that goes to zero as the number of repetitions grows. In a practical experiment, each quantum gate will be affected by a given error . By piling up sequences of such gates, for small gate errors , the total protocol will be affected by a global error, which is approximately the sum . This model for error accumulation has been proved in recent experiments (33); (32).

The effects of digitization in the simulation of LGT can be observed in Fig. 3. We have numerically integrated a Schrödinger equation regulated by the Hamiltonian in Eq. (8), choosing as initial state one of the states shown in Fig. 2a, for different simulated phases  (39). We compute both the evolution with the ideal Hamiltonian and with a digital sequence, in which we act with a single monomial at a time, repeating the protocol times. As a figure of merit, we compute the relative deviation from the ideal value of the gauge invariant, where stand for average values computed on the ideal (digitally simulated) state. The deviation goes to zero for large and small phases, defining surfaces at given error tolerance in the - space.

To simulate the pure-gauge interaction in a single triangular plaquette, we first consider a setup in which six tunable-coupling transmon qubits are coupled to a single microwave resonator (55); (56). Each tunable-coupling qubit is built using three superconducting islands, connected by two SQUID loops. Acting on these loops with magnetic fluxes, one can modify the coupling of the qubits with the resonator, without changing their transition frequencies. For further clarifications of the experimental setup involving tunable-coupling transmon qubits, see the Supplemental Material (39). By threading with magnetic fluxes at high frequencies, one can drive simultaneous red and blue detuned sidebands, and perform collective gates (57). Each many-body operator can be realized as a sequence of collective and single-qubit gates,

 UN=ei(π/4)∑i

where the indices and run from to . In this way, one can implement a generic (up to local rotations) -body evolution operator , with , using collective gates, and a number of single qubit gates which is upper bounded by , counting the single qubit gate in Eq. (9) and the rotations necessary to map to any -body operator. The simulation for one digital step will amount to implementing collective gates and a number of single qubit gates which is upper bounded by  (39).

We consider now the architecture of Fig. 1b, in which six Xmon qubits in a triangular geometry are capacitively coupled with an additional central ancillary qubit (58). In this case, the collective interactions can be decomposed and performed with pairwise -phase gates, using the central ancillary qubit to mediate non-nearest interactions. In this way, the quantum simulation of one digital step of the Hamiltonian in Eq. (8) will amount to realize  -Phase gates and a number of single-qubit rotations which is upper bounded by . For additional details on the experimental setup and gate counts, see (39).

To estimate the effects of experimental imprefections on the simulations, we plot in Fig. 4 the dynamics of a single plaquette, starting from a state. The observed periodic dynamical oscillations between the two paired states, depicted in Fig. 2b, are signaled by the periodic behavior of the overlap . The simulation is run for and steps, necessary to observe one flip between the two states and one full oscillation, respectively. On top of these oscillations, we plot fidelity bands that estimate the fidelity cap due to accumulated gate errors in the digital protocol, in the case of the two presented setups in Fig. 1, with collective and two-qubit gates. In order to see coherent oscillations, the fidelity of the two-body gates has to be one order of magnitude better with respect to the collective ones, with an error window of approximately in the collective gate case and a for the -phase gate setup. We assume a factor for the error loss on a single qubit gate. From the plots, it can be seen that, while the digital fidelity increases, the fidelity bands due to experimental errors broaden, affecting the simulation. The intersection of the digital fidelity line with the bands marks a regime dominated by experimental imperfections, for small , and one in which the error of the digital expansion prevails, for large .

The proposed simulation can be extended to larger lattices, where one could analyze dynamical properties of the confinement-deconfinement phase transition, which can be traced back to the breaking of the global center symmetry (59) in groups with a non-trivial center group, like the case. Moreover, such a LGT simulator can be used to perform quench experiments in a ladder configuration, analyzing breaking of gluon strings between pairs of particles and anti-particles (60); (61); (62).

In order to scale the quantum simulation to large qubit lattices, one has to consider that the accumulated gate error does not depend on the size of the lattice (33); (32), and only on the number of gates. The quantum resources necessary to run the simulation scale polynomially in the lattice size, and sub-polynomially in the digital error (63), making the whole protocol efficient. With a current total number of gates that exceeds 1000 (32); (34), we believe that in the near future simulations of LGT on large scales will be feasible. Furthermore, the possibility of performing quantum error correction for digital quantum simulations may drastically increase its effectiveness (64).

AM and LL acknowledge useful discussions with Rami Barends. ER thanks D. Banerjee, T. Calarco, M. Dalmonte, S. Montangero, U.-J. Wiese, P. Zoller for discussions, in the framework of quantum technologies for lattice gauge theories. We acknowledge financial support from Basque Government Grants IT472-10 and IT559-10, Spanish MINECO FIS2012-36673-C03-02, Ramón y Cajal Grant RYC-2012-11391, UPV/EHU Project No. EHUA14/04, UPV/EHU UFI 11/55, PROMISCE and SCALEQIT European projects.

## Appendix A Supplemental Material

In this supplemental material, we perform additional analysis that support the results obtained in the main article.

### a.1 Explicit Su(2) gauge operators in terms of Schwinger modes

In this section, we explicitly derive the link operators in terms of Schwinger bosons. Following (65), we define the operators , as right and left generators, acting on the finite-dimensional Hilbert space of the link,

 Ra(→x,^μ)=∑{α,β}∈{↑,↓}c†αR(→x,^μ)σaαβ2cβR(→x,^μ),La(→x,^μ)=∑{α,β}∈{↑,↓}c†αL(→x,^μ)σaαβ2cβL(→x,^μ), (10)

where are bosons that implement the Schwinger representation of the algebra, and the Pauli matrices

 σ(0)=(1001); σ(1)=(0110); σ(2)=(0−ii0); σ(3)=(100−1). (11)

The bosonic operators act on two different sites on the link between the two adjacent sites and and obey the commutation rules

 [Ra,Rb]=i∑cϵabcRc,  [Ra,Lb]=0,  [La,Lb]=i∑cϵabcLc, (12)

taking into account that they commute on different links (i.e., different , ). With the Schwinger representation, there are two multiplets with well-defined commutation relations

 [cαL,La]=∑γσaαγ2cγL,⎡⎣⎛⎝∑βσyαβc†βL⎞⎠,La⎤⎦=∑γσaαγ2⎛⎝∑βσyγβc†βL⎞⎠,[c†αR,Ra]=−∑γc†γRσaγα2,⎡⎣⎛⎝∑βcβRσyβα⎞⎠,Ra⎤⎦=−∑γ⎛⎝∑βcβLσyβγ⎞⎠σaγα2. (13)

Given these two multiplets, we can define the link operator as

 Uαβ=cαLc†βR+∑μνσyαμc†μLcνRσyνβ. (14)

We stress that these quantum link operators, or gauge fields, are associated with links between lattice points, transporting color information between the two points they connect. The two indices of the link operator are identified with the two ends of the link . The left (right) index is associated with the beginning (end) of the link as shown in Fig. 5. The gauge transformation acts on according to

 U(→x,^μ)=ei∑aθa(→x)σa2U(→x,^μ)e−i∑aθa(→x+^μ)σa2. (15)

Since in general and are different, our gauge-invariant equations will require invariance under right multiplication and left multiplication separately. An arbitrary gauge transformation can be built from individual gauge transformations at the point of the lattice. Therefore, it is sufficient to consider just a gauge transformation at the lattice site . Every link emanating from the point is affected by the gauge transformation generated by

 Ga(→x)=∑^μLa(→x−^μ,^μ)+Ra(→x,^μ), (16)

and the total color carried by a link is just the difference: .

The sum over non-Abelian electric fields emanating from a single point of the lattice is the lattice analog of the electromagnetic divergence of the electric field, . However, because the electric field itself varies from to along the link, there is an additional contribution given by the total color carried by the link. Thus, the generator is recast into: . The gauge invariance of the physical sector is defined by , i.e., at every vertex the state transforms locally as a singlet of . If we identify the local color density , then is the usual Gauss’ law. These are properties of non-Abelian lattice gauge theories, independent of their representation via quantum link models.

### a.2 Local Hilbert space of Su(2) quantum links

The local degrees of freedom in a link are encoded in the quantum state of four different bosonic modes . Different representations of this local Hilbert space are given by the occupation of these modes. The Hilbert space with one excitation in total is four-dimensional, and it can be decomposed in a spin degree of freedom with the symmetry group and a position degree of freedom.

A gauge invariant interaction commutes with the local generator of the at any vertex. As a consequence, the total angular momentum at any vertex of the lattice is a conserved quantity of the dynamics. In the simplest example of two quantum links, i.e. just one vertex, there are only two possible gauge sectors, the singlet sector where the total angular momentum around the vertex is zero, and the triplet sector where the total angular momentum is equal to one.

### a.3 Reduction to qubits - Triangular lattice

In the main text, we have considered the local four-dimensional Hilbert space spanned by the occupation of four different modes . By choosing the subspace with one excitation per link, we can faithfully represent its quantum state with two qubits, the position and spin qubits, . The excitations of these qubits are related to the left and right modes through and .

With these definitions, the left and right operators defined in the main text are recast into

 Unknown environment '% (17)

with , . The Hamiltonian for one plaquette can then be written down as

 HT=−JTr[U(→x,^μ)U(→x+^μ,^ν)U(→x+^μ+^ν,−^μ−^ν)]. (18)

In order to identify all the terms of the trace one can use the condition . The condition is met if and only if

 α=β=γ=0, α≠β≠γ≠0, α=0;β=γ≠0, β=0;α=γ≠0, γ=0;β=α≠0. (19)

The first condition gives term, the second terms, and terms for the last three equations, for a total of terms

 Tr[U12U23U31]=14{Γ012Γ023Γ031+∑abcϵabcΓa12Γb23Γc31−∑a[Γ012Γa23Γa31+Γa12Γ023Γa31+Γa12Γa23Γ031]}, (20)

where we have used the subscripts to represent the links , respectively. The pure-gauge Hamiltonian can then be written in terms of matrices as

 H=−J2{Γ012Γ023Γ031+∑abcϵabcΓa12Γb23Γc31−∑a[Γ012Γa23Γa31+Γa12Γ023Γa31+Γa12Γa23Γ031]}. (21)

Then, one can explicitly write this Hamiltonian in terms of up to six-body interactions between qubits, for instance

 Γ012Γ023Γ031=σxpos12σxpos23σxpos31Γa12Γb23Γc31=σypos12σypos23σypos31σam12σbm23σcm31∑aΓ012Γa23Γa31=σxpos12σypos23σypos31∑aσam23σam31. (22)

### a.4 Experimental setups

In this section we provide a brief overview of the experimental setups considered for the implementation of the proposed simulator. The first setup considered in the main text is comprised of six transmon qubits provided with tunable coupling to a single microwave resonator. A transmon qubit is a special instance of a Cooper pair box (66), operating at specific regimes.

A Cooper pair box is made of a piece of a superconducting island connected to a voltage source through a capacitor, and grounded via a Josephson junction. The temperature of these superconducting systems is low enough so that the quantum state of the setup can be defined by the number of Cooper pairs on the superconducting island. In order to add tunability and control to this quantum system, the single Josephson junction of the Cooper pair box is substituted with a SQUID loop, threaed with magnetic fluxes. These superconducting systems can be coupled to superconducting microwave resonators, which can be used to induce, via electromagnetic interactions, tunnelling of Cooper pairs in and out the superconducting island, as well as reading out the state of the system (67); (68). The Hamiltonian of the setup comprising a Cooper pair box coupled to a microwave resonator reads

 HS=4EC(n−ng)2−EJcosϕ+ωra†a+2βeVrmsn(a+a†). (23)

Here we have defined the quantized charge on the superconducting island, which is proportional to the number of Cooper pairs, the offset charge and the quantized flux that threads the SQUID loop. The () are bosonic operators that act upon the resonator field, of frequency . stands for the charging energy of the superconducting island, while is the Josephson energy of the dc-SQUID loop. The coupling coefficient is a function of the circuit capacitances, is the root mean square voltage of the resonator, and is the electron charge.

One of the main problems that arise in using Cooper pair boxes as controlled quantum systems is the high sensitivity of their quantum state to noise on the offset charge . By using large capacitive pads shunting the superconducting island, one can achieve insensitivity with respect to this noise. These Cooper pair boxes provided with large capacitive pads are known as transmon qubits (69). In this parameter regime, the Hamiltonian of Eq. (23) models a transmon qubit. It can be rewritten in a diagonal basis for the part which is not coupled to the resonator,

 HS→2∑i=0ωi|i⟩⟨i|+ωra†a+2∑i=0gi,i+1(|i⟩⟨i+1|+|i+1⟩⟨i|)(a+a†), (24)

where we have included the first three levels of this device. Higher levels are usually excluded from the quantum dynamics, as being more energetic and hard to access. The first three levels can be approximately modeled as the first three levels of an anharmonic oscillator, with a relative anharmonicity factor of . If is big enough, one can consider the transmon qubit device a two-level system (69).

The setups envisioned in this work use two variants of the transmon design. Each of the six devices sketched in Fig. 1a in the main text is made of three superconducting islands, connected via two SQUID loops, marked with different colors in the figure. The magnetic fluxes acting on these loops modify the quantum level structure of the system, and an appropriate manipulation of these fluxes in time allows for a fast modulation of the effective coupling to the resonator (see Eq. 24). This modulation does not affect the bare energy levels of the device (70); (71) and gives the possibility of performing red and blue sidebands upon the coupled qubit-resonator system, which, in turns, allows to generate many-body interactions (72).

The setup considered in Fig. 1b in the main text describes transmon qubits with improved design, the so called Xmon qubits (73). In this case, the shape of the capacitance pads is designed in order to improve coherence times of the device, by reducing coupling to material defects and radiative losses. The interactions among these qubits are mediated by direct capacitive couplings, not making use of a common superconducting resonator. For this reason, the many-body interactions are reduced in this case to nearest-neighbor two-qubit gates.

### a.5 Many-body interactions and gauge invariance for different sectors

In order to simulate the Hamiltonian associated with the non-Abelian model in Eq. (21), one can simulate stepwise the single many-body monomials that compose the total interaction. Each of them has the form of an -qubit string term

 HN=σi1σj2⋯σkN, (25)

where the set of superindices can assume values . Up to a maximum of local rotations, this general interaction can be mapped back to, e.g., . In what follows, we will show how to obtain this effective interaction in two different scenarios. In one case, one can use a set of tunable-coupling transmon qubits coupled to a single microwave resonator. In the second scenario, a set of charge-like qubits, e.g., transmon or Xmon qubits, are coupled to an additional ancillary qubit, for a total of qubits.

We consider the sequence of collective and single qubit gates (74); (72), as

 US(t)=ei(π/4)∑i

and obtain the set of interactions

 US(t)= e−igtσz1σx2⋯σxN,N=4n−1, US(t)= eigtσz1σx2⋯σxN,N=4n+1, US(t)= e−igtσy1σx2⋯σxN,N=4n, US(t)= eigtσy1σx2⋯σxN,N=4n−2. (27)

The collective interaction in Eq. (26) can be realized within a set of three-island transmon qubits, by threading the two SQUID loops embedded in each device with fast-tunable magnetic fluxes, tuned at appropriate frequencies (72). The general -body term of Eq. (25) can be simulated with this approach with a total of collective gates, and a number of single qubit gates which is upper bounded by . The Hamiltonian in Eq. (21) is composed of monomials, three-body term, six-body terms, and five-body terms. Therefore, the simulation for one digital step will amount to implement collective gates and a number of single-qubit gates which is bounded by .

One can implement the interaction in a set of capacitively-coupled charge-like qubits, by decomposing it into a set of -phase and single-qubit gates. First, let us consider the sequence of gates in Eq. (26) between qubits and an additional ancillary one, labeled by . We consider also that the gates are rotated about the axis of an angle . The sequence of gates in Eq. (26) can then be rewritten in the following way

 US,A(t)=ei(π/4)σzA∑Nj=1σzje−igtσxAe−i(π/4)σzA∑Nj=1σzj. (28)

This sequence of gates is equivalent to the following interactions

 US,A(t)= eigtσxAσz1⋯σzN,N=4n−1, US,A(t)= e−igtσxAσz1⋯σzN,N=4n+1, US,A(t)= e−igtσyAσz1⋯σzN,N=4n, US,A(t)= eigtσyAσz1⋯σzN,N=4n−2. (29)

If the ancillary qubit is initialized in an eigenstate of (), (), then the effective dynamics on the -qubit system would be equivalent to the interaction . The interaction can be realized as a sequence of gates between the ancillary qubit and all the other ones, . Each of these two-body gates can be realized as a sequence of two local -rotations and a -phase gate between the qubits and , . In fact, defining

 eiϕ/2e−i(ϕ/2)σz=(100eiϕ),CP1i(ϕ)=⎛⎜ ⎜ ⎜⎝100001000010000e−i2ϕ⎞⎟ ⎟ ⎟⎠, (30)

one has that

 e−i(ϕ/2)σz1e−i(ϕ/2)σziCP1i(ϕ)=ei(ϕ/2)σz1σzi. (31)

In this way, an arbitrary -body term can be encoded, up to a number of local rotations bounded by , with a sequence of -Phase gates and a number of single-qubit gates bounded by . Thus the quantum simulation of one digital step of the Hamiltonian in Eq. (21) will amount to realize -Phase gates and a number of single-qubit rotations which is bounded by .

The -Phase gate between two qubits can be performed by adiabatically changing the frequency of one qubit, bringing the level close to the anharmonic excitated level , to induce state-dependent phase accumulations. Unwanted errors due to leakage to the state are exponentially suppressed in the gate duration. For finite gate times errors can be minimized upon optimization in the Fourier spectrum of the control signal (75). This optimization has been carried out for gates in (76). For a generic quantum simulation involving -Phase gates, one should define the set of angles to be used in the experiment, and optimize over the duration and Fourier spectrum of the gate, to avoid unwanted leakage. This has been done in (77) by increasing the gate time from ns to ns, not observing significant state leakages up to phases of rads. Additional errors on the fidelity of the gate come from standard decoherence of the qubits, and control errors, see Ref. (76) for a discussion of the error budget.

Finally, to complete the information given in the main text, we show in this section, in Fig. 6, additional plots that complement Fig. 3 in the main text. The simulated gauge invariance is shown to converge for the singlet and triplet configurations on a single triangular plaquette.

The proposed experiment can be extended to larger simulated lattices, which are implementable by considering setups of superconducting qubits with enhanced connectivity. The required connectivity can be achieved by using several microwave resonators, acting as non-local quantum buses, or by considering superconducting pads that can provide with multiple capacitive couplings, as schematized in Fig. 7. For instance, we can envision, as shown in the main text, a ladder architecture that can support a singlet gauge invariant subspace, where non-Abelian string dynamics could be studied.

### a.6 Scaling of errors for larger lattices

It has been shown (78) that the scaling of the number of steps needed to achieve a digital error in a Suzuki-Lie-Trotter decomposition is almost linear in the simulated phase of the dynamics and sub-polinomial in , according to the formula

 Ns≤2m52k(m||H||t)1+1/2kϵ1/2k. (32)

In the formula is the number of terms that compose a generic Hamiltonian , and is the fractal degree of approximation (79), which we set to . For one plaquette one has , three-body term, six-body terms, and five-body terms, all present in the system with the energy . Considering that the norm of an arbitrary product of Pauli matrices is , , one has, using a triangular inequality, that the norm of the Hamiltonian for one plaquette is bounded by . As a consequence, the number of steps necessary to implement a dynamics with error for the time on a lattice of plaquettes is bounded by

 Ns≤2300N(N|J|t)3/2√ϵ. (33)

Notice that this formula gives a very loose bound which is not dependent on the underlying structure of the Hamiltonian. Nevertheless, it shows that the scaling of the number of steps is efficiently solvable by a quantum simulator, not being exponential in the resources used. The exact scaling of the number of steps cannot be known (up to an exact solution of the problem, which is not addressable numerically for larger lattices), however it is useful to extract the dependence from the upper-bound formula in Eq. (33). Assuming independent experimental errors on the single Trotter steps, one can relate the total number of steps to the required experimental fidelity for the implementation of a single Hamiltonian term. We impose an accumulated experimental error of the same order of magnitude of the Trotter error . This enforces the error on the single step to be of the order . Besides these issues on error scaling and accumulation, the accesible simulation time has to be smaller than and of the single qubits, defining a system of necessary conditions for the simulation of a lattice of plaquettes

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩Ns≤2300N(N|J|t)3/2√ϵϵs∼ϵNst≪T1,T2. (34)

### a.7 Matter-field Su(2) Gauge Interactions

In the main text we have focused on pure gauge dynamics with plaquette interactions, where usual perturbative approaches do not get energy scales which are big enough. In this section, we point out the basic elements needed for the simulation of the simplest gauge model coupled to a matter field in -dimensions. The basic ingredients and the validity of the different energy-scale regimes have been already discussed in the simulation of Abelian gauge models in trapped ions and superconducting circuits (80); (81). This section can be seen as a generalization to non-Abelian setups.

The ideal interaction Hamiltonian that we would like to simulate is given by

 Hint=∑ib†(i)U(i,i+1)b(i+1)+H.c.=∑iαβb†α(i){[cαL(i)c†βR(i+1)]+∑μν[σyαμc†μL(i)cνR(i+1)σyνβ]}bβ(i+1)+H.c., (35)

where and are hard-core bosonic modes, with , and .

In terms of these modes, the former Hamiltonian has two local (gauge) symmetries:

• gauge invariance around every site with generators: