DROID: discrete-time simulation for ring-oscillator-based Ising design

TL;DR

DROID is a fast and accurate simulation method for Ising machines, which solve complex problems using oscillator networks.

Contribution

DROID introduces an event-driven simulation technique for CMOS Ising machines that is significantly faster than traditional methods.

Findings

DROID is nearly four orders of magnitude faster than HSPICE simulations for oscillator arrays.

It achieves two orders of magnitude speedup over commercial fast SPICE solvers.

DROID produces solution distributions similar to actual hardware.

Abstract

Many combinatorial problems can be mapped to Ising machines, i.e., networks of coupled oscillators that settle to a minimum-energy ground state, from which the problem solution is inferred. This work proposes DROID, a novel event-driven method for simulating the evolution of a CMOS Ising machine to its ground state. The approach is accurate under general delay-phase relations that include the effects of the transistor nonlinearities and is computationally efficient. On a realistic-size all-to-all coupled ring oscillator array, DROID is nearly four orders of magnitude faster than a traditional HSPICE simulation and two orders of magnitude faster than a commercial fast SPICE solver in predicting the evolution of a coupled oscillator system and is demonstrated to attain a similar distribution of solutions as the hardware.

Figures

Period waveforms for a 5 $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document}$ 5 A2A array from (a) DROID, (b) HSPICE, and (c) GenAdler for the same initial conditions of the ROs.

Distribution of DROID runtimes for various A2A array sizes.

The evolution of RO-phases from random to binarized phases for a particular problem; phases are represented by the color bar at right.

Distributions of the Hamiltonian from 100 samples of a problem of density 0.2 obtained from (a) DROID and (b) the hardware. The X-axis has been normalized to the best Hamiltonian value obtained for the problem.

Three five-stage ROs, with a positive and a negative coupling. The green stage is the reference, and odd (even) stages are shown in red (blue).

The relative arrival times on nets *a* (green waveform) and *b* (blue waveform), top, impact the transition delay on net *c* (middle). The transition delay as a function of phase difference, and the interaction window *W* are illustrated at the bottom.

An illustration of the A2A concept through a small three-RO structure, showing shorting cells *S* on the diagonal, programmable off-diagonal coupling cells *C*, and enable cells *EN*.

The handling of an event on *net_a* is influenced by the knowledge of events on nets in some neighborhood around it, as shown on the top. CASE 1 shows the scenario when an event on the other input (*net_b*) of the same instance is known. CASE 2 shows the scenario where LOOK_BACK is invoked to find an event on *net_c* that can cause an event on *net_b* that might lie within the interaction window of the event at *net_a*.

Period waveforms for a 5 $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document}$ 5 A2A array from (a) DROID, (b) HSPICE, and (c) GenAdler for the same initial conditions of the ROs.

The phase differences across various edges are shown above. The phase differences $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{AB}$$\end{document}$ (red) and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{BC}$$\end{d

A simulator step with the sorted queue Q and the Net2Event map: PROCESS_EVENT consumes one or more events and generates future events.

The handling of an event on *net_a* is influenced by the knowledge of events on nets in some neighborhood around it, as shown on the top. CASE 1 shows the scenario when an event on the other input (*net_b*) of the same instance is known. CASE 2 shows the scenario where LOOK_BACK is invoked to find an event on *net_c* that can cause an event on *net_b* that might lie within the interaction window of the event at *net_a*.

Full text

Introduction

A critical computational domain for hardware accelerators is the area of solving combinatorial optimization problems (COPs) that are NP-complete or NP-hard—e.g., the traveling salesman, satisfiability, and knapsack problems. Today, such problems are solved on classical computers using heuristics with no optimality guarantees, or approximation algorithms with loose optimality bounds.

Ising computation is a promising emerging computational model for solving COPs. Ising machines are inspired by the work of Ernst Ising, who proposed a formulation based on binary states called spins, with allowable values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1$$\end{document} , to explain ferromagnetism. Such systems have a natural tendency to find a ground state with a configuration of spins that minimizes energy. Ising computation maps discrete combinatorial optimization problems to this paradigm. Under a linear transformation, Boolean quadratic unconstrained binary optimization (QUBO) problems can be formulated as two-body Ising interactions. Karp’s list of 21 NP-complete problems are shown to have an Ising formulation^1^, and many other problems can also be formulated in this way.

Recently, there has been great interest in building Ising hardware accelerators, realizing spins using superconducting loops in D-Wave machines^2,3^; optical parametric oscillators in coherent Ising machines^4,5^, ring oscillators (ROs) in CMOS-based Ising machines^6–8^, voltage-controlled nanomagnet stacks^9^, and memory cells in SRAM-based engines^10,11^. Oscillator-based methods use the phenomenon of synchronization, whereby a system of coupled oscillators with similar frequencies, converge to a common frequency and fixed phase difference through injection locking. The dynamics of coupled oscillators have been studied as early as 1663, when Huygens noticed the synchronization of pendulums connected to a common bar. Several works simulate synchronization behavior through differential equations: Adler derived a closed-form expression for locking in LC oscillators^12^; Winfree explored weak interactions of periodic behavior in biological rhythms; Kuramoto^13^ studied chemical oscillations under sinusoidal interactions. We study injection locking in an analog system that uses CMOS ROs.

Unlike platforms using exotic futuristic technologies, CMOS RO-based Ising machines use a mainstream semiconductor technology that is scalable, compact, economically and reliably mass-manufacturable today, and can operate at room temperature instead of requiring expensive high-power mK-level refrigeration schemes. The synchronization of RO-based Ising machines can be simulated using SPICE, but this is computationally intensive and does not scale well. Simulators for oscillator-based Ising machines are based on analytical solutions to the generalized Adler equation^14^ and the generalized Kuramoto equation^6^. A prior event-driven approach^15^ fast-forwards through multiple RO cycles until the phase difference between some pair of ROs crosses an integer multiple of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} : this is considered to be an event. If the phase difference remains close to an integer multiple of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} for some iterations, the associated coupling is removed from the system and a phase merging scheme is used to lock the phases of these oscillators henceforth. A hardware realization of a generalized Kuramoto equation solver has also been demonstrated^16^.

Prior methods have several limitations. First, they represent the phase of each oscillator by the phase of a single reference stage. However, the phase differences at specific coupling sites between two oscillators may differ from the differences in their reference phases. Second, methods that use phase merging^15^ can be misleading: the phase of an RO can diverge even after it appears to come close to another RO phase. This work proposes DROID (Discrete-Time Simulation for Ring-Oscillator-Based Ising Design), a method for simulating RO-based Ising machines, that overcomes the above limitations. Its contributions are as follows:

• We show that for coupled RO systems, prior continuous-time (CT) simulation abstractions, such as the generalized Adler formulation^14^, are abstractions of a discrete-event simulation, operating under restrictive assumptions that allow closed-form solutions, including assumptions of infinitesimal changes (section “Discrete-time vs. continuous-time simulation of coupled-RO systems”). Our approach removes these restrictions and uses lookup-table-based functions, characterized using HSPICE.
• Unlike prior methods that work in the continuous domain, we develop a discrete-time event-driven simulation methodology (section “Simulating the A2A array”) to predict the behavior of coupled RO systems; this method is inspired by timing analysis methods that are widely used for digital circuits, which achieve acceptable accuracy at a fraction of the runtime of HSPICE. Our approach is event-driven, where an event is defined with fine granularity, associated with a coupled transition between two oscillators.
• Our approach is 125 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} –7441 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} faster than HSPICE and 11 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} –109 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} faster than a commercial fast SPICE solver at similar accuracy, with larger speedups for larger systems. We match the distribution of our solutions, across 250 problems of various oscillator coupling densities, 100 samples per problem, and multiple initial conditions, against a CMOS RO-based Ising hardware solver^8^, and show that the distance between distributions, is small. The paper is organized as follows. Section “CMOS-based coupled-oscillator systems” summarizes the concepts that guide this work. Section “A silicon-proven all-to-all coupled-RO Ising machine” describes an all-to-all-connected Ising hardware accelerator that serves as our hardware testcase. Sections  “Discrete-time vs. continuous-time simulation of coupled-RO systems” and “Limitations of the continuous-time approximation” then analyze the relationship between discrete-time and continuous-time simulation of coupled oscillator systems. We describe our event-driven simulation scheme for the Ising hardware in section “Simulating the A2A array” and show our simulator results in section “Results and discussion”, finally concluding the paper in section “Conclusion”.

CMOS-based coupled-oscillator systems

The Ising model

The Ising formulation of a COP minimizes the following objective function, referred to as a Hamiltonian:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} H(\textbf{s}) = \textstyle -\sum _{i = 1}^{N} \sum _{j=1}^{N} J_{ij} s_i s_j - \sum _{i=1}^{N} h_i s_i, \end{aligned}$$\end{document}

In the magnetics domain, this models the energy of a system of N spins; spin is an intrinsic property associated with a subatomic particle, atom, or molecule, and can take on a value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+1$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1$$\end{document} . The Hamiltonian is the energy of a system of spins as a function of their interactions ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij} s_i s_j$$\end{document} ) and the effect of external magnetic fields on individual spins ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_i s_i$$\end{document} ). A physical Ising machine settles to a ground state of low-energy states favored by nature, thus minimizing the Hamiltonian. Therefore, by suitably mapping a COP to the weights \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_i$$\end{document} , an Ising machine can solve a COP formulated as a Hamiltonian.

Two spins \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_j$$\end{document} are in-phase if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i = s_j$$\end{document} , and out-of-phase otherwise. From (1), a positive [negative] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij}$$\end{document} encourages \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_j$$\end{document} to be in-phase [out-of-phase], and a positive [negative] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_i$$\end{document} pushes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i$$\end{document} to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+1$$\end{document} [ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1$$\end{document} ].

CMOS-ring-oscillator-based Ising machines

Core principle. Figure 1 shows a CMOS-based Ising machine with three coupled ROs. Each RO has an identical number of stages; each stage in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i\textrm{th}$$\end{document} oscillator, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} , is an inverter, except for one NAND stage with an enable signal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {en}_i$$\end{document} , to start the oscillator.

We denote the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k^\textrm{th}$$\end{document} stage of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_{i}^{k}$$\end{document} , with a phase, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _i^k$$\end{document} , given by the arrival time of the rising edge at the stage output. The phase \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _i$$\end{document} of an oscillator is defined as the phase, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _i^0$$\end{document} , at the output of its zeroth stage (shown in green in the figure). The phase at stage \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i^k$$\end{document} is found by adding to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _i$$\end{document} the sum of the stage delays from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i^0$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i^k$$\end{document} . The time between two consecutive rising edges at the output of a stage is the period of the RO. We denote the nominal period, i.e., the period of each uncoupled, freely oscillating RO, as T. When the oscillators are coupled together, they synchronize to the same period, which may be different from T.

We designate one of the oscillators as a reference oscillator, always setting its spin to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+1$$\end{document} ; without loss of generality, we refer to this as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_0$$\end{document} . To assign a spin value to an oscillator, its phase is compared with that of the reference oscillator: an oscillator that is in-phase with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_0$$\end{document} is said to have a spin of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+1$$\end{document} , and one that is out-of-phase has a spin of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1$$\end{document} .

Ising machines employ weak coupling^13^, where the delay change due to coupling is small compared to the nominal period T. A non-zero coupling coefficient, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij}$$\end{document} , in the Ising model is realized by coupling \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_j$$\end{document} . One way to couple a pair of ROs is by connecting the outputs of two corresponding stages in each RO by a resistor; other coupling schemes may also be used^17^. For resistive coupling, the coupling strength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij}$$\end{document} is determined by the values of the resistors between inverters in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_j$$\end{document} ; section “A silicon-proven all-to-all coupled-RO Ising machine” will describe a circuit that implements multiple programmable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij}$$\end{document} values.

We refer to inverters that are at an even-parity [odd-parity] distance from the reference stage as even [odd] stages, shown in blue [red] in Fig. 1. Coupling between two same-parity stages in different ROs is referred to as positive coupling ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_1$$\end{document} in Fig. 1), while coupling between opposite-parity stages is termed negative coupling ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_2$$\end{document} in Fig. 1). Positive coupling encourages the ROs to be in-phase, while negative coupling encourages the ROs to be out-of-phase. In Fig. 1, the net a driven by stage \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_0^1$$\end{document} is coupled to net b at the output of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_1^1$$\end{document} by a resistor; this causes the stage delays to change from their nominal values:

• When a rises, if b is low (i.e., yet to rise), it opposes the rise transition and the delay of a is increased by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _1^a$$\end{document} relative to the uncoupled case. On the other hand, if b is high (i.e., it has already risen), it aids the rise and reduces the delay of a by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _2^a$$\end{document} . In both cases, the rising edge of a is brought closer to the rising edge of b, reducing the phase difference between the signals.
• Similarly, when a falls, if b is high (i.e., yet to fall), the delay is increased by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _3^a$$\end{document} ; if b is low (i.e., has already fallen), its delay reduces by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _4^a$$\end{document} . In each case, the falling edges of a and b are brought closer, bringing the oscillators closer to phase-locking.

Net b behaves analogously, with delay shifts of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _1^b, \cdots , \delta _4^b$$\end{document} .

Consider the coupled oscillator system of Fig. 1, but with the coupling between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_2$$\end{document} removed, i.e., only \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_0^1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_1^1$$\end{document} are coupled. Figure 2a shows the effect of coupling over multiple cycles, where the waveforms at a and b begin with a phase difference, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{12}$$\end{document} . In the uncoupled case, this phase difference is unchanged, but under coupling, as a slows down and b speeds up, the phase difference decreases and the edges align. The inset in Fig. 2b shows the first cycle, i.e., the green region in Fig. 2a. Due to coupling, the first rising edge of a is delayed by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _1^a$$\end{document} while that of b arrives earlier by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _2^b$$\end{document} . Similarly, the falling edge at a is delayed by an additional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _3^a$$\end{document} while that at b is sped up by an additional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _4^b$$\end{document} . Thus, the phase difference is reduced by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\delta _1^a + \delta _2^b + \delta _3^a + \delta _4^b)$$\end{document} after the first cycle. Subsequently, as long as transitions in a are completed before those in b, a continues to be delayed by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _1^a+\delta _3^a$$\end{document} and b continues to speed up by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _2^b+\delta _4^b$$\end{document} every cycle, and the phase difference at the end of k cycles reduces to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{12} - k(\delta _1^a + \delta _2^b + \delta _3^a + \delta _4^b)$$\end{document} until the phases align.

Coupled oscillator systems. We abstract a general coupled oscillator system with a graph \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G=(V,E)$$\end{document} , where the vertex set V is the set of oscillators, and each element \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_{ij}$$\end{document} of the edge set E corresponds to a pair of oscillators \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_j$$\end{document} with an edge weight corresponding to the coupling strength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij}$$\end{document} . We denote the change in delay of the coupled stage of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} in each cycle, under a phase difference of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij}$$\end{document} , by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{J_{ij}}(\phi _{ij})$$\end{document} ; for the example in Fig. 2b, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{J_{12}}(\phi _{12}) = (\delta _1^a + \delta _3^a)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{J_{12}}(\phi _{21}) = (\delta _2^b + \delta _4^b)$$\end{document} . The net change in the period of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k^\textrm{th}$$\end{document} cycle, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathscr {D}_i^k$$\end{document} , is the sum of changes in the delay of each stage, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathscr {D}_i^k = \textstyle \sum _{(i,j) \in E} f_{J_{ij}}(\phi _{ij}^{k})$$\end{document} .

Synchronization. A pair of coupled oscillators that have the same period will have a constant phase difference. Conversely, when all phase differences are constant, it implies that all pairs of coupled oscillators have the same period, a phenomenon referred to as synchronization. In Fig. 2a, since the coupling is purely positive, all phases align upon synchronization. As a result, all stage delays go back to their nominal uncoupled values and the period reverts to the nominal period, T.

Practical considerations for CMOS-based coupled oscillator systems

Delay change as a function of coupling location At a given coupling location between two ROs, the delays from the reference stage of each RO to their coupling stages may be different. As a result, the phase difference at the coupling site may be different from the phase difference of the reference ROs. For example, in Fig. 1, the path from the reference stage to the coupling site for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{12}$$\end{document} involves two inverter delays in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_1$$\end{document} , but three in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_2$$\end{document} . Therefore, the phase difference at this coupling site is not the same as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{12}$$\end{document} , the phase difference between their reference stages. Prior simulators have not considered this issue.

Moreover, the stage delay depends on the presence or absence of coupling: ROs with more couplings may have different delays to a coupling site than ROs with fewer couplings. Thus, as the phase changes along an RO, delay shifts from earlier stages affect the phase difference at later stages, and using an identical \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{i}$$\end{document} at all coupling stages leads to inaccuracies. Many events within a cycle interact subtly to produce a total delay shift for each RO. Such interactions within a cycle are considered in our framework.

Delay change as a function of arrival time difference The speedup or slowdown in stage delay due to coupling depends on the arrival time differences between the signals on the coupled nets. Figure 3 illustrates this trend at five different relative signal arrival times, (P1, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cdots$$\end{document} , P5), on nets a (blue) and b (green) in Fig. 1. Near P1 and P5, the difference in arrival times is large enough that the transition on a does not overlap with a transition on b, and the change in delays is constant as a is effectively stable throughout the rise of b. At P2 and P4, when the edges are closer, the opposing or assisting transistor in the other RO is no longer completely on and it sees a reduced gate-source voltage, which reduces its effect on the delay. If two rising edges with identical transition times are exactly aligned, at P3, no current flows through the resistor, and the transitions do not affect each other. Thus, we see that there is a window around each edge beyond which the arrival of another edge causes a constant stage delay shift, but within this window, the delay shift is a function of the phase difference. This window of width W is called the interaction window and extends on either side of a transition as seen by the highlighted orange box in Fig. 3.

The delay at the output of a stage is also a function of the transition time of the signal at its input^18^. The transition time at net b is also seen to show a similar trend as shown in the bottom panel of Fig. 3, which has an effect on the delays of subsequent stages.

A silicon-proven all-to-all coupled-RO Ising machine

The design of a coupled-RO Ising machine is particularly tricky because of the need to ensure the uniformity of coupling weights between each pair of ROs. The shift in each transition time in the RO depends on RC parasitics associated with the coupling mechanism, which depend on the precise layout. It has been shown that through regular matched layouts, an all-to-all (A2A) coupled RO array with uniform coupling coefficients can be designed^8^. This A2A design is silicon-proven and will be used in our evaluations. Note that other designs with planar (hexagonal^19^ and King’s graph^7^) coupling have been proposed, but we focus on the A2A testcase because of its compactness and greater flexibility: in particular, an A2A design with N coupled ROs is equivalent to a planar hexagonal/King’s graph array with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim N^2$$\end{document} coupled ROs^20,21^ due to the need for planar arrays to replicate spins during minor embedding^8^. This family of A2A arrays has been applied to solve problems ranging from max-cut^8^ to maximal independent set^17^ to Boolean satisfiability^22^.

For illustration, a simplified schematic, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij} \in \{-1,0,+1\}$$\end{document} , for a three-RO system is presented in Fig. 4. Each RO is a combination of a vertical and a horizontal RO that are strongly coupled so as to implement the same spin, and have enable cells EN (similar to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$RO_i^0$$\end{document} in Fig. 1) outside the array. Strong couplings are implemented as shorts (upper inset) in the shorting cells S, placed at the diagonals (i, i). Each coupling cell C at off-diagonal location (i, j) has programmable coupling between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_j$$\end{document} . A coupling cell has two switches (lower inset) to enable either a positive or a negative coupling. The simplified figure shows possible couplings of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{-1, 0, +1\},$$\end{document} but coupling of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm 7$$\end{document} is demonstrated in silicon^8^. Using couplings at both (i, j) and (j, i), each programmable to integer values up to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm C_{max}$$\end{document} , the coupling coefficient of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i s_j$$\end{document} can implement integer \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij} \in [-2C_{max}, +2C_{max}]$$\end{document} .

Discrete-time vs. continuous-time simulation of coupled-RO systems

Traditional continuous-time formulations. The behavior of an LC oscillator, when injected with a sinusoidal signal, was described by Adler’s equation^12^; a slightly different equation is used by Kuramoto^13^. To extend this beyond a single coupling and sinusoidal signals, the generalized Adler (GenAdler) equation for a network of N coupled oscillators was shown^14^ to have the form:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{d \phi _{i}(t)}{dt} = (\omega _i - \omega ^*) + \omega _i \textstyle \sum _{j=1,j \ne i}^{N}c_{ij} (\phi _{ij} (t)) \end{aligned}$$\end{document}

Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij}(t) = \phi _{i}(t) - \phi _{j}(t)$$\end{document} is the difference between the phases of oscillators i and j, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _i$$\end{document} is the frequency of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i^\textrm{th}$$\end{document} oscillator, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega ^*$$\end{document} is the central frequency of the network, and for oscillators i and j, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_{ij}(.)$$\end{document} is a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\pi$$\end{document} -periodic function that represents the coupling-induced delay shift in each RO cycle. Prior methods^14^ abstract \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_{ij}$$\end{document} as a well-behaved function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij}$$\end{document} , the phase difference of the coupled ring oscillators; Fig. 5a shows an example HSPICE-characterized function showing the RO period shift against the phase difference.

Discrete-time formulation. The continuous formulation models a discrete-event system composed of a sequence of coupling events. We now examine the limitations of the continuous-time GenAdler model, as well as those of the phase delay shift model.

Figure 5a shows a model for the delay shift, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{J_{ij}}(\phi _{ij})$$\end{document} of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_j$$\end{document} at an example coupling strength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij} = 1$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij}$$\end{document} is the phase difference between the two ROs. The delay shifts of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_j$$\end{document} are shown by the red and green dotted lines, and the relative phase shift, i.e., the difference between these delay shifts, is shown by the solid blue line. It is seen that when the edges align to be in-phase or out-of-phase (i.e., at 0 or T/2), the relative phase shift is zero.

We model the coupled system using a sequence of discrete events that add up to a shift in an RO clock period at the end of a cycle. We denote the phase, period, and frequency as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{i}^{k}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{i}^{k}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{i}^{k}$$\end{document} , respectively, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} at the end of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k^\textrm{th}$$\end{document} cycle. We use a datum oscillator frequency, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega ^*$$\end{document} , which may be the frequency corresponding to the period T of each uncoupled oscillator.

The phase shift of each oscillator from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{i}^{k}$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{i}^{k+1}$$\end{document} during the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k+1)\textrm{th}$$\end{document} cycle is caused by two factors:

1. Frequency drift with respect to the reference oscillator:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta \phi _{1i}^{k+1} = \phi _{1i}^{k+1} - \phi _{1i}^{k} = (\omega _{i}^{k} - \omega ^*) T_{i}^{k} \end{aligned}$$\end{document}

2. Phase/frequency drift due to coupling to other ROs:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta \phi _{2i}^{k+1} = \phi _{2i}^{k+1} - \phi _{2i}^{k} = \textstyle \sum _{(i,j)\in E} f_{J_{ij}} (\phi _{ij}^{k}) \end{aligned}$$\end{document}

where E is the set of edges in the coupling graph (section “CMOS-ring-oscillator-based Ising machines”). The net phase shift in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k+1)^\textrm{th}$$\end{document} cycle is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta \phi _{i}^{k+1} =\Delta \phi _{1i}^{k+1} + \Delta \phi _{2i}^{k+1} =(\omega _{i}^{k} - \omega ^*) T_{i}^{k} + \textstyle \sum _{(i,j)\in E} f_{J_{ij}} (\phi _{ij}^{k}) \end{aligned}$$\end{document}

At synchronization, the clock frequency, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _i^k$$\end{document} is the same for all oscillators. Thus, in writing the relative phase difference between coupled oscillators \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_j$$\end{document} (note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ji} = -\phi _{ij}$$\end{document} ), their corresponding first terms in (5) cancel, and we have:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta \phi _{i}^{k+1} - \Delta \phi _{j}^{k+1} = \textstyle \sum _{(i,j)\in E} f_{J_{ij}} (\phi _{ij}^{k}) - f_{J_{ij}}(-\phi _{ij}^{k}) = 0 \end{aligned}$$\end{document}

The last equality arises because when the phases are locked, the difference between the delay shifts of locked oscillators is zero.

Relation between the continuous- and discrete-time formulations. Unlike coupled sinusoidal oscillators, as modeled by the GenAdler formulation, where the signal value changes throughout a cycle, for an RO system, the coupling component of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\phi _{i}(t)/dt$$\end{document} changes during signal transitions but not when signals are stable at logic 0 or logic 1. Under infinitesimal phase changes per cycle, the derivative can be approximated if the net phase change over m cycles is small. If the period of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} in cycle k is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{i}^{k},$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{d\phi _{i}(t)}{dt} \approx \frac{1}{m} \frac{\sum _{k=1}^m \Delta \phi _{i}^{k}}{\Delta t}, \hbox { where } \Delta t = \textstyle \sum _{k=1}^m T_{i}^{k} \end{aligned}$$\end{document}

From (5), the phase change in time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t$$\end{document} (m cycles of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_i$$\end{document} ) is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta \phi _{i} = \textstyle \sum _{k=1}^{m} \Delta \phi _{i}^{k} =\textstyle \sum _{k=1}^m \left( \omega _{i}^{k} - \omega ^* \right) T_i^k + \textstyle \sum _{(i,j)\in E} \textstyle \sum _{k=1}^{m} f_{J_{ij}} (\phi _{ij}^{k}) \end{aligned}$$\end{document}

The delay shifts over the m cycles must be assumed to be small, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{i}^{k} \approx \omega _i \; \forall \; k = 1, \cdots , m$$\end{document} . Under this assumption,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sum _{k=1}^{m} (\omega _{i}^{k} T_{i}^{k} - \omega ^*T_{i}^{k}) \approx (\omega _{i} - \omega ^*) \sum _{k=1}^{m} T_{i}^{k} = (\omega _{i} - \omega ^*) \Delta t \end{aligned}$$\end{document}

Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_i$$\end{document} be the average value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{i}^{k}$$\end{document} . Then

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} T_i = \frac{1}{m} {\textstyle \sum _{k=1}^{m} T_{i}^{k}}= \frac{\Delta t}{m} \Rightarrow m = \frac{\Delta t}{T_i} = \left( \frac{\omega _i}{2\pi } \right) \Delta t \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _i {\mathop {=}\limits ^{\Delta }} 2\pi /T_i$$\end{document} . Under small phase changes over m cycles, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l = 1, 2, \cdots , k$$\end{document} , we write \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{J_{ij}} (\phi _{ij}^{l}) \approx f_{J_{ij}} (\phi _{ij}^{1}) + S_{ij}(\phi _{ij}^{l}-\phi _{ij}^{1})$$\end{document} as a linear approximation, where the sensitivity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{ij} = \left. \left( \partial f_{J_{ij}}/\partial \phi _{ij} \right) \right| _{\phi _{ij}^{1}}$$\end{document} . Therefore,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textstyle \sum _{k=1}^{m} f_{J_{ij}} (\phi _{ij}^{l}) = m \left[ f_{J_{ij}} (\phi _{ij}^{1}) + S_{ij} \left( \frac{\textstyle \sum _{k=1}^{m}\phi _{ij}^{k}}{m} - \phi _{ij}^{1} \right) \right] = \omega _i \left[ \frac{f_{J_{ij}} (\phi _{ij})}{2\pi } \right] \Delta t \hspace{4mm} (\hbox {using}\, (10)) \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij} = \sum _{k=1}^{m}\phi _{ij}^{k}/m$$\end{document} is the mean value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij}^{k}$$\end{document} over m cycles. From (8), (9), and (11), since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t$$\end{document} is assumed to be small,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial \phi _i}{\partial t}&\approx \frac{\Delta \phi _i }{\Delta t} = (\omega _{i} - \omega ^*) + \omega _i \sum _{(i,j)\in E}\frac{f_{J_{ij}} (\phi _{ij})}{2\pi } \end{aligned}$$\end{document}

Setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_{ij} = f_{J_{ij}} (\phi _{ij})/(2\pi )$$\end{document} , this is the GenAdler equation, (2).

Limitations of the continuous-time approximation

The continuous-time approach is effective in matching a coupled CMOS RO system when:

• the phase difference between any pair of oscillators is independent of the location of coupling, and
• the phase difference between any pair of oscillators within a cycle is independent of their coupling to other oscillators.

From section “Discrete-time vs. continuous-time simulation of coupled-RO systems”, the function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_{ij}(.)$$\end{document} is crucial to the correctness of the model. For CMOS ROs, coupling is expressed through complex MOS models (e.g., BSIM4, BSIM-CMG), the mapping from the system to the coefficient is nontrivial.

To describe the impact of assumption (a), we present an example that shows that the phase difference at a coupling cell in an A2A array depends on the location of the cell within the array and the magnitude of coupling in previous stages of the ROs. Our example in Fig. 6a shows a section of an A2A array, where the labeled nets are inputs to a set of coupled cells in the array. The horizontal RO, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_Z$$\end{document} , runs through the lines z, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z'$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z''$$\end{document} . This oscillator couples with the vertical oscillators, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_W$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_X$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_Y$$\end{document} , at the stages with inputs w, x, and y, respectively, through tiles that implement the coupling coefficients, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ZW} = 6$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ZX} = 6$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ZY} = -6$$\end{document} . The locations w, x, y, and z correspond to reference phases of respective ROs. Given a set of initial reference phases for each RO, we show the transitions at various locations in Fig. 6b. As mentioned in section “CMOS-ring-oscillator-based Ising machines”, the phase difference at the coupling cell sites may differ from the phase difference at the reference stages of their oscillators. For instance, the phase difference at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ZX}$$\end{document} is 2.6 ps, corresponding to the difference in arrival times of transitions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z'$$\end{document} and x, which is different from the difference or 107.4 ps in the arrival times of transitions at z and x. As can be seen from Fig. 3, an inaccuracy of this magnitude (comparable to window W) will impact delay calculation for a coupling cell, also affecting the phase difference at the next cell.

Next, we consider the impact of assumption (b) alone, using corrected arrival times to eliminate the contribution of errors from assumption (a). As discussed in section “Practical considerations for CMOS-based coupled oscillator systems”, RO delays can vary with the arrival time difference, and these delays are also subtly impacted by changes in the signal transition time at the RO input. We determine the coupling between oscillators in Fig. 6c for a case where the arrival times of transitions at w, x, and y are identical, but the transition at z is slightly delayed, compared to the corresponding values in Fig. 6b. This small change has the effect of changing the coupling delay and transition time at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z'$$\end{document} , with a ripple effect on the timing of transitions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z'$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z''$$\end{document} , caused by the coupling between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_Z$$\end{document} and other oscillators: it can be seen that the small shift of < 20 ps at z shifts the transition at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z'$$\end{document} by > 50 ps, and at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z''$$\end{document} by > 100 ps. These magnified shifts arise because the arrival time at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z'$$\end{document} depends on the magnitude of coupling, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ZW}$$\end{document} , and that at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z''$$\end{document} depends on both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ZW}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ZY}$$\end{document} . Thus, it is not just the number of stages from the reference to the coupling stage that affects the delay shift; the magnitude of coupling in previous stages, and the precise timing relationship between waveforms in those stages, also affect phase differences at a later stage.

The GenAdler formulation in (2) makes assumptions (a) and (b), and simply uses the phase difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij}$$\end{document} between the reference stages of oscillators i and j. In section “Simulating the A2A array”, we present a fine-grained event-driven approach to overcome these limitations.

We know of only one prior event-driven approach^15^, but it uses a fundamentally different definition of events from ours, and speeds up the generalized Kuramoto simulation. This method inherits the assumptions as GenAdler, and hence, limitations (a) and (b). Its speedup mechanism determines whether, at any time, two or more ROs achieve a phase difference that is an integral multiple of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} radians: if so, it assumes that these oscillators will remain permanently phase-locked from that time onwards. Through this assumption, the number of variables is reduced as the simulation proceeds, reducing its computational cost.

We show a counterexample, based on HSPICE simulations, that illustrates that such preliminary phase-locking assumptions^15^ can be incorrect. Consider a system comprising four ROs, denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_A$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_B$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_C,$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_D$$\end{document} . Figure 7 shows the phase differences in radians between coupled RO pairs as the system evolves in time. The phase difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{AB}$$\end{document} (red) remains close to 0 from 20 ns to 400 ns and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{BC}$$\end{document} (blue) remains close to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\pi$$\end{document} from 20 ns to 400 ns, as shown in the highlighted green box. In the interval [20 ns, 400 ns], it appears as if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_A$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_B$$\end{document} are locked in-phase while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_B$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_C$$\end{document} are locked out-of-phase, but oscillator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_D$$\end{document} is not yet settled as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{AD}$$\end{document} (green) continues to change. As shown in the example in Fig. 6, couplings in earlier stages affect delays within the same cycle, and it is a result of this effect that the changes in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{AD}$$\end{document} become more dramatic around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=$$\end{document} 400 ns, it causes phase differences at other coupling cells to change. The net effect of these changes on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_A$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_B$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_C$$\end{document} is that they leave their seemingly phase-locked relationships as the system evolves, and settle to a different equilibrium at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=$$\end{document} 900 ns, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_A$$\end{document} is out-of-phase with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_B$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_D$$\end{document} , and in-phase with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_C$$\end{document} . If the phases of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_A$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_B$$\end{document} (or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_B$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_C$$\end{document} ) were merged into a single phase based on their behavior between 20 and 400 ns, this equilibrium stage would not be captured by the simulation.

Simulating the A2A array

Capturing the timing information

We define an event as a rise or fall transition in a digital signal at the input of a logic gate in an RO, which can cause the opposite transition at its output. We characterize each transition by its arrival time, transition time (rise/fall time), and whether the signal is rising or falling. Our approach is motivated by timing analysis in CMOS digital design, where a timing arc is used to propagate an event at the input of a gate to an event at its output. Cell timing information, i.e., the input-to-output delay and the output transition time, is captured in lookup tables as functions of the output load and the transition time of the input signal.

The invariant during timing analysis is the computation of the arrival time and transition time at a node. Given an event at the input of a cell, characterized by these two values, the timing information of the cell can be used to generate an output event(s) and their arrival time(s) and transition time(s). These events are expressed at the input of another cell, and the process continues.

As mentioned in section “A silicon-proven all-to-all coupled-RO Ising machine”, the A2A array has three types of cells: enable, coupling, and shorting. The simulator works with a timing view of these cells, and in the remainder of this section, we discuss this timing abstraction, using the notation in Fig. 4.

The enable cell has two inputs, one enable signal, enable, and another from the RO itself (inp). The cell is modeled by a timing arc from inp to outp. Since the load is the same for all enable cells, a one-dimensional table, characterized using HSPICE, is used to represent the cell rise delay as a function of the input transition time; a similar table characterizes the fall delay.

The coupling cell has four inputs ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{in}^{f}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{in}^{f}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{in}^{r}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{in}^{r}$$\end{document} ) and four outputs ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{out}^{f}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{out}^{f}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{out}^{r}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{out}^{r}$$\end{document} ), where the symbols \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v$$\end{document} correspond to the horizontal and vertical ROs, and the superscripts f and r represent the forward and reverse path, respectively, through the cell. The horizontal timing arc ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{in}^{f}$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{out}^{f}$$\end{document} ) of the horizontal RO interacts with the vertical timing arc ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{in}^{f}$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{out}^{f}$$\end{document} ) of the vertical RO within a window when the cell implements a non-zero coupling coefficient; otherwise the horizontal and vertical paths through the cell do not interact. To represent timing on the forward path, we use HSPICE-characterized three-dimensional tables for the delay and output transition times, indexed by the transition times of the two inputs, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{in}^{f}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{in}^{f}$$\end{document} , and the difference between the arrival times of the two input events, which ranges from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-W$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+W$$\end{document} . The precise value of the interaction window width, W, defined in section “Practical considerations for CMOS-based coupled oscillator systems”, is determined from HSPICE simulations. Each input may rise or fall, and the four resulting combinations of the transition types imply that we require four tables per coupling value. As a coupling cell implements \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2C_{max} + 1$$\end{document} levels, a total of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4(2C_{max} + 1)$$\end{document} three-dimensional tables are required. The timing arcs for the return paths ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{in}^{r}$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{out}^{r}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{in}^{r}$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{out}^{r}$$\end{document} ) do not interact as they are not coupled in the A2A architecture of section “A silicon-proven all-to-all coupled-RO Ising machine”. Therefore, the events on these arcs can be processed independently of each other, and one-dimensional tables will suffice as in the case of the enable cell.

The shorting cell has four inputs and four outputs that are labeled in the same way as the coupling cell, and the difference is that the coupling between the horizontal and vertical oscillators here is a short circuit. Since both the horizontal and vertical oscillators that meet at a shorting cell (i, i) are enabled by the same enable signal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$en_i$$\end{document} , any phase difference between them is a result of differences in coupling delays between the horizontal and vertical oscillators. Since Ising hardware uses weak coupling, these differential delays constitute a small fraction of the period. As a result, the arrival of a rising transition on the vertical RO will not be so severely delayed that it interacts with the falling transition of the horizontal RO at a shorting cell. Therefore, the lookup tables that capture rise-fall and fall-rise interactions are unnecessary, and two lookup tables suffice for shorting cells.

Overview of the event-driven simulator

The simulator requires the following inputs: (1) a timing file, with the characterized lookup tables and the interaction window (section “Capturing the timing information”); (2) the circuit netlist, a file that hierarchically captures the connections between devices and circuits; (3) a problem-specific coupling matrix, which maps the coupling coefficients of the Hamiltonian to the coupling cells, and is used to select the appropriate lookup table during simulation; (4) a maximum simulation time, which specifies the total simulation time; and (5) a tolerance value used to check for RO synchronization (section “CMOS-ring-oscillator-based Ising machines”). We use the following data structures:

• Event: an object that records an event, recording the net name, arrival time, transition time, and transition type (rising or falling).
• Q: a queue that sorts events by their arrival time, with the earliest occurring event at the head.
• Net2Event: a map with a net name as the key, pointing to an event at that net.
• PendingTrigger: a map with a net name as the key, pointing to a pending event with insufficient information for processing.

The simulator outputs are the map Net2Event and spin_vals, the set of spins that optimize the Hamiltonian for the coupling matrix.

An overview of the simulator is shown in Fig. 8, listing the event objects, the map Net2Event, and the scheduled events in queue Q. The simulator algorithm, described by pseudocode in the Supplementary Materials section, consists of the following steps:

Step 1: Initialize Initial events at the enable cells are scheduled to start the ROs. The queue, Q, and the map, Net2Event, are populated to reflect these events, and PendingTrigger is initialized to an empty map.

Step 2: Pop and process an event The earliest occurring event E is popped from Q. The event is passed to the PROCESS_EVENT function which generates new events that result from E. Consider an event that occurs at the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^f_{in}$$\end{document} pin of a coupling cell and the map Net2Event contains another event that occurs at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h^f_{in}$$\end{document} of the same cell. Then, PROCESS_EVENT will operate on these two events to generate events on output pins \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^f_{out}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h^f_{out}$$\end{document} of the coupling cell. We describe the PROCESS_EVENT function in section “PROCESS_EVENT: processing an event from the queue”.

Step 3: Check timeout and synchronization criteria The timeout criterion is met if the earliest event scheduled in Q has exceeded the maximum simulation time. The synchronization criterion, as defined in section “CMOS-ring-oscillator-based Ising machines”, is met when the periods of all coupled ROs are within the specified tolerance. We terminate the simulation when either of the above criteria is met.

Step 4: Assign spin values At the end of the simulation, the RO phases are translated to spin values, assigning a spin of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+1$$\end{document} to the reference RO. The phase difference between the RO in the A2A array and every other RO in the array is determined: if this phase difference is closer to 0 than it is to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} , a spin value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+1$$\end{document} is assigned to the RO, otherwise, we assign a spin value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1$$\end{document} .

PROCESS_EVENT: processing an event from the queue

We describe the intuition behind PROCESS_EVENT using an example to convey the complexities of looking forwards and backwards in time within the interaction window W; the pseudocode for PROCESS_EVENT is provided in our Supplementary Materials section.

Example 1: Figure 9 shows a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5 \times 5$$\end{document} A2A array, and focuses on three coupling cells within the array, as shown in the inset. The figure depicts two separate scenarios, CASE 1 and CASE 2, that will be used as examples in this subsection. We assume that the timing file specifies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W=$$\end{document} 75 ps for both examples. Consider the situation shown in CASE 1 of Fig. 9 where PROCESS_EVENT is called on a rising transition at net_a, which arrives at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=$$\end{document} 500 ps and has a transition time of 40 ps. The Net2Event map shows a rising transition on net_b at 520 ps, with a transition time of 35 ps. As the arrival time difference of the events is 20 ps which is less than the window, these events interact.

The output events are calculated using the three-dimensional lookup table mentioned in section “Capturing the timing information”. Note that if the event at net_b were to arrive at 580 ps instead of 520 ps, it would not interact with the event at net_a. In such a scenario, the event at net_a would be processed as a non-interacting event. The event(s) generated from processing net_a are inserted into Q, and any key-value pairs in Net2Event associated with net_a and any interacting event are removed. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square$$\end{document}

The above example considers events already in the Net2Event map, but the process could be complicated by as-yet-unprocessed events that could interact with a transition under consideration. For example, if net_b is not a key in Net2Event, LOOK_BACK is used to examine the predecessors of net_b to determine whether any upcoming event might interact with the event on net_a. We illustrate this with an example of a call to LOOK_BACK; the pseudocode for LOOK_BACK is in the Supplementary Materials section.

Example 2: Consider CASE 2 in Fig. 9 with events at net_a and net_c in Q. To process the event at net_a which arrives at 500 ps, an interacting event on net_b should arrive in the window (425 ps, 575 ps); there is no event in Net2Event with the key net_b. Thus, PROCESS_EVENT invokes LOOK_BACK with the arguments (net_b, (425 ps, 575 ps), 425 ps, Net2Event). The predecessor of net_b is net_c. Assume for this example, that the minimum and maximum delays of the coupling cell obtained from the timing file are 60 ps and 70 ps, respectively. An event that occurs on net_c can occur as early as 355 ps to incur the maximum delay of 70 ps and still generate an event on net_b in the required window. Similarly, an event on net_c can occur as late as 515 ps and incur the minimum delay of 60 ps to generate an interacting event on net_b. Thus, the window of arrival for an event on net_c is (355 ps, 515 ps).

Since Net2Event contains an event on net_c within this window, LOOK_BACK returns true. In PROCESS_EVENT, we stall the processing of net_a until net_b is scheduled, by adding the event to the map PendingTrigger with a key net_b. When the event at net_c is processed and it generates another at net_b, the pending event on net_a will be added back to the queue. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square$$\end{document}

Results and discussion

We have implemented DROID in Python 3.8.8 on a 64-bit RHEL 7.9 server with a 2.2GHz Intel^®^ Xeon^®^ Silver 4114 CPU. We used HSPICE and a commercial fast SPICE tool for circuit simulation using models from a commercial low-power 65nm technology. The characterization to build timing models for the cells, as described in section “Capturing the timing information”, took roughly 14 hours. This is a one-time step at a technology node, and its results can be reused across all designs and all simulations in the technology node. Therefore, like standard-cell library characterization, this amortization makes the runtime acceptable.

Comparison of DROID with GenAdler and HSPICE Figure 10 shows period waveforms of five ROs in a 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 5 A2A array as predicted by DROID, GenAdler, and HSPICE, for the same initial conditions, showing the evolution of RO states. The GenAdler simulation is based on available code^6^ and uses the tanh approximation for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{J_{ij}}$$\end{document} . In all plots, the period of each of the five ROs, labeled \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_0$$\end{document} through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {RO}_4$$\end{document} , evolves over the simulation until it reaches its final value. As pointed out in section “CMOS-ring-oscillator-based Ising machines”, the system tries to synchronize to a common period after which the phase differences between ROs become constant. It can be seen that the period waveform from DROID in Fig. 10a matches that from HSPICE (Fig. 10b), but the waveform from GenAdler (Fig. 10c) is noticeably different from HSPICE. Inaccuracies in the GenAdler waveforms are attributed to the limitations discussed in section “Limitations of the continuous-time approximation”.

Runtime trend as a function of array size We apply DROID to the silicon-proven A2A-coupled RO array^8^ (section “A silicon-proven all-to-all coupled-RO Ising machine”). Table 1 compares the runtimes of DROID, fast SPICE, and HSPICE for various array sizes, simulated up to 100 ns, for a single initial condition, over three array sizes, 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 5, 20 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 20, and 50 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 50, for random Ising problems. HSPICE is prohibitively slow and is about an order of magnitude slower than fast SPICE. DROID is two orders of magnitude faster than the fast SPICE simulator on a 50 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 50 array and is computationally efficient.

The analysis above considers a single graph topology and a single initial condition on the graph. We examine the scalability of DROID by determining the runtime for simulating the solution of Ising problems on A2A arrays of sizes up to 70 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 70, across various random graphs, and for multiple initial conditions. Recall from section “A silicon-proven all-to-all coupled-RO Ising machine”, that an A2A array of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \times N$$\end{document} is equivalent to a hexagonal/King’s graph array with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^2$$\end{document} coupled ROs^20,21^; therefore, a 70 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 70 A2A array has equivalent computation power to a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document} 4900-spin planar array, and avoids the problems of weakened spins caused by minor embedding in planar RO arrays. The random graphs are characterized by the value of graph density, defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2|E|/(|V||V|-1)$$\end{document} , where |V| is the number of vertices (i.e., ROs in the array) and |E| is edge count: the density varies from 0 for a graph with no edges to 1 for a complete graph.

DROID was exercised on arrays of various sizes: 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 5, 10 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 10, 20 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 20, 32 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 32, 50 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 50, and 70 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 70. For each size, we generate 80 random problems, 20 problems per graph density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\in \{0.4, 0.6, 0.8, 1.0\}$$\end{document} , with coupling values from a uniform probability distribution over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[-\,7, +\,7]-\{0\}$$\end{document} . Each problem was simulated with different initial conditions up to synchronization with a maximum simulation time of 250 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} s. The distribution of DROID runtimes that attained synchronization for each array size is shown in Fig. 11. Each box shows the lower quartile, median, and upper quartile of the runtimes over more than 1000 simulations for each array size. The median and upper quartiles of the synchronized runs for a 50 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 50 array are 10 minutes and 19.1 min, respectively, while for a 70 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 70 array these numbers are 33.1 min and 63.9 min, respectively. Such an experiment is not possible on the fast SPICE simulator, since even a single run takes over half an hour to simulate a 50 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 50 array for 100 ns, and simulating the array 1000 times for many microseconds would take more than 2 months.

Figure 12 shows the evolution of phases for a representative problem in a 50-spin A2A system, over multiple snapshots of time, starting from random phases. The blue and red colors correspond to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm 1$$\end{document} spins. As the state of the array evolves, the spins go through the continuum between the two states, shown by the color bar, until they converge to the final state, seen after 700 cycles.

Comparison with hardware Since the A2A Ising hardware platform does not allow us to probe terminals to observe the states of the ROs, we compare the synchronized state from the hardware with that from a software simulation. The A2A hardware has a 50 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 50 array^8^, using one RO from the array as a sampling clock and another as the reference spin. An arbitrary 48-spin Hamiltonian needs \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(48+1)$$\end{document} A2A ROs after accounting for the reference spin to implement the linear terms. Therefore, as a test set, we generate a set of 48-spin problems for five random values of graph density, {0.2, 0.4, 0.6, 0.8, 1.0}, with 50 problems for each density. As before, the coupling values for edges are chosen uniformly in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[-7, +7]-\{0\}$$\end{document} . For each problem, 100 samples obtained from the hardware are compared with DROID under randomized initial conditions. We compare the distributions of Hamiltonian (1), for the sample of random coupling weights, for the solutions from the hardware solver and the simulator.

The distributions of the optimized Hamiltonian from DROID and the hardware solver, for a specific problem instance with a density of 0.2, are shown in Fig. 13. The x-axis is normalized to the optimum Hamiltonian for the problem, and the bin size is 5% of the optimum Hamiltonian value, i.e., 0.05 units on this normalized axis. Visually, the histograms are seen to be very similar.

Quantitatively, we use Earth Mover Distance (EMD)^23^ to measure the similarity of the two distributions. The EMD quantifies the amount of work required to transform one distribution to the other and is robust to small changes to the bin size. Intuitively, if we view one distribution as a mass of earth spread in space and the other as a collection of holes in the very same space, the EMD measures the least amount of work needed to fill the holes with earth; if the two distributions are identical, the work is zero^23^. A large EMD value suggests dissimilarity of distributions as more work is required to transform one distribution to the other, and therefore small values of EMD are desirable. The EMD between the distributions in Fig. 13 is 0.009: this can be thought of as moving 9 of 100 samples over a distance of 10% of the x-axis, which indicates that the distributions are close.

Table 2 summarizes the EMD between the distributions across problems from the test set described earlier. The mean, minimum, maximum, and standard deviation of the EMD values for each graph density value correspond to a row in the table. The small EMD values indicate that over the test set, the distribution of solutions from DROID resembles that from the hardware.

Conclusion

This paper has presented DROID, an event-driven simulator for RO-based Ising machines. It is analytically shown that its formulation is more general than conventional continuous-time models, and that its event-driven nature captures effects that have not been addressed in previous analytical solvers. DROID is shown to be nearly two orders of magnitude faster than a commercial fast SPICE simulator, at similar accuracy on a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$70 \times 70$$\end{document} array, which has the computational power of a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document} 4900-spin planar array. Since A2A arrays are size-limited to about this dimension (larger problems can be solved using A2A arrays using decomposition^22^) due to noise limitations, DROID has demonstrated its utility for simulating arrays of practical size. This speedup enables us to simulate the all-to-all RO array with multiple initial conditions and obtain distributions of solutions that simulate the hardware.

Fig. 1. Three five-stage ROs, with a positive and a negative coupling. The green stage is the reference, and odd (even) stages are shown in red (blue).

Fig. 2(a) Waveforms for nets a and b in two ROs, when uncoupled, and with coupling enabled at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document} . (b) Detailed view of the green box in (a), showing the effect of coupling. The period of waveform at a increases while that at b decreases, reducing the phase difference.

Fig. 3. The relative arrival times on nets a (green waveform) and b (blue waveform), top, impact the transition delay on net c (middle). The transition delay as a function of phase difference, and the interaction window W are illustrated at the bottom.

Fig. 4. An illustration of the A2A concept through a small three-RO structure, showing shorting cells S on the diagonal, programmable off-diagonal coupling cells C, and enable cells EN.

Fig. 5(a) The results for an example characterization setup, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{ij}=1$$\end{document} , showing the delay shift \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{J_{ij}}$$\end{document} as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij}$$\end{document} . (b) HSPICE characterization of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{J_{ij}}$$\end{document} and its tanh approximation as a function of phase difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij}$$\end{document} .

Fig. 6(a) A section of the A2A array with three coupling cells. (b) Phase difference at coupling cells depends on the delay difference at the inputs of the coupling cell and varies with location. (c) The delay introduced by previous coupling cells affects the phase differences at later stages within the same cycle.

Fig. 7. The phase differences across various edges are shown above. The phase differences \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{AB}$$\end{document} (red) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{BC}$$\end{document} (blue) stay close to 0 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} between 20 and 400 ns giving the impression of being phase-locked. The phase differences deviate from these seemingly phase-locked positions beyond 400 ns.

Fig. 8A simulator step with the sorted queue Q and the Net2Event map: PROCESS_EVENT consumes one or more events and generates future events.

Fig. 9. The handling of an event on net_a is influenced by the knowledge of events on nets in some neighborhood around it, as shown on the top. CASE 1 shows the scenario when an event on the other input (net_b) of the same instance is known. CASE 2 shows the scenario where LOOK_BACK is invoked to find an event on net_c that can cause an event on net_b that might lie within the interaction window of the event at net_a.

Table 1. Simulation runtimes for various array sizes.Array size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5 \times 5$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20 \times 20$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$50 \times 50$$\end{document}

MOSFETs310052,000328,000DROID2.1 s3.5 s7.9 sFast SPICE Level 6 (1 core)34 s (16 $\documentclass[12pt]{minimal}

			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\times$$\end{document}$ )283 s (81 $\documentclass[12pt]{minimal}
			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\times$$\end{document}$ )0.68 h (310 $\documentclass[12pt]{minimal}
			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\times$$\end{document}$ )Fast SPICE Level 6 (8 cores)24 s (11 $\documentclass[12pt]{minimal}
			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\times$$\end{document}$ )111 s (32 $\documentclass[12pt]{minimal}
			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\times$$\end{document}$ )0.24 h (109 $\documentclass[12pt]{minimal}
			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\times$$\end{document}$ )HSPICE Level 5262 s (125 $\documentclass[12pt]{minimal}
			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\times$$\end{document}$ )1.04 h (1072 $\documentclass[12pt]{minimal}
			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\times$$\end{document}$ )16.33 h (7441 $\documentclass[12pt]{minimal}
			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\times$$\end{document}$ )

Fig. 10. Period waveforms for a 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 5 A2A array from (a) DROID, (b) HSPICE, and (c) GenAdler for the same initial conditions of the ROs.

Fig. 11. Distribution of DROID runtimes for various A2A array sizes.

Fig. 12. The evolution of RO-phases from random to binarized phases for a particular problem; phases are represented by the color bar at right.

Fig. 13. Distributions of the Hamiltonian from 100 samples of a problem of density 0.2 obtained from (a) DROID and (b) the hardware. The X-axis has been normalized to the best Hamiltonian value obtained for the problem.

Table 2EMD between distributions from DROID and A2A Ising chip for 5 graph densities, 50 problems/density, 100 samples/problem.DensityEMDMeanMinimumMaximumS.D.0.20.0210.0050.0730.0160.40.0210.0040.0540.0120.60.0140.0020.0360.0090.80.0170.0030.0400.0081.00.0160.0010.0440.009Average0.0180.011

Supplementary Information

Supplementary Information.