Hardware implementation of quantum stabilizers in superconducting circuits

TL;DR

This paper demonstrates a hardware implementation of quantum stabilizers directly in superconducting circuits, enabling more robust quantum error correction through engineered Hamiltonians.

Contribution

It introduces a novel superconducting circuit design that encodes stabilizers directly into the Hamiltonian, moving beyond software-based implementations.

Findings

Energy band dispersion softening observed with increasing frustrated elements

Close agreement between experimental results and numerical modeling

Potential for improved quantum error correction hardware

Abstract

Stabilizer operations are at the heart of quantum error correction and are typically implemented in software-controlled entangling gates and measurements of groups of qubits. Alternatively, qubits can be designed so that the Hamiltonian corresponds directly to a stabilizer for protecting quantum information. We demonstrate such a hardware implementation of stabilizers in a superconducting circuit composed of chains of $\pi$-periodic Josephson elements. With local on-chip flux- and charge-biasing, we observe a softening of the energy band dispersion with respect to flux that is exponential in the number of frustrated plaquette elements, in close agreement with our numerical modeling.

Tables

Table 1. Table 1: Inductance matrix (top) from InductEx simulations of device layout and (bottom) extracted from measurements of two-dimensional flux scans of readout cavity.

Simulated Inductance Matrix (pH)
	PB01	PB12	PB23	PB30
Plaq1	0.59	0.76	-0.17	0.11
Plaq2	0.15	-0.69	-0.55	0.24
Plaq3	0.07	-0.02	0.60	0.76
Extracted Inductance Matrix (pH)
	PB01	PB12	PB23	PB30
Plaq1	0.64	0.66	-0.15	0.05
Plaq2	0.20	-0.66	-0.54	0.15
Plaq3	0.13	-0.24	0.67	0.59

Table 2. Table 2: Simulated and extracted capacitance matrix.

Simulated capacitance matrix (aF)
	${\mathrm{CB}}_{\mathrm{sh}}$	CB1	CB2
${\mathrm{Isl}}_{\mathrm{sh}}$	57	419	353
Isl1	0	45	12
Isl2	0	27	88
Extracted capacitance matrix (aF)
	${\mathrm{CB}}_{\mathrm{sh}}$	CB1	CB2
${\mathrm{Isl}}_{\mathrm{sh}}$	57	501	327
Isl1	0	35	8
Isl2	0	73	120

Table 3. Table 3: Single frustration fit parameters and estimated parameters from design and fabrication tests.

	$E_J$ (K)	$E_C$ (K)	$E_L$ (K)	$E_{CL}$ (K)	$C_{\mathrm{sh}}$ (fF)	$\alpha$	$L_{\mathrm{factor}}$
Plaquette 1	1.65	3.65	1.12	5.60	1160	0.03	1.1
Plaquette 2	1.65	3.67	1.11	6.36	1190	0.02	1.1
Plaquette 3	1.97	4.00	1.27	6.66	1440	0.04	0.91
Estimated parameters	1.45	3.82	1.39	6.46	1000	0.02	1.0

Table 4. Table 4: Double frustration fit parameters and estimated parameters. *The estimated C isl subscript 𝐶 isl C_{\rm isl} for plaquette (12), (23), and (13) double frustration are 0.81, 5.0, 5.8 fF.

	$E_J$ (K)	$E_C$ (K)	$E_L$ (K)	$E_{CL}$ (K)	$C_{\mathrm{sh}}$ (fF)	$\alpha$	$C_{\mathrm{isl}}$ (fF)	$L_{\mathrm{factor}}$
Plaquette (12)	1.75	3.54	1.20	6.34	1240	0.03	1.5	0.99
Plaquette (23)	1.76	3.53	0.900	7.40	1290	0.04	5.7	1.0
Plaquette (13)	1.73	3.48	0.903	6.62	1310	0.03	8.1	0.98
Estimated parameters	1.45	3.82	1.39	6.46	1000	0.02	Vary*	1.0

Table 5. Table 5: Fit parameter error estimation based on (a) plasmon transition, (b) heavy fluxon transition, (c) light fluxon transition, and (d) anticrossing at flux points indicated in Fig.  S20 . (e) Intersection of fitted errors estimated from four types of transitions in (a-d).

Parameter uncertainty estimates within a 10% range based on:
$E_J$ (K)	$E_C$ (K)	$E_L$ (K)	$E_{CL}$ (K)	$C_{\mathrm{sh}}$ (fF)	$\alpha$	$C_{\mathrm{isl}}$ (fF)	$C_{\mathrm{int}}$ (fF)	$L_{\mathrm{factor}}$
(a) Plasmon transition ($|0_{ES}⟩\to |3_{ES}⟩$ transition at 20 m${\Phi }_0$)
1.65 $\sim$ 1.85	3.00 $\sim$ 4.00	0.770 $\sim$ 1.35	3.50 $\sim$ 15.0	1100 $\sim$ 1400	0.0 $\sim$ 0.1	0.200 $\sim$ 6.00	0.800 $\sim$ 2.10	0.83 $\sim$ 1.1
(b) Heavy fluxon transition ($|0_{ES}⟩\to |2_{OS}⟩$ transition at 6.5 m${\Phi }_0$)
1.73 $\sim$ 1.77	3.30 $\sim$ 3.70	1.10 $\sim$ 1.30	5.50 $\sim$ 8.00	1050 $\sim$ 1300	0.00 $\sim$ 0.04	1.26 $\sim$ 2.09	0.880 $\sim$ 1.13	0.83 $\sim$ 1.1
(c) Light fluxon transition ($|0_{ES}⟩\to |0_{EA}⟩$ transition at 0 m${\Phi }_0$)
1.73 $\sim$ 1.77	3.45 $\sim$ 3.62	1.15 $\sim$ 1.22	5.40 $\sim$ 7.80	1000 $\sim$ 1500	0.00 $\sim$ 0.05	1.48 $\sim$ 1.66	0.970 $\sim$ 1.06	0.71 $\sim$ 1.1
(d) Anticrossing between the $|1_{ES}⟩\to |4_{ES}⟩$ transition and the $|1_{ES}⟩\to |2_{OS}⟩$ transition at 12 m${\Phi }_0$
1.71 $\sim$ 1.76	3.50 $\sim$ 3.85	1.14 $\sim$ 1.38	5.80 $\sim$ 12.5	1000 $\sim$ 1500	0.02 $\sim$ 0.06	0.980 $\sim$ 1.67	0.880 $\sim$ 1.01	0.91 $\sim$ 1.2
(e) Intersection of fitted errors estimated from four types of transitions
1.73 $\sim$ 1.76	3.50 $\sim$ 3.62	1.15 $\sim$ 1.22	5.80 $\sim$ 7.80	1100 $\sim$ 1300	0.02 $\sim$ 0.04	1.48 $\sim$ 1.66	0.970 $\sim$ 1.01	0.91 $\sim$ 1.1

Table 6. Table 6: T 2 Φ superscript subscript 𝑇 2 Φ T_{2}^{\Phi} calculation from slope and curvature at 10  μ ​ Φ 0 𝜇 subscript Φ 0 \mu\Phi_{0} away from frustration for a protected qubit with the optimal parameters presented in Sec.  XI. . Λ T 2 superscript Λ subscript 𝑇 2 \Lambda^{T_{2}} is calculated from the ratio of T 2 Φ superscript subscript 𝑇 2 Φ T_{2}^{\Phi} values from the modeled energy levels, as described in Sec.  XI. ; Λ Δ SA superscript Λ subscript Δ SA \Lambda^{\Delta_{\rm SA}} is computed from Δ SA subscript Δ SA \Delta_{\rm SA} and h Z subscript ℎ 𝑍 h_{Z} following Eqs. ( S11 , S14 , S15 , S16 ).

	${\Delta }_{\mathrm{EO}}$ at 10 $\mu {\Phi }_0$ (MHz)	${\Delta }_{\mathrm{SA}}$ (GHz)	Slope at 10 $\mu {\Phi }_0$ (MHz/$m{\Phi }_0$)	Curvature at 10 $\mu {\Phi }_0$ (MHz/$m{\Phi }_0^2$)	$T_2^{\Phi }\left(\mu s\right)$	${\Lambda }_{\mathrm{DS}}^{T_2}$	${\Lambda }_{\mathrm{DS}}^{{\Delta }_{\mathrm{SA}}}$	${\Lambda }_{\mathrm{TD}}^{T_2}$	${\Lambda }_{\mathrm{TD}}^{{\Delta }_{\mathrm{SA}}}$	${\Lambda }_{\mathrm{TS}}^{T_2}$	${\Lambda }_{\mathrm{TS}}^{{\Delta }_{\mathrm{SA}}}$
plaquette 2	$3.7$	-	380	-	1.7	-	-	-	-	-	-
plaquette (12)	$6.7$	1.5	1.8	$1.8\times {10}^2$	340	200	200	180	150	-	-
plaquette (23)	$5.0$	0.87	3.2	$3.0\times {10}^2$	190	110	120	320	270	-	-
plaquette (123)	$1.5\times {10}^{-2}$	-	$1.0\times {10}^{-2}$	$1.2\times {10}^1$	$6.0\times {10}^4$	-	-	-	-	$3.5\times {10}^4$	$3.1\times {10}^4$

Equations

$$\vec{\Phi }=\mathrm{L}\vec{I}+\vec{x},$$

$$1/T_1=\text{absolutevalue}{⟨{\psi }_{\mathrm{odd}}\mid \hat{N}\mid {\psi }_{\mathrm{even}}⟩}^2{\Delta }_{\mathrm{EO}}\tan \delta .$$

$$\frac{1}{T_2^{\Phi }}\approx \frac{A_{\Phi }}{h}\times {\frac{\partial E}{\partial \Phi }|}_{\Phi ={\Phi }_{\mathrm{noise}}}$$

$$\frac{1}{T_2^Q}\approx A_Q\times \frac{{\Delta }_{\mathrm{EO}}}{e},$$

$$\begin{array}{rl}H& =-\frac{k\Delta {\Phi }_1}{2}{Z}_1-\frac{k\Delta {\Phi }_2}{2}{Z}_2-\frac{{\Delta }_{\mathrm{SA}}}{2}{X}_1{X}_2.\end{array}$$

$$\begin{array}{rlr}& E_1=-\sqrt{{\left(\frac{k\Delta {\Phi }_1}{2}+\frac{k\Delta {\Phi }_2}{2}\right)}^2+{\left(\frac{{\Delta }_{\mathrm{SA}}}{2}\right)}^2},& E_2=\sqrt{{\left(\frac{k\Delta {\Phi }_1}{2}+\frac{k\Delta {\Phi }_2}{2}\right)}^2+{\left(\frac{{\Delta }_{\mathrm{SA}}}{2}\right)}^2},\\ & E_3=-\sqrt{{\left(\frac{k\Delta {\Phi }_1}{2}-\frac{k\Delta {\Phi }_2}{2}\right)}^2+{\left(\frac{{\Delta }_{\mathrm{SA}}}{2}\right)}^2},& E_4=\sqrt{{\left(\frac{k\Delta {\Phi }_1}{2}-\frac{k\Delta {\Phi }_2}{2}\right)}^2+{\left(\frac{{\Delta }_{\mathrm{SA}}}{2}\right)}^2}.\end{array}$$

$$\begin{array}{rlrlr}& E_1=-|\frac{k\Delta {\Phi }_1}{2}+\frac{k\Delta {\Phi }_2}{2}|,& E_2=|\frac{k\Delta {\Phi }_1}{2}+\frac{k\Delta {\Phi }_2}{2}|,& E_3=-|\frac{k\Delta {\Phi }_1}{2}-\frac{k\Delta {\Phi }_2}{2}|,& E_4=|\frac{k\Delta {\Phi }_1}{2}+\frac{k\Delta {\Phi }_2}{2}|.\end{array}$$

$${\left|\frac{\partial \Delta E}{\partial \Delta \Phi }\right|}_{\mathrm{S}}=k,$$

$$\begin{array}{rlrl}{E}_1& =-\sqrt{(k\Delta \Phi {)}^2+{\left(\frac{{\Delta }_{\mathrm{SA}}}{2}\right)}^2},& E_2& =\sqrt{(k\Delta \Phi {)}^2+{\left(\frac{{\Delta }_{\mathrm{SA}}}{2}\right)}^2},\\ E_3& =-\frac{{\Delta }_{\mathrm{SA}}}{2},& E_4& =\frac{{\Delta }_{\mathrm{SA}}}{2}.\end{array}$$

$${\left|\frac{\partial \Delta E}{\partial \Delta \Phi }\right|}_{\mathrm{D}}=\frac{2k^2\Delta {\Phi }_{\mathrm{noise}}}{{\Delta }_{\mathrm{SA}}},$$

$${\Lambda }_{\mathrm{DS}}={\left|\frac{\partial \Delta E}{\partial \Delta \Phi }\right|}_{\mathrm{S}}/{\left|\frac{\partial \Delta E}{\partial \Delta \Phi }\right|}_{\mathrm{D}}=k\times \frac{{\Delta }_{\mathrm{SA}}}{2k^2\Delta {\Phi }_{\mathrm{noise}}}=\frac{{\Delta }_{\mathrm{SA}}}{2k\Delta {\Phi }_{\mathrm{noise}}}.$$

$$H=-\frac{k\Delta {\Phi }_1}{2}{Z}_1-\frac{k\Delta {\Phi }_2}{2}{Z}_2-\frac{k\Delta {\Phi }_3}{2}{Z}_3-\frac{{\Delta }_{\mathrm{SA}}^{(12)}}{2}{X}_1{X}_2-\frac{{\Delta }_{\mathrm{SA}}^{(23)}}{2}{X}_2{X}_3,$$

$$\text{absolutevalue}{\frac{\partial \Delta E}{\partial \Delta \Phi }}_{\mathrm{T}}=\frac{3k^3\Delta {\Phi }_{\mathrm{noise}}^2}{{\Delta }_{\mathrm{SA}}^{(12)}{\Delta }_{\mathrm{SA}}^{(23)}},$$

$${\Lambda }_{\mathrm{TD}}^{(12)}=\frac{2{\Delta }_{\mathrm{SA}}^{(23)}}{3k\Delta {\Phi }_{\mathrm{noise}}},$$

$${\Lambda }_{\mathrm{TD}}^{(23)}=\frac{2{\Delta }_{\mathrm{SA}}^{(12)}}{3k\Delta {\Phi }_{\mathrm{noise}}},$$

$${\Lambda }_{\mathrm{TS}}=\frac{{\Delta }_{\mathrm{SA}}^{(12)}{\Delta }_{\mathrm{SA}}^{(23)}}{3k^2\Delta {\Phi }_{\mathrm{noise}}^2}=\frac{3}{4}{\Lambda }_{\mathrm{TD}}^{(12)}{\Lambda }_{\mathrm{TD}}^{(23)}.$$

Full text

Protecting fragile information in quantum processors requires some form of quantum error correction (QEC). With typical “software” QEC techniques such as the surface code [1], stabilizing a single logical qubit requires many physical qubits, each of which is typically implemented as a weakly nonlinear oscillator. Error correction and computation is achieved by a string of operations and measurements that allow identification of bit-flip and phase-flip errors. An alternative is to implement quantum stabilizers directly in hardware. Here, error correction arises from the natural quantum dynamics, reducing the need for repeated entangling gates, measurements, and a multitude of control lines and complex classical control hardware. In this approach, the highly non-trivial Hamiltonian results in a tiny protected subspace within a huge Hilbert space.

Both approaches can be characterized by the error suppression factor $\Lambda$, the rate at which the logical error decreases with system size. The long time required by each round of software error correction for current transmon qubit arrays implies that $\Lambda$ is only marginally greater than one [2]. In this work, we experimentally demonstrate the potential to achieve much larger $\Lambda\gtrsim$100 with the Hamiltonian approach. The price that one pays is the appearance of relatively low energy modes with gaps $\lesssim 1\rm{\;GHz}$ that make initialization challenging; these gaps can be made higher through parameter optimization. Before building a scalable logical qubit with hardware QEC, it is crucial to demonstrate the effectiveness of protection based on Hamiltonian engineering as system size increases. In this manuscript, we observe and quantify the stabilizing interaction Hamiltonian between unprotected elements. We perform spectroscopic measurements with local flux control and observe signatures of stabilizer terms in the Hamiltonian. Specifically, we find a progressive flattening of the energy bands with respect to flux as system size increases, consistent with linear flux dispersion for a system size of one, quadratic for two, and cubic for three. In addition, we observe a characteristic periodic modulation with offset charge as we tune between regimes with different levels of protection.

A variety of qubit designs with intrinsic protection against decoherence have been studied previously [3, 4], including the 0$-\pi$ qubit [5, 6, 7], the two-Cooper-pair tunneling qubit [8], the bifluxon qubit [9], and rhombi arrays [10, 11, 12]. In this last work, previous devices had limited symmetry due to the inability to tune each element to the optimum flux independently; in addition, the devices were sensitive to offset charge fluctuations on internal nodes in each element, and the suppression of tunneling between the logical states was limited. Similar to previous protected qubit designs, our device is based on $\pi$-periodic Josephson elements [13], for which the Josephson energy is proportional to $\cos 2\varphi$, where $\varphi$ is the superconducting phase difference across the element. Here, charge transport consists of coherent tunneling of $4e$, as opposed to $2e$ for a conventional junction. We implement each element as a plaquette formed from a dc Superconducting QUantum Interference Device (SQUID), consisting of two conventional Josephson junctions and a non-negligible loop inductance. When flux-biased at frustration, $\Phi_{0}/2$ ($\Phi_{0}\equiv h/2e$), the first harmonic of the Josephson energy (proportional to $\cos\varphi$) vanishes. This leaves a second order term $E_{2}\cos 2\varphi$, with sequential minima separated by $\pi$; $E_{2}$ depends on the Josephson energy of the individual junctions $E_{J}$ and the energy of the SQUID inductance $E_{L}$ (Supplement [14], Sec. XI); $\varphi$ is thus a compact variable residing on a circle. Biasing below (above) $\Phi_{0}/2$ raises (lowers) the $\pi$ wells relative to the [math] wells; for flux bias at 0$\,{\rm mod}\,\Phi_{0}$, the potential becomes proportional to $\cos\varphi$. A small asymmetry between the two junctions has a similar but less severe effect on the $\cos 2\varphi$ potential compared to a small flux deviation from frustration (Supplement [14], Sec. I).

For a single frustrated plaquette with a large capacitive shunt $C_{\rm sh}$ [Fig. 1(a)], tunneling between the ground states in the $0,\pi$ wells is suppressed. In the phase basis, wavefunctions localized in the $0,\pi$ wells are thus disjoint and well protected against bit-flip errors. At the same time, the wavefunctions are spread out in the charge basis, corresponding for the 0($\pi$) states to superpositions of even (odd) multiples of Cooper pairs on the logical island where the plaquette connects to $C_{\rm sh}$. For bias away from frustration, the energy levels disperse linearly [Fig. 1(d)], with no protection against phase flips due to flux noise.

We next consider concatenation of multiple plaquettes while maintaining the large shunt $C_{\rm sh}$ across the array. At double frustration, when two plaquettes are simultaneously biased to $\Phi_{0}/2$, there are four minima in the two-dimensional surface defined by the phase drops across each plaquette: $00,\pi\pi,0\pi,\pi 0$. For the two-plaquette circuit this has the topology of a torus, since $\varphi$ for each plaquette is a compact variable with $2\pi$ periodicity [Fig. 1(f,g)]. If the capacitance of the intermediate island between plaquettes to ground $C_{\rm isl}$ is sufficiently small, with charging energy $E_{C}^{\rm isl}=(2e)^{2}/2C_{\rm isl}>E_{J}$, quantum fluctuations of the island phase cause hybridization along the direction between wells of the same parity; that is, $00$ will hybridize with $\pi\pi$ and $0\pi$ with $\pi 0$. Levels with the same parity develop a splitting near double frustration, with ground states corresponding to the symmetric superpositions $00+\pi\pi$ ($0\pi+\pi 0$) for even (odd) parity. Excited states are given by the antisymmetric superpositions 00-$\pi\pi$ ($0\pi$-$\pi 0$) for even (odd) parity; these states are separated by an energy $\Delta_{\rm SA}$ from the symmetric ground state of the same parity. The hybridized ground state wavefunctions of opposite parity are the logical states for the device [Fig. 1(h)] and form interlocking rings on the torus [Fig. 1(g)]. Due to delocalization and intertwining of the hybridized ground state wavefunctions, local perturbations affect the logical states symmetrically. Larger $E_{C}^{\rm isl}$ increases $\Delta_{\rm SA}$ and further flattens the bands [Fig. 1(i)], thus protecting against dephasing from flux noise.

Treating each plaquette as a spin-1/2 particle, the $\Delta_{\rm SA}$ splitting corresponds to an $XX$ stabilizer term in the Hamiltonian of frustrated plaquettes $i,j$: $H_{XX}=-(\Delta_{\rm SA}^{(ij)}/2)X_{i}X_{j}$, where $X_{i}$ is the Pauli $\sigma_{x}$ matrix for plaquette $i$. The error suppression factor $\Lambda$ can be approximated as the ratio of $\Delta_{\rm SA}^{(ij)}$ to twice the scale $h_{Z}$ of dephasing fluctuations for single plaquette $i$, $\delta H(t)=h_{Z}(t)Z_{i}$, which, for this device, will be dominated by flux noise (Supplement [14], Sec. XII). $C_{\rm sh}$ still suppresses tunneling between logical states of opposite parity, protecting against bit-flip errors.

In our experiments, we target a three-plaquette circuit with $E_{J}\sim E_{L}\sim 1.5\,{\rm K}$ ($k_{B}$=1), where $E_{L}$ is the energy $\left(\Phi_{0}/2\pi\right)^{2}/L$ of the inductance $L$ on each plaquette arm. We aim for a charging energy of each plaquette junction $E_{C}=(2e)^{2}/2C_{j}\sim 3.5\,{\rm K}$, where $C_{j}$ is the junction capacitance. These values can be achieved with conventional Al-AlOx-Al junctions. We implement the inductors with chains of large-area junctions, similar to fluxonium [17], thus eliminating charge fluctuations on the internal nodes between each small junction and inductor within a plaquette. The shunt capacitor $C_{\rm sh}$=1.2 pF is capacitively coupled to a resonator. There are four flux-bias lines, each of which couples strongly to one or two plaquettes. There are three charge-bias lines: one to the logical island that forms $C_{\rm sh}$, and one to each intermediate island between plaquettes (Supplement [14], Sec. II-V).

For device tune-up, we scan various pairs of flux-bias lines while monitoring the dispersive shift of the readout resonator. Each blue line in Fig. 2(a,b) corresponds to one plaquette passing through frustration. A crossing of two (three) lines indicates double (triple) frustration. The spacing between parallel sets of lines defines the period $\Phi_{0}$. We fit the slopes and spacing of the lines to extract the inductance matrix mapping bias levels on each flux line to net flux coupled to each plaquette (Supplement [14], Sec. VI). By inverting this matrix, we determine bias parameters for moving along arbitrary flux vectors.

We next map out the flux dispersion of the level transitions for different frustration conditions. With our ability to adjust the various plaquette fluxes independently using local flux-biasing, we maintain some plaquettes at unfrustration (0$\,{\rm mod}\,\Phi_{0}$), where the plaquette behaves like a conventional Josephson element, while we scan the flux of other plaquettes near frustration. In Fig. 3, we consider the expected level structure and define the types of possible transitions. We refer to transitions between levels in the same well as plasmons; transitions between different wells are referred to as heavy fluxons because of the vanishingly small gap associated with the corresponding anticrossing, a consequence of the large effective mass from $C_{\rm sh}$. Transitions between hybridized levels of the same parity but opposite symmetry, for example, $00$+$\pi\pi$ to $00$-$\pi\pi$, disperse sharply with flux; these are known as light fluxons due to the low effective mass in the $\varphi_{2}=-\varphi_{1}$ direction from the smallness of $C_{\rm isl}$.

To perform spectroscopy, we drive a microwave probe tone into the charge bias line coupled to $C_{\rm sh}$ while monitoring the cavity dispersive shift. Near single frustration, we initialize in the $\pi$ well prior to each spectroscopy pulse by setting the bias to $0.1\,\Phi_{0}$ from frustration, thus moving out of the protected space; we then quickly ramp the bias to the measurement point and apply spectroscopy and readout pulses (Supplement [14], Sec. VII). In Fig. 4(a), we show single-frustration measurements for plaquette 2. Features that disperse gradually correspond to plasmons within the $\pi$ well where the qubit is initialized. We continue to observe transitions out of the $\pi$ well even when the device is biased past frustration, where the $\pi$ well is higher in energy than the 0 well, due to suppressed tunneling between states of opposite parity. In addition to the 0-1, 0-2, and 0-3 transitions, we observe transitions out of excited states in the well, such as 1-2, 1-3, and 1-4, and even 2-3 and 2-4, due to insufficient cooling into the ground state of the $\pi$ well. Because of the spurious excitations to multiple levels, we are unable to apply initialization techniques that are commonly used for other low-gap qubits, such as heavy fluxonium [31, 32]. Nevertheless, we observe only weak transitions out of the 0 well, indicating that we are predominantly preparing the circuit in the $\pi$ well. In addition to the plasmons, we also observe heavy fluxons that disperse linearly with flux, which arise from transitions between various levels in the $\pi$ and 0 wells, where the barrier to tunneling is small because the initial state is an excited level or the wells are tilted by the flux bias; note that we do not observe the heavy fluxon between the protected ground states in the 0 and $\pi$ wells, which are the logical levels. We observe similar behavior for plaquettes 1 and 3 (Supplement [14], Sec. X).

The curves included in Fig. 4(a) are generated from detailed numerical modeling of the device energy levels (Supplement [14], Sec. IX). With the ability to calculate the level spectrum, we adjust the circuit parameters to fit the measured transitions from the spectroscopic data (Supplement [14], Sec. X). We observe excellent agreement, even capturing splittings that result when a fluxon crosses a plasmon due to resonant tunnel coupling between aligned levels in the 0 and $\pi$ wells. In addition, these splittings depend on the offset charge on the $C_{\rm sh}$ island [Fig. 4(e)] due to Aharonov-Casher (A-C) interference [33, 34] between tunneling paths clockwise (CW) or counterclockwise (CCW) in the $\cos 2\varphi$ potential [Fig. 1(c)] (Supplement [14], Sec. VIII). At single frustration, as expected, the heavy fluxon dispersion is linear down to zero energy, thus offering no protection against flux noise.

Upon tuning to double frustration, we observe a qualitatively different behavior. We initialize in the $\pi\pi$ well of the two-plaquette potential, then quickly ramp near double frustration. We scan both plaquette fluxes in tandem along the direction between the regimes with a global potential minimum at $\pi\pi$ and 00 and passing through double frustration. Spectroscopy at plaquette (12) double frustration shows plasmons similar to the single frustration measurements [Fig. 4(b)]. However, unlike single frustration, where suppressed tunneling between the $0,\pi$ wells allows the device to remain in the $\pi$ well even after the flux is ramped well past frustration, at double frustration, the large symmetric-antisymmetric gap $\Delta_{\rm SA}^{(12)}$ causes an adiabatic transition from $\pi\pi$ to 00 upon passing through double frustration. At higher frequencies, we observe steeply dispersing light fluxons, with the minimum at double frustration corresponding to $\Delta_{\rm SA}^{(12)}$ from hybridization of the 00 and $\pi\pi$ wells. For scans along the odd-parity flux direction, or if the circuit is initialized in an odd-parity well and scanned in the even-parity flux direction, the spectral features become swapped [Fig. 1(i), Supplement [14], Sec. X.C].

As with spectroscopy at single frustration, we include curves for the various transitions from numerical modeling and fitting for double frustration [Fig. 4(b)]. Here, the larger Hilbert space requires a significant increase in computational resources. Our modeled transition curves agree well with the measured spectroscopy, capturing both the plasmons and heavy fluxons. We are unable to directly drive a microwave transition between the logical states in the $00$+$\pi\pi$ and $0\pi$+$\pi 0$ wells due to the vanishing matrix element, the basis of protection. However, the increasing flatness of the higher fluxon transitions as one moves lower in the spectrum indicates that the logical levels will be the flattest. This can also be seen in the blue modeled curves near the bottom of the figure highlighting the dispersion of the logical level transition, which exhibits quadratic curvature. Additionally, our modeling captures the light fluxons to the antisymmetric levels.

The effectiveness of concatenation depends on $C_{\rm isl}$ of the intermediate island between the two frustrated plaquettes. For plaquette (12) double frustration, $\Delta_{\rm SA}^{(12)}$ is 2.7 GHz. At plaquette (23) double frustration, which involves a significantly larger $C_{\rm isl}$ because of the orientation of the plaquette 2 inductors, we observe a smaller $\Delta_{\rm SA}^{(23)}$ and a correspondingly larger curvature of the heavy fluxon transition. $\Delta_{\rm SA}^{(13)}$ is even smaller because of the excess capacitance to ground of the unfrustrated plaquette 2 (Supplement [14], Sec. X). Figure 4(h) shows the variation of $\Delta_{\rm SA}$ with $E_{C}^{\rm isl}$, including measured values of $\Delta_{\rm SA}^{(ij)}$ for each combination of double frustration, as well as numerically modeled values. For a typical flux noise level, $h_{Z}$ for these plaquettes will be $\sim$2 MHz, which, when combined with the measured $\Delta_{\rm SA}^{(12)}$, is consistent with $\Lambda\sim 700$. Note that this is an extracted parameter characterizing protection in one channel: dephasing. The complete $\Lambda$-parameter for a logical qubit must be derived from the scaling of $T_{1}$ and $T_{2}$ with system size, which is beyond the scope of this manuscript. Nonetheless, $\Lambda$ can also be expressed as the ratio of $T_{2}$ for a higher degree of frustration relative to $T_{2}$ at single frustration (Supplement [14], Sec. XII).

In addition to the symmetric/antisymmetric gap, another characteristic of the stabilizer term is the periodic modulation of $\Delta_{\rm SA}^{(ij)}$ with offset charge on the intermediate island between plaquettes $i$ and $j$. Destructive A-C interference of tunneling paths in the CW and CCW directions on the constant-parity circles for double frustration [Fig. 1(g)] causes $\Delta_{\rm SA}^{(ij)}$ to vanish for island offset charge near $e\,{\rm mod}\,2e$. We observe periodic modulation with charge bias to the islands with a spectroscopy pulse on the 0-1 transition [Fig. 4(f)]. While the island offset charge is stable on timescales up to one hour, it is critical there are no jumps to near $e\,{\rm mod}\,2e$. Thus, it is important to actively stabilize these offset charges through periodic calibrations (Supplement [14], Sec. VIII, IX).

By simultaneously frustrating all plaquettes, we measure spectroscopy near triple frustration [Fig. 4(c)]. In this case, we are unable to numerically fit the level spectrum since the Hilbert space size becomes prohibitively large. Nonetheless, we are able to compute the spectrum using parameter values from previous fits to double and single frustration, although the calculation takes several weeks to complete. We obtain reasonable agreement with the measurements, although the spectral features are more challenging to resolve compared to other degrees of frustration; the higher transitions are off by $\sim$5-10%, which is not unreasonable considering the circuit complexity and intertwined wavefunctions, given limitations on the number of quantum states needed for the computation to converge. Around 1.5 GHz, we observe a prominent central flat feature of width $\sim$7 m$\Phi_{0}$ around the 0-3 transition, which is uncharacteristic for parabolic, let alone linear, dispersion; below this, the 0-1 transition around 0.6 GHz is similarly flat. The transition between the logical states, which cannot be directly driven due to protection of these states from the environment, will be comparably flat (Supplement [14], Sec. X.D). Also, the light fluxon transitions are qualitatively different compared to double frustration. We additionally observe charge modulation with two different periods and slopes corresponding to separate tuning of offset charge on each intermediate island [Fig. 4(g)], characteristic of a Hamiltonian with two stabilizer terms: $H_{XX}=-(\Delta_{\rm SA}^{(12)}/2)X_{1}X_{2}-(\Delta_{\rm SA}^{(23)}/2)X_{2}X_{3}$. For our present device $\Delta_{\rm SA}^{(23)}$ is smaller than $\Delta_{\rm SA}^{(12)}$ due to excess ground capacitance from plaquette 2, resulting in the logical level dispersion at triple frustration being only marginally flatter than at double frustration [Fig. 4(d)] (Supplement [14], Sec. XI).

While our present device successfully demonstrates the implementation of stabilizer terms in hardware, development of protected qubits based on hybridized ground states of opposite parity requires larger gaps to the excited states. This, in conjunction with weaker radiative coupling to parasitic high-frequency modes from a more compact $C_{\rm sh}$, perhaps achieved using a parallel-plate rather than planar design, will avoid spurious excitations to multiple excited levels that complicate the initialization process for our present device. A device with higher excited-state energies that can be operated in the qubit regime requires larger $E_{J}$, ideally at least 3 K. We must also maintain even larger $E_{C}$ to have large $\Delta_{\rm SA}$ at double frustration with the resulting flat dispersion. For a qubit with these improved parameters subject to typical flux- and charge-noise levels, optimistic but feasible junction asymmetries, and dielectric loss from a parallel-plate $C_{\rm sh}$, we project $\Lambda\gtrsim$100, corresponding to $T_{1}\gg 1\,{\rm s}$ and $T_{2}\sim 60\,{\rm ms}$ (Supplement [14], Sec. XI), well beyond current state-of-the-art superconducting qubits.

This work is supported by the U.S. Government under ARO grant W911NF-18-1-0106. Fabrication was performed in part at the Cornell NanoScale Facility, a member of the National Nanotechnology Coordinated Infrastructure (NNCI), which is supported by the National Science Foundation (Grant NNCI-2025233). Portions of this work were supported by the National Science Foundation, Quantum Leap Challenge Institute for Hybrid Quantum Architectures and Networks, Grant No. 2016136.

I. $\pi$-periodic Josephson elements from dc SQUIDs

In our device, we implement each $\pi$-periodic Josephson element with a plaquette formed from a dc Superconducting QUantum Interference Device (SQUID), consisting of two conventional Josephson junctions and a non-negligible loop inductance [Fig. S1(a)]. Each junction has a critical current $I_{0}$ and $E_{J}=\Phi_{0}I_{0}/2\pi$; the inductance in each arm of the SQUID $L$ is related to the inductive energy $E_{L}=(\Phi_{0}/2\pi)^{2}/L$. In order to understand the origin of the $\cos 2\varphi$ potential, we consider the two-dimensional potential energy landscape as a function of the two junction phases, $\delta_{1}$ and $\delta_{2}$, which is determined by $E_{J}$, $E_{L}$, and the external flux bias $\Phi_{\rm ex}$ [1]. For now, we consider symmetric plaquettes where both junction critical currents are identical; later in this section we will consider the effects of junction asymmetry. Following convention for dc SQUIDs we plot the potential energy in terms of the common-mode and differential phase variables: $\delta_{p}=(\delta_{1}+\delta_{2})/2$ and $\delta_{m}=(\delta_{2}-\delta_{1})/2$. The phase dependence of the Josephson energy for each junction results in a 2D washboard pattern of potential minima. At the same time, the inductive energy associated with circulating currents flowing through the inductors corresponds to a parabolic sheet with its minimum along a line running parallel to $\delta_{p}$. Changing $\Phi_{\rm ex}$ shifts where the minimum of this inductive parabolic sheet falls with respect to the minima of the Josephson washboard, and thus determines the pattern of the global minima in the potential.

For a flux bias at unfrustration $\Phi_{\rm ex}=0\,{\rm mod}\Phi_{0}$, the minima are centered on $\delta_{m}=0$ and are spaced by $2\pi$ in $\delta_{p}$ [Fig. S1(b)]. Along $\delta_{m}$, there is only the one minimum at $\delta_{m}=0$ [Fig. S1(d)], corresponding to no circulating current around the SQUID loop. Along $\delta_{p}$ for $\delta_{m}=0$, the potential follows a $\cos\delta_{p}$ dependence. Thus, at unfrustration, the plaquette behaves like a single Josephson junction with critical current $2I_{0}$. When flux biased at $\Phi_{0}/2$, the plaquette exhibits a staggered pattern of energy minima about a line along $\delta_{p}$ for $\delta_{m}=\pi/2$ [Fig. S1(c)]. Figure S1(e) shows a linecut along a line between two adjacent minima as a function of $\delta_{m}$; the two minima correspond to opposite senses of circulating current around the plaquette loop, similar to a flux qubit [2] or fluxonium [3]. However, unlike these other qubits, these plaquettes also have another independent phase degree of freedom from $\delta_{p}$, which corresponds to the phase drop across the plaquette. Along $\delta_{p}$, the potential is simply $E_{2}\cos 2\varphi$, with sequential minima separated by $\pi$ [Fig. S1(g)], where the energy scale $E_{2}$ depends on the Josephson energies of the individual Josephson junctions $E_{J}$ and the inductive energy of the SQUID loop inductance $E_{L}$.

While the behavior described here is generic for any dc SQUID, achieving a $\cos 2\varphi$ potential at frustration with a significant barrier height $E_{2}$ requires a sufficiently large ratio $E_{J}/E_{L}$. In the conventional language of dc SQUIDs, screening effects are characterized by the parameter $\beta_{L}=2LI_{0}/\Phi_{0}=E_{J}/\pi E_{L}$. For SQUIDs in the limit $\beta_{L}\rightarrow 0$ and perfect symmetry, the critical current of the SQUID will modulate to zero at frustration. For such a device, not only is the first-order Josephson energy suppressed, but $E_{2}$ will be vanishingly small as well, and thus not support bound states in a $\cos 2\varphi$ potential. In order to have a significant $E_{2}$, $E_{J}/E_{L}$ must be of order unity. The dc SQUID in Fig. S1 has $E_{J}/E_{L}=\pi$ to highlight the development of the $\pi$-periodicity at frustration.

We next consider deviations from this ideal $\pi$-periodic plaquette behavior. With the flux bias moved below (above) frustration ($\Phi_{0}/2$), the $\pi$ wells are raised above (below) the [math] wells [Fig. S2(c)]. To account for asymmetries between the two junctions in a plaquette we define $\alpha=\left(E_{J2}-E_{J1}\right)/\left(E_{J2}+E_{J1}\right)$, where $E_{J1}\,(E_{J2})$ is the Josephson energy of the left (right) junction. With a non-zero $\alpha$, the common-mode potential along $\delta_{p}$ for $\Phi_{\rm ex}=\Phi_{0}/2$ has equal minima for the 0 and $\pi$ wells, but now the barrier heights between wells become asymmetric.

II. Device Fabrication

This device was fabricated on a high resistivity ($\geq$10 k$\Omega$-cm) silicon wafer that was given a standard RCA clean followed by an etch step in a buffered-2% per volume HF bath to remove native oxides immediately before loading into a vacuum chamber for the base-layer metal deposition. The base layer of 60-nm thick niobium is sputter-deposited and is then coated with DSK101-4 anti-reflective-coating (ARC) and DUV210-0.6 photoresist before performing deep-UV photolithography on a photostepper to define the ground plane, feedline, resonator, flux/charge bias lines, and the logical islands. The exposed wafer is then baked at 135∘C for 90 seconds, developed with AZ 726 MiF, briefly cleaned with an ARC etch to remove any remaining unwanted ARC and then dry etched using a BCl3, Cl2, and Ar in an inductively coupled plasma etcher. The wafer is then subject to another buffered HF dip to remove any further oxides that may have formed on the surface of the remaining niobium.

The next set of lithography steps creates ground straps that connect ground planes on either side of the flux, charge, and feedlines. The first step uses lift-off resist LOR3A and then DUV210-0.6 photoresist to expose a region underneath the intended ground straps where we deposit SiO2 to function as an insulating dielectric support for the aluminum ground straps to follow. The SiO2 is evaporated in an electron beam evaporator at a rate  3.5 ${\rm\AA}$/s until 100 nm is deposited. The wafer is then placed in 1165 Remover (N-Methly-2-pyrrolidone (NMP)) at 65∘C to lift off the excess SiO2 and resist and then another clean bath of NMP at 65∘C for further liftoff. The wafer is then sonicated for 10 seconds to remove any final remaining resist and SiO2. The second layer of the ground strap process is exposed in the same way, using LOR3A and DUV photoresist, but this time the pattern lies over the existing SiO2 and extends further so that once developed, there is an exposed region of the niobium ground plane for the aluminum to contact. The wafer is baked again and developed, and the ground straps are then deposited by electron beam evaporation of aluminum (100 nm thick). The wafer is once again subject to NMP to remove the remaining resist and excess aluminum.

Once clean, the wafer is put through a light oxygen plasma resist strip before a bilayer resist stack of MMA/PMMA is spun for electron beam lithography to define the Josephson junctions. The Al-AlOx-Al junctions are written at 100 keV to form a standard double-angle evaporation airbridge pattern. Following development, there is a brief ion mill step before the first electrode is deposited by electron beam evaporation. The bottom (top) electrode is 40 (80) nm thick. Once the junctions are deposited, the wafer is covered in S1813 photoresist and then diced to (6.25 mm)2 chips. After the dicing, the aluminum metallization is lifted off and the chips are then cleaned with a UV/ozone process before measurement.

III. Device Layout

In order to allow for local flux-biasing of the different plaquettes and charge-biasing of the various superconducting islands, our device incorporates a series of on-chip bias lines, indicated in Fig. S3. The heart of the device contains a chain of three plaquettes, each with two small Josephson junctions (130 nm $\times$110 nm) and two junction-chain inductors (seventeen 140 nm $\times$1070 nm junctions in series). As discussed in the main paper, minimizing $C_{\rm isl}$ for each intermediate island between two adjacent plaquettes is critical for successful concatenation. Thus, ideally the four Josephson junctions in two adjacent plaquettes will all be located near the island between the plaquettes so that the junction electrode that is closest to the island will be as short as possible and contribute a minimal amount of excess capacitance to ground. However, in a chain of three plaquettes, this is only possible for one of the two intermediate islands. The other island will necessarily have to be connected to the two inductors for one of the plaquettes, and the capacitance to ground for these inductors will enhance the effective island capacitance. In addition, the 3-plaquette chain has dummy plaquettes at either end, which have the same geometry as the other plaquettes, but the small junctions and inductor-chain junctions are intentionally shorted out. The dummy plaquettes are included to symmetrize the geometry and minimize the inductive coupling of the on-chip flux-bias lines to the $LC$ mode of oscillation of the plaquette chain, sometimes referred to as the $M^{\prime}$ coupling, as defined in Ref. [4].

There are four on-chip flux-bias lines for controlling the flux bias to each of the three plaquettes, with the labeling as described in the main paper. Each flux-bias line has a coplanar geometry and splits into a T-shaped path adjacent to the plaquette chain, with the two ends of the T connected directly to the ground plane. In order to have a well-defined path for the return currents, and to suppress slot-line modes between different portions of the ground plane, we fabricated superconducting ground straps across each flux-bias line in multiple locations. In addition to the flux-bias lines, we also have three charge-bias lines for tuning the offset charge to the shunt capacitor electrode and each of the two intermediate islands between pairs of plaquettes. These charge-bias lines are isolated from ground, but also include similar ground straps to the flux-bias lines.

Our design also includes a pair of series dc SQUIDs between the plaquette chain and $C_{\rm sh}$ that could be used for gate operations in a future implementation of a protected qubit based on concatenated $\pi$-periodic plaquettes. For the experiments presented here, this SQUID switch, which has separate flux-bias lines from the plaquettes, was not used and the two loops of the SQUID switch were maintained at a flux bias of 0$\,{\rm mod}\,\Phi_{0}$ throughout the experiment. At this bias point, the SQUIDs behave primarily as superconducting shorts, although we must still account for the nonlinearity of the SQUID junctions in modeling the energy levels for our device.

The target shunt capacitance, $C_{\rm sh}\sim 1200\,{\rm fF}$ for our present device is rather large compared to more conventional superconducting qubits. Nonetheless, in the present experiment, we implemented $C_{\rm sh}$ with a planar superconducting Nb electrode with a small gap to the ground plane around the perimeter. For measuring our device, we have a coplanar waveguide (CPW) readout resonator with a fundamental resonance at 4.7 GHz. This is a 1/4-wave resonator with one end inductively coupled to a CPW feedline that is connected to our measurement circuitry; the other end of the resonator has a coupling capacitance $C_{\rm c}=44\,{\rm fF}$ to our device.

The majority of our device is patterned in Nb, including the ground plane, bias lines, readout resonator, and shunt capacacitor. All Josephson junctions are fabricated from a standard Al-AlOx-Al double-angle shadow-evaporation process. As an initial attempt at superconductor gap engineering for reducing quasiparticle poisoning of the plaquette chain, we include two patches of Al for suppressing the Nb gap underneath – one patch is at the joint between the plaquette chain and the ground plane; the other patch is between the plaquette chain and shunt capacitor.

IV. Device and Measurement Setup

Measurements are performed on a cryogen-free dilution refrigerator running at a temperature below 15 mK. The device chip is wire-bonded into a machined Al sample box that is mounted on a cold-finger attached to the mixing chamber stage and surrounded by a Cryoperm magnetic shield. The detailed configuration of our cabling, attenuation, filtering, and shielding inside the refrigerator, and the room-temperature electronics hardware for control and readout, is shown in Fig. S4.

V. Device Parameters

Establishing clear stabilizer behavior at double frustration requires plaquettes with a dominant $\pi$-periodic potential and large quantum fluctuations in the direction of constant $\varphi_{1}$+$\varphi_{2}$ in the space of common-mode phases across each plaquette. The $\pi$-periodicity comes from a dc SQUID consisting of two conventional Josephson junctions and a non-negligible loop inductance. We implement inductors in each plaquette with chains of large-area Josephson junctions, similar to typical fluxonium designs [3]. The inductive energy of the junction chain can be extracted with $E_{L}=(\Phi_{0}\Delta/4e)/R_{n}^{L}$, where $R_{n}^{L}$ is the junction chain resistance at room temperature. To have large quantum fluctuations in the direction of constant $\varphi_{1}+\varphi_{2}$ for effective hybridization between the two plaquettes, we need large $E_{C}$ and $E_{C}^{\rm isl}$ compared to the barrier height, which determines the coupling between the 00 and $\pi\pi$ wells and the 0$\pi$ and $\pi$0 wells. For our device, we target $E_{J}\sim 1.5\,\rm{K}$, $E_{L}\sim 1.5\,\rm{K}$ and $E_{C}\sim 3.5\,\rm{K}$ ($k_{B}=1$). For a junction with large $E_{J}$ and $E_{C}$, if the junction plasma frequency $\omega_{p}=\sqrt{2E_{J}E_{C}}/\hbar$ approaches $2\Delta$ of the junction electrodes, the junction acquires an extra capacitance from quasiparticles on either side of the junction. This specific electronic capacitance can be expressed as $C_{\rm elec}^{\rm sp}=3\hbar eJ_{C}/16\Delta^{2}$ [5], where $J_{C}$ is the critical current density of the junction and $\Delta$ is the superconducting gap. Our target $E_{J}$ is $\sim$1.5 K and junction area is  $110\,{\rm nm}\times 130\,{\rm nm}$, and the corresponding $J_{C}\sim$4 $\mu\rm{A}/{\mu\rm{m}}^{2}$ and $C_{\rm elec}\sim 0.3$ fF. Our estimated specific geometric capacitance is $\sim$50 fF/$\mu{\rm m}^{2}$, so $C_{\rm geo}\sim 0.7$ fF. The total capacitance of the junction is $C_{J}=C_{\rm elec}+C_{\rm geo}=1$ fF, and $E_{C}\sim$4 K. Junctions of this size are close to the lower limit where we can maintain reasonably small junction asymmetry with our fabrication. Thus, making smaller junctions to reduce $C_{\rm geo}$ is not practical. The simulated geometric charging energy of each intermediate island to ground $E_{C}^{\rm isl}=(2e)^{2}/2C_{\rm isl}$ is $\sim 4.6\,\rm{K}$ for the island between plaquettes 1 and 2, and $\sim 0.74\,\rm{K}$ for the island between plaquettes 2 and 3. The $E_{C}^{\rm isl}$ is significantly smaller for the island between plaquettes 2 and 3 because the inductor junction chains in plaquette 2 contribute to the capacitance to ground of the intermediate island. The intermediate island also has capacitance to ground through the junction capacitors [Fig. S10], thus $E_{C}$ of each of the four junctions in the two plaquettes reduces the total charging energy of the intermediate island below $E_{C}^{\rm isl}$. Thus, minimizing the capacitance of these junctions, hence targeting large $E_{C}$, is crucial for strong hybridization between plaquettes.

To estimate $E_{2}$ for this device, we model a circuit that embeds a plaquette in an rf SQUID, vary the flux across the rf SQUID loop, calculate the energy levels, and obtain the Fourier components for the lowest energy level. The $E_{2}$ value then corresponds to the Fourier component for the $\cos 2\varphi$ term. The extracted $E_{2}$ for this circuit is $\sim 0.05$ K. For effective concatenation, both $E_{C}$ and $E_{C}^{\rm isl}$ need to be large compared to $E_{2}$. For this target value of $E_{2}$, we require a rather large shunt capacitor with $C_{\rm sh}$=1200 fF in order to suppress single Cooper pair tunneling on/off the logical island. The charging energy for this shunt capacitor is $E_{C}^{\rm sh}=0.003$ K, and $E_{C}^{\rm sh}/E_{2}=0.06$, so the coupling between the even- and odd-parity states will be suppressed.

When there are small asymmetries in the circuit, particularly between the $E_{J}$ values of the two junctions in a plaquette, the even- and odd-parity states experience slightly different potentials and the computational states do not have their minimum gap exactly at frustration. Based on our room-temperature measurements of the resistance of nominally identical junctions, as well as low-temperature measurements of the critical current modulation depth for dc SQUIDs fabricated with junctions that are identical to those in the plaquettes, we expect $\alpha\sim 0.02$.

VI. flux scans: calibrating inductance matrix

With the double SQUID switch and three plaquettes, the device has a total of five flux-tunable loops. The six on-chip bias lines allow us to tune the flux in each of the loops independently, provided we account for the different mutual inductances between the various bias lines and flux-tunable loops. As described in the main manuscript, we map out the flux-bias parameter space by performing two-dimensional scans of the dispersive shift of the readout cavity for different pairs of flux-bias lines. When the flux through one of the loops approaches frustration, the resonance frequency of the cavity will decrease in response to the transitions of the plaquette circuit shifting to lower frequencies. We measure transmission through the feedline at a fixed cavity frequency near the resonance when one of the loops is frustrated. This results in high transmission when the plaquettes are away from frustration, while near frustration, we are driving on resonance and get low transmission.

Following a series of two-dimensional scans of various combinations of pairs of flux-bias lines, as in Fig. 2(a,b) in the main manuscript, we fit the slopes and periods of the frustration lines, then calculate the mutual inductance matrix from the following relation:

[TABLE]

where $\vec{\Phi}$ is a length-3 vector of the plaquette fluxes, $\vec{x}$ is a length-3 vector of the flux offsets to each plaquette at zero bias, $\vec{I}$ is a length-4 vector of the bias currents in the four plaquette flux-bias lines, and $\mathbf{\rm{L}}$ is a $3\times 4$ matrix of the mutual inductances. The flux offsets at zero bias are due to small background magnetic flux that gets trapped in place when the ground plane goes superconducting during the initial cooldown of the device. These flux offsets can be stable for weeks at a time, although small changes that necessitate recalibration can occur occasionally.

For our spectroscopy measurements at different degrees of frustration for the plaquettes, we must be able to control the fluxes to an accuracy better than 1 m$\Phi_{0}$. The resolution of these two-dimensional scans over multiple $\Phi_{0}$ is not sufficient to determine the flux offsets and mutuals to this level. To achieve this, we zoom in near one of the double frustration points with finer voltage steps on the two flux-bias lines [Fig. S5(a)]. Here, we see fine structure in the feedline transmission that is symmetric around frustration that comes from higher energy levels of the device crossing the cavity when a plaquette is tuned near frustration. We calculate the currents through the flux lines that the applied voltages create by analyzing the resistor network that is formed by the attenuators on the line. From the slopes and offsets relative to the symmetry points in these high-resolution flux scans, we extract the locations of double frustration with high accuracy; in addition, we refine the calculation of the flux periods and slopes in order to compute the mutual inductances with the required precision. Table 1 shows the extracted mutual inductance matrix from our measurements, along with a comparison to the inductance matrix obtained from simulating the layout with the numerical software package InductEx [6].

Using the experimental mutual inductance matrix and vector of offset fluxes, we can apply combinations of currents in the flux-bias lines to cancel out the various crosstalk fluxes and take steps in the pure flux direction for any plaquette or combination of plaquettes. Thus, we are able to scan along arbitrary vectors in the three-dimensional flux space for the three plaquettes. Figure  S5(b) is another high-resolution scan near plaquette (12) double frustration, but now the fluxes have been orthogonalized and the axes step through the pure fluxes through plaquettes 1 and 2 while the flux in plaquette 3 is maintained at unfrustration.

VII. Spectroscopy measurements

For spectroscopy measurements at single frustration, before each spectroscopy pulse, we initialize the circuit in the $\pi$ well by setting the plaquette flux bias to 0.1 $\Phi_{0}$ away from frustration, while maintaining the other two plaquettes at unfrustration. We then ramp the flux bias to each $\Delta\Phi$ coordinate on the flux axis using a gaussian edge with a 167 ns standard deviation, idle for 5 $\mu$s, then apply a 5 $\mu$s spectroscopy pulse to the $C_{\rm sh}$ charge-bias line followed by a 5 $\mu$s cavity readout pulse [Fig. S6(a)]. For measurements at double or triple frustration, we perform a similar initialization sequence, but in the $\pi\pi$ ($\pi\pi\pi$) well for double (triple) frustration.

We choose the 0.1 $\Phi_{0}$ initialization point so that there is a single well for the system to relax into. Initialization points further from frustration would also produce a deep single well, but the larger flux amplitude would enhance flux distortions on the trajectory back near frustration for the spectroscopy measurements. We determine the 30-$\mu$s initialization time following measurements where we vary this wait time. For wait times much less than 30 $\mu$s, we observe significant excitations out of the 0 well in subsequent spectroscopy, indicating that the system hasn’t fully reset into the $\pi$ well. Waiting longer than 30 $\mu$s doesn’t provide any further benefit for initializing the system.

We use a 167-ns gaussian edge pulse shape so that we are moving sufficiently fast to be non-adiabatic, at least for measurements near single frustration, but not so fast that there are significant fourier components near the qubit transition frequency that could cause spurious excitations. With this particular edge time, following initialization in the $\pi$ well, if we ramp to a flux past frustration where the $\pi$ well becomes metastable, we still observe transitions out of the $\pi$ well. The 5-$\mu$s idle time before the spectroscopy pulse is applied provides time for the flux to settle. Flux distortions are commonly observed in low-temperature measurements with fast flux pulses, with various possible causes, including impedance mismatches on the line and eddy currents in the normal copper traces in the sample box [7]. In principle, it is possible to measure these distortions and compensate for them by applying a pre-distortion to the pulse waveform [8]. For the measurements presented here, the short idle time is sufficient for the flux to settle, as determined by varying this time; for short idle times, the frequencies of the spectroscopy features drift with respect to flux, but by 5 $\mu$s these settle to an asymptotic level.

Following the spectroscopy pulse, for scans near plaquette 1 or 2 single frustration, the flux bias is then brought with a square pulse to a common readout point 10 m$\Phi_{0}$ to the right of frustration. For scans near plaquette 3 single frustration, we read out at the same flux point as the spectroscopy because we need to overlap the spectroscopy pulse with the readout pulse in this case. This is likely due to a shorter $T_{1}$ lifetime for the plasmon states for plaquette 3 compared to that for plaquettes 1 and 2. For all the spectroscopy scans, we plot the quadrature distance between the heterodyne measurement of transmission through the feedline at the cavity resonance for the readout flux bias with and without a 5-$\mu$s spectroscopy pulse.

Figure S6(b) shows an example of a spectroscopy measurement at single frustration for plaquette 2. The features that disperse gradually with flux correspond to the plasmon excitations within the $\pi$ well where the qubit is initialized, as described in Fig. 3 in the main paper; in addition to the 0-1, 0-2, and 0-3 transitions, we also observe transitions out of excited plasmon levels, indicating that the device is not fully initialized into the ground state of the $\pi$ well. In addition, we observe heavy fluxon transitions that disperse linearly with flux, and with a much steeper slope than the plasmons, that arise from transitions between levels in the $\pi$ and 0 wells. We observe qualitatively similar behavior for single frustration of plaquettes 1 and 3, as can be seen in the spectroscopy plots in Fig. S15.

For spectroscopy measurements at double and triple frustration, we add an extra step to stabilize the offset charge on the intermediate island(s) between the frustrated plaquettes. Details on this procedure are described in the next section.

VIII. Offset charge scans

As described in the main paper, Aharonov-Casher interference of the CW and CCW tunneling paths at various degrees of frustration results in a periodic modulation of the energy-level structure with respect to the offset charge on the $C_{\rm sh}$ island and the two intermediate islands between pairs of plaquettes. The modulation with offset charge on the $C_{\rm sh}$ island is difficult to observe directly in spectroscopy because the tunnel splittings for the low-lying levels are small. However, levels near the top of the barrier, which are also close to the readout cavity, exhibit a large modulation that leads to a significant periodic charge tuning of the readout cavity dispersive shift. The large physical footprint of the $C_{\rm sh}$ island results in a large effective charge sensing area, so that offset charge jumps occur on a timescale of a few minutes, as presented in Fig. 4(e) of the main paper.

For measuring the modulation with respect to offset charge on the intermediate islands between plaquettes, we perform the spectroscopy sequence at the various combinations of double and triple frustration while applying a spectroscopy pulse at the 0-1 transition frequency for each particular frustration point and scanning the two island charge biases [Fig. S7]. Each of these two charge lines, CB1 and CB2, couples to the intermediate island adjacent to it, but there is also non-negligible crosstalk to the other intermediate island. Thus, in general, we observe a periodic modulation due to the charge sensitivity of the relevant $\Delta_{\rm SA}^{(ij)}$ for the particular frustration point, and these modulation features have a slope in the two-dimensional charge-bias space that depends on the capacitance between each bias line and the intermediate island(s) between the particular pair(s) of frustrated plaquettes. At plaquette (12) double frustration, the modulation is significantly faster with respect to CB1 because the capacitive crosstalk between CB2 and the intermediate island between plaquettes 1 and 2 is relatively weak [Fig. S7(a)]. By contrast, at plaquette (23) double frustration, the modulation is faster with respect to both CB1 and CB2 since the junction-chain inductors of plaquette 2 contribute to the intermediate island capacitance between plaquettes 2 and 3 and enhance its capacitance to both charge-bias lines [Fig. S7(b)]. The modulation is even faster for plaquette (13) double frustration, since now the effective intermediate island includes all of plaquette 2, which is unfrustrated, thus enhancing the capacitance to both charge-bias lines [Fig. S7(c)]. At triple frustration, we observe a double charge modulation, with one set of nearly vertical features corresponding to the modulation with respect to the offset charge on the intermediate island between plaquettes 1 and 2, and faster, more diagonal features from the modulation with respect to the offset charge on the intermediate island between plaquettes 2 and 3 [Fig. S7(d)]. From the slope and period of these various modulation features, we can extract the capacitance matrix between the charge bias lines and the plaquette islands (Table 2). These capacitances agree reasonably well with the Q3D numerical simulations of our device geometry [9].

While the offset charge jumps on the $C_{\rm sh}$ island occur every few minutes, we expect the offset charge jumps on the intermediate islands between plaquettes to be less frequent because of the much smaller charge sensing areas [10, 11]. In order to monitor offset charge jumps on both intermediate islands nearly simultaneously, we first scan the offset charge bias to island 1 (between plaquettes 1 and 2) at plaquettes (12) double frustratation while plaquette 3 is biased 50 m$\Phi_{0}$ away from frustration; we then shift the fluxes slightly and scan the offset charge bias to island 2 (between plaquettes 2 and 3) at plaquette (23) double frustration while plaquette 1 is biased 50 m$\Phi_{0}$ from frustration. We alternate back and forth between these two scans repeatedly over 11 hours (Fig.  S8). As expected, given the smaller charge sensing areas, the offset charge on the intermediate islands jumps less frequently than on the $C_{\rm sh}$ island. The island between plaquettes 1 and 2, which has the smallest charge sensing area, is the most stable, with roughly 1 hour between large offset charge jumps.

Despite the relative stability of the offset charge on the intermediate islands, we still need to actively stabilize the charge for long spectroscopy scans vs. flux at double and triple frustration, such as Fig. 4(b,c) in the main paper, to correct for occasional offset charge jumps. Thus, approximately every twenty minutes we interrupt the spectroscopy sequence to run a one-dimensional scan of the relevant intermediate island offset charge bias(es) while applying a spectroscopy pulse at the corresponding 0-1 transition frequency; each scan takes about 30 seconds. We then fit the resulting modulation signal with a cosine function to determine the appropriate adjustment to the charge bias to apply to maintain a constant total intermediate island offset charge.

In addition to the random offset charge jumps due to charge dynamics in the qubit environment, for example, from the impact of high energy particles [11], the various islands of our device are also subject to quasiparticle (QP) poisoning when a QP tunnels on or off the island. Because of the large physical footprint for the $C_{\rm sh}$ island, there will be spurious antenna resonances, as described in Refs. [12, 13], at frequencies extending from above the Al superconducting gap to below 100 GHz that couple resonantly to stray photons in the device environment and generate QPs at the plaquette junctions. In Fig. S9, we show measurements of spectroscopy near the 0-3 transition as a function of frequency and intermediate island offset charge near plaquette (12) and (23) double frustration. In both cases, we observe a periodic charge modulation of the transition, but with two bands that are offset by $e$, indicating QP poisoning on the intermediate island between the plaquettes on a timescale faster than the spectroscopy measurement. Such QP poisoning will need to be significantly suppressed in future devices for the successful implementation of a protected qubit.

IX. Modeling of energy levels

Modeling multi-plaquette devices is challenging – our 3-plaquette chip with a SQUID switch has eleven phase degrees of freedom, and the size of the truncated Hilbert space is $\sim 2\times 10^{5}$. Instead of choosing generalized coordinates manually, we use the SuperQuantPackage [14] to model the energy level spectra of the devices.

The SuperQuantPackage software framework was developed by Andrey Klots with the supervision of Lev Ioffe. This package is capable of modeling the energy spectrum of superconducting circuits with arbitrary configurations of Josephson junctions, capacitors, and inductors. The original flux coordinates of the nodes undergo a linear transformation that splits them into two classes: oscillator-like coordinates for which we choose a harmonic oscillator basis, and charge coordinates that correspond to clusters of nodes with quantized net charge and for which a natural charge basis is used. This automatically diagonalizes the inductive and capacitive parts of the Hamiltonian. At the same time, Josephson terms of the Hamiltonian assume a relatively simple and sparse form. Automatic assignment of physically meaningful coordinates does not require labor-intensive manual symmetry analysis of each configuration of all studied complex circuits. Meanwhile, it allows for efficient diagonalization of the Hamiltonians and relatively quick numerical convergence.

Despite this optimization of the numerics, modeling the full circuit of our most complex devices [Fig. S10(a)] would require at least several months on the most powerful processors available to our research group. Thus, we must devise strategies for simplifying the modeled circuit to make the calculation practical. A plaquette or SQUID biased at unfrustration behaves like a superconducting inductive short with an effective shunt capacitance. For example, when modeling single frustration, we simplify the circuit to a single plaquette connected in series with a $LC$ resonator that represents the other unfrustrated plaquettes and SQUID-switch elements [Fig. S10(b)]. The $LC$ resonator inductance and capacitance are shown in Fig. S10(b) as $L_{\rm extra}$ and $C_{\rm extra}$. As shown in Fig. S10, each plaquette contains two arms, each having one Josephson junction in series with a linear inductor. The Josephson junction is characterized by $E_{J}$ and $E_{C}$, where $E_{J}$ is the average energy of the Josephson junctions and $E_{C}$ is the charging energy set by the junction capacitance. The linear inductor is characterized by $E_{L}$ and $E_{CL}$, where $E_{L}$ is the average inductive energy of the junction-chain inductor and $E_{CL}$ is the charging energy across the junction-chain inductor. From our fabrication uniformity tests, our nominally identical junctions exhibit a spread in $E_{J}$ of a few percent. We account for this asymmetry between the two junctions in a plaquette with the parameter $\alpha=(E_{JL}-E_{JR})/(E_{JL}+E_{JR})$, with $E_{JL}$ and $E_{JR}$ the Josephson energy of the left and right junction, respectively. Each arm has capacitance to ground, and this is characterized by $C_{\rm int}$. $C_{\rm sh}$ is the capacitance of the shunt. We introduce the parameter $L_{\rm factor}$ to account for variations in $L_{\rm extra}$ due to small flux offsets in the bias of the nominally unfrustrated plaquettes or SQUIDs. After this, we can input the circuit elements in the SuperQuantPackage.

The matrix is typically quite sparse, thus we can use the scipy.sparse.linalg.eigsh() function to find the eigenvalues and eigenvectors efficiently. This function only requires calculating the first few lowest eigenvalues, so the calculation speeds up. Typically, we only require the first $\sim$16-32 eigenvalues, but when we need to calculate the transitions involving the readout cavity, we need to calculate the first 40 eigenvalues.

The next step is finding the minimum number of states for each coordinate. We start by using three states for each cyclic coordinate and one state for each oscillator coordinate. We then vary the number of states from 1-20 for each coordinate, while tracking how the transition frequencies change. When the transition frequencies change by less than 5%, we choose the corresponding number of states for that particular coordinate for the next iteration. Using the new number of states, we repeat the same procedure until the process converges. In Fig. S11, we show the convergence for each of the coordinates as a function of the number of states for double frustration.

We model single frustration by considering a single plaquette connected in series with an $LC$ circuit [Fig. S10(b)]. For double frustration, we model the unfrustrated elements as a single $LC$ circuit in series with the two plaquettes that are modeled near double frustration [Fig. S10(c)].

A. Double plaquette flux dispersion vs. intermediate island capacitance

We modeled double frustration flux dispersion with different intermediate island capacitance, and we present some of these results in Fig. S12. When $C_{\rm isl}=1$ fF, the wavefunction is hybridized strongly between the $00$ and $\pi\pi$ wells, resulting in a large $\Delta_{\rm SA}\approx 2.9\,{\rm~{}GHz}$ and rather flat ground-state energy band [Fig. S12(a)]. The ground state energy band is relatively flat near frustration. When $C_{\rm isl}=$ 5 fF, the effective mass along this direction is larger, but there is still somewhat effective hybridization between the $00$ and $\pi\pi$ well [Fig. S12(b)]. In this case, $\Delta_{\rm SA}\approx 1.3~{}\text{GHz}$, and the energy band curvature is larger. When $C_{\rm isl}=$ 10 fF, the hybridization is significantly weaker [Fig. S12(c)]. $\Delta_{\rm SA}\approx 0.5~{}\text{GHz}$ and the antisymmetric level is now lower than the first excited plasmon state. Also, the energy bands have a nearly linear dispersion near double frustration. When $C_{\rm isl}=$ 50 fF, the effective mass is so large that all four wells are nearly independent with vanishing coupling between them [Fig. S12(d)]; $\Delta_{\rm SA}\approx 0~{}\text{GHz}$ and the flux dispersion near double frustration is essentially linear.

B. Double plaquette flux dispersion for different $\alpha$

As shown in Fig. S14, when $\alpha=0$, the even-parity energy level and odd-parity energy level experience the same potential, and the energy levels are nearly degenerate, with the minimum separation at exact frustration. When $\alpha\neq 0$ the barrier heights between the $00$, $0\pi$, $\pi 0$, and $\pi\pi$ wells become different. This leads to the even-parity energy levels and odd-parity energy levels experiencing different potentials, which causes the levels to cross on either side of frustration. The crossing points move further apart for larger $\alpha$.

C. Double plaquette charge dispersion at double frustration

With effective hybridization at double frustration, the splitting between symmetric and antisymmetric levels exhibits Aharonov-Casher interference, based on the offset charge bias $Q_{\rm isl}$ of the intermediate island between the frustrated plaquettes (Fig. S13). When $\alpha=0$, the symmetric/antisymmetric energy levels for both even and odd parity have $4e$ periodicity. When the symmetric and antisymmetric energy levels cross at $1e\,$mod$2e$, the gap closes [Fig. S13(b)], because the even- (odd-) parity wavefunctions both experience a $\cos 2\varphi$ potential. When $\alpha$ of plaquette 1 is 0.01 and $\alpha$ of plaquette 2 is 0.03 [Fig. S13(c)], both of the gaps for even- and odd-parity states are not fully closed at $1e\,$mod$2e$ due to incomplete destructive interference. When we bias the island charge at $0e\,$mod$2e$, the transition between the even- and odd-parity logical states is first-order insensitive to charge noise on the intermediate island.

D. Structureless plaquette model

Modeling a fully-structured three-plaquette circuit is computationally expensive, with 11 nodes [Fig. S10(d)], and each node requires several charge states. The matrix size is 189,000$\times$189,000, and thus requires $\sim$300 GB of RAM and takes weeks to calculate the energy levels, even with processors with 40 cores. As an alternative, we can use the stuctureless plaquette model to approximate the full-structure plaquette model.

We first connect one arm of the plaquette to form a loop, then vary the flux in this loop to obtain the potential of this arm. We next extract the Fourier components of this potential. In the structureless plaquette model, we replace the Josephson potential with the potential extracted from one arm of the plaquette. Thus, we do not need the linear inductors in the circuit, which reduces the number of effective nodes from 11 to 5, and the matrix size is now 7,000$\times$7,000. We add a renormalization factor in the junction capacitance to simulate the effect of higher internal levels. The need for this renormalization stems from the fact that the gap separating the potential energy of a single plaquette from the plaquette’s higher internal states is created by a massive Dirac-like Hamiltonian, which makes the flux particle effectively chiral. This chirality impedes tunneling, but cannot be implicitly accounted for in the structureless model. We mimic that by replacing this chiral property of a flux particle by increasing its effective mass by a renormalization factor that is fitted numerically by comparing structureless and full-structure plaquettes. For high internal excited state, the phase particle is less chiral and the renormalization factor is $\sim 1.5$ and low-lying internal excited states make the flux particle more chiral and yield a renormalization factor from $\sim 2$ to $\sim 3.5$ in the most unfavorable cases. We find the structureless plaquette model has a good agreement with the full-structure plaquette near frustration. We thus use the strucureless plaquette model for the remainder of this section.

E. Triple-frustration modeling

We model the triple-frustration flux dispersion for a simultaneous scan of the flux bias to each plaquette along a line from 000 to $\pi\pi\pi$ in Fig. S21. In Fig. S21(a), at tripled frustration, the computational states are a superposition of even-parity wells ($000$, $0\pi\pi$, $\pi 0\pi$, $\pi\pi 0$) and a superposition of odd-parity wells ($\pi\pi\pi$, $00\pi$, $0\pi 0$, $\pi 00$). In the ideal case of symmetric plaquettes and small flux offsets, the nearly-degenerate computational states are $(\left|000\right\rangle+\left|0\pi\pi\right\rangle+\left|\pi 0\pi\right\rangle+\left|\pi\pi 0\right\rangle)\pm(\left|\pi\pi\pi\right\rangle+\left|\pi 00\right\rangle+\left|0\pi 0\right\rangle+\left|00\pi\right\rangle)$ that correspond to even and odd charge parities. However, as we go away from the protected regime, the computational states turn into even ($\left|000\right\rangle+\left|0\pi\pi\right\rangle+\left|\pi 0\pi\right\rangle+\left|\pi\pi 0\right\rangle$) and odd ($\left|\pi\pi\pi\right\rangle+\left|\pi 00\right\rangle+\left|0\pi 0\right\rangle+\left|00\pi\right\rangle$) flux states. A few $\rm m\Phi_{0}$ away from triple frustration, the lowest energy level corresponds to the wavefunction localized only in the 000 or $\pi\pi\pi$ well. The second lowest levels correspond to a superposition of $00\pi$, $0\pi 0$, $\pi 00$ on the left and $0\pi\pi$, $\pi 0\pi$, $\pi\pi 0$ on the right. The energy levels for the 0 and 1 logical states both have negative curvature with respect to flux, so the flux dispersion of the 0-1 transition is flatter compared to double frustration, thus further enhancing the protection against flux noise.

X. Fitting of energy-level spectra

A. General fitting strategy

For extracting the center frequency of each transition at every flux point, the complexity of the level spectrum for our device makes it not practical to implement an automated routine that fits a standard curve to each feature in the spectroscopy data. Instead, we identify each transition feature in the data manually, then extract the maximum of the spectroscopy signal for each of these features. We then check the correspondence of the extracted transitions with our numerical model of the energy-level spectrum. We use the initial estimates for the various device parameters to calculate the energy levels, which generally match the data qualitatively. From this, we can identify most of the transitions. During the first round of fitting, where we adjust the circuit parameters about their estimated values to match with the transition data, described in detail below, we only fit the transitions that we have identified correctly in the initial step. Following this stage, the fitted energy levels typically match with the data quite well and we can identify more transitions. We then use the newly identified transitions to further refine the fit.

After extracting the plasmon, heavy fluxon, and light fluxon transitions, as well as the anticrossings from the spectroscopy flux- and charge-dependence data, we use the single (double) plaquette model for the energy-level spectra described in Sec. IX. to fit the single- (double-) frustration spectroscopy data. At single frustration, we fit $E_{J}$, $E_{C}$, $E_{L}$, $E_{CL}$, $C_{\rm sh}$, $\alpha$, and $L_{\rm factor}$ using the model shown in Fig. S10(b). We fix $C_{\rm int}$ to be 1 fF, which is estimated from numerical modeling with Q3D and a theoretical estimation of the effect of the junction chain capacitance to ground. We introduce the parameter $L_{\rm factor}$ to account for variations in $L_{\rm extra}$ due to small flux offsets in the bias of the nominally unfrustrated plaquettes or SQUIDs. At double frustration, we fit the same parameter set as in the single frustration case, but with the addition of $C_{\rm isl}$. Similarly to the single frustration case, we fix $C_{\rm int}$ to be 1 fF. We assume the two plaquettes at frustration share the same set of parameters, because the actual parameters between the two plaquettes are typically only different by a few percent based on our test structures during the device fabrication. This allows us to reduce the number of fitting parameters from 15 to 8, and thus makes the fitting more practical.

The cost function of our fitting procedure is $\sum_{n}W_{n}\Delta f_{n}^{2}$, where $\Delta f_{n}$ is the difference between the modeled and experimental frequencies for transition $n$, and $W_{n}$ is the weight that we assign to transition $n$. The goal of the fitting process is to minimize the cost function and find the parameter set that has less than a 10% difference between the modeled transitions and the experimental transitions. We use the scipy.optimize.minimize function in Python to do the fitting. We have 7 and 8 parameters for single- and double-frustration fitting, respectively. We find that the Nelder-Mead method performs better for this fitting than gradient descent methods in terms of avoiding local minima.

With this high-dimensional fit, we need to choose the initial parameters carefully. We use the initial $E_{J}$ and $E_{L}$ values calculated from the Ambegaokar-Baratoff relation using the on-chip test junction resistances. The initial $\alpha$ of the junctions is estimated from our test chips that each contain 6 identical junctions. As mentioned earlier, we define the charging energy as $E_{C}\equiv(2e)^{2}/2C$. The initial $E_{C}$ and $E_{CL}$ values are calculated from the the relevant junction areas measured with scanning electron microscopy with a total specific capacitance 70 fF/$\mu$m2. The initial $C_{\rm sh}$ is estimated from Q3D simulation. The initial $L_{\rm factor}$ is set to 1 because our unfrustrated plaquettes are nominally biased at unfrustration. We choose the initial simplex for the minimization so that it covers the possible range for each parameter, which is typically $\pm 5\%$ to $\pm 30\%$ of the initial values.

B. Single frustration fitting

We fit the single frustration data with our single plaquette model [Fig. S10(b)]. We put equal weight on different transitions by setting $W_{n}=1$ for the fitting. The fit runs on a computer with a 12-core processor and takes $\sim$1 day and 500 iterations to converge. In Fig. S15(c), we show the fitting of plaquette 2 single frustration. The red lines correspond to the fitted plasmon transitions, and the blue lines correspond to the fitted heavy fluxon transitions. The dotted lines are transitions out of the $\ket{0_{\pi}}$ state, corresponding to the 0 level of the $\pi$ well. The dash-dotted lines are transitions out of the $\ket{1_{\pi}}$ state. The dashed lines are transitions out of the $\ket{2_{\pi}}$ state. The transitions match with the data within 10% error. In Fig. S15(a), we show the modeled energy levels using the fitting parameters and indicate some example plasmon and fluxon transitions out of the $\ket{0_{\pi}}$, $\ket{1_{\pi}}$ and $\ket{2_{\pi}}$ states.

In Fig. S15(b-d), we show the single frustration fitting of plaquette 1, 2, 3 single frustration. Plaquette 1 behaves similarly to plaquette 2, and thus the fitted transitions and parameter values are similar. Plaquette 3 single frustration behaves somewhat differently, and the fitted energy levels and parameters differ by a larger amount compared to plaquettes 1 and 2. The fitted parameters are listed in Table 3 and they are within 20% of the parameters that we estimate from the design and fabrication tests, although a few of the parameters for plaquette 3 have a slightly larger discrepancy. The spectroscopy measurements at plaquette 3 single frustration are not as clean as for plaquettes 1 and 2 single frustration, thus potentially accounting for the larger variation with the estimated values.

C. Double frustration fitting

We next use the model in Fig. S10(c) to fit the double frustration data. Because the Hilbert space is $\sim$11 times larger at double frustration using this model, we use a 48-core processor to do the fitting. This process takes between 4-7 days and $\sim$300 iterations to converge. The anticrossings between different transitions are important features to fit because they determine the coupling between the computational states, so we put $\sim$20 times more weight for the regions in the spectroscopy data that exhibit significant anticrossings. We also simultaneously fit the corresponding charge modulation data and we only use the minimum and maximum of the charge modulation data for fitting. In order to compensate for the relatively small number of charge modulation data points, we put $\sim$50 times more weight for these features in the fitting.

In Fig. S16(b), we show the plaquette (12) double frustration data and fitted transitions. The red lines are the fitted plasmon transitions, the blue lines are the fitted heavy fluxon transitions, and the purple lines are the fitted light fluxon transitions. The dotted lines are transitions out of the $\ket{0_{ES}}$ state, where $E$ corresponds to the even-parity hybridized well between plaquettes 1 and 2, $S$ corresponds to the symmetric hybridized energy level of plaquettes 1 and 2, and 0 corresponds to the lowest energy level with these conditions. The dash-dotted lines are transitions out of the $\ket{1_{ES}}$ state and the dashed lines are transitions out of the $\ket{2_{ES}}$ state. The red solid line is $\ket{0_{EA}}\rightarrow\ket{1_{EA}}$, which is the transition out of the 0 state of the antisymmetric energy levels in the even-parity wells to the 1 state of the antisymmetric energy levels in the even-parity wells. We see this transition because there is fast quasiparticle poisoning on the intermediate islands that is faster than our measurement timescale. When we prepare the qubit in the $\ket{0_{ES}}$ state, the fast quasiparticle poisoning closes and opens the symmetric and antisymmetric gap randomly, which allows the system to occasionally transfer population from the $\ket{0_{ES}}$ to $\ket{0_{EA}}$ states, thus leaving population in the excited antisymmetric state. This results in the transition indicated by the solid red line in Fig. S16(b). In Fig. S9(b), we show the plaquette (12) double frustration charge modulation data of the $\ket{0_{ES}}\rightarrow\ket{3_{ES}}$ transition at 17 m$\Phi_{0}$. We can clearly see two quasiparticle bands, which we indicate by the red and blue dotted lines for the fitted transitions. We show the flux dependence of the fitted energy levels and corresponding transitions in Fig. S16(a). In Fig. S9(c), we show the plaquette (23) double frustration charge modulation data of the $\ket{0_{ES}}\rightarrow\ket{3_{ES}}$ transition at 11 m$\Phi_{0}$.

In Fig. S16(b-d), we show the fit results for the flux spectroscopy at plaquette (12), (23), and (13) double frustration. The $\Delta_{\rm SA}$ for (12), (23), and (13) double frustration are $\sim$ 2.7, 1.0, and 0.5 GHz, respectively, as expected for a decreasing $\Delta_{\rm SA}$ and progressively weaker hybridization for a larger intermediate island capacitance to ground. The fitted curves capture the transitions, anticrossings, and charge modulation to within 10%. The fitted parameters are shown in Table 4 and are in reasonable agreement with our estimated parameters.

In Fig. S16(c,d), the solid black lines between -20 to -10 m$\Phi_{0}$ in the plaquette (23) and (13) double frustration plots correspond to transitions involving the readout cavity: $\ket{0_{ES,n}}\rightarrow\ket{1_{EA,n-1}}$, where $n$ and $n-1$ indicate photon number in the cavity.

Near double frustration, as described in the schematic plots in Fig. 1(k) in the main paper, the parity of the initial state and the direction of the scan through the 2D flux-bias space determines which levels will disperse with respect to flux and which will be flat. For example, preparation in an even-parity state, say, $\pi\pi$, followed by a scan in the even-parity flux-bias direction, so, going from a flux bias where the $\pi\pi$ well is the global minimum to a flux where the 00 well is the global minimum, should result in the even-parity levels dispersing with flux while the odd-parity levels should be flat. This behavior should swap for the opposite parity of the initialized state and scan direction.

To examine this difference, we measure single-frequency spectroscopy scans in the 2D flux-bias space near plaquette (12) double frustration for both even- and odd-parity initializations [Fig. S17]. For each point in Fig. S17(a), we initialize in the $\pi\pi$ well near double frustration, then bring the flux bias to the point indicated in the 2D flux space, apply a spectroscopy pulse at 720 MHz, and then bring the flux bias to double frustration to read out [Fig. S17(a)]. The axes of the plot correspond to the even- and odd-parity flux directions. The black dots along the 00 to $\pi\pi$ direction indicate the points where we see transitions cross 720 MHz in our simulation of the spectrum [Fig S17(b)]. These points correspond to the transitions that are visible in the spectroscopy scan in Fig. 4(b) of the main paper. To identify the somewhat hyperbolic features away from the even-parity flux-bias axis in the odd-parity flux-bias direction, we consider the modeled level spectrum for a preparation in the even-parity well, with a scan along the odd-parity flux-bias direction [Fig S17(c)]. For flux-bias points of 10 m$\Phi_{0}$ or more away from frustration, there are many transitions out of the flat even-parity level up to the dispersing odd-parity levels. The transitions that match the 720 MHz spectroscopy frequency are labeled by red dots. If the odd-parity flux-bias scan is done for a nonzero offset along the even-parity flux-bias axis, the points that correspond to this same spectroscopy frequency shift in position, as shown again by the red dots.

If we now initialize in an odd-parity well, the $\pi 0$ well in this case, then repeat the same measurement process [Fig. S17(d)], we see that all of the features from the scan with even-parity initialization repeat, but are now rotated by 90 degrees. We again perform the same numerical modeling process for the level structure, but now with the $\pi 0$ initial state and the parity of the simulated flux-bias axes swapped [Fig. S17(e,f)]. Again, the features visible in this 2D spectroscopy scan match up with the various transitions that we identify in the modeled level spectrum.

D. Triple frustration modeling

We model the triple frustration data with the circuit in Fig. S10(d). The Hilbert space of this triple frustration model is $\sim$100 times larger than for the double frustration model. Modeling one flux point takes $\sim$10 days, so it is impossible to model several flux points and fit to the triple frustration data. We model the triple frustration energy levels by simultaneously modeling different flux points on four virtual machines in parallel, which have 48-, 24-, 12- and 12-core processors, respectively. In each of the virtual machines, we use the multiprocess function in Python to model different flux points in parallel so that it uses all the computational power in that virtual machine. Because fitting is not practical here, we use the fitted parameters from the single- and double-frustration modeling, and only adjust the $L_{\rm factor}$ to account for the SQUIDs not being at exact unfrustration. Figure 4(c) in the main paper shows the triple frustration flux spectroscopy data and modeled transitions. The red lines are the modeled plasmon transitions, the blue line is the modeled heavy fluxon transition, and the purple line is the modeled light fluxon transition. The dotted lines are transitions out of the 0 state of the symmetric energy levels in the even-parity wells. Based on our modeling, we find the plasmon and light fluxon transitions are in good agreement between our modeled curves and the triple frustration data. In general, we observe that the transitions at triple frustration have an even flatter dispersion with respect to flux compared to double frustration. As shown in Fig. S18, within $\sim\pm 2$ m$\Phi_{0}$ of triple frustration, the heavy fluxon transitions and plasmon transitions are difficult to distinguish, as the even-parity and odd-parity energy levels are likewise nearly on top of each other. Because of this, the dispersion of the various heavy fluxon transitions, $\ket{0_{E}}\rightarrow\ket{0_{O}},\ket{1_{O}},\ket{2_{O}},\ket{3_{O}}$, is comparable and rather flat.

Figure S19 shows the light fluxon transition plus cavity photon, which is visible near the edges of our triple frustration spectroscopy data, and the comparison with the modeled level spectrum. The difference between the levels indicated by arrows minus the readout cavity frequency matches these particular features in our data.

E. Fit parameter error estimation

In general, each type of transition has a different sensitivity to the various fitting parameters for the device. In order to estimate the uncertainty in each fitted parameter, we compute the energy level spectrum while varying each parameter one at a time and keeping the other parameters at their best-fit values. We thus find the range over which each parameter can be varied while keeping the transition frequencies within 10% of the measured values. In Fig. S20, we show the four transitions we choose to estimate the fitted errors. In Table 5(a), we list the fitted errors obtained with this method for the $\ket{0_{ES}}\rightarrow\ket{3_{ES}}$ plasmon transition at 20 m$\Phi_{0}$. From this table, we see the plasmon transition is sensitive to $E_{J}$, $E_{C}$, $E_{L}$, $C_{\rm sh}$ and $L_{\rm factor}$, so we are confident of the parameters extracted from fitting the plasmon transitions. In Table 5(b), we list the fitted errors for the $\ket{0_{ES}}\rightarrow\ket{2_{OS}}$ heavy fluxon transition at 6.5 m$\Phi_{0}$. The heavy fluxon transition is especially sensitive to $E_{J}$, $E_{C}$, while moderately sensitive to the rest of the parameters, so by fitting to the heavy fluxon transition, we have high confidence in the $E_{J}$, $E_{C}$ values, with moderate confidence in the rest of the parameters. In Table 5(c), we list the fitted errors for the $\ket{0_{ES}}\rightarrow\ket{0_{EA}}$ light fluxon transition at 0 m$\Phi_{0}$. The light fluxon transition is quite sensitive to $E_{J}$, $E_{C}$, $E_{L}$, $E_{CL}$, $C_{\rm isl}$ and $C_{\rm int}$, so these fitted parameters extracted from fitting the light fluxon transition have small errors. In Table 5(d), we list the fitted errors for the anticrossing between the $\ket{1_{ES}}\rightarrow\ket{4_{ES}}$ transition and the $\ket{1_{ES}}\rightarrow\ket{2_{OS}}$ transition at 12 m$\Phi_{0}$. The anticrossing is highly sensitive to $E_{J}$, $E_{C}$, $E_{L}$, so the errors for these fitted parameters by fitting the anticrossing are quite small. Our fitting method fits all four types of transitions, so we are confident that the parameters of the actual chip are within the intersection of the estimated errors extracted from these four types of transitions, as shown in Table 5(e).

XI. Prospects for implementing protected qubits

Our experiment presented in the main paper demonstrates successful concatenation of $\pi$-periodic plaquettes through measurements of large $\Delta_{\rm SA}^{(ij)}$ and offset-charge tuning, characteristic of strong hybridization of frustrated plaquettes. In order to implement a protected qubit based on this approach, we need a device that robustly maintains nearly degenerate computational states while pushing the remaining fluxon and plasmon energy levels significantly upwards. This means that for qubit design, we want to gain control over the three major energy scales: (1) Heavy fluxon gap $\varDelta_{\rm EO}$, that is also the computational gap that protects from $T_{1}$-processes and charge noise, needs to be kept small; (2) Light fluxon gap $\varDelta_{\rm SA}$, that determines effectiveness of concatenation and protects from flux noise and $T_{2}$-processes, needs to be kept large; (3) Plasmon gap $\omega_{\text{pl}}$ that determines energies of unwanted low-lying plasmon states that deteriorate initialization fidelity and facilitate thermal excitations, needs to be kept high. For determining these energy scales, there are three primary device parameters that can be adjusted in the design: Shunt capacitance $C_{\rm sh}$; Effective intermediate island charging energy $E_{C}^{\rm isl}$, that is determined by both the geometric ground capacitance of the intermediate island and the effective screened capacitances of the Josephson junctions that also indirectly couple the intermediate island to ground; $E_{2}$ of an individual plaquette, that is mainly determined by a combination of $E_{L},E_{J},E_{C}$. The table below demonstrates how the three energy gaps change as the described device parameters increase. The entries of the table describe whether or not the gaps change in the desired direction, and whether they change exponentially (exp) or (sub)polynomially (poly). Each entry is also accompanied by a brief explanation.

[TABLE]

Note that each of these device parameters improves some energy scales, while deteriorating others. Thus, the goal of optimization is not to just maximize some circuit parameters, but to compensate the negative effects by tuning other parameters. For example, increasing $E_{C}^{\rm isl}$ significantly improves concatenation (light fluxon modes), but negatively, albeit weakly, affects the heavy fluxon and plasmon modes. Therefore, it is always beneficial to increase $E_{C}^{\rm isl}$, while keeping in mind the presence of mild adverse effects. On the other hand, increasing $C_{\rm sh}$ and $E_{2}$ is exponentially beneficial for the computational gap, but (sub)polynomially adversely affects concatenation (light fluxon) and the plasmon modes. It is encouraging that adverse effects of increasing $C_{\rm sh}$ and $E_{2}$ are smaller than the benefits. However, these adverse effects are not negligible and need to be treated very carefully. In addition to these device parameters that can be controlled through design, there are also uncontrolled ones, such as the flux noise amplitude and junction asymmetry $\alpha$. However, improved concatenation protects the device against these types of asymmetries. This means that decreasing $\alpha$ would also increase $T_{2}$ polynomially – as a power of the number of plaquettes minus one.

One set of parameters that provides a compromise between the competing circuit parameters needed for achieving a qubit with excellent coherence is the following: $E_{J}=3$ K, $E_{L}=2$ K, $E_{C}=5$ K, $C_{\rm sh}=1000$ fF, $C_{\rm isl1}=1$ fF, $C_{\rm isl2}=3$ fF, $\alpha=0.01$. Unfortunately, achieving these values of $E_{J}$ and $E_{C}$ with conventional Al electrodes is not possible due to the small superconducting gap and the electronic capacitance arising when the junction plasma frequency approaches the gap [5]. Thus, protected qubits incorporating this stabilizer mechanism will need junctions fabricated from a larger gap superconductor. Reducing $\alpha$ from $\sim$0.03 in our experimental device to 0.01 will be challenging, but should be achievable with improved control over the electron-beam lithography and thin-film growth for forming the junctions. Reducing the intermediate island capacitances from our experimental device in the main paper to the values listed above requires minimizing the size of the intermediate islands as well as the length of the junction chain inductors, but the levels listed here are feasible. An alternate approach for satisfying the competing constraints involving parallel chains of concatenated plaquettes will be discussed at the end of this section.

In addition to optimizing the device parameters as described above, the device layout must have weaker radiative coupling to parasitic high-frequency modes to avoid spurious excitations and fast quasiparticle poisoning of the intermediate islands. Achieving this in practice is challenging because of the large $C_{\rm sh}$. On our present device, the planar capacitor layout for $C_{\rm sh}$ required an exceedingly large physical footprint, resulting in significant low-frequency parasitic modes that couple strongly to the qubit environment. For a future protected qubit implementation, switching to a parallel-plate capacitor for $C_{\rm sh}$ would allow for a much more compact structure, thus pushing parasitic antenna modes to much higher frequencies. For a conventional superconducting qubit, such as a transmon, the microwave losses for the deposited amorphous dielectric layer in such a capacitor would result in quite poor $T_{1}$ performance. However, for a qubit formed from concatenated $\pi$-periodic plaquettes with a large $C_{\rm sh}$, the transition matrix elements between the logical qubit states are quite small. When combined with the vanishingly small value for $\Delta_{\rm EO}$, this results is an extremely long $T_{1}$, even with the large loss tangent for the deposited dielectric of $C_{\rm sh}$, as will be discussed below in our treatment of the projected coherence times.

In order to estimate the coherence for a protected qubit based on this design with the parameters described earlier in this section, we first compute the level spectrum for such a device at triple frustration (Fig. S21). With these $E_{J}$, $E_{L}$, and $C_{\rm sh}$ values, the separation between the logical qubit states is quite flat with respect to flux, with a minimum separation ($\Delta_{\rm EO}$) of $15$ kHz. The lowest plasmon transition is also pushed up to nearly 1 GHz above the ground-state doublet. At plaquette (12) and (23) double frustration, the symmetric-antisymmetric gaps are $\Delta_{\rm SA}^{(12)}=1.5$ GHz and $\Delta_{\rm SA}^{(23)}=0.9$ GHz, and so are comparable with the lowest plasmon transition.

Upon computing the level spectrum near triple frustration for this set of parameters, we can estimate the qubit coherence. For obtaining $T_{1}$, we first calculate the transition matrix element for the charge operator $\hat{N}$ between the even- and odd-parity ground states ($\psi_{\rm even},\psi_{\rm odd}$, respectively) at triple frustration for the above parameters: $\absolutevalue{\langle\psi_{\rm odd}\mid\hat{N}\mid\psi_{\rm even}\rangle}\sim 1\times 10^{-5}$. $T_{1}$ can then be expressed as

[TABLE]

Assuming we fabricate a parallel-plate shunt capacitor using electron-beam evaporated SiO2, which is compatible with our current device fabrication process, the corresponding loss tangent is $\sim$1/300 [15], resulting in $T_{1}\sim 1\times 10^{8}$ s. In general, for our architecture, the $T_{1}$ budget is quite high because it is determined by exponentially weak and easily suppressable interwell tunneling. The much bigger challenge is creating flat bands to enhance $T_{2}$.

For computing $T_{2}$ for the device described above, we consider the flux dispersion of the energy bands from our numerical modeling. Exactly at triple frustration, the slope with respect to flux bias vanishes and the curvature is $\sim 12$ MHz/${\rm m}\Phi_{0}^{2}$. Assuming the electronics used for supplying the flux-bias currents is capable of sufficiently fine steps, the flux resolution of our biasing will ultimately be limited by the rms flux noise level $\Phi_{\rm noise}$. We calculate the $1/f$-noise-limited $T_{2}$ with the formula:

[TABLE]

where the flux noise is assumed to follow the power spectrum $A_{\Phi}^{2}/f$, we consider a Ramsey pulse performed immediately after calibration. Taking a flux-noise power spectral density of $2\,\mu\Phi_{0}/\sqrt{\rm{Hz}}$ at 1 Hz and integrating over 10 decades yields an rms flux noise level of 10 $\mu\Phi_{0}$. Accounting for these various factors, we compute a flux-noise limited dephasing time $T_{2}^{\Phi}$ for single frustration, plaquette (12) double frustration, plaquettte (23) double frustration, and triple frustration: $1.7~{}\mu$s, $340~{}\mu$s, $190~{}\mu$s, $6.0\times 10^{4}$ $\mu$s, respectively.

We also consider dephasing due to charge noise at triple frustration by computing the level of fluctuations in the splitting between the computational levels. We estimate the lower limit of $T_{2}^{Q}$ with

[TABLE]

where we assume a charge noise power spectrum $A_{Q}^{2}/f$ and we again, like in the case of flux noise, assume that the Ramsey pulse is performed immediately after calibration. For a charge-noise spectral density of $2\times 10^{-2}e/\sqrt{{\rm Hz}}$ at 1 Hz [16, 10], we estimate a charge-noise limited dephasing time, $T_{2}^{Q}>10^{7}$ $\mu$s. This timescale decreases for larger $\Delta_{\rm EO}$, for example, with a smaller $C_{\rm sh}$. Thus, reducing $C_{\rm sh}$ in an attempt to increase $\omega_{\rm pl}$, will eventually cause dephasing due to charge noise to become non-negligible as $\Delta_{\rm EO}$ increases.

In addition to the optimization of our multi-plaquette chain discussed above, there are alternative device geometries that can mitigate the adverse effects of decreasing $C_{\rm sh}$. For example, connecting multiple plaquette chains in parallel, as shown in Fig. S22, results in a linear increase in the overall system $E_{2}$, thus raising the plasmon modes without the need for increasing $E_{2}$ of the individual plaquettes, which is helpful for maintaining strong concatenation. This allows for the use of a smaller $C_{\rm sh}$ for raising the plasmon transition frequencies without sacrificing an increase in $\Delta_{\rm EO}$. Of course, such an approach with two 3-plaquette chains will require twice as many flux bias controls compared to a single chain, plus one more flux bias for the loop formed by the two chains.

XII. $\Lambda$ calculation

As described in the main paper, $\Lambda$ characterizes the rate at which the logical error decreases with system size. $\Lambda$ can be expressed either as the ratio of $T_{2}$ between the larger and smaller system sizes, or, in the specific case of a qubit based on concatenated plaquettes, can also be determined by $\Delta_{\rm SA}$ and $h_{Z}$, the scale of dephasing fluctuations for a single plaquette. In this section, we show that these two approaches are essentially equivalent. In addition, we compute $\Lambda$ for different degrees of frustration for a qubit with the optimal parameters described in the previous section. For a given level of flux noise, $T_{2}^{\Phi}$ will be inversely proportional to the slope of the qubit transition energy with respect to flux. Thus, computing $\Lambda$ from the ratio of the $T_{2}$ values will be equivalent to the inverse ratio of the energy band slopes for the two different degrees of frustration.

For the semi-quantitative description of the system at double frustration presented in the main paper, we can treat each plaquette as a spin-1/2 particle with an $XX$ stabilizer term. The energy of each plaquette is $H_{i}=-(k\Delta\Phi_{i}/{2})Z_{i}$, where $k$ is the slope of the transition energy with respect to flux, and $\Delta\Phi_{i}$ is the flux offset from frustration for plaquette $i$. In this model, the double-plaquette Hamiltonian is then:

[TABLE]

Here, we assume the two plaquettes have the same $k$. The eigenvalues are

[TABLE]

The double-plaquette circuit typically has $k\sim 300$ MHz$/{\rm m}\Phi_{0}$ and $\Delta_{\rm SA}\sim 1$ GHz. Near single frustration for one of the plaquettes, the other plaquette will have $\Delta\Phi$ near 0.5 $\Phi_{0}$, thus $(k\Delta\Phi_{1}/2+k\Delta\Phi_{2}/2)^{2}\gg(\Delta_{\rm SA}/2)^{2}$. In this limit, we obtain:

[TABLE]

Thus, the transition energy between the lowest two energy levels is $\Delta E=k\Delta\Phi$ and the corresponding slope is

[TABLE]

where S indicates single frustration.

Near double frustration, the slope is highest in the $\Delta\Phi_{1}=\Delta\Phi_{2}$ or $\Delta\Phi_{1}=-\Delta\Phi_{2}$ directions. Taking $\Delta\Phi_{1}=\Delta\Phi_{2}=\Delta\Phi$, the energy levels are

[TABLE]

The transition energy between the lowest two energy levels is $\Delta_{\rm SA}/2-\sqrt{(k\Delta\Phi)^{2}+(\Delta_{\rm SA}/2)^{2}}$, and the corresponding slope with respect to $\Delta\Phi$ is $-k^{2}\Delta\Phi/\sqrt{(k\Delta\Phi)^{2}+(\Delta_{\rm SA}/2)^{2}}$. For a $1/f$ flux noise level of 2 $\mu\Phi_{0}/\sqrt{{\rm Hz}}$ at 1 Hz, the rms flux noise amplitude is $\Phi_{\rm noise}\sim$10 $\mu\Phi_{0}$, and we take this to be the minimum precision with which one can set the flux bias. For typical double-plaquette parameters, $k\sim 300$ MHz$/{\rm m}\Phi_{0}$ and $\Delta_{\rm SA}\sim 1$ GHz. In this limit, $\Delta_{\rm SA}/2\gg|(k\Delta\Phi_{1}/2+k\Delta\Phi_{2}/2)|$, and the slope of the transition energy with respect to flux at double frustration is

[TABLE]

where D indicates double frustration. $\Lambda_{\rm DS}$, the decrease in the logical error rate upon going from single to double frustration, is then

[TABLE]

Thus, in this approximation, $\Lambda$ defined as the ratio of $T_{2}$ at double and single frustration is equal to the ratio of $\Delta_{\rm SA}$ to $2h_{Z}$, where the factor of 2 accounts for the two plaquettes.

For triple frustration, the Hamiltonian for the coupled spin-1/2 model of a 3-plaquette chain is given by

[TABLE]

assuming all three plaquettes have the same single-frustration slope $k$. Finding the eigenvalues of this Hamiltonian and expanding them to third order in the limit of $\Delta_{\rm SA}^{(12)}/2\gg|k\Delta\Phi_{1}/2+k\Delta\Phi_{2}/2|$ and $\Delta_{\rm SA}^{(23)}/2\gg|k\Delta\Phi_{2}/2+k\Delta\Phi_{3}/2|$, we find the transition energy of the lowest two levels is $\Delta E=(k\Delta\Phi)^{3}/\Delta_{\rm SA}^{(12)}\Delta_{\rm SA}^{(23)}$. The corresponding slope with respect to flux is

[TABLE]

where $T$ indicates triple frustration.

We can then express $\Lambda_{\rm TD}$, the reduction in logical error rate upon going from double to triple frustration by taking the ratio of the expressions in Eq. (S11) and Eq. (S13). Thus, for plaquette (12) double frustration, we obtain

[TABLE]

and for plaquette (23) double frustration, we obtain

[TABLE]

Then, $\Lambda_{\rm TS}$, the reduction in logical error rate upon going from single to triple frustration, is given by

[TABLE]

In Table 6, we list the various energy scales and energy band parameters for the device with optimal parameters, described in Sec. XI., as well as the $T_{2}^{\Phi}$ values calculated from the modeled energy levels. We also list the $\Lambda$ values for the various degrees of frustration, calculated both from the ratio of $T_{2}^{\Phi}$ values from the modeled energy levels and from $\Delta_{\rm SA}$ and $h_{Z}$. We see that both approaches for computing $\Lambda$ are in agreement, consistent with our analysis described above.