Superconducting circuit protected by two-Cooper-pair tunneling

We present a protected superconducting qubit based on an effective circuit element that only allows pairs of Cooper pairs to tunnel. These dynamics give rise to a nearly degenerate ground state manifold indexed by the parity of tunneled Cooper pairs. We show that, when the circuit element is shunted by a large capacitance, this manifold can be used as a logical qubit that we expect to be insensitive to multiple relaxation and dephasing mechanisms.


A. Motivation
Superconducting circuits are widely recognized as a powerful potential platform for quantum computation and now stand at the frontier of quantum error correction [1]. Future progress will likely stem from two complementary strategies: (i) active error correction characterized by measurement-based [2][3][4] and autonomous [5][6][7][8] stabilization, and (ii) passive error correction characterized by protected qubits [9, and references therein]. We address strategy (ii) in this article by designing an experimentally accessible protected qubit.
The transmon qubit [10] has proven to be a remarkably successful prototype of a protected qubit. Its circuit (depicted in Fig. 1a) contains a Josephson junction whose potential energy is U = −E J cos ϕ, with E J being the tunneling energy and ϕ being the superconducting phase across the junction. The large shunt capacitance imposes a large effective mass for the analogous "particle in a box," confining the low-energy wavefunctions near ϕ = 0, the only minimum for ϕ ∈ (−π, π). This confinement suppresses the susceptibility of the qubit to offset charge noise and renders the energy spectrum approximately harmonic with level spacing much smaller than E J (see Fig. 1b). On the other hand, circuit elements with degenerate phase states that only allow tunneling of pairs of Cooper pairs, meaning their potential energy is U = −E J cos 2ϕ, have been developed in recent years as a building block for topologically protected qubits [11,12]. In this article, we propose a transmon-like qubit with additional protection from environmental noise by combining the large shunt capacitance of the transmon with such a cos 2ϕ circuit element (the cross-hatched box in Fig. 1c). In this case, the wavefunctions are localized near ϕ = 0, π (see Fig. 1d), resulting in a nearly degenerate harmonic level arrangement. While the detrimental effects of offset charge noise are similarly suppressed in this circuit, sensitivity of the qubit to other decoherence * Current Address: QUANTIC Team, INRIA Paris, 2 Rue Simone IFF, 75012 Paris, France; william.smith@inria.fr † michel.devoret@yale.edu mechanisms is also reduced, owing to the conservation of Cooper pair number parity.
In this article, we introduce a few-body transmontype qubit where the charge carriers are exclusively pairs of Cooper pairs. Our central result is that there exists an experimentally attainable parameter regime for which conservative predictions of relaxation and dephasing times exceed 1 ms, i.e. an order of magnitude higher than those of typical transmons, given the same environmental noise [13,14]. In the remainder of Sec. I, we describe a toy model for the protected qubit. We proceed by analytically and numerically examining the Hamiltonian for the full superconducting circuit in Sec. II. Our attention in Sec. III then turns to properties of the ground state manifold, which we envision using as a protected qubit. A brief discussion about the concept of protection and examples of protected qubits, as well as our perspectives on readout and control, follows in Sec. IV. Finally, we summarize our results in Sec. V.

B. cos 2ϕ element
We first examine the advantages of the ideal circuit in Fig. 1c as a protected qubit. This circuit can be viewed as a Josephson-junction-like element (the cross-hatched box) shunted by a capacitance. Pairs of Cooper pairs are the only charge excitations permitted to tunnel through this element [11]. In the Cooper pair number basis, the potential energy assumes the form where E J is the effective tunneling energy of the process. This expression follows from the conjugacy relation [ϕ, N ] = i, where N is the number of Cooper pairs that have tunneled. The invariance of the potential under translations in ϕ by multiples of π implies that halffluxons are able to traverse the element. The shunt capacitance and other charging effects produce a quadratic kinetic energy, yielding the Hamiltonian where E C is the charging energy and N g is the offset charge. This offset charge has been introduced due to the periodicity of the Hamiltonian in ϕ, which reflects the presence of a superconducting island in the circuit (as colored in Fig. 1c).
Since the circuit element only allows pairs of Cooper pairs to tunnel, the parity of the number of Cooper pairs that have tunneled is preserved under the action of the Hamiltonian. This property leads to two nearly degenerate ground states |+ and |− , which only consist of even and odd Cooper pair number states, respectively [12]. Since these states have no overlap in charge space (equivalently, they have opposite periodicity in phase space-see Fig. 1d), we have −|O|+ ≈ 0 for any sufficiently local operator O. Furthermore, the states [15] | = 1 √ 2 (|+ ± |− ) are respectively localized near ϕ = 0, π (see Fig. 1d). Because these states have suppressed overlap in phase space for large E J /E C (i.e. they are roughly inversely periodic in charge space), we have |O| ≈ 0 for similarly local O. These are precisely the conditions for simultaneously suppressing spurious transitions and phase changes between the states [16] [resembling a Gottesman-Kitaev-Preskill (GKP) encoding on a circle [17]].
More concretely, the ground state splitting obeys for large E J /E C (see App. A). The two ground state energies oscillate out of phase with one another in N g .
Moreover, this shows that the splitting, as well as the charge dispersion, is exponentially suppressed in E J /E C . Thus, the role of the shunt capacitance is to decrease the charging energy E C and hence mitigate offset charge noise, much like in the transmon qubit [10].

A. Hamiltonian
We now detail the superconducting circuit necessary to practically implement the sought-after cos 2ϕ Josephson element. This circuit, depicted in Fig. 2a, is composed of two identical arms, each containing a Josephson junction in series with superinductances [18,19], arranged in parallel [20,21] and shunted by a large capacitance. These superinductances are split in half and placed on either side of the respective Josephson junctions to avoid large capacitances to ground. Kirchoff's current law allows us to treat the phases across the superinductances in series as equal. When the external magnetic flux threading the inductive loop reaches half of a flux quantum, i.e. when  2. (a) Reduced electrical circuit for the physical protected qubit. When ϕext = π, the two Josephson junctions and superinductances collectively behave as the cross-hatched element. The superconducting island is indicated by color. (b) Contour plot of the potential energy U in Eq. 2 in the ϕ1ϕ2-plane at ϕext = π for θ = 0. The numerically computed instanton trajectory between adjacent potential minima is overlaid in black. Importantly, this trajectory closely resembles a sequence of straight lines. ϕ ext = π, the loop approximates the cross-hatched box in Fig. 1c. At this particular bias point, a Cooper pair can only tunnel through one of the Josephson junctions if it is accompanied by another Cooper pair tunneling through the other junction (in either direction). Conversely, a fluxon traversing a single Josephson junction corresponds to a half-fluxon traversing the whole element.
We consider the symmetric and antisymmetric combinations of superconducting phase coordinates and their conjugate charges {n, N, η} in numbers of Cooper pairs. Note that the prefactor in the definition of ϕ is chosen to bring the coordinate into agreement with the phase drop across the inductive loop in the limit that θ vanishes. The Hamiltonian reads where C and J are the single junction charging and tunneling energies, 2 L is the inductive energy of each superinductance, x ≡ C J /C shunt is the ratio of the junction capacitance to the shunt capacitance, and N g is the offset charge on the superconducting island (see Fig. 2a). From this expression, it is clear that this circuit has three strongly coupled modes. The φ mode is flux dependent and is strongly and nonlinearly coupled, via the Josephson junctions, to the ϕ mode. The ϕ mode is offset charge dependent and strongly but linearly capacitively coupled to the θ mode. Our analysis is concerned with the fluxonium-like [18] parameter regime J L and J C . In particular, the parameters chosen for numerical simulations are listed in Tab. I.

B. Semiclassical theory
In order to gain insight into the structure of the Hamiltonian in Eq. 2, we briefly revisit its potential energy U in the ϕ 1 ϕ 2 -plane, which is plotted in Fig. 2b for θ = 0. The cosine terms in the potential form a two-dimensional "egg carton" of wells. The minimum of the quadratic term in the potential occurs at ϕ 1 + ϕ 2 = ϕ ext , which generally falls between adjacent diagonal ridges of cosine wells. At the special value of ϕ ext = π, these two ridges are degenerate. Near this value of the external flux, we consider the path of the system between neighboring potential minima by numerically solving for a two-dimensional instanton [22] trajectory (see App. B).
We then constrain the system to this maximally probable tunneling path. For the parameters mentioned above, this path is well-described by where z ≡ L / J is a small parameter. Plugging this expression into the Hamiltonian in Eq. 2, approximating the resulting one-dimensional potential by the first few terms in its Fourier series, and Taylor expanding about (4) to leading order. Here, φ ext = ϕ ext − 4π round ϕext 4π and we have discarded terms higher than the second harmonic or O( L ) (see App. B for additional terms). This treatment exposes the 'cos 2ϕ nature' of the potential at ϕ ext = π, where the cos ϕ term vanishes. By comparison to Eq. 1, we see that the added complication is that the ϕ mode is strongly coupled to the θ mode. The resulting hybridization is a central ingredient to understanding properties of the system beyond the ground state manifold (see Sec. II D).
We comment that this approximation neglects quantum fluctuations that are perpendicular to the path in Eq. 3. This is consequently a semiclassical approximation: we have minimized the energy of the system with respect to the dynamical coordinate orthogonal to the trajectory [22]. Moreover, from Fig. 2b, it is clear that the approximation we have made is that fluxons traverse a single Josephson junction at a time.

C. Energy spectrum
From numerical diagonalization of Eq. 2 (see App. C), we obtain the dependence of the energy levels on external flux as shown in Fig. 3a [23]. At ϕ ext = π (the dashed line in Fig. 3a), the spectrum resembles a doubled harmonic oscillator with energy √ 16x L C . Once ϕ ext deviates from π, half of the energy levels increase in energy linearly with slope ∼ 32 3 L [24]. The other half of the energy levels form a flux-independent harmonic ladder.
We can understand this level structure as that of two emergent modes. The first mode is flux-dependent and its excitations correspond to the number of magnetic flux vortices, or fluxons, enclosed by the inductive loop in 2a. Each plasmon involves the two superinductances (energy 2 L in parallel) and the shunt capacitance (energy x C ), and hence has energy √ 16x L C . Hereafter, we refer to these modes as the "fluxon mode" and the "plasmon mode," respectively [25]. Additionally, we assign the labels |m • to the lowest-energy states, where m denotes the number of plasmons and • (or • ) denotes the presence (or absence) of a fluxon excitation relative to the ground state [26]. Note that the fluxon placeholder indices are defined by for ϕ ext mod 2π < π −/+, for ϕ ext mod 2π = π / , for ϕ ext mod 2π > π , where |m± = |m ⊗ 1 √ 2 (| ±| ) (as in Sec. I B) and the index represents the persistent current direction. This labeling serves the purpose of assigning quantum numbers consistently for all external flux values (see the colors in Fig. 3a), except the particular case ϕ ext mod 2π = 0 that we do not focus on.
To explain the energy level structure at ϕ ext = π, we factor Eq. 4 into the form [27] The parity sectors are governed by the Hamiltonians In the above, k + = 0 and k − = 1 whileφ andÑ should be viewed as conjugate operators corresponding to pairs of Cooper pairs. Note that this expression is an exact reformulation of Eq. 4. In the limit that J C , we may Taylor expand the potential aboutφ = 0 and retain the first few terms. This step discards the effects of the offset charge, rendering H + and H − identical. We then decompose the linear part of the Hamiltonian into normal modes, which yields to leading order in z. The corresponding transformation is given by Eq. 5 reveals two weakly anharmonic modes: the plasmon mode at the low frequency of √ 16x L C /h and a junction self-resonant mode at the high frequency of √ 8 J C /h. These modes are coupled by a quartic nonlinearity, which has the primary effect of inducing a Kerr shift on the junction self-resonant mode. At frequencies lower than √ 8 J C /h, the energy level structure is that of a two-fold degenerate harmonic oscillator, in agreement with the simulated energy spectrum at ϕ ext = π (see the dashed line in Fig. 3a).

D. Wavefunctions
Calculated charge wavefunctions N |ψ and phase wavefunctions ϕ, φ|ψ , obtained from numerical diag- The charge wavefunctions are grid states with Fockstate envelopes [17,28]. For fluxon excitation index +/−, these grid states are superpositions of even/odd Cooper pair number states [29]. Additionally, m corresponds to the order of the Fock state envelope. Note that a logical qubit encoded in |0+ and |0− is protected from spurious transitions except those mediated by operators that flip Cooper pair parity.
On the other hand, the phase wavefunctions are approximately Fock states localized within the potential energy wells (see Sec. II B) [30]. The fluxon index +/− denotes whether the state |m± is a symmetric (bonding) or antisymmetric (antibonding) superposition of states localized within opposite ridges of potential wells. These ridges correspond to persistent currents of opposite chirality, and hence also to the absence/presence of a fluxon in the inductive loop of the circuit [31]. In this picture, m refers to the Fock order of the localized states. Finally, operators that flip Cooper pair parity correspond to odd functions of φ or functions of ϕ with period an odd division of 2π, which can be seen from Fig. 3c to mediate the transition |m+ ↔ |m− .

III. QUBIT
In Sec. II, we analyzed the multi-mode Hamiltonian describing the superconducting circuit in Fig. 2a. Numerical diagonalization of this Hamiltonian showed the emergence of a linear plasmon mode and a nonlinear fluxon mode [32]. In this section, we consider the properties of the logical qubit formed by {|0+ , |0− }, the two lowest-energy eigenstates at ϕ ext = π, which generalizes to {|0 • , |0 • } away from ϕ ext = π.

A. Matrix elements
To better elucidate which types of operators can and cannot induce transitions between the two states of the qubit, we examine the relevant matrix elements corresponding to capacitive and inductive coupling. This discussion is particularly relevant to understanding the expected dominant loss mechanisms (Sec. III C) and designing a measurement and control apparatus that does not directly couple to the qubit (Sec. IV C).
For capacitive coupling, a generic voltage V couples to the superconducting island of the circuit in Fig. 2a via a gate capacitance C g and will append the term to the Hamiltonian in Eq. 2, in addition to dressing the capacitance matrix. This voltage may be a degree of freedom of another mode in the embedding circuit, a noise source, or an ac drive. We therefore see that the susceptibility of undergoing a transition from the ground state, due to capacitive coupling to the qubit island, is directly related to the matrix element ψ|η|0 • . For inductive coupling, a generic current I couples to the circuit via an inductance L s shared with the inductive loop, which adds the term to the Hamiltonian in Eq. 2. Here, φ 0 = /2e is the reduced magnetic flux quantum and L is the superinductance in each arm of the qubit (i.e. L = φ 2 0 /L). Like the voltage source, this current may represent an internal or environmental degree of freedom. We see that the susceptibility of undergoing a transition from the ground state, due to inductive coupling to the inductive loop, is related to the matrix element ψ|φ|0 • .
Limiting the Hilbert space to the six lowest-energy eigenstates {|m± : m = 0, 1, 2}, we numerically compute the normalized matrix elements Results are plotted in Fig. 4. Note that ψ |O ψ | 2 = 1 and |O ψ | 2 > 0, so we may reasonably consider these as transition probabilities via O. We see from Fig. 4a that transitions mediated by capacitive coupling to the qubit island are only allowed from |0 • to |1 • . These selection rules result from the decoupling of the even and odd Cooper pair number parity manifolds in Eq. 5. Most importantly, transitions between qubit states are forbidden, meaning capacitive coupling offers a promising ingredient for indirect qubit measurement and control (see Sec. IV C). Conversely, inductive coupling to the inductive loop of the qubit permits transitions between |0 • and |0 • in the vicinity of ϕ ext = π, as shown in Fig. 4b. This effect arises because the operator φ induces transitions between the Cooper pair parity manifolds, as can be seen from the Fourier series for Eq. 3. As a consequence, we expect that relaxation of the qubit will be primarily due to inductive loss in the superinductances (see Sec. III C).

B. Disorder
A highly symmetric superconducting circuit is usually fragile in view of unavoidable fabrication imperfections [31]. The symmetry of our circuit involving the two inductive arms in Fig. 2a may be broken in three parameters: the Josephson energies of the junctions, the capacitances of the junctions, or the superinductances. To analyze these effects, we numerically diagonalize Eq. 2 and examine the energy splitting ∆E [33] as well as the charge dispersion = max Ng ∆E−min Ng ∆E of the {|0+ , |0− } manifold at ϕ ext = π. A dimensionless quantity δ ∈ [0, 1) is introduced to parameterize the extent of asymmetry in all three cases, and the δ-dependence of the energies ∆E and is studied.

Disorder in J
We model disorder in the Josephson energies of the junctions by allowing the values of J to deviate. We therefore set the left and right junction tunneling energies to (1 ± δ J ) J , respectively, where δ J is the aforementioned asymmetry parameter. The Hamiltonian in Eq. 2 is perturbed by the term See Fig. 5 for a plot of ∆E and as a function of δ J . The important feature in these plots is that the charge dispersion decreases exponentially while the splitting increases exponentially with δ J [34]. These features arise from the effective Hamiltonian in Eq. 4 being accompanied by 2πperiodic terms in the presence of disorder. In this case of disorder in J , the approximations in Sec. II B lead to H ≈ H eff with which evidently permits the tunneling of single Cooper pairs across the element. The resulting qubit retains characteristics of the symmetric circuit as well as the asymmetric/transmon-like circuit.

Disorder in C
Analogous to Sec. III B 1, we set the left and right junction charging energies to C /(1 ± δ C ), respectively. Aside from dressing the capacitance matrix, the Hamiltonian in Eq. 2 inherits the term , whose effect on the qubit manifold is plotted in Fig. 5.
This form for capacitive disorder is assumed due to the fact that, if δ J = δ C ≡ δ A , then the product J C is kept constant for both junctions under the effects of disorder. This corresponds to the physical case where the junction plasma frequencies are fixed by oxidation, but their areas differ due to fabrication imperfections. Here, the junction areas are (1 ± δ A )A because A ∝ J / C . The consequences of area disorder are plotted in Fig. 5.

Disorder in L
Following the same procedure, we set the left and right superinductance inductive energies to L /(1±δ L ), respectively. This form is taken in order to fix the total linear inductance in the loop. Aside from dressing the induc-tance matrix, the Hamiltonian in Eq. 2 is perturbed by Note in Fig. 5 that the charge dispersion and energy splitting follow the same general trend for inductive disorder as for the other three. The key difference is that the charge dispersion decreases more quickly than for any other form of disorder. Oppositely, the splitting is initially the same as for area disorder, but the slope decreases in δ L . We conclude that disorder allows us to engineer a circuit with a sufficiently non-degenerate ground state manifold whose charge dispersion is largely suppressed. For reasons that will become clear in Sec. III D, these features are extremely valuable for designing a qubit that is protected from dephasing.

C. Relaxation
We model loss due to an arbitrary channel using Fermi's Golden Rule. This gives us the relaxation rate, including both emission and absorption, of In this expression, O is the operator of the circuit coupling to a noisy bath variable E(t), the noise spectral density of which is given by S EE (ω). Note that ∆E = ∆ω.
In this subsection, we calculate the relaxation rates for the expected dominant loss mechanisms of the qubit based on numerical diagonalization of the full Hamiltonian in Eq. 2. We perform this calculation for various degrees of inductive disorder δ L (see Fig. 5). We consider four possible loss channels for the qubit: capacitive loss, inductive loss, Purcell loss, and quasiparticle tunneling [36]. Capacitive loss involves dielectric dissipation in the Josephson junction capacitances. In this case, we have O = −2e 8 C φ i for the charge across the i-th junction and E = V for a bath voltage with Here, C J = e 2 /2 C is the junction capacitance, T is the temperature, and Q cap (ω) is a frequency-dependent quality factor [40] with nominal value Q cap ∼ 1×10 6 [35] (see App. D for more details). Depending on the specific implementation, inductive loss may occur within the superinductances via quasiparticle tunneling [36]. This situation can be modeled by taking O = φ 0 φ i for the flux across the i-th superinductance and E = I for a bath current with In this expression, L i = (1 ± δ L )L is the i-th superinductance and Q ind (ω) is a frequency-dependent quality factor with nominal value Q ind ∼ 500 × 10 6 [36]. As shown in App. D, this frequency dependence is included to extrapolate to small qubit transition frequencies.
Loss due to the Purcell effect should mainly arise from coupling of the qubit to the plasmon mode, which we model as dielectric loss in the shunt capacitance. As such, we have O = −2eη and E = V for a bath voltage with the same noise spectral density as Eq. 8, but with C J replaced by C shunt .
Finally, quasiparticle tunneling is expected to occur across either Josephson junction and contribute to qubit relaxation [36]. In this case, we have O = 2φ 0 sin ϕi 2 for the i-th junction and the noise spectral density with Re Y qp (ω) being the dissipative part of the Josephson junction admittance [41] (see App. D for the explicit form). Note that this admittance scales linearly in both J and the quasiparticle density x qp defined relative to that of Cooper pairs. The calculated relaxation rates and the corresponding components of Eq. 7 are shown in Tab. II for all four loss channels at ϕ ext = π. A key feature is the complete absence of quasiparticle loss, as in the fluxonium qubit [36]. Additionally, we see that the asymmetric qubit is marginally less susceptible to inductive loss than the symmetric qubit. This improvement comes at the cost of the susceptibility to capacitive and Purcell loss. However, we emphasize that the lifetimes shown Tab. II are conservative estimates that demonstrate T 1 1 ms, which is at least competitive with state-of-the-art qubit implementations [13,28,36,[42][43][44].

D. Pure dephasing
In this subsection, we examine the dependence of the qubit transition energy ∆E on various system parameters λ corresponding to different dephasing mechanisms. We estimate the dephasing times due to offset charge noise, external flux noise, photon shot noise in the plasmon mode, and critical current noise. The noise spectral densities are all assumed to be 1/f and hence given by S λλ (ω) = 2πA λ /|ω|, where √ A λ is the amplitude, with the sole exception of shot noise whose spectral density is assumed to be Lorentzian.

Charge noise
In the case of perfect symmetry, the charge dispersion is identically mapped to ∆E. As a consequence, resilience to dephasing from offset charge noise demands a high degree of degeneracy, making experimental implementation difficult. Fortunately, this issue can be avoided by introducing inductive disorder into the system (see Sec. III B 3). In the slow-varying limit for charge   noise, the pure dephasing time T φ is bounded by the solution to where e is Euler's number and J 0 is the Bessel function of the first kind [10]. In Tab. II, we see that this yields a strict bound on the decoherence time in the case of perfect symmetry, which is greatly alleviated in the presence of inductive disorder. Note that this estimate does not explicitly involve A Ng because the slow-varying limit has been taken where the offset charge assumes a random value for each measurement, but does not fluctuate within a given measurement.

Flux noise
At ϕ ext = π, the qubit manifold is only susceptible to external flux noise to second order. However, one can easily obtain [45] the relation ∂ 2 ∆E/∂ϕ 2 ext ∝ 1/∆E (at ϕ ext = π), which shows that insensitivity to flux-noiseinduced dephasing requires sufficiently weak degeneracy, at odds with the requirement for resilience to charge noise. Once more, introducing inductive disorder avoids this issue. Neglecting components that depend on details of the implementation, the pure dephasing rate [46] is bounded by where A ϕext /2π ∼ 3×10 −6 is the amplitude of the noise spectral density [38,47] for ϕ ext . As in the case of charge noise, this places a stringent limit on the decoherence time of the qubit in the absence of inductive disorder, as we see in Tab. II.

Shot noise
We expect the dominant contribution to dephasing due to thermal photons to arise from coupling to the plasmon mode. This mode has an angular frequency ω p ≈ √ 16x C L / and an occupation n p with mean given by n th = 1/(e ωp/kBT −1). To this end, we place a bound on the pure dephasing time T φ using the expression [13] 1 where χ is the dispersive shift of the qubit on the plasmon mode and κ = ω p /Q cap (ω p ) is the linewidth of the plasmon mode. Note that, as in Sec. III C, we use the frequency-dependent dielectric quality factor Q cap (ω) with nominal value Q cap ∼ 1 × 10 6 [35]. In all cases, this form of noise is not expected to limit the decoherence time of the qubit (see Tab. II).

Critical current noise
The final dephasing mechanism we investigate is critical current noise, or fluctuations in the Josephson energy J . As in the case of flux noise, we discard factors that are sensitive to experimental details and take the dephasing rate to be bounded by [46] 1 where A J ∼ 5 × 10 −7 J is the amplitude of the noise spectral density [39] for J . As in the case of photon shot noise, we do not expect critical current noise to place a strict bound on the decoherence time (see Tab. II).
At this point, it is clear that inductive asymmetry constitutes a necessary ingredient in the protection of this qubit. In light of Eq. 6, we attribute this to the resulting hybridization of the φ and θ modes, which are not directly coupled in the symmetric case (see Eq. 2). The fluxon transition between the qubit states inherits some character of the plasmon transition, thereby breaking the correspondence between and ∆E and reducing the flux matrix element 0−|φ|0+ .

A. Protection
We have proposed a superconducting circuit whose ground state manifold can be made protected, at the Hamiltonian level, from multiple common sources of noise. Here we describe some of the other strategies for achieving protection. The simplest manifestations involve improving quality factors associated with coupling to different thermal baths [42,48,49] and reducing noise spectral densities for different Hamiltonian parameters [13,50,51]. At a higher level, a popular approach to suppressing qubit relaxation has been to localize wavefunctions in disparate regions of phase space to lessen transition matrix elements [28,30,36,52]. On the other hand, delocalization of the same wavefunctions has been shown to mitigate dephasing effects by reducing qubit sensitivity to Hamiltonian parameters [10,18,43,[53][54][55]. Superconducting circuits with multiple degrees of freedom whose qubit wavefunctions are both localized and delocalized combine both of these approaches. In these circuits, quantum information is diffused among constituent local degrees of freedom, providing protection from local perturbations. The many-body limit promises topological protection [11,16,20,56], in which global operators are necessary to manipulate logical qubits, but the few-body case described here offers an experimentally realistic approximation [12,21,28,31,44,47,57]. In the following subsection, we outline the key differences between our circuit and similar proposals.
B. Comparison to other protected qubits

Rhombus qubit
The circuit in Fig. 2a bears resemblance to the rhombus [58], with two central differences. First, in the rhombus, the superinductances are replaced by single Josephson junctions, or arrays of a few larger junctions [9,28]. This changes the parameter regime of the circuit from L J to L ∼ J , decreasing the amplitude of the cos 2ϕ term in the Hamiltonian. Second, the shunt capacitance is replaced by a gate capacitance to a voltage source [28]. This is akin to substituting an electrostatic gate for a shunt capacitance to obtain the Cooper pair box [54] from the transmon [10]; overall suppression of the charge dispersion is traded for the ability to bias the circuit at its charge sweet spot.
Finally, the rhombus by itself is not designed to be a protected qubit. Rather, when multiple rhombi are arranged into a one-dimensional chain (or a twodimensional fabric), the ground states are eigenstates of a nonlocal operator, which provides topological protection [9,11,16,56]. On one hand, our qubit does not require such scaling to achieve protection. On the other hand, the protection we predict is inherently susceptible to local perturbations (see Secs. III C, III D).

0-π qubit
The 0-π qubit [21] also has a similar superconducting circuit to that in Fig. 2a, but with three essential distinctions. First, the pairs of superinductances on each arm of the circuit are combined, altering the embedding capacitance matrix, in particular the capacitances to ground. Second, there is the addition of a second large capacitance shunting the inductive loop between its two horizontally oriented nodes. When this second capacitance is precisely C shunt , this permits the exact decoupling of the θ mode from the ϕ mode in Eq. 2. Third, the 0-π qubit is operated in a parameter regime where L ∼ 0.01 J , as opposed to our circuit, where L ∼ 0.1 J [31]. This additional order of magnitude in the superinductance may prove marginally less accessible experimentally [18,19].
Notably, the inductive loop in the 0-π qubit is threaded with ϕ ext ≈ 0 at its working point instead of ϕ ext = π. This leads to a substantial change in the physics of the ground state manifold; |0+ and |0− are approximately localized in distinct potential wells [31]. In our case, these states are approximately the symmetric and antisymmetric superpositions of the localized wavefunctions, which are themselves localized in distinct Cooper pair number parities.

C. Readout and control
Protected qubits face the serious obstacle of realizing state manipulation and measurement while remaining sufficiently isolated from their environments to preserve their coherence. We envision performing readout and control using an ancillary mode structure, which enables cascaded dispersive readout and Raman indirect transitions, as outlined below.
For readout, we aim to exploit the sizable native dispersive shift of the qubit on the plasmon mode, χ/2π ∼ −20 MHz (see Sec. III D 3). Unfortunately, the small anharmonicity, which is at most of order 10 MHz, of the plasmon mode makes readout of the plasmon mode using a linear mode (e.g. a microwave cavity) difficult. As a remedy, we propose introducing an ancillary anharmonic mode by which to measure the plasmon state. This scheme involves incident microwave tones at the transition frequencies of the ancilla mode, the plasmon mode, and the qubit; and hence constitutes cascaded dispersive readout.
The above-mentioned ancillary anharmonic mode will also be useful for control of the protected qubit, where direct transitions, mediated by capacitive coupling to the qubit superconducting island, are forbidden (see Sec. III A). Note that this coupling method is chosen to minimize the contributions to decoherence resulting from the introduction of the readout/control circuit. For manipulation, we therefore propose transitioning indirectly through the |1+ or |1− state. In this case, the central challenge is to engineer selection rules to allow the |0± ↔ |1± transitions; however, this has recently been demonstrated in a fluxonium qubit coupled to its readout resonator using Superconducting Nonlinear Asymmetric Inductive eLements (SNAILs) [59,60]. We propose using Raman indirect transitions to avoid decoherence of the intermediate |1± state.

V. CONCLUSION
In summary, we have designed a few-body superconducting circuit in which the charge carriers are wellapproximated by pairs of Cooper pairs at a particular bias point. The Josephson tunneling element that supports these charge carriers is characterized by a cos 2ϕ term in the Hamiltonian, whose emergence we have shown analytically. Our numerical simulations supplement these arguments and demonstrate protection against a variety of common relaxation and dephasing sources. We find that this protection is substantially enhanced in the presence of disorder. Finally, we compared our circuit to similar proposals and offered our perspectives on readout and control.
As a final remark, we comment that engineering a circuit whose potential energy is dominated by a cos 2ϕ term opens the door to more exotic designs with potential energies of the form cos µϕ, with µ ∈ N, which could be obtained by introducing additional loops in the circuit. These could be tremendously valuable for quantum simulation, realizing nearly degenerate ground state manifolds with greater multiplicities, or performing degeneracy-preserving measurements of photon number parity [61].
Using the trial function g ≡ e −iNgϕ ψ, this equation becomes which is the Mathieu equation [10,62]. The k-th eigenenergy is given by where a denotes the Mathieu characteristic value and is an integer that sorts the eigenvalues. The crucial difference between this solution and that for the Cooper pair box is that the characteristic exponent 2(N g + ) is replaced by N g + [10]. This gives rise to the nearly degenerate ground states in the Hamiltonian in Eq. 1.
In particular, for E J E C , we can use the asymptotic form of the Mathieu characteristic values to find with k being the charge dispersion of the k-th level. Importantly, for even values of k, we have E k = E k+1 ≈ ( √ 32E J E C − 2E C )k (neglecting constants and nonlinear terms) and k = − k+1 . Additionally, the expression for the ground state splitting in Sec. I B is justified by and the fact that ∆E = E 1 − E 0 = 0 cos(πN g ).

Appendix B: Instanton treatment
The average tunneling trajectory between potential minima for the Hamiltonian in Eq. 2, upon constraining the dynamics to the ϕ 1 ϕ 2 -plane, can be found by examining an equivalent classical problem. This is done by inverting the potential and solving the classical equations of motion for a particle that starts at one maxima and ends at the other, each with zero kinetic energy [22]. Consequently, this classical trajectory requires an infinite amount of time. The corresponding differential equations, are solved numerically. Carrying out the procedure in Sec. II B, but retaining terms through the fourth harmonic and O(z 2 ), results in the extended version of the effective Hamiltonian which clearly reduces to Eq. 4. At ϕ ext = π, we see that the dominant correction term is of the form cos 4ϕ. Moreover, we observe that the odd Fourier terms exactly vanish for the symmetric circuit at this flux bias, while the coefficients for the even Fourier terms decay roughly in powers of z.

Appendix C: Numerical diagonalization
For numerical diagonalization of either Eq. 2 or Eq. 4, we employ a charge basis for the ϕ mode in order to efficiently capture the dependence of the Hamiltonian on N g [63]. This is due to the exact periodicity of the Hamiltonian in ϕ, which inhibits the use of a harmonic oscillator basis to capture the multi-well physics apparent in Fig. 2b. For the θ and φ modes, we use harmonic oscillator bases. Explicitly, the full Hamiltonian in Eq. 2 can be written as after discarding constants, where a † /a and b † /b are bosonic creation/annihilation operators for the φ and θ modes, respectively. Additionally, the zero point fluctuation amplitudes φ zpf = 8 C L 1/4 and η zpf = 1 have been introduced for convenience. The diagonalization basis is then with |N being the Cooper pair number eigenstates of the operator N and |p /|q being the photon number eigenstates of a † a/b † b. For numerical convergence, the dimensions N 0 7, p 0 7, and q 0 30 are used, depending on the desired truncation accuracy. The remaining operators are written as in this basis, while the final cosine term in Eq. C1 has an exact matrix expression [23], which we do not repeat here. Although the above discussion applies to the case of perfect symmetry, the discussion in Sec. III B allows a straightforward extension to the disordered case. Finally, we comment that, despite the high dimensionality of the matrix in Eq. C1, diagonalization can be made efficient by taking advantage of its sparsity. Now we detail the method used to compute the wavefunctions N |ψ and ϕ, φ|ψ (up to an overall phase) for |ψ = |m± plotted in Figs. 3b, 3c. First, the coordinate space wavefunction ϕ, φ, θ|ψ is constructed using ϕ|N = e −iN ϕ [63] and the expressions for φ|p and θ|q as harmonic oscillator wavefunctions: ϕ, φ, θ|ψ = where the complex coefficients c N pq are obtained by diagonalization (and consistent choice of overall phase). Then, the θ-dependence is removed by projection onto θ = 0 and normalization. This amounts to computing ϕ, φ|ψ = ϕ, φ, 0|ψ which is plotted in Fig. 3c. Next, the φ-dependence is removed by invoking Eq. 3. In this step, we project φ onto φ(ϕ) and normalize to obtain ϕ|ψ = ϕ, φ(ϕ)|ψ  δ L , we plot in Fig. 6 the results for the entire span of δ L . In addition, in Fig. 6c we also plot the expected relaxation, pure dephasing, and decoherence times combining the effects of all considered channels. Here, the overall decoherence time is defined using the relation 1/T 2 = 1/2T 1 + 1/T φ . We also note that, for the estimate of the dephasing rate due to charge noise, we do not use the extrapolation shown in Fig. 5a for the cases of disorder in J , C , and A . Numerical instabilities are not encounted for disorder in L until δ L ≈ 0.95 (which is outside the plotted range for both Fig. 5a and Fig. 6b).
In the remainder of this appendix, we describe the assumed frequency dependence of the three quality factors used for the relaxation time estimates in Sec. III C. First, for the dielectric quality factor, we use the form [36,40] Q cap (ω) = Q cap 2π × 6 GHz |ω| 0.7 so that the nominal value Q cap ∼ 1 × 10 6 corresponds to measurements performed at the resonant frequency of 6 GHz [35].
Second, for quasiparticle loss, the dissipative part of the admittance for a Josephson junction with tunneling energy J is given [41] to be Re Y qp (ω) = 2 π where R K = h/e 2 is the resistance quantum, ∆ is the superconducting gap, and K 0 is the modified Bessel function of the second kind. This expression holds under the assumption that the quasiparticle bath is in thermal equilibrium at temperatures T ∆/k B , which may not be correct, depending on the implementation [64]. As we have seen, the matrix elements for quasiparticle loss vanish for this circuit, so this admittance has little bearing on our relaxation time estimates.
On the other hand, inductive loss has been suggested to occur due to quasiparticle tunneling across Josephson junctions that compose superinductances [36]. If we expect the frequency-dependence of the dissipative part of the admittance to agree between quasiparticle and inductive loss, then we arrive at so that the nominal value Q ind ∼ 500 × 10 6 corresponds to measurements performed at the resonant frequency of 0.5 GHz.