Efficient one- and two-qubit pulsed gates for an oscillator stabilized Josephson qubit

We present theoretical schemes for performing high-fidelity one- and two-qubit pulsed gates for a superconducting flux qubit. The"IBM qubit"consists of three Josephson junctions, three loops, and a superconducting transmission line. Assuming a fixed inductive qubit-qubit coupling, we show that the effective qubit-qubit interaction is tunable by changing the applied fluxes, and can be made negligible, allowing one to perform high fidelity single qubit gates. Our schemes are tailored to alleviate errors due to 1/f noise; we find gates with only 1% loss of fidelity due to this source, for pulse times in the range of 20-30ns for one-qubit gates (Z rotations, Hadamard), and 60ns for a two-qubit gate (controlled-Z). Our relaxation and dephasing time estimates indicate a comparable loss of fidelity from this source. The control of leakage plays an important role in the design of our shaped pulses, preventing shorter pulse times. However, we have found that imprecision in the control of the quantum phase plays the major role in the limitation of the fidelity of our gates.


Introduction
Superconducting circuits containing Josephson junctions [1,2] are widely recognized to be promising systems for the physical implementation of quantum bits. These systems can be made and operated using well-established experimental techniques, and they have the clear potential, in principle, for scalability. Important experimental milestones in coupled superconducting qubits, including the observation of two-qubit gates [3,4] and the measurement of entanglement [5], have already been reached. But because superconducting qubits are condensed matter systems, it has been a hard task to isolate them from their environment. Because of that, these systems tend to suffer from short coherence times, which has imposed serious limitations on achieving very high fidelity gates.
In this paper, we report a theoretical study of a universal set of one-and two-qubit gates implemented using only shaped dc flux pulses for an oscillator stabilized flux qubit. We introduce a simplified but accurate model to describe the dynamics of the lowest states of the qubit as a function of the external control parameters. This model provides a simpler way to analyze the physics of the problem, helping to make it easy to see what operations are necessary for performing the desired quantum gates. As we shall see, there is a trade-off between the speed of a gate and the amount of leakage produced by it (leakage = evolution of states out of the 0-1 computational basis). Smart choices for the shape of the pulses are required to keep leakage at a tolerable level, such that the gate fidelities are not compromised by this process [6,7]. We have been able to keep leakage at the 0.12% level for gate times of the order of a few tens of nanoseconds.
The other important goal of the present work is to consider the effect of low-frequency noise during the gate operation and in the memory state. When designing quantum gates, we search for the best pulse paths such that the loss of fidelity due to the low-frequency noise is minimized. We find that for all gates of interest, the loss of fidelity due to low-frequency noise is never greater than 1%. Achieving this low level of infidelity requires a careful use of symmetries (so that errors accumulated in the first half of a pulse can be cancelled out in the second, for 3 example) and of 'sweet spots' [8,9] (points in the control parameter space that are first-order insensitive to fluctuations).
The last part of this paper analyzes a pulsed two-qubit gate, the controlled phase gate. An interesting feature of the coupling used in this gate scheme is that, although the physical inductive coupling between the qubits is assumed fixed, the effective qubit-qubit interaction is tunable as a function of the control parameters on both qubits, due to the change of character of the bare qubit states as the control parameters are varied. In fact, the effective qubit-qubit interaction can be made negligible in large regions of the flux space, allowing us to perform high-fidelity one-qubit gates for the two-qubit system without extremely stringent control over electrical cross-talk.
The outline of this paper is as follows. In sections 2 and 3, we discuss the physics of our system, and present the simplified model used to simulate the dynamics of the lowestlying levels of an oscillator-stabilized flux qubit. The appendix carries out a detailed derivation of the system Hamiltonian, and the regime of validity of the simplified model is discussed. Section 4 presents the different schemes designed to perform the following basic single-qubit operations: measurement in the standard basis (0/1), measurement in the conjugate basis (+/−), Z -rotation gates and the Hadamard gate. The fidelity of all these gates is given as a function of unwanted shifts from the optimal point of operation in flux and time synchronization of the pulses. In addition, we characterize the noise through the operator-sum representation of the system superoperator. In section 5, the two-qubit system is analyzed. The form of and the reasons for an effective tunable interaction are discussed. Then, a gate in the equivalence class of the controlled-Z gate is proposed, and its fidelity and a characterization of the nature of the noise are presented. Finally, section 6 gives some conclusions.

System Hamiltonian
Our analysis will be focused on the IBM qubit [10,11]. This device, shown in figure 1, consists of a bare qubit, which is a type of flux qubit containing three Josephson junctions and three loops, and a high-quality superconducting transmission line [9,11,12]. The bare qubit are subject to external control via flux lines which change the total magnetic fluxes threading the loops. As previously reported [13], and summarized in the appendix, the bare qubit has a gradiometric structure, and hence its behavior, to good approximation (see the appendix for the validity of this assumption), is only a function of the difference of the magnetic flux in the two large loops, which we denote as ≡ − p . Whenever is an integer multiple of the flux quantum 0 = h/2e, the system potential has a perfectly symmetric structure. These lines in the flux space are referred to as the 'S lines' (S for symmetric).
The bare qubit has two control parameters that define its flux space: one, , causes departures from the parity symmetry just described, and the other, the so-called control flux˜ c , changes the structure of the qubit potential, varying the height of the potential barrier between the two classically stable states. As˜ c is changed, the system potential passes from a doublewell to a single-well structure. This change in the form of the potential as a function of a control parameter produces an essential feature in the quantum behavior of the system: a well-defined change in the character of the qubit eigenstates. While for the double-well regime the qubit eigenstates are almost localized orbitals in the left and right wells, as the single-well regime is approached the qubit eigenbasis goes over to delocalized states, being close to states that are symmetric and antisymmetric with respect to reflection around the midpoint of the potential. In order to operate the qubit, flux lines (highlighted in black tracks) are used to change the total fluxes threading the loops. Readout SQUIDs (two such structures are shown on the left-hand side of the picture) perform the measurement of the state of the qubit. (c) Proposed scheme of the two-qubit system. The qubit-qubit interaction is assumed to arise via the indicated mutual inductance between the two big loops [14]. This results in a qubit-qubit interaction Hamiltonian of the formσ z ⊗σ z . Physical parameters for this qubit: the capacitances and critical currents of the junctions are assumed to be C =10 fF and I c = 1.3 µA; L T = 5.6 nH, L 1 = 32 pH and L 3 = 680 pH are the transmission line, the small and big loop inductances, respectively. The mutual inductances between qubit-transmission line, small loop-control flux line and big loop-bias flux line are respectively 200, 0.8 and 0.5 pH. The transmission line is designed to have a fundamental mode frequency of ω T = 2π × 3.1 GHz. Finally, the qubit-qubit mutual inductance is M = 12 pH.
This change of the character of the states happens over a very small interval of˜ c , because of the exponential increase in the amplitude in tunneling between the wells as the barrier between the two wells is decreased. This transition interval is referred to as the 'portal' [10,11], and pulsing the flux parameters through this portal can create some of the elementary actions needed for the construction of quantum gates. For example, a non-adiabatic pulse through this region can create superpositions between states |0 and |1 .
The fundamental mode of the open-ended transmission line acts as a harmonic oscillator of frequency ω T coupled to the bare qubit. The presence of that structure modifies the quantum behavior of the qubit when the energy splitting of the ground and first excited states of the bare qubit is comparable to or larger thanhω T . In this regime, the two lowest-lying states of 5 the system both have the bare qubit in its ground state; they differ only in their transmission line quantum number. By tuning the energy splitting of the ground and first excited states, one can move information stored in the bare qubit to the transmission line and vice versa. When the qubit has been transferred to this transmission-line embodiment, we say that it is parked. Parking results in a very useful stabilization of the 0-1 frequency as a function of changes in {˜ c , }. In addition, the quantum coherence times while parking are seen to reach several microseconds [15]. Thus, the parking regime will be used as the memory state of the qubit; it will stay in this state when it is awaiting operation. Far away from parking, when the 0-1 energy gap is much smaller thanhω T , the transmission line does not play any role in the dynamics of the 0-1 states, which are just those of the bare qubit. Since the coherence times here are expected to be of the order of only tens of nanoseconds, this regime should be avoided; this part of the parameter space will only be used for state measurement.
So far, we have given a qualitative description of the qubit dynamics. This description is made quantitative with the methodology introduced by Burkard, Koch and DiVincenzo (BKD) [16]. Using the network graph theory, BKD developed a universal method for analyzing any lumped element electrical circuit containing Josephson junctions. The result of this theory is a mapping of the circuit dynamics to that of a massive particle in a potential, whose mass tensor and degrees of freedom are associated with the system capacitances. The system Hamiltonian in this formulation is given by Here L −1 J ;i ≡ (2π/ 0 )I c;i , where I c;i is the critical current of junction i. The diagonal matrix C contains the capacitances of the system. The topology of the circuit is encoded in the matrices M 0 ,N andS (see the appendix for a more detailed presentation of the BKD theory). The first term of the potential equation (2) represents the energy due to the presence of the Josephson junctions. The second term is associated with the inductive energies of each branch of the circuit. The last two terms take into account coupling to external sources of magnetic fluxΦ x and current sources I B . Quantization of the system is introduced by imposing the canonical commutation relation for the variables of charge, Q C , and phase, ϕ: The analysis of the system potential equation (2) for our qubit reveals that, instead of working with the real applied fluxesΦ x ≡ {˜ c , , p }, it is more convenient to introduce 'nonorthogonal' flux coordinates ≡ − p and c ≡˜ c + (where L 1 stands for the small loop inductance, L 3 for the big loop inductance, and M 15 and M 35 for the the small-big loop and two big loops mutual inductances, respectively). In fact, as presented in the appendix, the potential has a definite symmetry as a function of { c , }, rather than {˜ c , }. Since the potential symmetry is a key feature of the system, and we will explore it 6 when designing our gates, an accurate calculation must consider that fact, even though the term proportional to represents a small correction to the real applied flux˜ c . From here onwards, we will only refer to the pair of effective fluxes { c , }.
Since our qubit has four capacitances (three associated with the Josephson junctions and one representing the fundamental mode of the transmission line), the BKD theory leads us to a four-dimensional (4D) potential, for which direct calculations and analysis are more difficult. Thus, we follow the procedure of [13] to reduce the system dimensionality to two, one representing the bare qubit and the other representing the transmission line. The procedure involves organizing the degrees of freedom into 'fast' coordinates, in which the potential rises very steeply (such that the system dynamics are frozen into ground state along these direction), and ones that are 'slow'. Then, using a Born-Oppenheimer approach, the fast coordinates are traced out, resulting in small modifications to the remaining slow-coordinate potential energy. Figure 2 shows the first seven levels of the system spectrum, calculated following the above described steps, as a function of c , for three different values of the bias flux . Note that, by convention, the ground level is always at E = 0. It is clear that the equally spaced harmonic oscillator levels cut through this spectrum for all values of c and ; since the transmission line sees the control flux only via interaction with the bare qubit, its energies are very stable except in the vicinity of energy crossings. Observe that the lowest-lying states in the 'parking' regime, at high values of c , involve only transmission line quantum numbers (i.e. they are the states of a harmonic oscillator). This is true because the energy splitting of the eigenstates of the bare qubit becomes much larger than thehω T energy splitting of the transmission line states.
In addition, because of the interaction between the bare qubit and the transmission line, another structure present for each is an observed avoided crossing gap between the bare qubit and the transmission line states, occurring close to c = 1.45 0 . As one can see, the bias flux plays an important role for small values of c , where the system potential has a double-well structure, and the bare qubit energy splitting is smaller thanhω T . The first plot presents the case on the S line, = 0. Since at the S line the system potential is symmetric, the qubit states are symmetric and antisymmetric superpositions of the degenerate localized orbitals of the left and right wells of the potential. As we move away from the S line ( = 0), the symmetry of the potential is broken and the localized left and right orbitals are no longer degenerate. This explains the appearance of the gap between the ground and first excited states in the second and third plots at small values of c . One important feature of the system potential is a symmetry in going from + to − (see the appendix); under this transformation, the energies of the left and right states are interchanged. Finally, in the third plot it can be seen that the first excited state stabilizes at the transmission line frequency for c < 1.43 0 . This occurs because for most values of c the bias term in the Hamiltonian is large enough by itself to make the bare qubit energy splitting larger thanhω T (but not much larger). As a result, for small c the first excited state corresponds to one excitation of the transmission line, and the second excited state to the excited state of the bare qubit. The qubit is encoded using the lowest two eigenstates. For large values of c , the energy splitting of the bare qubit states is much larger thanhω T , so that the lowest-lying qubit states have purely harmonic oscillator character. (a) The level structure on the S line. In this case, the system potential has a symmetric double-well structure for small values of control flux, consequently the ground and first excited states are nearly degenerate there. (b) and (c) present cases away from the S line. Here, the states |L and |R are no longer degenerate and the system has a gap between the ground and first excited states for small values of control flux.

Four-level model
Although it would be possible, with considerable computational effort, to design and simulate the desired quantum gates by direct evaluation of the time-dependent Schrödinger equation for the circuit Hamiltonian equation (1), a considerable economy is achieved by introducing a simplified model that correctly describes the lowest eigenstates, since we will be only interested in their dynamics for our quantum gates. We start the derivation of such a simplified model by using the first-order perturbation theory to treat the bare qubit coupling with the flux line and with its transmission line. The resulting Hamiltonian due to that approximation has the form whereσ i are the Pauli matrices,σ z |L ≡ −|L ,σ z |R ≡ +|R , and {â † ,â} are canonical bosonic creation and annihilation operators. The states |L and |R represent the localized orbitals found in the left and right potential wells. Because of the change of the state character as a function of c , the amplitude of tunneling , the bias term coefficient b and the qubit-transmission line coupling g also become control flux-dependent. Figure 3 presents their behavior as a function of c . As one can see, the tunneling amplitude becomes negligible for small values of c . This happens because, in this regime, the very large barrier between the two wells ( 100 GHz) in the double-well potential exponentially suppresses the tunneling between their lowest states. However, as the value of c increases, the barrier decreases and the two minima become closer. Consequently, the amplitude of tunneling rapidly increases, and the left and right orbitals start to become more and more delocalized. Around c ≈ 1.45 0 the barrier vanishes rather abruptly and the bare qubit potential enters a single-well regime. The region around c ≈ 1. Once the bare qubit reaches the single-well regime, those coefficients become less sensitive to changes in c , since the nature of the states barely changes as one passes through this regime.
In addition, the inset shows that the amplitude of tunneling can rapidly reach values of several tens of GHz in the parking regime, which are, in general, much larger than the typical value of the transmission line energy splittinghω T ; hence, the lowest-lying states in the parking regime involve only transmission line states.
Unfortunately, even though the Hamiltonian equation (4) already represents an important simplification for the system description, since it has reduced the bare qubit Hilbert space to that of a two-state system, we cannot obtain an analytical solution for its eigenstates/eigenvalues. Thus, we still have to deal with an infinite set of states due to the harmonic oscillator. As we shall see, for the purposes of our work, we must have not only the correct dynamics of the states |0 and |1 , but also an excellent agreement for the minimal gap between the computational basis {|0 , |1 } and the rest of the spectrum of the system. Because of that, the truncation of the harmonic oscillator Hilbert space at its first excited state does not work as a fair approximation of our system Hamiltonian. As presented in figure 4, the simple truncation of the harmonic oscillator Hilbert space fails to give a good description of the gap between the |1 and |2 states. This happens because at the regime of small values of c the function g is appreciable. Hence, the shift applied to the oscillator due to the interaction with the bare qubit becomes important. Consequently, the most adequate description of the system is obtained using the representation of shifted harmonic oscillator states, where the truncation of the harmonic oscillator Hilbert space should give better results.
In fact, a more careful inspection of the Hamiltonian equation (4) reveals that its harmonic oscillator part has the form of a well-known shifted oscillator. Thereby, changing to the representation of the shifted harmonic oscillator states should give us a better picture of the system dynamics, and the appropriate basis in order to perform further truncations of the Hilbert space. Nevertheless, the shift applied to the harmonic oscillator is dependent on the bare qubit state, as can be seen by the term g( c )(â +â † )σ z of equation (4). This important feature of the system leads us to introduce the unitary transformationD(s,σ z ) ≡ e (sâ † −s * â )σ z as the conditional displacement operator of the system. Observe that the operatorD does not commute with the spin part of the Hamiltonian equation (4). Thus, changing the Hamiltonian representation to the shifted harmonic oscillator states, H →D † HD, and then performing the truncation of the harmonic oscillator Hilbert space at its first excited state, we arrive at the following simplified four-level model: where we have defined the operatorsĉ ≡ |0 HO 1 HO | andĉ † ≡ |1 HO 0 HO |, in which |0 HO and |1 HO represent the ground and first excited harmonic oscillator states, respectively. During the procedure described above, we have introduced an ad hoc parameter, the shift parameter s, which can be used to parameterize the Hamiltonian equation (5) in order to obtain the closest level dynamics as possible to that determined by the Hamiltonian equation (1).
A natural choice for the parameter s would be the standard value s = −g( c )/hω T , which leads to the cancellation of the last two terms of equation (5), and the Hamiltonian diagonalization when ≈ 0. However, since the shift imposed on the harmonic oscillator is conditioned on that of the bare qubit state, and the form chosen for s does not take into account the change of character of the bare qubit state, one should expect it would fail when the regime of high tunneling amplitude is reached. Indeed, as can be seen in figure 4, this choice for the parameter s only gives the correct description when the tunneling amplitude is  (5), using the shifts s = −g/hω T (trianglesymbolized curve) and s = −g/ ω 2 T + 2 (dashed red curve), equation (6), for three different values of the bias flux . The simple truncation fails when g becomes appreciable (small values of c ). In this regime, the shifted harmonic oscillator states are the preferred system representation. The four-level model using the shift s = −g/hω T does not provide the parking harmonic stabilization, since it does not consider the change of character of the bare qubit. At last, the four-level model using the conditional shift equation (6) gives a very fair approximation of the lowest-lying levels for all regimes of the bias and control flux.
negligible. In fact, the analytical solution of equation (5) using the above shift reveals that, at the limit → ∞, one should expect ω 01 → ∞, which is in complete opposition to the parking stabilization expected for the system eigenstates.
Therefore, a smart choice of the parameter s has to consider the fact that the harmonic oscillator shift is conditioned to the bare qubit state, and that its state character changes as a function of c . As we already know, at the limit → 0, the parameter s should asymptotically reach the value s → −g( c )/hω T , since it correctly decouples the spin and harmonic oscillator degrees of freedom in the Hamiltonian equation (5). In the other limit, → ∞, one should expect observing no shifts to be imposed on the harmonic oscillator, since the system states would be frozen at the bare qubit ground state. Thus, we expect a reasonable interpolation for the parameter s between the two regimes to be In fact, as presented by the dashed curve in figure 4, the form equation (6) for the parameter s gives a very fair approximation for the lowest level dynamics, in particular the energy splittings ω 01 and ω 12 . In addition, from the analytical solution of equation (5) using equation (6), we can check that the ω 01 has the correct harmonic oscillator stabilization, i.e. → ∞ ⇒ ω 01 → ω T .
It is worth pointing out that, although the model equation (5) does not give the correct parking stabilization for the system second excited state |2 , since this state involves the second excited state of the transmission line, for the purposes of our work, it turns out that is not a limitation for the model. Indeed, the lack of parking stabilization for the second excited state |2 would lead to underestimations of the transitions between the states |1 and |2 , and consequently wrong leakage estimations. However, as we are going to assume ω T of the order of several GHz, the energy splitting ω 12 for the parking regime is big enough to suppress those transitions during our dc shaped pulse operations. Furthermore, since at the parking regime the ground state, and the first and second excited states have the same nature (i.e. they are the states of the harmonic oscillator), observing the induced 0-1 transitions due to the dc pulse operations also gives a measurement of the expected 1-2 transitions in this regime. Thus, we end with a very controllable 4D Hilbert space model, spanned by the ground and first excited states of the bare qubit and transmission line, that accurately mimics the main features of the system dynamics determined by the exact circuit Hamiltonian model equation (1). From now onwards, we will only use the Hamiltonian equation (5) with (6), when designing and simulating the quantum gates.
Completing the system description, figure 5 presents the plots for the frequency difference between the ground and first excited states, ω 01 , and the frequency difference between the first two excited states, ω 12 , as a function of { c , }. As mentioned before, the ω 12 plot quantifies the minimal gap between the computational basis {|0 , |1 } and the rest of the spectrum of the system, giving a measure of how likely the system is to undergo leakage during the sweep of a pulsed gate. The ω 01 plot shows how fast the relative phase of the two computational basis states advances if the system is held at some flux values. We envision that there will be a very stable master clock at the transmission-line frequency, and that errors in the accumulated phase difference are unlikely if the system is held in parking, where its phase advance is synchronous with this clock. But when ω 01 departs far from ω T , phase accumulation with respect to the reference is very fast, and we expect (and our numerical studies confirm) that the system is much more susceptible to phase errors in this regime. These plots provide a good 'map' for designing the one-qubit gates, since they indicate the regimes where one should expect an appreciable amount of leakage and very fast relative phase accumulation. The ω 12 plot quantifies the minimal gap between the computational basis {|0 , |1 } and the rest of the spectrum of the system, giving a limit on the rate at which the system can evolve without appreciable leakage. (b) The ω 01 plot indicates the phase accumulation rate with respect to the reference. In regions of very fast rates, a very precise control of the external applied fluxes is required in order to avoid phase noise. Illustrated with dots is the path in the flux space used to implement the Hadamard gate. The number of dots indicates the time spent when passing through that region.

One-qubit gates
Before we move on to the discussion of each gate individually, it is worth summarizing some features present in all of them (including the two-qubit gate). First, the memory and measurement points must be defined. The parking point is taken to be at c = 1.6 0 , and the measurement point to be at c = 1.4 0 . These choices are optimal, in the following sense: for the parking point, the important feature to be considered is how much the ω 01 frequency changes as a function of external controls. Since we envision the qubit spending long times at this position, it is imperative that the frequency deviation δω 01 has very small values for reasonable flux shifts due to the noise. We calculate that the frequency sensitivity is δω 01 /δ c ≈ 10 MHz/ 0 at the chosen memory point, which will give acceptably small memory error. In addition, inaccuracies in the transmission line fundamental model frequency can be corrected using a π -pulse stabilization scheme [17]. This operation (not analyzed in detail here) consists of initially applying a π-rotation gate; the system is then left to evolve in a free evolution for the same amount of time previously spent at the memory state, and finally a new π -rotation is applied. As a result, undesired Z -phase accumulated at the parking point ends as a global phase of the system state.
For the measurement point, we must make sure the left and right states are experimentally distinguishable, and that, since the potential has a periodic structure [13], the potential barrier between different pairs of minima is high enough to avoid tunneling between their states. At c = 1.4 0 , we found that the barrier between the two principal minima is of the order of 10 THz, and the barrier to other minima is even higher (about 20 THz).
Another common feature of all our shaped pulses is our design methodology. We shape our dc pulses using simple sums of tanh functions. As we shall see, our gates can be divided into several parts, and with each part is associated a function of the form δ tanh((t − t max )/τ ), where δ, τ and t max determine the flux excursion, the maximum rate and its time position, respectively. Thus, with these parameters, one can adjust the rate and the flux excursion of each part of the gate, in order to optimize the gate. It turns out that our designed gates do not require maximum flux slew rate and a bandwidth higher than 7 × 10 6 0 s −1 and 1 GHz, respectively. We measure the gates' fidelity using the entanglement fidelity [18] where A Q ≡ U † ideal U real (δ c , δ , δt), with U ideal and U real (δ c , δ , δt) representing the ideal gate and the final achieved transformation, respectively. ρ Q is an equal distribution of the computational basis states |0 and |1 : ρ Q ≡ 1 2 |0 0| + 1 2 |1 1|. Finally, we estimate the probability of leakage using the projection of an arbitrary evolved state, where P j represents the projection operator in the state | j . Figures 6-8 show the proposed measurement gates, phase gate and the Hadamard gate, whose matrix representations in the {|0 , |1 } basis are respectively given by where θ 01 and θ +− are arbitrary phases, whose values are not relevant for the gate implementation, provided that the gates equations (8) and (9) The gate consists of two different regimes: the first, t 7 ns, is designed to be an adiabatic process, thus minimizing leakage when passing through the avoided crossing. The other, t > 7 ns, evolves the system in a non-adiabatic process through the portal, leading to the desired transitions between the |0 and |1 states. The applied flux is maintained constant during the gate, and close to the S line: = 30 µ 0 . Insets: the gate fidelity as a function of unwanted shifts in c and from the optimal point of operation. The main source of the loss of fidelity for the 0/1 measurement is related to 0/1 transitions, and for the +/− measurement it is related to the phase noise. Note the difference in the scales. equation (10). Finally, since the states |0 and e iθ 0 |0 , and |1 and e iθ 1 |1 are physically identical, the convention for the global phases θ 0 and θ 1 is chosen such that the physical implementation of the Hadamard gate has the matrix representation given by equation (11). Once fixed, the phase convention is maintained for all other gate implementations.
In addition, figures 6-8 also present the gate fidelities as a function of unwanted constant shifts from the optimal point of operation in the whole pulse profile. We model the effect of 1/ f flux noise by an ensemble of random constant shifts of this sort. Our assumption of constancy is accurate for those components of the 1/ f noise at frequencies below the inverse of the gate time, around 100 MHz. The 1/ f noise at this frequency and higher is not accurately modeled in this way, but at these frequencies we believe that other sources of high-frequency noise (that is, white noise from the resistances in our circuit) become more important than 1/ f noise. In our work, these noise sources are modeled separately using a quantum bath; see [13] for these calculations.
Our estimates for T 1 and T 2 [16] times indicate that we can expect a very long coherence time in parking, O(1s), and a very short dephasing time, 10 ns, at the measurement position. In the portal region, we calculate coherence times of the order of hundreds of microseconds. As we shall see, except in the measurement processes, our gates are designed to not go lower than the upper limit of the portal. Thus, in principle, for gates of duration of 20-30 ns, the fidelity should not be compromised by decoherence much more strongly than by the imperfections explored in these figures: unwanted shifts in the applied fluxes due to low-frequency noise, and temporal shifts between c and pulses.
Those shifts are assumed uncorrelated and constant during a single gate operation, but random from one 'shot' to another. We model the noise as a normalized Gaussian probability distribution, µ(x) ≡ 1 N e −x 2 /2σ 2 , which we assume to have (at 1 Hz) 6 µ 0 (flux shifts) and 6 ps (time shifts) as its root mean square (rms) deviation, σ , such that approximately 90% of the distribution is found between the values ±10 µ 0 (flux shifts) and ±10 ps (time shifts); we believe that this quality of control will be readily achievable in the laboratory in the coming years. Indeed, recent experimental results [19] indicate that our assumption of having 6 µ 0 as rms deviation of the 1/ f flux noise at 1 Hz is already achievable for Josephson qubits.
We note that by using the 1/ f rms amplitude at 1 Hz in this model, we are implicitly assuming that the 1/ f noise is cut off below this frequency. In fact, the 1/ f spectrum goes much lower (at least three orders of magnitude lower in [19] and in other similar previous studies); but we may assume that in quantum computer operation, qubits are frequently taken offline and recalibrated. Assuming that this recalibration takes place once per second gives the 1 Hz cut-off. Our assumption of Gaussian statistics merely embodies the expectation that the 1/ f noise arises as the summed effect of many independent fluctuators; this is also borne out by the traces of [19]. Table 1 summarizes the average fidelity obtained for the gates discussed considering each channel noise separately: F = dδxµ(δx)F(δx), where F is given by equation (7). As one can see, the minimal expected fidelity found was 99.47% for the Hadamard gate.  For the phase, Hadamard and controlled-Z gates, we present the operator-sum representation [20] of the system superoperator where the operators {E k } are the operation elements for the quantum operation E. Having the set of operators {E k } permits the characterization of the noise present during the system evolution.
We find that the model noise considered in this work indicates that the effective noise for the physical implementation of our universal set of quantum gates is heavily biased, i.e., we have found that the effect of phase noise is at least one order of magnitude higher than bit-flip noise and leakage processes. New strategies for fault-tolerant computation have recently been worked out which take advantage of this sort of biasing in the noise [21].

Measurement gates
We present two distinct measurement gates, both very useful for effective quantum error correction and universal quantum computation. As one can see from their matrix representations equations (8) and (9), by measurement gates we mean the unitary operations that are used to prepare the state for the final projection process (measurement). The first gate is a measurement in the standard 0/1 basis. That is, starting from parking, we are to distinguish whether the system is in its ground state (0 quanta in the transmission line mode) or its first excited state (1 quantum in the transmission line mode). This gate, shown as the first plot of figure 6, works by performing an adiabatic evolution of the states |0 and |1 from the parking point to the measurement point, maintaining the bias flux at a constant value off the S line ( = 0). So, through this adiabatic transformation, the qubit states evolve from the configuration |S, 0 (ground state) and |S, 1 (first excited state), to |L , 0 and |R, 0 (where the first label corresponds to the bare qubit states, and the second to those of the transmission line). The matrix representation equation (8) shows that we end with an irrelevant relative phase θ 01 between the |0 and |1 states. The sources of loss of fidelity are leakage at the avoided crossing gap, and 0/1 transitions at the portal. For the path shown, the probabilities of observing leakage and 0/1 transitions are 3 × 10 −4 and 6 × 10 −3 %, respectively. As a result, the expected total net gate fidelity was found: dδ c dδ µ(δ c )µ(δ )F(δ c , δ ) ∼ 99.99%. The second measurement gate is in the conjugate basis. It evolves equal superpositions of |0 and |1 states at the parking point, |± = 1 √ 2 (|0 ± |1 ), to the final states |L , 0 and |R, 0 , respectively, at the measurement point. Since we know an adiabatic evolution would preserve the amplitude of probability of each state in the superposition, we have to design a non-adiabatic pulse in order to implement this +/− measurement. Indeed, by passing at an appropriate rate through the portal region, we achieve the desired transformation. The second plot of figure 6 presents the proposed gate. It consists of two distinct parts: the first occurs up to t 7 ns, when an adiabatic c pulse is applied to the qubit. This part is done slowly to minimize the leakage when the system passes through the avoided crossing gap; it is also tuned so that the correct relative phase is accumulated between the |0 and |1 states, with respect to a reference phase. The second part of the pulse, t 7 ns, causes the qubit to undergo a non-adiabatic evolution through the portal, producing the transitions needed to implement the gate. The bias flux is maintained constant during the whole gate, and it is held as close as is practical to the S linein our calculation we take = 30 µ 0 . The matrix representation equation (9) reveals the ideal transformation, which ends with an irrelevant relative phase between the ground and first excited states. For the path proposed, the probability of leakage is found to be 0.03%. Because the qubit has to stay, for an appreciable amount of time in a region of very fast relative phase accumulation, the main type of noise is a phase noise. The expected total net fidelity for this gate is: dδ c dδ µ(δ c )µ(δ )F(δ c , δ ) ∼ 99.8%.

Phase gate
The phase gate, figure 7, is the simplest operation of our set. In order to accumulate the desired phase, i.e. θ z ( t) ≡ t 0 ω 01 (t) dt, it is sufficient to adiabatically bring the system from the parking point to a position where the ω 01 frequency deviates a few hundreds of MHz from the reference frame; then one just has to wait the appropriate amount of time, corresponding to the desired phase accumulation (see the inset of figure 7). Since this position can be reached without passing through the avoided crossing, the leakage probability is extremely low; we calculate a leakage probability of 10 −7 %. Since we remain essentially within a very stable parking regime, the gate, for any value of θ z , is very insensitivity to fluctuations of both the applied fluxes. Thus, the gate can be performed with a very high total net fidelity: dδ c dδ µ(δ c )µ(δ )F(δ c , δ ) ∼ 99.999%. Since the operation is designed to perform an adiabatic evolution of the system, the 0/1 transitions are almost completely suppressed. Consequently, the main source of (small) fidelity loss can be characterized as a phase noise. Indeed, if we take a look at the operator-sum representation of the system superoperator, equation (12), we obtain the following decomposition: E 0 ≈ 0.99999σ z , E 1 ≈ 0.0031, E 2 = 0 and E 3 = 0, which clearly indicates that we have just one form of noise.

The Hadamard gate
Completing our universal set of one-qubit gates is the Hadamard gate, figure 8. Its implementation is more complex than the previous gates, since it requires the synchronization of the two parameters, c and . As one can see, the gate can be basically divided into three parts: the first, t 10 ns, evolves the qubit in an adiabatic process from the memory state to the portal. During this process, only control flux is changed; once the start of the portal is reached, the second part of the gate is applied. It consists of a non-adiabatic pulse performed by the bias flux . This pulse is responsible for very quickly moving the qubit from one side of the S line to another. This process, 10 ns t 14 ns, is designed to create the desired superpositions between the states |0 and |1 . In order to require reasonable rise-times for this pulse, we have to be as far as possible from the S line. Nevertheless, this excursion in the bias flux is limited due to the very small gap ω 12 region seen in the flux space, see figure 5. Finally, the last step, t 14 ns, adiabatically brings the qubit back to its initial position-again, only c is changed in this regime.
As highlighted in figure 5, the path taken in the flux space due to this fluxes c (t) and (t) composition has the form of a 'U'-shape. This is a very useful feature, since working on both sides of the S line, one can expect to correct in the second half of the pulse some of the errors accumulated during the first. Indeed, because of the symmetry around the S line, in the first order of approximation, the errors due to the shifts in the bias flux, δ , should be canceled out, since at this order ω 01 (± + δ ) = ω 01 (| |) ± αδ . However, for this to be true, a precise synchronization of the c (t) and (t) pulses is required, in order to obtain a symmetric path. For the paths in figure 5, the number of dots indicates the rate as the system evolves: the fewer the dots, the faster the time flux rate; one can clearly see the other strategies used to optimize the gate: at the parking regime, we can perform a very fast evolution, since the ω 01 and ω 12 are very big (∼3.1 GHz); once we have reached the avoided crossing position, in order to avoid leakage, we slow down the evolution. However, this regime also coincides with a fast phase accumulation rate, so we try to not spend too much time at this region. Thus, we find a trade-off that has to be worked out during the optimization of the gate. For the path shown, we expect to observe a probability of leakage of 0.12%.
The second plot of figure 8 shows the fidelity as a function of δ c , δ , and the flux desynchronization δt. It turns out that the gate proposed is very insensitive to fluctuations in δ c . This occurs because we have designed the gate to explore a 'sweet' spot in control flux. As one can see in figure 2, the region close to c = 1.447 0 presents a first-order insensitive point to fluctuations of c for several values of . Thus, since this point is at the portal, we have chosen it to 'sit' the qubit while we perform the bias pulse. The average fidelity considering each noise channel separately is presented in table 1. Because the δ and δt fluctuations essentially have the same effect, changing the relative phase of the two computational basis states, their loss of fidelity is found to be very similar. Thus, we see that the phase noise again plays a major role in the loss of fidelity. The complete noise characterization is obtained through the operator-sum representation of the system superoperator: The net fidelity considering the effects of all noise channels together is given by: dδ c dδ dδtµ(δ c )µ(δ )µ(δt)F(δ c , δ , δt) ∼ 99.46%.

The two-qubit system
The two-qubit system we have used to simulate the two-qubit gate is sketched in figure 1. This layout preserves the same structures present for the one-qubit system, i.e. the transmission line, the readout SQUIDs and the flux lines; in addition, it has a qubit-qubit interaction that is assumed to arise due to the mutual inductance between the two big loops [14] (other qubitqubit coupling implementations using tunable interactions are demonstrated in [22,23]). The qubit-qubit mutual inductance is considered to be a small parameter, such that a first-order perturbation theory is expected to give a very fair description of the system dynamics. Thus, following the procedure adopted to derive equation (5), we obtained the two-qubit system Hamiltonian given by (see the appendix) H A,B represent the single-qubit Hamiltonians, equation (5), of the qubits A and B, respectively. The qubits are assumed to have identical bare qubits, but with different fundamental-mode transmission lines, ω A T = ω B T . This choice was made in order to avoid possible double excitation of the transmission lines due to the transfer of one quantum of energy from one transmission line to another. In addition, because we would like to use the same master clock to track the dynamics of both qubits, we had to impose a rational relation between ω A T and ω B T . In our calculations, we have assumed that ω A T /ω B T = 3/4. The system interaction Hamiltonian H I is given by equation (14). Its terms may be understood as arising due to 'classical' and 'quantum' magnetic field components, in the following sense: for the first two terms, we clearly can identify them having the form of a bias field B applied to each qubit (compare with the bias term of equation (4)). As expected, the bias field B applied to one qubit, let us say qubit i, depends on the control parameter of the other, qubit j. This happens because the external circulating current of qubit j creates a magnetic field, B ( j c ), seen by qubit i. Nevertheless, for a given j c , the qubit i feels the same bias field B ( j c ) whatever the quantum state of qubit j is. Thus, B should be seen as simply the result of Faraday's law applied to a classical circuit: the field j c induces a persistent current in the circuit j, which, in turn, creates a magnetic field B ( j c ) seen by the other qubit. Not surprisingly, the first two terms of H I cannot generate qubit-qubit entanglement. The dot-dashed curve in figure 9 shows B as a function of the control flux. Because of the presence of the Josephson junctions (nonlinear inductances), B has nonlinear behavior for some values of c (double-well regime). As a consequence of the existence of the magnetic field B , for the two-qubit system, the S line of qubit i is shifted to the lines i = k 0 ± B ( j c ) (where k is an integer, and the choice of the sign depends on which loop the flux B is threaded through). Thus, we refer to those 'new' qubit symmetric lines as the S lines. This is a very important effect and must be taken into account when performing one-qubit operations for the two-qubit system.
The last term of equation (14) describes a magnetic field seen by one qubit, determined by the quantum state of the other. The result is a qubit-qubit coupling of the formσ A z ⊗σ B z , which is responsible for generating qubit-qubit entanglement in the system. As one can see in figure 9, even though the physical qubit-qubit inductive coupling is assumed fixed, the effective interaction J is tunable as a function of both control flux parameters. In fact, as shown in the second plot of figure 9, the coupling J never reaches values higher than 450 MHz when one qubit is maintained deeply in the single-well regime (solid curve), rather than 5 GHz observed when both qubits are at the measurement point (dashed curve). This is a direct manifestation of the change of the system state character as a function of control parameter.
In addition, because the transmission line of qubit A sees the qubit B only via interaction with its bare qubit, and the interaction between the bare qubits is in fact small, it turns out that the transmission line states are very stable under the changes of control parameters of the other qubit. Indeed, since the qubit-qubit interaction equation (14) is capable neither of moving the qubits from the single-to the double-well potential regime (for that, terms of the formσ x are needed) nor of changing the quantum number of the transmission lines, once one qubit is parked, the two-qubit system eigenstates stay 'frozen' at the subspace in which the parked qubit Table 2. The expected one-qubit gate fidelities, as a function of unwanted shifts in δ c , δ and δt, when performed on the two-qubit system using the same schemes described in section 4. We observe a small difference between the oneand two-qubit system fidelities only for the +/− measurement gate. eigenstates are close to the states |S, 0 and |S, 1 . Thus, parking one qubit imposes a selection rule on the system that leads to a further reduction, below the already small value of J , of the effective qubit-qubit coupling. In fact, our calculations of the interaction matrix elements where the states i, j, l, m are the 0-1 one-qubit states) show that they never reach values higher than 3 MHz when one of the qubits is parked.
As a result, parking leads to a very effective way to decouple the qubits, thus allowing one to perform one-qubit gates for the two-qubit system discussed. By simply taking into account of the fact that the qubit S lines in the flux space are moved to the S lines, we can use exactly the same shaped pulse schemes presented previously for the one-qubit system, in order to perform the universal set of one-qubit gates. Indeed, our simulations showed that one can expect the same probability of leakage previously reported, and virtually the same expected fidelities, see table 2.

Controlled-Z gate
The two-qubit gate proposed is a gate in the equivalence class [24] of the controlled-Z gate, i.e. it only differs from the controlled-Z gate by local operations to the qubits A and B. The designed gate, the first plot of figure 10, involves an adiabatic evolution of the system states, in order to accumulate the correct relative phases. Both c and of qubits A and B are used to perform the operation. In order to obtain the strongest qubit-qubit interaction, and thus the shortest gate possible, the control fluxes of both qubits are changed simultaneously, A c (t) = B c (t). Also the bias fluxes change identically and, as for the Hadamard scheme, we work on both sides of the S line. The desired unitary transformation has the matrix representation in the two-qubit eigenstate basis {|0 , |1 , |2 , |3 } given by The relative phase is the pertinent parameter of this gate. Since the local invariants [24] of the gate U Z Z are determined by and those of the controlled-Z gate are given by G 1 = 0 and G 2 = 1, we end with the specific condition θ Z Z = π(1 + 2k) for the relative phase.
Since the evolution is done adiabatically, we can write the relative phase as where E i is the instantaneous eigenenergy of the state i. Thus, it is evident that θ Z Z only has non-negligible values when the term Jσ A z ⊗σ B z of equation (14) becomes appreciable. However, as shown in figure 10, to reach that region we have to pass through a regime of a very small gap (∼300 MHz) between the |3 and |4 states (minimum gap between the computational basis {|0 , |1 , |2 , |3 } and the rest of the spectrum of the system). Because of that, our shaped dc gate has a total duration of 60 ns in order to avoid leakage during the evolution. In addition, in order to explore 'sweet' spots of operations, the designed c (t) pulse is not completely flat when passing through the minimum gap region (see figures 10(b) and (c)).
The gate is performed in the following way: starting from the memory state, both control fluxes are changed to the region around c ≈ 1.453 0 . Then, by changing the bias fluxes A,B , the qubits are slowly moved through the region of minimum gap ω 34 , passing from one side of the S line to the other. At this stage, the parameter θ Z Z becomes appreciable (hundreds of MHz). The gate is designed to spend the right amount of time, such that the final evolution will satisfy the necessary condition for θ Z Z . Once the other side of the S line is reached, the qubits are brought back to the memory state. As observed for the Hadamard gate, the two-qubit gate proposed also has a 'U'-shape in the flux space.
The probability of leakage due to this process is ∼ 0.04%. Since the gaps between the computational states are bigger than 700 MHz, unwanted transitions are strongly suppressed. Therefore, phase noise is the main source of the loss of fidelity. As for one-bit gates, this phase noise arises from physical 1/ f noise and pulse timing variations; we assume their magnitudes to be the same as for the one-qubit case, with no correlations between the two qubits. We then calculate (assuming, again, Gaussian noise statistics) the operator-sum representation of the system superoperator: where The total net fidelity considering the effects of all noise channels is given by dδ c dδ dδ tµ(δ c )µ(δ )µ(δ t)F(δ c , δ , δ t) ∼ 99.67%.
The energy splitting ω 34 gives a limit on the rate at which the system can evolve without appreciable leakage. The θ Z Z plot gives the phase accumulation rate of the principal parameter of the gate, equation (16). The path (a) taken in the flux space is illustrated with dots in (b) and (c). The number of dots indicates the time spent when passing through that region.

Conclusion
Obviously, it is hoped that the results of this paper will be of direct relevance for experiments on qubits being performed at IBM, as well as being of general relevance to the problem of precision quantum control in any flux qubit system. The noise mechanisms analyzed here-magnetic 1/ f noise, Johnson noise of circuit resistances and timing errors in pulse channels-are a mixture of fundamental and practical sources of error that are currently understood in the laboratory. Thus, the quantitative fact that the values of the gate infidelity (one minus the fidelity) are at the 1% level and below, is the major result of this paper.
One might ask, given that this paper's main claim to importance is quantitative, whether the approximations made in the quantum model are actually adequate to give sub-1% accuracy. Perhaps the fundamental circuit Hamiltonian may be accepted as being very accurate; but can it be that the sequence of subsequent approximations, namely (i) the representation of the transmission line by a single oscillator mode, (ii) the dimension-reduction resulting from the Born-Oppenheimer approximation, and (iii) the final spectral truncation, involving an interpolated oscillator-displacement parameter, resulting in our four-level model, have the desired sub-1% accuracy? We can see, for example, in figure 4, that the spectral truncation results in changes in the absolute value of the energy eigenvalues, in some places by as much as 7%.
Despite this, we believe that our infidelity numbers are nevertheless quite sound, so that if we calculate an infidelity here of 1%, it may be, in a more accurate calculation, actually equal to 0.9 or 1.1%, but not much different from that. We have a few reasons for saying this. First, for our gate designs, it is generally not necessary that an energy gap be exactly a certain value at a particular point on the flux axis; it is more important that, somewhere in the near vicinity of a particular flux value, the energy gap has a certain size. This will be true even for a model Hamiltonian of moderate accuracy.
Traced to a deeper level, the successful functioning of our gates depends primarily on the following general features of the model. (i) Since many of the evolutions we consider are adiabatic, it is important merely that some integrated properties of the energy eigenvalues over some paths in parameter space be correct. This is fairly easily achieved by a model of moderate accuracy. (ii) In a few cases, basis-changing, non-adiabatic evolutions are important. These are achieved in just three ways: by passing through the portal, by crossing the S line, and (unintentionally) by moving the system in and out of parking. Our model obviously contains all of these features, with trends in the size of gaps that certainly track those of an exact calculation very closely. Furthermore, the general structure of the low-lying eigenvectors is captured in our truncation, meaning that trends in the matrix elements that determine the magnitude of the Landau-Zener tunneling effects are also well represented.
So, we conclude that our qubit, as currently understood, should be capable of gate operation at the 1% noise level. Unfortunately, in the lab, in any given day (or month) there are 'bugs' in the experiment that cause the qubit, often for unexplained reasons, to have fidelities much worse than what we estimate here. But we nevertheless hope that these calculations will have significant practical value. Our gate set is universal, in one of two different ways, in fact: the 0/1 measurement/preparation gates, Hadamard, controlled-Z , and one-qubit phase gates form one well-known universal set [20]; it is less well-known that, with more overhead, Hadamard can be replaced by the +/− measurement/preparation gates [25].
So, can a 'debugged' IBM qubit be used soon for universal quantum computation? The answer is, in our opinion, ultimately yes. The answer would certainly be no if the noise threshold for fault-tolerant quantum computation were in the neighborhood of the oft-quoted value of 10 −5 [26]- [28]. It is not inconceivable for the experiment to get to these values someday, since we find that the infidelities decrease much faster than linearly with the assumed noise levels.
(To get to 10 −5 we would need to get to the very daunting levels of 100 n 0 at 1 Hz for the 1/ f noise amplitudes and 100 fs for timing accuracies; there is optimism that both of these numbers are ultimately attainable.) Fortunately, while 10 −5 was the threshold as it was understood 10 years ago, much recent work shows that with good designs, much higher thresholds are possible [29,30]. According to [30], 1% is in fact on the high end of the noise levels for which fault tolerance may be possible.
But the final answer to the question of whether fault tolerance is possible will require resolving a number of further uncertainties. One effect not included here is the fact that noise will be correlated between qubits that are physically nearby, due to ordinary electrical cross-talk, say. This is certainly deleterious for noise thresholds. However, a great virtue of our parking scheme is that, while parked, qubits are quite immune to electrical noise from nearby sources. Another point that might count in our favor is that, as seen in our operator-sum expansions, the noise has a definite structure (i.e. it is mostly phase noise in most cases), while most analyses of noise thresholds have assumed worst-case, structureless noise. On the other hand, these special noise sources, to the extent that they arise from the physical 1/ f noise, are substantially correlated in time. While fault tolerance is known to survive in the presence of such 'non-Markovian' sources [31], it is not clear what such an effect does to the numerical values of the threshold. Finally, there are purely 'architectural' considerations, e.g. the theoretical analyses assumed that a qubit can be moved to the proximity of any other one without cost. This will surely not be true in any real Josephson architecture.
To summarize, we have demonstrated that using only low-bandwidth electrical pulses, a universal set of high-fidelity quantum gates can be achieved for an oscillator-stabilized Josephson qubit, with a required 'clock time' around 60 ns, and achievable fidelities above 99% per gate operation. Essential to these results are the existence of 'parking' made possible by coupling of the flux qubit to a transmission line, the possibility of mostly adiabatic control, and the availability of two robust basis-changing effects obtained by crossing the 'portal' and the 'S line'. Time will tell if these physical elements make it possible, in the face of the many sources of noise that are present in the solid-state environment, to do large-scale quantum computation.

Appendix A. Hamiltonian derivation
In this appendix, we supply derivations of the Hamiltonians presented in the main text, equations (4) and (14). The main assumption in the derivation is the use of the first-order perturbation theory to treat terms of the system potential. As presented in figure 4 and previously discussed, performing the procedure described here, we could arrive at a very good description for the lowest states of the IBM qubit.

A.1. Brief review of BKD
BKD [16] have introduced a universal method for analyzing any electrical circuit containing Josephson junctions, provided only that its elements can be represented by lumped elements. The methodology can be summarized in a few steps: a network graph is written for the circuit. A network graph is simply a drawing of the circuit where each two-terminal element (inductor, capacitor, Josephson junction, etc) is represented as an oriented labeled branch connecting two nodes. Then, a tree of the network graph (a subgraph that does not contain any loops) is chosen using the following criteria: the tree has to contain all of the capacitors in the system, no resistors 26 or external impedances, no current sources, no Josephson junctions and as few linear inductors as possible. The branches that do not belong to the tree are called chords.
Associated with the tree chosen, there are the so-called sub-loop matrices F X Y , where the label X represents the tree elements (X = C for capacitors and X = K for the tree inductors) and Y the chord branches (Y = J for Josephson junctions, Y = L for linear inductors, Y = R for shunt resistors, Y = Z for external impedances and Y = B for bias current sources). The subloop matrices have entries −1, 0 or 1, and give information about the interconnections in the circuit, determining which tree branches X are present in which loop defined by the chords Y . The entries in F X Y matrix are found as follows: if the corresponding tree element X i does not belong to the loop defined by the corresponding chord Y j , its entry F i, j X Y is zero. If it belongs to the loop but has opposite orientation to Y i , we have F i, j X Y = +1. Finally, if it belongs and has the same orientation, we have F i, j X Y = −1. These steps give an algorithm to encode the topology of the system in a matrix representation. In addition, the formalism assumes that all capacitors should be considered to be in parallel with a Josephson junction, even if it is one with zero critical current.
The physics of the circuit is introduced by imposing the Kirchhoff laws at each node of the network, and defining the electrical characteristics of each branch type as Here, the diagonal matrix I c contains the critical currents of the junctions, and sin ϕ is the vector (sinϕ 1 , sinϕ 2 , . . . , sinϕ N J ). The phase ϕ i represents the superconducting phase difference across the junction i. The linear capacitors are described by equation (A.2), with their capacitance values given by the entries of the diagonal matrix C. The junction resistors are assumed to follow Ohm's law, equation (A.3), where R is the real and diagonal shunt resistance matrix. The impedances are described by equation (A.4), which relates the Fourier transforms of the current and voltage. Z(ω) is the diagonal impedance matrix. Since linear inductors can be tree branches as well as chords, we have to distinguish them to apply the network graph theory. Therefore, the inductance matrix must be organized in the block form shown in equation (A.5), where L is the inductance matrix of the tree inductors, L K is that for the chord inductors, and L L K represents the mutual inductances between the tree and chord inductors. BKD arrived at the system Hamiltonian where the external applied fluxes and current are respectively represented byΦ  (66) in [16]. Thus, BKD maps the circuit dynamics to that of a massive particle in a potential, whose masses and degrees of freedom are associated with the system capacitances. The quantization of the system is introduced by imposing the canonical commutation relation for the variables of charge, Q C , and phase ϕ: [ 0 2π ϕ i , Q C; j ] = ihδ i j .

A.2. IBM qubit network graph
The L matrices are denoted: In addition, the capacitances and critical currents of the junctions are taken to be C = 10 fF and I c = 1.3 µA.

A.3. Bare qubit Hamiltonian
Our aim in this section is to present the steps used to derive the approximate Hamiltonian for the bare qubit. In order to reach the appropriate conditions to perform a first-order perturbation theory to treat the coupling between the qubit and the bias flux line, we firstly have to identify the 'slow' and 'fast' coordinates of the system. The slow coordinate is that which joins the two minima when the system presents a double-well structure, or that which has the slowest curvature in the single-well regime. The fast coordinates are those in which the potential rises very steeply, such that the system dynamics is frozen into the ground state along these directions.
The procedure to identify the slow and fast degrees of freedom starts from the exact system potential equation (2) and a unitary transformation R that diagonalizes M 0 (M 0 is a real and symmetric matrix [16]). Using R we can decouple the quadratic form of the potential by transforming the system coordinate into the new coordinates n ≡ Rϕ. Thus, from equation (2), we obtain (ignoring the term due to external sources of current) where we have defined the diagonal matrix λ ≡ RM 0 R T , whose entries are the exact eigenvalues Because of L 1 −2M 15 2(L 3 +M 35 −M 15 ) ≈ 0.024 and the fact that we expect to work with magnetic fluxes of the order of only a few 0 s, we have that, to good approximation, the system potential is only a function of the difference of the magnetic flux in the two large loops, defined as . This property is a direct manifestation of fact that the designed bare qubit has a 'gradiometer' structure. Even though this correction is expected to be very small, we shall consider it during the whole procedure we are going to describe here. The new coordinate system allows us to easily identify, by inspection of each term of equation (A.13), the following symmetry in the potential U n,˜ c , , p = U (n 1 , n 2 , n 3 , c , ) = U (−n 1 , −n 2 , n 3 , c , − ) .
(A. 16) This shows that the system potential presents a definite symmetry in the plane determined by the directions {n 1 , n 2 }, when the system goes from + to − (if we simultaneously adjust the control flux˜ c to maintain the new flux coordinate c , equation (A.15), constant). In addition, equation (A.16) also permits us to determine the conditions for which the system presents the same physical potential U (n 1 , n 2 , n 3 , c , ) = U n 1 , n 2 , n 3 , c , + γ , where γ is an irrelevant constant. In fact, we find that whenever the bias flux difference = − is an integer multiple of the flux quantum, let us say = (k 1 − k 2 − 2k 3 ) 0 (each k i is any integer), equation However, since λ 2 /λ 3 ≈ 0.059 and L −1 j;i /λ 3 ≈ 0.061 (i.e. the potential is much steeper in the direction s 3 than in the others), we have the global minimum of the potential occurring at the position s 3 ≈ 0. As a result, the bare qubit low-level dynamics can be considered frozen at the plane s 3 ≈ 0.
Another important result is derived from the fact that the low-level system dynamics is only governed by the degrees of freedom s 1 and s 2 : the system potential has a symmetry line (S line) at = 0 in the flux space. Indeed, as one can see, when = 0, U (s, c , ) has a perfect symmetry around the origin in the {s 1 , s 2 } plane. Because of that, the low-level system wavefunctions have definite parity symmetry at = 0. In addition, because the same physical potential is found when is changed by an integer multiple of 0 , an S line occurs whenever is an integer multiple of 0 .
Unfortunately, because λ 2 ∼ L −1 j;i , it is not simple, as it was for the direction s 3 , to determine the soft and fast directions in the plane defined by s 3 ≈ 0. However, from the symmetry of the S line one can see that, if the potential has a minimum in s min 1 = 0 and s min 2 = 0, the potential presents a symmetric double-well structure (with maximum at s 1 = s 2 = 0) in the direction connecting the minimum points. This direction is expected to be the slow coordinate of the system, while its perpendicular direction should determine the fast degree of freedom. Consequently, by finding the minima positions in the plane {s 1 , s 2 }, one can determine those directions. Following this procedure, we perform one more rotation, q 1 = s 1 cos θ − s 2 sin θ, q 2 = s 1 sin θ + s 2 cos θ and q 3 = s 3 , such that the direction q 1 connects the minima. When the potential does not have a double-well structure, we define the last rotation so that the potential curvature in the direction q 1 is the smallest.
Thus, after appropriate transformations, we end with phase and flux coordinate systems in which the symmetries presented by the system are much more clearly stated. In addition, the slow, q 1 , and fast, {q 2 , q 3 }, degrees of freedom of the system are also identified. As already mentioned, because the potential is very steep in the fast directions, the low-level system dynamics is frozen into the ground state along these directions. Following the Born-Oppenheimer approach developed in [13], these directions can be traced out and their effects incorporated as small corrections to the remaining slow-coordinate potential energy.
At last, following all the steps described above, we are now in a position to construct term by term Hamiltonian equation (5): • No bias flux line ( = 0): As the system potential has a perfect symmetry around the origin, the bare qubit wavefunctions also have a parity symmetry. Thus, the ground and first excited states are expected to be symmetric and antisymmetric with respect to reflection around the origin. As described in the main text, the system potential is found in a doublewell structure for small values of c . There, the bare qubit eigenstates can be understood as symmetrical and antisymmetrical equal superpositions of the classical orbital states of the left and right wells. For the single-well regime, c 1.5 0 , the bare qubit states are associated with those of a harmonic oscillator. Thus, in this representation, the Hamiltonian of the bare qubit without bias flux line can be written as − 1 2 ( c )σ x , where ( c ) (shown in figure 3) is the ground and first excited state energy splitting.
• The bias flux term: As previously discussed, the bias flux is responsible for breaking the perfect potential symmetry observed at the S line. However, the system presents, as stated in equation (A.16), another important symmetry in going from one side of the S line to another (i.e. passing from − to + ): the system potential at < 0 is identical to that in > 0 by a mirror reflection (a π -rotation for the full multi-dimensional potential) at the origin. In the double-well regime that corresponds to the interchange of the left and right states under the mentioned transformation. The bias term of the Hamiltonian equation (5) is obtained applying a first-order perturbation theory to treat the bare qubitbias flux line coupling. Indeed, because the potential term representing their coupling, U q ≈ q 1 π √ (2/3)λ 2 sin θ, is very small, for typical values of bias flux used in this work ( = O(m 0 )), compared to the other potential terms, the first-order approximation can be done in a very controllable way. Because of the symmetry of the bare qubit eigenstates, the interaction matrix elements i |U q | j (where = S, A are the symmetrical and antisymmetrical bare qubit eigenstates) are nonzero only for those connecting the ground and excited states. Thus, the bias term has the form: 1 2

A.4. Bare qubit-transmission line coupling
Another structure present in the IBM qubit is a high-quality superconducting transmission line (figure 1). We model the fundamental mode of the transmission line as a simple LC circuit of definite frequency ω T coupled to the bare qubit. We assume that its characteristic impedance and inductance are given by Z 0 = √ L T /C T = 110 and L T = 5.6 nH, respectively, and the bare qubit-transmission line mutual inductance equals M qT = 200 pH. In order to find the slow coordinate directions, the pertinent terms of the interaction potential are given by As one can see, the linear terms in equation (A.23) are due to the 'frozen' values of the fastest degree of freedom of each qubits. Moreover, the linear term of the slow coordinate of qubit i arises due to the frozen coordinate of qubit j and vice versa. The interaction matrix elements due to these terms are exactly that of the one-qubit bias term, leading to a Hamiltonian form: Finally, the term q A 1 q B 1 of equation (A.23) only connects states of different parity in both qubits. Thus, this term leads to the form: