{'fraction_non_alphanumeric': 0.04975124378109453, 'fraction_numerical': 0.017011715615471032, 'mean_word_length': 4.31084865457202, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum annealing with a network of all-to-all connected, two-photon driven Kerr nonlinear oscillators

Quantum annealing aims to solve combinatorial optimization problems mapped on to Ising interactions between quantum spins. A critical factor that limits the success of a quantum annealer is its sensitivity to noise, and intensive research is consequently focussed towards developing noise-resilient annealers. Here we propose a new paradigm for quantum annealing with a scalable network of all-to-all connected, two-photon driven Kerr-nonlinear resonators. Each of these resonators encode an Ising spin in a robust degenerate subspace formed by two coherent states of opposite phases. The fully-connected optimization problem is mapped onto local fields driving the resonators, which are themselves connected by local four-body interactions. We describe an adiabatic annealing protocol in this system and analyze its performance in the presence of photon loss. Numerical simulations indicate substantial resilience to this noise channel, making it a promising platform for implementing a large scale quantum Ising machine. Finally, we propose a realistic implementation of this scheme in circuit QED.


INTRODUCTION
Many 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. In general, finding the ground state of the Ising Hamiltonian, referred to as the Ising problem, is an NP-hard problem [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 slowly evolved from the nondegenerate 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 behavior 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 onto two levels of a quantum system, i.e. a qubit, and the optimization problem is encoded in the interactions between these qubits. Adiabatic optimization with a variety of physical systems 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 other physical systems with improved resilience to noise. A general Ising problem is defined on a fully connected graph of Ising spins, and efficient embedding of large problems with such long-range interactions in physical systems with local connectivity is a challenge. 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. The physical spins represent the relative configuration of a pair of logical spins and 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 problem independent four-body coupling to enforce the constraints. While this scheme exhibits some intrinsic fault tolerance to weakly correlated errors [15], better decoding strategies are required to enhance its performance in the presence of correlated spin-flip noise [16]. Nevertheless, the simple design requiring only precise control of local fields makes it attractive for scaling to large problem sizes.
In search of a physical platform for quantum annealing that is both scalable and has adequate robustness to noise, we propose to encode the Ising problem in a network of two-photon driven Kerr nonlinear resonators (KNRs). In our scheme, a single Ising spin is mapped onto 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 [17]. A similar encoding has been used to implement classical Ising machines in other systems [18][19][20][21]. In contrast, we propose to realize quantum adiabatic algorithms by encoding a quantum spin in the quasiorthogonal coherent states. The dominant source of error in this system is single-photon loss from the resonators. However, since a coherent state is invariant under the action of the photon jump operator, the encoded Ising spin is stabilized against bit flips. We describe a circuit arXiv:1609.07117v1 [quant-ph] 22 Sep 2016 QED implementation of local magnetic fields and fourbody coupling between such resonators to build a quantum annealing platform, where a fully connected graph of Ising spins is embedded using the LHZ scheme. The adiabatic optimization is carried out by initializing the resonators in the network to vacuum, and varying only single-site drives to evolve them to the correct ground state of the embedded Ising problem. Interestingly, as we demonstrate numerically with a single driven KNR, the probability for the system to jump from the instantaneous ground state to one of the excited states during the adiabatic protocol due to photon loss is greatly suppressed, when compared to conventional qubit implementations with equal noise strengths. This resilience to the detrimental effects of photon loss, combined with easy state initialization and final state detection by homodyne measurement of the resonators' field amplitudes, opens the door to realizing a large scale quantum annealer with favorable noise resistance.

Quantum annealing is executed by evolving a system of N spins under a time-dependent Hamiltonian
whereĤ i is the initial trivial Hamiltonian whose ground state is known, andĤ p is the final Hamiltonian at t = τ which encodes an Ising spin problem: Here,σ z,i = |1 1| − |0 0| is the Pauli-z matrix for the i th spin and J i,j is the interaction strength between the i th and j th spin. For simplicity, we have assumed a linear time dependence of the Hamiltonian in Eq. (1), but more complex annealing schedules can be used. Crucially, the initial and final Hamiltonian do not commute [Ĥ i ,Ĥ p ] = 0. The system, initialized to the ground state ofĤ i , adiabatically evolves to the ground state of the problem Hamiltonian,Ĥ p , at time t = τ if ∆ min τ 1, where ∆ min is the minimum energy gap [6]. In most existing schemes, the binary states of a qubit represent an Ising spin and the initial Hamiltonian is given by local transverse magnetic fields [22]. Here we show how quantum adiabatic annealing can be implemented by mapping the spin state {|0 , |1 } on two coherent states {|−α , |α } of a KNR. We, moreover, show that the time-dependent Hamiltonian can be implemented through simple single-site tunable drive fields.
Single spin in a two-photon driven KNR The Hamiltonian of a two-photon driven KNR in the frame rotating at the frequency of the drive is given by,Ĥ 0 = −Kâ †2â2 + E p (â †2 +â 2 ), where K is the Kerrnonlinearity and E p is the strength of the two-photon drive. The coherent states |±α 0 are eigenstates of the photon annihilation operatorâ|±α 0 = ±α 0 |±α 0 and are stabilized in such a resonator with α 0 = E p /K [17] (see also Methods). Intuitively, this is seen from the metapotential obtained by replacing the operatorsâ and a † with the complex classical variables x + iy and x − iy in the expression forĤ 0 [23]. As shown in Fig. 1(a), this metapotential features an inverted double well with two peaks of equal height at (± E p /K, 0). These are the two stable points in the metapotential (see Supplementary Note 1). This is consistent with the quantum picture according to which the coherent states |±α 0 are two degenerate eigenstates ofĤ 0 with eigenenergy E 2 p /K [17] (see also Methods). Having a well-defined two-state subspace, we choose to encode an Ising spin {|0 , |1 } in the stable states {|−α 0 , |α 0 }. Importantly, if the rate of single photon loss (κ) is small (κ 8E p ), this mapping remains robust against single-photon loss from the resonator [17]. Moreover, the photon jump operatorâ leaves the coherent states invariantâ|0/1 = ±α 0 |0/1 . As a result, if the amplitude α 0 is large such that 0 /1|â|1/0 = ∓α 0 e −2|α0| 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 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 strength E 0 ,Ĥ p = −Kâ †2â2 +E p (â †2 +â 2 )+E 0 (â † +â). The metapotential for this Hamiltonian, illustrated in Fig. 1(b) for small E 0 , is an asymmetric inverted double well with peaks of unequal heights at (±α 0 , 0). Depending on if E 0 > 0 or E 0 < 0, the peak at (−α 0 , 0) is lower than the one at (α 0 , 0) or vice versa. These two states remain stable, but have different energies, indicating that the small single-photon field induces an effective magnetic field on the Ising spins {|0 , |1 }. Indeed, a full quantum analysis shows that if |E 0 | 4K|α 0 | 3 , then |±α 0 remain the eigenstates ofĤ p but their degeneracy is lifted by 4E 0 α 0 [17]. In other words, in the spin subspace,Ĥ p can be expressed asĤ p = 2E 0 α 0σz + const., withσ z = |1 1 | − |0 0 |, which is the required problem Hamiltonian for a single spin in a magnetic field. Note that, if E 0 increases then the eigenstates can deviate from coherent states (see Supplementary Note 2). Choosing |E 0 | 4K|α 0 | 3 ensures that {|0 , |1 } are indeed coherent states, so that 0 /1|â|1/0 ∼ 0 and the encoded subspace is well protected from the photon loss channel.
In correspondence with Eq. (1), we require an initial Hamiltonian which does not commute with the final problem Hamiltonian and has a simple non-degenerate ground state. This is achieved by introducing a finite detuning δ 0 > 0 between the drives and resonator. In a frame rotating at the frequency of the drives, the initial Hamiltonian is chosen asĤ i = δ 0â †â − Kâ †2â2 with δ 0 < K. In this frame, the ground and first excited states are the vacuum |0 and single-photon Fock state |n = 1 respectively, which are separated by an energy gap δ 0 . If a single-photon is lost from the resonator, the excited state |n = 1 decays to the ground state |0 which, on the other hand, is invariant to photon loss. Since it is simple to prepare in the superconducting circuit realizations 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 single-photon drive strengths and detuning so thatĤ 1 (t) = (1 − t/τ )Ĥ i + (t/τ )Ĥ p realizing Eq. (1) for a single-spin (see Methods). Note that the form of the HamiltonianĤ 1 (t) conveniently ensures that the nonlinear Kerr term is time-independent so that one only needs to vary the drives. By adiabatically controlling the frequency and amplitude of the drives it is possible to evolve the state of the KNR from the vacuum |0 at t = 0, to the ground state of a single Ising spin in a magnetic field at t = τ . Figure 2(a) shows the change of the energy landscape in time found by numerically diagonalizing the instantaneous Hamiltonian,Ĥ 1 (t) for E p = 4K, α 0 = 2, E 0 = 0.2K and δ 0 = 0.2K. The minimum energy gap ∆ min at the avoided level crossing is also indicated. As illustrated by the plots of the Wigner functions in the inset of Fig. 2(a), a resonator initialized to the vacuum state at t = 0 evolves through highly non-classical and non-Gaussian states, towards the ground state |0 at t = τ , with τ ∼ 30/∆ min in this example. If, on the other hand, the KNR is initialized to single-photon Fock state at t = 0 then it evolves to the first excited state |1 at t = τ . The average probability to reach the correct ground state is 99.9% for both E 0 > 0 and E 0 < 0. The 0.1% probability of erroneously ending in the excited state arises from non-adiabatic errors and can be decreased by increasing the evolution time. For example, if τ = 60/∆ min then the success probability increases to 99.99%.

Effect of single-photon loss
An appealing feature of this implementation is that, at the start of the adiabatic protocol at t = 0, the ground state (vacuum) is invariant under single-photon loss. Similarly, at the end of the adiabatic protocol at t = τ , irrespective of the problem Hamiltonian (i.e., E 0 > 0 or E 0 < 0) the ground state (coherent states |0 or |1 ) is also invariant under single-photon loss. It follows that towards the beginning and end of the protocol, photon loss will not induce any errors. Moreover, we find that, even at intermediate times 0 < t < τ , the ground state ofĤ 1 (t) remains largely unaffected by photon loss. This can be understood intuitively from the distortion of the metapotential, as shown in Fig. 2(b) for the example depicted in Fig. 2(a) at t = 0.25τ . The metapotential still shows two peaks, however, 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. Quantitatively, the effect of single-photon loss is seen by numerically evaluating [24,25] the transition matrix elements ψ g (t)|â|ψ e (t) , ψ e (t)|â|ψ g (t) for the duration of the protocol, where |ψ g (t) and |ψ e (t) are the ground and excited state ofĤ 1 (t) respectively. As shown in Fig. 2(c) the transition from the ground to excited state is greatly suppressed throughout the whole adiabatic evolution. This asymmetry in the transition rates distinguishes the adiabatic protocol described here 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
Consider the problem of two interacting spins, which can be embedded in a system of two linearly coupled KNRs, each driven by a two-photon driveĤ p = Transition matrix elements between the ground |ψg(t) and excited states |ψe(t) in the event of a photon jump during the adiabatic protocol.
Here J 1,2 is the strength 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 asĤ p = 4J 1,2 |α 0 | 2σ z,1σz,2 + const., which is the required problem Hamiltonian [17]. The nature of the interaction, that is, ferromagnetic or anti-ferromagnetic, is encoded in the phase of the coupling J 1,2 < 0 or J 1,2 > 0, respectively. For the initial Hamiltonian, we takeĤ i = (1), the full time-dependent Hamiltonian for the twospin problem isĤ 2 (t) = (1 − t/τ )Ĥ i + (t/τ )Ĥ p . 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 during the adiabatic evolution.
tion is indeed restricted to this computational subspace.
In the presence of single-photon loss, the coherence between the states is reduced. However the success probability (see Methods) to solve the Ising problem, which depends only on the diagonal elements of the density matrix (e.g. 0 ,1|ρ(τ )|0,1 ) remains high. For instance, with a large loss rate κ = 50/τ , the fidelity with respect to the superposition state N (|0,0 + |1,1 ) or N (|0,1 + |1,0 ) decreases to 37.6%, but the average success probability of the algorithm is 75.2%.
To characterize the effect of noise, a useful figure of merit is the ratio of the minimum energy gap to the loss rate (∆ min /κ). The dependence of the average success probability on this ratio is presented in Fig. 3 when the algorithm is implemented using KNRs (red squares) with single-photon loss (κ) or qubits (blue squares) with pure dephasing (γ φ ). The success probability is averaged over all instances of the problem of two coupled spins (i.e., ferromagnetic and anti-ferromagnetic). All points are generated by varying the loss rates, while keeping ∆ min and τ = 20/∆ min fixed. In the presence of pure dephasing, the success probability with qubits saturates to 50% for large γ φ . 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, with KNRs, in the presence of photon loss the rate at which the instantaneous ground state jumps to the excited state (∝ ψ e (t)|â|ψ g (t) ) is small compared to the rate at which the instantaneous excited state jumps to the ground state (∝ ψ g (t)|â|ψ e (t) ). For example, the success probability is ∼ 75% even when ∆ min /κ ∼ 1. This shows that the algorithm implemented with twophoton driven KNRs has superior performance compared to that implemented with qubits in the presence of equal strength noise.
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 nonlocal Hamiltonian is a challenging hardware problem, but it may be solved by embedding schemes that map the Ising problem on a graph with only local interactions. As mentioned earlier, the LHZ scheme [14] is one such technique in which the relative configuration of pairs of N logical spins is mapped on M = N (N − 1)/2 physical spins. A pair of logical spins in which both the spins are aligned |1, 1 or |0, 0 (or anti-aligned |1, 0 or |0, 1 ) 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 spins J k (k = 1, ..M ). For a consistent mapping, M − N + 1 energy penalties in the form of four-body coupling are introduced which enforce an even number of spin-flips around any closed loop in the logical spins. It was shown in Ref. [14] that a fully connected graph can then be encoded in a planar architecture with only local connectivity. The problem Hamiltonian in the physical spin basis becomesĤ LHZ,N p = M k=1 J kσz,k − i,j,k,l Cσ z,iσz,jσz,kσz,l where i, j, k, l denotes the nearest-neighbor spins enforcing the constraint.
We now describe a circuit QED platform implementing the LHZ scheme by embedding the physical spins in the eigenbasis {|0 , |1 } of two-photon driven KNRs. Such a resonator can be realized as a microwave resonator terminated by a flux-pumped SQUID. The non-linear inductance of the SQUIDs induces a Kerr non-linearity, and a two-photon drive is introduced by flux-pumping at twice the resonator frequency. This is the exact same setup as is used for Josephson parametric amplifiers (JPAs), and we will therefore refer to this implementation of a Kerr nonlinear resonator as a JPA in the following [17,[27][28][29]. We envision the quantum annealing platform to be built with a group of four JPAs of frequencies ω r,i (i = 1, 2, 3, 4) coupled to a single Josephson Junction (JJ) as shown in Fig. 4(a). To obtain the time-dependent twophoton 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 the resonator frequency, ω p,i (t) 2ω r,k (see Methods and Supplementary Note 4). Additional single-photon drives whose amplitude and frequency can be varied in time are also applied to each of the JPAs to provide the effective local magnetic field. Local four-body couplings are realized by the nonlinear inductance of the central JJ, see Supplementary Note 4. Choosing ω p,1 (t) + ω p,2 (t) = ω p,3 (t) + ω p,4 (t) and taking the resonators to be detuned from each other, the central JJ induces a coupling of the form −C(â † 1â † 2â 3â4 + 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. Other circuits with tunable four-body coupling are possible, Supplementary Note 6. This 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. 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 . As an illustration, Fig. 4(b) depicts all the possible interactions in an Ising problem with N = 5 logical spins and Fig. 4(c) shows the corresponding triangular network of coupled KNRs. A final necessary component for a quantum annealing architecture is readout of the state of the physical spin. Here, this is realized by homodyne detection which can resolve the two coherent states | ± α 0 allowing 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, H p = J k,j=1,2,3σ z,kσz,j with J > 0. 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 3 JPAs) and one physical spin fixed to up state (in our case a JPA initialized to the stable eigenstate |1 ). Since, The anti-ferromagnetic coupling between the logical spins is represented by the single-photon drives on each JPA with amplitude J > 0. At t = 0, the ground state of this Hamiltonian is the vacuum |0, 0, 0 . If the four-body coupling is weak then the problem Hamiltonian can be expressed asĤ p = 2Jα 0 3 k=1σ z,i − 2C|α 0 | 4σ z,1σz,2σz,3σz,4 +const., with α 0 = E p /K. This realizes the required problem Hamiltonian in the LHZ scheme.
To illustrate the performance of this protocol, we numerically simulate the evolution subjected to the Hamiltonian in Eq. (2) with the three resonators initialized to vacuum and the fourth to the state |1 . With E p = 2K, α 0 = √ 2, J = 0.095K, C = 0.05K, τ = 40/∆ min and κ = 0, we find that the success probability to reach the ground state is 99.3%. The reduction in fidelity arises FIG. 5. Success probability for the frustrated threespin problem. Probability of successfully finding the ground state of a frustrated three-spin Ising problem by implementing the adiabatic algorithm on a plaquette of four KNRs with single photon loss (green squares) for E 0 p = 2K, δ0 = 0.45K, C = 0.05K, J = 0.095K. The success probability for an implementation with qubits with pure dephasing rate γ φ is also shown (red squares). The two cases are designed to have identical ∆min and computation time τ = 40/∆min. The quality factor Q = ωr/κ is indicated on the top axis for a KNR of frequency ωr/2π = 5 GHz.
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 is 99.98% indicating that the final state is indeed restricted to this subspace. Figure 5 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 the adiabatic protocol with JPAs (or two-photon driven KNRs) has superior performance in the presence of equal strength noise.
This example of a simple frustrated three-spin problem demonstrates the performance of a single plaquette.
Embedding of large Ising problems would require more plaquettes connected together as shown in Fig. 4. Nevertheless, even in 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 , |1 . Notably, the fully connected Ising problem is realized in the frame rotating at the drive frequencies, which is typically of the order of 5-10 GHz in experimental implementations. The thermal fluctuations around these frequencies are negligible in superconducting circuits operating at 10-30 mK. As a result thermal excitations in the reservoir are unlikely to drive the system out of the ground state even if the minimum energy gap is small.

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 local interactions. The 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 through non-Gaussian states, with negative-valued Wigner functions. How this behaviour persists in the presence of photon loss as the problem sizes are scaled up is an interesting question, as negativity of the Wigner function is directly related to classical non-simulability [30][31][32].
An intriguing question is how the continuous variable nature of our system influences the annealer's computational capabilities when compared to a more conventional approach based on two-level systems evolving under a transverse field Ising Hamiltonian, i.e., where Eq. (1) is built from H i = iσ x,i and H p = i,j J i,jσz,iσz,j [9,33]. For instance, we showed how the nature of the of 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, potentially leading to computational speedups [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â † iâ iâ † jâ j (see Supplementary Note 5). 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.
Another appealing feature that motivates further study into the complexity of our protocol is that the time-dependent Hamiltonian we use is generically non-stoquastic in the number basis. A stoquastic Hamiltonian by definition only has real, non-positive off-diagonal entries [36], and the significance of this is that 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 is stoquastic. In contrast, our Hamiltonian has off-diagonal terms k J k (â † k +â k ) in the LHZ embedding (or ij J i,jâiâ † j + h.c. if this embedding is not used) with problem dependent signs (note that a simple diagonalization does not solve the problem due to the presence of quartic terms). The same is true if one considers matrix elements in the over-complete coherent state basis, α|Ĥ(t)|β . It therefore does not appear straightforward to adapt quantum Monte Carlo techniques to this system.
Ultimately, further investigation into the performance of our adiabatic protocol on larger problem size 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.

Eigen-subspace of a two-photon driven KNR:
Following Ref. [17], the Hamiltonian of the two-photon driven KNR can be expressed aŝ This form form makes it clear that the two coherent states |±α 0 with α 0 = E p /K, which are the eigenstates of the annihilation operatorâ, are also degenerate eigenstates of Eq. (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 timedependent Hamiltonians needed for the adiabatic evolution. As an illustration, consider the example of a twophoton driven KNR with additional single-photon drive whose Hamiltonian is written in the laboratory frame aŝ Here, ω r is the fixed KNR frequency and ω p (t) is the timedependent two-photon drive frequency. The frequency of the single-photon drive, of amplitude E 0 (t), is chosen to be ω p (t)/2 such that it is on resonance with the twophoton drive. In a rotating frame defined by the unitary transformationÛ = exp[iω p (t)tâ †â /2], this Hamiltonian readŝ Choosing the time dependence of the frequency as ω p (t) = 2ω r − 2δ 0 (1 − t/2τ ), and drive strengths as E p (t) = E p t/τ and E 0 (t) = E 0 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 is (see Supplementary Note) where ω r,k are the fixed resonator frequencies, ω p,k (t) are 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 ω p,k (t), ω 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 the master equationρ = −[Ĥ 2 (t),ρ] + κD[â 1 ] + κD[â 2 ], where the photon loss is accounted for by the Lindbladian D[ [24,25]. It is important to keep in mind that even though the energy gap is small in the rotating frame, the KNRs laboratory frame frequencies ω r,k are by for the largest energy scale. As a result, this standard quantum optics master equation correctly describes damping in this system [37]. Moreover, because we are working with KNR frequencies in the GHz range, as is typical in supercondcuting circuits, thermal fluctuations are negligible. From his master equation, the success probability can be evaluated as the probability of occupation of the correct ground state at the final time t = τ , that is, 0 ,1|ρ(τ )|0,1 + 1 ,0|ρ(τ )|1,0 and 0 ,0|ρ(τ )|0,0 + 1 ,1|ρ(τ )|1,1 for E 0 > 0 and E 0 < 0, respectively.

SUPPLEMENTARY NOTE 1: STABILITY ANALYSIS OF THE META-POTENTIAL FOR A SINGLE TWO-PHOTON DRIVEN KNR
The coherent states |±α 0 , where α 0 = E p /K, are eigenstates of the two-photon driven KNR Hamiltonian [1] H 0 = −Kâ †2â2 + E p (â †2 +â 2 ). (S1) Following Ref. [1], we will outline a proof of this statement that will serve as the foundation for the analysis in the Supplementary Notes 2 and 3. We start by applying a displacement transformation D(α) = exp(αâ † − αâ) toĤ 0 , so that In the above expression, we have dropped the constant term E(α) = −K|α| 4 + E * p α 2 + E p α * 2 which represents a shift in energy. We choose α such that the coefficient ofâ † satisfies −2Kα 2 α * + 2E p α * = 0, which is equivalent to finding the turning points of the metapotential of Fig. 1 of the manuscript. This equation has three solutions: (0, 0), (±α 0 , 0) corresponding to the dip and two peaks of the inverted double-well metapotential. At (0, 0), the Hamiltonian in Eq. (S3) represents a resonantly driven parametric amplifier. Therefore, in the absence of losses, large fluctuations make the system unstable around (0, 0) [2]. On the other hand, at (±α 0 , 0), the Hamiltonian Eq. (S3) takes the form With this normally ordered form, we immediately conclude that the vacuum |0 is an eigenstate ofĤ 0 in the displaced frame. It immediately follows that the coherent states |±α 0 are the eigenstates of the HamiltonianĤ 0 in the non-displaced frame. Furthermore, these are degenerate eigenstates as E(α 0 ) = E(−α 0 ).

SUPPLEMENTARY NOTE 2: THE EIGEN-SUBSPACE IN THE PRESENCE OF SINGLE-PHOTON DRIVE
The Hamiltonian of a two-photon driven KNR with an additional single-photon drive is given bŷ Under a displacement transformation D(α) this Hamiltonian readŝ where we have again dropped the constant term E(α) = −K|α| 4 + E p (α 2 + α * 2 ) + E 0 (α + α * ) representing a shift in energy. Following the same steps as above, the coefficient of theâ † terms vanish if For small E 0 , this equation has three solutions of the form (±α 0 + , 0), ( , 0) with = E 0 /4E p . In practice, we assume that the amplitude of the single-photon drive is small compared to that of the two-photon drive, E 0 E p , so that → 0. If the condition Eq. (S7) is satisfied, then the Hamiltonian in the displaced frame reduces tô Following Supplementary Note 1, |0 is an eigenstate ofĤ except for the first term which represents an off-resonant parametric drive of strength |E 0 /2α * | detuned by 4K|α| 2 . If |E 0 /α| 4K|α| 2 , fluctuations around α are small and |0 remains an eigenstate in the displaced frame. Of the three solutions to Eq. (S7), only (±α 0 + , 0) satisfy the condition |E 0 /(±α 0 + )| 4K| ± α 0 + | 2 . The third solution ( , 0) is unstable because of the large quantum fluctuations around this point. As a result, in the non-displaced frame, the eigenstates of the system are |α 0 + and |−α 0 + , where is a small correction. From the expression for E(α), it also clear that the degeneracy between the eigenstates |α 0 + and |−α 0 + is lifted by an amount E(α 0 ) − E(−α 0 ) = 4E 0 α 0 . Increasing E 0 , leads to squeezing due to the first term of Eq. (S8). As a result, the states are no longer coherent or, in other words, they are no longer the eigenstates of the photon annihilation operatorâ and hence are no longer protected against the single-photon loss channel. To illustrate this effect quantitatively, we numerically diagonalize the Hamiltonian in Eq. (S5) to evaluate the ground state for a fixed two-photon drive strength E p = 2K and variable strength single-photon drive amplitude. We numerically solve the master equation for the Hamiltonian in Eq. (S5), ρ = −i[Ĥ,ρ] + κ[âρâ † − (â †âρ +ρâ †â )/2] with the system initialized to the ground state at t = 0 and κ = 0.2K [3,4]. Finally, from these results, we compute the probability of the system to remain in the ground state ψ g |ρ|ψ g after time t, presented in Supplementary Figure 1. As seen from these numerical results, the ground state is invariant to single-photon loss for small E 0 , confirming our prediction from the simple theoretical analysis.

SUPPLEMENTARY NOTE 3: EFFECT OF SINGLE PHOTON LOSS DURING THE ANNEALING PROTOCOL FOR A SINGLE DRIVEN KNR
In this section, we provide a simple theoretical analysis to understand the behaviour of the instantaneous ground and excited state during the annealing protocol with a single driven KNR in the presence of single-photon loss. The time-dependent Hamiltonian of this system is given bŷ In the presence of single-photon loss, the dynamics of the system is described by the master equationρ = −i[Ĥ 1,eff (t)ρ− ρĤ † 1,eff (t)] + κâρâ † , where κ is rate of photon loss andĤ 1,eff (t) =Ĥ 1 (t) − κâ †â /2. Following Supplementary Notes 1 and 2, we apply a displacement transformation D(α) toĤ 1,eff (t) For convenience, the variables have been redefined as δ(t) = δ 0 (1 − t/τ ), E p (t) = E p t/τ and E 0 (t) = E 0 t/τ . Again, we take α to satisfy such that the effective Hamiltonian readŝ Eq. (S11) admits three solutions: ∼ (0, 0), (±α 0 , 0), with α 0 = [2E p (t) + δ(t)]/2K. Repeating the procedure described in Supplementary Note 1 and 2 we find, at short times when δ(t) 2E p (t), that |0 is approximately the lower energy, stable eigenstate of the Hamiltonian. As the evolution proceeds, the strength of the detuning and two-photon drive is modified as δ(t) 2E p (t). In this case (0, 0) is not a stable eigenstate. On the other hand, at (±α 0 , 0) the Hamiltonian readŝ (S13) In the absence of the first term, |0 would be the eigenstate of the above Hamiltonian and hence the coherent states |±α 0 would be the eigenstates in the non-displaced frame. The first term represents a parametric drive whose amplitude has a magnitude ζ = | − (δ(t) − iκ/2)(α 0 /2α * 0 ) − E 0 /α * 0 | and is detuned by | − 4K|α 0 | 2â †â + δ(t)â †â − iκ/2|. In other words, the effect of detuning, photon-loss and single-photon drive is to squeeze the fluctuations around (±α 0 , 0). When δ(t) < 2E p (t), κ < 8E p (t) and E(t) < 4K|α 0 | 3 then squeezing is negligible and the eigenstates in the non-displaced frame are approximately coherent states |±α 0 . If E 0 > 0, then (−α 0 , 0) corresponds to the lower energy state, and, on the other hand, if E 0 < 0, (α 0 , 0) corresponds to the lower energy state. Note that the strength of squeezing for the lower energy state |ζ| = | − [δ(t) − |E(t)|/ 2E p (t) + δ(t)/2K] − iκ/2| is smaller than that for higher energy state |ζ| = | − [δ(t) + |E(t)|/ 2E p (t) + δ(t)/2K] − iκ/2| and hence, the deviation of lower energy state from the coherent state is smaller than the deviation of the higher energy state from the coherent state. This forms the basis for the observation that during the adiabatic evolution, the lower energy state which is the computational ground state is more stable against photon-jump operation. There are eight fully-connected Ising problems for three logical spins with weights J ij = ±1: one where all spins are ferromagnetically coupled, one where all are anti-ferromagnetically coupled, three problems where two spins are ferromagnetically coupled, while the coupling to the third is anti-ferromagnetic and three problems where two spins are anti-ferromagnetically coupled, while the coupling to the third is ferromagnetic. All these problems can be encoded on a single plaquette by appropriately choosing the sign of the single photon drive. For example, as shown in the manuscript, the problem of all the logical spins coupled anti-ferromagnetically is embedded on the KNRs by applying single-photon drives such that J > 0. Similarly, to embed a problem where all logical spins are coupled ferromagnetically, all the applied single-photon drives must be such that J < 0. In order to embed the problem correctly in the LHZ scheme, the four-body coupling must be larger than local magnetic fields [5] and in our case, this implies C|α 0 | 4 > J|α 0 |. The adiabatic algorithm, as described in the manuscript, can be applied once the problem is correctly embedded. the dependence of the average success probability for all fully connected problems on a single plaquette when the adiabatic protocol is implemented using KNRs (green squares) characterized by single photon loss rate κ and qubits (red squares) characterized by the dephasing rate γ φ . The total computation time τ is the same for both cases.
For example, when all the logical spins are coupled anti-ferromagnetically, the ground state is sixfold degenerate: |0, 0, 1 , |0, 1, 0 , |1, 0, 0 , |1, 1, 0 , |1, 0, 1 and |0, 1, 1 . In the basis of physical spins, which map the relative configuration of the logical spins, the ground state is threefold degenerate |0,0,1 , |0,1,0 and |1,0,0 . Supplementary Figure 2(a) shows the evolution of the energy spectrum during the adiabatic protocol with E p = 2K and C = 0.05K. The single photon drive to all the resonators is J = 0.095K. The energy is measured from the ground state and the system starts from a non-degenerate ground state at t = 0. As expected, at t = τ , the ground state becomes threefold degenerate.
If the logical spins are coupled ferromagnetically, then there are two possible ground states: |0, 0, 0 and |1, 1, 1 , which means that in the physical spin basis there exists single non-degenerate ground state: |1,1,1 . Supplementary  Figure 2(b) shows the evolution of the energy spectrum during the adiabatic protocol with E p = 2K and C = 0.05K but in this case the single photon drive is J = −0.095K on all the resonators. The system starts from a non-degenerate ground state at t = 0 and indeed, at t = τ , the ground state remains non-degenerate.
Similarly, it is possible to encode all eight problems on the plaquette and we find that in the absence of single-photon loss the average success probability to find the correct ground state with the adiabatic algorithm is 99.7% in time τ = 200/K. Supplementary figure 3 shows the dependence of the success probability on the rate of single-photon loss. It also presents the success probability of the protocol implemented using qubits with dephasing rate γ φ . The time-dependent Hamiltonian for qubits is designed to have the same minimum energy gap as that with the JPAs and the duration of the protocol is also chosen to be the same. Clearly the adiabatic protocol with the JPAs outperform that with qubits in the presence of equal strength noise.

SUPPLEMENTARY NOTE 4: PHYSICAL REALIZATION OF A PLAQUETTE
To implement the LHZ scheme [5] in our proposed setup, Josephson parametric amplifiers (JPAs) are arranged in a triangular lattice and coupled together a single coupling Josephson junction (JJ). As described in the main text, a single plaquette is comprised of four JPAs and the coupling JJ. The JPA is realized by embedding a SQUID in a resonator and modulating the flux through the SQUID at a frequency close to twice the frequency of the resonator [1,[6][7][8].
Recalling the Hamiltonian of the adiabatic protocol in the manuscript (Eq. (2) and Methods section in the manuscript), the flux modulation frequency is linearly varied from ω p,k (0) = 2ω r,k − 2δ 0 at t = 0 to ω p,k (τ ) = 2ω r,k − δ 0 at t = τ . With an additional single-photon drive, the Hamiltonian of the k th JPA can then be written aŝ H JPA,k = ω r,kâ † kâ k − Kâ †2 kâ 2 k + J(t)[e −iωp,k(t)t/2â † + e iωp,k(t)t/2â ] + E p (t)[e −iωp,k(t)tâ †2 + e iωp,k(t)tâ2 ]. (S14) The Kerr nonlinearity K and frequency ω r,k are determined by the charging and Josephson energy of the junctions in the SQUID. By choosing appropriate parameters it is therefore possible to realize JPAs with different frequencies.
The two-photon drive E p depends on the magnitude of the flux modulation. Consider Supplementary Figure 4 which shows four JPAs of different frequencies (indicated by different colors) capacitively coupled to a single JJ, which, as mentioned above, corresponds to a single plaquette and, thus, is the building block for the LHZ scheme. The