## Abstract

Quantum phenomena have the potential to speed up the solution of hard optimization problems. For example, quantum annealing, based on the quantum tunneling effect, has recently been shown to scale exponentially better with system size than classical simulated annealing. However, current realizations of quantum annealers with superconducting qubits face two major challenges. First, the connectivity between the qubits is limited, excluding many optimization problems from a direct implementation. Second, decoherence degrades the success probability of the optimization. We address both of these shortcomings and propose an architecture in which the qubits are robustly encoded in continuous variable degrees of freedom. By leveraging the phenomenon of flux quantization, all-to-all connectivity with sufficient tunability to implement many relevant optimization problems is obtained without overhead. Furthermore, we demonstrate the robustness of this architecture by simulating the optimal solution of a small instance of the nondeterministic polynomial-time hard (NP-hard) and fully connected number partitioning problem in the presence of dissipation.

- Quantum optimization
- superconducting circuits
- all-to-all connectivity
- decoherence
- continuous variable quantum computation
- nonlinear quantum optics

## INTRODUCTION

Optimization problems are ubiquitous in nature and throughout human activities such as computational biology, combinatorial chemistry, or corporate planning. Consider, for example, the following: Given a set of assets with fixed values, is it possible to share them fairly between two parties? This decision problem is an instance of what is known as the number partitioning problem (NPP), which is nondeterministic polynomial-time hard (NP-hard) (*1*, *2*). This means that any known exact algorithm will take a time exponential in the number of assets to solve at least some instances of this problem. It turns out that this problem, like many others (*3*), is isomorphic to finding the ground-state configuration of an antiferromagnetic long-range Ising model with free energy

Finding a ground state means finding a configuration of spins {*s*_{1}, *s*_{2}, …, *s*_{N}}, which minimizes this energy. To map the NPP onto an Ising problem, for example, we set *J*_{ij} ∝ *n*_{i}*n*_{j}, where *n*_{i} denotes the value of the *i*th asset. The Ising spin *s*_{i} labels which of the two sets in a given partition of the set of assets , the asset *i* is in. If *s*_{i} = +1, then the asset is in *F*, whereas if *s*_{i} = −1, then the asset is in the complement . A fair partition exists if and only if the corresponding ground-state energy is 0; since then, . If no fair partition exists, minimizing the energy will yield the least unfair partition.

The NPP is also noteworthy as an example of a class of optimization problems with dense connectivity graphs that require only limited tunability (*3*). Because the free energy can be expressed as the square of a linear combination of the Ising spins, only tunable parameters are necessary as opposed to the total number of pairwise connections. Furthermore, for the NPP, approximate solutions obtained via classical heuristics, such as simulated annealing (*4*), can be very poor (*5*), making it an ideal test bed for quantum algorithms. The optimization landscapes of hard problems are typically highly nonconvex, and heuristic approaches tend to get trapped in local minima. In 1998, Kadowaki and Nishimori (*6*) introduced the idea that the phenomenon of quantum tunneling could help to escape from local potential minima. This insight and subsequent works (*7*–*9*) have led to high interest in quantum annealing both in academia and in the private sector (*10*–*14*). Although a genuine quantum speedup over the best known classical heuristic algorithms remains to be demonstrated, there are strong indications that this could be achieved in the near future. In particular, for problem instances with high and narrow barriers, the d-wave quantum annealer (QA) has been shown to succeed exponentially faster than thermally activated, simulated annealing (*15*).

### Two main challenges

Many hard optimization problems map onto Ising models with dense connectivity graphs. However, the d-wave architecture has limited connectivity (*15*). This is due to the fact that the interconnections between qubits are realized by physical coupler devices. Because the number of connections grows quadratically with the number of qubits, this approach quickly represents an intractable design challenge. To circumvent this connectivity problem, the standard approach is to use a minor embedding scheme (*16*, *17*) to map a fully connected graph onto the sparse physical graph. However, this comes at the cost of a substantial overhead in the number of physical qubits. More recently, another mapping was proposed, where, instead of the *N* bits of the Ising model, the *N*(*N* – 1)/2 binary pairings are encoded in a two-dimensional lattice of qubits with four-qubit nearest-neighbor interactions (*18*). Implementations of this idea with superconducting circuits have been proposed by Leib *et al*. (*19*) and Chancellor *et al*. (*20*). Also, in this case, a substantial overhead as compared to a direct implementation remains. In contrast, here, we propose a superconducting continuous variable Ising machine (CVIM) with full connectivity and zero overhead, paving the way to experimentally demonstrate quantum supremacy (*21*).

Ideally, quantum annealing relies on the coherent evolution of the ground state of an isolated quantum system, the Hamiltonian of which is varied in time adiabatically (*22*). Typically, the Ising spin *s*_{i} is encoded in the two states of a qubit, and the Hamiltonian is of the form (we set )(1)

The value of the control parameter ε(*t*) is initially chosen such that |ε(*t* = 0)| ≫ max |*J*_{ij}|, and the system is prepared in the ground state of *H*_{A}(0). is then gradually reduced to 0 such that . If the adiabatic condition is satisfied (*22*), the system at time *T* is in the ground state of the quantum Ising model. However, real-world QAs are open quantum systems. It was recognized early on that one of the potential strengths of quantum annealing is its relative robustness to certain types of errors affecting the underlying qubits (*23*–*27*). More precisely, if decoherence takes place only in the instantaneous energy eigenbasis of (1), then it does not decrease the success probability of the optimization as compared with the coherent limit (*23*). However, in general, decoherence takes place in different channels depending on the physical hardware. For example, dephasing errors in the Ising basis, described by the random action of operators, are deleterious to the success of the optimization, because the error operators do not commute with the transverse part of the Hamiltonian (*1*) and hence lead to transitions out of the ground state in the initial stages of the annealing process (*13*). In contrast, we show here that quantum annealing with our continuous variable system is remarkably robust to decoherence. A key result of the present work is that if the annealing rate is below the problem-specific adiabatic threshold but larger than the dissipation rate, then a single run of the continuous variable optimizer will succeed with probability > 0.5 for ramp rates smaller than but arbitrarily close to the threshold value.

## RESULTS

### Continuous variable Ising machine

The binary Ising spin variable *s*_{i} is encoded into the quantized phase of a Kerr parametric oscillator (KPO) (see the Supplementary Materials) (*28*–*32*) above threshold, which can take on two values: 0, corresponding to *s*_{i} = +1, or π, corresponding to *s*_{i} = −1. The dynamics of this system is described by the Hamiltonian

Here, *a* and *a*^{†} are the bosonic annihilation and creation operators, respectively. *K* > 0 is the strength of the Kerr nonlinearity, and ε is the strength of a two-photon drive. If Δ < 0, then at ε = 0, the vacuum is the ground state of the system. As the drive strength is increased, the system undergoes a bifurcation at the threshold value ε_{th} = |Δ| into a superposition of coherent states , with , also called a cat state.

In the studies by Goto (*30*, *31*), a system of *N* KPOs coupled via a term of the form was considered. It was shown, using perturbation theory, that as the two-photon drive strength of each KPO is varied from ε = 0 to ε ≫ |Δ|, the multimode vacuum is adiabatically connected to a multimode cat state of the form(2)where *s*_{1}, …, *s*_{N} ∈ {−1,1} are such that the Ising energy − ∑_{n,m}*J*_{nm}*s*_{n}*s*_{m} is minimized. This system thus presents the opportunity to encode an Ising optimization problem in the adiabatic dynamics of a continuous variable quantum system. To build some intuition, we first consider the simple case of two identical coupled KPOs and denote their operators with *a* and *b*. The mutual coupling then has the form *J*(*a*^{†}*b* + *b*^{†}*a*). In the weakly nonlinear limit *K* ≪ 4|*J*|, the Hamiltonian of this system can be conveniently written in the basis of the symmetric and the antisymmetric modes as (see Materials and Methods)

This corresponds to two KPOs with different frequencies, coupled via a cross-Kerr term, −*Kd*^{†}*dc*^{†}*c*. If Δ ± *J* < 0, then the two-mode vacuum is the ground state at ε = 0. Because of the different frequencies, the thresholds of the two modes are shifted to |Δ + *J*| for the symmetric mode and to |Δ − *J*| for the antisymmetric mode. Hence, as the two-photon drive strength of the two KPOs is increased, the soft mode (*d* if *J* > 0 and *c* if *J* < 0) undergoes a bifurcation before the hard mode (*c* if *J* > 0 and *d* if *J* < 0). Once the soft mode starts to bifurcate and becomes populated with photons, the cross-Kerr coupling makes the hard mode even harder, pushing its bifurcation threshold further away. Hence, the Kerr nonlinearity provides a stabilizing feedback mechanism. This is schematically illustrated in Fig. 1 (A to C).

### Physical implementation

We next turn to the physical implementation of such a machine with superconducting circuits. Prototypes of similar bifurcation-based coherent Ising machines have been built with optical systems (*33*–*36*). However, up until now, no physical realization that implements a dense connectivity graph in a scalable fashion has been proposed. As compared with optical systems, superconducting Josephson circuits offer the crucial advantage of stronger nonlinearities that have recently reached the quantum regime, where the nonlinear frequency shifts are larger than the resonance linewidths (*37*–*39*).

A single KPO can be engineered with superconducting circuits by modulating the flux through a split Josephson junction (Fig. 2) at close to twice its natural resonance frequency (see the Supplementary Materials) (*32*). The Kerr nonlinearity is provided by the Josephson potential expanded to the fourth order in a regime where the Josephson energy exceeds the charging energy.

Although several superconducting KPOs can be coupled inductively, this typically yields only a short-range interaction between nearest neighbors. However, long-range coupling between all pairs of *N* oscillators can be obtained if the oscillators are connected in series and shunted by an inductive element as shown in Fig. 2. Flux quantization imposes a constraint on the sum of all phase drops across the KPOs ({ϕ_{n}}_{n=1,…,N}) and across the shunt (ϕ_{0}) such that . The inductive energy of the shuntthen immediately yields an interaction term between all KPOs. Note that we use units such that the flux quantum . After including the capacitive energies and proceeding with standard circuit quantization (*40*), , where *Z*_{n} denotes the mode impedance, we obtain the effective Hamiltonian in a frame rotating with half the ac flux modulation frequencies Ω_{n} = 2(ω_{n} − Δ) (see the Supplementary Materials)(3)

Here, is the natural frequency of the *n*th KPO expressed in terms of the charging and Josephson energies (*40*) of the Josephson junctions and and the dc part of the flux bias . The two-photon drive strength is , where and is the strength of the ac flux modulation at frequency Ω_{n}. We have assumed |Δ| ≪ ω_{n} and applied the rotating wave approximation to suppress fast rotating terms. In the regime , the strength of the Kerr nonlinearity of oscillator *n* is , and finally, the interaction strength between oscillators *n* and *m* is .

To achieve all-to-all coupling, the *N* KPOs must be made resonant with each other. This can be achieved by tuning the *N* dc flux biases such that ω_{1} ≃ ω_{2} ≃ … ≃ ω_{N}. Note that if the shunt is a conventional inductor, the coupling strengths so far are all positive, implying that only ferromagnetic instances of the Ising model can be accessed. Whereas some nontrivial optimization problems can be mapped onto ferromagnetic Ising models with inhomogeneous longitudinal fields (*13*), others such as the NPP require antiferromagnetic couplings. Moreover, antiferromagnetic couplings can give rise to frustration that is intimately related to spin glass physics (*10*, *41*). Antiferromagnetic coupling can be achieved by substituting the shunt inductor with a large-area Josephson junction and by biasing the loop created by the *N* KPOs and the shunt with half a flux quantum (see Materials and Methods). Under the condition that , where denotes the Josephson energy of the shunt junction, the latter effectively acts as a negative inductor and the coupling matrix elements become (see the Supplementary Materials). We remark that a π-junction (*42*–*44*) shunt provides an alternative to realize antiferromagnetic couplings. Tunability of the matrix elements is enabled by connecting a tunable capacitor in parallel with the split junction of each KPO. Although tunable high-*Q* capacitors in the microwave regime are not yet part of the standard toolbox of circuit quantum electrodynamics, their development is an active area of research (*45*, *46*). This yields tunable parameters, which are sufficient for a number of relevant optimization problems with dense connectivity graphs (*3*). Furthermore, circuits with multiply connected topology can be used to extend tunability (see the Supplementary Materials).

### Robustness to dissipation

The dominant decoherence source in this system is the energy dissipation that occurs because of the internal losses in the device or via the capacitively coupled readout lines (see Fig. 2 and the Supplementary Materials). We describe this dissipation by including photon losses with rate κ using a standard Lindblad master equation(4)

Hereand *H*_{N} is given by Eq. 3. This way of writing the master equation emphasizes the two different aspects of photon loss: The first term on the right-hand side represents a nonunitary but deterministic evolution of the state, whereas the second term on the right-hand side represents the stochastic, that is, nondeterministic, “jump” action of an annihilation operator on the state (*47*).

A key property of the continuous variable Ising encoding is its robustness to photon loss. This robustness can be illustrated with the simple case of two coupled identical KPOs (Fig. 1). The dissipation modifies the thresholds of the soft and hard modes as(see the Supplementary Materials). If the drive strength is varied adiabatically, then under the action of *H*_{NH}, the wave function of the system splits deterministically into an equal superposition between the two maxima of the potential to remain in the instantaneous steady state (*48*). Note that this evolution preserves the purity of the state. Then, in the original basis, the two-mode vacuum evolves into a two-mode cat state with even photon number parity (Fig. 3A)

Far above threshold, the amplitude can be found from the solution of the semiclassical equations of motion (see Materials and Methods)(5)with

Let us now turn to the jump part of the dissipative process (see Eq. 4). Below threshold, the average photon number is low (thin solid curves in Fig. 3). Therefore, the probability of a photon loss event below threshold is strongly suppressed, and consequently, the impact of dissipation on the annealing process is reduced. Above threshold, the jumps lead to the random switching of the photon number parity (Fig. 3B). In the ensemble-averaged picture, this results in dephasing of the pure states into mixtures of two-mode coherent states, that is, (, as well as . Crucially though, photon loss above threshold does not corrupt the Ising spin correlations, that is, for *J* > 0 (*J* < 0), the two oscillators still oscillate in-phase (with opposite phase) (Figs. 1D and 3). This robustness of the annealing process to dissipation is to be contrasted with a conventional discrete qubit-based implementation of a QA, where initially the state of the annealer is a fragile coherent superposition of all possible spin configurations and the impact of qubit dephasing in the Ising basis before the avoided crossing strongly reduces the population of the ground state (see the Supplementary Materials) (*13*).

### Application to number partitioning

To demonstrate the capabilities of our proposed device, we have simulated an instance of the NPP with *N* = 4 oscillators using a Monte Carlo quantum trajectory algorithm (*49*). The NPP is defined by the ordered set *A* = {4,5,6,7}, for which a fair partition is *F* = {4,7} and . To encode the corresponding Ising problem into our circuit, we set the coupling matrix elements as *J*_{ij} = − *J*_{0}*A*_{i}*A*_{j} and *J*_{ii} = 0. Here, *J*_{0} > 0 is a scale factor that leaves the Ising problem invariant and that is convenient to satisfy the physical constraints on the coupling strengths. In section S6, we provide realistic circuit parameter choices compatible with state-of-the-art technology. In the simulations, we choose *J*_{0} = 1/42. The two Ising spin configurations that satisfy this NPP are {*s*_{1} = *s*_{4} = 1, *s*_{2} = *s*_{3} = −1}, and the configuration with all spins flipped. Figure 4A shows the success probability as a function of (*d*ε/*dt*)^{−1} = *T*/ε_{MAX} and the photon loss rate *κ*. Figure 4B shows the corresponding average number of jump events. Each point in both figures represents an average over 40 trajectories. A particular run of the optimizer is deemed successful if the phase correlation between each pair of oscillators, as measured by , is equal to *s*_{i}*s*_{j}, at the final drive strength value ε_{MAX}. For comparison, Fig. 4C shows the success probability, obtained for the same NPP, when simulating a standard QA consisting of four qubits that are subject to dephasing with rate γ in the basis that diagonalizes the Ising terms. The optimization in the CVIM is more robust. Success probabilities above 0.5 are obtained even when the average number of photons lost during the annealing is larger than 1 (see Fig. 4B and the Supplementary Materials). In contrast, the optimization with the spin-based QA fails already after a single dephasing event (see Fig. 4D and the Supplementary Materials). This can be understood by observing that, in the initial state of the annealer, where the spins are polarized along for *i* ∈ {1,2,3,4}, a dephasing error described by induces transitions out of the ground state to excited states that are not adiabatically connected with the solution of the optimization problem (see the Supplementary Materials).

## DISCUSSION

In conclusion, we have proposed and investigated, both analytically and numerically, the implementation of a quantum Ising optimization machine with superconducting circuits, which addresses two of the major challenges in the field. First, we show that flux quantization enables the realization of all-to-all connectivity among the Ising spins, without overhead; that is, *N* oscillators are sufficient to encode *N* Ising spins with full connectivity. Second, we show that, with this continuous variable encoding of the Ising minimization problem, quantum annealing succeeds with high probability in a dissipative regime with high error rate, where conventional discrete qubit-based quantum annealing breaks down. Our results open up new perspectives for quantum optimization. After the present work was made public, a related but independent work was reported (*50*).

## MATERIALS AND METHODS

### Two coupled KPOs

The Hamiltonian of two coupled KPOs based on the symmetric/antisymmetric modes *d* and *c* reads

Given that the frequencies of the symmetric/antisymmetric modes are Δ ± *J*, the terms on the second line rotate with frequency *4J* and can be neglected in rotating wave approximation.

### Semiclassical equations of motion

The semiclassical equations of motion are obtained by replacing the quantum operators *a* and *b* by complex functions in the Heisenberg-Langevin equations of motion. They read

A linear stability analysis of the solutions of these equations is provided in the Supplementary Materials.

### Antiferromagnetic coupling

The potential energy of the system with a large-area Josephson junction shunt and a half-flux quantum flux bias, Φ_{e} = Φ_{0}/2, is

For small oscillations, the potential minimum satisfies the transcendental equation (see the Supplementary Materials)

When , the only solution of this equation is φ_{n} = 0. Expanding the shunt potential around this classical minimum and quantizing yields the antiferromagnetic interaction term given in the main text (see the Supplementary Materials).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/4/e1602273/DC1

section S1. Derivation of the circuit Hamiltonian

section S2. Readout and dissipation

section S3. Two coupled parametric oscillators

section S4. Robustness to decoherence

section S5. Quantum trajectory comparison between CVIM and QA

section S6. Parameter estimates

section S7. Typical quantum trajectories for the NPP

fig. S1. Illustration of a superconducting quantum interference device–based KPO circuit.

fig. S2. Schematics of a multiply connected circuit with logarithmically extended tunability.

fig. S3. Highest part of the eigenspectrum of two coupled parametric oscillators in the rotating frame.

fig. S4. Populations of the seven highest eigenstates as a function of the normalized ramp rate.

fig. S5. Stability diagram for the symmetric mode when *J* > 0.

fig. S6. Logarithmic negativity of the two-mode Gaussian state during annealing.

fig. S7. Conditional success probability conditioned on the total number of jumps.

fig. S8. Conditional success probability conditioned on the total number of jumps (*N* = 6) and conditioned on the number of jumps in a given oscillator.

fig. S9. Effect of dephasing on the discrete qubit QA.

fig. S10. Ratio of ac amplitudes as a function of .

fig. S11. Conditional expectation values of the interoscillator phase correlations for .

fig. S12. Conditional expectation values of the interoscillator phase correlations for .

fig. S13. Conditional expectation values of the interoscillator phase correlations for .

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**The numerical calculations were performed in a parallel computing environment at sciCORE (www.scicore.unibas.ch/) scientific computing core facility at the University of Basel.

**Funding:**S.E.N. and R.P.T. acknowledge financial support from the Swiss National Science Foundation. N.L. was financially supported by the National Centre of Competence in Research Quantum Science and Technology.

**Author contributions:**S.E.N. developed the concepts, carried out the calculations, wrote the numerical code and the manuscript, and contributed to the interpretation of the results. R.P.T. contributed to the calculations, the interpretation of the results, and the manuscript. N.L. contributed to the numerics, the interpretation of the results, and the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2017, The Authors