High-accuracy Ising machine using Kerr-nonlinear parametric oscillators with local four-body interactions

A two-dimensional array of Kerr-nonlinear parametric oscillators (KPOs) with local four-body interactions is a promising candidate for realizing an Ising machine with all-to-all spin couplings, based on adiabatic quantum computation in the Lechner-Hauke-Zoller (LHZ) scheme. However its performance has been evaluated for only a few KPOs. By numerically simulating more KPOs, here we show that the performance can be dramatically improved by reducing inhomogeneity in photon numbers induced by the four-body interactions. The discrepancies of the photon numbers can be corrected by tuning the detunings of KPOs depending on their positions, without monitoring their states during adiabatic time evolution. This correction can be used independent of the number of KPOs in the LHZ scheme and thus can be applied to large-scale implementation.


INTRODUCTION
Hardware devices for solving the Ising problem [1] have attracted attention because this problem can represent various combinatorial optimization problems [2]. Such Ising machines are thus expected to be applied to realworld problems, such as integrated circuit design [3], computational biology [4,5], and financial portfolio management [6]. The Ising machines have been developed in several implementations including quantum annealing [7,8] or adiabatic quantum computation [9][10][11] with superconducting circuits [12], and coherent Ising machines with laser pulses [13][14][15][16]. In these devices, quantum effects such as quantum tunneling and quantum fluctuations are expected to be exploited for solving the Ising problem. Classical Ising machines have also been implemented with digital circuits based on simulated annealing [17][18][19] and simulated bifurcation [20][21][22].
Another adiabatic quantum computation for the Ising problem has been proposed using Kerr-nonlinear parametric oscillators (KPOs) [23,24]. A KPO is a parametrically driven oscillator with Kerr nonlinearity, and exhibits a bifurcation [25]. The bifurcation allows KPOs to represent binary Ising spins. Furthermore, a KPO can generate non-classical states such as a quantum superposition of coherent states, known as a Schrödinger cat state, via a quantum adiabatic bifurcation [23,24]. This Ising machine has been referred to as a quantum bifurcation machine (QbM) [24,26]. As driven states in KPOs are used, QbMs can contrast with quantum annealers operated in thermal equilibrium [12,27].
QbMs for the Ising problem with all-to-all spin couplings, which allow us to solve a wider class of real-world problems, have been proposed [38][39][40]. The architecture in Ref. [39] is based on the Lechner-Hauke-Zoller (LHZ) scheme [41][42][43][44][45][46], where the all-to-all spin couplings can be realized by a two-dimensional array of KPOs with local four-body interactions. This architecture is a promising candidate for large-scale implementation of the LHZ scheme, because in the case of QbM, the four-body interaction is realized by four-wave mixing in a single Josephson junction [39]. This simple four-body interaction can be an advantage over quantum annealers in the LHZ scheme with flux qubits [47] or transmon qubits [48], where a four-body interaction may need complicated circuits with multiple ancilla qubits.
In the previous numerical studies, however, to our knowledge the number of KPOs in the LHZ scheme has been limited up to three [49], and the performance of its large-scale implementations has not been revealed. In this paper, we investigate adiabatic quantum computation using a KPO network in the LHZ scheme (LHZ-QbM) with a larger number of KPOs. By a variational method, we first find that mean photon numbers in the KPOs become inhomogeneous, i.e. unequal, owing to asymmetry in the LHZ scheme. This inhomogeneity can degrade accuracy in solving the Ising problem. We then propose a method to reduce the inhomogeneity by scheduling the detunings of KPOs without monitoring their quantum states during adiabatic time evolution. Finally we numerically demonstrate that this method can improve the solution accuracy dramatically. This method is applicable to arbitrary numbers of KPOs in the LHZ-QbM and therefore to its large-scale implementation.

LHZ-QbM
We first briefly explain the Ising problem and the LHZ scheme, and then introduce an LHZ-QbM. In the Ising problem [1], a dimensionless Ising energy, In the LHZ scheme [41], the product of two Ising spins s i s j is mapped to a variables k = ±1, which we call an LHZ spin. This mapping reduces the two-body interaction J ij to a local external field J k , while increasing the number of the spins to L = N (N − 1)/2. The additional degrees of freedom are removed by four-body constraints on neighboring LHZ spins:s kslsmsn = 1. In Fig. 1, which shows the lattice structure of the LHZ scheme for N = 4 (L = 6), a four-body constraint is imposed on four LHZ spins connected by a square. The lower edge of the lattice is terminated by ancillary LHZ spins fixed to 1.
To satisfy the constraints, terms proportional tõ s kslsmsn are added to the Ising energy, resulting in the LHZ energy, E LHZ = − L k=1 J ksk − C k,l,m,n s kslsmsn , to be minimized. Here the summation in the second term is over all the constraints. The second term means four-body interactions with their strength C. With sufficiently large C, the four-body constraints are satisfied by {s k } minimizing E LHZ (see Appendix A).
We search for the solution of the Ising problem by embedding the LHZ energy into a network of L KPOs described by a Hamiltonian H = H KPO + H LHZ , where a k and a † k are, respectively, the annihilation and creation operators for the kth KPO, and , K, p, and ∆ k are the reduced Planck constant, the Kerr coefficient, the pump amplitude, and the detuning frequency for the kth KPO, respectively. This expression is obtained in a frame rotating at half the pump frequency, ω p /2, of the parametric drive and in the rotating-wave approximation [23,25,39]. In this work, we assume that ∆ k can be controlled individually.
H LHZ corresponds to the LHZ energy. In Eq. (2), we introduce parameters, ξ and A, to scale H LHZ and the term including {J k }, respectively, where ξ has the dimension of frequency and A is dimensionless. H LHZ is small compared with H KPO , ξ K, for accurate embedding [23,39]. A is chosen to reproduce E LHZ (see Methods). The first and second terms in H LHZ physically mean local external drives at ω p /2 [26,28,39] and four-body interactions, respectively. Since the ancillary LHZ spins on the lower edge of the lattice are fixed to 1, as mentioned above, the corresponding KPOs can be replaced by classical driving [39]. Thus their annihilation and creation operators in H LHZ are replaced by classical driving amplitudes denoted by β (see Appendix B). A solution of the Ising problem is found from the ground state of the LHZ-QbM after adiabatic time evolution via bifurcations [23,24,39]. The sign of a quadrature amplitude x k = a k + a † k /2 provides the LHZ spiñ s k [23,39], and the LHZ spins are transformed to Ising spins, giving the solution to the Ising problem. Since the LHZ spins represent the relative signs of pairs of the Ising spins, a configuration {s k } corresponds to both {s i } and {−s i } simultaneously, which give an identical Ising energy. Thus in mapping the LHZ spins to the Ising spins, we can fix one of the Ising spins, e.g., s 1 = 1. The probabilities for the LHZ spins, hence for the Ising spins, can be evaluated from the wave function of the ground state.
In simulating a larger number of KPOs than before, we found that computational costs can be larger for the computation of the probabilities for the spins than that of time evolution. These large costs come from the continuous degrees of freedom of the KPOs, which are finally projected to the binary variables. We succeeded in calculating these probabilities within a reasonable computational time by noting that the Ising spins can be determined from only (N − 1) KPOs [41] and thus tracing out the others (see Methods). For instance, only the KPOs 1, 2, and 3 are readout in the case where N = 4, as shown in Fig. 1.
Before showing simulation results, we evaluate the ground state by an analytical method, in which we find inhomogeneity in the KPOs. Since this inhomogeneity will degrade the accuracy of the solutions, we propose a method to correct the inhomogeneity.

Inhomogeneous photon numbers in the LHZ-QbM
We examine the mean photon numbers in the ground state of H by using a variational method. Here, we assume a small ξ satisfying ξ/K 1. Although the terms yielded by the four-body interactions in Eq. (2) seem intractable, these terms are simplified when the four-body constraints hold (see Methods). As a result, the mean photon number in the kth KPO can be approximated by where α 0 = p −∆ /K 1/2 , and∆ is the average of ∆ k over the KPOs. z k is the number of the four-body interactions connected to the kth KPO. The first term  (3), the terms proportional to J ksk and Cz k are caused by the local external drives and the four-body interactions, respectively. Because the term proportional to Cz k typically becomes dominant under the condition for satisfying the four-body constraints, we focus on this term in the following. In general, z k takes the values of 1, 2, 3, and 4 as can be seen from Fig. 1. Hence the variation in the mean photon number due to the fourbody interactions can become smaller for KPOs near to the edge of the network while larger around its center.
This inhomogeneity due to the four-body interactions effectively changes the local drive in each KPO and consequently prevents the local external drives in Eq. (2) from accurately implementing {J k }, which results in low solution accuracy. The solution accuracy is thus expected to be improved by reducing the inhomogeneity in the photon numbers. Here we propose a method to reduce the inhomogeneity without measurement.

Correcting the inhomogeneous photon numbers
We reduce the inhomogeneity in the mean photon numbers by setting ∆ k as wherez is the average of z k . Thus, within the first order approximation in ξ/K, α 2 0 in Eq. (3) can be replaced by (p − ∆)/K. Substituting Eq. (4) into the first term in Eq. (3), we obtain homogeneous mean photon numbers: where the second term in Eq. (4) has canceled the term proportional to C in Eq. (3), and the term proportional to J ksk has been dropped.
The correction by Eq. (4) is possible without knowing the solutions to the Ising problems, because the term proportional to C in Eq. (3) is independent of {J k } or {s k }. The parameters are therefore scheduled in advance, and no measurement of the states is necessary during the time evolution.
In the following, we simulate the time evolution of the LHZ-QbM shown in Fig. 1 by numerically solving the Schrödinger equation with H. (See Methods for details.) N = 4, hence L = 6, is chosen because this network is the smallest among those where the KPOs are unequal. We first apply this LHZ-QbM to the Ising problem with uniform interactions to check the inhomogeneity, and then with random interactions in order to evaluate the performance.

Simulation results for uniform antiferromagnetic interactions
Here we solve the Ising problem with uniform all-to-all connected antiferromagnetic interactions [14], which is expressed by the same negative J for all {J k }. In this work, we normalize the {J k } such that L k=1 |J k | = 1.    Figure 2a shows the mean photon number in each KPO after time evolution, where we compare the results for uniform ∆ (without the correction) and modified ∆ k in Eq. (4) (with the correction). For uniform ∆, the mean photon numbers increase depending on z k (z 1 = z 3 = z 6 = 1, z 4 = z 5 = 2, and z 2 = 3 as can be seen from Fig. 1), which are consistent with Eq. (3). Since {J k } is uniform, this inhomogeneity originates from the four-body interactions.
For modified ∆ k , on the other hand, all the mean photon numbers are nearly equal to 3, which coincides with the value predicted by Eq. (5). This result shows that the photon numbers can be made almost homogeneous by using ∆ k in Eq. (4) as expected. Figure 2b shows the probabilities for the Ising spin configurations. Although the probabilities are finite only for the three degenerate ground states, namely two spins in 1 and the others −1, the probabilities for these three configurations are not equal. In particular, uniform ∆ results in the much higher probability for {+ + −−} than the other configurations. The reason for this can be explained as follows. The configuration {+ + −−} is mapped to the LHZ spins ofs 1 =s 3 = 1 and the others in −1. While the negative J k favorss k = −1, the four-body constraints lead tos 1 =s 3 = 1. Here J 1 and J 3 become effectively small because of the relatively small mean photon numbers in KPOs 1 and 3, as shown in Fig. 2a, which is caused by the four-body interactions. Thus the probability for this configuration becomes particularly high compared to the others.
When ∆ k in Eq. (4) are used, the probabilities are similar for two of the ground states, {+ − −+} and {+ + −−}, as shown in Fig. 2b. Thus, the bias to one ground state {+ + −−} observed for uniform ∆ are suppressed by reducing the inhomogeneity in the photon numbers.

Simulation results for random interactions
Finally we evaluate the performance of the LHZ-QbM for the Ising problem with all-to-all connected random interactions, generating {J ij } uniformly from {−1, −0.99, −0.98, · · · , 1} for 100 instances [23,24,26], and normalizing them such that L k=1 |J k | = 1, as mentioned above. The performance is measured by a success probability P s and a residual energy E res . The success probability is defined as a probability for obtaining {s i } minimizing E Ising . The residual energy is given by the expectation value of E Ising subtracted by the minimum E Ising [23,24,50]. A lower residual energy indicates higher accuracy, and takes its minimum of 0 when the success probability is 1. Figures 3a and 3b show the distributions of P s and E res without and with the correction, respectively. In each case, the values of (C, ξ/K) are set to (0.3, 0.3) and (0.4, 0.6), in order to maximize the success probabilities averaged over the instances (see Appendix D). Without the correction, for several instances the success probabilities are nearly 0. With the correction, in contrast, the success probabilities are at least higher than 0.3 and close to 1 for most instances. With the nonzero success probability, repeated use of this LHZ-QbM gives a correct solution with probability rapidly approaching 1. The residual energies are also substantially lowered by the correction, which means that obtained solutions become more accurate.
In Fig. 3c we show improvement rates of a failure probability P f = 1 − P s and E res in each instance, where an improvement rate is defined as the ratio of the value without the correction to that with it. Figure 3c shows that both P f and E res are improved by one or two orders of magnitude in many instances, with the factors of up to P

DISCUSSION
In the LHZ-QbM, i.e. a KPO network for adiabatic quantum computation in the LHZ scheme, we have shown that inhomogeneous photon numbers due to four-body interactions degrade its performance, and furthermore that the inhomogeneity can be suppressed by controlling detunings, which can dramatically improve the performance. This method does not need to refer to the states of the KPOs during adiabatic time evolution, offering a simple operation. This method can be used regardless of the number of KPOs in an LHZ-QbM, thus allowing its large-scale implementation.
While we have modified only the detunings in the present work, similar corrections are possible by setting Kerr coefficients or pump amplitudes. In implementations with superconducting circuits, these three parameters can be controlled experimentally through parameters characterizing these circuits [24,39].
In the present study, we have assumed sufficiently small H LHZ compared to H KPO , and based on this assumption, we have determined the detunings to correct the inhomogeneity. From the simulation results, however, we have found that the optimal value, ξ/K = 0.6, is rather large (see Appendix D). This large optimal ξ/K may be because the time evolution becomes more adiabatic, that is, the larger H LHZ widens gaps between the energy levels of ground states and excited states, and prevents transitions between these states during the time evolution [11]. These results imply further possible improvement by increasing ξ/K together with more elaborate control of detunings than the present analytic ones.

Simulation of the LHZ-QbM
Each KPO is initialized to a vacuum, and the parameters in H are set such that the vacuum is its ground state. For the initial parameters p/K = A = β = 0 and ∆/K > 0, the vacuum is the ground state if ξC/K ≤ 1, which is valid within an estimation by a variational method [23,26] (see Appendix C). Time evolution is then started, and the parameters are slowly varied to induce bifurcations. During the evolution, the state is expected to remain in the ground state. Figure 4 shows the parameters as a function of time [23,26], where the time interval for the evolution is T = 500/K. To reproduce E LHZ , A and β for t T are chosen such that A α 3 and β α, where α = [(p − ∆)/K] 1/2 .
The time evolution is simulated by numerically solving the Schrödinger equation, where H and the state |ψ are represented in the photon number basis with the largest The others, K, ξ, and C, are set to be constant [24,39], and positive K is assumed [23,24]. (For negative K, the same results can be obtained by changing the signs of p, ∆, and ξ [23,24]. Negative K has been used in Ref. [39].) photon number 12 truncating the Hilbert space [23,26]. After the time evolution, the probabilities for the LHZ spins can be formulated using eigenstates of x k as [23,26,51] P (s 1 , · · · ,s L ) = Tr |ψ ψ| However we found that their calculation becomes the largest in the simulation for large L. Then, noting that (N − 1) LHZ spins can determine the Ising spins [41], we trace out the other KPOs in Eq. (6) without corresponding M k (s k ) and obtain the probabilities within reasonable computational costs.

Inhomogeneity in the KPOs in the LHZ-QbM
We evaluate the ground state using a variational method and assuming the bifurcations, small H LHZ , and the fourbody constraints. We employ |α 1 · · ·|α L as a trial wave function [23,26], and minimize the expectation value H with respect to {α k , α * k }, where α * k is the complex conjugate of α k . The following nonlinear equations are then yielded By assuming sufficiently small ξ, Eq. (8) is approximately solved within the first order in ξ as where α 0k = [(p − ∆ k ) /K] 1/2 , and α 0ksk is the solution of Eq. (8) for ξ = 0 after the bifurcation p − ∆ k > 0.
Here ∆ k −∆ is assumed to be the first order in ξ [which is valid in Eq. (4)].
The term proportional to C in Eq. (9) can be simplified under the four-body constraintss kslsmsn = 1 as follows. Sinces lsmsn =s k holds, the summation in Eq. (9) reduces to the number of connected four-body interactions z k , resulting in This expression becomes accurate for ξ/K 2/(|J k | + Cz k ) ∼ 1. Equation (10) shows that H LHZ affects the amplitude |α k |. Its square |α k | 2 gives the mean photon number in Eq. (3) in this approximation.

This work was supported by JST ERATO (Grant No. JPMJER1601).
Appendix A: Condition on C for satisfying the four-body constraints We derive a formula that expresses the condition on C for satisfying the four-body constraints for {s k } minimizing E LHZ . We then show that C > 1 is a sufficient condition under the normalization condition LHZ , is expressed as is the spin configuration minimizing E LHZ with b broken constraints. By using these notations, the condition for satisfying the constraints is written as E (0) LHZ for all b ≥ 1. Thus, using Eq. (A1), we have the following inequality: Note that the right-hand side in Eq. (A2) offers a tight lower bound for C to satisfy the constraints for an arbitrary instance {J k }. In the case of the uniform antiferromagnetic interaction in the main text (L = 6 and J k = −1/6), the values of − k for b = 0, 1, 2, 3 are −1/3, −2/3, −1, and −2/3, respectively, and thus Eq. (A2) gives C > 1/6. Next, we roughly evaluate a condition on C valid for any instance {J k }. The right-hand side of Eq. (A2) can be estimated as where we have used b ≥ 1, s ≤ 2, and the normalization condition L k=1 |J k | = 1. Equation (A3) shows that the right-hand side of Eq. (A2) is at most 1, indicating that C > 1 is a sufficient condition for satisfying the four-body constraints.  where the products are over k = 2, 4, 5, 6. If K ≥ ξC, the quantity in the square bracket is always positive and thus α 2 α 4 α 5 α 6 must be zero. α 2 α 4 α 5 α 6 = 0 and Eqs. (C3)-(C6) result in α 2 = α 4 = α 5 = α 6 = 0. Around α k = 0, H is approximated by H 6 k=1 ∆ |α k | 2 , which shows that α k = 0 corresponds to the minimum. This result means that if ξC/K ≤ 1 then the initial ground state is the vacuum within the estimation by the variational method. Figures 5a (without the correction) and 5b (with the correction) show the success probabilities averaged over the 100 random instances in the main text as functions of C and ξ. By the proposed correction, the highest average success probability is increased from 89.5% at (C, ξ/K) = (0.3, 0.3) to 97.1% at (0.4, 0.6), which are used for Fig. 3. The average residual energies are shown in Figs. 5c (without the correction) and 5d (with the correction), which are decreased by one order of magnitude by the correction.
The dependences of the average success probabilities on C are understood as follows. (Those of the residual energies can be understood in the same way.) We first explain the case without the correction shown in Fig. 5a. As a function of C, the average success probability shows a maximum. The low success probabilities for C 1 are because the four-body constraints are frequently broken. With increasing C, the inhomogeneity in the photon numbers due to the four-body interactions is increased and degrades the success probabilities. Thus the average success probabilities are the highest at an intermediate C.
For the case with the correction shown in Fig. 5b, the average success probability also shows a maximum as a function of C. For small C, the success probabilities are low, again, owing to violated four-body constraints as mentioned above. For intermediate C, the correction can suppress the effects of the inhomogeneity in the photon numbers, unlike the above case, leading to higher performance. For larger C, however, the success probabilities suddenly decrease, because the estimation given by Eq. (3) is no longer valid for these (C, ξ/K), and the resultant large modulation of detunings according to Eq. (4) rather degrades the performance.