Oscillator Ising Machines (OIMs)

• Oscillator Ising Machines are systems where binary variables are encoded in the phases (0 and π) of nonlinear oscillators.
• They employ programmable pairwise coupling and sub-harmonic injection locking to map the oscillator network to the Ising Hamiltonian for optimization tasks.
• OIMs have been demonstrated on various hardware platforms, achieving rapid convergence and improved performance on benchmarks such as MAX-CUT.

Oscillator Ising Machines (OIMs) are Ising machines in which binary variables are represented by the phases of coupled self-sustaining nonlinear oscillators rather than by static bistable devices. Their central construction is to combine programmable pairwise coupling with a binarizing perturbation—typically second-harmonic or sub-harmonic injection locking—so that the network evolves on a continuous phase manifold but settles into two phase wells, usually $0$ and $\pi$, that encode spins $s_i \in \{+1,-1\}$. In the canonical formulation, the network admits an energy or Lyapunov function that coincides with the Ising Hamiltonian on the binary-phase manifold, so low-energy steady states of the oscillator system implement low-energy, and often ground-state, solutions of combinatorial optimization problems such as MAX-CUT (Wang et al., 2017, Wang et al., 2019).

1. Formal definition and Ising correspondence

The optimization target is the Ising Hamiltonian. Across the OIM literature, sign conventions vary, but the standard zero-field form is written as

$$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$

with symmetric couplings $J_{ij}$ and optional local fields $h_i$ (Vosoughi, 2020). Equivalent formulations absorb the linear term into an augmented spin vector by introducing a fixed spin $s_{n+1}=+1$, or write the quadratic term with the opposite sign and reinterpret $J_{ij}$ accordingly (Wang et al., 2017).

In OIMs, each spin is encoded by oscillator phase. The binary mapping most commonly used is

$$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$

with $\phi_i \in \{0,\pi\}$ corresponding to $\pi$0 and $\pi$1, respectively (Wang et al., 2017, Wang et al., 2019). In reference-based formulations, the phase is measured relative to a designated oscillator,

$\pi$2

so that in-phase and anti-phase states relative to the reference encode the two spin values (Roy et al., 5 Feb 2025).

The mapping from combinatorial problems to couplings is direct for MAX-CUT. For a weighted graph with edge weights $\pi$3, one common convention sets $\pi$4, so that opposite spins lower the Ising energy and correspond to cut edges (Roy et al., 5 Feb 2025, Wang et al., 2019). Under an alternative convention used in some circuit papers, $\pi$5 on edges and the sign is absorbed into the oscillator coupling term (Vosoughi, 2020). This suggests that, in OIM practice, problem embedding is less about a unique algebraic convention than about maintaining consistency between the physical coupling sign and the discrete cost function.

A key convenience for hardware mappings is that local fields can be realized either by coupling to a fixed reference spin or by adding a first-harmonic bias term proportional to $\pi$6 (Wang et al., 2017). That observation underlies later extensions in which OIMs are used not only for zero-field MAX-CUT instances but also for bias-driven energy-based learning and probabilistic sampling (Gower, 14 Oct 2025, Ekanayake et al., 21 Aug 2025).

2. Phase dynamics, SHIL/SHI, and energy descent

The canonical OIM dynamics are Kuramoto/Adler-type phase equations augmented by a binarization term. A representative form is

$\pi$7

where $\pi$8 models frequency mismatch, $\pi$9 sets mutual coupling strength, $s_i \in \{+1,-1\}$0 is the sub-harmonic or second-harmonic injection strength, and $s_i \in \{+1,-1\}$1 represents Gaussian white noise (Wang et al., 2017). Closely related formulations write

$s_i \in \{+1,-1\}$2

or add bias terms $s_i \in \{+1,-1\}$3 for unary costs (Wang et al., 2019, Gower, 14 Oct 2025).

The term $s_i \in \{+1,-1\}$4 is the defining OIM regularizer. Under injection near twice the oscillator frequency, it creates two stable phase equilibria separated by $s_i \in \{+1,-1\}$5, so the continuous phase state is driven onto a bistable manifold (Wang et al., 2017, Vosoughi, 2020). Terminology varies across papers: some describe this as sub-harmonic injection locking (SHIL), others as second-harmonic injection (SHI), but the dynamical role is the same—phase binarization into $s_i \in \{+1,-1\}$6 (Farasat et al., 26 Aug 2025).

For the sinusoidal coupling model with symmetric weights, the network is gradient-like and admits a Lyapunov function. One common expression is

$s_i \in \{+1,-1\}$7

with

$s_i \in \{+1,-1\}$8

in the deterministic phase-only model (Wang et al., 2019). Restricting to binary phases $s_i \in \{+1,-1\}$9 gives $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$0 and $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$1, so $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$2 reduces to the Ising quadratic form up to a constant offset (Wang et al., 2017).

Operationally, reported OIM schedules initialize phases randomly, apply a global SYNC signal near $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$3, ramp mutual coupling $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$4 or $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$5, and maintain controlled noise to permit barrier crossing before final readout via $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$6 (Wang et al., 2017, Wang et al., 2019). Later work generalizes the interaction beyond pure sinusoids by using the Generalized Adler Equation or by shaping oscillator perturbation projection vectors (PPVs) and waveforms; smooth square- and triangle-like interactions were reported to improve performance on large MAX-CUT benchmarks (Wang et al., 2017).

Although phase-only models dominate the early OIM literature, more recent phase–amplitude formulations replace scalar phase variables by complex oscillator states $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$7 and derive a real-valued energy $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$8 through conjugate-symmetrized polynomial interactions. That construction extends OIMs to arbitrary-order cost functions and preserves monotonic energy decrease under ideal Wirtinger-gradient dynamics (Sun et al., 1 Apr 2025).

3. Stability structure, bias, and theoretical controversies

The existence of an Ising-like Lyapunov function does not by itself determine which equilibria are dynamically selected. A control-theoretic analysis showed that OIM functionality depends on the relative magnitude of the second-harmonic forcing $$H(\mathbf{s}) = - \sum_{i<j} J_{ij}\, s_i s_j - \sum_i h_i\, s_i, \qquad s_i \in \{\pm 1\},$$9 and the coupling $J_{ij}$0: if $J_{ij}$1 is too small, phases need not binarize; if it is too large, many suboptimal local minima become stable (Bashar et al., 2023). In a 20-node, 152-edge example with $J_{ij}$2, values $J_{ij}$3 destabilized even globally optimal phase configurations, $J_{ij}$4 stabilized only a subset of global minima, and $J_{ij}$5 stabilized all global minima but also many local minima, increasing trapping (Bashar et al., 2023).

A complementary note formalized the relation between linearization and energy curvature: for the standard phase-only OIM,

$J_{ij}$6

so positive definiteness of the Hessian of the energy at a fixed point is equivalent to negative definiteness of the Jacobian of the dynamics at that point (Bashar et al., 2023). This result does not remove the bias problem; it clarifies that local asymptotic stability and strict local minimality are the same statement for gradient-form OIMs.

Several later papers sharpened the point that not all binary equilibria are computationally equal. One study showed analytically that the Jacobian spectrum of an OIM depends explicitly on the spin configuration through the signed interaction structure, so higher-energy states can be selectively destabilized, whereas in bistable latch Ising machines all discrete Ising configurations have identical linear stability (Hasan et al., 6 Mar 2026). Another study recast stability in terms of a signed graph Laplacian and showed that heterogeneous regularization parameters $J_{ij}$7 can bias asymptotic stability toward lower-energy states; in the frustration-free case, global minimizers are stable for any $J_{ij}$8, while sufficiently small $J_{ij}$9 destabilize nonoptimal configurations (Allibhoy et al., 6 May 2025).

The binarizing SHI term is itself double-edged. A 2025 analysis introduced “spin freezing” to describe the regime in which the SHI restoring torque prevents an oscillator phase from crossing the neutral thresholds $h_i$0 or $h_i$1, even when such a spin flip would reduce the Ising energy (Farasat et al., 26 Aug 2025). For the stochastic model

$h_i$2

the paper derived a local condition for an unfrozen spin,

$h_i$3

and reported that initializing all oscillators at $h_i$4 or $h_i$5 delayed freezing and consistently yielded higher-quality cuts than random initialization on Erdős–Rényi MAX-CUT instances (Farasat et al., 26 Aug 2025).

An even more fundamental qualification concerns the gradient-flow assumption itself. A 2026 analysis showed that gradient-flow dynamics are not generic across hardware-realistic oscillators; they require a harmonic-by-harmonic quadrature relation between the oscillator waveform and its infinitesimal phase response curve (iPRC) (Hasan et al., 17 May 2026). Harmonic misalignment introduces even components into the pairwise interaction function, producing non-conservative phase dynamics with nonzero curl. The paper quantified this departure through

$h_i$6

and reported substantial non-gradient contributions in representative ring-oscillator models (Hasan et al., 17 May 2026). This suggests that the classical OIM Lyapunov picture is exact for some oscillator classes and parameterizations, but not universal across all implementations.

4. Hardware realizations and architectural variants

The original OIM proposal emphasized that “almost any nonlinear self-sustaining oscillators” can, in principle, implement the paradigm, including CMOS LC oscillators, CMOS ring oscillators, MEMS resonators such as Resonant Body Transistors, and spin-torque nano-oscillators (Wang et al., 2017). Subsequent work broadened this into a hardware taxonomy spanning breadboard demonstrators, PCB systems, analog CMOS macros, injection-distribution networks, on-chip SHIL generators, GPU simulations, and fixed-point ASIC emulations.

Implementation class	Representative reported details	Reported outcome
CMOS LC and early prototypes	six coupled $h_i$7 LC oscillators in ngspice-26; four-oscillator breadboard at $h_i$8; 8/32/64/240-spin prototypes with programmable couplings	correct small MAX-CUT solutions; OIM240 delivered solutions every $h_i$9 at $s_{n+1}=+1$0
Injection-distribution optimized CMOS OIM	$s_{n+1}=+1$1 coupled LC oscillators at $s_{n+1}=+1$2 in $s_{n+1}=+1$3 PTM bulk CMOS with centralized vs distributed current mirrors	phase-lock time reduced from $s_{n+1}=+1$4 to $s_{n+1}=+1$5; $s_{n+1}=+1$6 speed improvement
Off-the-shelf analog experimental OIM	phase-shift oscillators at approximately $s_{n+1}=+1$7, $s_{n+1}=+1$8 matrix of 10-bit digital potentiometers, SHIL around $s_{n+1}=+1$9 and $J_{ij}$0	convergence in a few oscillation periods on 4-node and 8-node MAX-CUT instances
On-chip SHIL and digital emulation	ROA-SHIL at $J_{ij}$1 for a 324-node ROSC OIM; $J_{ij}$2 65 nm fixed-point ASIC at $J_{ij}$3	ROA-SHIL preserved $J_{ij}$4 to $J_{ij}$5 accuracy under PVT variation; ASIC solved 400-node tasks in $J_{ij}$6

The earliest reported circuit demonstrations used coupled LC oscillators with resistive interconnects whose conductances encoded $J_{ij}$7. The 2017 study reported six coupled oscillators at $J_{ij}$8 converging in $J_{ij}$9, and a four-oscillator breadboard built from ALD1106/7 CMOS devices, $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$0 inductors, $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$1 capacitors, and six $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$2 potentiometers (Wang et al., 2017). A 2019 late-breaking paper expanded this hardware line to 8-, 32-, 64-, and 240-oscillator systems, with the 240-spin prototype using 200 AD5206 digital potentiometers for 1200 programmable couplings and producing solutions every $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$3 after a $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$4 synchronization window (Wang et al., 2019).

Injection delivery emerged as a distinct architectural issue. In a 10-oscillator LC OIM simulated in Virtuoso Cadence and ADS Keysight with a $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$5 PTM model, centralized SHIL routing introduced unintended coupling terms through the shared injection network, whereas distributed current mirrors restored the intended Kuramoto-plus-injection dynamics (Vosoughi, 2020). At total injection current $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$6, the reported phase-lock time dropped from $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$7 to $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$8, the average phase-lock error decreased by $$s_i = \cos \phi_i \quad\text{or}\quad s_i = \operatorname{sign}(\cos \phi_i),$$9, and total power increased by only $\phi_i \in \{0,\pi\}$0 to $\phi_i \in \{0,\pi\}$1 (Vosoughi, 2020).

A separate experimental line used phase-shift oscillators rather than LC tanks. The 2025 “Experiments with an oscillator based Ising machine” paper built 4- and 8-oscillator systems from off-the-shelf op-amp phase-shift oscillators, digital potentiometer matrices, analog summers, multipliers, and integrators, with approximate resonance at $\phi_i \in \{0,\pi\}$2 and SHIL around $\phi_i \in \{0,\pi\}$3 (Roy et al., 5 Feb 2025). That work stressed accessibility and rapid convergence—typically a few oscillation periods—while also documenting practical issues such as cross-coupling via shared power rails, threshold tuning in multiplier–integrator phase detectors, and the scaling burden of all-to-all coupling, which grows as $\phi_i \in \{0,\pi\}$4 (Roy et al., 5 Feb 2025).

On-chip SHIL generation and distribution became a separate scaling problem. A 2026 study proposed rotary oscillator array (ROA) bricks based on rotary traveling-wave oscillators to generate a coherent $\phi_i \in \{0,\pi\}$5 SHIL signal for a ring-oscillator OIM with $\phi_i \in \{0,\pi\}$6 (Sica et al., 18 May 2026). In a 324-node MAX-CUT benchmark, ROA-SHIL preserved binarization and $\phi_i \in \{0,\pi\}$7 to $\phi_i \in \{0,\pi\}$8 accuracy under voltage, temperature, and process variation, whereas a distributed ring-oscillator SHIL source failed to maintain injection locking under several PVT corners (Sica et al., 18 May 2026).

Parallel to analog hardware, digital emulations increasingly adopted OIM dynamics rather than conventional annealing. One 2026 ASIC emulated OIM/OPM phase updates in 8-bit fixed point on a $\phi_i \in \{0,\pi\}$9 processing-element array with king’s-graph connectivity, achieving $\pi$00 post-layout frequency, about $\pi$01 active power, and $\pi$02 solution time for 1000 Euler steps (Gonul et al., 15 Apr 2026).

5. Benchmarks, applications, and functional extensions

MAX-CUT has been the dominant OIM benchmark from the outset. The 2017 paper reported reliable convergence to the global optimum on a fully connected 6-node random MAX-CUT instance in both phase macromodel and LC-oscillator SPICE simulations, then evaluated the large G22 instance with $\pi$03 and $\pi$04 edges (Wang et al., 2017). Over 100 runs with random initialization, ramped coupling, and a nontrivial $\pi$05 schedule, the reported mean cut was 13253 and the best cut was 13305, compared with coherent Ising machine results of mean 13248 and best 13313; removing noise lowered the mean to 13203, and removing SYNC lowered it further to 13050 (Wang et al., 2017).

A broader 2019 arXiv study ran OIMs on all 54 G-set MAX-CUT instances, sizes 800 to 3000, and reported matching best-known cut values for 38 of 54 instances while improving 17 beyond prior published results (Wang et al., 2019). A contemporaneous late-breaking results paper, using 100 random initializations and fixed 1000-cycle simulations, reported matching 21 of 54 and improving 17 of 54 previously published optimal solutions, while also presenting hardware energy histograms for a 240-spin programmable OIM on random Ising instances (Wang et al., 2019).

Graph coloring was already present in the early literature through a 4-coloring formulation of the United States map using four spins per vertex (Wang et al., 2019). Later work generalized the oscillator formalism from Ising to Potts variables. The 2026 ASIC-emulated OIM/OPM reported 97–100% maximum accuracy on 400-node unweighted and weighted MAX-CUT problems and a mean accuracy of 92.02% with 4.22% standard deviation on a 400-node graph 3-coloring problem (Gonul et al., 15 Apr 2026). A GPU-based simulator extended the same line to GSET and SATLIB, reporting up to 11295x speed-up over CPUs with up to 99.27% accuracy on large MAX-CUT instances and at least 95–99.5% constraint satisfaction on reported 3-coloring benchmarks (Gonul et al., 28 May 2025).

Higher-order optimization objectives moved OIMs beyond pairwise Ising formulations. A 2025 phase–amplitude OIM introduced a real-valued complex energy for polynomial interactions and evaluated 3-SAT instances from the SATLIB uf-20, uf-50, uf-75, uf-100, and uf-150 sets, totaling 1700 instances (Sun et al., 1 Apr 2025). The proposed model solved a higher percentage of instances than the prior complex-valued Hopf baseline and achieved lower final energies more frequently, while preserving monotonic energy descent under ideal Wirtinger-gradient dynamics (Sun et al., 1 Apr 2025).

Functional extensions now reach beyond deterministic combinatorial optimization. One 2025 paper showed that OIMs can be configured as p-bit engines by combining first-harmonic injection with SHI, thereby implementing a stochastic update law of the form

$\pi$06

and demonstrated a 15-node, 59-edge MAX-CUT instance reaching the global optimum cut of 39 after stochastic sampling events improved on an initial deterministic OIM solution (Ekanayake et al., 21 Aug 2025).

Another extension uses OIMs as trainable energy-based models. Equilibrium Propagation was mapped directly onto OIM hardware by identifying model couplings, biases, and synchronization fields with the OIM parameters $\pi$07, $\pi$08, and $\pi$09 (Gower, 4 May 2025, Gower, 14 Oct 2025). The reported results include $\pi$10 test accuracy on MNIST without hardware modification, robustness to moderate phase noise and 10-bit parameter quantization, and in a later study $\pi$11 on MNIST and $\pi$12 on Fashion-MNIST (Gower, 4 May 2025, Gower, 14 Oct 2025).

A different learning-oriented direction trains OIM couplings to assign dynamic stability rather than merely low energy. The Hamiltonian-Regularized Eigenvalue Contrastive Method (HRECM) exploits the fact that all $\pi$13 equilibria are structurally stable while their dynamic stability depends on the coupling matrix. In this formulation, desired binary patterns are made asymptotically stable and undesired ones unstable, yielding Hopfield-like associative memory (Cheng et al., 18 Jul 2025).

6. Limitations, open questions, and research directions

Several limitations recur across the literature. The first is scaling of programmability. Analog all-to-all coupling grows as $\pi$14, which is explicitly visible in small experimental systems: 12 programmable couplings for 4 oscillators and 56 for 8 oscillators in the phase-shift prototype, with the extrapolation to $\pi$15 oscillators implying $\pi$16 weights (Roy et al., 5 Feb 2025). Dense programmability also competes with routing, cross-talk suppression, delay management, and SHIL distribution, which motivated both distributed current-mirror injection in 2020 and ROA-brick SHIL in 2026 (Vosoughi, 2020, Sica et al., 18 May 2026).

A second limitation is model fidelity. The standard phase-only OIM assumes symmetric couplings, negligible amplitude dynamics, and interaction functions that are effectively odd in phase difference. That assumption is challenged both by the amplitude-aware polynomial model, which explicitly adds amplitude dynamics and higher-order interactions, and by the harmonic-misalignment analysis, which argues that hardware-realistic oscillators can violate gradient-flow conditions altogether (Sun et al., 1 Apr 2025, Hasan et al., 17 May 2026). This suggests that future OIM theory may need nonequilibrium tools in addition to Lyapunov analysis.

A third limitation is the trade-off between exploration and binarization. SHI is necessary for clean phase readout and Ising mapping, but it can also induce premature locking, spin freezing, and enlarged basins for suboptimal states (Farasat et al., 26 Aug 2025, Bashar et al., 2023). Reported mitigations include ramping $\pi$17, injecting controlled noise, introducing heterogeneous $\pi$18, using deterministic initializations such as $\pi$19 or $\pi$20, and periodically turning synchronization on and off (Wang et al., 2019, Allibhoy et al., 6 May 2025, Farasat et al., 26 Aug 2025).

A fourth issue is graph dependence. One 2025 study introduced the Dynamical Ising Machine (DIM), which minimizes the same Ising Hamiltonian as an OIM but uses phase-sum rather than phase-difference coupling and exhibits a markedly different bifurcation and stability structure (Ekanayake et al., 21 Mar 2025). The reported results on random graphs showed that OIM and DIM success frequencies depend on graph topology, and that running multiple dynamical systems with different stability structures can reduce sensitivity to the input graph (Ekanayake et al., 21 Mar 2025). A plausible implication is that “physics-inspired optimization” may ultimately be an ensemble design problem rather than a single-model problem.

Open directions in the OIM literature are correspondingly broad: device-specific PPV and waveform engineering; integrated implementations across CMOS rings, LC tanks, MEMS, and spin-torque devices; higher-order and multi-valued phase states; adaptive or node-wise SHIL; delayed coupling; broader benchmark comparisons; associative-memory training; and learning-centric co-design in which OIMs act simultaneously as optimizers, samplers, and trainable energy-based processors (Wang et al., 2017, Sun et al., 1 Apr 2025, Cheng et al., 18 Jul 2025, Gower, 14 Oct 2025). What remains stable across these variants is the defining OIM idea: encode discrete variables in oscillator phases, shape the phase landscape so that binary or multi-valued equilibria correspond to computational states, and exploit the collective synchronization dynamics of the network as the computation itself.