Mean-Field Oscillator Ising Machines: Gradient Flows and Classification of Limit Solutions

Abstract: Oscillator Ising Machines (OIMs) have emerged as promising computational architectures for approximating solutions to combinatorial optimization problems. We derive and analyze the mean-field limit of an OIM model and show that it inherits the gradient-flow structure of the finite-dimensional dynamics. We identify conditions under which this mean-field evolution admits an Eulerian formulation as a gradient flow on the Wasserstein space of probability measures, and contrast this with a Lagrangian formulation which is always available. The gradient-flow structure strongly constrains the long-time dynamics and enables a complete classification of limit solutions and their stability in the symmetric case. In particular, all limit solutions are fixed points whose phases cluster into at most four groups, and for almost all parameter values, only binarized fixed points -- those with clusters at $0$ and/or $π$ -- can be stable. Since binarized states are exactly those for which a feasible solution to the original problem can be read out, this shows that feasible solutions can almost always be recovered. We provide tight bounds on the parameter thresholds for which fixed points in this binarized family are stable, thereby identifying the threshold for binarization in this model. We also present numerical evidence that the mean-field model correctly predicts behavioral regimes in large random networks, including Erdős-Rényi networks.

• The paper derives Lagrangian $L^2$ and Eulerian 2-Wasserstein gradient-flow models, proving global well-posedness and energy dissipation for the mean-field dynamics.
• The paper classifies every equilibrium in the totally symmetric case as an atomic phase distribution with at most four clusters and provides an algorithm for enumerating them.
• The paper proves that, except at $K/K_s=-2$, only binarized states can be directionally stable, with nontrivial two-cluster stability requiring $-2<K/K_s<2/|2m-1|$, a result supported by simulations on four random-network classes.

Overview and motivation

Oscillator Ising Machines (OIMs) solve combinatorial optimization problems by embedding the discrete Ising Hamiltonian into a continuous surrogate over coupled Kuramoto oscillators with second-harmonic injection locking, then running gradient descent dynamics that can be implemented directly in analog hardware such as FPGAs or optical lasers. The paper under review develops a rigorous mean-field theory for this architecture. Its central contributions are: (i) a formal derivation of Lagrangian (phase-function) and Eulerian (phase-density) mean-field models, both of which inherit the gradient-flow structure of the finite-dimensional system; (ii) a complete classification of limit solutions in the totally symmetric case; and (iii) tight stability bounds showing that only binarized fixed points—those whose phases cluster at $0$ and/or $\pi$—can be stable for almost all parameter values.

The practical significance of the last point is immediate: binarized states are exactly those from which a feasible spin configuration $\sigma_i = \cos(\theta_i)$ can be read out for the original problem. The analysis therefore establishes that feasible solutions are recoverable from the mean-field dynamics almost always, and identifies the parameter threshold governing when nontrivial binarized states exist.

Mean-field models

Starting from finite-$N$ energy $$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$ with dense-graph limits described by graphons $K(x,y)$, the Lagrangian model replaces indices with continuum variables:

$$\partial_t \theta(x) = \int K(x,y)\sin(\theta(x)-\theta(y))\,\mathrm{d}y - K_s\sin(2\theta(x)).$$

In the totally symmetric case $K(x,y)\equiv K$, relabeling invariance makes the index set extraneous, and the state pushes forward to a phase density $\rho = \theta_\#\mathscr{L}$ on the torus. The resulting Eulerian model is a continuity equation

$$\partial_t\rho = -\nabla\cdot\big(\rho\,[K(\sin * \rho) - K_s\sin(2\,\cdot)]\big),$$

with energy expressible through the first order parameter $\pi$0. The pushforward map is many-to-one (phase functions related by measure-preserving relabelings yield the same density), so the Eulerian description is well-defined precisely when the system is relabeling-invariant.

Gradient-flow structure

Both models are gradient flows, but with respect to different geometries. The Lagrangian dynamics are the gradient flow of the energy with respect to the $\pi$1 Riemannian metric on the Hilbert manifold of phase functions; since $\pi$2 bounded implies the energy is $\pi$3 with globally Lipschitz gradient, the flow is globally well-posed and satisfies an Energy Dissipation Equality $\pi$4.

The Eulerian model is a gradient flow with respect to the 2-Wasserstein metric, which the authors identify as exactly the quotient of the $\pi$5 metric modulo measure-preserving relabelings—a structural explanation for Wasserstein's ubiquity in collective systems where labels are extraneous. However, the Wasserstein space has a boundary (atomic measures) where the Riemannian structure degenerates, and all fixed points of interest turn out to be atomic. Neither working on the interior nor using metric-space gradient-flow theory provides a usable second-order stability calculus at these boundary points, which is why the stability analysis is carried out in the Lagrangian setting.

A subtle regularity issue shapes the entire treatment: the infinite-dimensional energy is merely $\pi$6, so positivity of the Hessian does not imply that a fixed point is a true local minimum. Instead it yields a directional local minimum, necessary but not sufficient for asymptotic stability in the limit. Notably, the binarized fixed points form a continuous one-parameter family of strong directional local minima that are not local minima—an instructive counterexample to finite-dimensional intuition. Because first- and second-order information transfers between finite and infinite dimensions while higher-order information does not, Hessian conditions still correctly predict observed behavior in large finite systems.

In the constant-$\pi$7 case, monotone rearrangement plus the comparison principle bounds total variation by $\pi$8, giving precompactness via Helly's selection principle. LaSalle–Krasovskii then shows every trajectory's accumulation set is a connected union of fixed points at a single energy level. Full convergence of trajectories to individual fixed points remains open: neither isolation of fixed points nor a Łojasiewicz inequality holds in the limit (fixed-point sets are generically two-dimensional). A conditional route exists through the order parameter—the authors prove $\pi$9 converges if and only if $\sigma_i = \cos(\theta_i)$0 converges—but the available estimate gives only $\sigma_i = \cos(\theta_i)$1, whereas $\sigma_i = \cos(\theta_i)$2 integrability would be needed.

Classification of equilibria

Fixed points require $\sigma_i = \cos(\theta_i)$3 supported on the zero set of the velocity field $\sigma_i = \cos(\theta_i)$4, together with self-consistency of $\sigma_i = \cos(\theta_i)$5. Rewriting $\sigma_i = \cos(\theta_i)$6 as $\sigma_i = \cos(\theta_i)$7 with $\sigma_i = \cos(\theta_i)$8, its zeros are the intersections of a rectangular hyperbola and the unit circle; Bézout's theorem yields at most four zeros, and the hyperbola's passage through the origin guarantees at least two. Every equilibrium is therefore atomic with at most four Dirac masses, with masses given by convex coefficients expressing $\sigma_i = \cos(\theta_i)$9 in the zero set. This yields a concrete five-step algorithm enumerating all equilibria from the two-parameter sweep $N$0.

Stability and the binarization threshold

The stability analysis proceeds by testing the second-variation quadratic form against progressively richer perturbation families. Two results complete the classification:

• Necessary condition: for almost all values of $N$1, any directional local minimum must be valued in $N$2. The single exceptional case occurs at $N$3 with two equal masses at symmetric phases $N$4, where the Jacobian has a zero eigenvalue and the equilibria appear neutrally stable.
• Sufficient condition: for two-cluster binarized states with mass fraction $N$5 at phase $N$6, the following bounds hold, and they are tight (necessary when inequalities are made non-strict):

Mass fraction	Stability condition
$N$7	$N$8
$N$9	$$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$0
otherwise	$$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$1

These bounds identify the threshold for binarization at $$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$2: below it, only single clusters at $$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$3 or $$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$4 are stable, so no nontrivial solution structure survives. The proofs rest on sharp Cauchy–Schwarz-type bounds on the quadratic form; notably, the naive bound for negative $$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$5 in the general-$$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$6 case is not tight, and a more careful variational argument (minimizing over mass-splitting perturbations) recovers the exact constant.

Numerical validation

Simulations on $$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$7 random networks—constant-weight and uniformly weighted all-to-all graphs, binary and weighted Erdős–Rényi graphs—perturb binarized states and record the settled mass fraction against $$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$8. The observed stability region agrees closely with the theoretical wedge $$E = \tfrac{1}{2}\sum_{ij} A_{ij}\cos(\theta_i - \theta_j) - \tfrac{K_s}{2}\sum_i \cos(2\theta_i)$$9 across all four network classes, supporting the claim that the mean-field model predicts behavioral regimes of large homogeneous random graphs. This agreement is expected to be good precisely because only first- and second-order information transfers across the finite-to-infinite-dimensional limit.

Limitations and open questions

Several restrictions bound the scope of the results. The classification applies only to the totally symmetric case $K(x,y)$0; general graphons, sparse-graph limits such as lattices, small-world networks, and power-law networks are not treated. Convergence of trajectories to individual fixed points in the infinite-dimensional model remains open, as does the genericity of convergence to asymptotically stable states (stable-manifold theory in infinite dimensions was not developed here). The Eulerian/Wasserstein framework lacks a second-order calculus at atomic measures, leaving open whether a boundary-appropriate stability theory could be formulated intrinsically. Finally, the numerical validation covers homogeneous random graphs at moderate size; heterogeneous topologies may deviate from the constant-$K(x,y)$1 heuristic of averaging interaction strengths.

Conclusion

This paper places OIMs within the optimal-transport/gradient-flow framework, deriving Lagrangian ($K(x,y)$2) and Eulerian (2-Wasserstein) mean-field models that preserve the dissipative structure of the finite system. In the symmetric case it delivers a complete equilibrium classification—at most four atomic clusters—and proves that, for almost all coupling ratios, only binarized fixed points can be directionally stable, with tight thresholds $K(x,y)$3. These results both explain empirically observed binarization behavior and delimit the parameter regime in which the hardware architecture can produce meaningful solutions, while leaving convergence questions and nonsymmetric graphons as the principal open problems.