Quantum annealing with all-to-all connected nonlinear oscillators

Quantum annealing aims at solving combinatorial optimization problems mapped to Ising interactions between quantum spins. Here, with the objective of developing a noise-resilient annealer, we propose a paradigm for quantum annealing with a scalable network of two-photon-driven Kerr-nonlinear resonators. Each resonator encodes an Ising spin in a robust degenerate subspace formed by two coherent states of opposite phases. A fully connected optimization problem is mapped to local fields driving the resonators, which are connected with only local four-body interactions. We describe an adiabatic annealing protocol in this system and analyse its performance in the presence of photon loss. Numerical simulations indicate substantial resilience to this noise channel, leading to a high success probability for quantum annealing. Finally, we propose a realistic circuit QED implementation of this promising platform for implementing a large-scale quantum Ising machine.

M any hard combinatorial optimization problems arising in diverse areas such as physics, chemistry, biology and social science [1][2][3][4] can be mapped onto finding the ground state of an Ising Hamiltonian. This problem, referred to as the Ising problem, is in general NP-hard 5 . Quantum annealing, based on adiabatic quantum computing (AQC) 6,7 , aims to find solutions to the Ising problem, with the hope of a significant speedup over classical algorithms. In AQC, a system is evolved slowly from the non-degenerate ground state of a trivial initial Hamiltonian to that of a final Hamiltonian encoding a computational problem. During the time-evolution, the energy spectrum of the system changes and, for the adiabatic condition to be satisfied, the evolution must be slow compared to the inverse minimum energy gap between the instantaneous ground state and the excited states. The scaling behaviour of the gap with problem size, thus, determines the efficiency of the adiabatic annealing algorithm.
In order to perform quantum annealing, the Ising spins are mapped to two levels of a quantum system, that is, a qubit, and the optimization problem is encoded in the interactions between these qubits. Adiabatic optimization with a variety of physical realizations such as nuclear magnetic spins 8 and superconducting qubits 9,10 has been demonstrated. However, despite great efforts, whether these systems are able to solve large problems in the presence of noise remains an open question 11 . As a consequence, it is imperative to search for implementations with improved resilience to noise.
A general Ising problem is defined on a fully connected graph of Ising spins. However, efficient embedding of large problems with such long-range interactions is a challenge because physical systems more naturally realize local connectivity. In one approach, a fully connected graph of Ising spins is embedded in a so-called Chimera graph 12,13 . Alternatively, a more recent embedding scheme was proposed by Lechner, Hauke and Zoller (LHZ) 14 in which N logical Ising spins are encoded in M ¼ N(N À 1)/2 physical spins with M À N þ 1 constraints. Each physical spin represents the relative configuration of a pair of logical spins. An all-to-all connected Ising problem in the logical spins is realized by mapping the logical couplings onto local fields acting on the physical spins and a problemindependent four-body coupling to enforce the constraints. This simple design requires only precise control of local fields, making it attractive for scaling to large problem sizes.
In the following, we present a physical platform for quantum annealing that is both scalable and shows robustness to noise; here we propose to encode the Ising problem in a network of two-photon-driven Kerr-nonlinear resonators (KNR). In our scheme, a single Ising spin is mapped to two coherent states with opposite phases, which constitute a two-fold degenerate eigenspace of the two-photon-driven KNR in the rotating frame of the drive 15 . Here we propose to realize quantum adiabatic algorithms by encoding a quantum spin in quasi-orthogonal coherent states. The dominant source of error in this system is single-photon loss from the resonators. However, since coherent states are invariant under the action of the photon jump operator, the encoded Ising spin is stabilized against bit flips. We describe a circuit QED implementation of a quantum annealing platform, where a fully connected graph of Ising spins is embedded using the LHZ scheme, relying on effective local magnetic fields and four-body coupling between KNRs. The adiabatic optimization is carried out by initializing the resonators to vacuum, and varying only single-site drives to adiabatically evolve the system to the ground state of the embedded Ising problem. This realization allows us to encode arbitrary Ising problems with no restriction on connectivity, or on the signs and amplitudes of the spin-spin couplings.
Encoding Ising spins in the phase of coherent states has previously been explored in the context of classical Ising machines [16][17][18][19][20][21] . The quantum case was considered in ref. 22. However, this previous study focussed on idealized quantum systems without noise analysis and did not consider practical implementations of these ideas. In contrast, our analysis considers the performance of quantum annealing in the presence of single-photon loss, by far the dominant loss mechanism. Crucially, we numerically demonstrate that the probability for the system to jump from the instantaneous ground state to one of its excited states due to photon loss during the adiabatic protocol is greatly suppressed as compared to conventional qubit implementations with equal noise strengths. This resilience to the detrimental effects of photon loss leads to high success probabilities in finding the optimal solution to optimization problems mapped on two-photon-driven KNRs. This noise resilience, in combination with simple initialization and final state detection by homodyne measurement of the resonators' field amplitudes, opens the door to realizing a large-scale quantum annealer with favourable noise resistance.

Results
Adiabatic protocol for quantum annealing. Quantum annealing is executed with the time-dependent Hamiltonian whereĤ i is the initial trivial Hamiltonian whose ground state is known, andĤ p is the final Hamiltonian at t ¼ t which encodes an Ising spin problem: h j is the Pauli-z matrix for the ith spin and J i,j is the interaction strength between the ith and jth spin. Crucially, the initial and final Hamiltonian do not commute. For simplicity, we have assumed a linear time dependence, but more complex annealing schedules can be used. The system, initialized to the ground state ofĤ i , adiabatically evolves to the ground state of the problem Hamiltonian,Ĥ p , at time t¼t ) 1=D min , where D min is the minimum energy gap 6 .
Single spin in a two-photon-driven Kerr-nonlinear resonator. The Hamiltonian of a two-photon driven KNR in a frame rotating at the drive frequency is given byĤ 0 ¼ À Kâ y2â2 þ E pâ y2 þâ 2 À Á , where K is the Kerr-nonlinearity and E p the strength of the two-photon drive. In a KNR, the coherent states AE a 0 j i, which are eigenstates of the photon annihilation operator, are stabilized by the two-photon drive with a 0 ¼ ffiffiffiffiffiffiffiffiffiffiffi 15). This statement can be visualized more intuitively by considering the metapotential obtained by replacing the operatorsâ andâ y with the complex classical variables x þ iy and x À iy in the expression forĤ 0 (ref. 23). As shown in Fig. 1a, this metapotential is an inverted double well with two peaks of equal height at (±a 0 , 0), corresponding to two stable points (see Supplementary Note 1). This is consistent with the quantum picture according to which the coherent states AE a 0 j iare two degenerate eigenstates of H 0 with eigenenergy E 2 p =K (ref. 15) (see Methods). Taking advantage of this well-defined two-state subspace, we choose to encode an Ising spin 0 j i; 1 j i f gin the stable states À a 0 j i; a 0 j i f g . Importantly, this mapping is robust against single-photon loss from the resonator when the rate of single-photon loss is small k ( 8E p , a condition than can readily be realized in superconducting circuits 15 . Moreover, the photon jump operator a leaves the coherent states invariantâ 0= 1 j i¼ AE a 0 0= 1 j i. As a result, if the amplitude a 0 is large such that 0= 1 h jâ 1= 0 j i¼ Ç a 0 e À 2 a0 j j 2 $ 0, a single-photon loss does not lead to a spin-flip error.
Having defined the spin subspace, we now discuss the realization of a problem Hamiltonian in this system. As an illustrative example, we first address the trivial problem of finding the ground state of a single spin in a magnetic field. Consider the Hamiltonian of a two-photon-driven KNR with an additional weak single-photon drive of amplitude E 0 , As illustrated in Fig. 1b for small E 0 , under this drive, the two peaks located at (±a 0 , 0) in the metapotential associated toĤ p are now of unequal height with the peak at ( À a 0 , 0) lower than the one at (a 0 , 0) if E 0 40, and vice versa for E 0 o0. These two states remain stable, but have different energies, indicating that the singlephoton drive induces an effective magnetic field on the Ising spins 0 j i; 1 j i f g. Indeed, a full quantum analysis shows that if E 0 j j ( 4K a 0 j j 3 , then AE a 0 j i remain the eigenstates ofĤ p but their degeneracy is lifted by 4E 0 a 0 (ref. 15). In other words, in the spin subspace,Ĥ p can be expressed asĤ p ¼2E 0 a 0 s z þ const:, witĥ s z ¼ 1 j i 1 h j À 0 j i 0 h j. This is the required problem Hamiltonian for a single spin in a magnetic field. This simple observation will play an essential role in the implementation of the LHZ scheme discussed below. Note that for larger E 0 , the eigenstates can deviate from coherent states (see Supplementary Note 2). Choosing E 0 j j ( 4K a 0 j j 3 , however, ensures that 0 j i; 1 j i f g are indeed coherent states to an excellent approximation, such that 0= 1 h jâ 1= 0 j i$0 and the encoded subspace remains well protected from the photon loss channel. Following equation (1), we require an initial Hamiltonian which does not commute with the final problem Hamiltonian and which has a simple non-degenerate ground state. This is achieved by introducing a finite detuning d 0 40 between the drives and resonator frequency. In a frame rotating at the frequency of the drives, the initial Hamiltonian is chosen asĤ i ¼ d 0â yâ À Kâ y2â2 with d 0 oK. This choice of initial Hamiltonian generates large phase fluctuations that helps maximize quantum tunnelling to states with well-defined phase at the final stages of the adiabatic evolution. In this frame, the ground and first excited states are the vacuum 0 j i and singlephoton Fock state n¼1 j i, respectively, and which are separated by an energy gap d 0 . If a photon is lost from the resonator, the excited state n¼1 j i decays to the ground state 0 j i which, on the other hand, is invariant under photon loss. Since it is simple to prepare in the superconducting circuit implementation that we consider below, the vacuum state is a natural choice for the initial state.
The time-dependent Hamiltonian required for the adiabatic computation can be realized by slowly varying the two-and singlephoton drive strengths and detuning so thatĤ 1 ðtÞ¼ 1 À t=t ð ÞĤ i þ t=t ð ÞĤ p , realizing equation (1) for a single-spin. Note that the form ofĤ 1 ðtÞ conveniently ensures that the nonlinear Kerr term is time-independent. The time-dependent detuning is achieved by tuning the single-and two-photon drive frequencies (see Methods). By adiabatically controlling the frequency and amplitude of the drives it is possible to evolve the state of the KNR from the vacuum 0 j i at t ¼ 0, to the ground state of a single Ising spin in a magnetic field at t ¼ t. Figure 2a shows the change of the energy landscape throughout this evolution found by numerically diagonalizing the instantaneous HamiltonianĤ 1 ðtÞ for E p ¼4K, a 0 ¼ 2, E 0 ¼0:2K and d 0 ¼ 0.2K. The minimum energy gap D min is indicated. As illustrated by the Wigner functions in Fig. 2b, a resonator initialized to the vacuum state at t ¼ 0 evolves through highly non-classical and non-Gaussian states towards the ground state 0 j i at t ¼ t, with t $ 30=D min in this example. If, on the other hand, the KNR is initialized to the single-photon Fock state at t ¼ 0, then it evolves to the first excited state 1 j i at t ¼ t. The average probability to reach the correct ground state is 99.9% for both E 0 40 and E 0 o0. The 0.1% probability of erroneously ending in the excited state arises from non-adiabatic errors and can be reduced by increasing the evolution time. For example, for t ¼ 60/D min we find a success probability of 99.99%.
Effect of single-photon loss. An appealing feature of this implementation is that, at the start of the adiabatic protocol, the ground (vacuum) state is invariant under single-photon loss. Similarly, at the end of the adiabatic protocol at t ¼ t, irrespective of the problem Hamiltonian (that is, E 0 40 or E 0 o0) the ground state (coherent states 0 j i or 1 j i) is also invariant under singlephoton loss. It follows that towards the beginning and end of the protocol, photon loss will not induce errors. Moreover, we find that, even at intermediate times 0otot, the ground state of H 1 ðtÞ remains largely unaffected by photon loss. This can be understood intuitively from the distortion of the metapotential, as shown in Fig. 2c at t ¼ 0.2t for the same parameters as Fig. 2a. While the metapotential still shows two peaks, the region around the lower peak (corresponding to the ground state) is a circle whereas that around the higher peak (corresponding to the excited state) is deformed. This suggests that the ground state is closer to a coherent state and, therefore, more robust to photon loss than the excited state (see also Supplementary Note 3). Quantitatively, the effect of single-photon loss is seen by numerically evaluating 24,25 the transition matrix elements hc g ðtÞjâ c e ðtÞ j i, c e ðtÞ h jâjc g ðtÞi for the duration of the protocol, where jc g ðtÞi and c e ðtÞ j i are the ground and excited state of H 1 ðtÞ respectively. As shown in Fig. 2d, the transition from the ground to excited state is greatly suppressed throughout the whole adiabatic evolution. This asymmetry in the transition rates distinguishes AQC with two-photon-driven KNRs from implementations with qubits 26 , something that will be made even clearer below with examples.
Two coupled spins with driven KNRs. Before going to larger lattices, consider the problem of two interacting spins embedded in a system of two linearly coupled KNRs, each driven by a two-photon drive and given by the Hamiltonian Here, J 1,2 is the amplitude of the single-photon exchange coupling and, for simplicity, the two resonators are assumed to have identical parameters. For small J 1,2 , this Hamiltonian can be expressed in Figure 1 | Contour plot of the metapotential. Metapotential corresponding where K is the Kerr-nonlinearity, E p and E 0 are the strengths of the two-photon and singlephoton drive respectively, with E p ¼4K and (a) and, following equation (1), the full timedependent Hamiltonian for the two-spin problem iŝ Although it is possible to tune these parameters in time, with the above form ofĤ 2 ðtÞ, both the linear coupling and the Kerr-nonlinearity are fixed throughout the adiabatic evolution.
The ground state ofĤ 2 ð0Þ is the vacuum state if the initial detuning is greater than the single-photon exchange rate, d 0 4J 1,2 . On the other hand, at t ¼ t, the two degenerate ground states for anti-ferromagnetic (ferromagnetic) coupling are 0; 1 Accordingly, numerical simulations with both resonators initialized to vacuum show the coupled system to reach the entangled state N 0; 1 j iþ 1; 0 j i ð Þ and N 0; 0 j iþ 1; 1 j i ð Þ , under anti-ferromagnetic and ferromagnetic coupling, respectively. In these expressions, N ¼1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi , the state fidelity is 99.9%. Moreover, the probability that the system is in any one of the states 0= 1; 0= 1 j iis 99.99%, showing that the evolution is to a very good approximation restricted to this computational subspace.
While the states used in this encoding are tolerant to photon loss, coherence between a superposition of those states is reduced. However the success probability (see Methods) in solving the Ising problem remains high as it depends only on the diagonal elements of the density matrix (for example, 0; 1 h jrðtÞ 0; 1 j i). As an illustration, with the large loss rate k ¼ 50/t, while the fidelity of the final state to the desired superposition N 0; 0 j iþ 1; 1 j i ð Þor N 0; 1 j iþ 1; 0 j i ð Þ decreases to 37.6%, the average success probability of the algorithm is 75.2%.
To characterize the effect of noise, a useful figure of merit is the ratio D min /k of the minimum energy gap to the loss rate. The dependence of the average success probability on this ratio is presented in Fig. 3 for the algorithm implemented using KNRs (green squares) with single-photon loss k or qubits (red squares) with pure dephasing g f . In practice, the average success probability is computed by varying the loss rates at fixed D min and t ¼ 20/D min , and is averaged over all instances of the problem of two coupled spins (that is, ferromagnetic and anti-ferromagnetic). In the presence of pure dephasing, the success probability with qubits saturates to 50% for large g f . This is a consequence of the fact that the steady state of the qubits is an equal weight classical mixture of all possible computational states. On the other hand, for KNRs with a finite k, the rate at which the instantaneous ground state jumps to the excited state (/ c e ðtÞ h jâjc g ðtÞi) is small compared to the rate at which the instantaneous excited state jumps to the ground state (/ hc g ðtÞjâ c e ðtÞ j i). As a result, even with large single-photon loss rate, for example D min =k $ 1, the success probability is B75%. Consequently, in the presence of equivalent strength noise, a two-photon-driven KNR implementation of AQC has superior performance compared to a qubit implementation (see Methods for details).
All-to-all connected Ising problem with the LHZ scheme. The above scheme can be scaled up with pairwise linear couplings in a network of KNRs, while still requiring only single-site drives. However, unlike the above one-and two-spin examples, most optimization problems of interest require controllable long-range interactions between a large number of Ising spins. Realizing such highly non-local Hamiltonian is a challenging hardware problem, but it may be solved by embeddings such as the LHZ scheme 14 that map the Ising problem on a graph with local interactions only. In this approach, the relative configuration of pairs of N logical spins is mapped to M ¼ N(N À 1)/2 physical spins. A pair of logical spins, in which both spins are aligned |0, 0i or |1, 1i (or anti-aligned |0, 1i or |1, 0i), is mapped on the two levels of the physical spin. The coupling between the logical pairs J i,j (i ¼ 1, ..N) is encoded in local magnetic fields on the physical We now describe a circuit QED platform implementing the LHZ scheme by embedding the physical spins in the eigenbasis 0 j i; 1 j i f gof two-photon-driven KNRs. In practice, KNRs can be realized as a superconducting microwave resonator terminated by a flux-pumped SQUID. The non-linear inductance of the SQUIDs induces a Kerr-nonlinearity, and a two-photon drive is introduced by flux-pumping at twice the resonator frequency. This is the exact same setup that is used to realize Josephson parametric amplifiers (JPAs), and we will therefore refer to this implementation of a KNR as a JPA in the following 15,[27][28][29] . We envision a quantum annealing platform to be built with groups of four JPAs of frequencies o r,i (i ¼ 1, 2, 3, 4) interacting via a single Josephson junction (JJ) as illustrated in Fig. 4a. In the figure, the four different frequencies of the resonators are indicated by four different colours. To realize a time-dependent two-photon drive, the SQUID loop of each JPA is driven by a flux pump with tunable amplitude and frequency. The pump frequency is varied close to twice the resonator frequency, o p;k ðtÞ ' 2o r;i (see Methods). Additional single-photon drives, whose amplitude and frequency can be varied in time, are also applied to each of the JPAs and play the role of effective local magnetic fields. Local four-body couplings are realized through the nonlinear inductance of the central JJ (see Supplementary Note 5). Choosing o p,k (t) þ o p,l (t) ¼ o p,m (t) þ o p,n (t) and taking the resonators to be detuned from each other, the central JJ induces a coupling of the form À Cðâ y kâ y lâ mân þ h:c:Þ in the instantaneous rotating frame of the two-photon drives. This four-body interaction is an always-on coupling and its strength C is determined by the JJ nonlinearity. Such a group of four JPAs, which we will refer to as a plaquette, is the central building block of our architecture and can be scaled in the form of the triangular lattice required to implement the LHZ scheme 14 . Note that while JPAs within a plaquette have different frequencies, only four distinct JPA frequencies are required for the entire lattice as illustrated in Fig. 4c. Lastly, the LHZ scheme also requires additional N À 2 physical spins at the boundary that are fixed to the up state and which are implemented in our scheme as JPAs stabilized in the eigenstate 1 j i by two-photon drives. As an illustration, Fig. 4b depicts all the possible interactions in an Ising problem with N ¼ 5 logical spins and Fig. 4c shows the corresponding triangular network of coupled JPAs. To implement the adiabatic protocol for a general N-spin Ising problem with the triangular network of M JPAs, the timedependent Hamiltonian in a frame where each of the JPAs rotate at the instantaneous drive frequency can be written aŝ wherê As mentioned above, M ¼ N(N À 1)/2 while J k is the singlephoton drive which induces the local effective magnetic field on the kth resonator and C is the local four-body coupling between  Loss-rate dependence of the success probability for the two-spin adiabatic algorithm in a system of two-photon-driven KNRs with single-photon loss k (green squares) and qubits with pure dephasing at rate g f (red squares Here, this is realized by standard homodyne detection which can resolve the two coherent states AE a 0 j iallowing the determination of the ground state configuration of the spins. In order to demonstrate the adiabatic algorithm for a non-trivial case, we embed on a plaquette a simple three-spin frustrated Ising problem, in which the spins are anti-ferromagnetically coupled to each other,Ĥ p ¼J P k;j¼1;2;3ŝ z;kŝz;j with J40. This Hamiltonian has six degenerate ground states in the logical spin basis. Following the LHZ approach, a mapping of N ¼ 3 logical spins requires M ¼ 3 physical spins (in our case three JPAs) and one physical spin fixed to up state (in our case a JPA initialized to the stable eigenstate 1 j i). Since the physical spins 0 j i; 1 j i f g encoded in the JPAs constitute the relative alignment of the logical spins, there are three possible solutions in this basis: 1; 0; 0 j i, 0; 1; 0 j i and 0; 0; 1 j i. The time-dependent Hamiltonian for the adiabatic protocol then follows from equation (2) with M ¼ 3. The anti-ferromagnetic coupling between the logical spins is represented by the single-photon drives on each JPA with amplitude J k ¼ J40. At t ¼ 0, the ground state of this Hamiltonian is the vacuum 0; 0; 0 j i. For appropriate magnitude of the four-body coupling (see Supplementary Note 4), the problem Hamiltonian can be expressed asĤ p ¼Jð4a 0 E 2 p =KÞ P 3 k¼1 s z;i À 2C a 0 j j 4 s z;1 s z;2 s z;3 s z;4 þ const: . This realizes the required problem Hamiltonian in the LHZ schemeĤ LHZ;N¼3 p .
To illustrate the performance of this protocol, we numerically simulate the evolution subjected to the Hamiltonian of equation (2) with the three resonators initialized to vacuum and the fourth initialized to the state 1 and k ¼ 0, we find that the success probability to reach the ground state to be 99.3%. The reduction in fidelity arises from the non-adiabatic errors. The probability for the system to be in one of the states 0= 1; 0= 1; 0= 1; 0= 1 j i is 99.98% indicating that, with high accuracy, the final state is restricted to this subspace. Figure 5a shows the dependence of the success probability on single-photon loss rate (green). It also presents the success probability when the algorithm is implemented with qubits (red) subjected to dephasing noise (see Methods). Again, we find that, in the presence of equal strength noise, the adiabatic protocol with JPAs (or two-photon-driven KNRs) has superior performance with respect to qubits. Figure 5b also shows the success probability for the same problem but without using the LHZ embedding, that is, when the three KNRs (green) or qubits (red) are directly coupled to each other via a two-body interaction of the formĤ p ¼J P k;j¼1;2;3ŝ z;kŝz;j with J40 (see Methods). As with embedding, the success probability with KNRs is higher than with qubits for equal strength noise. For the particular example considered here, the degeneracy of the ground state is higher in the un-embedded problem (six) than for the embedded problem (three). As a result, the likelihood to remain in one of the ground states increases and, in the presence of noise, the un-embedded problem performs slightly better than the embedded problem. These examples of simple frustrated three-spin problems demonstrate the performance of a single plaquette. Embedding of large Ising problems requires more plaquettes connected together as shown in Fig. 4. Even in such a larger lattice, each JPA is connected to only four other JPAs, making it likely that the final state remains restricted to the encoded subspace spanned by the states 0 j i, 1 j i.

Discussion
We have introduced an adiabatic protocol performing quantum annealing with all-to-all connected Ising spins in a network of non-linear resonators with only local interactions. We have analysed the performances of this annealer in the presence of single-photon loss and shown that the success probability is considerably higher compared to qubits with same amount of loss. Although the implementation of the LHZ scheme has been explored here, other embeddings schemes such as minor embedding could be realized by taking advantage of single-photon exchange and the corresponding two-body couplings that it results in. A distinguishing feature of our scheme is that the spins are encoded in continuous-variable states of resonator fields. The restriction to two approximately orthogonal coherent states only happens in the late stage of the adiabatic evolution, and in general each site must be treated as a continuous variable system displaying rich physics, exemplified by non-Gaussian states, with negative-valued Wigner functions. Because the negativity of the Wigner function is directly related to classical non-simulability [30][31][32] , how this behaviour persists in the presence of photon loss with increasing problem sizes is an interesting question. Another promising avenue to explore is how the continuous variable nature of our system influences the annealer's The success probability for an implementation with qubits with pure dephasing rate g f is also shown (red squares). The two cases are designed to have identical D min and computation time t ¼ 40/D min . The quality factor Q ¼ o r /k is indicated on the top axis for a KNR of frequency o r /2p ¼ 5 GHz. (b) Without encoding: Probability of successfully finding the ground state of a frustrated three-spin Ising problem by implementing the adiabatic algorithm on three directly coupled KNRs with single-photon loss (green squares) for E 0 p ¼2K, d 0 ¼ 0.45K, J k,j ¼ 0.095K for k,i ¼ 1, 2, 3. Note that the local drive J in the embedded problem is same as the coupling J k,j in the un-embedded one and the minimum energy gap in the un-embedded problem is twice that of the embedded problem. The success probability for an implementation with qubits without encoding and with pure dephasing is also shown (red squares).
computational capabilities as compared to more conventional approaches based on two-level systems evolving under a transverse field Ising Hamiltonian, that is, where equation (1) is built from H i ¼ P iŝ x;i and H p ¼ P ij J i;jŝz;iŝz;j (refs 9,33). For instance, we showed how the nature of the quantum fluctuations around the instantaneous ground and excited states leads to increased stability of the ground state. As the size of the system increases, these continuous variable states might alter the nature of phase transitions during the adiabatic evolution, something which may lead to changes in computational power 34,35 . It is also worth pointing out that our circuit QED implementation easily allows for adding correlated phase fluctuations given by interaction terms likeâ y iâ iâ y jâ j (see Supplementary Note 6). These terms do not affect the energy spectrum of the encoded problem Hamiltonian, but may modify the scaling of the minimal gap during the annealing protocol.
Yet another appealing feature that motivates further study is that the time-dependent Hamiltonian of KNRs is generically non-stoquastic in the number basis. A stoquastic Hamiltonian by definition only has real, non-positive off-diagonal entries 36 , and the Hamiltonians in this class are directly amenable to quantum Monte Carlo simulations (stoquastic Hamiltonians do not have the so-called 'sign problem'). As an example, the transverse field Ising Hamiltonian, which is the focus of much current experimental efforts 9,33 , is stoquastic. In contrast, the Hamiltonian of our system has off-diagonal terms P k J k ðâ y k þâ k Þ in the LHZ embedding (or P ij J i;jâiâ y j þ h:c: if this embedding is not used) with problem-dependent signs (note that simply mappingâ k ! Àâ k does not solve the problem due to the presence of the quartic terms equation (2)). The same is true if one considers matrix elements in the over-complete basis of coherent states. Non-stoquasticity has been linked to exponential speedups in quantum annealing 37 , and is widely believed to be necessary to gain more than constant speedup over classical devices 38 .
Ultimately, further investigation into the performance of our adiabatic protocol on larger problem sizes is warranted. Currently, the large Hilbert space size prevents numerically exact simulations with more than a few resonators. Nonetheless, the results here strongly suggest that the adiabatic protocol with two-photon-driven KNRs has excellent resistance to photon loss and thermal noise. Together with the highly non-classical physics displayed during the adiabatic evolution, this motivates the realization of a robust, scalable quantum Ising machine based on this architecture. After completing this work, we became aware of an alternative approach to quantum annealing with Kerr parametric oscillator 39 .

Methods
Eigen-subspace of a two-photon-driven KNR. Following ref. 15, the Hamiltonian of the two-photon-driven KNR can be expressed aŝ This form makes it clear that the two coherent states AE ffiffiffiffiffiffiffiffiffiffiffi E p =K p , which are the eigenstates of the annihilation operatorâ, are also degenerate eigenstates of equation (4) with energy E 2 p =K.
Time-dependent Hamiltonian in the instantaneous rotating frame. We describe the required time-dependence of the amplitude and frequency of the drives to obtain the time-dependent Hamiltonians needed for the adiabatic protocol. As an illustration, consider the example of a two-photon-driven KNR with additional single-photon drive whose Hamiltonian is written in the laboratory frame asĤ Here, o r is the fixed KNR frequency and o p (t) is the time-dependent two-photon drive frequency. The frequency of the single-photon drive, of amplitude E 0 ðtÞ, is chosen to be o p (t)/2 such that it is on resonance with the two-photon drive. In a rotating frame defined by the unitary transformationÛ¼ exp io p ðtÞtâ yâ =2 Â Ã , this Hamiltonian readsĤ 1 ðtÞ¼ÛðtÞ yĤ Lab ðtÞÛðtÞ À i _ UðtÞ yÛ ðtÞ; Choosing the time dependence of the drive frequency as o p (t) ¼ 2o r À 2d 0 (1 À t/2t), and the drive strengths as E p ðtÞ¼E p t=t and E 0 ðtÞ¼E 0 t=t, the above Hamiltonian simplifies tô This has the standard form of a linear interpolation between an initial Hamiltonian and a problem Hamiltonian that is required to implement the adiabatic protocol. As a second illustration, the time-dependent Hamiltonian for finding the ground state of a frustrated three-spin problem embedded on a plaquette iŝ where o r,k are the fixed resonator frequencies and o p,k (t) the time-dependent two-photon drive frequencies. The resonators labelled k ¼ 1, 2 and 3 are driven by time-dependent two-photon and single-photon drives of strengths E p ðtÞ, J(t) and frequency o p,k (t), o p,k (t)/2, respectively. On the other hand, the frequency and strength of the two-photon drive on the k ¼ 4 resonator is fixed.
Estimation of success probability:. To estimate the success probability of the adiabatic algorithm with KNRs, as shown by the green squares in Fig. 3, we numerically simulate 24,25 the master equation _ r¼ ÀĤ 2 ðtÞ;r where photon loss is accounted for by the Lindbladian Dâ i ½ ¼â irâ y i À ðâ y iâir þrâ y iâi Þ=2. It is important to keep in mind that even though the energy gap is small in the rotating frame, the KNRs laboratory frame frequencies o r,k are by far the largest energy scale. As a result, the above standard quantum optics master equation correctly describes damping in this system 40 Here,ŝ z;i andŝ x;i are Pauli operators in the computational basis formed by the ground g j i and excited state e j i of the ith qubit. In these simulations, the qubits are initialized to the ground state of the initial transverse field, and the success probability (red squares in Fig. 3) is computed as the probability of occupation of the correct ground state at t ¼ t, that is, g; e h jrðtÞ g; e j iþ e; g h jrðtÞ e; g j i and g; g h jrðtÞ g; g j iþ e; e h jrðtÞ e; e j i for J40 and Jo0, respectively. Finally, to obtain the data for the resonators in Fig. 5 The success probability is computed as the probability of occupation of the correct ground state at t ¼ t, that is, 0; 1; 0 h jrðtÞ 0; 1; 0 j iþ 1; 0; 0 h jrðtÞ 1; 0; 0 j iþ 0; 0; 1 h jrðtÞ 0; 0; 1 j i(green squares in Fig. 5) and g; e; g h jrðtÞ g; e; g j iþ e; g; g h jrðtÞ e; g; g j iþ g; g; e h jrðtÞ g; g; e j i(red squares in Fig. 5).
Two paradigms of quantum annealing. The KNR-based quantum annealer proposed in this paper is based on an adiabatic non-equilibrium evolution, where the system is subject to driven and dissipative processes. This is in stark contrast to the conventional approach to quantum annealing, where a quantum system is at all times in thermal equilibrium at very low temperature such that it stays close to the ground state, as Hamiltonian parameters are adiabatically varied 1 . It is crucial to understand the very different roles played by the bath, modelling the environment of the annealer, in these two different approaches. In the conventional approach, as long as the temperature is sufficiently low compared to the energy gap, the thermal population of the first excited state is negligible and the system is effectively in its ground state. In fact, for large gaps, the coupling to the environment typically helps the annealer by constantly cooling the system towards its ground state. On the other hand, the bath becomes detrimental and will typically lead to large errors as soon as D $ k B T. Since the gap decreases exponentially with problem size for hard problems, this is a major roadblock for conventional quantum annealing. Even for easier problems when the gap closes polynomially, it quickly becomes extremely challenging to go to large system sizes.
The type of non-equilibrium quantum annealer considered in the present paper overcomes this roadblock, but trades the difficulty for a related but different challenge. The crucial point is that although the gap is small in the rotating frame where the annealing schedule is realized, the system still probes the environment at a very high frequency. All coupling and interaction terms in the Hamiltonian are effectively small perturbations of the resonator's bare HamiltonianĤ 0 ¼ o râ yâ , such that the energy cost of adding a thermal photon is to a very good approximation ' o r , giving a negligible thermal population, N th ¼ exp( À ' o r / k B T), for typical frequencies in the 5-15 GHz range and temperatures TB10 mK. This is also the justification for the master equations used when computing the success probabilities Figs 3 and 5.
Although thermal noise is no longer a bottleneck for this type of non-equilibrium quantum annealing, another challenge now arises. Since the system is not in equilibrium, the eigenstates of the rotating frame Hamiltonian are not global eigenstates of the total system, including the bath, and the interaction with the bath therefore does not generically drive the system towards the rotating-frame ground state, even at zero temperature (see Supplementary Note 3). This leads to local dephasing noise for the KNR implementation due to resonator photon loss. We emphasize that when comparing to a qubit implementation in Figs 3 and 5, we are comparing to an analogous implementation where the qubits are also only subjected to local dephasing noise, as opposed to thermal noise due to a small gap. This allows a fair comparison of two different physical systems used to realize the Ising spins under equal noise strength, and applies for example to realizations of the type proposed in ref. 26. How non-equilibrium quantum annealing compares to conventional equilibrium quantum annealing more generally is to the best of our knowledge an open problem, and an interesting and important avenue for future research.
Realization of four-body coupling. The physical realization of the four-body coupling is described here and more details can be found in Supplementary Notes 6 and 8. The photon annihilation operators of the four KNRs each of which are driven with a two-photon drive are denoted byâ i , with i ¼ 1, 2, 3, 4. These resonators are capacitively coupled to a central JJ described by the annihilation operatorâ c and of energy À E J cos f c =f 0 ð Þâ y c þâ c À Á È É . In this expression, E J is the Josephson energy, f 0 ¼ ' /2e is the reduced flux quantum and f c is the standard deviation of the zero-point flux fluctuation for the junction mode. The coupling strength g i between the resonators and the junction is smaller than the detuning between them D i , g i ( D i . In this dispersive, limit the mode of the junction becomes dressed with the resonator modesâ c !â c À g 1 =D 1 ð Þ a 1 À g 2 =D 2 ð Þ a 2 þ g 3 =D 3 ð Þ a 3 þ g 4 =D 4 ð Þ a 4 and the Josephson energy becomes The fourth-order expansion of the cosine leads to a coupling À Câ y 1â y 2â3â4 þ h:c, where C ¼ E J f 4 c =f 4 0 À Á g 1 g 2 g 3 g 4 =D 1 D 2 D 3 D 4 . In the rotating frame of the drive, this coupling becomes resonant when the frequencies are chosen such that o p,1 (t) þ o p,2 (t) ¼ o p,3 (t) þ o p,4 (t). In addition to the above four-body coupling, cross-Kerr terms of the formâ y iâiâ y jâj are also resonant. As shown in Supplementary Note 7, these terms do not affect the success probability of the algorithm. The strength of the coupling can be estimated with typical parameters E J /2p ¼ 600 GHz, f c ¼ 0.12f 0 , g i =D i $ 0:12, resulting in C/2p ¼ 63 KHz. For a typical strength of Kerr-nonlinearity K/2p ¼ 600 KHz, this leads to C=K $ 0:1.
Data availability. The data that support the findings of this study are available from the corresponding author on reasonable request.