Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network

The dynamics of nonlinear systems qualitatively change depending on their parameters, which is called bifurcation. A quantum-mechanical nonlinear oscillator can yield a quantum superposition of two oscillation states, known as a Schrödinger cat state, via quantum adiabatic evolution through its bifurcation point. Here we propose a quantum computer comprising such quantum nonlinear oscillators, instead of quantum bits, to solve hard combinatorial optimization problems. The nonlinear oscillator network finds optimal solutions via quantum adiabatic evolution, where nonlinear terms are increased slowly, in contrast to conventional adiabatic quantum computation or quantum annealing, where quantum fluctuation terms are decreased slowly. As a result of numerical simulations, it is concluded that quantum superposition and quantum fluctuation work effectively to find optimal solutions. It is also notable that the present computer is analogous to neural computers, which are also networks of nonlinear components. Thus, the present scheme will open new possibilities for quantum computation, nonlinear science, and artificial intelligence.

driving). This model is not only the simplest one for the present purpose but also physically feasible. Promising candidates for implementing this model are superconducting microwave resonators with Josephson junctions 3,10-13 , nanoelectromechanical systems 14,15 , and carbon nanotubes 16 (also see ref. 2).
In a frame rotating at half the pump frequency of the parametric drive and in the rotating-wave approximation, its Hamiltonian is given by 2 K a a  p a  a  2  2  1   1  2 2  2  2 where a and a † are the annihilation and creation operators for quanta of the oscillator (the quanta are, e.g., photons for electromagnetic resonators or phonons for mechanical oscillators), Δ is the detuning of the oscillator eigenfrequency from half the pump frequency, K is the Kerr coefficient for the Kerr effect, and p is the pump amplitude for the parametric drive. Hereafter, we assume that K and Δ are positive constants and p is a nonnegative control parameter. When K is negative, similar discussion is straightforward. A cat state is generated via quantum adiabatic evolution as follows. Initially, the state of the oscillator and the pump amplitude are set to the "vacuum" 0 and zero, respectively. Note that 0 is the ground state of H 1 with p = 0. As p is increased slowly, the system follows adiabatically the ground state of H 1 (t). When p becomes much larger than Δ , the final state becomes approximately the ground state of H 1 with Δ = 0, which is a superposition of two coherent states ± / p K . ( has been ignored assuming sufficiently large p. This result can be understood qualitatively from a viewpoint of classical nonlinear dynamics as follows. A classical model for the above system is defined by 2 2 where x and y are real variables corresponding to the Hermitian operators ( + )/ † a a 2 and ( − )/ † a a i 2 , respectively, and the dots denote differentiation with respect to time t. These equations are derived by replacing a with x + iy in the Heisenberg equations of motion for a with H 1 . This model exhibits a bifurcation at p = Δ as follows. When p ≤ Δ , the origin is a single fixed point, which is stable. (Fixed points are defined by = =   x y 0 1 .) When p > Δ , the origin becomes an unstable fixed point and two stable ones are created, the positions of which are (± ( − ∆)/ , ) p K 0 . The dependence of the fixed points on p is depicted in Fig. 1a,b as the bold lines, where the solid and broken lines correspond to the stable and unstable fixed points, respectively. (In Fig. 1, Δ is set to K.) Such figures are called bifurcation diagrams 1 . The oscillating thin curves in Fig. 1a,b are obtained by numerically solving Eqs (2) and (3), where p(t) is increased linearly from p(0) = 0 to ( / ) = p K K 500 5 and the initial condition is set as ( ) = .
x 0 01 and ( ) = y 0 0. It is found that the state changes along one of the stable branches. This result is explained by the adiabatic theorem in classical mechanics 18 (see Supplementary Information for details). Thus, the cat-state generation described above is interpreted as a superposition of two coherent states corresponding to the two classical stable branches. In other words, while the classical system chooses one of the two branches (which branch the system will choose may be unpredictable because of chaos), the quantum system can follow both the branches "simultaneously" as a superposition state. To emphasize this nonclassical feature of a quantum nonlinear oscillator, we refer to such an intriguing process as a quantum-mechanical bifurcation.
We performed numerical simulation of the cat-state generation. In this simulation, the Schrödinger equation with H 1 is numerically solved, where the Hilbert space is truncated at a "photon" number of 20, the initial state is set to 0 , and p is increased linearly from ( ) = p 0 0 to ( / ) = p K K 500 5 . Figure 1c,d show the Wigner function ( , ) W x y at = . ∆ p 0 9 and = ∆ p 5 , respectively. Here the Wigner function is a quasiprobability distribution for x and y 3,10,17 (also see Methods). In Fig. 1c On the other hand, the interference fringe between the two peaks in Fig. 1d means that the two oscillation states are superposed 17 , that is, a cat state is generated successfully. The negative values of ( , ) W x y in Fig. 1d indicate the nonclassicality of the cat state. To realize the cat-state generation, we require a large Kerr coefficient compared to a loss rate. While this requirement is too stringent for optical and mechanical systems, superconducting circuits with Josephson junctions have already achieved this situation 11,12 . Thus, superconducting systems are the most promising for implementing the present scheme.
Adiabatic quantum computation with a KPO network. Next, our quantum computer with KPOs is described. If we have N independent KPOs, we will obtain a superposition of 2 N oscillation states via the quantum-mechanical bifurcation described above. To exploit the exponentially large number of states for solving combinatorial optimization problems, we couple the KPOs to one another appropriately depending on given problems.
The combinatorial optimization problem studied here is the Ising problem: given a dimensionless Ising energy give the same Ising energy, and therefore there are always two solutions for each problem. This Ising problem is extremely hard unless the coupling topology is too simple; more precisely, it is known to be non-deterministic polynomial-time hard (NP-hard) in computational complexity theory 19 . Recently, machines specially designed for the Ising problem have attracted much attention [20][21][22][23][24][25][26][27] .
For the above problem, we couple N KPOs as follows: where ( ) H i 1 is the Hamiltonian for the ith KPO of the form of Eq. (1) with an individually controllable detuning Δ i and ξ 0 is a positive constant with the dimension of frequency. Note that the coupling Hamiltonian describes standard linear couplings, and therefore is physically feasible. It is also notable that H is symmetric under simultaneous parity inversion defined as → − a a i i for all i simultaneously. In the following, we show that the KPO network can solve the Ising problem via quantum adiabatic evolution. To use a quantum adiabatic evolution for finding a configuration minimizing E Ising , the initial state 0 should be the ground state of H with p = 0. This condition can be satisfied by setting the detunings such that the following matrix M becomes positive semidefinite (see Supplementary Information for  By increasing p slowly, we obtain the ground state of H with large p assuming that the so-called adiabatic condition 5,9,28 is satisfied.
When p becomes much larger than ∆ i and ξ | | , J i j 0 , the nonlinear terms in H are dominant, the ground state of which is 2 N -fold degenerate and the eigenspace is spanned by the tensor products of ± / p K . By the perturbation theory to the lowest order 28 , the correction to the energy of a tensor product where / − / p K p K has been ignored assuming sufficiently large p. Note that the first term in the right-hand side of Eq. (8) is independent of s { } i and the second one is proportional to the Ising energy. Consequently, the ground state is two-fold degenerate and the eigenspace is spanned by are the two solutions of the Ising problem. The degeneracy comes from the simultaneous parity symmetry of H. Taking the simultaneous parity symmetry of H into account, the ground state obtained as a final state in the above adiabatic evolution is given by Thus, it turns out that we can find a solution of the Ising problem by measuring the amplitudes of the KPOs and identifying their signs with the Ising spins.
The above degeneracy means that the energy gap between the ground state and the first excited state will vanish, which seems problematic for the adiabatic approach. However, this causes no problem for the following two reasons. First, the transition between the two states is prohibited by the simultaneous parity symmetry of H. (The ground and first excited states have even and odd parities, respectively.) Second, even if the transition occurs by some accidental errors, we can find a solution correctly because the first excited state is also a superposition of (the same as ψ | 〉 f except for a negative relative phase). In the above discussion, we have shown that the KPO network can solve the Ising problem by assuming the adiabatic condition. From the adiabatic condition, it turns out that the speed of adiabatic quantum computation is limited by the minimum energy gap between the ground and excited states 5,9 . Thus, the study of the energy level structure is important for evaluating and improving the computational speed [29][30][31] . In the present scheme, however, the energy level structure may be much more complex than the ones in conventional adiabatic quantum computation, because the fundamental component of our quantum computer is a nonlinear oscillator described by an infinite-dimensional Hilbert space, not a simple two-level system (qubit). In the present work, we instead performed numerical simulations, the results of which support the above discussion (see below).
It is also notable that an entangled cat state given by Eq. (9) is generated as a result of the quantum computation. We confirmed this fact by numerical simulation in the case of two spins (see Supplementary Information). Thus, the present scheme also provides a method for the generation of the intriguing states via quantum adiabatic evolution.
Numerical simulation results. Finally, we present numerical simulation results of the quantum computation for four-spin problems, which are more difficult than two-and three-spin ones in the sense that, in the four-spin case, there may be not only frustration but also a nonglobal local minimum. In these simulations, the Schrödinger equation with H in Eq. (5) is numerically solved, where the Hilbert space is truncated at a "photon" number of 18 for each KPO, the initial state is set to 0 , ξ = . K 0 5 0 , the detunings are set as Eq. (7), and p is increased linearly from ( ) = p 0 0 to ( / ) = p K K 700 7 . We generated 1000 instances of the problem with the coupling coefficients chosen randomly in the range −1 to 1. We estimated the success probability and the residual energy for each instance (see Methods for details), where the residual energy is defined as the difference between the Ising energy obtained by simulation and its minimum value 8 . The histograms of the success probabilities and the residual energies are shown in Fig. 2a,b, respectively. In Fig. 2, we treat two configurations s { } i and −s { } i as a pair because these give the same value of E Ising . We also simulated a classical model for the quantum computation, the equations of which are derived in a similar manner to Eqs (2) and (3) (see Supplementary Information for details). From the results for a single KPO, the comparison between the quantum and classical models may be helpful for understanding the simulation results. For each instance, we repeated the classical simulation 10 3 times, setting the initial values of x i and y i to random numbers in the range − 10 −6 to 10 −6 . The success probability and the residual energy for each instance are estimated by taking averages. The results are shown in Fig. 2c,d.

Discussion
Here we discuss the results shown in Fig. 2.
First, it is notable that the classical model can find optimal solutions with high probability. This result comes from the fact that the classical model can approximately solve the Ising problem (see Supplementary  Information). The high success probability for the classical model means that the approximation is fairly good.
(Thus, this model may provide a new approach to combinatorial optimization problems with nonlinear dynamics. Although the study of the classical model, which may exhibit chaotic behaviors, is interesting, this is beyond the scope of the present work.) Next, it is clear that the quantum model can achieve higher performance than the classical one. Since the differences between the two models may be mainly quantum superposition and quantum fluctuation, the high performance may come from these quantum effects. To examine this point, we look into one of the most difficult instances, for which the classical model almost always fails.
The time evolutions of the probabilities for the spin configurations in the quantum and classical models are shown in Fig. 3a,b, respectively. Figure 3c shows  In this problem, there are two local minima. The classical model is trapped around the nonglobal local minimum (Fig. 3b). On the other hand, in the quantum model, a superposition of the two local minima arises from quantum fluctuations, and finally the probability for the global minimum converges to unity via quantum adiabatic evolution (Fig. 3a). From this result, we conclude that quantum superposition and quantum fluctuation work effectively to find optimal solutions in the present computation. (Although a more detailed comparison between the two models is interesting, which may become a new theme in the field of quantum chaos, this is beyond the scope of the present work.) While most previous computational models for quantum computation use qubits as fundamental components 5,32-34 , the present model uses quantum nonlinear oscillators. Moreover, this harnesses continuous degrees of freedom of a nonlinear oscillator network, that is, its quantum-mechanical bifurcation based on quantum adiabatic evolution for quantum computation for the first time. It is also notable that the present model is analogous to a neural network 35 . Thus, the present scheme will open a new paradigm in the fields of quantum computation, nonlinear science, and artificial intelligence.