Demonstration of nonstoquastic Hamiltonian in coupled superconducting flux qubits

TL;DR

This paper demonstrates a nonstoquastic Hamiltonian in coupled superconducting flux qubits by coupling via charge and flux degrees of freedom, potentially improving quantum annealing and simulation capabilities.

Contribution

It introduces a novel coupling scheme for flux qubits that achieves a nonstoquastic Hamiltonian, expanding the possibilities for quantum annealing and simulation.

Findings

Charge coupling manifests as YY interaction

Destructive interference observed in quantum oscillations

Hamiltonian shown to be nonstoquastic over a wide parameter range

Abstract

Quantum annealing (QA) is a heuristic algorithm for finding low-energy configurations of a system, with applications in optimization, machine learning, and quantum simulation. Up to now, all implementations of QA have been limited to qubits coupled via a single degree of freedom. This gives rise to a stoquastic Hamiltonian that has no sign problem in quantum Monte Carlo (QMC) simulations. In this paper, we report implementation and measurements of two superconducting flux qubits coupled via two canonically conjugate degrees of freedom (charge and flux) to achieve a nonstoquastic Hamiltonian. Such coupling can enhance performance of QA processors, extend the range of quantum simulations. We perform microwave spectroscopy to extract circuit parameters and show that the charge coupling manifests itself as a YY interaction in the computational basis. We observe destructive interference in quantum coherent oscillations between the computational basis states of the two-qubit system. Finally, we show that the extracted Hamiltonian is nonstoquastic over a wide range of parameters.

Equations

$$H=-\frac{1}{2}\sum _i{\Delta }_i{\sigma }_i^x+\sum _i{h}_i{\sigma }_i^z+\sum _{i<j}{J}_{ij}{\sigma }_i^z{\sigma }_j^z,$$

$$H$$

$$H=-\frac{{\Delta }_1}{2}{\sigma }_1^x-\frac{{\Delta }_2}{2}{\sigma }_2^x+h_1{\sigma }_1^z+h_2{\sigma }_2^z+\sum _{\alpha ,\beta }{J}_{\alpha \beta }{\sigma }_1^{\alpha }{\sigma }_2^{\beta },$$

$$H=H_q+H_{\mathrm{cjj}}-E_J\cos \left(\frac{\pi {\Phi }_{\mathrm{cjj}}}{{\Phi }_0}\right)\cos \left(2\pi \frac{{\Phi }_q}{{\Phi }_0}\right).$$

$$H_{\kappa }=\frac{Q_{\kappa }^2}{2C_{\kappa }}+\frac{({\Phi }_{\kappa }-{\Phi }_{\kappa }^x{)}^2}{2L_{\kappa }},\kappa =\{q,\mathrm{cjj}\}.$$

$$U_{\mathrm{eff}}({\Phi }_q)={\min }_{{\Phi }_{\mathrm{cjj}}}U({\Phi }_q,{\Phi }_{\mathrm{cjj}}),$$

$$U({\Phi }_q,{\Phi }_{\mathrm{cjj}})=\sum _{\kappa =q,\mathrm{cjj}}\frac{({\Phi }_{\kappa }-{\Phi }_{\kappa }^x{)}^2}{2L_{\kappa }}-$$

$$E_J\cos \left(\frac{\pi {\Phi }_{\mathrm{cjj}}}{{\Phi }_0}\right)\cos \left(2\pi \frac{{\Phi }_q}{{\Phi }_0}\right).$$

$$H=\frac{Q_q^2}{2C_q}+U_{\mathrm{eff}}({\Phi }_q).$$

$$U_{\mathrm{eff}}({\Phi }_q)=\frac{({\Phi }_q-{\Phi }_q^x{)}^2}{2L_q}-E_J({\Phi }_{\mathrm{cjj}}^x)\cos \left(2\pi \frac{{\Phi }_q}{{\Phi }_0}\right),$$

$$E_J({\Phi }_{\mathrm{cjj}}^x)=\frac{{\Phi }_0{I}_c}{2\pi }\cos \left(\pi \frac{{\Phi }_{\mathrm{cjj}}^x}{{\Phi }_0}\right).$$

$$U_M=\frac{M_{12}}{L_1{L}_2}({\Phi }_{q,1}-{\Phi }_{q,1}^x)({\Phi }_{q,2}-{\Phi }_{q,2}^x),$$

$$K=\sum _i\frac{C_i{\dot{\Phi }}_{q,i}^2}{2}+\frac{C_{12}({\dot{\Phi }}_{q,2}-{\dot{\Phi }}_{q,1}{)}^2}{2},$$

$$\mathcal{L}=K-\sum _{i=1}^2{U}_{\mathrm{eff},i}({\Phi }_{q,i})-U_{12},$$

$$Q_1=(C_1+C_{12}){\dot{\Phi }}_{q,1}-C_{12}{\dot{\Phi }}_{q,2},$$

$$Q_2=(C_2+C_{12}){\dot{\Phi }}_{q,2}-C_{12}{\dot{\Phi }}_{q,1}.$$

$$H=\sum _{i=1}^2{H}_i+U_M+U_C.$$

$$H_i=\frac{Q_i^2}{2{\tilde{C}}_i}+U_{\mathrm{eff},i}({\Phi }_{q,i})$$

$${\tilde{C}}_1=C_1+\frac{C_{12}{C}_2}{C_2+C_{12}},$$

$${\tilde{C}}_2=C_2+\frac{C_{12}{C}_1}{C_1+C_{12}}.$$

$$U_C=\frac{C_{12}{Q}_1{Q}_2}{C_1{C}_2+(C_1+C_2)C_{12}}.$$

$$H_i=\sum _{\mu =1}^{N_i}{\epsilon }_{\mu }^i|{\chi }_{\mu }^i⟩⟨{\chi }_{\mu }^i|.$$

$$H=\sum _{\mu =1}^{N_1}{\epsilon }_{\mu }^1|{\chi }_{\mu }^1⟩⟨{\chi }_{\mu }^1|+\sum _{\nu =1}^{N_2}{\epsilon }_{\nu }^2|{\chi }_{\nu }^2⟩⟨{\chi }_{\mu }^2|+U_M+U_C.$$

$$|{\eta }_a⟩=\sum _{\mu }^{N_1}\sum _{\nu }^{N_2}{c}_{\mu \nu }^a|{\chi }_{\mu }^1\otimes {\chi }_{\nu }^2⟩,$$

$$c_{\mu \nu }^a=⟨{\chi }_{\mu }^1\otimes {\chi }_{\nu }^2|{\eta }_a⟩.$$

$$H=\sum _{a=1}^N{\epsilon }_a|{\eta }_a⟩⟨{\eta }_a|,$$

$$|\downarrow {⟩}_i$$

$$|\uparrow {⟩}_i$$

$$|{\eta }_a⟩\simeq \frac{1}{N_a}\sum _{\mu =1}^2\sum _{\nu =1}^2{c}_{\mu \nu }^a|{\chi }_{\mu }^1\otimes {\chi }_{\nu }^2⟩,$$

$$N_a=\sqrt{\sum _{\mu =1}^2\sum _{\nu =1}^2|c_{\mu \nu }^a{|}^2}.$$

$$|1⟩=|{\downarrow }_1\otimes {\downarrow }_2⟩,|2⟩=|{\downarrow }_1\otimes {\uparrow }_2⟩,$$

$$|3⟩=|{\uparrow }_1\otimes {\downarrow }_2⟩,|4⟩=|{\uparrow }_1\otimes {\uparrow }_2⟩.$$

$$\displaystyle H=\left[\begin{array}[]{cccc}J_{zz}-h_{1}-h_{2}&-\Delta_{2}/2&-\Delta_{1}/2&J_{xx}-J_{yy}\\ -\Delta_{2}/2&-J_{zz}-h_{1}+h_{2}&J_{xx}+J_{yy}&-\Delta_{1}/2\\ -\Delta_{1}/2&J_{xx}+J_{yy}&-J_{zz}+h_{1}-h_{2}&-\Delta_{2}/2\\ J_{xx}-J_{yy}&-\Delta_{1}/2&-\Delta_{2}/2&J_{zz}+h_{1}+h_{2}\end{array}\right],$$

Full text

I Introduction

Early generation processors Johnson11 for Hamiltonian-based quantum computation were designed to perform quantum annealing (QA), a heuristic algorithm for finding low-energy configurations of a system Kadowaki98 ; Farhi01 ; Santoro02 . With additional features, other applications such as quantum simulation Harris18 ; King18 and machine learning Mott17 ; Amin18 have become possible. Currently available large-scale quantum processors are made from a network of radio-frequency superconducting quantum interference devices (rf-SQUIDs) Johnson11 ; Harris09 ; Harris10 . Interaction between pairs of devices is realized through tunable magnetic coupling of their flux degrees of freedom. The low-energy dynamics of individual rf-SQUIDs are effectively captured with their two lowest-energy eigenstates, allowing one to approximate rf-SQUIDs as qubits, described by Pauli matrices $\sigma^{x,y,z}$. The computational basis states $|{\uparrow}\rangle$ and $|{\downarrow}\rangle$ (eigenfunctions of $\sigma^{z}$) correspond to directions of persistent current in the body of the rf-SQUID. This network implements the transverse field Ising model Hamiltonian:

[TABLE]

where $\Delta_{i}$ and $h_{i}$ are tunneling amplitude and energy bias of qubit $i$, respectively, and $J_{ij}$ is the magnetic coupling strength between qubits $i$ and $j$. Quantum annealing is performed by adjusting $\Delta_{i}\gg h_{i},J_{ij}$ at the beginning of the annealing process and gradually evolving until $\Delta_{i}\ll h_{i},J_{ij}$ at the end.

One important property of Hamiltonian (1) is that it is stoquastic, meaning that its equilibrium properties can be simulated using stochastic algorithms such as quantum Monte Carlo (QMC) Isakov2016 ; Amin17 ; Albash2018 . These algorithms typically operate in a local basis, in which all basis vectors are product states of individual qubits. Otherwise, representation of the basis vectors requires exponential resources. A Hamiltonian is stoquastic if there exists a local basis such that all its off-diagonal elements are real and nonpositive Bravyi08 ; Lidar18 ; Klassen18 . A positive off-diagonal element would cause negative transition probabilities, which cannot be simulated by stochastic processes. This issue is called sign problem in QMC Loh1990 . If a Hamiltonian is nonstoquastic, there should exist no local basis in which all off-diagonal elements are negative or zero. For some stoquastic Hamiltonians, finding the local transformation that cures the sign problem is by itself intractable Lidar18 . Therefore, proving that a Hamiltonian is nonstoquastic is also intractable. To achieve nonstoquasticity, additional interactions such as $\sigma_{i}^{x}\sigma_{j}^{x}$ or $\sigma_{i}^{y}\sigma_{j}^{y}$ are needed in Hamiltonian (1). Such interactions can enhance performance of QA Hormozi17 ; Nishimori17 , extend the range of quantum simulations Love14 , provide a path towards annealing-based universal quantum computation Aharonov ; Biamonte ; Jordan , and enable error supression in quantum annealing Marvian17 ; Jiang17 . Nonstoquasticity may also be achieved via nonadiabatic evolution Vinci .

II Hamiltonian

To implement a nonstoquastic Hamiltonian, we use two rf-SQUIDs coupled both inductively, through a tunable mutual inductance $M_{12}$ Harris09 , and capacitively, through a fixed capacitance $C_{12}$ as shown in Fig. 1(a). The Hamiltonian of the system is

[TABLE]

where $Q_{i}$ and $\Phi_{\text{q},i}$ are charge and flux variables that satisfy the commutation relation $[\Phi_{\text{q},i},Q_{j}]=i\hbar\,\delta_{ij}$, $\Phi_{\text{q},i}^{x}$ and $\Phi_{{\rm cjj},i}^{x}$ are external fluxes, and $\Phi_{0}=\pi\hbar/e$ is the flux quantum. Each rf-SQUID is characterized by its capacitance $C_{i}$, inductance $L_{i}$, and tunable Josephson energy $E_{Ji}(\Phi_{{\rm cjj},i}^{x})\approx(\Phi_{0}/2\pi)I_{ci}\cos(\pi\Phi_{{\rm cjj},i}^{x}/\Phi_{0})$ Harris10 (for a more detailed description see appendix A.1). The tunable mutual inductance $M_{12}$ is adjusted with the coupler control bias $\Phi_{\rm co}^{x}$. The renormalized capacitances are defined as $\widetilde{C}_{1(2)}=C_{1(2)}+C_{12}C_{2(1)}/(C_{2(1)}+C_{12})$. Capacitive couping has been proposed and analyzed theoretically for flux qubits in Levitov01 ; Bruder05 ; Satoh15 ; Kerman18 and is commonly used in transmon qubits google ; Kounalakis18 .

The potential energy of each rf-SQUID can have a double-well shape (see the two-qubit potential in Fig. 1(b)). The barrier height of the potential is controlled by $\Phi_{{\rm cjj},i}^{x}$, which tunes the tunneling amplitude $\Delta_{i}$, but also changes the persistent current. The potential is monostable when $\Phi_{{\rm cjj},i}^{x}=$ 0.5 $\Phi_{0}$. The flux bias $\Phi_{\text{q},i}^{x}$ adjusts the tilt of the potential, setting $h_{i}$. Both $\Phi_{\text{q},i}^{x}$ and $\Phi_{{\rm cjj},i}^{x}$ are controlled by high bandwidth coaxial lines, allowing for microwave operation and fast quench of the coherent dynamics. At the end of quench, which involves raising the tunneling barrier rapidly, the persistent current is measured via a shift register coupled to a microwave resonant readout Berkley2010 ; Whittaker16 .

II.1 Parameter extraction

Before characterizing the circuit parameters, we calibrate the tunable magnetic coupler, which provides the function $M_{12}(\Phi_{\rm co}^{x})$ as discussed in Ref. Harris09 . For the rest of the manuscript, we treat the coupler as a simple tunable mutual inductance, assuming dynamics of the coupler are significantly faster than single- and coupled-qubit dynamics. Next, we measure the persistent current of each qubit for a range of $\Phi_{\rm cjj}^{x}$ in a regime where its tunneling amplitude is negligible and the other qubit is kept monostable. We fit these measurements to a classical rf-SQUID model Harris10 and obtain $I_{c,1}=$3.227\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}, $L_{1}=$231.9\text{\,}\mathrm{p}\mathrm{H}, $I_{c,2}=$3.157\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}, and $L_{2}=$239.0\text{\,}\mathrm{p}\mathrm{H}.

To extract the remainder of the circuit parameters, we perform microwave spectroscopy on the single- and two-qubit systems Quintana17 . Applying a fixed $\Phi^{x}_{\rm cjj,2}$ to the second qubit, we sweep the barrier height of the first qubit, controlled by the external flux bias $\Phi^{x}_{\rm cjj,1}$. At every $\Phi^{x}_{\rm cjj,1}$, a long microwave pulse ($1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$) is applied to the first qubit to excite the two-qubit system. The energy eigenstates are read out by applying an adiabatic tilt to both qubits to transform the energy eigenstates into persistent current states, followed by a quench to freeze the dynamics of both qubits before readout. The pulse sequence is shown in Fig. 2(a). The excited state population of each qubit serves as a signal for detecting the energy spectrum of the coupled system. We collect effective single-qubit data by removing the potential barrier of the other qubit, making it monostable, and perform two-qubit spectroscopy for various $\Phi^{x}_{\rm cjj,2}$. In both sets of experiments we keep $\Phi_{q,1}^{x}=\Phi_{q,2}^{x}=0$, $M_{\text{12}}=0$ pH, while the capacitive coupling is always present. Jointly fitting the effective single- and two-qubit spectroscopy data to the coupled rf-SQUID model, Eq. (2), we obtain the rest of the circuit parameters, $C_{1}=119.5$ fF, $C_{2}=116.4$ fF, and $C_{12}=132$ fF.

Figures 2(b)-(d) show two-qubit spectroscopy data along with the numerical fit using the rf-SQUID model (2). One can see a clear $\Delta_{2}$ dependent anticrossing that is in good agreement with simulations (solid lines). Without any type of coupling, the spectral lines representing the first excited states of noninteracting qubits would cross as shown by the dashed lines 111Although the capacitive coupling is off, the capacitive loading on the individual qubits are kept in the simulation.. The anticrossing is therefore a signature of capacitive coupling (at $M_{12}=0$) and its energy gap is a measure of the coupling strength. The extracted anticrossing gaps of 0.77, 1.14, and 1.78 GHz at $\Delta_{2}=1.5$, 2.0, and 3.0 GHz, respectively, suggest a strong capacitive coupling.

II.2 Reduction to qubit model

We now reduce the continuous rf-SQUID model to a two-state (qubit) model, relevant for quantum computation. The flux degree of freedom is described by $\sigma^{z}_{i}$, with $\sigma_{1}^{z}\sigma_{2}^{z}$ representing the inductive coupling. The charge operator $Q_{i}{=}-i\hbar\frac{\partial}{\partial\Phi_{i}}$, on the other hand, is related to $\sigma^{y}_{i}$ since both are complex in the computational basis. Thus, the electrostatic coupling between the rf-SQUIDs gives rise to a $\sigma_{1}^{y}\sigma_{2}^{y}$ interaction. This term describes direct two-qubit cotunneling. In addition, a $\sigma_{1}^{x}\sigma_{2}^{x}$ term is obtained reflecting cotunneling mediated by the high energy states of the rf-SQUIDs. An effective two-qubit Hamiltonian can be represented as

[TABLE]

where $\alpha,\beta=\{x,y,z\}$. Since the continuous Hamiltonian (2) is real, the reduced Hamiltonian will also be real in the chosen basis, therefore, $J_{\alpha y}=J_{y\alpha}=0$, for $\alpha\neq y$. Parameters of Hamiltonian (3) are derived from the rf-SQUID model through a reduction techniques described in appendix A.2. Circuit parameters of the rf-SQUID model are extracted by fitting to experimental data obtained from uncoupled ($M_{12}{=}0$) rf-SQUIDs.

III Experimental results

In this section, we show that the reduced Hamiltonian (3) can explain experimental observations with no further fitting parameters.

III.1 Microwave spectroscopy

Figures 3(a)-(c) depict two-qubit spectroscopy at nonzero energy bias ($h_{i}\neq 0$). Solid white lines in these panels correspond to numerical simulations obtained using Hamiltonian (3) with no fitting parameters. The presence of a nonzero longitudinal field is necessary for nonstoquasticity, as there exists a unitary transformation that can remove positive off-diagonal elements at $h_{i}=0$ (see appendix A.3). The experimental parameters in Fig. 3(a) are the same as those in Fig. 2(b), except for the flux bias $\Phi_{q,i}^{x}$ determining $h_{i}$. This nonzero bias manifests itself in Fig. 3(a) as an avoided level crossing between the second and third excited states. Figure 3(b) shows the energy spectrum as a function of $\Phi_{q,1}^{x}$ while $\Phi_{q,2}^{x}=0.1$ m$\Phi_{0}$ and $M_{12}=0.55$ pH. The top two energy levels cross when $\Phi_{q,1}^{x}=-\Phi_{q,2}^{x}$. Energy spectrum as a function of $M_{12}$ is presented in Fig. 3(c). As in Fig. 3(a), the avoided level crossing observed at $M_{12}=0.55$ pH is a result of nonzero energy bias ($h_{i}\approx 0.15$ GHz at this point). Zero bias simulations are shown by the dashed lines. One can clearly see that the capacitive coupling introduces an asymmetry between the AFM and FM sides of the magnetic coupling. Without the coupling capacitor, the energy spectrum is expected to be symmetric about $M_{12}=0$. The theoretical (solid) lines in Fig. 3(c) are produced using the reduced Hamiltonian (3) with coupling parameters shown in Fig. 3(d). In Fig. 3(a)-(c), theory agrees well with experiment, with no additional fitting.

III.2 Coherent Oscillations

Finally, we measure quantum coherent oscillations in the coupled two-qubit system. Qubits are initialized in a computational basis state by applying a strong flux bias $\Phi_{q,i}^{x}$. The coherent oscillations are induced by quickly (within 200 ps) pulsing down the barriers of both qubits simultaneously. Before lowering the barriers, the flux bias $\Phi_{q,i}^{x}$ on each qubit is changed from its value at preparation pulse to its final value. Since the computational basis states are not the eigenstates of the total Hamiltonian, the system undergoes coherent oscillations between these states. After some dwell time $\tau$, the qubits are simultaneously quenched by rapidly raising their energy barriers via $\Phi^{x}_{\rm cjj,i}$, followed by qubit state readout. We repeat this process for a range of dwell times $\tau$ and coupling strengths $M_{12}$.

In Figs. 4(a),(b), we show the measured state population $P_{\uparrow\downarrow}$ and $P_{\downarrow\downarrow}$ for qubits initially prepared in $|{\downarrow\uparrow}\rangle$ and $|{\uparrow\uparrow}\rangle$ configurations, respectively, for $h_{i}=0$. The agreement between time domain experiments and numerical simulations, shown in Fig. 4(d),(e), justifies the two-qubit Hamiltonian (3) as a valid description of the circuit in Fig. 1(a). Deviation between theory and experiment can be attributed to decoherence, which is absent in simulations. Single qubit measurements reveal relaxation and dephasing times of $T_{1}=17$ ns and $T_{2}=16$ ns, respectively. These timescales are consistent with measurements of isolated flux qubits manufactured in the same stack but without the capacitive coupling.

Features of the energy spectrum as a function of $M_{12}$ (Fig. 3(c)) are reflected in coherent oscillations (Figs. 4(c),(f)). The initial configuration, $|{\uparrow\uparrow}\rangle$, has significant overlap with the second and third excited states in the AFM region ($M_{12}>0$). The slow oscillation frequency on the right half of figures is therefore related to the gap between these two states. The minimum gap at $M_{12}\approx$ 0.55 pH in Fig. 3(c) corresponds to the maximum slowdown at the same point in Figs. 4c,f. When $h_{1}=h_{2}=0$ this gap vanishes (see dashed lines in Fig. 3(c)), nullifying oscillations as seen in Fig. 4(b) and 4(e). This is a result of the destructive interference between the two-qubit cotunneling channels due to $\sigma_{1}^{x}\sigma_{2}^{x}$ and $\sigma_{1}^{y}\sigma_{2}^{y}$ terms and the indirect tunneling channels through sequential single qubit flips caused by $\sigma_{1}^{x}$ and $\sigma_{2}^{x}$ terms in Hamiltonian (3), as schematically illustrated in Fig. 1(b). Moreover, we see an additional signature of $J_{yy}$, since without such coupling Fig. 4(a) would be a reflection of Fig. 4(b) with respect to $M_{12}=0$ (see appendix D).

IV Stoquasticity

While existence of $J_{xx}$ or $J_{yy}$ coupling is a necessary condition for nonstoquasticity, it is not sufficient. One needs to show that the sign problem survives all local transformations. While this is intractable at large scales Lidar18 , it is feasible for two qubits Klassen18 . To demonstrate this, we extract coefficients in Hamiltonian (3) which describes Fig. 3(c). Figure 3(d) plots interaction parameters except for $J_{xz}$ and $J_{zx}$, which are negligibly small. Other parameters are provided in appendix A.2. We see that the electrostatic coupling between the rf-SQUIDs gives rise to a pronounced $J_{yy}$ that is almost constant over the whole range of $M_{12}$. The two-qubit cotunneling mediated by the higher energy states of the rf-SQUIDs leads to a $\sigma_{1}^{x}\sigma_{2}^{x}$ coupling with a coefficient $J_{xx}$. Both $J_{xx}$ and $J_{zz}$ depend on $M_{12}$ and therefore cannot be tuned independently. The magnitude of $J_{xx}$ is comparable to $J_{yy}$ for large antiferromagnetic coupling. Considering rotations in $x$-$z$ plane, both $J_{xx}$ and $J_{zz}$ terms can reduce the nonstoquastic contribution of $J_{yy}$. The Hamiltonian becomes stoquastic if either of them exceeds $J_{yy}$ in magnitude, as highlighted by the shaded area in the figure. Applying all possible local unitary transformations outlined in Ref. Klassen18 , we confirm nonstoquasticity in the unshaded region. Note that $\Delta$, $J_{xx}$ and $J_{yy}$ depend on the barrier heights, with $\Delta$ and $J_{yy}$ being proportional to each other (see appendix A.3). As a result, the width of the nonstoquastic region changes with $\Phi_{{\rm cjj},i}^{x}$, which controls the qubit tunneling amplitudes $\Delta$.

V Conclusion

In summary, we have fabricated two superconducting flux qubits coupled both inductively and capacitively. Starting from an rf-SQUID model, with experimentally extracted circuit parameters, we have obtained a reduced two-qubit Hamiltonian with $\sigma_{1}^{x}\sigma_{2}^{x}$, $\sigma_{1}^{y}\sigma_{2}^{y}$, and $\sigma_{1}^{z}\sigma_{2}^{z}$ interactions. We show that the reduced Hamiltonian can explain spectroscopy and coherent oscillation experiments. Considering all local transformations, we prove that the Hamiltonian is nonstoquastic in a wide range of parameters. Equilibrium statistics of such nonstoquastic quantum processors cannot be simulated by QMC due to the sign problem. Implementation of couplings via conjugate degreed of freedom such as charge and flux is an important step towards the development of universal quantum annealers Aharonov ; Biamonte ; Jordan . Our implementation is based on current, scalable superconducting fabrication technology that is ready to be expanded to a large number of qubits. Finally, we should mention that for large scale systems having $\sigma^{y}_{i}\sigma^{y}_{j}$ terms can make QMC simulation intractable even if the Hamiltonian is stoquastic. This is because finding a local transformation to cure sign problem at large scale is by itself intractable Lidar18 .

Acknowledgements

We are grateful to D. Lidar and M. Marvian for pointing out unitary transformations that cure sign problem at zero bias. We also thank J. Biamonte, D. DiVincenzo, I. Hen, J. Klassen, P. Love, J. Raymond, P. Saint-Jean, and B. Terhal for fruitful discussions, and F. Hanington and A. King for carefully reading the manuscript.

Appendix A System Hamiltonian

Our goal here is to derive a Hamiltonian of two rf-SQUIDs that have both, inductive and capacitive, couplings.

A.1 rf-SQUID with a cjj-loop

We begin with a description of a single rf-SQUID having a symmetric compound Josephson junction (CJJ) loop, as it is outlined in details in Ref. Harris10 (see Fig. 1 in the main text). The Hamiltonian of this system has three components:

[TABLE]

Here the Hamiltonians $H_{q}$ and $H_{\rm cjj}$ describe the main body of the qubit and the cjj-loop, respectively. They are expressed as

[TABLE]

The third term in Eq. (4) introduces an interaction between these two loops. The fluxes in the main body of the qubit and in the cjj-loop are denoted as $\Phi_{q}$ and $\Phi_{\rm cjj}$, their inductances are $L_{q}$ and $L_{\rm cjj}$, respectively. The external fluxes applied to the main body of the qubit and to the cjj-loop are $\Phi^{x}_{q}$ and $\Phi^{x}_{\rm cjj}.$ The cjj-loop has two symmetric branches $a$ and $b$. Each of them has one Josephson junction, with capacitance $C_{a}(C_{b})$ and with a persistent current $I_{a}(I_{b})$. The total current flowing through both branches of the cjj-loop is denoted as $I_{c}=I_{a}+I_{b}$, with the total Josephson energy $E_{J}=\frac{\Phi_{0}I_{c}}{2\pi},$ where $\Phi_{0}=\pi\hbar/e$ is the flux quantum, $e$ is the electron charge. Effective capacitances of the main qubit and the cjj-loop are $C_{q}=C_{a}+C_{b}$ and $C_{\rm cjj}=C_{a}C_{b}/(C_{a}+C_{b}),$ respectively. The charge, $Q_{\kappa}=-i\hbar\frac{\partial}{\partial\Phi_{\kappa}}$, and flux, $\Phi_{\kappa}$, are conjugate operators obeying the commutator: $[\Phi_{\kappa},Q_{\kappa}]_{-}=i\hbar.$

For qubits 1 and 2 studied in the main text of the paper we have $L_{\rm cjj,1}=17.0$ pH and $L_{\rm cjj,2}=17.2$ pH, while the body inductances are $L_{1}=231.9$ pH and $L_{2}=239.0$ pH. When $L_{\rm cjj}\ll L_{q},$ the dynamics of the fast degrees of freedom described by the operator $\Phi_{\rm cjj}$ is determined by the ground state of the Hamiltonian $H_{\rm cjj}$ (see Ref. Harris10 ). This ground state adiabatically follows the flux degrees of freedom, $\Phi_{q},$ in the main body of the rf-SQUID. One can then neglect the kinetic part of the Hamiltonian $H_{\rm cjj}$ and define an effective potential of the main loop by finding the minimum potential energy of the system (4) for each given flux $\Phi_{q}$:

[TABLE]

where

[TABLE]

This leads to the following effective Hamiltonian for the rf-SQUID,

[TABLE]

As a rough approximation, the effective potential energy of the rf-SQUID can be written as

[TABLE]

with a tunable Josephson energy

[TABLE]

However, in our simulations, we do not use the above approximation but actually do the minimization with respect to $\Phi_{\rm cjj}$ as given in Eq. (6) and described in detail in Ref. Harris10 .

We now analyze two rf-SQUIDs connected by a mutual inductance $M_{12}$ and by a capacitor $C_{12}$ as depicted in Fig. 1a of the main text. The inductive coupling between the main loops of the SQUIDs is given by the formula

[TABLE]

where $L_{i}\equiv L_{qi}$ ($i=1,2).$

The kinetic energy of two electrostatically-coupled SQUIDs, with capacitances $C_{i}\equiv C_{qi},$ has the form

[TABLE]

where we use the relation $V_{i}=\dot{\Phi}_{q,i},$ between a voltage $V_{i}$ on the $i$-junction and the flux $\Phi_{q,i}.$ Charge $Q_{i}=\partial{\cal L}/\partial\dot{\Phi}_{q,i}$ of the $i-$qubit is defined as a derivative of the two-qubit Lagrangian

[TABLE]

where $U_{{\rm eff},i}(\Phi_{q,i})$ is the effective potential energy of the $i-$qubit (see Eqs. (6) and (9)). Using Eqs. (12, 13), we obtain the relation between charge and time derivatives of flux,

[TABLE]

We can now write the total Hamiltonian of two coupled rf-SQUIDs, $H=\sum_{i}\dot{\Phi}_{q,i}Q_{i}-{\cal L}$, as

[TABLE]

Here

[TABLE]

is the Hamiltonian of $i$-qubit. The qubits are now characterized by the loaded capacitances:

[TABLE]

The inductive interaction between qubits is given by (11). The capacitive coupling is determined by the potential

[TABLE]

At large scales, the capacitative network renormalizes the capacitances of each rf-SQUID. In addition, capacitative loading should reduce sensitivity of the qubits to charge fluctuations, thus making this device more immune to charge noise.

A.2 Reduction approach

Diagonalizing the single-SQUID Hamiltonian $H_{i}$ (16) we obtain a set of energy eigenstates, $|{\chi_{\mu}^{i}}\rangle$, and eigenenergies, $\varepsilon_{\mu}^{i},$ such that

[TABLE]

For each SQUID we take into account a large number, $N_{i}\gg 1,$ of the energy eigenstates and write the Hamiltonian (15) of two coupled rf-SQUIDs in the basis formed by direct products $|{\chi_{\mu}^{1}\otimes\chi_{\nu}^{2}}\rangle\equiv|{\chi_{\mu}^{1}}\rangle\otimes|{\chi_{\nu}^{2}}\rangle$,

[TABLE]

Here charge and flux operators should be also written in the $|{\chi_{\mu}^{1}\otimes\chi_{\nu}^{2}}\rangle$ basis. The eigenstates, $|{\eta_{a}}\rangle,$ of the Hamiltonian (A.2) become

[TABLE]

where the amplitudes are given as

[TABLE]

The two-SQUID Hamiltonian,

[TABLE]

is characterized by the energy spectrum $\varepsilon_{a},$ with a total number of levels $N=N_{1}N_{2}.$

Working with continuous models becomes computationally challenging beyond a small number of coupled rf-SQUIDs. With the goal to use rf-SQUIDs as qubits, one needs to reduce the continuous Hamiltonian to a discrete (qubit) Hamiltonian. For uncoupled SQUIDs, we choose the following superpositions of two lowest-energy states with the mixing angle $\theta_{i}$:

[TABLE]

The basis states $|{\downarrow_{i}}\rangle$ and $|{\uparrow_{i}}\rangle$ correspond to the left and right circulating currents, or, equivalently, to the left and right sides of the SQUID potential well. The mixing angle $\theta_{i}$ is chosen to maximize the left-well population in the state $|{\downarrow_{i}}\rangle$ and, thus, the right-well population when the SQUID is in the state $|{\uparrow_{i}}\rangle.$ The interaction Hamiltonian mixes the states $|{\downarrow_{i}}\rangle,|{\uparrow_{i}}\rangle$ with higher energy states of the individual rf-SQUIDs that one needs to take into account for a correct description of the coupled system.

Since we are interested in the low energy spectrum of the coupled system, the number of eigenstates taken into account in Eqs. (19,A.2) can be truncated to just two states for each rf-SQUID. With this truncation, four eigenvectors, $|{\eta_{a}}\rangle$ ($a=1,\ldots,4$), of the two-SQUID system are approximated as

[TABLE]

where the amplitudes $c^{a}_{\mu\nu}$ are given in Eq. (22), and the normalization coefficient is calculated as

[TABLE]

We apply the Gram-Schmidt procedure to the four states given in Eq. (25) to obtain the orthonormalized set of the two-qubit basis states obeying the relation: $\langle\eta_{a}|\eta_{b}\rangle=\delta_{ab}.$ This reduction approach only works in the limit where $N_{a}\approx 1$ for $a=1,\ldots,4$. In this case the most of the population of two-SQUID system is distributed over the tensor products $|{\chi_{\mu}^{1}\otimes\chi_{\nu}^{2}}\rangle$ of two lowest eigenstates of the isolated rf-SQUIDs ($\mu,\nu=1,2$).

In order to derive the reduced Hamiltonian of two coupled SQUIDs we start with a Hamiltonian (23) truncated to the four lowest-energy states. In the energy basis this Hamiltonian is described by the diagonal matrix: $H={\rm diag}(\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},\varepsilon_{4}).$ The Hamiltonian $H$ can be written in the $\chi-$basis formed by the four vectors $|{\chi_{\mu}^{1}\otimes\chi_{\nu}^{2}}\rangle$ ($\mu,\nu=1,2)$ by applying the rotation matrix $R_{1}=(|{\eta_{1}}\rangle\;|{\eta_{2}}\rangle\;|{\eta_{3}}\rangle\;|{\eta_{4}}\rangle)$ such that $H_{\chi}=R_{1}HR_{1}^{T}$. Finally, the unitary matrix $R_{2}=(|{1}\rangle\,|{2}\rangle\,|{3}\rangle\,|{4}\rangle)$ rotates the Hamiltonian into the computational basis, formed by the vectors:

[TABLE]

After this rotation the Hamiltonian $H=R_{2}H_{\chi}R_{2}^{T}$ can be represented by the 4x4 matrix,

[TABLE]

or by the Hamiltonian of two coupled spins (qubits):

[TABLE]

where $\alpha,\beta={x,y,z}$. We use the following representation of the Pauli matrices:

[TABLE]

The parameters of the Hamiltonian (32) are extracted by comparing the operator $H=R_{2}H_{\chi}R_{2}^{T}$ with the matrix (31). Note that $J_{yx}=J_{xy}=J_{yz}=J_{yz}=0,$ since the Hamiltonian is real and the coefficients $J_{xz}$ and $J_{zx}$ are negligible: $J_{xz},J_{zx}\ll\Delta_{1},\Delta_{2}.$ The non-negligible parameters of the Hamiltonian $H$ describing interactions between the qubits are shown in Fig. 3(d) of the main text.

Figure 5(a) shows how the Hamiltonian parameters change during the annealing which is performed by changing $\Phi^{x}_{\rm cjj,i}$. All off-diagonal terms vanish as $\Phi^{x}_{\rm cjj,i}$ varies from right to left as the annealing proceeds. Linear dependence of $J_{xx}$ and $J_{yy}$ on $\Delta_{i}$ is evident from Fig. 5(b).

One can equivalently use the projection technique of Ref. Amin12S to derive a reduced Hamiltonian for two interacting subsystems, namely our SQUIDs. The reduced Hamiltonian obtained with the projection approach agrees with the two-qubit Hamiltonian given by Eqs. (31) and (32).

A.3 Unitary transformations and stoquasticity

Hamiltonian (31) can have positive off-diagonal matrix elements under several conditions but not all of those make the Hamiltonian nonstoquastic. One has to take into account all local transformations of the basis before deciding about nonstoquasticity of the Hamiltonian. Using the unitary transformation $H\to U^{\dagger}HU$, with $U=\otimes_{i}(\sigma^{z}_{i})^{(1-{\rm sign}\,\Delta_{i})/2}$, one can make all $\Delta_{i}\sigma^{x}_{i}$ terms in the Hamiltonian nonpositive regardless of the sign of $\Delta_{i}$. Hamiltonian (31) has positive matrix elements if $|J_{yy}|>|J_{xx}|$, since matrix elements $H_{14}$ and $H_{23}$ have opposite signs, due to $J_{yy}$, making one of them positive. However, this condition is not sufficient for nonstoquasticity. To demonstrate this fact, we apply a Hadamard rotation on both qubits turning $x-$axis to $z-$axis and vice versa, so that $\sigma^{x}_{i}=\tau^{z}_{i},\;\sigma^{y}_{i}=\tau^{y}_{i},$ and $\sigma^{z}_{i}=-\tau^{x}_{i}.$ Here $\tau^{x}_{i},\tau^{y}_{i},\tau^{z}_{i}$ are Pauli matrices of $i$-qubit in the rotated frame. After the Hadamard rotation the two-qubit Hamiltonian has the form:

[TABLE]

This Hamiltonian has positive matrix elements at $|J_{yy}|>|J_{zz}|.$ At non-zero biases, $h_{i}\neq 0,$ (and at non-zero $\Delta_{i}$) only rotations in the $x$-$z$ plane are allowed Klassen18 since rotations in $x$-$y$ and $z$-$y$ planes introduce complex matrix elements in the Hamiltonian. In Fig. 3(d) of the main text, we search over all possible local rotations outlined in detail in Ref. Klassen18 and find, in particular, that at the borders of the nonstoquastic range the $\sigma^{z}\sigma^{z}$ coupling is of order of the $\sigma^{y}\sigma^{y}$ interaction strength: $|J_{zz}|\sim|J_{yy}|.$ Figure 5(c) is an extension to this figure, displaying nonstoquastic region as a function of $M_{12}$ and $\Phi^{x}_{cjj,i}$, which varies during annealing. It follows from Fig. 5(a) and (c) that during the critical region of annealing, when off-diagonal elements of the Hamiltonian are present, the Hamiltonian can be nonstoquastic depending on the magnitude of the coupling.

The matrix (38) also provides the intuition necessary to see effects of $J_{xx}$ and $h_{1},h_{2}$ terms on experimentally measured spectra. When $J_{zz}=0$ and $h_{1}=h_{2}=0$, the Hamiltonian is block-diagonal as there is no interaction between aligned and anti-aligned states. In the absence of the $J_{xx}$ term, the eigenvalues of the Hamiltonian becomes $\pm\sqrt{J_{yy}^{2}+\left(\frac{\Delta_{1}\pm\Delta_{2}}{2}\right)^{2}}$. Without the $J_{xx}$ term, the two highest eigenvalues can only cross if $\Delta_{1}\Delta_{2}=0.$ However we see a clear level crossing in spectroscopy vs. $\Phi_{\text{cjj},1}^{x}$ plots in the middle row of Fig. 9 where $\Delta_{i}\neq 0$. This is a clear signature of a $J_{xx}$ type interaction. Furthermore, the second and third excited states, that cross in presence of $J_{xx}$ belong to two separate blocks of the Hamiltonian. When $h_{i}\neq 0$ the Hamiltonian is no longer block diagonal and the two highest excited states start interacting. Hence, the avoided crossings visible in energy vs $\Phi_{\text{cjj},1}^{x}$ and $M_{12}$ in Fig. 3(a),(c) of the main text are clear signatures of non-zero $h_{1,2}$.

Appendix B Qubit and coupler parameter characterization

B.1 Coupler characterization

We follow Ref. Harris09 for the characterization of the tunable magnetic coupler. In Fig. 6, we show the measured coupler $M_{12}$ versus the coupler flux bias $\Phi_{co}^{x}$. The data is fitted to a classical model:

[TABLE]

with $M_{co,q}$ the mutual inductance between the qubit and the coupler, $L_{co}$ the coupler inductance, $\beta=2\pi L_{co}I_{c,co}/\Phi_{0}$, $M_{12}^{(0)}$ the stray mutual inductance between the qubits, and $\varphi_{co}^{x}=2\pi\Phi_{co}^{x}$ the normalized external bias of the coupler. From the fitting, we obtained $M_{co,q}^{2}/L_{co}=10.77$ pH, $\beta=1.416$, and $M_{eff}^{(0)}=1.848$ pH. Using Eq. (39) and the above fitting parameters, we can set the coupler to any coupling strength within the range $|M_{12}|\leq 8.145$ pH.

B.2 Quasi-static qubit characterization

When the tunneling barrier is high, where the single-qubit tunneling is largely suppressed, the properties of a flux qubit can be described by a classical model Harris10 . In this regime, we measure the qubit persistent current versus the flux bias that control the barrier height $\Phi_{{\rm cjj}}^{x}$ across three flux quantum $\Phi_{0}$ for each qubit individually. Fits to the classical model (see Fig. 7) yield the following qubit parameters: $I_{c,1}=3.227~{}\mu$A, $I_{c,2}=3.157~{}\mu$A, $L_{1}=231.6~{}$pH, $L_{2}=239.0~{}$pH.

B.3 Single-qubit spectroscopy

The qubit energy eigenstates can be characterized by performing microwave spectroscopy. In Fig. 8(a), we show the pulse sequence for measuring the qubit spectroscopy and the effective qubit potential at each segment of the pulse sequence. The qubit is first initialized in its ground state by adiabatically preparing a symmetric double-well potential with a relatively low tunneling barrier. In this limit, the qubit frequency, which is set by the tunneling amplitude, is much larger than $1$GHz $(\Delta\gg k_{B}T)$. Then, we apply a microwave pulse $\Phi_{\text{q}}^{x}$ to excite the qubit to its excited state. We sweep the frequency of the microwave pulse from 0.5GHz to 8GHz to probe all excited states in that range. After the microwave pulse, an adiabatic tilt followed by a quench, that increases the barrier height, is applied to the qubit to project the qubit ground and excited states to the clockwise and counter-clockwise persistent current states for readout. The same pulse sequence was used Ref. Quintana17 . In this experiment, we use a rise time of 1 ns for both the adiabatic tilt and the quench on the barrier height. Throughout the single qubit spectroscopy experiments, we keep the magnetic coupling strength at $M_{12}=0$ and the other qubit at $\Phi_{\text{cjj}}^{x}=0.5\Phi_{0}$. In Fig. 8(b), we show the spectroscopy versus $\Phi_{\text{cjj,1}}^{x}$ for Qubit 1. The extracted qubit frequency from this figure along with the fit to rf-SQUID numerical model is shown in Fig. 8(c).

Effective single-qubit $\Delta$ corresponds to the energy gap between the ground state and the first excited state in the single-qubit spectroscopy measurements. Once we have an accurate model that predicts the location of the first excited state, the appropriate $\Phi_{\text{cjj},i}^{x}$ that needs to be applied to set the qubit to a given $\Delta$ is determined using the model extracted above. Note that the effective single-qubit $\Delta$ is influenced by the capacitive coupling that is always present. As a result the individual qubit dynamics are always affected by the presence of the secondary qubit even if it is in monostable state ($\Phi_{\text{cjj}}^{x}=0.5\Phi_{0}$).

B.4 Two-qubit spectroscopy

The coupled system energies can be probed by microwave spectroscopy in a way that is similar to single-qubit spectroscopy. In the main text we show the pulse sequence of the two-qubit spectroscopy. Although a microwave pulse is only applied to Qubit 1, the energy eigenstates of both qubits are read out by simultaneous adiabatic tilts followed by quench. We plot the spectroscopy, as the excited state population of Qubit 1, versus $\Phi_{\text{cjj},1}^{x}$ that controls the barrier height of Qubit 1. The barrier height of the second qubit is set to an effective single-qubit tunneling amplitude of $\Delta_{2}$. Using the previously extracted circuit parameters and initializing the fitting procedure by the design values for the unknown parameters, we extract $C_{1}=119.5$ fF, $C_{2}=116.4$ fF, and $C_{12}=132.0$ fF and $L_{\text{cjj},1(2)}=17.016(17.175)$pH. At this point, we have all the circuit parameters of this system as summarized in the following table:

[TABLE]

We also perform similar two-qubit spectroscopy for various $\Delta_{2}$ and $M_{12}$ at $h_{1}=h_{2}=0$ and compare the system energy spectra with numerical simulation using previously extracted parameters. In Fig. 9, we show the two-qubit spectroscopy data and calculated energy spectra of the system for all the combination of $\Delta_{2}/h=1.5$ GHz, 2.0 GHz, and 3.0 GHz and $M_{12}=0$ pH and $\pm 2$ pH. Very good agreement between the numerical model and the data is achieved. In all figures the solid (dashed) lines represent excitations from the ground (first excited) state.

Appendix C Pulse distortion compensation for coherent oscillations

Short duration pulse ($\leq 10$ ns) distortion imposes a great limitation on the fidelity of coherent qubit operation. In our experiment, the qubit control involves applying fast pulses to lower and rise the tunneling barrier of the flux qubits. These pulses have short rise and fall times ($\approx 200$ ps). Here, we discuss our method of measuring and correcting pulse distortion in-situ.

We first measure the single-qubit coherent oscillations at $\Delta/h=5$ GHz, where the qubit population in the computational basis is supposed to oscillate at the frequency $\Delta$. Pulse distortion is mainly caused by reflections. As a result, the $\Phi_{\text{cjj}}^{x}(t)$ signal reaching the qubit deviates from the ideal square pulse which may distort the frequency of the coherent oscillations from the target $\Delta$ at a given time $\tau$. To measure the distorted pulse in the time domain, we slice the coherent oscillation data into small time windows, typically with a length of 400 ps that contains more than one period of oscillation. We then extract the coherent oscillation frequency $\Omega(\tau)$ in each slice starting at time $\tau$ and treat this frequency as a measurement of the instantaneous $\Delta_{m}(\tau)=\Omega(\tau)$. The pulse distortion at any given time is then calculated as $\tau$ by $\delta\Phi_{\text{cjj}}^{x}(\tau)=(\Delta_{m}(\tau)-\Delta)/\frac{\partial\Delta}{\partial\Phi_{\text{cjj}}^{x}}$ where $\frac{\partial\Delta}{\partial\Phi_{\text{cjj}}^{x}}$ can be evaluated numerically using the qubit model. To correct for the pulse distortion, we simply apply the first order correction to the applied pulse with $\Phi_{\text{cjj,corr}}^{x}(\tau)=\Phi_{\text{cjj}}^{x}(\tau)+\delta\Phi_{\text{cjj}}^{x}(\tau)$ at the same time $\tau$. As the above correction at time $\tau$ may lead to distortions at times $t>\tau$, we iterate the entire measurement and correction procedure until the corrected pulse converges. In practice, this pulse distortion correction procedure converges within 5 iterations. In Fig. 10, we show the coherent oscillation data and the extracted instantaneous tunneling amplitude $\Delta_{m}$ with uncorrected pulse (blue) and corrected(green) pulses $\Phi_{\text{cjj}}^{x}$, shown respectively. The measured tunneling amplitude $\Delta_{m}/h$ is fixed to the target value 5 GHz after the pulse distortion compensation is applied, which significantly increases the fidelity of the coherent oscillation.

Appendix D Coherent oscillation protocol

The pulse sequence and the effective two-qubit potential at each pulse segment of a coherent oscillation protocol is shown in Fig. 11. The qubits are initialized in two steps. A large flux-bias, $\Phi_{\text{q,prep}}^{x}$, is applied to tilt the potential in its monostable state. Next, the tunneling strength is reduced, keeping the potential tilted. Once both qubits are prepared in a computational state controlled by the signs of $\Phi_{\text{q,prep}}$ pulses, the tilts are removed. The coherent oscillations are induced by reducing the barriers of both qubits simultaneously to allow quantum fluctuations drive the coherent dynamics of the system. As the computational basis states are not the eigenstates of the total Hamiltonian, the system undergoes coherent oscillations between the computational basis states with near-degenerate Ising energies. After some time $\tau$, the states are read out by measuring qubit-persistent currents after simultaneously quenching both qubits. The rise and fall times of the pulse, about 200 ps, are much faster than the dynamics of the qubits to snapshot the qubit population at the end of evolution. We repeat this process for a range of dwell times $\tau$ and magnetic coupling strengths. This protocol works equally well for a single-qubit if the qubit potential is kept at its monostable state.

In coherent oscillation experiments, the system is prepared and measured in one of the four computational basis states $|\!\uparrow\uparrow\rangle$, $|\!\uparrow\downarrow\rangle$, $|\!\downarrow\uparrow\rangle$, and $|\!\downarrow\downarrow\rangle$. In Fig. 11, we show that the coherent oscillations are indeed mostly between the two low-energy configurations at strong magnetic couplings. At strong FM coupling with $M_{12}=-2$ pH, the coherent oscillations occur between $|\!\uparrow\uparrow\rangle$ and $|\!\downarrow\downarrow\rangle$, whereas the coherent oscillations occur between $|\!\uparrow\downarrow\rangle$ and $|\!\downarrow\uparrow\rangle$ at strong AFM coupling with $M_{12}=2$ pH.

Appendix E Coherence

The phase coherence time $T_{2}$ between the energy eigenstates of each qubit can be characterized by measuring single-qubit coherent oscillations. As in Eq. (A.2), at zero bias we define the computational basis as $|{\uparrow}\rangle{=}\frac{|{\chi_{1}}\rangle{+}|{\chi_{2}}\rangle}{\sqrt{2}}$ and $|{\downarrow}\rangle{=}\frac{|{\chi_{1}}\rangle{-}|{\chi_{2}}\rangle}{\sqrt{2}}$ for each qubit. In Fig. 12, we show the coherent oscillations of Qubit 1 at $\Delta_{1}=2$ GHz. A fit to the decay envelop of the oscillations yields $T_{2}=16.2$ ns. Energy relaxation time $T_{1}$ at the same operation point is characterized by measuring the qubit excited state population versus the delay time between initialization and readout. By fitting the excited-state population decay to an exponential function, we extract $T_{1}=17.5$ ns. We perform the same characterization on Qubit 2 and obtain $T_{2}=14.3$ ns and $T_{1}=17.4$ ns at $\Delta_{2}=2$ GHz. We have also fabricated uncoupled qubits with the same parameters and observed the same coherence parameters, confirming that the capacitive coupling does not degrade qubit coherence.