Robust quantum optimizer with full connectivity

See allHide authors and affiliations

Science Advances  07 Apr 2017:
Vol. 3, no. 4, e1602273
DOI: 10.1126/sciadv.1602273


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


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 energyEmbedded Image

Finding a ground state means finding a configuration of spins {s1, s2, …, sN}, which minimizes this energy. To map the NPP onto an Ising problem, for example, we set Jijninj, where ni denotes the value of the ith asset. The Ising spin si labels which of the two sets in a given partition of the set of assets Embedded Image, the asset i is in. If si = +1, then the asset is in F, whereas if si = −1, then the asset is in the complement Embedded Image. A fair partition exists if and only if the corresponding ground-state energy is 0; since then, Embedded Image. 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 Embedded Image tunable parameters are necessary as opposed to the Embedded Image 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 (79) have led to high interest in quantum annealing both in academia and in the private sector (1014). 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 si is encoded in the two states of a qubit, and the Hamiltonian is of the form (we set Embedded Image)Embedded Image(1)

The value of the control parameter ε(t) is initially chosen such that |ε(t = 0)| ≫ max |Jij|, and the system is prepared in the ground state of HA(0). is then gradually reduced to 0 such that Embedded Image. 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 (2327). 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 Embedded Image 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.


Continuous variable Ising machine

The binary Ising spin variable si is encoded into the quantized phase of a Kerr parametric oscillator (KPO) (see the Supplementary Materials) (2832) above threshold, which can take on two values: 0, corresponding to si = +1, or π, corresponding to si = −1. The dynamics of this system is described by the HamiltonianEmbedded Image

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 Embedded Image, with Embedded Image, also called a cat state.

In the studies by Goto (30, 31), a system of N KPOs coupled via a term of the form Embedded Image 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 Embedded Image is adiabatically connected to a multimode cat state of the formEmbedded Image(2)where s1, …, sN ∈ {−1,1} are such that the Ising energy − ∑n,mJnmsnsm 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(ab + ba). In the weakly nonlinear limit K ≪ 4|J|, the Hamiltonian of this system can be conveniently written in the basis of the symmetric Embedded Image and the antisymmetric Embedded Image modes as (see Materials and Methods)Embedded Image

This corresponds to two KPOs with different frequencies, coupled via a cross-Kerr term, −Kddcc. If Δ ± J < 0, then the two-mode vacuum Embedded Image 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).

Fig. 1 Illustration of the bifurcation-based annealing process in two antiferromagnetically coupled KPOs.

(A to C) Evolution of the two-dimensional potential landscape with increasing two-photon drive strength ε. The negative curvature is a consequence of the negative sign of the Josephson-induced Kerr nonlinearity. The antiferromagnetic case (J < 0), where the antisymmetric mode ab is softer, and thus has a lower bifurcation threshold, than the symmetric mode a + b, is shown. Consequently, the system evolves from the vacuum Embedded Image at ε = 0 to the two-mode cat state Embedded Image at large ε, as shown in (D), where we plot the relevant part of the eigenspectrum in the rotating frame (see the Supplementary Materials). Photon loss events predominantly take place above threshold and induce transitions between the cat states Embedded Image and Embedded Image, as indicated by the red arrows in (D).

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 (3336). 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 (3739).

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.

Fig. 2 Schematics of the proposed superconducting CVIM.

It consists of a chain of KPOs shunted by an effective inductor, Leff. Each KPO is realized by a capacitively shunted split Josephson junction threaded by an ac flux bias at twice the resonance frequency. The Ising spin variables are encoded in the quantized oscillation phases (either 0 or π) of the KPOs above threshold. The inductive shunt induces all-to-all coupling between the KPOs. To obtain antiferromagnetic coupling, a large-area Josephson junction can be used as a shunt together with a flux bias of Φe = Φ0/2. Homodyne readout of the oscillator phases is enabled via capacitively coupled transmission lines (see the Supplementary Materials).

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 Embedded Image. The inductive energy of the shuntEmbedded Imagethen immediately yields an interaction term between all KPOs. Note that we use units such that the flux quantum Embedded Image. After including the capacitive energies and proceeding with standard circuit quantization (40), Embedded Image, where Zn 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)Embedded Image(3)

Here, Embedded Image is the natural frequency of the nth KPO expressed in terms of the charging and Josephson energies (40) of the Josephson junctions Embedded Image and Embedded Image and the dc part of the flux bias Embedded Image. The two-photon drive strength is Embedded Image, where Embedded Image and Embedded Image 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 Embedded Image, the strength of the Kerr nonlinearity of oscillator n is Embedded Image, and finally, the interaction strength between oscillators n and m is Embedded Image.

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 Embedded Image 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 Embedded Image, where Embedded Image denotes the Josephson energy of the shunt junction, the latter effectively acts as a negative inductor and the coupling matrix elements become Embedded Image (see the Supplementary Materials). We remark that a π-junction (4244) 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 Embedded Image 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 equationEmbedded Image(4)

HereEmbedded Imageand HN 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 asEmbedded Image(see the Supplementary Materials). If the drive strength is varied adiabatically, then under the action of HNH, 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)Embedded Image

Fig. 3 Comparison between coherent and dissipative quantum annealing for two antiferromagnetically coupled equal KPOs.

The quantum trajectory Embedded Image is obtained by numerically solving the stochastic Schrödinger equation with the Hamiltonian (3) and the photon loss rate κ. The fidelities with respect to the vacuum Embedded Image as well as the three states Embedded Image and Embedded Image, are shown. The amplitude α is given by Eq. 5. (A) Without dissipation: κ = 0. The system evolves from the vacuum Embedded Image (dashed blue line) at t = 0 to the even parity cat state Embedded Image (full red line) at t = T. The latter state encodes the ground state of the corresponding antiferromagnetic Ising model (J < 0). The population of the odd photon number parity state Embedded Image remains 0 (dashed purple line). The bifurcation dynamics is clearly visible as a kink of the population of Embedded Image, when the drive strength reaches the threshold value Embedded Image (vertical thin dashed black line). (B) With dissipation: κ = 0.01 MHz. A quantum trajectory with six jumps obtained from a Monte Carlo simulation of the dissipative dynamics is shown. A photon loss event induces a transition between the even and odd photon number parity cat states. However, note that both Embedded Image and Embedded Image correctly encode the antiferromagnetic Ising spin correlations. Also, note the absence of jumps below threshold, where the average photon number (thin magenta line) is close to 0. The parameter values used in both simulations are as follows: Δ = −1 MHz, J = −0.5 MHz, K = 0.7 MHz, T = 400 μs, εMAX = 2.0 MHz, and ε(t) = εMAX(t/T).

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

Let us now turn to the jump part of the dissipative process (see Eq. 4). Below threshold, the average photon number Embedded Image 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, (Embedded Image, as well as Embedded Image. 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 Embedded Image. To encode the corresponding Ising problem into our circuit, we set the coupling matrix elements as Jij = − J0AiAj and Jii = 0. Here, J0 > 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 J0 = 1/42. The two Ising spin configurations that satisfy this NPP are {s1 = s4 = 1, s2 = s3 = −1}, and the configuration with all spins flipped. Figure 4A shows the success probability as a function of (dε/dt)−1 = TMAX 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 Embedded Image, is equal to sisj, 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 Embedded Image for i ∈ {1,2,3,4}, a dephasing error described by Embedded Image 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).

Fig. 4 Comparison of performance between the CVIM (A and B) and a standard discrete qubit-based QA (C and D).

(A and C) Success probability of the NPP with set A = {4,5,6,7} as a function of the inverse ramp rate (dε/dt)− 1 and the photon loss rate κ in (A) or dephasing rate γ in (C). (B and D) Corresponding mean number of jump events [photon loss for (B) and dephasing events for (D)]. Although the success probability of the QA drops sharply already after a single (on average) dephasing event [see region within the white contour in (D)], the CVIM still succeeds with probability > 0.5 even when more than one photon has been lost [see region within the white contour in (B)]. Also, the adiabatic ramp rate threshold for the CVIM is substantially higher than that for the QA. Finally, the success probability for the CVIM is typically above that for a random guess, ~ 1/8 = 0.125, in the entire region shown above the adiabatic threshold, whereas the success probability of the QA quickly drops below the random guess value delimited by thin solid black contour lines in (A) and (C). The parameter values for the CVIM simulations in (A) and (B) are as follows: Δ = −1.5 MHz, K = 0.6 MHz, εMAX = 2 MHz, and ε(t) = εMAX(t/T) for T in the range (0, 200) μs. The parameters for the QA simulations in (C) and (D) are as follows: ε(t) = εMAX(1 − t/T) with εMAX = 6 MHz and T in the range (0, 600) μs. The dashed black curves in (A) and (C) indicate the points where κT = 1 and γT = 1. Each point in all figures corresponds to an average over 40 trajectories.


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).


Two coupled KPOs

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

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 readEmbedded Image

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, isEmbedded Image

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

When Embedded Image, 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 material for this article is available at

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 Embedded Image.

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

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

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

References (5158)

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.


Acknowledgments: The numerical calculations were performed in a parallel computing environment at sciCORE ( 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.

Stay Connected to Science Advances

Navigate This Article