Superconducting Quantum Simulator for Topological Order and the Toric Code

Topological order is now being established as a central criterion for characterizing and classifying ground states of condensed matter systems and complements categorizations based on symmetries. Fractional quantum Hall systems and quantum spin liquids are receiving substantial interest because of their intriguing quantum correlations, their exotic excitations and prospects for protecting stored quantum information against errors. Here we show that the Hamiltonian of the central model of this class of systems, the Toric Code, can be directly implemented as an analog quantum simulator in lattices of superconducting circuits. The four-body interactions, which lie at its heart, are in our concept realized via Superconducting Quantum Interference Devices (SQUIDs) that are driven by a suitably oscillating flux bias. All physical qubits and coupling SQUIDs can be individually controlled with high precision. Topologically ordered states can be prepared via an adiabatic ramp of the stabilizer interactions. Strings of qubit operators, including the stabilizers and correlations along non-contractible loops, can be read out via a capacitive coupling to read-out resonators. Moreover, the available single qubit operations allow to create and propagate elementary excitations of the Toric Code and to verify their fractional statistics. The architecture we propose allows to implement a large variety of many-body interactions and thus provides a versatile analog quantum simulator for topological order and lattice gauge theories.

Topological order is now being established as a central criterion for characterizing and classifying ground states of condensed matter systems and complements categorizations based on symmetries. Fractional quantum Hall systems and quantum spin liquids are receiving substantial interest because of their intriguing quantum correlations, their exotic excitations and prospects for protecting stored quantum information against errors. Here we show that the Hamiltonian of the central model of this class of systems, the Toric Code, can be directly implemented as an analog quantum simulator in lattices of superconducting circuits. The four-body interactions, which lie at its heart, are in our concept realized via Superconducting Quantum Interference Devices (SQUIDs) that are driven by a suitably oscillating flux bias. All physical qubits and coupling SQUIDs can be individually controlled with high precision. Topologically ordered states can be prepared via an adiabatic ramp of the stabilizer interactions. Strings of qubit operators, including the stabilizers and correlations along non-contractible loops, can be read out via a capacitive coupling to read-out resonators. Moreover, the available single qubit operations allow to create and propagate elementary excitations of the Toric Code and to verify their fractional statistics. The architecture we propose allows to implement a large variety of many-body interactions and thus provides a versatile analog quantum simulator for topological order and lattice gauge theories.
Topological phases of quantum matter [1, 2] are of scientific interest for their intriguing quantum correlation properties, exotic excitations with fractional and non-Abelian quantum statistics, and their prospects for physically protecting quantum information against errors [3].
The model which takes the central role in the discussion of topological order is the Toric Code [4][5][6], which is an example of a Z 2 lattice gauge theory. It consists of a two-dimensional spin lattice with quasi local fourbody interactions between the spins and features 4 g degenerate, topologically ordered ground states for a lattice on a surface of genus g. The topological order of these 4 g ground states shows up in their correlation properties. The topologically ordered states are locally indistinguishable (the reduced density matrices of a single spin are identical for all of them) and only show differences for global properties (correlations along non-contractible loops). Moreover local perturbations cannot transform these states into one another and can therefore only excite higher energy states that are separated by a finite energy gap.
In equilibrium, the Toric Code ground states are therefore protected against local perturbations that are small compared to the gap, which renders them ideal candidates for storing quantum information [4]. Yet for a selfcorrecting quantum memory, the mobility of excitations needs to be suppressed as well, which so far has only been shown to be achievable in four or more lattice dimensions [7,8]. Excitations above the ground states of the Toric Code are generated by any string of spin rotations with open boundaries. The study of these excitations, called anyons, is of great interest as the neither show bosonic nor fermionic but fractional statistics [4,5].
In two dimensions, the Toric Code can be visualized as a square lattice where the spin-1/2 degrees of freedom are placed on the edges, see Fig. 1a. Stabilizer operators can be defined for this lattice as products of σ x -operators for all four spins around a star (orange diamond in Fig. 1a), i.e. A s = j∈star(s) σ x j , and products of σ y -operators for all four spins around a plaquette (green diamond in Fig. 1a), i.e. B p = j∈plaq(p) σ y j . The Hamiltonian for this system reads, where J s (J p ) is the strength of the star (plaquette) interactions and the sums s ( p ) run over all stars (plaquettes) in the lattice [9]. Despite the intriguing perspectives for quantum computation and for exploring strongly correlated condensed matter physics that a realization of the Toric Code Hamiltonian would offer, progress towards its implementation has to date been inhibited by the difficulties in realising clean and strong four-body interactions while at the same time suppressing two-body interactions to a sufficient extent. Realisations of topologically encoded qubits have been reported in superconducting qubits [10,11], photons [12] and trapped ions [13]. The Toric Code is also related to a model of geometrically frustrated spins on a honeycomb lattice [14]. Yet, in contrast to the Toric Code itself, this Honeycomb model only features two-body interactions, which can naturally occur in spin lattices, and evidence for its realization in iridates [15] and ruthenium-based materials [16] has been found in neutron scattering experiments.
There is moreover an important difference in the use arXiv:1608.04565v2 [quant-ph] 17 Mar 2017 and purpose of an implementation of the Hamiltonian (1) compared to error detection via stabilizer measurements as pursued in [12,[17][18][19][20]. In particular, two-body interactions are all that is needed for arbitrary stabiliser measurements, while lattice gauge theories [21][22][23][24] and Hamiltonians such as (1) require k-body (with k > 2) effective interactions. We note that the realisation of k-body interactions is also important in quantum annealing architectures, where implementations have been proposed utilising local interactions and many-body constraints [25][26][27]. We anticipate that the approach we describe here may also be relevant for quantum annealing, however our focus in this paper is on the realisation of Hamiltonians supporting topological order.
Superconducting circuits have recently made tremendous advances in realizing engineered quantum dynamics for quantum simulation [28][29][30][31][32] and quantum information processing [17][18][19][33][34][35][36][37]. Their quantum features originate from integrated Josephson junctions that behave like ideal nonlinear inductors. They are typically built as SQUIDs that can be tuned and driven via magnetic fluxes. Currently superconducting circuits are assembled in ever larger connected networks [17,18,[34][35][36]38]. Here we make use of these exquisite properties of superconducting circuits for a scheme implementing the Hamiltonian of the two-dimensional Toric Code in a direct way by creating four-body interactions of the form A s = j∈star(s) σ x j and B p = j∈plaq(p) σ y j between the respective qubits and fully suppressing all two-body interactions.
Superconducting circuits are an ideal platform for experimentally exploring topological order since high precision control and measurements of both, individual qubits as well as multi-qubit clusters have been demonstrated [17,30,[39][40][41] . Importantly, this platform is also suitable for realizing periodic boundary conditions, since conductors can cross each other via air-bridges [42,43]. This property, which is essential for exploring topological order, is so far not accessible in other platforms unless one resorts to digitized approaches.
Our approach is designed to keep the necessary circuitry as simple as possible and only employs standard superconducting qubits that are coupled via dc-SQUIDs [44,45]. To generate the desired interactions we make use of the nonlinear nature of the coupling SQUIDs and the possibilities to bias them with constant and oscillating magnetic fluxes. In particular the four-body stabilizer interactions of the Toric Code are generated by driving the SQUIDs with suitably oscillating magnetic fluxes. They can therefore be switched on and off at desired times or even be smoothly tuned in strength. In this way the Toric Code Hamiltonian is realized in a frame where the qubits rotate at their transition frequencies. Importantly, the transformation between the rotating and the laboratory frame is composed of single qubit unitaries only, so that the entanglement properties of the states are identical in both frames. Yet for implementing the Toric Code Hamiltonian (1), considering a rotating frame has the pivotal advantages that the effective transition frequencies of the individual spins become arbitrarily small and that the approach becomes largely insensitive to disorder in the lattice. In contrast to earlier approaches to multibody interactions based on gate sequences [46][47][48][49], the four-body interactions in our approach are directly imple-mented as an analog quantum simulator and are therefore time-independent, which avoids any errors routed in the Trotter sequencing. Moreover, our approach suppresses two-body interactions to an extent that they are negligible compared to the four-body interactions.
A minimal two-dimensional lattice covering the surface of a torus requires eight qubits connected via superconducting wires such that the lattice has periodic boundary conditions in both directions and is thus realizable with available technology [17,18,34]. We show that the topologically ordered states of the Toric Code can be prepared via an adiabatic sweep from a product state for experimentally realizable parameters and that measurements of individual qubits and qubit-qubit correlations can be performed to demonstrate the topological order of these states. Moreover the robustness of the topologically ordered states against local perturbation that would change the energy of the system by a limited amount can be verified experimentally.

I. THE PHYSICAL LATTICE
For realizing the Toric Code lattice, star and plaquette interactions, indicated by orange and green diamonds in Fig. 1a, are both realized via the circuit shown in Fig. 1b but with different external fluxes Φ ext applied to the coupling SQUID. The lattice resulting from the orange and green diamonds in Fig. 1a with spins on the lattice vertices rather than edges is a convenient representation introduced in [50], see also [51,52]. The small black squares in Figs. 1a and 1b indicate spins formed by the two lowest energy eigenstates of superconducting qubits. For the transition frequencies between these two lowest states, we assume a periodic pattern in both lattice directions such that no nearest neighbors and no next nearest neighbors of qubits have the same transition frequency and hence all qubits in each four-spin cell have different transition frequencies, see Sec. III for details.
For each cell the four qubits on its corners are connected by a dc-SQUID formed by two identical Josephson junctions with Josephson energy E J , see Fig 1b (an asymmetry between the junctions causes a negligible perturbation only). The dc-SQUID is threaded by a time dependent external flux, Φ ext and its ports are connected to the four qubits via four inductors, each with inductance L. There are thus auxiliary nodes at both ports of the SQUID. As we show in Sec. II, these auxiliary nodes (a; n, m) and (b; n, m) can be eliminated from the description while, together with the flux bias Φ ext on the SQUID, they mediate the desired four-body interactions of the stabilizers.

A. Hamiltonian of the physical lattice
We describe the lattice in terms of the phases ϕ j and their conjugate momenta π j for each node, which obey the commutation relations [ϕ j , π l ] = i δ j,l . The nodes are labelled by composite indices j = (n, m), where n labels the diagonals extending from top left to bottom right and m labels the qubits along each diagonal in Fig. 1a.
Each qubit is characterized by a Josephson energy E Jq;j and a charging energy E Cq;j = e 2 /(2C q;j ), where e is the elementary charge and C q;j the qubit's capacitance. Including a contribution 4E L ϕ 2 j that arises from the coupling to its four adjacent auxiliary nodes via the inductances L, the Hamiltonian of the qubit at node j reads, where E L = φ 2 0 /(2L) and φ 0 = /(2e) is the rescaled quantum of flux. We consider transmon qubits [17,18,29,30,34,35,53] because of their favorable coherence properties and moderate charging energies, c.f. appendix D 1.
For the description of the dc-SQUIDs it is useful to consider the modes ϕ ±;j = ϕ a;j ± ϕ b;j , where ϕ a/b;j is the phase at the auxiliary node (a/b; n, m). The modes ϕ +;j behave as harmonic oscillators described by the Hamiltonians H +;j = 4E C+ −2 π 2 +;j + E L ϕ 2 +;j , where π ±;j are the canonically conjugate momenta to ϕ ±;j and E C+ = e 2 /C g , with C g the capacitance between each of the nodes (a; j) or (b; j) and ground. The modes ϕ −;j in turn are influenced by the dc-SQUID. Their Hamiltonians read (3) where ϕ ext = Φ ext /φ 0 is the phase associated to the external flux through the SQUID loop, E J the Josephson energy of each junction in the SQUID and E C− = e 2 /2C J with C J the capacitance of one junction in the SQUID, including a shunt capacitance parallel to the junction.
We choose the indices j = (n, m) for the auxiliary nodes such that node (a; n, m) couples to the qubits at nodes (n, m) and (n, m + 1) whereas node (b; n, m) couples to the qubits at nodes (n + 1, m) and (n + 1, m + 1), see Fig. 1b. The interaction terms between the qubit variables ϕ j and the SQUID modes ϕ ±;j due to the inductances L thus read, where ϕ q±;j = ϕ n,m +ϕ n,m+1 ±ϕ n+1,m ±ϕ n+1,m+1 . Here, inductive couplings via mutual inductances are an alternative option [54]. The full Hamiltonian of the lattice is therefore given by A detailed derivation of the Hamiltonian H, starting from a Lagrangian for the circuit, is provided in appendix A. We now turn to explain how the four-body interactions of the stabilizers emerge in our architecture.

II. STABILIZER OPERATORS
The guiding idea for our approach comes from the following approximate consideration (We present a quantitatively precise derivation in Sec. II A). In each single cell, as sketched in Fig. 1b, the auxiliary nodes (a; n, m) and (b; n, m) have very small capacitances with respect to ground and therefore cannot accumulate charge. As a consequence the node phases ϕ a;n,m and ϕ b;n,m are not dynamical degrees of freedom but instantly follow the phases at nodes (n, m) and (n, m + 1) or (n + 1, m) and (n + 1, m + 1) as given by Kirchoff's current law saying that the sum of all currents into and out of such a node must vanish. Therefore the Josephson energy of the coupling SQUID, 2E J cos(ϕ ext /2) cos(ϕ a;n,m − ϕ b;n,m ), becomes a nonlinear function of the node phases ϕ n,m , ϕ n,m+1 , ϕ n+1,m and ϕ n+1,m+1 , which contains four-body terms.
The external flux applied through the SQUID is then chosen such that it triggers the desired four-body interactions forming the stabilizers. To this end it is composed of a constant part with associated phase ϕ dc and an oscillating part associated with phase amplitude ϕ ac , where the oscillation frequencies of F (t) are linear combinations of the qubit transition frequencies such that the desired four-body spin interactions of the Hamiltonian (1) are enabled and all other interaction terms are suppressed. Moreover we choose ϕ dc = π to suppress the critical current I c = 2E J φ −1 0 cos(ϕ dc /2) of the dc-SQUID [55] and thus unwanted interactions of excitations between qubits at opposite sides of it, c.f. Fig. 1b. Assuming furthermore that the oscillating part of the flux bias is small compared to the constant part, we approximate cos(ϕ ext /2) ≈ −ϕ ac F (t)/2.
The above consideration uses approximations and only serves us as a guide. To obtain a quantitatively precise picture, we here take the small but finite capacitances of the Josephson junctions in the SQUID into account and derive the stabilizer interactions in a fully quantum mechanical approach. Our derivation makes use of methods for generating effective four-body Hamiltonians in the low energy sector of a two-body Hamiltonian [14,[56][57][58] by employing a Schrieffer-Wolff transformation [59] and combines these with multi-tone driving in order to single out specific many-body interactions while fully suppressing all two-body interactions.
Further intuition for our approach can be obtained from the macroscopic analogy explained in Fig. 2. In short, the plasma frequency of the coupling SQUIDs differs form all transition frequencies of the qubits they couple to as well as from all frequencies contained in the oscillating flux that is threaded through their loops. They thus decouple from the qubits but mediate the desired four-body interactions via the drives applied to them. Further details of this analogy are explained in the caption of Fig. 2. We now turn to present the quantitatively precise derivation. The qubits are nonlinear oscillators so that one can think of them as masses (corresponding to the qubits' capacitances) attached to anharmonic springs (corresponding to the qubits Josephson energy). We here sketch four nonlinear oscillators representing four qubits as green spheres connected by black (anharmonic) springs to the support frame. Note that all four qubits have different transition frequencies so that all four spheres representing them here would have different masses and their springs different spring constants. For vanishing critical current the relevant SQUID modes ϕ−;j, behave as harmonic oscillators and in our macroscopic analogy we can represent them by a mass (here sketched as a red sphere, corresponding to the capacitance 2CJ ) that is attached via four springs (sketched in orange, corresponding to the inductors L) to the four masses representing the neighboring qubits. For vanishing critical current, the only potential energy of the modes ϕ−;j comes from their inductive coupling to the neighboring qubits, and hence the mass representing the SQUID mode (red sphere) has no spring connecting it to a support frame. We now imagine that we shake the mass representing the SQUID mode with a driving term that is proportional to the fourth power of the mass' position variable, ϕ 4 −;j . This corresponds to the oscillating flux bias ∝ ϕacF (t). If we shake this mass at a frequency ω d that is distinct from its eigenfrequency ω−;j but equals the sum of the oscillation frequencies of all four neighboring masses (neighboring qubits), we drive a term proportional to q † 1 q † 2 q † 3 q † 4 , where q † j creates an excitation in qubit j. If, on the other hand we shake the SQUID-mass at a frequency ω d = ω1 + ω2 − ω3 − ω4 we drive a term proportional to q † 1 q † 2 q3q4. The selective driving of specific four-body terms, as we consider it here, of course relies on choosing the transition frequencies of all four qubits to be different. Our approach considers a superposition of various drive frequencies such that the sum of the triggered four-body terms adds up to the desired stabilizer interaction, see text.
A. Elimination of the SQUID degrees of freedom For the parameters of our architecture, the oscillation frequencies of the SQUID modes ϕ ±;j differ from the transition frequencies of the qubits and the frequencies of the time dependent fluxes threaded through the SQUID loops by an amount that exceeds their mutual couplings by more than an order of magnitude. These modes will thus remain in their ground states to very high precision. To eliminate them from the description and obtain an effective Hamiltonian for the qubits only, it is useful to consider the Schrieffer-Wolff transformation which eliminates the qubit-SQUID interactions described in Eq. (4). The effective Hamiltonian for the qubits is then obtained by projecting the resulting Hamiltonian onto the ground states of the SQUID modes ϕ ±;j . The details of this derivation are discussed in appendix B. For the transformed HamiltonianH, see Eq. (7), we consider terms up to fourth order in S to get the leading contributions of the terms which are in resonance with the oscillating part of ϕ ext (t), c.f. Eq. (6) . On the other hand, a truncation of the expansion is here well justified since both fields ϕ q±;j and π ±;j have low enough amplitudes, see Eq. (B16), so that higher order terms can safely be ignored. The transformed Hamiltonian can be written asH where H 0q = (N,M ) j=(1,1) H q;j with H q;j as in Eq.
(2) describes the individual qubits and interactions between neighboring qubits that are ineffective since their transition frequencies differ by an amount that exceeds the relevant coupling strength by more than an order of magnitude, see appendix D 1. The remaining two terms H S and H qS in Eq. (8) describe the SQUIDs and the residual couplings between the qubits and the SQUIDs, see Eqs. (B7) and (B8) in appendix B. Since the oscillation frequencies of the SQUID modes are chosen such that they cannot be excited by the applied driving fields or via their coupling to the qubits, an effective Hamiltonian, describing only the qubits, can be obtained by projecting H qS onto the ground states of the SQUID modes, As the SQUID modes remain in their ground states, their Hamiltonian H S may be discarded. The projected Hamiltonian H (0) qS contains two-body interactions and four-body interactions between the qubits, with χ = 0 −;n,m | cos(ϕ −;n,m )|0 −;n,m , where |0 −;n,m is the ground state of the mode ϕ −;n,m . We emphasize here that our approach requires nonlinear coupling circuits, such as the considered SQUIDs. The terms in the second line of Eq. (11), which give rise to the desired four-body interactions, only emerge because the Hamiltonian of the coupling SQUIDs contains terms that are at least fourth power in ϕ −;n,m , see appendix B.
The oscillation frequencies contained in the external flux, i.e. in ϕ ext , can be chosen such that selected terms in the second line of Eq. (11) are enabled as they become resonant. For suitable amplitudes of the frequency components of ϕ ext , these terms add up to the desired fourbody stabilizer interactions whereas all other interactions contained in Eq. (11) are for our choice of parameters strongly suppressed in a rotating wave approximation.
Before discussing this point in more detail in the next section, let us here emphasize the versatility of our approach: With suitably chosen Fourier components of ϕ ext almost any two-, three-or four-body interaction can be generated. The only limitation here is that terms with the same rotation frequency can not be distinguished, e.g. one can not choose between generating σ z j σ z k or σ z j σ z l for j = k = l which are both non-rotating. On the other hand, straightforward modifications of the circuit allow to implement higher than four-body (e.g. five-and sixbody) interactions.

B. Effective Hamiltonian of one cell
We now focus on the interactions between the four qubits in one cell only, see Fig. 1b. To simplify the notation we label these four qubits clockwise with indices 1,2,3 and 4, i.e. (n, m) → 1, (n + 1, m) → 2, (n + 1, m + 1) → 3 and (n, m + 1) → 4. We thus divide one term, say for (n, m) = (n 0 , m 0 ), in the sum of Eq. (11) into the following parts, where contains the four-body interactions of the stabilizers, contains all terms that are quadratic in the node phases ϕ j and contains all terms that are fourth-order in the node phases ϕ j except for the four-body interactions. These are terms where at least two of the indices i, j, k and l are the same as indicated by the notation (i,j,k,l) . Details of the derivation of the Hamiltonian (12) are given in appendix C, where we also state the exact form of the coefficients C ijkl . The Hamiltonian (12) can be written in second quantized form using ϕ j = ϕ j (q j + q † j ), where q † j (q j ) creates (annihiliates) a bosonic excitation at node j with energy ω j = 8E C;j (E Jq;j + 4E L ) and zero point fluctuation amplitude ϕ j = [2E C;j /(E Jq;j + 4E L )] 1/4 . Note that the effective Josephson energy E Jq;j +4E L includes here contributions of the attached inductances, see appendix B 1 for details. Due to the qubit nonlinearities we can restrict the description to the two lowest energy levels of each node by replacing q j → σ − j (q † j → σ + j ) and get for the four-body interactions, Here, the sum a,b,c,d∈{+,−} runs over all 16 possible strings of σ + and σ − operators, which are listed in table II in appendix C.
The explicit form of the Hamiltonians H 2 in second quantized form is given in Eq. (D15). The terms contained in H 2 only lead to frequency shifts for the qubits and negligible two-body interactions. The frequency shifts can be compensated for by modified frequencies for the driving fields and the residual two-body interactions can be made two orders of magnitude smaller than the targeted four-body interactions, see Sec. II C and appendix D. Moreover the terms contained in H 3 only lead to corrections that are at least an order of magnitude smaller than those of H 2 and can safely be discarded, see appendix D 3. Here, we therefore focus on the part H 1 which gives rise to the dominant terms in the form of four-body stabilizer interactions.

C. Triggering four-body interactions
In a rotating frame where each qubit rotates at its transition frequency, σ + , the strings σ a 1 σ b 2 σ c 3 σ d 4 that appear in Eq. (13) rotate at the frequencies ω a,b,c,d = aω 1 + bω 2 + cω 3 + dω 4 . Hence, if we thread a magnetic flux, which oscillates at the frequency ω a,b,c,d , through the SQUID loop, the corresponding operator σ a 1 σ b 2 σ c 3 σ d 4 acquires a pre-factor with oscillation frequency ω a,b,c,d , which renders the term non-oscillating and hence enables this particular four-body interaction. Importantly, all qubits of a cell have different transition frequencies such that a drive with frequency ω a,b,c,d only enables the particular interaction σ a 1 σ b 2 σ c 3 σ d 4 whereas all other interactions are suppressed by their fast oscillating prefactors.
As the star and plaquette interactions J s σ x 1 σ x 2 σ x 3 σ x 4 and J p σ y 1 σ y 2 σ y 3 σ y 4 can be written as linear combinations of operator products σ a 1 σ b 2 σ c 3 σ d 4 , their generation requires magnetic fluxes that are a superposition of oscillations at the respective frequencies. To enable the specific four-body interactions of the Toric Code we thus consider, where The driving field F s generates a star interaction A s and F p a plaquette interaction B p . Applying F s or F p to a cell, its Hamiltonian (13) thus reads, where r.t. refers to the remaining fast rotating terms, and Hence star and plaquette interactions are realized with the same coupling circuit (shown in figure 1b) and only differ by the relative sign of the F X and F Y contributions, i.e. a relative phase of π in the respective Fourier components. Consequently, a driving field with F s (F p ) is applied to the coupling circuit of every orange (green) diamond in figure 1a (see also figure 3) and the resulting checkerboard pattern of driving fields generates the Toric Code Hamiltonian of equation (1). Note that the signs and magnitudes of the interactions J s and J p can for each cell be tuned independently via the respective choices for the drive amplitudes ϕ ac .
Besides the stabilizer interactions, there are further terms in the Hamiltonian (8). As we show in appendix D, all these other interactions are strongly suppressed by their fast oscillating prefactors (the oscillation frequencies of the prefactors exceed the interaction strengths by more than an order of magnitude in all cases). Their effect can be estimated using time dependent perturbation theory. We find for our parameters that they lead to frequency shifts for the individual qubits, which can be compensated for by a modification of the frequencies of the driving field ϕ ext . For our choice of qubit transition frequencies, where no nearest neighbors or next nearest neighbors of qubits have the same transition frequencies, the leading corrections to the Toric Code Hamiltonian of Eq. (1), besides local frequency shifts, are two-body interactions between qubits that are further apart (i.e. neither nearest neighbors nor next neighrest neighbors). For suitable parameters as discussed in Sec. III, these interactions are about two orders of magnitude weaker than the stabilizers and thus even weaker than dissipation processes in the device, see appendix D.
In this context, we also note that the leading contribution to the perturbations of the Hamiltonian (1) is proportional to the charging energy E Cq;j of the employed qubits, c.f. appendix D. Hence the considered transmon qubits are a well suited choice for our aims.

D. Including dissipative processes
As we have shown, the dynamics of the architecture we envision can be described by the Hamiltonian (1) in a rotating frame with respect to H 0 = j ω j σ + j σ − j , where we assume that all frequency shifts discussed in the previous section (see also appendix D) have been absorbed into a redefinition of the ω j . In an experimentally realistic scenario, the qubits will however be affected by dissipation. The full dynamics of the circuit can therefore be described by the Markovian master equation, where ρ is the state of the qubit excitations, H is as in In the following we turn to discuss realistic experimental parameters that lead to an effective Toric Code Hamiltonian, the preparation and verification of topological order in experiments and protocols for probing the anyonic statistics of excitations in a realization of the Toric Code with our approach. In all these discussions, we fully take the relaxation and dephasing processes described by Eq. (17) into account.

III. TOWARDS AN EXPERIMENTAL REALISATION
An experimental implementation of our approach is feasible with existing technology as it is based on transmon qubits which are the current state of the art in many laboratories. Moreover, periodic driving of a coupling element at frequencies comparable to the qubit transition frequencies is a well established technique [41,60]. The minimal lattice of a Toric Code with periodic boundary conditions in both lattice directions is the eight-qubit lattice sketched in Fig. 3. An important benefit of our rotating frame approach is its inherent insensitivity to disorder due to largely spaced qubit transition frequencies. The parameters we discuss in the sequel do thus Lattice for minimal realization on a torus with 8 physical qubits including 4 star interactions (orange diamonds) and 4 plaquette interactions (green diamonds). The qubits are indicated by small colored squares, where the coloring encodes their transition frequencies and the numbering indicates the periodic boundary conditions. Here, 8 Fig. 3. These frequencies are assumed to already include the frequency shifts calculated in appendix D.
not need to be realized exactly but only with moderate precision.

A. Parameters
A working example for this eight-qubit realization with experimentally realistic parameters is provided by the qubit transition frequencies given in table I, which are realized in qubits with Josephson energies in the range 4.9 GHz ≤ E Jq /h ≤ 94 GHz (h is Planck's constant) and capacitances in the range 52 fF ≤ C q;j ≤ 77 fF. Note that the frequencies are chosen such that neighboring qubits are detuned from each other by several GHz whereas the transition frequencies of next nearest neighbor qubits differ by ∼ 100 MHz. A practical advantage of our approach with respect to frequency crowding is that it only requires a small set of continuous single frequency driving fields, whereas the frequency width of the pulses in digitized schemes grows inversely proportional to their durations.
For the coupling SQUIDs, the Josephson energy of each junction should be E J /h ≈ 10 GHz and the associated capacitance 3.47 fF. Since high transparency Josephson junctions typically have capacitances of 1 fF per 50 GHz of E J , each SQUID is shunted by a capacitance of a few fF. For these parameters the modes ϕ −;j oscillate at ω − /2π = 8.45 GHz, whereas the oscillation frequency of the modes ϕ +;j is much higher since C g C J . Furthermore the inductances coupling the qubits to the SQUIDs need to be chosen large enough so that the resulting coupling term remains weak compared to the qubit nonlinearities. This requires inductances of about 100 nH which can be built as superinductors in the form of Josephson junction arrays [61,62] or with high kinetic inductance superconducting nanowires [63]. We note that our concept is compatible with nonlinearities Josephson junction arrays may have. The amplitude of the drive is ϕ ac = 0.1 rad corresponding to an oscillating flux bias of amplitude 0.1 × Φ 0 at the SQUID. These parameters give rise to a strength for star and plaquette interactions of J s = J p = h × 1 MHz, which exceeds currently achievable dissipation rates of transmon qubits significantly [32,[64][65][66][67]. The strength of the stabilizer interactions can be enhanced further by increasing the Josephson energies of the SQUIDs E J and/or the amplitudes of the oscillating drives ϕ ac .
For the above parameters, the total oscillation frequency, i.e. the sum of the frequencies of the qubit operators and the frequency of the external flux, is for any remaining two-body interaction term at least 10-30 times larger than the corresponding coupling strength.
To leading order, these terms give rise to frequency shifts of the qubits which can be compensated for by a modification of the drive frequencies. All interaction terms contained in higher order correction are negligible compared to the four-body interactions forming the stabilizers. We note that the experimental verification of the topological phase is further facilitated by its robustness against local perturbations, such as residual two-body interactions or shifts of the qubit transition frequencies, c.f. [68,69].
As an alternative to our approach, one could aim at implementing the Honeycomb model [14,70], which allows to obtain the Toric Code via a 4th order perturbative expansion in the limit where one of the spin-spin couplings strongly exceeds the two others [5]. This strategy would however require twice the number of qubits as our approach. For typical qubit-qubit couplings of ∼100 MHz, it would moreover only lead to stabilizer interactions that are two orders of magnitude weaker than in our approach.

B. Control and read-out
In common with all two-dimensional lattices, the architecture we envision requires an increasing number of control lines for the qubits, SQUIDs and the read-out as one enhances the lattice size. Whereas solutions to this complexity challenge are already being developed, e.g. by placing the control lines on a second layer [35,71,72], a minimal realization of the Hamiltonian (1) requires only eight qubits that form four stars and four plaquettes, see Fig. 3. Moreover, these eight physical qubits are arranged on a 4 × 2 lattice so that every qubit is located on its boundary, facilitating the control access.
Controlling single qubit states and two qubit correlation measurements would require a charge line and a readout resonator for each qubit. Our scheme, where a single cell features four qubits with different transition frequencies, has the beneficial property that four readout resonators can be connected to a single transmission line and the qubit states read out using frequency division multiplexing [73][74][75]. Such a joint readout is capable of measuring the state of all four qubits and qubit pairs even when their transition frequencies are separated by several GHz. An improvement over the readout through a common transmission line would be to employ a common Purcell filter [76,77] with two-modes, one mode placed at 3.5 GHz and the second at 12.5 GHz. Such scheme would not only allow a joint readout of four qubits but also protect them against the Purcell decay. Frequency division multiplexing can also be used for qubit excitation where a single charge line can be split or routed to four qubits in the cell. If a common transmission line is used for the readout the same line can also be used as a feedline for qubit state preparation [18].

C. Scalability
Importantly, our scheme is scalable as it also allows to realize the Toric Code on larger lattices. As already mentioned in sections I and II C, the requirement for the choice of the qubit transition frequencies is that no nearest neighbor qubits and no next nearest neighbor qubits have the same transition frequencies. The resulting pattern of transition frequencies for a large lattice is sketched in Fig. 4. It involves 16 qubit transition frequencies and can be realized with circuit parameters similar to the ones given above. These 16 qubit transition frequencies are periodically repeated in both lattice directions and therefore, in principle, sufficient for building any lattice size.
Larger and more realistic Toric Code realizations will obviously suffer from increased complexity. This could be drastically simplified when adding on-chip RF switches [78] for routing microwave signals to different cells and using selective broadcasting technique which would also reduce the number of microwave devices such as signal generators and arbitrary waveform generators [79].

IV. ADIABATIC PREPARATION OF TOPOLOGICAL ORDER
On the surface of a torus, i.e. for a two dimensional lattice with periodic boundary conditions in both lattice directions, the Hamiltonian (1) has four topologically ordered, degenerate eigenstates |ψ a (a = 1, 2, 3, 4), which are joint eigenstates of the stabilizers A s and B p with Here four principal transition frequencies are needed. These differ by a few GHz and are indicated by the filling colors brown, blue, red and purple. To suppress two-body interactions between next nearest neighbor qubits, each pair of these needs to be detuned from each other by ∼ 100 MHz. The resulting frequency differences are indicated by the colored frames around the respective squares. The gray diamond highlights the frequency set which is periodically repeated in both lattice directions. eigenvalue +1, see appendix E for a possible choice of the |ψ a for an eight-qubit lattice. A state in this manifold of topologically ordered states can either be prepared via measurements of the stabilizer operators or via an adiabatic sweep that starts from a product state [80]. Measurements of stabilizers can in our architecture be performed via a joint dispersive readout of all qubits in a star or plaquette, see Sec. V.
For the adiabatic sweep, one starts in our realisation with the oscillating driving fields turned off, ϕ ac = 0 and the qubit transition frequencies blue-detuned by a detuning ∆ with respect to their final values ω j . In the cryogenic environment of an experiment, this initializes the system in a product state, where all qubits are in their ground state, |ψ 0 = j |0 j . Then the amplitudes of the external oscillating fluxes, ϕ ac , are gradually increased from zero to their final values ϕ ac,f , i.e. ϕ ac (t) = λ(t) ϕ ac,f , where λ(t = 0) = 0, λ(T S ) = 1, and T S is the duration of the sweep. The qubit transition frequencies are decreased to ω j , i.e. ω j (t) − ω j = ∆[1 − λ(t)] > 0, with ω j (0) − ω j = ∆ and ω j (T S ) − ω j = 0. For this sweep, the Hamiltonian of the circuit can be written in the rotating frame as, Hamiltonian is H(λ = 0) = j ∆σ z j /2 and the unique initial ground state is the vacuum with no qubit excitations (all spins down), i.e. |ψ 0 = j |0 j . The final Hamiltonian is the targeted Toric Code Hamiltonian given in Eq. (1). Whereas the initial detuning ensures that the state |ψ 0 is also an eigenstate in the rotating frame for J s = J p = 0, the sweep will drive the system adiabatically into the manifold of protected states, provided the parameters are tuned sufficiently slow.
The spectrum of the minimal eight-qubit lattice, c.f. Fig. 3, is shown in Fig. 5 as a function of the tuning parameter λ. There is one unique ground state at λ = 0 which eventually becomes degenerate with three other states to form the degenerate four-dimensional manifold of topologically ordered ground states at λ = 1. In finite size lattices, there is however, for any value of λ, a finite energy separation of the order of J s or J p between the ground state and all states which do not end up in the topologically ordered manifold at λ = 1. We have numerically confirmed this for 8, 12 and 16 qubits. This spectral gap makes the sweep robust because it is not important which state in this manifold is prepared. All states in the manifold and any linear combination of them are equivalent. Therefore we use the probability of ending up in the topologically ordered manifold, F = 4 a=1 ψ a |ρ(T S )|ψ a , as a criterion for a successful preparation. Moreover, in the absence of dissipation one will end up in a pure state since a unitary evolution cannot change the purity of a state. To quantify how unavoidable dissipative processes in a real experiment would affect the preparation we also compute the purity of the final state P = Tr(ρ 2 ). Figure 6a shows the fidelity F for being in the topologically ordered qubit subspace as a function of the sweep time T S and the initial detuning ∆. Here ρ(t) is the actual state of the system in the rotating frame and the |ψ a (a = 1, . . . , 4) form the subspace of degenerate pro-tected ground states. To test that the prepared state is pure, we also plot the purity P as a function of T S and ∆ in Fig. 6b. We choose J s = J p = h×2 MHz, see Sec. II D, relaxation and dephasing times T 1 = T 2 = 50 µs, and a ramp profile λ(t) = (1 + ζ)(t/T S ) 6 /[1 + ζ(t/T S ) 6 ] with ζ = 100. For these values, a protected state can successfully be prepared.
For larger lattices one may also follow the scheme proposed in [80] which does not rely on a finite energy separation between the ground states and the higher excited states that do not become topologically ordered in the sweep. In our case, the scheme of [80] proceeds as follows. All qubits, which are initially in their states |0 j , are first transformed into states |+ j with σ x j |+ j = +1 |+ j via a π/2-pulse. Since the resulting state j |+ j is an eigenstate of all star stabilizers A s , these may then be switched on abruptly together with a term −B x j σ x j which stabilizes the state j |+ j . Subsequently the plaquette interactions are turned on and the term −B x j σ x j is turned off in a sweep that is adiabatic since the relevant transition matrix elements of the instantaneous Hamiltonian always vanish [80].
Once they are prepared, it is of course desirable to stabilize the topologically ordered states as long as possible. Several approaches to this task have been put forward, see [81,82] and references therein, which are compatible with the rotating frame implementation we consider. Most approaches consider a shadow lattice of auxiliary qubits that are coupled to the primary lattice. In our architecture, such auxiliary systems may not be required as the coupling SQUIDs could be employed for this task.

V. MEASURING TOPOLOGICAL CORRELATIONS
Our scheme allows to experimentally explore and verify topological order. To this end one prepares the lattice in one of the topologically ordered states |ψ a and measures the state of individual physical qubits as well as a correlation such as l∈v σ x l for a non-contractible path v along one of the lattice directions, which is closed in a loop to implement the periodic boundary conditions. For the minimal eight qubit realization this is merely a measurement of a two-spin correlation, i.e. σ x m σ x n . Noncontractible loops of σ y operations, such as l∈v σ y l , transform the topologically ordered states |ψ a into one another [5]. Here, these can be realized via strings of π-pulses applied to the qubits along the loop. After applying a string of π-pulses along one lattice direction, which transforms |ψ a into a different topologically correlated state |ψ b , one again measures individual spins as well as the correlation l∈v σ x l . Whereas the states of the individual spins should be the same for |ψ a and |ψ b because the topologically correlated states are locally indistinguishable, the results for the correlations should be distinct.
The property of only being distinguishable via spin-spin correlations, but not by single-spin quantities, is not unique to topologically ordered states but also holds for Bell states. In contrast, a property that is unique to topologically ordered states is that local perturbations cannot transform one topologically ordered state |ψ a into another one |ψ b . This crucial property can be experimentally tested with our approach.
In the rotating frame we consider, single-spin perturbations on lattice site j can be taken to read, where oscillating prefactors of σ α j have been included into the time dependence of g(t). In this rotating frame, the excited states of the Toric Code Hamiltonian (1) are separated from the degenerate topologically ordered states by an energy 4J (assuming J s = J p = J), irrespective of the lattice size. Therefore, the perturbations H pert,j (t) cannot transform a topologically ordered state into a state outside the topologically ordered manifold provided |g(t)| 4J. Yet because of the topological order, no local perturbation is capable of transforming a topologically ordered state |ψ a into another one |ψ b . As a consequence, provided their strength fulfills the condition |g(t)| 4J, no local perturbation of the form (19) can affect a topologically ordered state |ψ a . This robustness can be tested experimentally in our approach by verifying that the prepared topologically ordered states are insensitive to any local operation, provided it is weak enough.
However, as indicated above, the topologically ordered states are in our rotating frame implementation not robust against local dissipation because the coupling to an environment can change the energy in the Toric Code lattice by an amount that exceeds 4J.
This experimental verification of topological correlations can even be performed for the minimal lattice consisting of eight qubits, where it only requires measuring individual qubits and two-qubit correlations. An example of possible outcomes for these measurements are given in table III and discussed in appendix E.
In recent years several strategies for state tomography in superconducting circuits have been demonstrated experimentally [18,39,40,83]. These either proceed by reading out the qubit excitations via an auxiliary phasequbit or resonator. Both concepts are compatible with our circuit design. In particular correlations between multiple qubits, including the stabilizers A s and B p can be measured in a joint dispersive read-out via a common coplanar waveguide resonator [18,30,39,40] or several readout resonators connected to a single transmission line using frequency multiplexing [73][74][75].

VI. BRAIDING AND DETECTION OF ANYONS
Our approach is also ideally suited for probing the fractional statistics of excitations in the Toric Code. Once the lattice has been prepared in a topologically ordered state, for example via the sweep described in Sec. IV, the anyonic character of the excitations can be verified in an experiment using single qubit rotations and the measurement of one stabilizer only. We here discuss a protocol that has originally been proposed for cold atoms [84] but is very well suited for superconducting circuits.
For the prepared topologically ordered state |ψ , the protocol considers applying a single qubit rotation σ y to one of the qubits. This operation, which in our approach can be implemented via a microwave pulse in a charge line to the qubit, creates a pair of so called eparticles [5] on the neighboring vertices and puts the system into a state |ψ, e . Alternatively, half a σ y -rotation creates the state (|ψ + |ψ, e )/ √ 2. Similarly, with a σ x -rotation on another qubit, one creates a pair of mparticles. Then one can move one of these m-particles around one of the e-particles via successive σ x -rotations on the qubits along the chosen path and fuse both mparticles again. After this sequence, the system is in the state (|ψ − |ψ, e )/ √ 2 due to the fractional phase π acquired from braiding the e-and m-particles. This phase flip can for instance be detected with another √ σ y -operation that maps (|ψ − |ψ, e )/ √ 2 → |ψ and (|ψ + |ψ, e )/ √ 2 → |ψ, e and a subsequent measurement of the corresponding stabilizer A s . We note that in our approach all single qubit operations can be implemented via pulses applied through charge lines and the stabilizer measurement can be realized via the technique described in Sec. V.

VII. CONCLUSION AND OUTLOOK
We have introduced a scheme and architecture for a quantum simulator of topological order that implements the central model of this class of quantum many-body systems, the Toric Code, on a two dimensional superconducting circuit lattice. Our approach is based on a scheme that allows to realize a large variety of many-body interactions by coupling multiple qubits via dc SQUIDs that are driven by a suitably oscillating flux bias. The advanced status of experimental technology for superconducting circuits makes available an ample toolbox for exploring the physics of topological order with our scheme.
All physical qubits and coupling SQUIDs can be individually controlled with high precision, topologically ordered states can be prepared via an adiabatic ramp of the stabilizer interactions, multi-qubit correlations including stabilizers and correlations along non-contractible loops can be measured and anyons can be braided and their fractional statistics can be verified. Moreover, superconducting circuit technology allows for realizing periodic boundary conditions, a feature that not easily available in other platforms but crucial for exploring topological order. We thus expect our work to open up a realistic path towards an experimental investigations of this intriguing area of quantum many-body physics that at the same time will enable access for detailed measurements.
Our scheme is highly versatile and not restricted to implementations of the Toric Code only. Several directions for generalizations appear very intriguing at this stage.
(i) By tuning the time-independent flux bias of the coupling SQUIDs away from the operating point for the Toric Code, one can generate additional two-body interactions of the form σ z j σ z l between neighboring qubits. These can become stronger than the stabilizer interactions. The ability to tune the transition frequencies of the qubits makes our approach a very versatile quantum simulator for exploring quantum phase transitions between topologically ordered and other phases of twodimensional spin lattices.
(ii) In our concept, the many-body interactions mediated by SQUIDs are not limited to 4-body interactions. By connecting one port of a coupling SQUID to k l qubits and the other to k r qubits, k-body interactions for k = k l + k r can be implemented in kth order of the Schrieffer-Wolff expansion (7) together with a suitable choice of frequencies for the flux bias ϕ ext . Thus 3-body interactions, for example, that are an order of magnitude larger than the stabilizer interactions of equation (16) can be implemented. For the parameters considered in section III, also 6-body interactions that are stronger than dissipation would be feasible and even interactions among a larger number of qubits are within reach with some improvement of circuit parameters. These capabilities of our scheme open the door to realizations of higher dimensional Toric Codes [7,8] and other topological codes such as color codes [85]. A useful property here is that the implementation of the sum of a star and plaquette interaction in one cell would only require half the frequency components for the driving fields.
(iii) If two qubits that are connected to the same port of a coupling SQUID are identified with one lattice site of the simulated model, our approach can straight forwardly be generalized to the analogue quantum simulation of non-Abelian lattice gauge theories, where the engineered 4-body interactions play the role of the coupling matrix between the color degrees of freedom of the adjacent lattice sites.
(iv) Considering regimes, where the degrees of freedom of the coupling SQUID can no longer be adiabatically eliminated, one can modify our approach towards quantum simulations of dynamical gauge fields theories, where the coupling SQUID would represent the dynamical degrees of freedom of the gauge field. In this appendix we present a detailed derivation of the Hamiltonian in Eq. (5) that describes the considered circuit lattice. We start with the full Lagrangian describing the lattice depicted in Fig. 1, which  2E J cos( ϕ ext 2 ) cos(ϕ a;n,m − ϕ b;n,m ).
As explained in the main text, ϕ n,m is the phase at node (n, m), C q;n,m is the capacitance and E Jq;n,m the Josephson energy of the qubit at that node. ϕ a;n,m and ϕ b;n,m are the phases at the auxiliary nodes next to the SQUID located between nodes (n, m), (n, m + 1), (n + 1, m), and (n + 1, m + 1). The C g are the capacitances of these auxiliary nodes with respect to ground and C J is the capacitance of each junction in the SQUID including a shunt capacitance parallel to the junctions. E J is the Josephson energy of each junction in the SQUID.
As the Josephson junctions of the dc-SQUIDs only couple to the difference between the phases of the two adjacent auxiliary nodes, we define the modes, ϕ ±;n,m = ϕ a;n,m ± ϕ b;n,m , and describe the dc-SQUIDs in terms of ϕ +;n,m and ϕ −;n,m from now on. We now turn to derive the Hamiltonian associated to the Lagrangian in Eq. (A1). As there is neither capacitive coupling between the qubits in the lattice nor between the qubits and the coupling SQUIDs, the kinetic energy can be decomposed into the contribution of individual qubits and that of the individual coupling SQUIDs, is the capacitive matrix associated with one SQUID.
The Hamiltonian of the lattice is obtained in the standard way via a Legendre transform of the Lagrangian L in Eq. (A1) by defining the qubit and SQUID momenta, π n,m = ∂L ∂φ n,m , and π ±;n,m = ∂L ∂φ ±;n,m . (A5) The kinetic energy thus reads, where we have used that C g C J in inverting C s . We thus arrive at the Hamiltonian, where we have introduced the notation ϕ q±;n,m = ϕ n,m + ϕ n,m+1 ± ϕ n+1,m ± ϕ n+1,m+1 . (A8) is the inductive energy associated with an inductor L and the charging energies of the qubits and SQUID modes are where e is the elementary charge. The first line in Eq. (A7) describes the qubits, the second line the SQUIDs and the third line the interactions between qubits and SQUIDs via the inductances L.
The Hamiltonian (A7) is quantized in the standard way by imposing the canonical commutation relations [ϕ j , π l ] = i δ j,l and [ϕ ±;j , π ±;l ] = i δ j,l . (A9) The external flux applied through the SQUID loops and hence the associated phase is composed of a constant dcpart and a time dependent ac-part as given in Eq. (6). We choose ϕ dc = π and ϕ ac ϕ dc = π, and obtain to leading order in ϕ ac ,

Appendix B: Effective time-dependent Hamiltonian
In this appendix we present details of the derivation of the effective Hamiltonian for the qubits only, which will be dominated by the four-body interactions corresponding to the stabilizer operators in the Toric Code Hamiltonian, see Eq. (1).
As the modes ϕ ±;n,m of the SQUIDs oscillate at frequencies that strongly differ from the transition frequencies of the qubits and the oscillation frequencies of the applied flux biases, they remain in their ground states to a high degree of precision and can be adiabatically eliminated from the description. To this end we apply the Schrieffer-Wolf transformation introduced in Eq. where we have here introduced the abbreviations, π S±;n,m = π ±;n,m + π ±;n,m−1 ± π ±;n−1,m ± π ±;n−1,m−1 . Truncating the expansion (B2) at fourth order is a good approximation since the amplitudes of the terms ϕ n,m π ±;n,m / are much smaller than unity, see Eq. (B16) for details. Moreover, as we explain in the sequel, the frequencies contained in the flux bias φ ext that is applied to the SQUIDs can be chosen such that interactions contained in the fourth order terms [ The Hamiltonian of the individual qubits, H 0q , and the qubit-qubit interactions, H Iq , are given in Eqs. (8) and (9). We here repeat them for completeness, where again we have neglected terms ∝ π n,m π S+;n,m because |ω + − ω n,m | |ω − − ω n,m |, see Eq. (B15), which makes these terms highly off-resonant so that we only keep the dominant contributions ∝ π n,m π S−;n,m .

Second-quantized form
The free Hamiltonian of the qubits can be written in second-quantized form by defining the ladder operators q i and q † i via Here, q j (q † j ) is the annihilation (creation) operator for qubit j [j = (1, 1), . . . , (N, M )] and the effective Josephson energy of the qubits is E Jq;j + 4E L due to the contributions of the four adjacent inductors. In terms of these creation and annihilation operators, H 0q reads, where we have approximated cos(ϕ j ) by its expansion up to fourth order in ϕ j and the transition frequencies ω j and nonlinearities U j are given by For the SQUID modes we introduce the annihilation (creation) operators s ± (s † ± ) according to, , and get where and we have neglected the last term of Eq. (B7) since cos (ϕ ext /2) ≈ −(ϕ ac /2)F (t) and the frequencies contained in F (t) strongly differ from ω + and ω − . Note that C J C g and thusẼ C+ Ẽ C− , so that ω + ω − .

Accuracy of truncating the transformation
After having introduced the amplitudes ϕ j and ϕ ± in Eqs. (B9) and (B13), we can now estimate the amplitude of the Schrieffer-Wolff generator S, see Eqs. (7) or (B1), to check whether a truncation of the expansion (B2) at fourth order provides a good approximation. For the scenario with at most a single excitation in each qubit and vanishingly small excitation numbers in the SQUID modes, which is of interest here, the generator S has a magnitude given by 1 4 for each site. Typical parameters of our scheme lead to ϕ j 0.7 and ϕ ± 1.65 so that ϕ j /4ϕ ± ≈ 0.1 and the truncation of the expansion in Eq. (B2) is well justified.

Appendix C: Effective Hamiltonian in the rotating frame
To analyze which terms of the transformed Hamiltonian (8) respectively (B4) provide the dominant contributions, we now move to a rotating frame, where all qubits rotate at their transition frequencies ω n,m and the SQUID modes rotate at their oscillation frequencies ω ± . In transforming to this rotating frame, the raising and lowering operators transform as, q j → e −iωj t q j , and s ±;j → e iω±t s ±;j (C1) for j = (1, 1), . . . (N, M ).

Adiabatic elimination of the SQUID modes
For the parameters of our architecture, the SQUID frequencies ω ± differ strongly from the qubit transition frequencies ω j and the frequencies contained in the driving fields. The SQUID modes ϕ +;n,m and ϕ −;n,m therefore remain to a very good approximation in their ground states, |0 +;n,m and |0 −;n,m . They can thus be adiabatically eliminated from the description, where the dominant contributions come from the ϕ −;n,m modes because ω + ω − . The leading order of the adiabatic elimination is obtained by projecting H S and H qS onto the ground states of the SQUIDs, qS is a sum of contributions from each cell formed by four qubits at its corners which interact via one coupling SQUID.
When considering the terms of one cell only, as described in Sec. II A, and labelling its qubits by indices 1,2,3 and 4, i.e. (n, m) → 1, (n + 1, m) → 2, (n + 1, m + 1) → 3, (n, m + 1) → 4, one can divide the Hamiltonian H (0) qS for one cell into three parts as described in Eq. (12), where Here, H 1 contains the four-body interactions which will give rise to stabilizer interactions, H 2 contains all twobody interactions and H 3 all fourth-order interactions beyond the four-body interaction, i.e. terms where at least two of the indices i, j, k and l are the same as we indicate by the notation (ijkl) . The coefficients C ijkl have the values C ijkl = ±1/32 for terms of the from ϕ 2 j ϕ k ϕ l with j = k, j = l and k = l, C ijkl = ±1/64 for terms of the from ϕ 2 j ϕ 2 k with j = k, C ijkl = ±1/96 for terms of the from ϕ 3 j ϕ k with j = k and C ijkl = ±1/384 for terms of the from ϕ 4 j . Here, "+" signs apply whenever an even number of the indices i, j, k, l refer to the same side of the SQUID (i.e. the same diagonal) and "−" signs otherwise. Whereas we discuss perturbations originating from H 2 and H 3 in appendix D, we now turn to discuss H 1 in more detail.

Four-body interactions
Due to the nonlinearities of the qubits, we can here restrict the description to the two lowest energy levels of each node by replacing q j → σ − j (q † j → σ + j ). The Hamiltonian H 1 contains 16 terms of the form S 1 S 2 S 3 S 4 , where S j ∈ {e −iωj t σ − j , e iωj t σ + j } (see C6). Ignoring the oscillating pre-factor F (t), these terms rotate at the angular frequencies as listed in table II. Denoting the sum term frequency term frequency of the terms in the first column of table II by V 1 and the sum of the terms in the second column by V 2 it can readily be verified that Motivated by the specific form of the Hamiltonian H 1 , see Eq. (C6), we thus define two oscillating fields F X and F Y , F X (t) = f 1,1,1,1 + f 1,1,−1,−1 + f 1,−1,1,−1 + f 1,−1,−1,1 F Y (t) = f 1,1,1,−1 + f 1,1,−1,1 + f 1,−1,1,1 + f −1,1,1,1 (C7) with f a,b,c,d = cos[(aω 1 + bω 2 + cω 3 + dω 4 )t]. Each of the terms in F X (F Y ) brings two of the four-body processes in the first (second) column of table II into resonance (makes them non-rotating). Next to resonant terms, each field also generates rotating terms. The rotation frequency of these terms is the detuning between different components of F X and F Y . As the least detuning is in the GHz range (∼ 0.4 GHz) and the strength of four-body interaction is in the MHz range (∼ 1 MHz), the rotating terms can be neglected in a RWA. Hence the external oscillating field F s = F X + F Y gives rise to the star interaction A s and the external field F p = F X − F Y gives rise to the plaquette interaction B p , (C8) and the strength of the stabilizer interaction is Besides these stabilizer interactions, there are also twobody interactions between the qubits which may lead to corrections in the effective Hamiltonian (1). We now turn to estimate their effect.

Appendix D: Perturbations
In this section we estimate the effects of the remaining terms of the transformed HamiltonianH, see Eq. (8) or (B4). These include contributions contained in the Hamiltonians H Iq , see Eq. (B6), and H qS , see Eq. (B8). In our approach, the transition frequencies of the qubits ω n,m and the frequency spectrum of the oscillating fluxes ν k are chosen such the all undesired terms are highly offresonant and therefore strongly suppressed. They however lead to a frequency shift j and ∆ (4) j are given in Eqs. (D8), (D11),(D13) and (D18). These frequency shifts can easily be compensating for by modifying the frequency spectrum of the applied driving fields accordingly, so that we did not include them in our preceding analysis. They can be calculated via the following time dependent perturbation theory.
In the frame where all qubits rotate at their transition frequencies, the transformed HamiltonianH can be written asH where H TC is the Toric Code Hamiltonian as introduced in Eq.
Since ||ĥ α || ω α , the first term of H eff,2 and the first term of H eff,3 are at least two orders of magnitude smaller than H TC and we neglect them. We now turn to discuss the other perturbations to H TC one by one.

Corrections to the adiabatic elimination of the SQUIDs
Besides the terms contained in H where we have neglected the terms in Eq. (B8) that are proportional to ϕ 3 q−;n,m as these are a factor ϕ 2 n,m /24 1 smaller than the terms that are linear in ϕ q−;n,m . We now consider corrections to H T C coming from second order processes of ∆H qS .
Using the formulation in second quantization, see Eqs. (B9) and (B13), one finds for the first term of ∆H qS that the perturbation theory of Eq. (D3) is applicable if For this regime, the effective Hamiltonian H eff,1 according to Eq. (D3) for the first term of ∆H qS just contains a frequency shift for each qubit j = (n, m). Here, the prefactor 4 accounts for the fact that each qubit has 4 neighboring SQUIDs. For the set of parameters discussed in Sec. III, we find for the qubits with high transition frequencies E Cq;j /(hϕ j ϕ − ) ∼ 500 MHz leading to ∆ (1) j ∼ 200 MHz, i.e. a shift of 50 MHz from each coupling to a SQUID. For qubits with low transition frequencies the shifts are ∆ (1) j ∼ 100 MHz, i.e. a shift of 25 MHz from each coupling to a SQUID. Besides frequency shifts, the first term of ∆H qS can also lead to interactions between two qubits of one cell. These are ineffective provided any pair of qubits in a cell is sufficiently detuned from each other, which holds since any pair of qubits in one cell is at least detuned by 700 MHz from each other, see table I. Note that the detuning between two neighboring qubits on a diagonal n that would couple via two SQUIDs is significantly higher (∼ 10 GHz). The next higher order nonvanishing term of the expansion (D3) is here the second term of H eff,3 which is more than two orders of magnitude smaller than ∆ (1) j . In a similar way as second-order terms can lead to qubit-qubit interactions in one cell, fourth-order terms could mediate interactions between two qubits in neighboring cells. For our parameters, these interactions would have a strength of ∼ 0.5 MHz but are ineffective as qubits in neighboring cells are always detuned by at least 100 MHz.
Via an analogous calculation, one finds for the second term of ∆H qS that the perturbation theory of Eq. (D3) is applicable if and leads to the frequency shifts (taking into account that each qubit couples to 4 SQUIDs) ∆ (2) j = −2Ẽ for each qubit j, where the frequencies ν k that are contained in the driving field F (t) are defined in table II.
For the set of parameters as given in Sec. III, we find |∆ (2) j | 20 MHz. The next higher order non-vanishing term of the expansion (D3) is here the second term of H eff,3 which is more than two orders of magnitude smaller than ∆ (2) j .
Furthermore, there could be interactions between qubits in neighboring cells, say qubit (n, m) and qubit (n, m + 2), since qubits (n, m) and qubit (n, m + 1) couple via 2E L ϕ n,m ϕ n,m+1 and qubits (n, m + 1) and qubit (n, m + 2) couple in a second order process via their interaction with a common SQUID as discussed in Sec. D 1. This third-order interactions are strongly suppressed since their strength is ∼ 3 MHz and thus 30 times smaller than the detuning between the respective qubits.