An Ising Hamiltonian solver based on coupled stochastic phase-transition nano-oscillators

Combinatorial optimization problems belong to the non-deterministic polynomial time (NP)-hard complexity class, and their computational requirements scale exponentially with problem size. They can be mapped into the problem of finding the ground state of an Ising model, which describes a physical system with converging dynamics. Various platforms, including optical, electronic and quantum approaches, have been explored to accelerate the ground-state search, but improvements in energy efficiencies and computational abilities are still required. Here we report an Ising solver based on a network of electrically coupled phase-transition nano-oscillators (PTNOs) that form a continuous-time dynamical system (CTDS). The bi-stable phases of the injection-locked PTNOs act as artificial Ising spins and the stable points of the CTDS act as the ground-state solution of the problem. We experimentally show that a prototype with eight PTNOs can solve an NP-hard MaxCut problem with high probability of success (96% for 600 annealing cycles). We also show via numerical simulations that our Ising Hamiltonian solver can solve MaxCut problems of 100 nodes with energy efficiency of 1.3 × 107 solutions per second per watt, offering advantages over other approaches including memristor-based Hopfield networks, quantum annealers and photonic Ising solvers. An Ising solver that is based on a network of electrically coupled phase-transition nano-oscillators, which provides a continuous-time dynamical system, can be used to efficiently solve a non-deterministic polynomial time (NP)-hard MaxCut problem.

networks of degenerate optical parametric oscillators 14,15 .However, qubit-based quantum annealers have a high cost of operation and complexity because of their cryogenic environment.Further, although an optical coherent Ising machine has shown competitive performance compared with a quantum annealer 16 , it requires a long fibre ring cavity for implementing Ising spins through temporal multiplexing, as well as an extremely fast and power-hungry field-programmable gate array (FPGA) for implementing coupling in a measure-and-feedback scheme 15 .Digital CMOS annealers [11][12][13] rely on an external source for random number generation to introduce stochasticity, but maintaining true randomness is still technologically challenging and requires substantial post-processing.
In this Article, we report an Ising Hamiltonian solver based on coupled electronic phase-transition nano-oscillators (PTNOs) that form a continuous-time dynamical system (CTDS).The bi-stable phases of the PTNOs emulate the Ising spins and coupling between the oscillators emulates the Ising interaction, allowing an optimization problem to be mapped into the PTNO network by carefully choosing the coupling matrix.The dynamical evolution of this physical system reaches an energy minimization point, which represents the solution of the optimization problem.Such dynamics of physical systems have been explored before for a wide range of applications, including understanding neural activity [17][18][19] , achieving robotic locomotion control 20,21 and solving optimization problems using discrete-time Hopfield networks 22,23 .However, our PTNO-based Ising solver has an inherently distributed nature and continuous-time dynamics, allowing for completely synchronous updates of all the Ising spins.This property drastically lowers the time duration of each annealing cycle to one oscillation time period, reducing the time-to-solution value.
We use a PTNO-based CTDS prototype with eight oscillators to solve a benchmark NP-hard MaxCut problem.By increasing the strength of the injection-locking signal, PTNO can progressively converge to a more optimal solution, similar to performing classical annealing.We obtain a success probability of 30% when no annealing cycles are implemented and a success probability of 96% when over 600 annealing cycles are implemented.We also show in simulations that our Ising Hamiltonian solver can be used for dense MaxCut problems of up to 100 nodes.We calculate a high energy efficiency of 1.3 × 10 7 solutions per second per watt for MaxCut problems of 100 nodes, which translates to a fivefold improvement over the recently demonstrated memristor-based Hopfield network 23 and improvement of several orders of magnitude over other approaches, including central processing units (CPUs), graphics processing units (GPUs), quantum annealers and photonic Ising solvers.

PTNo-based cTDS as an Ising Hamiltonian solver
The combinatorial optimization problem is reformulated in terms of the Ising Hamiltonian H, which is defined by spin vector σ and symmetric coupling matrix J, and mapped onto the Ising solver, as shown in Fig. 1a,b.The Ising Hamiltonian with N discrete spins σ 1≤i≤N ∈ {−1, +1} N and symmetric coupling matrix J in the absence of an external magnetic field is given by J ij σ i σ j , where i and j are two spins.Each Ising spin represents a node in the graph and is emulated by an insulator-to-metal (IMT) phase-transition nano-oscillator.In our physical system, PTNO comprises a two-terminal phase-transition hysteretic device, based on vanadium dioxide (VO 2 ) as a prototypical IMT material 24 , in series with a transistor (Fig. 1c).The working principle of a VO 2 -based PTNO has been reported elsewhere 20,25 .Below the phase-transition temperature and under zero external electric field or current, VO 2 shows insulating behaviour.On application of an electric field across the two terminals of the device that forces current to flow through, the material undergoes an abrupt phase transition from the insulating to metallic state.The hysteretic phase transition is reflected as an abrupt hysteretic current-voltage (I-V) characteristic of the device.Pairing a conventional n-type metal-oxide-semiconductor (NMOS) transistor in series with VO 2 such that the load line passes through the unstable region of the I-V curve results in self-sustained oscillations, as shown in Fig. 1c.We used a supply voltage of V DD = 2.0 V and gate voltage V GS = 0.8 V in our experiments.We create a highly interconnected PTNO-based network where coupling matrix W is derived from adjacency or coupling matrix J of the Ising model, as shown in Fig. 1a,b.In this example of the Ising model with antiferromagnetic interactions, the presence of an edge between two nodes is denoted as J = −1 and is represented by coupling capacitance C C .
Figure 1d shows an overview of the experimental setup of the PTNO-based CTDS.The main computing kernel comprises eight PTNOs connected using coupling capacitance following coupling matrix W. The fabricated VO 2 device array (labelled as 1-8 in Fig. 1d) is connected with eight NMOS transistors in series to create eight PTNOs.To emulate artificial Ising spin σ, an external injection-locking signal S inj is applied to all the oscillators using injection capacitance C inj .As explained later, this creates bi-stable oscillator phase θ that can be used as state variables for computing.The continuous-time dynamics of the network are dictated by the internal energy of the system, or the Lyapunov function E(θ), which closely resembles the Ising Hamiltonian H.The continuous-time dynamics of the PTNO network and evolution of oscillator phases θ(t) with time t will then be determined by the inherent energy minimization property of the network and can be described by dθ(t)/ dt = −∂E(θ)/∂θ.The dynamics of the network continuously evolves in time so as to naturally minimize E(θ) and in the process minimizes the Ising Hamiltonian to reach the ground state, as shown in Fig. 1e.Once the network reaches an energy minima state, the output of the network is read out in terms of state variable θ and subsequently reformulated to provide the optimal solution of the original optimization problem.One formidable challenge in solving the combinatorial optimization problem is the inevitable increase in the complexity of the energy landscape with problem size.Specifically, the presence of a large number of local minima degrades the probability of reaching a global optimum.Stochasticity is a well-known technique utilized in Boltzmann machines and simulated annealing, to overcome the issue of getting trapped in local minima.In the latter case, a stochastic noise is used to perturb the state of the system and the magnitude of noise is slowly reduced over time (replicated as a change in the temperature parameter in the algorithm) as the system approaches the global optimum.In our system, we exploit the inherent stochasticity present in the IMT material 26 to escape the local minima as well as to introduce a novel way of gradually reducing the temporal fluctuations in the oscillator phases by increasing the strength of injection-locking signal S inj .This allows us to perform a classical annealing operation and obtain progressively better solutions, as schematically illustrated in Fig. 1e.

Artificial Ising spin
The binary degrees of freedom in the phase space of PTNOs arises from the second-harmonic injection-locking phenomenon.We apply a sinusoidal injection-locking signal S inj = V inj sin(2πf inj t) with f inj = 2f 0 as an external input to the oscillators across capacitor C inj = 20 pF, as shown in Fig. 1d.When V inj = 0, the oscillator is free running.Figure 2a(i) shows the voltage output waveform (V out ) measured over multiple runs.The corresponding phase of the oscillator, measured with respect to a reference sinusoidal signal of the same frequency f 0 , shows a constantly varying phase with a uniform probability distribution over the entire phase space, as shown in Fig. 2b(i).In contrast, when the oscillator is perturbed with S inj at the first harmonic, that is, f inj ≈ f 0 (also referred to as first-harmonic injection locking (FHIL)), V out shows a constant 80° phase-locking configuration with S inj , as shown in Fig. 2a(ii).The probability distribution of the phase, measured over multiple runs, shows a single Gaussian distribution, as shown in Fig. 2b(ii).Interestingly, when f inj ≈ 2f 0 , the oscillator waveform shows both in-phase (40°) and out-of-phase (220°) injection-locking configuration when measured over multiple runs, as shown in Fig. 2a(iii).The corresponding probability distribution shown in Fig. 2b(iii) exhibits a double Gaussian distribution highlighting an equiprobable and bi-stable phase portrait.This bi-stability provides the ideal means to encode the Ising spin in the electrical domain, where a phase of 40° represents up-spin, that is, σ = +1, and a phase of 220° represents down-spin, that is, σ = −1.
The continuous-time dynamics of phase difference θ between the oscillator and injection-locking signal can be described by a generalized version of Adler's equation (Gen-Adler) given by where n H is the nth harmonic of the IMT nano-oscillator, f 0 is the natural frequency of the oscillator, ϑ is the integration variable and The first term describes the frequency mismatch between V out and S inj , which contributes to phase slipping.The second term depends on the phase delay incurred due to the perturbation caused by S inj and is described in terms of the perturbation projection vector (PPV), ξ (Supplementary Notes 1 and 3 provide the details).The corresponding 'energy' function or Lyapunov function of the oscillator is given by where ϕ and ϑ are the integration variables.
The first energy term comes from the contribution of the frequency mismatch that creates an overall bias in the energy land-scape.However, since it is a linear term, it does not introduce any additional valleys or peaks in the energy landscape.The second term describes the interaction between the injection-locking signal and the oscillator.The corresponding probability distribution of the oscillator's phase can be calculated as P (θ i ) = e −E(θ i )/η Z , where Z = ∑ i e −E(θ i )/η is the partition function and η is analogous to the kT term in the Boltzmann distribution (k is the Boltzmann constant and T is temperature) and can be interpreted as a measure of the stochastic noise in the IMT oscillator.Expectedly, with zero injection locking, the 'energy' function of the oscillator stays flat, as shown in Fig. 2b, indicating a uniform distribution over the phase space as is experimentally obtained in Fig. 2a.For FHIL, as long as the relative frequency difference (f inj -f o ) is small compared with the second term in equation (2), dθ(t) dt = 0 exhibits one stable point at θ ≅ 0.4π, that is, one injection-locked equilibrium phase and the calculated E(θ) (assuming zero frequency mismatch) displays a single energy minimum.This results in a single Gaussian peak in the probability distribution of the oscillator's phase, as verified experimentally.In the case of second-harmonic injection locking (SHIL), dθ(t) dt exhibits two stable points at θ ≅ 0.2π and θ ≅ 1.2π, that is, the two equilibrium phases.The calculated E(θ) (assuming zero frequency mismatch) evolves into a double-well energy landscape that results in a double Gaussian distribution in the phase space, as faithfully reproduced in the measurements (Supplementary Notes 1 and 3 provide further details).The steady-state analysis of the Gen-Adler equation in the case of SHIL also predicts that in the presence of a frequency mismatch between the oscillator and the injection-locking signal, there exists a critical amplitude of injection signal V inj below which no stable solution exists.In our experiments, this critical amplitude is close to 1 V for a frequency mismatch of 0.1%, as observed in the measurements.Above this threshold, the bi-stability in the phase space begins to appear.Figure 2c shows the measured oscillator phase as a function of the amplitude of V inj .For very low V inj close to the critical value, the phase of the oscillator measured over multiple runs shows a wide distribution in the phase space, with the distribution narrowing as V inj increases.This behaviour can be understood by considering the perturbation in the energy landscape, as shown in Fig. 2e, for V inj = 1, 3 and 5 V.For V inj = 1 V, energy barrier E B separating the two stable equilibrium phases is low (E B ≈ 0.006 calculated from equation ( 2)).Hence, in the presence of stochastic noise, the oscillator's phase constantly fluctuates between the two stable phases, as shown in Fig. 2d.An increase in V inj to 3 and 5 V 0 increases barrier height E B to 0.012 and 0.020, respectively.This reduces the fluctuations in the measured oscillator's phase.This is reflected in the experimentally measured mean time between each phase flip, referred to as the dwell time τ dwell (analogous to Néel relaxation time for magnetization) as a function of V inj (Fig. 2f).The increase in τ dwell with increasing V inj accurately follows the Arrhenius's relation η , where τ 0 = 1 f 0 is the characteristic or attempt time (equal to the time period of the oscillator), α is the fitting parameter and η is the stochastic noise in the IMT oscillator.This characteristic of reduction in temporal fluctuations of the oscillator's phase with increasing amplitude of the injection-locking signal proves to be the key towards performing classical annealing in our hardware.

Replicating the interaction term in the Ising Hamiltonian
To implement a PTNO-based Ising solver that can replicate an artificial Ising spin system, the oscillators need to be connected to one another using coupling elements that emulate the ferromagnetic and antiferromagnetic nature of interaction.We first study the nature of interaction in a pairwise coupled oscillator system, as schematically shown in Fig. 3a, using capacitance C C as the coupling element.Figure 3b shows the experimentally measured phase distribution of the two oscillators using injection-locking capacitance C inj = 20 pF and coupling capacitance C C = 56 pF.The oscillators remain out of phase with each other and the two configurations, namely, (40°, 220°) and (220°, 40°), are two equally probable states.We measured the probability of the out-of-phase configuration for varying coupling strengths, as shown in Fig. 3c.The probability remains close to 0.5 for low coupling capacitance, implying that both in-phase and out-of-phase configurations are equally probable; the probability increases close to 1.0 for stronger coupling.We also compare our experimental results with experimentally calibrated PPV-based numerical simulations (Methods), as shown in Fig. 3c, showing excellent agreement.To understand the exact nature of capacitive coupling, we compare it with a two-spin Ising model with antiferromagnetic interaction (negative J), where the individual spins prefer to remain antiparallel with one another.In such a system, the probability of one of the possible configurations, namely, (↑↑, ↓↓, ↑↓, ↓↑), is determined by the Boltzmann distribution P (σ 1 , σ 2 ) = e −H(σ 1 ,σ 2 )/kT /Z , where σ 1 , σ 2 are the two spins, T is the temperature and e −H(σ′ 1 ,σ′ 2 )/kT is the partition function.With an antiferromagnetic interaction (negative J), the probability for antiparallel configuration (↑↓, ↓↑) increases on varying the interaction strength from J = 0 to J = −2.Thus, the antiferromagnetic interaction in an Ising Hamiltonian can be faithfully replicated using capacitive coupling in this IMT nano-oscillator-based system.Mathematically, the continuous-time dynamics of such a PTNO network can be further described by extending equations ( 2) and (3) to incorporate an additional coupling term and is given by where I osc,j = C C,j dVosc,j dt .The additional third term describes the coupling interaction energy between the pairs of oscillators.The corresponding 'total energy' function or the global Lyapunov function of the PTNO-based CTDS is then given by We use equation ( 4) to calculate the two-dimensional energy landscape for pairwise capacitively coupled IMT nano-oscillators, as shown in Fig. 3b.The calculated energy landscape exhibits four stable points or attractor states-two degenerate global minima in the out-of-phase configuration and two degenerate local minima in the in-phase configuration.By increasing the strength of capacitive coupling, the attractor states for the out-of-phase configuration become more prominent.
A similar investigation is performed for pairwise resistively coupled oscillators, as schematically shown in Fig. 3d.Contrary to the previous case, the measured oscillator phases as well as numerical simulations reveal a higher probability to be in phase with each other in either (40°, 40°) or (220°, 220°) configuration, as shown in Fig. 3e, for a coupling resistance of R C = 40 kΩ.To establish the nature of resistive coupling, the probability of in-phase configuration is measured as a function different R C values and compared with a two-spin Ising model with ferromagnetic interaction (positive J).Note that a lower R C value represents a higher coupling strength and hence higher J.The increase in the probability of in-phase configuration with decreasing R C , that is, increasing coupling strength, agrees well with the theoretical Ising model that predicts an increase in probability of parallel configuration (↑↑, ↓↓) for increasing J, as in Fig. 3f.We also use equation ( 4) with I osc,j = − V osc,j R C,j to calculate the two-dimensional energy landscape, as shown in Fig. 3e, The PTNOs are connected following the same adjacency (or connectivity) matrix of the graph using coupling capacitances.b, Schematic of the annealing schedule used in the experiment.A sinusoidal injection-locking signal at twice the oscillator frequency f 0 is applied with a linearly increasing amplitude over annealing time T anneal that corresponds to a linear annealing schedule.This is followed by phase readout time T readout .Thus, the total computation time T comp = T anneal + T readout .Annealing time T anneal is varied from 0 (corresponding to no annealing) to 660 oscillation cycles.c, Evolution of the phases of the eight oscillators, settling to either the in-phase or out-of-phase configuration with the injection-locking signal.d, The calculated temporal evolution of the Ising Hamiltonian shows energy minimization accompanied by an increase in the graph cut size.For the no annealing scheme, the network converges to a suboptimal solution with a higher energy, while annealing over 250 cycles allows the network to converge to the optimal solution with lower Ising energy and higher probability.Numerical simulations using the same annealing schemes show very good agreement with our experimental results.e, Experimental data and numerically simulated results for the success probability for different annealing times, showing a steady increase from 30% (with no annealing) to 96% (for over 600 cycles of linear annealing).
revealing two global and two local energy minima for the in-phase and out-of-phase configurations, respectively.This establishes that a ferromagnetic interaction can be replicated using resistive coupling in a PTNO-CTDS-based Ising solver.By increasing the strength of the resistive coupling, the attractor states for the in-phase configuration become more prominent.

experimental demonstration of Maxcut with annealing
We investigate the performance of the PTNO-CTDS-based Ising Hamiltonian solver on an NP-hard graph problem of MaxCut for an undirected and unweighted graph.The MaxCut problem statement is defined as follows: given an undirected graph G = (V, E) with V nodes and E non-negative weights on its edges, the problem requires partitioning G into two subsets, namely, W and X, such that the total weight on the edges connecting the two subsets is maximized.The MaxCut problem can be formulated into an equivalent Ising problem using antiferromagnetic interaction (J = −1) and we assume the linear Zeeman term to be zero.The cut size for a given spin configuration σ has a direct mapping to the Ising Hamiltonian H(σ), which is given by C As such, minimizing Ising Hamiltonian H maximizes cut-set C. We chose an undirected and unweighted Möbius ladder graph with eight nodes for our experiment, as shown in Fig. 4a.Phases of the oscillators θ are converted into Ising spins σ using discretization windows in the phase space (Methods).The PTNOs are connected following the adjacency (or connectivity) matrix of the given graph G using coupling capacitances of equal magnitude.The sinusoidal injection-locking signal at twice the oscillator frequency f 0 is applied across the injection capacitances.We implement a linear annealing schedule where the amplitude of the injection-locking signal V inj is linearly ramped from zero to a maximum peak-to-peak value of 10 V over annealing time T anneal .This is followed by the phase readout time T readout .Thus, the total computation time T comp = T anneal + T readout .To experimentally investigate the efficacy of annealing, we vary annealing time T anneal , as schematically shown in Fig. 4b.For all the cases, we keep T readout fixed at 100 oscillation cycles, while T anneal is varied from 0 (representing no annealing) to 660 oscillation cycles.
Figure 4c shows the evolution of the phases of PTNOs for a single run for T anneal = 3.7 ms.The equivalent number of oscillations is calculated as N osc = 250 cycles.The temporal evolution of the Ising energy and the resultant cut-set C is shown in Fig. 4d for the case of no annealing and annealing for 250 cycles.For the case of no annealing, the application of a high V inj immediately binarizes the phases of PTNOs.Hence, the corresponding temporal evolution of the Ising energy shows a steep descent.However, as highlighted in Fig. 2f, the high V inj results in a high dwell time that substantially reduces temporal fluctuations in the oscillator phases.Thus, with very little freedom to escape the local minima, the network converges to a suboptimal solution with a higher energy.On the other hand, when we linearly increase V inj over 250 cycles, the dwell time exponentially increases, as highlighted in Fig. 2f, and the temporal fluctuations in the oscillator phases gradually reduce.The network slowly performs energy minimization with more freedom to escape the local minima and converges to the optimal solution with a higher probability, as shown in Fig. 4d.Thus, we can perform classical annealing in our PTNO-based CTDS by controlling the temporal fluctuation in the oscillator phases, which is very similar to that of simulated annealing with a decaying noise (or temperature).To quantify the performance of the PTNO-based CTDS for varying d, Success probability of finding the MaxCut for random cubic graphs of varying size.We compare the simulation results for three different annealing cycles: 100, 500 and 1,000.e, Success probability of solving dense MaxCut problems with 50% connectivity for varying problem sizes.We consider three different annealing cycles, namely, 100, 500 and 1,000, for comparison.IQR, interquartile range.

Nature electroNics
annealing conditions, we run the network multiple times to calculate the success probability for finding the MaxCut on this graph instance.The success probability is defined as the ratio of the number of trials that returned the true ground-state energy to the total number of trials.To obtain the true ground-state energy for comparison, we run the same graph instance using the Biq Mac solver that executes an exact algorithm (branch-and-bound) on a digital CPU 27 and using the QUBO software called qbsolv (ref. 28) developed by D-Wave.Figure 4e shows the success probability increasing with the annealing cycles.The scenario of no annealing resulted in a success probability of 30% obtained experimentally, while increasing T anneal to over 600 cycles resulted in a success probability of 96%.It is to be noted that the measured performance of our Ising solver is limited by the parasitics introduced by the experimental setup.One major contributor is the parasitic coupling capacitance C C,W arising from the breadboard.The presence of a large C C,W (estimated to be around 22 pF in our experiment) introduces undesired coupling among the oscillators in addition to the intended coupling C C determined by the adjacency matrix and hence lowers the success probability (Supplementary Note 8 provides the details).To compare with the experimental results, we perform numerical simulations for an eight-oscillator network coupled using capacitances, as shown in Fig. 4a.We use a PPV-based framework with experimentally calibrated device and circuit parameters.The simulation details are delineated in Methods and Supplementary Notes 2 and 8.We also introduce the same annealing scheme as our experiments.
The success probability obtained from the simulations shows very good agreement with our experimental results, as shown in Fig. 4e.Overall, the experimental and simulation results validate our proposed methodology of progressively obtaining better solutions through annealing.The excellent agreement of our numerical simulations with the experimental results also enables us to use the simulation framework to predict the performance of PTNO-based CTDS for a larger problem size.

Scaling with problem size
Next, we use the experimentally calibrated PPV-based numerical simulation framework to explore and benchmark the performance of our PTNO-CTDS-based Ising solver for solving MaxCut with increasing problem size.It is to be noted that such analogue computing using CTDS exhibits an inevitable challenge arising from parasitics and variability as we scale up to a larger network size.As such, in our simulation framework, we incorporate the non-idealities such as interconnect/wire parasitics in terms of line-to-ground capacitance, line-to-line capacitance and frequency variability among the oscillators.Coupling capacitance C C is optimally chosen in our simulations such that the parasitic coupling capacitance remains an order or magnitude lower.This ensures that we obtain a high success probability (Supplementary Notes 6 and 8).The simulation details are highlighted in Methods and Supplementary Notes 6 and 7. We first start by investigating the performance of the Ising solver on a 100 node Möbius ladder graph, which is a regular cubic graph of degree three, as shown in Fig. 5a.We use a linear annealing scheme where the amplitude of SHIL is increased over 500 cycles.Figure 5b shows the temporal evolution of the oscillator phases governed by the process of energy minimization.Figure 5c shows the decrease in Ising energy accompanied by an increase in cut size as the network evolves towards the ground-state configuration.We run the simulation 100 times to calculate a success probability of 94%.Next, we extend the investigation to random cubic graphs of degree three for different problem sizes ranging from 10 to 100 nodes.Figure 5d shows the success probability with increasing problem size.We compare the simulation results for three different annealing cycles, namely, 100, 500 and 1,000.Evidently, for 100 node problems, the success probability reduces to 15% for a short annealing time of 100 cycles.Increasing the annealing time to 500 and 1,000 cycles boosts the success probability to 44% and 86%, respectively.Figure 5e shows the success probability for dense MaxCut problems (connectivity of 50%) as a function of problem size and different annealing cycles.For the 100 node graph, the success probability increases from 10% for a short annealing time of 100 cycles to 27% and 80% for 500 and 1,000 cycles, respectively.Overall, we see that for increasing problem size, it becomes difficult to find the global optimal solution.However, decreasing the reduction in success probability can be mitigated by increasing the number of cycles.Note that here we report the success probability for obtaining the absolute ground-state or optimum MaxCut value.
Relaxing the solution accuracy to 99% or 95% of the optimum MaxCut value incurs a much higher success probability.This means that the PTNO-based CTDS is capable of finding near-optimal

Performance comparison with other approaches
Table 1 shows the performance of the PTNO-CTDS-based Ising solver compared with other approaches for solving 100 node, random dense MaxCut problems.We highlight the relevant metrics, such as the time-to-solution value, energy-to-solution value, power dissipated and energy efficiency (calculated as solutions per second per watt).For comparison, we include five different approaches: (1) a well-known simulated annealing algorithm 8 running on an iMac computer with four 3.5 GHz Intel Core i5 processors, (2) a noisy mean-field annealing algorithm running on an NVIDIA GeForce GTX 1080 Ti GPU (ref. 29), (3) a D-Wave 2000Q quantum annealer containing 2,048 qubits (ref. 16), (4) a coherent Ising machine (CIM) based on an optical parametric oscillator with an FPGA feedback loop 14,15 and (5) a discrete-time memristor-based hybrid analogue-digital accelerator implementing a Hopfield neural network (mem-HNN) (ref. 23) (Supplementary Note 10 provides the performance comparison details).We calculate that our PTNO-based CTDS scheme has an energy efficiency of 1.3 × 10 7 solutions per second per watt, which is a fivefold improvement over the recently demonstrated mem-HNN and several orders of magnitude improvement over other candidates such as CPU and GPU, D-Wave and CIM.Such a performance gain can be attributed to (1) the inherent advantage provided by the CTDS approach that allows synchronized updates of all the oscillator phases, (2) short annealing times to obtain a high success probability and (3) low power dissipation of PTNOs.While the success probability obtained in this work is by utilizing a linear annealing scheme, we believe further improvement is possible in terms of exploring better annealing methodologies or modifying the energy function to avoid non-solution attractor states -based CTDS An Ising Hamiltonian solver based on coupled stochastic phase-transition nano-oscillators S. Dutta 1,5 , A. Khanna 1,5 , A. S. Assoa 2 , H. Paik 3 , D. G. Schlom 3 , Z. Toroczkai 4 , A. Raychowdhury 2 and S. Datta 1 ✉

Fig. 1 |
Fig. 1 | overview of PTNo-based cTDS as an Ising Hamiltonian solver.a,b, The combinatorial optimization problem is reformulated in terms of the IsingHamiltonian H defined by the spin vector σ and symmetric coupling matrix J, and mapped onto the Ising solver.Each Ising spin-representing a node in the graph-is represented by an IMT PTNO.The nano-oscillators are coupled to each other using passive elements such as capacitances.Coupling matrix W for the PTNO network is derived from adjacency or coupling matrix J of the Ising model.c, Top left: schematic of a PTNO consisting of a two-terminal phase-transition hysteretic device (VO 2 ) in series with a transistor.Top right: a false-coloured scanning electron microscopy image of the two-VO 2 device with a length of 100 nm.Bottom: as the load line of the series transistor passes through the unstable hysteresis region of the VO 2 device, self-sustained oscillations are created.V out is the output voltage from the oscillator node.L VO2 is the length of the fabricated device.d, Experimental setup of our PTNO-based CTDS.The main computing kernel comprises eight PTNOs connected using the coupling capacitance following coupling matrix W. The SHIL phenomenon is used to create bi-stable oscillator phases, emulating artificial Ising spin.e, The inherent stochasticity present in the PTNO along with a novel technique of gradually reducing the temporal fluctuations in the oscillator phases by increasing the strength of the injection-locking signal S inj (top) is utilized to perform classical annealing and obtain progressively better solutions (bottom).

Fig. 2 |
Fig.2| creating artificial Ising spin using SHIl.a, Three different injection-locking scenarios are shown.The no synchronization case refers to a free-running oscillator with a uniform oscillator phase in the phase space.For FHIL, the oscillator remains injection-locked with the input signal at the in-phase configuration.For SHIL, the oscillator phase gets binarized into the in-phase and out-of-phase configuration with equal probability.b, Measured probability distribution function (PDF) of the oscillator phases for the three scenarios along with the equivalent Ising energy for a single PTNO is shown depicting zero, one and two energy minima for stable phase locking.c, Increasing the amplitude (V inj ) of SHIL forces the oscillator into the bi-stable phase configuration.d, The increase in V inj also reduces the temporal fluctuations between the two stable phases in the presence of noise.The three curves correspond to V inj = 1 V (top), 3 V (middle) and 5 V (bottom).e, Modulation of the energy landscape with V inj.The three curves correspond to V inj = 1 V (top), 3 V (middle) and 5 V (bottom).Increasing V inj increases the energy barrier between the two bi-stable energy minima states.This tightens the phase distribution shown in c. f, The mean dwell time (τ dwell ) spent by an oscillator before hopping between the two phases is plotted against the applied injection-locking amplitude and is found to follow an Arrhenius-type law (ατ 0 e E B η ) with inherent stochasticity (η) playing the role of temperature.Here

Fig. 3 |
Fig. 3 | Replicating ferromagnetic and antiferromagnetic Ising interactions.a, Schematic of pairwise capacitively coupled PTNOs.b, Measured distribution of the phases of the two PTNOs, highlighting the preference of the oscillators to remain in a stable out-of-phase configuration.The calculated energy landscape for capacitive coupling also highlights the presence of global energy minima corresponding to the out-of-phase configuration.c, The measured probability of the out-of-phase configuration (P out-of-phase ) as a function of varying capacitive coupling strength establishes the equivalent antiferromagnetic nature of interaction when compared with a two-spin Ising model.The out-of-phase probability P ↑↓ is calculated from the Boltzmann distribution.d, Schematic of pairwise resistively coupled PTNOs.e, The measured phase distribution for resistive coupling along with the calculated energy landscape showing the preference for in-phase configuration.f, The measured probability of in-phase configuration (P in-phase ) versus the coupling strength highlights the ferromagnetic nature of interaction for resistively coupled oscillators.The in-phase probability P ↑↑ is calculated from the Boltzmann distribution.

Fig. 4 |
Fig.4| experimental demonstration of Maxcut and performance enhancement with annealing.a, An undirected and unweighted eight-node Möbius ladder graph used for investigating the MaxCut problem.The PTNOs are connected following the same adjacency (or connectivity) matrix of the graph using coupling capacitances.b, Schematic of the annealing schedule used in the experiment.A sinusoidal injection-locking signal at twice the oscillator frequency f 0 is applied with a linearly increasing amplitude over annealing time T anneal that corresponds to a linear annealing schedule.This is followed by phase readout time T readout .Thus, the total computation time T comp = T anneal + T readout .Annealing time T anneal is varied from 0 (corresponding to no annealing) to 660 oscillation cycles.c, Evolution of the phases of the eight oscillators, settling to either the in-phase or out-of-phase configuration with the injection-locking signal.d, The calculated temporal evolution of the Ising Hamiltonian shows energy minimization accompanied by an increase in the graph cut size.For the no annealing scheme, the network converges to a suboptimal solution with a higher energy, while annealing over 250 cycles allows the network to converge to the optimal solution with lower Ising energy and higher probability.Numerical simulations using the same annealing schemes show very good agreement with our experimental results.e, Experimental data and numerically simulated results for the success probability for different annealing times, showing a steady increase from 30% (with no annealing) to 96% (for over 600 cycles of linear annealing).

Fig. 5 |
Fig.5| Scaling with problem size.a, Schematic of a 100 node Möbius ladder graph used for numerical simulation.b, Evolution of the phases of the oscillators for a single run as a function of oscillation cycles.c, Increase in the graph cut size accompanied by a decrease in the equivalent Ising energy as the system evolves towards the ground-state configuration.The simulation was performed for 100 trials to calculate a success probability of 94%.The light coloured curves show the simulation results for 100 different trials and the dark coloured curve shows one representative simulation result.d, Success probability of finding the MaxCut for random cubic graphs of varying size.We compare the simulation results for three different annealing cycles: 100, 500 and 1,000.e, Success probability of solving dense MaxCut problems with 50% connectivity for varying problem sizes.We consider three different annealing cycles, namely, 100, 500 and 1,000, for comparison.IQR, interquartile range.

Fig. 6 |
Fig. 6 | Performance evaluation of PTNo-based Ising solver.a, The total time-to-solution value for solving dense MaxCut problems with 50% connectivity for different annealing times and increasing problem size.A non-trivial dependence of time to solution on the annealing time is seen, where shorter annealing time is preferred for smaller problems and the success probability remains insensitive to annealing time.Longer annealing is preferred for larger problems where the success probability dominates.b, Energy-to-solution value for solving dense MaxCut problems is plotted for different annealing times and for varying problem sizes showing a similar trade-off.

Table 1 | Performance comparison between PTNo-cTDS-based Ising solver and other state-of-the-art approaches for solving Maxcut on 100 node random cubic graphs
Simulated annealing algorithm implemented on a CPU using four Intel Core i5 processors each running at 3.5 GHz.b Hybrid scheme updates 10 nodes at a time.c Reported for 50 cycles (running sequentially).d Reported for 1,000 cycles (running sequentially).e 10 4 s for N = 55. a