A transmon quantum annealer: decomposing many-body Ising constraints into pair interactions

In PAQC, a system with a larger number of spins but only local interactions is used to mimic the spin-glass system with all-to-all Ising interactions at the end of the adiabatic sweep. This is possible because there exists an unambiguous mapping of eigenstates of the spin-glass system to the lowest-energy states of the PAQC system, which are separated by an energy gap from the remaining states of the enlarged Hilbert space of the PAQC system.

The mapping between states of the spin-glass system with all-to-all connectivity to the states of the PAQC system is accomplished as follows. The Ising interaction ${\sigma }_{z,i}{\sigma }_{z,j}$ has two two-fold degenerate eigenvalues: +1 and −1 corresponding to parallel ($| \uparrow \uparrow \rangle$ and $| \downarrow \downarrow \rangle$) and anti-parallel ($| \uparrow \downarrow \rangle$ and $| \downarrow \uparrow \rangle$) alignment of the spins. In the PAQC system, the state of this Ising interaction is represented by a single spin, say state $| \uparrow \rangle$ corresponds to a parallel configuration and $| \downarrow \rangle$ to a anti-parallel configuration of spins. This already provides us with an unambiguous mapping from states of the spin-glass with all-to-all connectivity with N spins to the states of the PAQC system with $K=N(N-1)/2$ spins. However, the larger Hilbert space of the PAQC system already prohibits an unambiguous mapping of all states to corresponding states of the spin-glass Hamilton operator. The effective reduction of the Hilbert space of the PAQC system is accomplished with constraints. Constraints are terms in the Hamilton operator that inflict energy penalties, larger than any other energy scale in the system, on configurations of the PAQC systems that do not have a counterpart in the spin-glass system. The PAQC system is a spatial arrangement of the spins that enables only local constraint, i.e. energy penalties involving only neighbouring spins (c.f. figure 2(b)) and [25]). It consists of a cubic lattice of spins with overall pyramid geometry where the constraints act locally on every unit cell of the cubic lattice. The Hamilton operator of PAQC has the form

$$\begin{eqnarray}&&{H}_{{\rm{PAQC}}}=A(t)\displaystyle \sum _{i}^{K}{\sigma }_{x}^{(i)}+B(t)\displaystyle \sum _{i}^{K}{J}_{i}{\sigma }_{z}^{(i)}+C(t)\displaystyle \sum _{l}^{L}{\sigma }_{z}^{(l,n)}{\sigma }_{z}^{(l,s)}{\sigma }_{z}^{(l,e)}{\sigma }_{z}^{(l,w)}.\end{eqnarray} \tag{ 2 }$$

Here, ${\sigma }_{\{x,y,z\}}^{(i)}$ are the Pauli operator spins in the PAQC system. The terms $A(t)$, $B(t)$ and $C(t)$, are three independent schedule functions that define the annealing protocol. $A(t)$, the schedule function of the initial Hamiltonian, is tuned from its maximal value ${A}_{{\rm{\max }}}$ to 0 during the adiabatic sweep. $B(t)$, the schedule function of the Hamiltonian that encodes the optimization problem, is tuned from 0 to its maximal value ${B}_{{\rm{\max }}}$ during the adiabatic sweep. $C(t)$ is the schedule function of the constraint term. $C(t)$ can be switched from 0 to ${C}_{{\rm{\max }}}$, with ${C}_{{\rm{\max }}}$ dominating every other energy scale in the system. However, PAQC also opens the possibility to keep the constraints 'always on'. The constraints are implemented with four-body Ising interactions involving the north $(l,n)$, east $(l,e)$, south $(l,s)$ and west $(l,w)$ spins of every plaquette l. They are constraints because at least at the end of the adiabatic sweep they are the dominating energy contribution and all dynamics is restricted to reside in their ggroundstate manifold. The scheme features error correction capabilities for uncorrelated noise models [32]. However, in [33], it is pointed that this does not correct for correlated errors as expected from a error model based on quantum Monte Carlo simulations, where minor embedding schemes show larger robustness.

In the following, we show how to decompose a general many-body Ising constraints to two-body Ising interactions with the help of ancilla spins.

The constraint decomposition follows a two-step process: (i) We present a general recursive technique to decompose many-body Ising constraints to three-body Ising constraints with ancilla spins and use this technique to decompose the four-body Ising constraint of PAQC to two three-body Ising constraints, and (ii) we show a decomposition of a three-body Ising constraint with two-body Ising interactions with one ancilla spin. This constraint decomposition technique follows a similar logic to the embedding described above. However, instead of embedding the whole eigensystem of an operator in another, we only need to embed the ggroundstate manifold of the constraint and make sure that the ground state manifold of the embedding constraint is still separated by a sufficient gap from the rest of the states.

Decomposing many- to three-body constraints. The many-body Ising interaction $C{\prod }_{i=1}^{M}{\sigma }_{z}^{(i)}$, has two multiply degenerate eigenvalues 1 and −1 corresponding to product states of an even or odd number of spins in the spin down state, respectively. In the following, we call them even- or odd-parity states, respectively. Depending on the sign of the constraint strength C, either the states with even or odd parity are the ground state manifold, i.e. are 'allowed' under the constraint. Given a set of qubits in eigenstates of their respective ${\sigma }_{z}$ operators, we could detect the parity of the state by either detecting the parity of the whole set or subdividing the set into two subsets, A (spins 1 through to k) and B (spins $k+1$ through to N), and detecting the parity of the two subsets. If the two subsets have the same parity, the parity of the set is even and vice versa. We can translate this principle into the domain of Ising interactions with the help of an ancilla spin $(a)$,

$$\begin{eqnarray}&&C\displaystyle \prod _{i=1}^{N}{\sigma }_{z}^{(i)}\equiv \pm | C| \left({\sigma }_{z}^{(a)}\displaystyle \prod _{i=1}^{k}{\sigma }_{z}^{(i)}-{\rm{sgn}}(C){\sigma }_{z}^{(a)}\displaystyle \prod _{i=k+1}^{N}{\sigma }_{z}^{(i)}\right).\end{eqnarray} \tag{ 3 }$$

The equivalence '$\equiv$' is here defined as the existence of a one-to-one mapping between the ground state manifolds in the sense that the states are identical for the non-ancilla spins. Additionally, all remaining states of the embedding constraint have to be separated by a gap from the ground state manifold. The overall sign of the embedding constraint $\pm | C|$ above is actually arbitrary because the states of the ground state manifolds for the system with ancilla are in both cases identical for the non-ancilla spins. The constraint decomposition technique is valid for odd $(C\gt 0)$ and even $(C\lt 0)$ parity constraints. The four-body Ising constraint of PAQC is a even party constraint, which is why we have chosen to further illustrate even parity decomposition. The reasoning for odd parity constraints parallels the reasoning given in the following.

The ground state manifold of the embedding constraint for an even parity constraint $(C\lt 0)$ is characterised by a odd parity for the subset A including the ancilla spin $(a)$ and the subset B including the same ancilla spin $(a)$. If the ancilla spin $(a)$ is in the spin-up state, then both subsets A and B have to have odd parity in order for the state to be in the ground state manifold. If the ancilla spin $(a)$ is in the spin down state, both subsets have to have even parity. The ancilla spin $(a)$ therefore 'communicates' the parity information between subsets. Therefore, the set of spins not including the ancilla spin has the same ground state manifold as for the many-body Ising interaction. The gap of the embedding constraint above, c.f. Equation (3), is actually of the same size as the gap of the constraints involving the subsets A and B by virtue of the non-perturbative nature of the decomposition. This opens the door to very strong effective Ising constraints.

To summarise, with the ancilla spin $(a)$, one can decompose a constraint of order N, i.e. involving N spins, into two constraints of order $k+1$ and $N-k+1$. This decomposition can be iterated recursively with both subsystems, thus one can decompose any constraint of arbitrary order N into a set of coupled three-body constraints. The recursive sequence of decomposition steps generates a tree structure which is a highly desirable feature for two-dimensional setups like the ones used in circuit quantum electrodynamics because qubits are connected without any crossings and no air-bridges are needed, not including control and readout circuitry.

Figure 3 depicts an example of this recursive algorithm for the decomposition of a five-body constraint. We start by splitting the five-body term into a right branch with subset $\{1,2\}$ and a left branch with subset $\{3,4,5\}$. In each branch, we add the ancilla spin $({a}_{1})$ resulting in two terms, a three-body constraint in the right branch and a four-body constraint in the left branch. The left branch is finished, as three-body constraints cannot be further decomposed with this scheme. We continue with the left branch and split the four-body constraint into two subsystems containing $\{4,5\}$ and $\{{a}_{1},3\}$, respectively. Adding the ancilla qubit $({a}_{2})$ on both sides results in a three-body constraint in all leafs of the tree which terminates the procedure. After joining the leafs, the final decomposed five-body constraint is depicted in figure 3(b).

Decomposing three- to two-body constraints. The second decomposition step (ii) aims at replacing the remaining three-body constraints by pair interactions including another ancilla spin. The above recursive decomposition cannot be applied to three-body terms, as one subbranch would again result in a three-body constraint. However, it is possible to construct a set of interactions that feature the same degenerate ground state as the three-body Ising constraint if we ignore the state of the ancilla spin for the moment [30, 31, 34]. There are actually many possible combinations of pair interactions that share the same degenerate ground state. We chose the particular solution, where (i) all interactions have the same sign and (ii) there are no crossings (depicted in figure 1(b))

$$\begin{eqnarray}{\sigma }_{z}^{(1)}{\sigma }_{z}^{(2)}{\sigma }_{z}^{(3)} & \equiv & {\sigma }_{z}^{(1)}{\sigma }_{z}^{(2)}+{\sigma }_{z}^{(2)}{\sigma }_{z}^{(3)}+{\sigma }_{z}^{(3)}{\sigma }_{z}^{(1)}\\ & & +\displaystyle \sum _{i=1}^{3}[2{\sigma }_{z}^{(i)}{\sigma }_{z}^{(a)}-{\sigma }_{z}^{(i)}]-2{\sigma }_{z}^{(a)}.\end{eqnarray} \tag{ 4 }$$

Here, the symbol '$\equiv$' quotes the definition of equivalence for constraints given above. Using equation (4) together with equation (6), the parity quantum annealing scheme [25] can be implemented with pair interactions only. The full layout of interactions and ancillas is shown in figure 1(c).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let us now apply the above algorithm to decompose all constraints in the PAQC scheme. The first step (i) is a trivial tree in this case. We simply split the four-body constraint into two three-body constraints,

$$\begin{eqnarray}&&-C{\sigma }_{z}^{(n)}{\sigma }_{z}^{(e)}{\sigma }_{z}^{(s)}{\sigma }_{z}^{(w)}\to -C({\sigma }_{z}^{(n)}{\sigma }_{z}^{(e)}{\sigma }_{z}^{(a)}+{\sigma }_{z}^{(a)}{\sigma }_{z}^{(s)}{\sigma }_{z}^{(w)}).\end{eqnarray} \tag{ 5 }$$

Plugging equation (5) into equation (2) results in the adiabatic protocol

$$\begin{eqnarray}H(t) & = & A(t)\displaystyle \sum _{i}^{K}{\sigma }_{x}^{(i)}+B(t)\displaystyle \sum _{i}^{K+M}{J}_{i}{\sigma }_{z}^{(i)}\\ & & +C(t)\displaystyle \sum _{l}^{L}({\sigma }_{z}^{(n)}{\sigma }_{z}^{(w)}{\sigma }_{z}^{(a)}+{\sigma }_{z}^{(a)}{\sigma }_{z}^{(e)}{\sigma }_{z}^{(s)}).\end{eqnarray} \tag{ 6 }$$

Here, we add one ancilla spin for each plaquette which makes a total of $M={(N-2)}^{2}/2$ spins as shown in figure 1. Each of the ancilla spins is shared by two three-body Ising interactions. Compared with the implementation with four-body Ising constraints, the number of constraints doubled in the above realisation with three-body Ising constraints. In the second step (ii), each three-body term in equation (6) is replaced by the right hand side of equation (4), which results in the final layout with pair interaction only depicted in figure 1(c). Note that this scheme may be optimised by combining pairs of spins (figure 1(c) (yellow)), which results in a layout with less spins but introduces crossings of interactions.

The main source of errors in our annealing scheme are Landau-Zener transitions at the minimal gap, as we are considering a coherent annealing sweep. Therefore, it is vital that the constraint decomposition technique presented above does not considerably decrease the minimal gap during the annealing sweep. In the following, we illustrate the differences in the instantaneous eigenenergies of one plaquette for the original PAQC scheme with the four-body Ising constraint (c.f. figure 1(a)) and the above presented decomposition with three ancillas and pair interactions (c.f. figure 1(b)). We start by diagonalizing the instantaneous Hamiltonian of a plaquette with the decomposed four-body Ising constraint. We consider a protocol that starts with a strictly local transverse field Hamiltonian, A = 1, B = 0 and C = 0 in equation (2). Then the local longitudinal fields and constraint terms are adiabatically switched on ($B\to 1$, $C\to 1$ ), while the the ${\sigma }_{x}$ term is switched adiabatically off ($A\to 0$).

In the absence of local programmable fields J_i, the result is a superposition of all eight constraint-satisfying configurations. The time-dependent spectrum of this ideal sweep is shown in figure 4(a). In the presence of local programable fields, the final ground state is the solution of the optimization problem (c.f. figure 4(b)). The local field terms J_i lift the degeneracy of the final state. All unphysical, or not constraint-satisfying states, are separated by a energy gap of size C at the end of the sweep, thereby validating our above decomposition of the four-body Ising constraint.

Zoom In Zoom Out Reset image size

Instead of adiabatically switching on the constraints, there is also the possibility of always-on constraints, $C(t)=\mathrm{const}.=1$, which circumvents difficulties in the implementation. Either there is an efficient protocol to generate the ground state of the initial Hamiltonian with A = 1 and C = 1, or one has dominating local transversal fields $A\gg C$ and initialises the system in the product state of all local transverse field Hamiltonians. The instantaneous eigenenergies of the two protocols, the ramping up of constraints and always-on constraints (A = C = 1 at the beginning of the sweep), are illustrated in figure 5 using the same parameters as in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To quantify the efficiency of the protocols, we compare the minimal gaps from the idealised model with four-body Ising interaction and the ancilla-based approach. Figure 6(a) shows a comparison of the minimal gap of a plaquette using N = 500 random instances with $| J| \lt 0.5$. The minimal gap in the proposed three-ancilla setup $({{\rm{\Delta }}}_{\min }^{a})$ is compared with the minimal gap (${{\rm{\Delta }}}_{\min }^{4{\rm{body}}}$) of the model with idealised four-body interaction. The gap in the idealised model is in general larger than in the ancilla-based implementation, but remarkably, for a large number of instances, they are almost identical. We further compared the gap in the ramp protocol (${{\rm{\Delta }}}_{\min }^{a}$) as in figure 6(a) with the gap $({{\rm{\Delta }}}_{\min }^{c})$ in the always-on protocol starting from a superposition state (c.f. figure 6(b)). The gap in the always-on protocol is smaller for most instances, with only a small systematic difference between the two in favour of the ramp protocol.

The Transmon was developed on the basis of charge qubits as a improvement on their coherence times. However, while the charge qubit provides a natural mapping to a spin Hamilton operator with longitudinal and transverse fields, the Transmon is more appropriately described as a harmonic oscillator with small non-linearity. In the following, we show how the longer coherence times of the Transmon are a direct consequence of the lost transverse field term and how to combine the transverse field term and the long coherence times in the following subsection.

A charge qubit is an arrangement of two superconducting islands that are connected by a Josephson junction. The state of the Cooper pair condensate on the two islands is completely described by a canonical conjugate pair of operators: the number of Cooper pairs N that have tunnelled through the Josephson junction starting from the electrical neutral state and the phase difference of the Cooper pair condensate between the two islands ϕ. The Hamilton operator of the charge qubit is given by the sum of the electrical energy stored in the capacitor (C) that is formed by the two islands and the inductive energy of the Josephson junction,

$$\begin{eqnarray}&&{H}_{{\rm{cq}}}=4{E}_{C}{(N-{n}_{g})}^{2}-{E}_{J}\cos (\phi ),\end{eqnarray} \tag{ 7 }$$

where ${E}_{C}={e}^{2}/2C$ is the charging energy with the elementary charge e, E_J the Josephson energy and ${n}_{g}={C}_{g}{V}_{g}/(2e)$ the Cooper pair equivalent of an external potential V_g with capacitance C_g. Here, n_g can be a externally applied classical field or the quantum electrical field of uncontrolled sources as well as quantum fields from coupled resonators or other charge qubits. If we ignore the Josephson energy for the moment, the states of the charge qubit as a function of the externally applied classical electric field n_g are parabolas with origin at integer n_g and degeneracies for half integer n_g. Each parabola is the energy of a state $| n\rangle$ with exactly n Cooper pairs. The degeneracy at half integer n_g is lifted by the introduction of the Josephson energy. This is possible because of the tunnelling of Cooper pairs between the islands through the Josephson junction. Around this avoided crossing in the spectrum, one can truncate the Hilbert space of the charge qubit to the two states of exactly defined Cooper pairs and thereby accomplish a mapping to an effective spin Hamilton operator with a longitudinal field given by the externally applied classical electrical field ${h}_{z}=2{E}_{C}(1-2{n}_{g})$ and transverse field given by the Josephson energy ${h}_{x}={E}_{J}/2$. The qubit states corresponding to spin 'up' and 'down' states are states of well-defined Cooper pairs $| n\rangle \to | \uparrow \rangle$ and $| n+1\rangle \to | \downarrow \rangle$ away from the avoided crossing. Exactly in the middle of the avoided crossing, the states are symmetric and anti-symmetric superpositions of states of well-defined charge, $(| n\rangle +| n+1)/\sqrt{2}\to | \uparrow \rangle$ and $(| n\rangle -| n+1)/\sqrt{2}\to | \downarrow \rangle$. Because the charge is not well-defined at the avoided crossing, the charge qubit shows increased resilience to external fluctuations in the electric field. The Transmon improves on this idea by generating a universal avoided crossing for a very strong Josephson tunneling which mixes all charge states $| n\rangle$, ${E}_{J}/{E}_{C}\gg 1$. In this parameter regime, the difference between the maximal and minimal energy as a function of n_g decreases exponentially as a function of ${E}_{J}/{E}_{C}$, while the non-linearity, i.e. the difference between the energy ground state to the first excited state ${E}_{1}-{E}_{0}$ and the first to the second excited ${E}_{2}-{E}_{1}$ state decreases polynomially. Here, E_n are the eigenenergies of the Transmon. The Transmon thereby combines the advantages of a universal sweet spot with sufficient non-linearity, necessary for individual addressing of its eigenstates with microwave drives.

The eigenstates of the Transmon are characterised by small zero-point fluctuations in the phase by virtue of ${E}_{J}\gg {E}_{C}$, which is why it is appropriate to truncate the Josephson energy that is given by the cosine of the phase to fourth order,

$$\begin{eqnarray}&&{H}_{{\rm{cq}}}=4{E}_{C}{(N-{n}_{g})}^{2}-{E}_{J}\cos (\phi )\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\approx \,4{E}_{C}{(N-{n}_{g})}^{2}-{E}_{J}+\displaystyle \frac{{\phi }^{2}}{2}-{E}_{J}\displaystyle \frac{{\phi }^{4}}{24}.\end{eqnarray} \tag{ 9 }$$

The first two terms of the truncated Hamilton operator can be identified by the Hamilton operator of a harmonic oscillator. The externally applied n_g represents a shift in space for the harmonic oscillator that does not change the eigenenergies which shows the existence of the universal 'sweet spot' [2]. We may therefore as well set ${n}_{g}\to 0$ in the following and introduce lowering and raising operators to describe the number of Cooper pairs and phase of the Transmon,

$$\begin{eqnarray}&&N=\displaystyle \frac{i}{2}{\left(\displaystyle \frac{{E}_{J}}{2{E}_{C}}\right)}^{\displaystyle \frac{1}{4}}(a-{a}^{\dagger })\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\phi ={\left(\displaystyle \frac{2{E}_{C}}{{E}_{J}}\right)}^{\displaystyle \frac{1}{4}}(a+{a}^{\dagger }).\end{eqnarray} \tag{ 11 }$$

The truncated Hamilton operator of the Transmon can be approximated with the help of a rotating wave approximation that transforms the non-linear fourth-order phase term to a Kerr non-linearity,

$$\begin{eqnarray}&&4{E}_{C}{(N-{n}_{g})}^{2}-{E}_{J}+\displaystyle \frac{{\phi }^{2}}{2}-{E}_{J}\displaystyle \frac{{\phi }^{4}}{24}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\approx (\sqrt{8{E}_{J}{E}_{C}}-{E}_{C}){a}^{\dagger }a-\displaystyle \frac{{E}_{C}}{2}{a}^{\dagger }{a}^{\dagger }{aa}={H}_{1T}.\end{eqnarray} \tag{ 13 }$$

Here, we shifted the ground state energy to zero. Here, $(\sqrt{8{E}_{J}{E}_{C}}-{E}_{C})$ is the Transmon energy including a renormalisation of the Transmon energy coming from the normal ordering procedure of the fourth-order phase term. The Transmon is therefore more appropriately described as a non-linear harmonic oscillator and as a consequence of the universal 'sweet spot' there is no mapping to a spin Hamilton operator with a transverse field. In the following subsection, we show how to reintroduce a transverse field in a rotating frame by a constant microwave drive.

A viable way to reintroduce a transverse field is to excite the Transmon with a constant microwave drive of frequency ${\omega }_{d}$, strength A, and Hamilton operator ${H}_{{\rm{drive}}}=A({{ae}}^{i{\omega }_{d}t}+{a}^{\dagger }{e}^{-i{\omega }_{d}t})$. In a frame rotating with this microwave drive, $U=\exp (-i{\omega }_{d}{{ta}}^{\dagger }a)$, the transverse field term is reintroduced,

$$\begin{eqnarray}&&{U}^{\dagger }({H}_{{\rm{1T}}}+{H}_{{\rm{drive}}})U-i{\rm{\hslash }}{U}^{\dagger }\dot{U}=2{{Ja}}^{\dagger }a-\displaystyle \frac{{E}_{C}}{2}{a}^{\dagger }{a}^{\dagger }{aa}+A(a+{a}^{\dagger }),\end{eqnarray} \tag{ 14 }$$

where $2J=\sqrt{8{E}_{C}{E}_{J}}-{\rm{\hslash }}{\omega }_{d}$ is the energy equivalent of the frequency difference between Transmon and microwave drive. Here, ${\rm{\hslash }}$ is the reduced Planck constant. The symbol J is chosen in anticipation for the use in the PAQC scheme as the on-site longitudinal field that encodes the coupling strength of the associated spin-glass annealer. The quantum annealing processor is operated in a regime where we can neglect occupation of states higher than the first excited state, which is ensured by choosing a driving strength smaller than the non-linearity $A\lt {E}_{C}$. Therefore, we may project the Hamilton operator to the qubit subspace to get the effective Hamilton operator in the qubit subspace,

$$\begin{eqnarray}&&{H}_{{\rm{1T,Qbit}}}=J{\sigma }_{z}+A{\sigma }_{x}.\end{eqnarray} \tag{ 15 }$$

To summarize, in a frame rotating with the microwave drive, the effective longitudinal field is defined by the energy equivalent of the frequency difference between the Transmon and the microwave drive J, which allows longitudinal fields smaller in magnitude than the Ising pair interaction presented below. Additionally, a transverse field given by the strength of the microwave drive A, is reintroduced.

Transmon qubits are naturally coupled via their electric field [1]. First, we shortly illustrate why this capacitive coupling is linear in the Transmon field operators while the desired Ising pair interaction is of fourth order. Afterward, we present how to implement the Ising pair interaction with the non-linearity provided by Josephson junctions.

For capacitively coupled Transmons, n_g, as given above in equation (8), of one Transmon $(a)$ is a function of the electric field of another Transmon $(b)$ and their mutual capacitance C_c. The interaction term in the Hamilton operator is

$$\begin{eqnarray}&&{H}_{{\rm{int,cap}}}=\displaystyle \frac{{(2e({N}_{a}-{N}_{b}))}^{2}}{2{C}_{c}}\approx {g}_{{\rm{cap}}}({a}^{\dagger }b+{{ab}}^{\dagger }).\end{eqnarray} \tag{ 16 }$$

where we assumed the validity of the Transmon approximation and performed a rotating wave approximation valid for small coupling strengths $g\lt {\omega }_{a}+{\omega }_{b}$. Here, ${\omega }_{a}$ and ${\omega }_{b}$ are the frequencies of the two Transmons, while a and b are their field operators. The natural capacitive coupling of two Transmons therefore provides us with an exchange interaction that is linear in the field operators of the two Transmons. However, the fundamental building block of the four-body Ising constraint that is needed for PAQC is, by virtue of the decomposition techniques described above (c.f. figure 1(c)), the Ising pair interaction ${\sigma }_{z}^{(a)}{\sigma }_{z}^{(b)}$. In terms of field operators of the Transmon, a and ${a}^{\dagger }$, every ${\sigma }_{z}$ operator is already a quadratic operator because $2({a}^{\dagger }a-1/2)$ corresponds to ${\sigma }_{z}$ in the qubit subspace and consequentially our desired Ising pair interaction term is of fourth order. Therefore, this interaction has to be implemented with the only non-linear element at our disposal for superconducting qubits, the Josephson junction. One possibility is to connect an island of the first Transmon to one of the islands of the other, which results in a coupling energy given by

$$\begin{eqnarray}&&{H}_{{\rm{int,ind}}}={E}_{{J}_{c}}\cos ({\phi }_{a}-{\phi }_{b})\approx {g}_{{\rm{ind}}}({a}^{\dagger }b+{{ab}}^{\dagger })+{g}_{{\rm{nonlin,ind}}}{(a+{a}^{\dagger }-b-{b}^{\dagger })}^{4}.\end{eqnarray} \tag{ 17 }$$

Here, we again assumed the validity of the Transmon approximation and performed a rotating wave approximation to get the linear exchange interaction. Additionally, we truncated the cosine of the Josephson junction at fourth order which is valid given the Transmon approximation. In addition to our desired fourth-order non-linear coupling, we get a linear exchange interaction. It is possible to cancel the linear terms provided by a Josephson junction with a parallel capacitive interaction [37, 38] because ${g}_{{\rm{cap}}}$ and ${g}_{{\rm{ind}}}$ have a different sign. However, this scheme limits the interaction energy to a small fraction of the Transmon non-linearity because the coupling capacitance has to be considerably smaller than the Transmon's capacitance to inhibit unwanted long-range interactions [38].

We present another possibility to couple Transmons with Josephson junctions in the form of a JRM [36] where the linear coupling vanishes for symmetry reasons. The interaction energies exceed the coupling energies attainable by the capacitively shunted Josephson junction scheme by at least an order of magnitude without introducing unwanted long-range interactions.

Four superconducting islands joined to a ring by identical Josephson junctions define a JRM (c.f. appendix A). Its energy is proportional to the product of the cosines of three orthogonal modes if the flux threaded through the JRM loop vanishes. The JRM's modes comprise two differential modes involving opposing superconducting islands ${\varphi }_{x}={\varphi }_{1}-{\varphi }_{3}$, ${\varphi }_{y}={\varphi }_{2}-{\varphi }_{4}$ and a third mode involving all islands ${\varphi }_{z}={\varphi }_{1}-{\varphi }_{2}+{\varphi }_{3}-{\varphi }_{4}$. Here, ${\varphi }_{i}={\int }_{-\infty }^{t}{V}_{i}{dt}$ is the flux variable defined as the time integral of the electrical potential V_i of island i. It is connected to the phase variable introduced above by ${\varphi }_{i}=\phi {\varphi }_{0}$ with the reduced magnetic flux quantum ${\varphi }_{0}={\rm{\hslash }}/(2e)$.

We associate the two differential modes ${\varphi }_{x}$ and ${\varphi }_{y}$ of the JRM with the modes that register the two qubits ${\varphi }_{a}$ and ${\varphi }_{b}$ by connecting them with conducting leads. For the sake of simplicity, we assume all JRM junctions to be equal with Josephson energy ${E}_{{\rm{jrm}}}$, although our setup tolerates the usual fabrication inaccuracies (c.f. appendix B). The qubit capacitances C_J are also assumed to be equal, although unequal qubit capacitances do not alter the main result. The qubit Josephson junctions can be implemented as direct current superconducting quantum interference devices (dc-SQUID) and are modelled as tunable Josephson energies E_Ja and E_Jb. Note that the Transmons might as well be implemented with fixed frequencies, provided that we can individually change the drive frequencies. They are the 'nobs' of the quantum annealing processor that need to be adjusted to encode the problem we want to solve. Contrary to the typical use case of the JRM with ${\varphi }_{{ext}}\ne 0$, we are interested in the case ${\varphi }_{{ext}}=0$. The two-Transmon Lagrangian with the corresponding JRM contribution reads (for a full derivation of the Lagrangian of the JRM c.f. 4),

$$\begin{eqnarray}&&{L}_{2T}^{\prime }=\displaystyle \frac{{C}_{J}}{2}{\dot{\varphi }}_{a}^{2}+{E}_{{Ja}}\cos \left(\displaystyle \frac{{\varphi }_{a}}{{\varphi }_{0}}\right)+\displaystyle \frac{{C}_{J}}{2}{\dot{\varphi }}_{b}^{2}+{E}_{{Jb}}\cos \left(\displaystyle \frac{{\varphi }_{b}}{{\varphi }_{0}}\right),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&+\,4{E}_{{\rm{jrm}}}\cos \left(\displaystyle \frac{{\varphi }_{a}}{2{\varphi }_{0}}\right)\cos \left(\displaystyle \frac{{\varphi }_{b}}{2{\varphi }_{0}}\right)\cos \left(\displaystyle \frac{{\varphi }_{z}}{2{\varphi }_{0}}\right).\end{eqnarray} \tag{ 19 }$$

Notice that the third mode ${\varphi }_{z}$ of the JRM does not possess any capacitive term. In a mechanical picture, this mode is a massless particle moving in a one-dimensional potential that depends on the position of other massive particles. Therefore, it will immediately adjust itself to the potential minimum. Even a small capacitance for the ${\varphi }_{z}$-mode which will be present in every setup does not change this qualitative picture. In analogy to the elimination of the coupling degree of freedom for the g-mon [39], we find

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\displaystyle \frac{\partial {L}_{2T}^{\prime }}{\partial {\dot{\varphi }}_{z}}=\displaystyle \frac{\partial {L}_{2T}^{\prime }}{\partial {\varphi }_{z}}=0,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\cos \left(\displaystyle \frac{{\varphi }_{a}}{2{\varphi }_{0}}\right)\cos \left(\displaystyle \frac{{\varphi }_{b}}{2{\varphi }_{0}}\right)\sin \left(\displaystyle \frac{{\varphi }_{z}}{2{\varphi }_{0}}\right)=0.\end{eqnarray} \tag{ 21 }$$

${\varphi }_{a}$ and ${\varphi }_{b}$ are the flux variables of Transmons, i.e. their zero-point fluctuations are small and centred around 0. Therefore, ${\varphi }_{z}=0$ is the only possible solution to the constraint given by the Euler–Lagrange equation for mode ${\varphi }_{z}$. We Legendre transform the resulting Lagrangian and quantise the theory to get the Hamilton operator,

$$\begin{eqnarray}&&{H}_{2T}^{\prime }=4{E}_{C}{N}_{a}^{2}-{E}_{{Ja}}\cos ({\phi }_{a})+4{E}_{C}{N}_{b}^{2}-{E}_{{Jb}}\cos ({\phi }_{b}),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&-\,4{E}_{{\rm{jrm}}}\cos \left(\displaystyle \frac{{\phi }_{a}}{2}\right)\cos \left(\displaystyle \frac{{\phi }_{b}}{2}\right),\end{eqnarray} \tag{ 23 }$$

where ${N}_{a/b}=C{\dot{\varphi }}_{a/b}/(2e)$ and ${\phi }_{a/b}={\varphi }_{a/b}/{\varphi }_{0}$. In preparation for the Transmon approximation, we introduce bosonic lowering and raising operators with the following dependency on the Cooper pair and phase operators,

$$\begin{eqnarray}&&{N}_{x}=\displaystyle \frac{i}{2}{\left(\displaystyle \frac{{E}_{{Jx}}+{E}_{{\rm{jrm}}}}{2{E}_{C}}\right)}^{\tfrac{1}{4}}(x-{x}^{\dagger }),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\phi }_{x}={\left(\displaystyle \frac{2{E}_{C}}{{E}_{{Jx}}+{E}_{{\rm{jrm}}}}\right)}^{\tfrac{1}{4}}(x+{x}^{\dagger }),\end{eqnarray} \tag{ 25 }$$

where $x\in \{a,b\}$. Truncating the Hamilton operator to quartic order results in the canonical Transmon approximation of a harmonic oscillator mode with a quartic non-linearity given by the charging energy of the Transmon,

$$\begin{eqnarray}&&{H}_{2T}^{\prime }\approx {H}_{2T}={E}_{a}{a}^{\dagger }a-\displaystyle \frac{{E}_{C}}{2}{a}^{\dagger }{a}^{\dagger }{aa}+{E}_{b}{b}^{\dagger }b-\displaystyle \frac{{E}_{C}}{2}{b}^{\dagger }{b}^{\dagger }{bb}-g{(a+{a}^{\dagger })}^{2}{(b+{b}^{\dagger })}^{2}\end{eqnarray} \tag{ 26 }$$

with

$$\begin{eqnarray}&&{E}_{x}=\sqrt{8{E}_{C}({E}_{{Jx}}+{E}_{{\rm{jrm}}})},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&g=\displaystyle \frac{{E}_{C}}{2}\displaystyle \frac{{E}_{{\rm{jrm}}}}{\sqrt{{E}_{{Ja}}+{E}_{{\rm{jrm}}}}\sqrt{{E}_{{Jb}}+{E}_{{\rm{jrm}}}}}.\end{eqnarray} \tag{ 28 }$$

As a last step, we add the individual microwave drives introduced above for both Transmons and transform into a frame rotating with the microwave drives with ${U}_{x}=\exp (-i{\omega }_{d,x}{{tx}}^{\dagger }x)$ for $x\in \{a,b\}$. Note, that the ${\sigma }_{z}$ coupling strength in equation (28) is of the order of the non-linearity of the Transmon. The validity of our projection to the individual qubit subspaces still only requires the driving strengths A_a and A_b to be small compared with the on-site non-linearity E_C. The effective Hamilton operator in the rotating frame projected to the qubit subspace is, therefore,

$$\begin{eqnarray}&&{H}_{{\rm{2T,Qbit}}}={A}_{a}{\sigma }_{x}^{(a)}+{A}_{b}{\sigma }_{x}^{(b)}+{J}_{a}{\sigma }_{z}^{(a)}+{J}_{b}{\sigma }_{z}^{(b)}-g{\sigma }_{z}^{(a)}{\sigma }_{z}^{(b)},\end{eqnarray} \tag{ 29 }$$

where ${J}_{x}={E}_{x}-\delta -{\rm{\hslash }}{\omega }_{d,x}$ for $x\in \{a,b\}$. With δ a renormalization of the Transmon qubit energy that originates in the normal ordering process of the qubit non-linearity as well as the non-linear Ising coupling.

With the longitudinal and transverse field terms and the Ising pair interaction, we have all the ingredients to build a Transmon quantum annealer with the PAQC scheme. Errors on the local field terms ${\sigma }_{z}$ and ${\sigma }_{x}$ are highly reduced in the rotating frame due to ns-precision of state-of-the-art microwave drives used in recent experiments [40]. As the optimization problem is encoded in the z-direction of the Hamiltonian, any error in the ${\sigma }_{z}$ terms must be smaller than the range of programmability $C\gg | J| \gg {\delta }_{z},{k}_{B}T$. Errors arise from the second-order processes neglected in the derivation of equation (29) that result in additional ${\sigma }_{x}{\sigma }_{x}$ couplings (c.f. appendix B). Let us also note additional opportunities for implementations with our scheme. The presented approach allows for an implementation of a Transmon quantum annealer with fixed frequency qubits [41] that show increased resilience with respect to environmental noise. By driving each Transmon individually ${\omega }_{d,a}\ne {\omega }_{d,b}$ and upon changing into the rotating frame, one gets an effective time-independent Hamilton operator. One might even conceive of a Transmon quantum annealer where each Transmon is replaced by a capacitor and the JRMs alone provide the inductive energy for the on-site Transmons, as the effective Josephson energy of both Transmons is the sum of the Josephson energy of the Transmon junction and the JRM junction equation (27).