Ising Hamiltonian Solver using Stochastic Phase-Transition Nano-Oscillators

Computationally hard problems, including combinatorial optimization, can be mapped onto the problem of finding the ground-state of an Ising Hamiltonian. Thus, physical annealing systems such as Ising machines can be employed to accelerate the search for their solution. Here, we leverage the continuous-time dynamics of a network of stochastic phase-transition nano-oscillators and construct an electronic Ising Hamiltonian solver operating at room temperature. We fabricate a prototype network of capacitively coupled eight SHIL oscillators that finds the ground state of an Ising Hamiltonian and solves the maximum-cut (Max-Cut) in non-planar graphs - a problem of non-deterministic polynomial time (NP) complexity - with high probability of success. Numerical simulations are utilized to find the Max-Cut in non-planar graphs as large as 800 nodes, demonstrating up to six orders of magnitude improvement in time-to-solution and nearly ten orders of magnitude improvement in energy-to-solution compared to standard heuristic algorithms running on a state-of-the-art digital computer.


INTRODUCTION
Combinatorial optimization is ubiquitous in various fields such as artificial intelligence, bioinformatics, drug discovery, cryptography, quantitative finance, operation research, resource allocation, trajectory and route planning. Such problems often belong to the NP-hard or NPcomplete complexity class, requiring computational resources (time or energy) that scale exponentially with the problem size. Exact methods such as branch-and-bound are often limited to problem sizes of only few hundred variables. Alternatively, approximate algorithms and heuristics such as relaxation to semidefinite program (1), simulated annealing (2) and breakout local search (3) are widely used in digital computers to find an optimal or near-optimal solution. However, even for moderately sized problem instances, the time to find a near-optimal solution can become prohibitively large. Hence, there is a growing interest towards finding alternative approaches that can solve large-scale problems efficiently. One promising approach is to build special purpose physical computing substrates or accelerators that can determine the ground state of an Ising model (4) or its equivalent Quadratic Unconstrained Binary Optimization (QUBO) problem (5) in a timescale much shorter than a digital computer, with less hardware resources. Various techniques have been developed to decompose large-scale problems into smaller QUBO problems that are tractable with today's digital accelerators (6,7). By running such special purpose machines in parallel one can reduce the correlation between consecutive samples and effectively search the solution space with multiple independent runs (8). Furthermore, hybrid approaches such as augmenting the fast approximate solution offered by the optimization solver with other metaheuristic local-neighborhood search algorithms such as Tabu search (9) can further improve the quality of the final solution.
The Ising model, describing the properties of spin glass, was put forward as a tool of statistical physics to explain ferromagnetism. The Ising Hamiltonian with discrete spins !"#"$ ∈ {−1, +1} $ and a symmetric coupling matrix in the absence of an external magnetic field is given by Finding the ground state of Ising model belongs to the NP-hard complexity class (10) and can be extended to other NP-hard and NP-complete problems including all of Karp's 21 NP-complete problems through polynomial time mapping (11). Recently, various schemes for building special non-von Neumann processors, called Ising machines have been devised using a variety of technology platforms. These include super-conducting qubit-based adiabatic quantum computing (AQC) and quantum annealing (12,13), digital and mixed-signal complementary metal-oxidesemiconductor (CMOS) annealers (14)(15)(16), memristive Boltzmann machine (17), memristor-based hardware accelerators for classical annealing (18) nanomagnetic arrays (19), coherent networks of degenerate optical parametric oscillators (20,21) and synchronized networks of electronic oscillators (22)(23)(24)(25). Qubit-based quantum annealers incur high cost and complexity arising from cryogenic operating temperature requirement. The optical coherent Ising machine has shown competitive performance compared to the quantum annealer (26). However, they utilize a kilometer long fiber ring cavity for implementing 2,000 spins using time-multiplexing and require extremely fast (and power hungry) field-programmable gate array (FPGA) for implementing coupling in a measure-and-feedback scheme (21). CMOS annealers relies on an external source for random number generator and finds it technologically challenging to maintain true randomness in CMOS implementation and requires significant post-processing.
To address the above shortcomings (such as , we devise and experimentally demonstrate a chipscale insulator-to-metal phase-transition (IMT) nano-oscillator-based Ising Hamiltonian solver operating entirely in the electronic domain at room-temperature. The hardware features the key ingredients necessary for implementation of the Ising Hamiltonian, namely representing the Ising spin and the coupling between the spins. We establish the mathematics of continuous-time dynamics of a network of coupled IMT oscillators and establish the analogy between the "energy" or Lyapunov function of such a dynamical system and the Ising Hamiltonian. Finally, we utilize the inherent stochasticity of the insulator to metal phase transition in IMT oscillators (27) and controllable crossover between unsynchronized and synchronized states to implement an annealing schedule that allows escaping the local minima. We use Max-Cut problem of NP-hard complexity as a benchmark task to quantify the performance of this Ising solver and quantify the advantages in terms of time-to-solution and energy-to-solution over classical digital machines.

Overview of an artificial Ising solver with IMT oscillator network
An overview of an IMT nano-oscillator-based Ising Hamiltonian solver system is shown in Fig.  1(a). The optimization problem is reformulated in terms of the Hamiltonian of an Ising spin glass defined by the spin vector and the symmetric coupling matrix and mapped onto the Ising solver. The main computing kernel comprises of a network of IMT nano-oscillators. Artificial Ising spins are generated by applying external injection locking signals to the oscillators, while the ferromagnetic or anti-ferromagnetic interaction of the Ising Hamiltonian is established using resistive and capacitive elements used for coupling the nano-oscillators as explained later. In this solver, the phases of the nano-oscillators ⃗ are used as state variables for computing. The continuous-time dynamics of the network is dictated by an "energy" or the Lyapunov function ( ⃗ ) that closely resembles the Ising Hamiltonian . The dynamics of the network evolves continuously in time, so as to naturally minimize ( ⃗ ), and in the process minimizes the Ising Hamiltonian to reach the ground state. Once the network reaches an energy minimum state, the output of the network is read out in terms of the state variables ⃗ and subsequently reformulated to provide the optimal solution of the original optimization problem. One daunting task in solving such combinatorial optimization problem is the increase in the complexity of the energy landscape with the 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 machine and simulated annealing, to help with 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 reduced slowly 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 IMT material (28) to escape the local minima. We introduce a novel way to modulate the energy landscape (or the energy barrier ' ) by tuning the amplitude of the injection locking signal that, in turn, controls the probability of stochastic transitions between the adjacent valleys. As such, this Ising Hamiltonian solver naturally implements classical annealing by exploiting the continuous time dynamics of the underlying physical system.

Artificial Ising spin
To emulate an Ising spin system, we require a coupled network of elements with binary degrees of freedom. We use a network of injection locked IMT nano-oscillators as an artificial Ising spin system. The IMT nano-oscillator comprises a two-terminal phase-transition hysteretic device in series with a transistor. We use Vanadium Dioxide (VO2) as a prototypical IMT material in our experiments. The device dimensions used in our experiments are 1 m in length and 2 m wide. Below the insulator to metal phase transition temperature and under zero external electric field or current, VO2 shows insulating behavior. Upon 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 insulating to metallic state. The hysteretic phase transition is reflected as an abrupt hysteretic current-voltage (I-V) characteristic of the device (See Materials and Methods and Supplementary Information section S2 for details). Pairing a conventional n-type metal-oxidesemiconductor (NMOS) transistor in series with the VO2 such that the load line passes through the unstable region of the I-V curve results in self-sustained oscillations. We used VDD = 2V and a gate voltage VGS = 0.9V in our experiments. The measured fundamental frequency of the oscillator The binary degrees of freedom in the phase space of an IMT nano-oscillator arises from the phenomenon of second harmonic injection locking. We apply a sinusoidal injection locking signal #)% = #)% (2 #)% ) with #)% = 2 ( as an external input to the oscillators across the capacitor #)% = 20 as shown in Fig. 1(a). When #)% = 0, the oscillator is free running. Fig. 1(b) shows the voltage output waveform *+, measured over multiple runs. The corresponding phase of the oscillator, measured with respect to a reference sinusoidal signal of same frequency ( , shows a constantly varying phase with a uniform probability distribution over the entire phase space as seen Fig. 1(c). In contrast, when the oscillator is perturbed with #)% at the first harmonic, i.e. #)% ≈ ( also referred to as first harmonic injection locking (FHIL), *+, shows a constant inphase locking configuration with #)% . The probability distribution of the phase, measured over multiple runs, shows a single gaussian distribution around 50 0 . Interestingly, when #)% ≈ 2 ( , the oscillator waveform shows both in-phase and out-of-phase injection locking configuration when measured over multiple runs. The corresponding probability distribution exhibits a double gaussian distribution highlighting an equiprobable and bi-stable phase portrait. This bi-stability provides an ideal means to encode the Ising spin in the electrical domain, where phase 40 0 represents up-spin, i.e., = +1 , and phase 220 0 represents down-spin, i.e., = −1. The dynamics of the phase difference between the oscillator and the injection locking signal is described by a generalized version of Adler's equation (Gen-Adler) given by where 1 is the n th harmonic of the IMT nano-oscillator and #)% 1 = 2 1 * #)% #)% #)% . The first term describes the frequency mismatch between *+, and #)% , which contributes to phase slipping. The second term depends on the phase delay incurred due to the perturbation caused by #)% and is described in terms of the perturbation-projection-vector (PPV), (see Materials and Methods section and Supplementary information section S4 for details). The corresponding "energy" or Lyapunov function of the oscillator is given by The first energy term comes from the contribution of the frequency mismatch that creates an overall bias in the energy landscape. However, this being a linear term does not introduce any new 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 over the phase space can be calculated as # is the partition function and is analogous to the term in the Boltzmann distribution 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. 1(d) indicating a uniform distribution over the phase space as is experimentally obtained in Fig. 1(c). For first harmonic injection locking (FHIL), as long as the relative frequency difference E #)% − * F is small compared to the second term in equation, -, exhibits two stable points at ≅ 0.2 and ≅ 1.2 , i.e. two equilibrium phases. The calculated ( ) (assuming zero frequency mismatch) evolves into a double well energy landscape that results in a double gaussian distribution in the phase space as reproduced faithfully in the measurements (see Materials and Methods for detail).
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 the injection signal #)% below which no stable solution exists. In our experiments, this critical amplitude is close to 1V for a frequency mismatch of 3% as observed in the measurement. Above this threshold, the bi-stability in the phase space begins to appear. Fig.  2(a) shows the measured oscillator phase upon varying #)% . For very low #)% close to the critical value, the phase of the oscillator measured over multiple runs shows a wide distribution in the phase space, while the distribution narrows as #)% increases. This behavior can be understood by considering the perturbation in the energy landscape as shown in Fig. 2(b) for #)% = 1V, 3V and 5V. For #)% = 1 , the energy barrier ' separating the two stable equilibrium phases is low (around 0.006 a.u. calculated from Eq. 3). Hence, in the presence of stochastic noise, the oscillator's phase fluctuates between the two stable phases as seen in the Fig. 2(c). An increase in #)% to 3V and 5V increases the barrier height ' to 0.012 a.u. and 0.02 a.u., 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 ;<4== (analogous to Neel relaxation time for magnetization) as a function of #)% as shown in Fig. 2

Replicating the interaction terms in the Ising Hamiltonian
To implement an asynchronous network of IMT nano-oscillators 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 shown schematically in Fig. 3(a) using capacitance ? as the coupling element. Fig. 3(b) shows the experimentally measured phase distribution of the two oscillators using an injection locking capacitance #)% = 20 and coupling capacitance ? = 56 . The oscillators remain out-of-phase with each other and the two configurations -(40 0 , 220 0 ) and (220 0 , 40 0 ) are two equally probable states. We measured the probability of out-of-phase configuration for varying coupling strength. The probability remains close to 0.5 for low coupling capacitance, meaning both in-phase and out-of-phase configurations are equally probable, and increases close to 1 for higher coupling. To understand the exact nature of capacitive coupling, we compare it with a 2-spin Ising model with antiferromagnetic interaction (negative ) where the individual spins prefer to remain anti-parallel with one another. In such a system, the probability of one of the possible configurations-(↑↑, ↓↓, ↑↓, ↓↑) is determined by the Boltzmann distribution ( ! , 2 ) = 61(@ + ,@ , )/BC ⁄ , where is temperature and = ∑ 61(@D + ,@D , )/BC @D + ,@D , is the partition function. With antiferromagnetic interaction (negative ), the probability for anti-parallel configuration (↑↓, ↓↑) increases upon varying the interaction strength from = 0 to −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 dynamics of such a coupled oscillator network is further described by extending Eq. 2 and Eq. 3 to incorporate an additional coupling term and is given by where *EF,% = ?,% -H -./ ,% -, . The additional third term describes the coupling interaction energy between pairs of oscillators. The corresponding "total energy" or the global Lyapunov function of the system is then given by We use Eq. 5 to calculate the two-dimensional energy landscape for pairwise capacitively coupled IMT nano-oscillators as shown in Fig. 3(d). The calculated energy exhibits the presence of two degenerate global minima at the out-of-phase configuration and two additional degenerate local minima at in-phase configuration. The corresponding probability distribution of the phases of the coupled oscillators over the two-dimensional phase space can be calculated as ( ! , 2 ) = is the partition function and accurately describes the experimentally obtained result in Fig. 3 Fig. 3(d) are two calculated phase trajectories of the oscillator network as it settles in the two equiprobable energy minima.

(b). Also shown in
A similar investigation is performed for pairwise resistively coupled oscillators as shown schematically in Fig. 3(e). Contrary to the previous case, the measured oscillator phases reveal a higher probability to be in-phase with each other in either (40 0 ,40 0 ) or (220 0 ,220 0 ) configuration as shown in Fig. 3(f) for a coupling resistance of ? = 40 Ω. To establish the nature of resistive coupling, the probability of in-phase configuration is measured as a function of varying ? and compared to a 2-spin Ising model with ferromagnetic interaction (positive ). Note that a lower ? represent a higher coupling strength and hence a higher . The increase in the probability of inphase configuration with decreasing ? , i.e. increasing coupling strength, agrees well with the theoretical Ising model that predicts an increase in probability of parallel configuration (↑↑, ↓↓) for increasing as shown in Fig. 3(g). We also use Eq. 5 with *EF,% = − to calculate the twodimensional energy landscape as shown in Fig. 3(h), revealing global and local energy minima for in-phase and out-of-phase configurations, respectively. The corresponding probability distribution of phases ( ! , 2 ), calculated from the energy landscape, is in excellent agreement with the experimental results, thus establishing that a ferromagnetic interaction can be replicated using resistive coupling in the IMT nano-oscillator based Ising solver. The two simulated trajectories of the oscillator network settling to the two equiprobable energy minima are also shown.

Experimental Demonstration of Solving Max-Cut
We investigate the performance of the IMT nano-oscillator-based Ising Hamiltonian solver on a NP-hard graph problem of Max-Cut for an undirected and unweighted graph. The Max-Cut problem statement is defined as: Given an undirected graph G = (V, E) with V nodes and E nonnegative weights on its edges, the problem requires partitioning G into two subsets W and X such that the total weight on the edges connecting the two subsets is maximized. The Max-Cut problem can be formulated into an equivalent Ising problem using antiferromagnetic interaction ( = −1) and we assume the linear Zeeman term to be zero (See Materials and Methods for details). The cut size for a given spin configuration ⃗ has a direct mapping to the Ising Hamiltonian ( ⃗), given by As such, minimizing the Ising Hamiltonian maximizes the cut-set . We chose an undirected and unweighted Mobius Ladder graph with 8 nodes as shown in Fig. 4(a) for our experiments. The phases of the oscillators ⃗ are converted to the Ising spins ⃗ using discretization windows in the phase space (see Materials and Methods for detail). The IMT nano-oscillators 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 ( is applied across the injection capacitances. We implement a linear annealing schedule where the amplitude #)% is linearly ramped from zero to a maximum of 10V peak-to-peak over an annealing time M))4M= . Fig. 4(b) shows, for a single run, the evolution of the phases of the oscillators for M))4M= = 4.6 . The equivalent number of oscillations is calculated as *EF = M))4M= ( , where ( is the operating oscillator frequency in our experiments. The temporal evolution of the "energy" of the network E ⃗ F, calculated from Eq. 5 using the measured oscillator phases, is shown in Fig.  4(c) for multiple runs. We see the network minimizing the energy and evolving towards the lowest energy state. Fig. 4(d) depicts the Ising solver maximizing the cut set while minimizing the equivalent Ising Hamiltonian . We run the Ising solver multiple times to calculate the success probability for finding the Max-Cut on this graph instance. The success probability is defined as the ratio of the number of runs that returned the true ground-state energy to the total number of runs. To obtain the true ground-state energy for comparison, we run the same graph instance on the BiqMac solver that runs an exact algorithm (branch-and-bound) on a digital CPU (29). Over 62 runs, we obtained the ground state 56 times resulting in a success probability of 90%. We also performed similar experiments for a shorter and longer anneal time of M))4M= = 1 and 9.6 , respectively. The success probability increased to 96% as the anneal time is increased as shown in Fig. 4(e). The distribution of the obtained cut set or the equivalent ground state energy is shown in the inset of Fig. 4(e). Besides finding the ground state for a large fraction of the runs, the IMT nano-oscillator Ising solver also finds the cut sets or the equivalent energy states that are in close proximity to the global minimum, thus providing acceptable approximate solutions. The dependence of the success probability on the annealing time can be understood by looking at the distribution for the number of oscillation cycles required to reach the ground state energy. As shown in Fig. 4(f), for most of the times, the ground state energy is reached within a period of 300 oscillations ( M))4M= = 4.5 ) giving us a high success probability of 90%. This explains the marked increase in the success probability when increasing M))4M= from 1ms to 4.5ms. A further increase in the anneal M))4M= allows the solver to also capture the worst-case scenarios where the number of required oscillation cycles reaches above 500 oscillations as seen in Fig. 4(f).

Scaling with Problem Size
Here, we investigate the performance of the IMT nano-oscillator-based Ising solver for problems of large size. We chose Mobius Ladder graphs up to 800 nodes for which the ground state energy and Max-Cut solution was first obtained using the BiqMac solver. For our Ising solver, we use a PPV-based numerical simulation methodology as delineated in the Materials and Methods section later. For obtaining the PPV function , we use a SPICE compatible macro-model of the IMT nano-oscillator (30) to quantitively match the dynamics of the oscillators under SHIL condition, and perform cycle accurate time domain simulations of the IMT oscillator using the Cadence Spectre circuit simulation framework (31). The details are described in the Materials and Methods section. We additionally optimize the values for the coupling capacitance ? and the oscillator jitter noise to obtain the best possible cut-set values. Fig. 5(a) shows the numerical simulation result of the IMT nano-oscillator-based Ising machine solving a 200 node Mobius ladder graph for a single run. The phases of the oscillators evolve as a function of time or oscillation cycles, governed by the process of energy minimization of the overall system. Fig. 5(b) shows the increase in the graph cut size accompanied by a decrease in the equivalent Ising energy as the system evolves towards the ground state configuration. We see that the ground state energy is reached within 400 cycles for most cases. To calculate the success probability of finding the ground state, we re-run the simulation 100 times for each graph size. Fig. 5(c) shows the success probability as a function of the graph size. We plot E+FF4EE,N% , which is the success probability of obtaining a solution greater than % of the optimal Max-Cut value in a single run, with x = 100%, 99% and 98%. It is interesting to see that, for larger graph sizes like N = 500, even though the solver returns the optimal solution only 5 times out of the total 100 runs ( E+FF4EE,!((% = 5%), the solver is able to find sub-optimal solutions with higher probability ( E+FF4EE,PP% = 37% and E+FF4EE,PQ% = 88%). This can be understood by looking at the distribution of the cut set solutions obtained for N = 60, 200 and 500 nodes as shown in Fig. 5

(d).
A key metric for benchmarking the performance of any Ising solver is the total computation time required to obtain at least one solution, which is either the ground state or within an accuracy range of 98% -99% of the ground state, with a probability of 99%. The total time to solution is calculated directly from Fig. 5 With increasing graph size, the success probability of a single run decreases exponentially, necessitating the solver to re-run the problem several times (either serially on the same solver or in a batch mode over multiple solvers) for ensuring a 99% cumulative success probability of obtaining at least one solution which is either the ground state or within an accuracy range of the ground state (see supplementary section S8 for details). As shown in Fig. 5(e), the total time-to-solution as a function of the graph size also follows an exponential nature ( R$ ) highlighting the NP-hard complexity of the problem. However, in absolute terms, the IMT nano-oscillator-based Ising solver, assumed to be running at 100 MHz in our simulations, is able to find the optimum Max-Cut solution of a 300 node Mobius ladder graph in just 134 s. For sub-optimal solutions with 98% accuracy, the total time to solution is 17 s. To highlight the advantage of IMT nano-oscillator based Ising solver, we calculate the time-to-solution required for running heuristic algorithms for the same Mobius ladder graphs up to 800 nodes as shown in Fig. 5(e). For instance, the BiqMac solver that runs the exact branch-and-bound algorithm on a CPU takes around 547s to find the optimal Max-Cut solution of a 300 node Mobius ladder, while a simulated annealing (SA) algorithm running on an iMac computer with 3.5 GHz Intel Core i5 processor take around 103s. Thus, our simulation results show over five orders of magnitude improvement in total time to solution over running an approximate algorithm such as SA and nearly six orders of magnitude improvement over the exact branch and bound algorithm, both running on a digital computer.
We also investigate another key metric for benchmarking the performance -the energy-to-solution for solving such graph problems. For the IMT nano-oscillator Ising solver, the average power consumption for the main compute kernel, i.e. the coupled oscillator network, is estimated to be around 10 per oscillator as obtained from circuit simulations. We incorporate an additional energy overhead of 20% that will be incurred from the peripheral circuits such as for reading out the phases of the oscillators (the details of the peripheral circuit design are beyond the scope of this work). The total consumed energy-to-solution is then estimated considering the total time to solution from Fig. 5(e). Fig. 5(f) illustrates the energy consumption of our Ising solver for obtaining solutions with different accuracies for 99% cumulative success probability. For comparison, we also estimate the consumed energy-to-solution from running the BiqMac solver and SA algorithm on digital computers. We conservatively budget an average CPU power consumption of 15W per core and 4 cores per CPU. From the simulation results, the IMT nanooscillator based Ising solver outperforms digital systems running heuristic algorithms by ten orders of magnitude.

DISCUSSION
The notion of solving hard optimization problems using the continuous dynamics of a physical system opens novel avenues of exploration of dynamical systems for compute applications. There is much enthusiasm in building special purpose machines (or accelerators) for solving graph problems belonging to the NP-hard and NP-complete complexity class as part of a strong push towards a heterogenous compute platform. There is a rapidly growing demand to analyze and uncover hidden relationships between similar or diverse datasets in real-time and service applications such as customer analytics, risk and compliance management, recommendation engines, route optimization, fraud detection, asset allocation and risk management. Hence, we are witnessing a resurgence in building dedicated optimization processing units (such as Ising Hamiltonian solvers) that can complement general-purpose central processing units (CPU) and graphics processing units (GPU). Ising solvers are gaining attraction in real-life as many relevant NP-hard and NP-complete problems can be reformulated into the problem of finding the ground state of an Ising model (11). In this work, we propose and demonstrate phase-transition nanooscillator based approach that enables a chip-scale, room temperature compatible, all-electrical, energy efficient and high-performance Ising Hamiltonian solver. While the concept of utilizing a network of coupled electronic oscillators for solving optimization problems have been initially proposed by Wang et. al. (32) and later also studied by Cho et. al. (24) and Ahmed et. al. (25), such hardware involved bulky LC oscillators and ring oscillators with latch-based coupling. In comparison, we demonstrate a novel hardware implementation of a capacitively coupled oscillator network using one-transistor and one-resistor (1T-1R) IMT nano-oscillator (see Supplementary information section S9 for details). The ability to harness the intrinsic stochasticity in phase-transition devices such as VO2 allows us to traverse complex energy landscapes and escape local minima to reach an optimum solution without requiring any additional random number generator circuit for generating noise.
The IMT nano-oscillator Ising solver provides several advantages when compared with other contemporary Ising machine proposals. For example, D-Wave's 2000Q quantum annealer containing 2,048 qubits can compute the Max-Cut of a 200 node cubic graph (degree = 3) in 11ms which is around 200x higher than the time required by our Ising solver (26). Moreover, D-Wave's 2000Q currently has a specific "chimera" coupling architecture with sparsely connected qubits which requires further minor embedding techniques such as native clique embeddings to find chimera subgraphs (33). D-Wave's quantum annealer requires around 25kW power to operate at cryogenic temperature. We expect this overhead from cooling to be constant with scaling up to thousands of nodes. The coherent Ising machine (CIM) based on optical parametric oscillator proves to be another promising candidate. Recently, a fully connected 2000 node CIM has been demonstrated. However, it requires one-kilometer long fiber cavity to accommodate the degenerate optical parametric oscillator (DOPO) pulses and incurs a cavity round trip time of 5 . This limits the time-to-solution to around 51ms to solve Max-Cut problem on a 200 node cubic graph which is around four orders of magnitude slower than that required in the IMT nano-oscillator Ising solver (26). Moreover, the CIM requires additional DOPO pulses and additional mechanisms for maintaining phase stabilization (20,21). One key advantage in our Ising solver is the ability to replicate bi-directional ferromagnetic and anti-ferromagnetic coupling using simple electrical elements and thus achieve highly parallelized connectivity. This is in contrast to an elaborate measurement-and-feedback schemes used in CIM for coupling the spins, requiring a 12-bit analogto-digital converter (ADC), multiple field-programmable-gate-arrays (FPGAs), and high speed transceivers to connect the FPGAs (21) and is inherently serial.
Besides D-Wave and CIM, several CMOS-based digital or mixed-signal hardware accelerators have been recently proposed. For example, the Ising chip by Hitachi uses CMOS static random access memory (SRAM) cells as spins while the coupling is realized using digital logic gates (15). As such, it provides a hardware implementation of the classic simulated annealing algorithm by performing deterministic on-chip computation paired with an external random number generator is simulate. Due to the non-von Neumann architecture, it exhibits a 50x lower energy-to solution over running a SG3 algorithm (34) on a CPU for finding Max-Cut of a 200 node random cubic graph. However, the energy dissipation is still orders of magnitude more than our Ising solver. Additionally, the analog or continuous-time dynamics of our solver provides a key potential advantage which lowers the time to solution compared to a CMOS annealer operating in discrete time. Along the lines of improving the energy efficiency of hardware accelerators, recently analog in-memory computing in resistive-RAM (RRAM)-based crossbar arrays has been proposed that provides four orders of magnitude improvement in energy over CPU (18). In comparison, the energy of our IMT nano-oscillator-based Ising solver still remains around three order of magnitude lower than the RRAM-based hardware, mainly due to the ultra-low power dissipation of the IMT oscillators. The time-to-solution (or cycles-to-solutions), however, remains almost similar for both our solver and the RRAM-based hardware accelerator.
This work implements a simple annealing schedule for finding the ground state of the Ising model. While the success probability decreases as the problem size increases using this simple annealing schedule as seen in Fig. 5(a), further improvement can be made in terms of investigating better annealing methodologies or modifying the energy function to avoid non-solution attractors that trap the system in local minima (35,36). Hybrid approaches can also be adopted to improve the quality of solution, such as augmenting the search of an Ising solver in the first phase with other metaheuristic local-neighborhood search algorithms such as Tabu search (9) in the second phase. Another avenue of future work is to introduce reconfigurability of the oscillator network such that we can solve a diverse set of combinatorial optimization problems. Our current implementation with passive capacitive elements requires human-in-the-loop to set the coupling strengths depending on the problem being solved. A promising pathway will be to investigate switched capacitor based coupling which will allow multiple levels of programmable ferromagnetic and anti-ferromagnetic type interaction between oscillators with a single control signal. This would not only automate coupling network configuration, but also open the possibility of exploring additional annealing schemes using time varying coupling terms to augment the search for the ground state solution. The current breadboard-based connections of (non-VO2) components such as series transistors, coupling elements and injection locking sources lowers the overall frequency of the oscillators due to large inherent parasitic capacitances and limits the ability to couple large number of oscillators. As such, the ultimate pathway towards realizing a large-scale phase-transition nanooscillator based Ising solver chip requires developing wafer-scale growth of high quality phasetransition oxides using CMOS compatible processes and subsequently building arrays of fully integrated IMT oscillators with fully programmable on-chip coupling. Finally, it is to noted that, while an all-to-all connected oscillator network may not be practically feasible to map a largescale real-life problem, various avenues of research remain open to investigate methodologies for decomposing such large-scale problems into smaller Ising problems that are fully compatible with this oscillatory network.

Sample preparation
10nm think Vanadium dioxide (VO2) is grown on a substrate of (001) TiO2 substrate using Veeco Gen10 molecular beam epitaxy (MBE) system. The widths of the two terminals are defined by dry etching with CF4. The device length is defined by depositing Pd/Au metal contacts using electron beam evaporation. The fabricated VO2 devices varied in length from 100nm to 1um with resulting insulator-to-metal transition threshold voltages ranging from 0.7V to 4V. All our experiments have been performed on devices with dimensions 1 × 2 . Fig. 1(a) shows the schematic of an IMT oscillator which is realized by connecting an n-channel MOSFET (ALD1103) transistor in series with the two-terminal VO2 device. A VDD of 2V is applied and the amplitude of the relaxation oscillations is ~1.7V. A gate voltage VGS = 0.9V is applied to the series transistor that set the oscillation frequency ( ≈ 100 and a full-with-halfmaxima (FWHM) = 200Hz. A schematic of the full experimental setup for the Ising solver is shown in supplementary figure Fig. S1. Eight VO2 devices places in the Keithley 4200-SCS probe station were connected using multi-contact probes. The VDD and the analog gate voltages VGS of the 8 series transistors are applied using a multichannel analog voltage card connected to a computer. The injection locking signal was applied to the 8 oscillators across injection locking capacitances #)% = 20

Experimental setup
using an external voltage generator. #)% was realized using discrete off-chip capacitors connected on a breadboard. The output voltage waveforms of the oscillators were measured using a multichannel digital oscilloscope. The coupling among the oscillators was realized using discrete off-chip coupling capacitances F = 56 connected on the breadboard.

Data-processing
The output voltage waveform of the oscillators was measured using a multichannel digital oscilloscope and subsequently analyzed and processed in MATLAB on a digital computer. The phase of an oscillator is calculated with reference to a reference sinusoidal signal with the same frequency ( as the oscillator and half of the injection locking signal. The phase is defined as the time difference the minima point of the discharging phase of the IMT oscillator and the minima of the sinusoidal signal divided by the time period of the oscillator (see supplementary figure Fig.  S5). The measured oscillator phases were converted to Ising spins using a discretization window in the phase space such that if 90 ( < < 270 ( , = 1, else = 0. Subsequently, the Ising energy and the cut set are calculated following Eq. 1 and Eq. 6 and = −1. We consider that the Ising solver has reached the minimum energy state if the phases of the oscillators and the calculated cut set does not change over 10 oscillation cycles after reaching the known final configuration.

Theory
The Adler's Equation describes the evolution of the phase difference of a free-running sinusoidal oscillator (or LC oscillator) with respect to a sinusoidal (AC) injection locking signal. The generalised Adler's Equation (37) given by Eq. 2 extends this to any non-linear oscillator with an arbitrary periodic injection locking signal. Analytical and numerical techniques to predict the locking behaviour for LC and ring oscillators have been shown (38) based on the Perturbation Projection Vector (PPV). The PPV tracks the excess phase delay generated in a free-running oscillator from a small current (or voltage) perturbation at the oscillating node. We numerically calculate the PPV (ξ( )) from the impulse response of the IMT oscillator through SPICE simulations (see Supplementary Section S4 for details). The phenomenon of injection locking involves a continuous perturbation from an input signal ( ( )) such as a sinusoidal signal with a fixed frequency. The excess phase delay of the oscillator due to this external perturbation signal can then be calculated by multiplying the PPV with the input signal. The excess phase delay generated over one time period is given by where ( ) is the initial phase difference between the oscillator and injection locking signal. The integration is performed in the phase space over one oscillation period. For FHIL, we can write the generalised Adler's equation describing the rate of change phase difference between the oscillator and signal as The steady-state analysis produces two solutions for ,one of which is stable. Note that for the producing a stable solution, the relative frequency difference between the oscillator and the injection locking signal must be smaller than the magnitude of ( ) (see supplementary section S4 for details). In the case of SHIL, two stable solutions are produced which are shown graphically in the supplementary Fig. S4 Thus, the "energy" of a single injection locked oscillator is given by The presence of one stable solution for FHIL and two stable solutions for SHIL is also predicted by the calculated energy landscape that shows a single energy well and double energy well for FHIL and SHIL, respectively ( Fig. 1(d)).
The gen-Adler equation is be further extended to describe the interaction between the pair of nonlinear oscillators which results in a coupled oscillator system. We can write a coupled pair of differential equations as Thus, overall, the dynamics of the network including injection locking and mutual coupling between oscillators can be written as and the corresponding "total energy" or the global Lyapunov function of the system is then given by

Large Scale Simulation of Oscillator Network
The dynamics of the IMT oscillator network was simulated by converting Eq. 16 into a stochastic differential equation and numerically solving it using the Euler-Maruyama method. The stochasticity of the IMT nano-oscillator is implemented by adding a Gaussian phase noise twice every oscillation cycle. This replicates the inherent stochasticity in abrupt transition from insulator to metallic phase and vice versa occurring twice in every oscillation cycle. The PPV function (ξ( )) and oscillator currents ( *EF ) in Eq. 16 was calculated though SPICE simulations using the Cadence Spectre circuit simulator (31). The hysteretical VO2 SPCIE model is based on a previous work (30) and the circuit schematic is shown in supplementary figure Fig. S3. In our simulations, the insulator-to-metal transition voltage VIMT was considered as 0.7V and a VDD = 1V was used in our simulations. The total capacitance between the VO2 device and ground was assumed to be U = 200 , while the metallic and insulating resistance of the device was 2kΩ and 100kΩ, respectively. The injection locking capacitance #)% = 1 was used. The resultant operating frequency of the IMT oscillator in our simulation was ( = 100 . Note that with the current experimental setup being dominated by large parasitic capacitances, the experimentally measured frequency is in kHz. Upon building integrated IMT oscillators with on-chip coupling, higher operation frequencies up to 100 MHz or above can be practically achieved (27). The time step of our simulation was fixed at 100 ps used. Each simulation run was performed for an anneal time of M))4M= =10 s which resulted in 1000 oscillation cycles. The simulation was further re-run 100 times to calculate the success probability. We found that the success probability of our Ising solver for finding the Max-Cut is sensitive to parameters like the coupling strength and the stochasticity or noise in the system. As such, we varied the strength of the capacitive coupling and the oscillator jitter to find an optimal coupling value and the oscillator jitter that can maximize the success probability (see supplementary information section S7 for details). For the rest of the scaling simulations reported in Fig 5, we use an optimized coupling capacitance ? = 30 and oscillator jitter of 0.5%.  Fig. 1

. Oscillator based Ising spin. (A)
A schematic of the Ising solver consisting of a network of stochastic phase transition nano-oscillators which are injection locked to an external signal at twice the network frequency to replicate artificial Ising spins. The adjacency matric of the graph to be solved is mapped to the connectivity matrix of the oscillators. The continuous-time dynamics of the network is dictated by an energy that closely resembles the Ising Hamiltonian . As such, the dynamics of the network evolves continuously in time, so as to naturally minimize energy and in the process minimizes the Ising Hamiltonian to reach the ground state. The stochasticity of the IMT oscillators is utilized to escape local minima. (B) Three different injection locking scenarios are shown. No synchronization case refers to free running oscillator with uniform of oscillator phase in the phase space. For first harmonic synchronization (FHIL), the oscillator remains injection locked with input signal at in-phase configuration. For second harmonic synchronization (SHIL), the oscillator phase gets binarized into in-phase and out-of-phase configuration with equal probability. (C) Measured distribution of the oscillator phases for the three scenarios. (D) The equivalent Ising energy for a single oscillator is shown depicting zero, one and two energy minima for stable phase locking, respectively.    The total computational time required to obtain a solution of a particular accuracy with a 99% cumulative success probability. Also shown for comparison, the computation time required by BiqMac solver for running an exact branch-andbound algorithm on a CPU and running a simulated annealing (SA) algorithm on a CPU. (F) Energy consumption for finding solutions with 99% cumulative success probability on Mobius ladder graphs for varying sizes. The energy consumption for BiqMac solver running an exact branch-and-bound algorithm on a CPU and SA algorithm running on a CPU are also estimated for comparison. The VO2 are connected using multi-contact probes. The injection locking signal is applied across the injection locking capacitances #)% using an arbitrary waveform generator. The VDD and the analog gate voltages VGS of the series transistors are applied using a multichannel analog voltage card connected to a computer. The output voltage waveforms of the oscillators were measured using a multichannel digital oscilloscope and analyzed on the computer. The coupling among the oscillators was realized using discrete off-chip coupling capacitances F connected on the breadboard. The perturbation projection vector (PPV) is the phase response of an oscillator to an impulse current input normalized to the amount of injected charge and the oscillator frequency. We setup a SPICE circuit simulation to numerically calculate the PPV. A single IMT nano-oscillator consisting of a 2-terminal hysteretic phase transition device (30) is placed in series with a nchannel MOSFET (Fig. S4.1(a)) to create an oscillating waveform as shown in Fig. S4.1(A). An ideal current source ( ) is connected to the oscillating node for producing a current impulse with a duration much smaller than the oscillator's time period. We sweep the delay of the perturbation over 100 steps covering one time period of the oscillation. For each delay, we compare the altered

S2. Characterization of VO2
waveform (or trajectory) to the original unperturbed trajectory and measure the time difference Δ between the two as shown in Fig. S4.1(B). The corresponding phase difference is measured as = 2 / V4W#*-. The PPV is then calculated as below : (θ) = θ ∫ ( ) Fig. S4.1(C) and (D) show the original waveform of the IMT nano-oscillator and its corresponding PPV. The excess phase generated over 1 oscillation cycle due to the injection locking signal is given by  Energy landscape due to capacitive coupling interaction between two IMT oscillators is similar to interaction energy of the anti-ferromagnetic coupled XY model with a minimum at ( ! − 2 = ). (D) SHIL introduces a local minimum at ( ! − 2 = 0) which was previously a maximum, however the shape of energy landscape around global minimum ( ! − 2 = ) remains the same as before. The XY model is a generalization of the Ising model where each spin is represented not by a discrete variable # ∈ {−1, +1}, but rather by continuous vector # = (sin( # ) , cos ( # )). The interaction energy between two spins in the XY model is given by This becomes equivalent to the Ising model when E # − % F ∈ {0, ± }, i.e. when the spins are completely parallel or antiparallel with each other. As shown in Fig. S4.3(B), the two-oscillator interaction energy landscape closely resembles the XY-model especially around the neighbourhood of the minima or maxima points, i.e. ( ! − 2 ) ∈ {0, ± } and become equal to the Ising model at E # − % F ∈ {0, ± }. Upon introducing the SHIL, the energy landscape is changed. As shown in figure Fig. S4.3(d), this introduces a local minimum at ! − 2 = 0. Hence, the shape of the energy landscape deviates The energy landscape of the oscillator network is sensitive to the choice of simulation parameters such as capacitive coupling and stochasticity in the system. In our simulations, we set the parameters such as U , #)% , Y and # as described in the Materials and Method section of the main text to achieve an operating frequency of ( = 100 . The amplitude of the injection locking signal is chosen to be sufficient enough to reliably binarize the oscillator phases. The interaction energy between the oscillators is modulated by the value of the coupling capacitance ? and we find the optimum value for this interaction by measuring the success probability of a 10 node mobius ladder Max-cut problem for different values of capacitance as shown in Fig.  S7(A). Lower coupling strength means that the interactions between the oscillators are too weak, resulting in the energy difference between the different Ising spin configurations being smaller, leading to a lower chance of the oscillators to find the minimum energy configuration and hence a lower success probability. Conversely, if the interaction energy is too strong compared to the injection locking strength, then the binarization of the oscillator phases is reduced resulting in a wider distribution in the phase space. This apparent destruction of the 'up' and 'down' spin nature of the phases in some simulation leads to a lower overall success probability. We found that a coupling capacitance strength of ? = 30 was an optimum value to achieve maximum success probability and is approximately 16% of the total VO2 capacitance U . We used this coupling capacitance value for all our simulations up to 200 oscillators. Stochasticity is the second parameter that can greatly affect the likelihood of the oscillator network to escape the local energy minima and reach the global optimum energy configuration. We introduce stochasticity in the form or oscillator jitter phase (or frequency) jitter. Too little stochasticity would mean the system gets stuck in local minima and is not able to jump out to reach the global minimum. On the other hand, a very large amount of intrinsic noise would cause the system to jump out of the global energy minimum configuration even at highest injection locking leading to lower probability of success when the phase is read-out. We varied the amount of noise for a 10-oscillator system and found that a 0.5% oscillator jitter was the ideal amount of noise to achieve maximum success probability while solving a mobius ladder Max-Cut graph (Fig S7(B)). We used this amount of oscillator jitter for all our simulations up to 800 oscillators. Since for larger graphs the success probability decreases exponentially, we need to re-run the problem several times (either serially on same solver or in parallel over multiple solvers) for ensuring 99% probability. The number of iterations required to obtain at least one solution with % accuracy and with a success probability of 99% is given by [log (1 − 0.99)] ologE1 − E+FF4EE,N% Fp ⁄ , where E+FF4EE,N% is the success probability. We see that the number of iterations required to obtain 99% success probability over multiple iterations on the same solver increases exponentially. Fig. S9. Comparison between LC oscillator, ring oscillator and IMT oscillator. Fig. S9 shows the circuit schematic of an LC oscillator (32), a seven-stage ring oscillator (25) and the phase-transition nano-oscillator. The inductor in an LC oscillator takes up most of the area, close to 0.1 mm 2 (40). The ring oscillator has a comparatively lower area. Ref. (25) reports a unit cell area of 950 2 . In contrast, our IMT nano-oscillator comprises a transistor and a scaled VO2 device which takes up orders of magnitude lower area.