2508.15234v2

# Bridging the Analog and the Probabilistic Computing Divide: Configuring Oscillator Ising Machines as P-bit Engines **Authors**: E.M.Hasantha Ekanayake, Nikhat Khan, Nikhil Shukla > University of Virginia, Charlottesville, VA, USA Abstract Oscillator Ising Machines (OIMs) and probabilistic bit (p-bit) platforms have emerged as promising non-Von Neumann paradigms for tackling hard computational problems. While OIMs realize gradient-flow dynamics, p-bit platforms operate through stochastic sampling. Although traditionally viewed as distinct approaches, this work presents a theoretically grounded framework for configuring OIMs as p-bit engines. We demonstrate that this functionality can be enabled through a novel interplay between first- and second harmonic injection to the oscillators. Our work identifies new synergies between the two methods and broadens the scope of applications for OIMs beyond combinatorial optimization problems to those that entail stochastic sampling. We further show that the proposed approach can be applied to other analog dynamical systems, such as the Dynamical Ising Machine. preprint: APS/123-QED I Introduction The endeavor to devise efficient solutions to complex computational problems has been a longstanding focus of science and technology research owing to its far-reaching implications for practical applications. A particularly promising direction is the design of special-purpose hardware that accelerates such tasks while improving energy efficiency, especially for problems that remain challenging for conventional digital platforms. Within this context, two emerging approaches have attracted significant attention: analog oscillator Ising machines (OIMs) and probabilistic bit (p-bit) based computing engines– the focus of the present work. OIMs, first proposed by Wang et al. [1] exploit the elegant equivalence between the ’energy function’ characterizing the dynamics of a network of coupled oscillators under second harmonic injection (SHI) and the Ising Hamiltonian. Minimizing the Ising Hamiltonian is an archetypal combinatorial optimization problem (COP) where the goal is to find spin configurations $s∈\{-1,+1\}$ that minimizes the Ising Hamiltonian given by $H=-\sum{J_{ij}s_{i}s_{j}}$ , where $J_{ij}$ is the interaction between spin $i$ and $j$ . OIMs realize gradient-flow dynamics on a continuous relaxation of the Ising model and have been extensively explored for solving combinatorial optimization problems (COPs) [2, 3]. They have been experimentally realized across diverse platforms, including optical [4, 5], acoustic [6], electronic [7, 8, 9, 10, 11, 12], spin-wave [13], and quantum systems [14]. These demonstrations have been complemented by extensive theoretical studies [15, 16, 17]. P-bit-based computing engines, on the other hand, are a complementary computing paradigm uniquely capable of Boltzmann sampling. A p-bit can be viewed as a tunable random number generator whose output probability depends on the synaptic input. At the network level, this yields a binary stochastic neural network (BSNN) [18], with spin updates governed by, $$ s_{i}^{+}=\text{sgn}\left[\tanh\left(\beta\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}s_{j}\right)-\mu\right] \tag{1} $$ where, $\mu$ is a random number typically selected from a uniform distribution between $[-1,1]$ , and $\beta$ is the equivalent of inverse temperature [19]. Furthermore, if we interpret the spins as the phases of oscillators—a perspective that is useful in the context of this work—then the state update rule can be expressed as, $$ s_{i}^{+}=\cos\phi_{i}=\text{sgn}\left[\tanh\left(\beta\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\cos\phi_{j}\right)-\mu\right] \tag{2} $$ where, $\phi∈\{0,\pi\}$ (wrapped phase form). In both Eqs. (1) and (2), the self bias term has not been considered although it is relatively straightforward to include it into the approach presented here. The intrinsic stochastic sampling capability of p-bit engines offers a complementary approach to OIMs in the context of solving COPs [20, 21, 22, 23, 24, 25, 26, 27]. In addition, the stochastic sampling can be exploited in other application such as probabilistic inference and learning. Establishing such a capability within OIMs—hitherto not demonstrated—could therefore significantly broaden their functional scope beyond combinatorial optimization. To date, OIMs and p-bit engines have largely been pursued as independent approaches, with limited exploration of their potential synergies [28, 29]. In this work, we bridge this gap by demonstrating the feasibility of configuring OIMs as p-bit engines, thereby establishing them as an alternative analog platform for probabilistic computing. <details> <summary>Figure1r.png Details</summary>  ### Visual Description ## Multi-Panel Scientific Diagram: Neuronal Dynamics and Energy Landscapes ### Overview The image presents a multi-panel scientific diagram analyzing neuronal dynamics through harmonic injections, phase trajectories, and energy landscapes. Key elements include: - First/Second Harmonic Injection (FHI/SHI) mechanisms - Phase-space trajectories with synaptic noise - Energy landscapes as functions of phase (ε) and gain (γ) - Stability analysis across parameter space ### Components/Axes #### Panel (a): Harmonic Injection Mechanism - **Diagram Elements**: - Circular trajectory labeled "Trajectory of Θ,ε" - Arrows indicating: - Noise: ±γ (orange) - Synaptic input: +γ (blue) - Labels: - FHI: First Harmonic Injection (green) - SHI: Second Harmonic Injection (orange) - Axes: Θ (horizontal), ε (vertical) #### Panel (b): Phase Trajectory Heatmap - **Axes**: - x-axis: γ (gain parameter, -0.5 to 0.5) - y-axis: ε(π) (phase, -0.5 to 0.5) - **Color Scale**: - Blue (-0.5) to Yellow (0.5) representing -∇E = ε - **Key Features**: - Vertical dashed line at γ = 0 - Yellow arrows marking: - ε = 0 (horizontal equilibrium) - -∇E = ε (diagonal equilibrium) - Ks = 0.15 (critical coupling parameter) #### Panel (c): 3D Energy Landscape - **Axes**: - x-axis: ε(π) (-0.5 to 0.5) - y-axis: γ (-0.2 to 0.2) - z-axis: Energy (a.u., -0.25 to 0.1) - **Color Gradient**: - Blue (low energy) to Yellow (high energy) - **Key Feature**: - Saddle-like structure with energy minima/maxima #### Panel (d): Energy vs Phase Curves - **Axes**: - x-axis: ε(π) (-0.5 to 0.5) - y-axis: Energy (a.u., -0.3 to 0.15) - **Data Series** (Ks = 0.15): - Orange: γ = 0 (baseline) - Red: γ = -0.1 (downward slope) - Purple: γ = -0.15 (steepest descent) - Green: γ = -0.2 (most negative slope) - Blue: γ = 0.1 (ascending slope) - Pink: γ = 0.15 (shallow ascent) - Green: γ = 0.2 (steepest ascent) ### Detailed Analysis #### Panel (a) - Circular trajectory shows phase-amplitude coupling - Noise (γ) modulates trajectory width - FHI/SHI labels suggest harmonic perturbation mechanisms #### Panel (b) - Phase trajectories cluster near: - ε = 0 (horizontal equilibrium) - -∇E = ε (diagonal stability boundary) - Ks = 0.15 indicates moderate coupling strength #### Panel (c) - Energy landscape reveals: - Double-well structure for γ < 0 - Monotonic increase for γ > 0 - Critical point at γ = 0, ε = 0 #### Panel (d) - Energy curves demonstrate: - γ < 0: Energy decreases with increasing ε - γ > 0: Energy increases with ε - Steeper slopes for |γ| > 0.15 ### Key Observations 1. **Bistability**: Panel (c) shows distinct energy minima for negative γ values 2. **Critical Coupling**: Ks = 0.15 appears in both (b) and (d), suggesting parameter consistency 3. **Slope Correlation**: Panel (d) confirms steeper energy gradients for |γ| > 0.15 4. **Equilibrium Points**: Panel (b) identifies two stable equilibria (ε=0 and -∇E=ε) ### Interpretation This diagram illustrates how synaptic noise (γ) and phase (ε) interact to shape neuronal energy landscapes. The FHI/SHI mechanism (a) likely represents input perturbations that drive the system toward different stability regimes. The heatmap (b) reveals how gain (γ) modulates phase trajectories, with Ks = 0.15 marking a critical coupling threshold. The 3D landscape (c) visualizes energy minima/maxima, while panel (d) quantifies how γ affects energy-phase relationships. Notably, negative γ values create bistable energy wells, suggesting potential for memory-like states in neuronal dynamics. The consistent Ks value across panels implies a unified parameter space for analyzing these interactions. </details> Figure 1: Dynamics of a harmonic oscillator under SHI (a) Schematic illustration of the signals required to program an harmonic oscillator as a p-bit. (b) Force field as a function of $\gamma$ ( $K_{s}=0.15$ , where $K_{s}$ is the strength of SHI.). (c) Corresponding energy landscape demonstrating its evolution with the synaptic input, $\gamma$ . (d) Specific cuts of the energy landscape at $\gamma=\{0,\,±\,0.1,\,±\,0.15,\,±\,0.2\}$ ( $K_{s}=0.15$ ). To establish the foundation of our oscillator-based p-bit engine, we begin by demonstrating two key properties of harmonic oscillators (the class of harmonic oscillators considered in this study) and their networks operating under SHI: (a) An oscillator subjected to first- and second harmonic injection can function as a binary stochastic neuron (BSN); and (b) A network of such coupled oscillators—specifically an OIM—can operate as a binary stochastic neural network (BSNN). These two properties form the conceptual and mathematical basis for designing an oscillator-based p-bit engine. II Configuring Oscillators as Stochastic Neurons To design a binary stochastic neuron (BSN) using an oscillator, we first examine the dynamics of a harmonic oscillator subjected to two external inputs (Fig. 1 (a)). (a) The first input is a signal with a frequency nearly equal to the oscillator’s natural frequency—a condition commonly referred to as injection locking. Within a specific locking range, this external signal entrains the oscillator, steering its output toward the phase and frequency of the injected signal, effectively synchronizing the oscillator’s behavior with the input. An energetics-based explanation of the relevant phase behavior is provided in Appendix A. We refer to this signal as the fundamental harmonic injection (FHI). All phase values are defined with respect to a common reference signal. (b) The second input is SHI, operating at twice the oscillator’s natural frequency. SHI drives the oscillator phase toward either $\phi=0$ or $\phi=\pi$ . The resulting phase dynamics of such an oscillator can be described by, $$ \begin{split}\frac{d\phi}{dt}&=-K_{c}\sin{\left(\phi-\theta\right)}-K_{s}\sin(2\phi)\\ \end{split} \tag{3} $$ where, $\phi$ denotes the output phase of the oscillator, while $\theta$ represents the phase offset of the FHI signal input. The parameters $K_{c}$ and $K_{s}$ are coupling constants of the FHI and SHI signals, with the first and the second terms on the RHS of Eq. (LABEL:neuron1) capturing the influence of the FHI and SHI, respectively. In this work, we will restrict our attention to cases where $\theta∈\{0,\frac{\pi}{2},\pi\}$ . Initially, we will focus on the analysis where $\theta∈\{0,\pi\}$ , with the case $\theta=\frac{\pi}{2}$ becoming relevant further on. Under this constraint, we can recast Eq. (LABEL:neuron1) as, $$ \begin{split}\frac{d\phi}{dt}=-\gamma\sin{\left(\phi\right)}-K_{s}\sin(2\phi)\end{split} \tag{4} $$ where $\gamma\,(=± K_{c})$ denotes the scaled synaptic input; $\gamma=K_{c}$ when $\theta=0$ and $\gamma=-K_{c}$ when $\theta=\pi$ . For simplicity, we will generally refer to $\gamma$ as the synaptic input in the following discussion. The details of the scaling factor will be elaborated in the subsequent section. Furthermore, for convenience, we perform a frame rotation by $0.5\pi$ , redefining the phase as $\phi=\frac{\pi}{2}+\epsilon$ . We also note that the relevant values of $\theta$ for this rotated frame shift to $\theta^{{}^{\prime}}∈\{-\frac{\pi}{2},0,\frac{\pi}{2}\}$ . Under this transformation, equation (4) becomes, $$ \frac{d\epsilon}{dt}=-\gamma\cos{\left(\epsilon\right)}+K_{s}\sin(2\epsilon) \tag{5} $$ We now analyze the properties and the behavior of Eq. (5). First, the fixed points of this system lie at $\epsilon_{1}^{*}=± 0.5\pi$ , and at values satisfying $\sin(\epsilon_{2}^{*})=\frac{\gamma}{2K_{s}}$ . A detailed stability analysis of all the fixed points has been presented in Appendix B. Next, we investigate the expected dynamical behavior of the system by analyzing the effective force field $(-∇ E=\dot{\epsilon})$ driving the phase evolution as a function of $\gamma$ , as shown in Fig. 1 (b). The dotted line in the figure represents the set of phase points where the phase velocity vanishes, i.e., $\frac{d\epsilon}{dt}=0$ , and is described by the equation: $$ -\gamma+2K_{s}\sin(\epsilon)=0 $$ This curve defines the nullcline of the system, separating regions of positive and negative phase flow. Interestingly, when the system is initialized at $\epsilon=0$ , the direction of phase evolution depends entirely on the sign of $\gamma$ . For $\gamma<0$ , the dynamics flow toward $\epsilon=\frac{\pi}{2}$ (corresponding to $\phi=\pi$ ). Conversely, for $\gamma>0$ , the dynamics flow toward $\epsilon=-\frac{\pi}{2}$ (i.e., $\phi=0$ ). Moreover, the magnitude of $\gamma$ determines the steepness of the phase flow at $\epsilon=0$ . At the critical point $\gamma=0$ , the direction of phase evolution becomes entirely stochastic in the presence of noise. Figure 1 (c) shows the evolution of the oscillator’s energy landscape with the synaptic input, $\gamma$ , with Fig. 1 (d) showing cuts at specific values of $\gamma=\{0,\,±\,0.1,\,±\,0.15,\,±\,0.2\}$ . In Fig. 1 (d), the relative symmetry about $\epsilon=0$ is evident across all cases, with the energy profiles for $\gamma=\{+0.1,\,+0.15,\,+0.2\}$ and $\gamma=\{-0.1,\,-0.15,\,-0.2\}$ appearing as mirror images, respectively, and the $\gamma=0$ case exhibiting a perfectly symmetric landscape. In alignment with prior work on magnetic tunnel junctions (MTJs)-based p-bits [30], we engineer our oscillator-based BSN to operate within the low energy barrier regime. Since the height of the barrier is controlled by $K_{s}$ (see Appendix C for details), this is achieved by tuning $K_{s}$ to a small value ( $K_{s}\ll 1$ ) during the sampling event. Thus, the oscillator-based approach can enable a BSN with a tunable barrier. From 1 (b) as well as from the form of Eq. (5), it can be observed that the system exhibits a unique symmetry at $\epsilon=0$ , within the domain $\epsilon∈\left[-\frac{\pi}{2},\frac{\pi}{2}\right]$ , whereby the magnitude of the (scaled) synaptic input has a symmetric effect for $+\gamma$ and $-\gamma$ , but induces phase flows in opposite directions. This symmetry plays a critical role in enabling the oscillator to function as a BSN, as it defines a neutral point from which the system can stochastically evolve toward one of two stable states depending on the sign of the synaptic input. As we will show later, this stochastic response is also highly non-linear. In practical settings, the oscillator can be driven to this neutral point $\epsilon=0$ using an FHI signal with a phase offset of $\theta^{{}^{\prime}}=0$ (we refer to this input as $\text{FHI}^{0}$ ), and without applying the SHI. Thus, operating the oscillator as a BSN entails the following steps: (i) Set oscillator phase to the neutral point using $\text{FHI}^{0}$ . (ii) Remove $\text{FHI}^{0}$ followed by application of the synaptic input (also an FHI signal when considering a single oscillator) and SHI signal (applied as a ramp) to evaluate the (stochastic) neuron’s state—this constitutes a sampling event. As we will demonstrate later, in the OIM network, such synaptic input is effectively generated by other connected oscillators within the network. III Dynamics of an Oscillator-based BSN To quantitatively analyze the stochastic nonlinear response of the oscillator-based BSN, we begin with the oscillator dynamics described in Eq. (5). We first define the updated spin state in terms of $\epsilon$ , which, in the context of the continuous time dynamics considered here, is given by, $$ s^{+}=\text{sgn}\left(\cos(\phi^{+})\right)=-\text{sgn}\left(\sin(\epsilon^{+})\right) $$ where, $\phi^{+}$ and $\epsilon^{+}$ refer to the phase of the system at a small time increment $\Delta t→ 0$ after the sampling has been initiated. We now evaluate the solution to the dynamics presented in Eq. (5). Although deriving an explicit solution is challenging, Eq. (5) has an implicit analytical solution given by, $$ \frac{\zeta(\epsilon)-\gamma\tanh^{-1}(\sin(\epsilon))}{\gamma^{2}-4K_{s}^{2}}=t \tag{6} $$ where, $$ \begin{split}\zeta(\epsilon)&=K_{s}\big(-2\log(\gamma-2K_{s}\sin(\epsilon))\\ &+\log(1-\sin(\epsilon))+\log(\sin(\epsilon)+1)\big)\\ \end{split} \tag{7} $$ For it’s applications as a BSN, we focus on the direction of the initial phase trajectory. To understand and evaluate the system dynamics at the sampling instant—defined as the moment when (FHI 0) is suppressed and the SHI and the synaptic input are asserted —we adopt the following approximation: (i) The dynamics are evaluated in the limit $t→ 0$ . (ii) We model the noise as Gaussian white noise with $\langle\eta(t)\rangle=0,\quad\langle\eta(t)\eta(t^{\prime})\rangle=2K_{n}\delta(t-t^{\prime}),$ with $K_{n}$ denoting the noise intensity i.e., $\eta(t)\,dt=\sqrt{2K_{n}}\,dW_{t},$ where $W_{t}$ is a Wiener process. Equivalently, over a finite timestep $\Delta t$ , $∈t_{t}^{t+\Delta t}\eta(\tau)\,d\tau\;\sim\;\mathcal{N}\!\big(0,\,2K_{n}\,\Delta t\big).$ (iii) At the onset of sampling ( $t→ 0$ ) , we assume $K_{s}→ 0$ such that $K_{s}\ll|\gamma|$ . This reflects the requirement for a low energy barrier in the probabilistic regime. Subsequently, $K_{s}$ must be ramped up for reasons discussed in the following section. Under these approximations, the oscillator phase can be expressed as: $$ \begin{split}\sin(\epsilon_{i})&\approx\tanh\left(-\gamma t+\epsilon^{\eta}(t)\right)\\ \\ \ &\approx\tanh\left(-\gamma t\right)+\epsilon^{\eta}(t).\text{sech}^{2}(\gamma t)\\ \\ \end{split} \tag{8} $$ Here, $\epsilon^{\eta}(t)\sim\mathcal{N}\!\big(0,\,2K_{n}\,t\big)$ , and is small such that the $\tanh(.)$ term can be linearized. Accordingly, at a small time increment $\Delta t→ 0$ after the onset of sampling at t=0, the updated state of the oscillator-based BSN can be described as, $$ \begin{split}s^{+}&\approx\text{sgn}\left[\tanh(\gamma\Delta t)-\epsilon^{\eta}(\Delta t).\text{sech}^{2}(\gamma\Delta t)\right]\\ \\ &\equiv\text{sgn}\left[\tanh(\gamma\Delta t)-\vartheta\right]\end{split} \tag{9} $$ where, $\vartheta\equiv\epsilon^{\eta}(\Delta t).\text{sech}^{2}(\gamma\Delta t)$ . Since $\text{sech}^{2}(\gamma\Delta t)∈(0,1]$ , and noise intensity is assumed to be small, $\vartheta$ has a high probability of being in the interval [-1,+1]. The infinitesimal time ( $\Delta t$ ) can be interpreted as the time duration required by the oscillator phase to reach a critical threshold magnitude $|\epsilon_{th}|$ , beyond which the oscillator phase cannot ‘reverse’ trajectory. $\Delta t$ can be tuned through properties of the SHI pulse, such as its slew rate. Additional details on $\Delta t$ are discussed in the following sections. The derivation of Eqs. (8) and (9) has been presented in Appendix D. Furthermore, while Eq. (8) addresses the regime where $\gamma\gg K_{s}$ , we also consider (in Appendix D) the complementary case where both $\gamma$ and $K_{s}$ are small and comparable. In this setting, the synaptic bias $\gamma$ and the stochastic perturbations act as competing drivers of the phase dynamics. A large $|\gamma|$ produces a predictable exponential drift, whereas strong noise leads to rapid amplification of fluctuations and broad dispersion of trajectories. The observed evolution is therefore governed by the balance between deterministic drive and stochastic forcing. We note that the actual synaptic input to the BSN would be applied as the voltage (or current) amplitude, $V_{\text{inj}}$ , and phase (in-phase or out-of-phase) of the injected signal, which relate to $\gamma$ and the coupling constant $K_{c}$ as $\gamma=±\,K_{c}\,≈±\,\,\frac{\omega_{0}}{2Q}·\frac{V_{\text{inj}}}{V_{\text{osc}}}$ [31, 32], where $Q$ is the quality factor of the oscillator, $\omega_{0}$ is the natural frequency, and $V_{\text{osc}}$ represents the amplitude of the oscillator, respectively. This formulation implies that $\gamma$ serves as a scaled representation of the synaptic input, as noted earlier. Furthermore, the effective inverse temperature is then given by $\beta_{\text{eff}}=\frac{\omega_{0}\Delta t}{2V_{\text{osc}}}.\frac{1}{Q}$ implying that $\beta_{\text{eff}}$ can be modulated using the quality factor, $Q$ , of the oscillator as well as other parameters such as the oscillation amplitude among others. This tunability provides a practical mechanism for controlling the stochastic behavior of the system. The updated state is then expressed as $$ \begin{split}s^{+}&\approx\text{sgn}\left[\tanh\left(\frac{\omega_{0}\Delta t}{2QV_{\text{osc}}}\cdot V_{\text{inj}}\right)-\vartheta\right]\\ \\ &\equiv\text{sgn}\left[\tanh\left(\beta_{\text{eff}}\cdot V_{\text{inj}}\right)-\vartheta\right]\end{split} \tag{10} $$ <details> <summary>Figure2r.png Details</summary>  ### Visual Description ## Chart/Diagram Type: Probability vs. Injected Voltage with Fit Parameters ### Overview The image contains two panels (a) and (b) depicting probability distributions as a function of injected voltage (V_inj). Both panels include: - Sigmoidal curves fitted with the equation: $ p = \frac{1 + \tanh(\beta_{\text{eff}} V_{\text{inj}})}{2} $ - Tables listing numerical parameters (K_n, β_eff) - Legends mapping colors to K_n values - Key annotations for maximum curvature (K_s,max) and fit parameters ### Components/Axes #### Panel (a) - **X-axis**: V_inj (ranging from -0.8 to 1.0) - **Y-axis**: Probability, p (0.0 to 1.0) - **Legend**: Located at top-right, colors correspond to K_n values (0.01, 0.05, 0.10, 0.15, 0.20) - **Fit Equation**: Annotated with $ K_{s,\text{max}} = 0.15 $ - **Table**: Embedded in the plot, showing K_n and β_eff values: | K_n | β_eff (from fit) | |-------|------------------| | 0.01 | 88.76 | | 0.05 | 12.631 | | 0.10 | 6.289 | | 0.15 | 3.837 | | 0.20 | 2.704 | #### Panel (b) - **X-axis**: V_inj (ranging from -1.0 to 1.0) - **Y-axis**: Probability, p (0.0 to 1.0) - **Legend**: Located at top-right, colors correspond to β_eff ≈ 8, 6, 4, 2 - **Fit Equation**: Annotated with $ K_{s,\text{max}} = 0.15 $, $ K_n = 0.15 $ - **Table**: Embedded in the plot, showing β_eff approximations: | β_eff ≈ | Color | |---------|--------| | 8 | Purple | | 6 | Orange | | 4 | Green | | 2 | Blue | ### Detailed Analysis #### Panel (a) - **Curves**: Five sigmoidal curves (dark blue to olive green) represent increasing K_n values. Higher K_n curves are steeper and shift rightward. - **Fit Equation**: The tanh function models probability as a function of V_inj, with β_eff inversely proportional to K_n (e.g., β_eff = 88.76 for K_n = 0.01 vs. 2.704 for K_n = 0.20). - **K_s,max**: Annotated at V_inj = -0.15, indicating the voltage at maximum curvature for all curves. #### Panel (b) - **Curves**: Four sigmoidal curves (purple to dark blue) represent decreasing β_eff values. Lower β_eff curves are shallower and shift leftward. - **Fit Equation**: Same tanh form as panel (a), but β_eff is fixed at ≈8, 6, 4, 2, with K_n = 0.15 constant. - **K_s,max**: Same value (0.15) as panel (a), but curves are less steep due to fixed K_n. ### Key Observations 1. **Inverse Relationship**: Higher K_n values (panel a) correlate with lower β_eff, reducing the sensitivity of probability to V_inj. 2. **Steepness vs. Spread**: Panel (a) shows steeper curves for smaller K_n, while panel (b) shows broader distributions for lower β_eff. 3. **Consistency**: Both panels use the same fit equation, but panel (b) fixes K_n = 0.15, altering the parameter space. 4. **Color-Legend Matching**: All curves in both panels align with their respective legend colors (e.g., dark blue = K_n = 0.20 in panel a; purple = β_eff ≈8 in panel b). ### Interpretation The data demonstrates how probability distributions of an injected voltage-dependent process are modulated by two parameters: 1. **K_n**: Controls the steepness of the sigmoidal curve. Smaller K_n values (e.g., 0.01) produce sharper transitions, while larger K_n (e.g., 0.20) yield broader distributions. 2. **β_eff**: Represents the effective voltage sensitivity. Higher β_eff (e.g., 88.76) amplifies the effect of V_inj, whereas lower β_eff (e.g., 2.704) dampens it. The fit equation $ p = \frac{1 + \tanh(\beta_{\text{eff}} V_{\text{inj}})}{2} $ suggests a logistic-like behavior, where probability transitions from 0 to 1 as V_inj crosses zero. The tables confirm that β_eff scales inversely with K_n (panel a) and is discretized in panel (b). The consistent K_s,max = 0.15 across panels implies a universal curvature point, independent of K_n or β_eff. This could reflect a physical constraint (e.g., maximum response rate) in the system being modeled. </details> Figure 2: Oscillator-based BSN. Firing probability (symbols) as a function of the synaptic input for: (a) varying levels of noise ( $K_{n}$ ). (b) different values of $\beta_{\text{eff}}$ derived for $K_{n}=0.15$ . We note that the $\beta_{\text{eff}}$ profile will change for a different $K_{n}$ . The lines in the plot indicate fits using the equation $p=\frac{1+\tanh(\beta_{\text{eff}}V_{\text{inj}})}{2}$ along with the calculated $\beta_{\text{eff}}$ ; $p:$ firing probability. All fits exhibit $R^{2}>0.999$ Equation (LABEL:final_BSN) showcases the BSN’s capability to perform Boltzmann sampling, and consequently, function as a p-bit when initialized at the critical phase point, $\epsilon=0$ . In the following section, we show that, in OIMs, $V_{\text{inj}}$ —and consequently $\gamma$ —is generated and regulated through the mutual feedback among the coupled oscillators. To validate the dynamics derived above, we simulate the oscillator-based BSN’s switching using a stochastic differential equation solver implemented in MATLAB ®. We first examine the evolution of stochastic behavior under varying noise levels. The phase is initialized at $\epsilon=0$ and allowed to evolve in the presence of different noise intensities and synaptic inputs. The firing probability is then estimated over 2000 such cycles. Figure 2 (a) presents the simulation results (symbols) for the output firing probability of the oscillator-based BSN as a function of $V_{\text{inj}}$ . These results are fitted using the function $p=\frac{1+\tanh(\beta_{\text{eff}}V_{\text{inj}})}{2}$ (solid lines), showing excellent agreement with the simulated data ( $R^{2}>0.999$ ); p is the probability of the neuron firing i.e., switching to $s=+1$ . In practical implementations, modulating external noise to control temperature and stochasticity may not be feasible. Instead, the effective inverse temperature—and thus the degree of stochastic behavior—can be tuned by adjusting $\beta_{\text{eff}}$ , which is achievable through modulation of the oscillator’s quality factor $Q$ . Figure 2 (b) illustrates the evolution of firing probability (symbols) for different values of $\beta_{\text{eff}}$ , where the switching probability again exhibits tanh(.) dependence on $V_{\text{inj}}$ as confirmed by the corresponding fitted curves (solid lines) with $R^{2}>0.999$ . A noise strength of $K_{n}=0.15$ was used in the simulation. These simulation results further support that the oscillator exhibits the characteristic behavior of a BSN capable of performing Boltzmann sampling. A key consideration in engineering the oscillator’s dynamics for BSN functionality is the relative magnitude of $\gamma$ and the SHI strength, $K_{s}$ . Realizing effective Boltzmann sampling behavior requires a small energy barrier, which corresponds to the regime $K_{s}→ 0$ such that $K_{s}\ll|\gamma|$ . If this condition is not satisfied, the system may still operate as a BSN; however, its dynamics may deviate from the Boltzmann sampling behavior. This is detailed in the analysis in Appendix D. In contrast, to preserve the oscillator’s phase trajectory after sampling and to enable reliable readout, a lower bound on $K_{s}$ must be satisfied. Specifically, this bound ensures that once the phase magnitude exceeds a certain threshold ( $|\epsilon_{\text{th}}|$ ), the phase continues to evolve in the same direction until it reaches the corresponding fixed point. While a detailed analytical derivation is provided in Appendix E, we offer here a qualitative explanation of the origin of this constraint by examining the dynamics described by Eq. (5). Within the interval $\epsilon∈(-\frac{\pi}{2},\frac{\pi}{2})$ , the cosine term satisfies $\cos(\epsilon)>0$ , and the synaptic input term $-\gamma\cos(\epsilon)$ therefore drives the phase evolution in the direction opposite to the sign of $\gamma$ . This implies that the synaptic input alone tends to push the phase toward the fixed point opposite to the sign of $\gamma$ , unless counteracted by the SHI term. If noise initially drives $\frac{d\epsilon}{dt}$ in the direction not favored by the synaptic input (Fig. 8 (b)), and $K_{s}$ is too small, the oscillator may reverse its trajectory to align with the trajectory favored by the synaptic input. In such cases, the final state may not reflect the oscillator’s initial trajectory or the intended output $s^{+}$ . The SHI term— specifically, the $K_{s}\sin(2\epsilon)$ component in Eq. (5), acts to reinforce the initial direction of $\frac{d\epsilon}{dt}$ generated by the stochastic sampling process, provided $K_{s}$ is sufficiently large. This reinforcement helps ensure that the phase continues toward the correct fixed point. The critical condition to ensure that an oscillator at $\epsilon=|\epsilon_{\text{th}}|$ maintains its trajectory is given by: $$ |\gamma|<2K_{s}\sin(|\epsilon_{\text{th}}|) $$ It is important to note that due to the presence of noise, this threshold is inherently probabilistic, resulting in a diffuse rather than deterministic boundary. This condition appears to contradict the earlier requirement concerning the relative magnitudes of $\gamma$ and $K_{s}$ . To reconcile these seemingly opposing constraints, we propose implementing the SHI input as a ramp signal with a carefully engineered slew rate. Specifically, the SHI can be initialized at a low amplitude—ensuring that $K_{s}\ll|\gamma|$ —to facilitate stochastic sampling in the low-barrier regime. Subsequently, the amplitude of the SHI is increased to reinforce the resulting phase trajectory. Additionally, as noted above the characteristics of the applied SHI input, particularly its slew rate, also modulate $\Delta t$ . A smaller slew rate increases $\Delta t$ , which in turn modifies the inverse temperature and also makes the system less stochastic (see Appendix F for details). The analysis presented in this section reveals that the competition between synaptic input and noise gives rise to a stochastic bifurcation, which constitutes the core operating principle of the proposed oscillator-based BSN. IV Configuring OIMs As Binary Stochastic Neural Networks Building upon the ability to configure an oscillator as a BSN capable of performing Boltzmann sampling, we now investigate the possibility of configuring an OIM as a p-bit engine, or in other words, a BSNN. The core premise of this idea is that after the system reaches steady state, i.e., $\epsilon∈\{-\frac{\pi}{2},\frac{\pi}{2}\}$ , the dynamics of a randomly sampled oscillator driven to the phase point $\epsilon=0$ (i.e., $\phi=\frac{\pi}{2}$ ) still approximate Boltzmann sampling. Steady state initialization at $\epsilon∈\{-\frac{\pi}{2},\frac{\pi}{2}\}$ can be achieved by applying an SHI signal prior to the onset of stochastic sampling. As noted earlier, in the OIM, the feedback from other coupled oscillators acts as the effective synaptic input, which subsequently modulates the oscillator’s stochastic dynamics. To establish this result, in the following sections, we divide the oscillators in the network into two categories and analyze their dynamics: (i) Randomly sampled oscillator $i$ initialized to $\mathbf{\epsilon_{i}=0}$ . The phase evolution of an oscillator in the OIM network can be described by the equation, $$ \frac{d\phi_{i}}{dt}=-K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\phi_{i}-\phi_{j})-K_{s}\sin(2\phi_{i})\\ $$ which, in the rotated frame of reference, can be expressed as, $$ \frac{d\epsilon_{i}}{dt}=-K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\epsilon_{i}-\epsilon_{j})+K_{s}\sin(2\epsilon_{i}) \tag{11} $$ As alluded to earlier, we assume that: (i) The selected oscillator is initialized at $\epsilon=0$ . This can be accomplished using $\text{FHI}^{0}$ signal with large amplitude. (ii) Since we begin performing stochastic sampling after the OIM network has achieved steady state, all other oscillators are at $\epsilon=±\frac{\pi}{2}\>\>\big(\phi∈\{0,\pi\}\big)$ . Furthermore, we will ensure that the oscillators maintain their state during the sampling event by applying a sufficiently large $K_{s}$ . This condition is required to ensure that the oscillator dynamics map to Gibbs sampling (see Appendix G). In the subsequent sections, we will discuss how these conditions can be implemented. Under the constraints outlined above, the phase dynamics of the selected oscillator simplify to: $$ \begin{split}\frac{d\epsilon_{i}}{dt}&=\left(K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\epsilon_{j})\right)\cos(\epsilon_{i})+K_{s}\sin(2\epsilon_{i})\\ \\ &\equiv-\gamma_{i}\cos(\epsilon_{i})+K_{s}\sin(2\epsilon_{i})\end{split} \tag{12} $$ where, $$ \gamma_{i}=-K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\epsilon_{j}) $$ We now show that $\gamma_{i}$ , as defined above, represents the synaptic input to oscillator $i$ from the other connected oscillators in the network. To establish this, we consider the relationship between $\epsilon$ and $\phi$ , and the fact that $s=\cos(\phi)$ when $\phi∈\{0,\pi\}$ : $$ \begin{split}\gamma_{i}=-K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\epsilon_{j})&=-K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin\left(\phi_{j}-\frac{\pi}{2}\right)\\ =K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\cos(\phi_{j})&=K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}s_{j}\end{split} $$ This establishes that $\gamma_{i}$ represents the net synaptic input received by oscillator $i$ from the connected oscillators under the conditions described above. Additionally, we note that the self-biasing term can be incorporated by injecting an FHI signal to oscillator with the strength and phase of the signal representing the self-bias input. The updated state of oscillator $i$ can then be expressed as, $$ \begin{split}s_{i}^{+}&\approx-\text{sgn}\left[\tanh(-\gamma_{i}\Delta t)+\epsilon^{\eta}(\Delta t).\text{sech}^{2}(\gamma\Delta t)\right]\\ \\ &\approx\text{sgn}\left[\tanh\left(K\Delta t\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}s_{j}\right)-\vartheta\right]\end{split} \tag{13} $$ which closely resembles the state update rule for p-bits (Eq. (1)), with the factor $\beta_{\text{eff}}=K\Delta t$ serving as the effective inverse temperature. As detailed in [1], the value of K depends on the perturbation projection vector function for the oscillator as well as the amplitude of the perturbation from the oscillators in the network, which can be tuned via the coupling element in the network. This equivalence implies that even in the OIM, the oscillator can approximate Boltzmann sampling thereby enabling the network to function as a p-bit platform. One of the critical requirements of realizing the above dynamics is to drive the randomly selected oscillator $i$ to $\epsilon_{i}=0\>\>\big(\phi_{i}=\frac{\pi}{2}\big)$ . This can be achieved by applying a large FHI signal, $\text{FHI}^{0}$ , while suppressing the SHI signal $(K_{s}=0)$ . The resulting dynamics can be described by: $$ \frac{d\epsilon_{i}}{dt}=-K_{c,i}\sin(\epsilon_{i})-K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\epsilon_{i}-\epsilon_{j}) \tag{14} $$ The largest magnitude of the second term is $D_{i}$ —the degree of the node (oscillator) $i$ in the network. By designing $K_{c,i}\gg KD_{i}$ ensures that the contributions of the second term are small. Consequently, the phase will be driven to $\epsilon_{i}≈ 0$ , thereby preparing it for the subsequent stochastic sampling event. (ii) Oscillators not being sampled Under steady state, such oscillators have phases $\epsilon=±\frac{\pi}{2}$ , and the goal is to maintain the configuration when oscillator $i$ is sampled. To achieve this, we apply a strong SHI signal. The corresponding dynamics for an oscillator $j$ that is not being sampled can be then described as follows: $$ \begin{split}\frac{d\epsilon_{j}}{dt}=-K\sum_{\begin{subarray}{c}k=1\\ k\neq j\end{subarray}}^{N}J_{jk}\sin(\epsilon_{j}-\epsilon_{k})+K_{s}\sin(2\epsilon_{j})\\ \\ \quad\forall j\in\{1,2,\dots,i-1,i+1,\dots,N\}\end{split} \tag{15} $$ Furthermore, $\sin(\epsilon_{j}-\epsilon_{k})≈ 0$ since $$ \epsilon_{k}\in\left\{-\frac{\pi}{2},\frac{\pi}{2}\right\}\quad\forall j,k\in\{1,2,\dots,i-1,i+1,\dots,N\} $$ Consequently, the dynamics can be reduced to, $$ \begin{split}\frac{d\epsilon_{j}}{dt}=-KJ_{ji}\sin(\epsilon_{j}-\epsilon_{i})+K_{s}\sin(2\epsilon_{j})\\ \\ \quad\forall j\in\{1,2,\dots,i-1,i+1,\dots,N\}\end{split} \tag{16} $$ The maximum magnitude of the first term is $|KJ_{ji}|$ . Thus, by using $2K_{s,j}\gg KJ_{ji}$ , the phase can be maintained at $\epsilon_{j}≈\{-\frac{\pi}{2},\frac{\pi}{2}\}$ . These assumptions are further validated through simulations presented in the following section. Based on the above analysis, algorithm 1 presents the scheme for configuring OIMs as p-bit engines. Although the analysis above was conducted in a rotated frame of reference, we present the operational scheme in terms of the original phase variable $\phi$ , to maintain consistency with the prevailing conventions in the OIM literature, where the spin states are defined as $\phi∈\{0,\pi\}$ . As previously noted, the relationship between the two frames is given by $\phi=\frac{\pi}{2}+\epsilon$ . Table 1: Operating OIM as a p-bit platform | 1: 2: 3: | Initialize all oscillators. Apply SHI signal to all oscillators. The OIM performs gradient descent and achieves a steady state characterized by $\phi∈\{0,\pi\}$ . while not converged or for a fixed number of iterations do | | --- | --- | | 4: | Randomly select an oscillator $i$ . | | {Step (i): Prepare oscillator for stochastic sampling} | | | 5: | Apply $\mathrm{FHI}^{0}$ signal to oscillator $i$ with appropriate amplitude and set SHI signal to oscillator $i$ to $0$ . | | {Step (ii): Initiate probabilistic evolution} | | | 6: | Reduce $\mathrm{FHI}^{0}$ signal to $0$ . | | 7: | Ramp SHI signal to oscillator $i$ . | | 8: | Let $\phi_{i}$ probabilistically relax to $\phi=0$ or $\phi=\pi$ based on synaptic feedback from connected nodes and intrinsic noise. | | 9: | end while | V Computing with OIMs configured as P-bit engines V.1 5-Node Adder <details> <summary>Figure_sampling.png Details</summary>  ### Visual Description ## Bar Chart: Boltzmann Law vs Measured (Oscillators) Probabilities ### Overview The chart compares two probability distributions: "Boltzmann law" (blue bars) and "Measured (oscillators)" (red bars) across discrete categories labeled [C_i B A S C_o]. The x-axis combines categorical labels and decimal values (0–31), while the y-axis represents probability (0–0.2). A table at the top maps categorical combinations to decimal values. ### Components/Axes - **X-axis**: - Labels: [C_i B A S C_o] (categorical) + decimal values (0, 6, 10, 13, 18, 21, 25, 31). - Categories correspond to combinations of binary variables (C_i, B, A, S, C_o) listed in the table. - **Y-axis**: - Label: "Probability" (0–0.2 in increments of 0.05). - **Legend**: - Position: Right of the chart. - Entries: - Blue: "Boltzmann law" - Red: "Measured (oscillators)" - **Table**: - Position: Top-left corner. - Columns: C_i, B, A, S, C_o, Decimal. - Rows: 8 combinations (e.g., 00000 → 0, 00110 → 6, etc.). ### Detailed Analysis - **Boltzmann Law (Blue Bars)**: - Values: - 0: ~0.12 - 6: ~0.13 - 10: ~0.12 - 13: ~0.13 - 18: ~0.12 - 21: ~0.11 - 25: ~0.13 - 31: ~0.12 - Trend: Relatively flat with minor fluctuations (max ~0.13, min ~0.11). - **Measured (Oscillators) (Red Bars)**: - Values: - 0: ~0.12 - 6: ~0.13 - 10: ~0.12 - 13: ~0.13 - 18: ~0.12 - 21: ~0.11 - 25: ~0.13 - 31: ~0.12 - Trend: Nearly identical to Boltzmann law, with slight deviations (e.g., 25: red bar slightly taller). - **Table Data**: - Decimal values map to binary combinations (e.g., 00000 → 0, 00110 → 6, 01010 → 10, etc.). ### Key Observations 1. **Alignment**: Both distributions show nearly identical probabilities across all categories, with red bars (measured) slightly exceeding blue bars (Boltzmann) at decimal 25 (~0.13 vs. ~0.12). 2. **Consistency**: Probabilities remain within a narrow range (0.11–0.13), suggesting minimal variability. 3. **Outliers**: No significant outliers; deviations are within measurement uncertainty. ### Interpretation The chart demonstrates strong agreement between theoretical predictions (Boltzmann law) and experimental measurements (oscillators), with discrepancies likely attributable to experimental noise or minor model approximations. The flat trend across categories implies the system’s behavior is stable under the tested conditions. The slight divergence at decimal 25 warrants further investigation to determine if it reflects a systematic error or a genuine deviation from the law. </details> Figure 3: Full adder Probability histogram measured using the oscillators (orange bars, obtained using $5× 10^{5}$ sweeps) compared with the target Boltzmann distribution (blue bars). Following the convention of Ref. [19], states are indexed by the decimal value of the binary word $[C_{\mathrm{in}}\;A\;B\;S\;C_{\mathrm{out}}]$ . The dominant peaks correspond to the valid entries of the full-adder truth table, as highlighted in the inset. The two distributions show excellent agreement, with a measured KL divergence of $7.68× 10^{-4}$ . Simulation parameters: $K=0.18$ ; $K_{s}(t)=K_{s,\text{max}}(1-e^{-\frac{t}{10^{-2}}})$ ; $K_{n}=0.1$ . To evaluate the sampling capability of our oscillator-based BSN, we benchmark its output distribution against the target Boltzmann distribution for a full adder, following the approach presented by Camsari et al. [19]. For this purpose, we construct the corresponding $14× 14$ adjacency matrix with the following assignment of nodes: 1–9 represent auxiliary and handle bits, 10 corresponds to $C_{\mathrm{in}}$ , 11 to $A$ , 12 to $B$ , 13 to the sum $S$ , and 14 to the carry-out $C_{\mathrm{out}}$ . The network is operated in the so-called truth table mode, where all input and output terminals are left floating. Under this condition, the distribution obtained using the oscillator-based BSN (using $5× 10^{5}$ sweeps) aligns closely with the target Boltzmann distribution, thereby confirming its ability to perform correct stochastic sampling. Here, we note that simulating the analog dynamics of all oscillators for $5× 10^{5}$ sweeps (corresponding to a total simulated time of $\sim 2× 10^{7}$ per oscillator) within the stochastic differential equation (SDE) framework (to account for noise) is computationally prohibitive. Therefore, we only simulate the dynamics of oscillator being sampled, while the remaining oscillators are held at their fixed phase values. However, in the following example—computing MaxCut—which requires a smaller number of sweeps, we simulate the entire oscillator network. Figure 3 presents the resulting probability histogram, showing an excellent match to the target Boltzmann distribution (KL divergence= $7.68× 10^{-4}$ ). Each state is labeled by the decimal value of the binary string $[\,C_{\mathrm{in}}\;\;A\;\;B\;\;S\;\;C_{\mathrm{out}}\,]$ and the inset highlights the valid truth-table states of the full adder, which appear as the tallest peaks in the distribution. V.2 Computing MaxCut <details> <summary>Figure3r.png Details</summary>  ### Visual Description ## Oscillator Sampling Sequence and Dynamics Analysis ### Overview The image presents a multi-part technical analysis of an oscillator system, combining sampling sequences, phase dynamics, and network partitioning. It includes three distinct components: (a) a sampling sequence with color-coded transitions, (b) phase dynamics over time, and (c) a network partitioning diagram. --- ### Components/Axes #### (a) Sampling Sequence - **Header**: "Oscillator sampling sequence:" in red text. - **Transition Arrows**: - Sequence: `5 → 6 → 15 → 7 → 1 → 3 → 9 → 4 → 11 → 8 → 14 → 12 → 13 → 10 → 2 → 14 → 9 → 12 → 8 → 5 → 11 → 3 → 15 → 10 → 13 → 6 → 4 → 7 → 2 → 1` - Arrows are blue, with numbers in black. - **Vertical Bars**: - **Top Set**: Labeled `FHI⁰ > 0, Kₛ = 0` (pink, black, blue, teal, purple, orange, brown, green, red, etc.). - **Bottom Set**: Labeled `FHI⁰ = 0, Kₛ > 0` (same color palette). - Colors are evenly distributed across the bars, with no explicit legend provided in the image. #### (b) Phase Dynamics (φ(π) vs. Time) - **Axes**: - **X-axis**: "Time" (0–80, linear scale). - **Y-axis**: "φ(π)" (range: -1 to 1, with a secondary axis labeled "40" at the bottom). - **Data Series**: - Multiple colored lines (pink, black, blue, teal, purple, orange, brown, green, red) representing different states or parameters. - Arrows point to specific time points (e.g., ~50, ~60, ~70) with annotations like "↓" and "↓". - **Legend**: Not explicitly visible in the image, but colors correspond to distinct data series. #### (c) Network Partitioning (OIM) - **Title**: "OIM" in red text. - **Graph**: - Nodes labeled 1–15, connected by blue lines. - A vertical green "cut" line at **time = 35**, intersecting nodes 1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15. - Nodes are densely interconnected, forming a complex network. - **Annotations**: - "Cut 35" labeled on the y-axis (left side). --- ### Detailed Analysis #### (a) Sampling Sequence - The sequence shows a cyclical pattern of transitions between 15 states (1–15), with some states revisited (e.g., 14 → 9 → 12 → 8 → 5 → 11 → 3 → 15). - The two sets of vertical bars likely represent two experimental conditions: - **FHI⁰ > 0, Kₛ = 0**: Higher variability in transitions (more color diversity). - **FHI⁰ = 0, Kₛ > 0**: More uniform transitions (consistent color distribution). #### (b) Phase Dynamics - **Trends**: - All lines oscillate around φ(π) = 0, with some dipping below -1 (e.g., purple line at ~50). - The orange line shows the most pronounced dips, while the black line remains relatively stable. - Arrows highlight critical transitions (e.g., a sharp drop at ~50 and a recovery at ~70). - **Uncertainty**: No numerical values provided; trends are inferred from visual inspection. #### (c) Network Partitioning - The "cut" at time = 35 divides the network into two subgraphs: - **Left Subgraph (time < 35)**: Nodes 1–15 are fully connected. - **Right Subgraph (time > 35)**: Nodes 1–15 remain connected, but the cut line suggests a threshold for partitioning. - **OIM Label**: Likely refers to "Optimal Information Measure" or a similar metric, indicating a critical point in the network's structure. --- ### Key Observations 1. **Sampling Sequence**: The oscillator cycles through 15 states with a repeating pattern, suggesting periodic or quasi-periodic behavior. 2. **Phase Dynamics**: The system exhibits oscillatory behavior with varying stability across parameters (FHI⁰, Kₛ). 3. **Network Partitioning**: The cut at time = 35 may represent a critical threshold for network connectivity or information flow. --- ### Interpretation - **System Behavior**: The oscillator's dynamics (φ(π)) are influenced by parameters FHI⁰ and Kₛ. The sampling sequence and phase dynamics suggest a balance between stability and variability. - **Network Criticality**: The OIM cut at time = 35 implies a phase transition or bifurcation in the network's structure, potentially marking a shift in system behavior. - **Uncertainty**: The absence of numerical data limits precise quantification of trends. Further analysis would require explicit values for φ(π) and transition probabilities. </details> Figure 4: Operating OIMs as p-bit platforms. (a) Randomly generated FHI sequence to the oscillators. The application of FHI to the oscillator is accompanied by the suppression of SHI, and vice-versa. (b) Phase response of the oscillators over time. (c) Evolution of the computed graph cut over time/iterations. With stochastic sampling, the system is able to reach the globally optimal solution (MaxCut =39). The red arrows highlight sampling events that improved cut. Simulation parameters: $K=1;K_{s,\text{max}}=2;\>\text{FHI}^{0}=50;\>K_{n}\text{ for sampled oscillator}=0.04$ ; ramping schedule of the SHI signal: $K_{s}(t)=K_{s,\text{max}}\big(1-e^{-\frac{t}{\tau}}\big)$ , where $K_{s,\text{max}}=2$ and $\tau=10^{-2}$ . We now evaluate the archetypal MaxCut problem using the OIM’s stochastic sampling mode. Computing the MaxCut of a graph is a NP-hard problem where the objective is to partition the nodes such that the weight of the edges shared among the two sets (i.e., intersect the cut) is maximized. The MaxCut problem directly maps to the solution of the corresponding anti-ferromagnetic Ising Hamiltonian i.e., $J_{ij}=-W_{ij}$ , where $\>W_{ij}$ represents the weight of the edges in the graph to be partitioned. While the MaxCut problem has been used as a representative example, the above approach can be applied to other COPs such as graph coloring and finding maximum independent sets, among others. We demonstrate the functionality using a randomly generated graph with 15 nodes and 59 edges. Figure 4 (a) illustrates the sequence of $\text{FHI}^{0}$ inputs applied to various randomly sampled oscillators in a sequential manner. The assertion of $\text{FHI}^{0}$ is accompanied by suppression of the SHI signal to that oscillator. The resulting phase dynamics of the oscillators are shown in Fig. 4 (b). It can be observed that during each sampling event, the sampled oscillator is first driven to a phase of $\phi_{i}=\frac{\pi}{2}$ by the application of the $\text{FHI}^{0}$ signal, while the phases of the other (non-sampled) oscillators remain at $\phi∈\{0,\pi\}$ . Subsequently, the $\text{FHI}^{0}$ signal is suppressed, and the corresponding SHI input is reasserted (not shown in Fig. 4 (a) for clarity). As observed in the dynamics, the sampled oscillator then stochastically relaxes toward either $\phi_{i}=0$ or $\phi_{i}=\pi$ , with the direction of phase relaxation determined by the competition between the synaptic input and the noise, as noted earlier. Arrows in Fig. 4 (b) indicate representative cases where the oscillator flips its state. It is noteworthy that in traditional OIMs, such transitions are likely to have a low probability of occurrence once all oscillator phases have settled to $\phi_{i}=0$ or $\phi_{i}=\pi$ since $\frac{d\phi}{dt}$ , which effectively represents the diving force ( $-∇ E=\frac{d\phi}{dt}$ ) on the oscillator phase is close to zero for all oscillators i.e., $\frac{d\phi}{dt}≈ 0$ . Figure 4 (c) shows the corresponding evolution of the graph cut. The red arrows in Figs. 4 (c) highlight the sampling events that lead to an increase in the graph cut allowing the system to compute the optimal MaxCut. Simulation results on 10 additional randomly generated graphs have been presented in Appendix H. <details> <summary>Figure_MaxCutBZ.png Details</summary>  ### Visual Description ## Scatter Plots: Relationship Between Energy, ACF, and β Parameter ### Overview The image contains two scatter plots labeled (a) and (b), analyzing relationships between variables with three distinct β (beta) parameter values. Both plots use color-coded data points (purple, orange, green) to represent different β values, with high R² values indicating strong linear or autocorrelation patterns. --- ### Components/Axes #### Plot (a) - **X-axis**: Energy (ranging from -20 to 10) - **Y-axis**: log(p) (ranging from -15 to 0) - **Legend**: - Purple: β = 1.145 (R² = 0.99) - Orange: β = 0.662 (R² = 0.977) - Green: β = 0.362 (R² = 0.99) - **Trend**: All three lines show a steep downward slope, with purple (β=1.145) being the steepest. #### Plot (b) - **X-axis**: lag (ranging from 0 to 30) - **Y-axis**: ACF (Autocorrelation Function, ranging from 0 to 1) - **Legend**: - Purple: β = 1.145 - Orange: β = 0.662 - Green: β = 0.362 - **Trend**: All three lines show a gradual decline, with purple (β=1.145) starting highest and green (β=0.362) ending lowest. --- ### Detailed Analysis #### Plot (a) - **Data Points**: - Purple (β=1.145): Points tightly clustered along a steep linear regression line (R²=0.99). - Orange (β=0.662): Slightly less steep slope (R²=0.977), with minor deviations. - Green (β=0.362): Flattest slope (R²=0.99), indicating near-perfect linear fit despite lower β. - **Spatial Grounding**: Legend positioned in the top-right corner of the plot. #### Plot (b) - **Data Points**: - Purple (β=1.145): Highest initial ACF values, decaying slowly with lag. - Orange (β=0.662): Intermediate decay rate. - Green (β=0.362): Fastest decay, reaching near-zero ACF at lower lags. - **Spatial Grounding**: Legend positioned in the top-right corner of the plot. --- ### Key Observations 1. **Plot (a)**: - All β values show strong linear relationships (R² ≥ 0.977). - Higher β values correspond to steeper slopes, suggesting a direct relationship between β and the rate of change in log(p) with energy. - Green (β=0.362) has the flattest slope despite the highest R², indicating minimal energy dependence. 2. **Plot (b)**: - ACF decays with increasing lag for all β values, consistent with diminishing autocorrelation over time/intervals. - Higher β values (e.g., purple) retain stronger initial autocorrelation, while lower β values (green) exhibit rapid decay. --- ### Interpretation - **Plot (a)** suggests that the parameter β modulates the sensitivity of log(p) to energy changes. Higher β values (e.g., 1.145) imply a stronger dependency, while lower β values (e.g., 0.362) indicate near-independence. The near-perfect R² values suggest the linear models are highly reliable. - **Plot (b)** reveals that β also influences the autocorrelation structure. Higher β values maintain stronger correlations over longer lags, potentially indicating persistent dependencies in the system. The decay patterns align with typical time-series behavior, where correlations weaken as temporal or spatial separation increases. - **Cross-Plot Insight**: The consistent β values across both plots imply a unified parameter governing both energy-response dynamics (Plot a) and temporal autocorrelation (Plot b). This could reflect a physical or statistical property (e.g., damping, memory) in the underlying system. </details> Figure 5: Boltzmann sampling behavior. (a) Plot of $\log$ (p) (p: probability) versus energy for varying noise intensities ( $K_{n}=0.05$ , 0.1, 0.15), corresponding to different effective temperatures. The extracted effective inverse temperatures ( $\beta$ ) and the quality of the linear fits ( $R^{2}$ ) are shown in the inset. (b) Autocorrelation function (ACF) of the system energy as a function of lag, exhibiting a decay toward zero, confirming that the stochastic dynamics effectively decorrelate successive configurations. While the above simulations demonstrate the application of the OIM’s stochastic sampling capability in solving COPs such as MaxCut, we also employ the same example to further corroborate their Boltzmann sampling behavior. Specifically, we perform simulations on the 15-node graph considered in Fig. 4 under varying noise intensities, which effectively correspond to different effective temperatures. Similar to the five-state adder analysis, we simulate only the dynamics of the oscillator being sampled, while the remaining oscillator phases are held fixed. Figure 5 (a) shows the resulting $\log$ (p) (p: probability) versus energy distributions obtained over 5000 sweeps. The observed linear dependence confirms that the sampled state probabilities follow the expected Boltzmann relation, with the slope yielding the effective inverse temperature ( $\beta$ ). Furthermore, Fig. 5 (b) presents the normalized autocorrelation function (ACF) [33] of the system energy as a function of the lag. In all cases, the autocorrelation decays toward zero with increasing delay, indicating that the stochastic dynamics efficiently decorrelate successive configurations. <details> <summary>Figure4r.png Details</summary>  ### Visual Description ## Line Graph with Step Function: Bifurcation Analysis ### Overview The image contains two panels: - **(a)** A complex multi-line graph showing fluctuating values of Φ(π) over time, with sharp transitions and plateaus. - **(b)** A step function graph depicting a binary "Cut" state over time, with a single drop and recovery. ### Components/Axes #### Panel (a): - **Y-axis**: Φ(π) (unitless, range: 0.0 to 1.0). - **X-axis**: Time (range: 0 to 300). - **Lines**: Multiple colored lines (red, blue, purple, orange, pink, teal) representing distinct data series. - **Legend**: Located on the right, but colors do not match the lines (e.g., red in legend does not correspond to red lines). - **Annotations**: - Red arrow pointing to a sharp drop in Φ(π) around time 200. - Text "Before Bifurcation" near the bottom left. #### Panel (b): - **Y-axis**: Cut (unitless, range: 35 to 40). - **X-axis**: Time (range: 0 to 300). - **Line**: Single green step function. - **Annotations**: - Red arrow pointing to a drop in "Cut" at time 10. - Text "DIM" near the end of the graph. ### Detailed Analysis #### Panel (a): - **Trends**: - All lines exhibit sharp transitions between plateaus (e.g., red line drops from ~0.8 to ~0.4 at time 200). - Lines overlap significantly, with no clear separation between data series. - Values oscillate between 0.0 and 1.0, with most lines clustering near 0.5–0.8. - **Key Data Points**: - Sharp drop in Φ(π) at time 200 (red arrow). - Multiple plateaus (e.g., ~0.6–0.8 between 50–150, ~0.4–0.6 between 150–250). #### Panel (b): - **Trends**: - "Cut" starts at 40, drops to 35 at time 10, remains constant until time 250, then jumps back to 40. - **Key Data Points**: - Drop at time 10 (red arrow). - Recovery at time 250. ### Key Observations 1. **Panel (a)**: - Multiple data series show synchronized fluctuations, suggesting correlated events. - The red arrow highlights a critical drop in Φ(π) at time 200, potentially indicating a bifurcation event. - Legend colors do not match lines, raising concerns about data integrity. 2. **Panel (b)**: - Binary "Cut" state changes abruptly at time 10 and 250, possibly representing a threshold-based system. - "DIM" annotation at time 250 may indicate a secondary event or state change. ### Interpretation - **Panel (a)** likely represents a system with multiple interacting variables (e.g., physical, chemical, or biological processes) exhibiting chaotic or threshold-driven behavior. The sharp drop at time 200 could signify a bifurcation point where the system transitions to a new state. - **Panel (b)** simplifies the system into a binary "Cut" state, which aligns with the drop in Φ(π) in Panel (a) at time 10. The "DIM" label at time 250 might correlate with the recovery in Panel (a), suggesting a delayed response or secondary effect. - **Discrepancies**: The mismatched legend in Panel (a) implies potential errors in data labeling or visualization. Further validation of the legend is critical for accurate interpretation. - **Relationship**: Both panels share the same time axis, implying temporal coupling. The "Before Bifurcation" text in Panel (a) suggests the data captures pre-critical behavior, with Panel (b) possibly representing a simplified proxy for the same phenomenon. </details> Figure 6: Operating DIMs as p-bit platforms. Evolution of (a) phases in the DIM model; and (b) corresponding graph cut over time. The same graph shown in Fig. 4 has been considered in the simulation. With stochastic sampling, the system is able to compute the MaxCut (=39) ( $K=1;K_{s,\text{max}}=4;\>\>\text{FHI}^{0}=50;\>\>K_{n}\text{ for sampled oscillator}=0.006$ ). VI Realizing P-bit Engines Using Other Analog Ising Machines Beyond conventional OIMs, the proposed sampling methodology exhibits potential for generalization to a broader class of analog dynamical systems. As a case in point, we evaluate the implementation of the proposed sampling mode in the Dynamical Ising Machine (DIM) recently introduced by the authors [34]. Unlike the traditional Kuramoto model, which relies on phase differences, the DIM employs additive phase interactions. The DIM dynamics can be described by: $$ \begin{split}\frac{d\phi_{i}}{dt}&=-K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\phi_{i}\bm{+}\phi_{j})-K_{s}\sin(2\phi_{i})\end{split} \tag{17} $$ Specifically, when $K_{s}$ is below a certain threshold, the system stabilizes at the trivial state $\phi=\frac{\pi}{2}$ . As $K_{s}$ increases beyond this threshold, the system undergoes a bifurcation, leading to the emergence of stable phase configurations at $\phi∈\{0,\pi\}^{N}$ , which can subsequently be mapped to a spin configuration. The phase dynamics of the DIM, without stochastic sampling, are presented in Appendix I for the graph considered in Fig. 4 . While a detailed analysis of the DIM dynamics has been presented in [34], it is worth noting that the system exhibits a pitchfork bifurcation, qualitatively similar to that observed in the other popular Ising machine models such as the simulated bifurcation machine (SBM) [35]. This similarity suggests the feasibility of performing stochastic sampling in a broad class of analog dynamical systems beyond the OIM. Similar to the OIM, the DIM dynamics in the rotated frame of reference are given by, $$ \begin{split}\frac{d\epsilon_{i}}{dt}&=K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\epsilon_{i}+\epsilon_{j})+K_{s}\sin(2\epsilon_{i})\\ \\ &=\cos(\epsilon_{i})\left(K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\epsilon_{j})\right)+K_{s}\sin(2\epsilon_{i})\\ \\ &\equiv-\gamma_{i}\cos(\epsilon_{i})\ +K_{s}\sin(2\epsilon_{i})\end{split} \tag{18} $$ where $\gamma_{i}$ has the same definition and meaning as that in the case of the OIM. Moreover, equation (18) is exactly the same as the corresponding equation (12) derived for the OIM. Therefore, by employing the same approach and constraints used for the OIM, the updated spin state using the DIM can be derived as follows: $$ \begin{split}s_{i}^{+}&\approx-\text{sgn}\left[\tanh(-\gamma_{i}\Delta t)+\epsilon^{\eta}(\Delta t).\text{sech}^{2}(\gamma\Delta t)\right]\\ \\ &\approx\text{sgn}\left[\tanh\left(K\Delta t\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}s_{j}\right)-\vartheta\right]\end{split} \tag{19} $$ Figure 6 presents a simulation demonstrating the operation of the DIM in sampling mode to compute the MaxCut of the same graph considered in Fig. 4. As shown in Fig. 6 (a), the phase dynamics initially exhibit a bifurcation as $K_{\text{s}}$ is ramped up from $K_{\text{s}}=0$ to $K_{\text{s}}=4$ (not shown in the figure). At this stage, however, the system has not yet reached the optimal solution. With the onset of stochastic sampling, the system begins to sample the solution space and eventually converges to the optimal MaxCut value of 39, as shown in Fig. 6 (b).We note that while both the models exhibit Gibbs sampling, the characteristics (e.g., effective temperature) obtained for a specific set of parameters ( $K,K_{\text{s}},K_{\text{n}}$ ) will depend on the choice of exact dynamics. VII Conclusion This work builds a conceptual bridge between two paradigms that have traditionally been regarded distinct: analog oscillator-based Ising machines (OIMs) and stochastic sampling-based p-bit engines. By leveraging the natural dynamics of coupled oscillators—specifically through the interplay of SHI and FHI signals—we demonstrate that analog OIMs can perform stochastic sampling without requiring explicit computation of energy functions or the synaptic feedback. An intrinsic feature of this approach is the initial phase evolution, during which the oscillator network naturally performs gradient descent that involves the oscillator phases evolving simultaneously. Starting from random initial conditions, the phases converge to discrete states ( $\phi_{i}≈ 0$ or $\phi_{i}≈\pi$ ), effectively settling into a local minimum of the Ising energy landscape. In certain applications such as combinatorial optimization, this simultaneous evolution has the potential to offer a potential speed-up, positioning the system in a low-energy configuration even before the onset of the sampling-mode operation. However, these promising features come with the trade-off of requiring physical connectivity among oscillators—digital p-bit designs are better suited to implement the interaction between p-bits. From an implementation standpoint, this makes the analog approach more appropriate for sparse architectures which aligns well with the operational regimes where traditional p-bit platforms are expected to perform well. Nevertheless, it is important to recognize several key factors that must be optimized to ensure reliable operation: (a) The SHI pulse characteristics (amplitude and slew rate) must be carefully designed, as they determine $\Delta t$ , which in turn governs the effective temperature, the effective noise, and the validity of the approximations used in Eq. (8) to achieve Gibbs sampling. In practice, a high slew-rate SHI signal (resulting in a small $\Delta t$ ) is generally desirable to simultaneously satisfy these requirements. (b) Beyond $\Delta t$ , the effective temperature also depends on the coupling strength, $K$ , among oscillators. Therefore, the coupling strength must be co-designed with $\Delta t$ to realize the desired temperature range, while maintaining frequency locking and satisfying the weak-coupling condition implicit in the dynamics considered here. Since $K$ itself is determined by factors such as the natural frequency, oscillation amplitude, and the magnitude of the physical quantity (e.g., coupling resistance) characterizing the coupling element, these parameters must also be carefully considered to ensure robust operation. (c) Finally, the present analysis does not account for parameter variations or natural frequency mismatches–effects that may play a critical role in determining whether robust operation can be achieved. These design aspects and challenges will need to be systematically investigated going forward. From an algorithmic standpoint, we note that since Gibbs sampling is inherently sequential, incorporating techniques such as the graph-coloring method proposed by Niazi et al. [36] which enables multiple spins to be updated simultaneously— need exploration to further enhance performance, and represents a compelling direction for future exploration. Looking ahead, the capability of OIMs—when configured as p-bit engines—to perform stochastic sampling opens promising avenues for neural network applications, including the training of restricted Boltzmann machines. The overlaps between the analog and probabilistic paradigms identified in this work motivate the extension of this framework to more generalized computational models, including higher-order Ising machines [37, 38], p-bits with more than two states [39, 22], as well as to other analog systems [35, 40, 41, 42, 43, 44] beyond those explored here. ACKNOWLEDGEMENTS We gratefully acknowledge Prof. Kerem Camsari for providing valuable insights on the sampling properties of p-bits. This material is based upon work supported in part by ARO award W911NF-24-1-0228 and National Science Foundation grants ( $\#$ 2422333, #2433871). APPENDIX A INJECTION LOCKING IN SINGLE OSCILLATOR We analyze the impact of FHI on oscillator dynamics from an energy-based perspective. For an oscillator driven by an FHI signal at its natural frequency but with a phase offset $\theta$ , the phase dynamics can be described using Adler’s equation [31, 32] as follows: $$ \frac{d\phi}{dt}=-K_{c}\sin(\phi-\theta) \tag{20} $$ The corresponding energy function that the above dynamics minimize can be expressed as: $$ E(\phi)=-K_{c}\cos(\phi-\theta) \tag{21} $$ For this system, the relationship $-\frac{dE}{d\phi}=\frac{d\phi}{dt}$ holds, implying that: $$ \frac{dE}{dt}=\frac{dE}{d\phi}\cdot\frac{d\phi}{dt}=-\left(\frac{d\phi}{dt}\right)^{2}\leq 0 $$ Thus, the energy monotonically decreases over time, and the system evolves toward a stable fixed point. The minimum of the energy function occurs at $\phi^{*}=\theta$ , where $\frac{dE}{dt}=0$ . Within the domain $\phi∈[\theta,\theta+2\pi)$ , the only other fixed point satisfying $\frac{dE}{dt}=\frac{d\phi}{dt}=0$ is at $\phi^{*}=\theta+\pi$ , which corresponds to a local maximum, and is therefore unstable. For all other values of $\phi$ in this domain, $\frac{dE}{dt}<0$ . Consequently, the oscillator phase converges to the stable fixed point $\phi=\theta$ , although perturbations may be necessary to prevent the dynamics from settling at the unstable fixed point $\phi^{*}=\theta+\pi$ . APPENDIX B STABILITY OF THE FIXED POINTS OF THE DYNAMICS We analyze the stability of the fixed points associated with the dynamics described in Eq. (5) of the main text. As previously discussed, the fixed points of the dynamics are given by, $\epsilon_{1}^{*}=±\frac{\pi}{2}$ and $\sin(\epsilon_{2}^{*})=\frac{\gamma_{i}}{2K_{s}}$ . To investigate the local stability of the system, we analyze the sign of the second derivative: $$ H(\epsilon_{i})=\frac{d^{2}\epsilon_{i}}{dt^{2}}=\gamma_{i}\sin(\epsilon_{i})+2K_{s}\cos(2\epsilon_{i}) $$ We seek conditions under which $H(\epsilon_{i})<0$ , indicating local stability. - Stability of $\epsilon_{1}^{*}=\frac{\pi}{2}$ At this point, $\sin(\epsilon_{1}^{*})=1$ , $\cos(2\epsilon_{1}^{*})=-1$ , so: $$ H(\epsilon_{1}^{*})=\gamma_{i}+(-2K_{s})=\gamma_{i}-2K_{s} $$ Hence, $H(\epsilon_{1}^{*})<0$ when $\gamma_{i}<2K_{s}$ . - Stability of $\epsilon_{1}^{*}=-\frac{\pi}{2}$ Here, $\sin(\epsilon_{1}^{*})=-1$ , $\cos(2\epsilon_{1}^{*})=-1$ , therefore: $$ H(\epsilon_{1}^{*})=-\gamma_{i}-2K_{s} $$ Thus, $H(\epsilon_{1}^{*})<0$ when $\gamma_{i}+2K_{s}>0$ . - Combined Condition: Both fixed points $\epsilon_{1}^{*}=±\frac{\pi}{2}$ are stable when: $$ \gamma_{i}^{2}<4K_{s}^{2} $$ - Stability of $\sin(\epsilon_{2}^{*})=\frac{\gamma_{i}}{2K_{s}}$ Substituting into $H(\epsilon_{i})$ , we find that: $$ H(\epsilon_{2}^{*})<0\quad\text{when}\quad\gamma_{i}^{2}>4K_{s}^{2} $$ APPENDIX C TUNING ENERGY BARRIER WITH SHI STRENGTH ( $K_{s}$ ) <details> <summary>Figure5r.png Details</summary>  ### Visual Description ## Line Graph: Energy vs. ε (π) with Varying K_s Values ### Overview The graph depicts three energy curves plotted against ε (π), with distinct peaks and troughs. Each curve corresponds to a different K_s value (0.10, 0.15, 0.20), represented by orange, green, and purple lines respectively. The y-axis measures energy in arbitrary units (a.u.), while the x-axis ranges from -0.5 to 0.5 in ε (π). A label "γ = 0" is positioned in the top-right corner. --- ### Components/Axes - **Y-Axis**: Labeled "Energy (a.u.)" with ticks at -0.1, 0.0, and 0.1. - **X-Axis**: Labeled "ε (π)" with ticks at -0.5, 0.0, and 0.5. - **Legend**: Located on the right side, associating: - Orange line → K_s = 0.10 - Green line → K_s = 0.15 - Purple line → K_s = 0.20 - **Additional Text**: "γ = 0" in the top-right corner. --- ### Detailed Analysis 1. **Curve Trends**: - **Purple (K_s = 0.20)**: Highest peak at ε = 0 (~0.08 a.u.), deepest trough at ε = ±0.5 (~-0.08 a.u.). - **Green (K_s = 0.15)**: Intermediate peak (~0.06 a.u.) and trough (~-0.06 a.u.). - **Orange (K_s = 0.10)**: Lowest peak (~0.04 a.u.) and trough (~-0.04 a.u.). - All curves are symmetric about ε = 0, forming parabolic-like shapes. 2. **Key Data Points**: - At ε = 0: - Purple: ~0.08 a.u. - Green: ~0.06 a.u. - Orange: ~0.04 a.u. - At ε = ±0.5: - Purple: ~-0.08 a.u. - Green: ~-0.06 a.u. - Orange: ~-0.04 a.u. --- ### Key Observations - **Symmetry**: All curves exhibit even symmetry (mirror-image behavior around ε = 0). - **K_s Dependency**: Higher K_s values correlate with increased energy magnitude (both peaks and troughs). - **Spacing**: Curves are closely spaced, suggesting a gradual relationship between K_s and energy. --- ### Interpretation The graph likely models a quantum mechanical or vibrational system where K_s represents a coupling constant or interaction strength. The symmetry implies the system has even parity (e.g., a harmonic oscillator or particle in a symmetric potential). The γ = 0 label suggests no damping or external perturbation is considered. The energy levels' dependence on K_s indicates that stronger interactions (higher K_s) result in more pronounced energy states, which could relate to binding energy, resonance frequencies, or stability thresholds in the system. </details> Figure 7: Tuning energy barrier with $\mathbf{K_{s}}$ . Evolution of the energy barrier with $K_{s}$ We analyze the impact of $K_{s}$ on the energy barrier. The function corresponding to the dynamics described by Eq. (5) is, $$ E=\gamma\sin(\epsilon)+\frac{1}{2}K_{s}\cos(2\epsilon) \tag{22} $$ Figure 7 shows the resulting energy landscape for different values of $K_{s}$ ( $\gamma=0$ ). The energy difference between the peak energy (at $\epsilon=0$ ) and either valley (at $\epsilon=±\frac{\pi}{2}$ ) is given by $(\Delta E)_{max}=|K_{s}|$ . APPENDIX D TEMPORAL EVOLUTION OF THE PHASE We now analyze the dynamics presented in Eq. (5). The implicit analytical solution is given by, $$ \frac{\zeta(\epsilon)-\gamma\tanh^{-1}(\sin(\epsilon))}{\gamma^{2}-4K_{s}^{2}}=t+C \tag{23} $$ where, $$ \begin{split}\zeta(\epsilon)&=K_{s}\big(-2\log(\gamma-2K_{s}\sin(\epsilon))\\ &+\log(1-\sin(\epsilon))+\log(\sin(\epsilon)+1)\big)\\ \end{split} \tag{24} $$ and $C$ is the constant of integration. $C=0$ since $\epsilon(t=0)=0$ , and $K_{s}=0$ at $t=0$ . We now analyze the dynamics under the constraints specified in the main text, which are also restated below for reference: (i) The dynamics are evaluated in the limit $t→ 0$ . (ii) We model the noise as Gaussian white noise with $\langle\eta(t)\rangle=0,\quad\langle\eta(t)\eta(t^{\prime})\rangle=2K_{n}\delta(t-t^{\prime}),$ with $K_{n}$ denoting the noise intensity i.e., $\eta(t)\,dt=\sqrt{2K_{n}}\,dW_{t},$ where $W_{t}$ is a Wiener process. Equivalently, over a finite timestep $\Delta t$ , $∈t_{t}^{t+\Delta t}\eta(\tau)\,d\tau\;\sim\;\mathcal{N}\!\big(0,\,2K_{n}\,\Delta t\big).$ (iii) At the onset of sampling ( $t→ 0$ ) , we assume $K_{s}→ 0$ such that $K_{s}\ll|\gamma|$ . This reflects the requirement for a low energy barrier in the probabilistic regime. Eq. (23) can be rearranged to yield, $$ \begin{split}\sin(\epsilon_{i})&=\tanh\left(-\gamma t+\frac{4K_{s}^{2}t+\zeta(\epsilon)}{\gamma}+\epsilon^{\eta}(t)\right)\\ \\ &=\tanh\left(-\gamma t+\frac{4K_{s}^{2}t+\zeta(\epsilon)}{\gamma}+\epsilon^{\eta}(t)\right)\end{split}\ \tag{25} $$ Here, $\epsilon^{\eta}(t)$ represents the perturbation induced by noise. We note that in the above analysis, we approximate the impact of noise as phase jitter; a more exhaustive treatment would model the phase dynamics as a stochastic differential equation. From the above equation, we group all the terms dependent on $K_{s}$ as, $$ \begin{split}&\Theta=\frac{4K_{s}^{2}t+\zeta(\epsilon)}{\gamma}\\ \\ \Rightarrow\quad&\sin(\epsilon_{i})=\tanh\left(-\gamma t+\Theta+\epsilon^{\eta}(t)\right)\end{split} $$ and evaluate $\Theta$ under the constraints stated above. We begin by simplifying the following terms: $\quad\log(1-\sin(\epsilon))+\log(1+\sin(\epsilon))=\log(\cos^{2}(\epsilon))$ . Thus, $\zeta(\epsilon)=2K_{s}\left(\log\left(\cos(\epsilon)\right)-\log\left(\gamma-2K_{s}\sin(\epsilon)\right)\right)$ Next, we apply the following approximations (applicable under the constraints stated above), $\bullet$ $\quad\cos(\epsilon)≈ 1-\frac{\epsilon^{2}}{2}\Rightarrow\log(\cos(\epsilon))≈-\frac{\epsilon^{2}}{2}$ $\bullet$ $\quad\sin(\epsilon)≈\epsilon$ $\bullet$ $\quad\log(\gamma-2K_{s}\sin(\epsilon))≈\log(|\gamma|)-\frac{2K_{s}\epsilon}{\gamma}$ Substituting these approximations into expression for $\zeta$ yields, $$ \zeta\approx 2K_{s}\left(\frac{-\epsilon^{2}}{2}-\log(|\gamma|)+\frac{2K_{s}\epsilon}{\gamma}\right) $$ Subsequently, $\Theta$ can be approximated as, $$ \begin{split}\Theta&=\frac{4K_{s}^{2}t+2K_{s}\left(\frac{-\epsilon^{2}}{2}-\log(|\gamma|)+\frac{2K_{s}\epsilon}{\gamma}\right)}{\gamma}\\ \\ &=\frac{4K_{s}^{2}t}{\gamma}-\frac{K_{s}\epsilon^{2}}{\gamma}-\frac{2K_{s}\log(|\gamma|)}{\gamma}+\frac{4K_{s}^{2}\epsilon}{\gamma^{2}}\end{split} $$ The above expression shows the leading-order behavior of terms dependent on $K_{s}$ in the phase behavior. Moreover, when $K_{s}→ 0$ such that $K_{s}\ll|\gamma|$ , $\Theta→ 0$ . Nevertheless it is important that the system parameters be carefully designed to ensure that the dynamics emulate Boltzmann sampling as close as possible. <details> <summary>Figure6r.png Details</summary>  ### Visual Description ## Line Graphs: ξ vs Time for Different K_s Values ### Overview Three panels (a, b, c) display the evolution of ξ over time (0–6) for distinct stiffness constants (K_s = 0, 2, 5). Each panel shows multiple colored lines representing different initial conditions or parameters, converging toward ξ = 0 (dashed gray line). ### Components/Axes - **X-axis**: Time (0–6), labeled "Time" in black. - **Y-axis**: ξ (ranging from -0.5 to 0.5), labeled "ξ" in black. - **Legend**: Located on the right side of each panel, mapping colors to K_s values (e.g., red = K_s=0, blue = K_s=2, etc.). - **Dashed Line**: Horizontal line at ξ = 0 across all panels. - **Annotations**: - Panel (b): "ε_initial = -0.05π" (blue arrow pointing to orange line). - Panel (a): "K_s=0" (red text). - Panel (b): "K_s=2" (red text). - Panel (c): "K_s=5" (red text). ### Detailed Analysis #### Panel (a): K_s=0 - **Lines**: 10+ colored lines (red, blue, green, purple, etc.) curve upward, approaching ξ=0 asymptotically. - **Trends**: All lines increase monotonically but never cross ξ=0. - **Key Data**: - Red line (bottom): Starts at ξ ≈ -0.5, reaches ξ ≈ -0.1 by t=6. - Blue line: Starts at ξ ≈ -0.4, reaches ξ ≈ -0.05 by t=6. - Green line: Starts at ξ ≈ -0.3, reaches ξ ≈ 0 by t=6. #### Panel (b): K_s=2 - **Lines**: 10+ colored lines start at ε_initial = -0.05π (≈ -0.157) and rise toward ξ=0. - **Trends**: Lines flatten faster than in (a), with sharper initial slopes. - **Key Data**: - Orange line (highlighted): Starts at ξ ≈ -0.157, reaches ξ ≈ -0.02 by t=2, then plateaus. - Blue line: Starts at ξ ≈ -0.157, reaches ξ ≈ 0 by t=4. - Green line: Starts at ξ ≈ -0.157, overshoots ξ=0 slightly before stabilizing. #### Panel (c): K_s=5 - **Lines**: 10+ colored lines remain nearly flat near ξ=0. - **Trends**: Minimal deviation from ξ=0 across all time. - **Key Data**: - All lines cluster tightly around ξ=0 (within ±0.01). - No significant movement observed. ### Key Observations 1. **Convergence Behavior**: - Higher K_s values (5 > 2 > 0) result in faster convergence to ξ=0. - At K_s=5, the system stabilizes immediately; at K_s=0, convergence is slowest. 2. **Initial Condition Impact**: - Panel (b) shows a deliberate perturbation (ε_initial = -0.05π), causing transient oscillations before stabilization. 3. **Color Consistency**: - Legend colors match line colors across panels (e.g., red = K_s=0 in all panels). ### Interpretation The graphs model a damped system approaching equilibrium (ξ=0) with stiffness-dependent dynamics: - **K_s=0**: No damping; system evolves slowly toward equilibrium. - **K_s=2**: Moderate damping; initial perturbations decay with oscillations. - **K_s=5**: Overdamped; system stabilizes instantly without overshoot. The ε_initial perturbation in (b) demonstrates how initial displacements are mitigated by increasing K_s. The dashed ξ=0 line serves as a universal reference for equilibrium across all conditions. </details> Figure 8: Role of SHI in maintaining phase trajectory. Phase deviation $\epsilon$ as a function of time for: (a) $K_{s}=0$ ; (b) $K_{s}=2$ ; (c) $K_{s}=5$ . The results in the illustrative example show that a critical SHI is needed to help the oscillator phase continue to evolve along the direction of its initial perturbation. $|\gamma_{1}|=1$ has been considered in this example. Under these conditions, Eq. (LABEL:appendix4-2) can be approximated as, $$ \begin{split}\sin(\epsilon_{i})&\approx\tanh\left(-\gamma t+\epsilon^{\eta}(t)\right)\\ \\ &\approx\tanh\left(-\gamma t\right)+\epsilon^{\eta}(t).\text{sech}^{2}(\gamma t)\\ \\ \end{split} $$ Here, $\epsilon^{\eta}(t)\sim\mathcal{N}\!\big(0,\,2K_{n}\,t\big)$ , and is small such that the $\tanh(.)$ term can be linearized. The noise term, $\epsilon^{\eta}(t).\text{sech}^{2}(\gamma t)$ , can be considered as a gaussian distribution representing white noise, and scaled by a function of the synaptic input—sech ${}^{2}(\gamma t)$ . The updated spin state, $s^{+}=-\text{sgn}\left(\sin(\epsilon^{+})\right)$ , at a short time instant $\Delta t→ 0$ after sampling has been initiated (at $t=0$ ), can be expressed as, $$ \begin{split}s^{+}&\approx\text{sgn}\left[\tanh(\gamma\Delta t)-\epsilon^{\eta}(\Delta t).\text{sech}^{2}(\gamma\Delta t)\right]\\ \\ &\equiv\text{sgn}\left[\tanh(\gamma\Delta t)-\vartheta\right]\end{split} $$ where, $\vartheta\equiv\epsilon^{\eta}(\Delta t).\text{sech}^{2}(\gamma\Delta t)$ . Since $\text{sech}^{2}(\gamma\Delta t)∈(0,1]$ , and noise power is assumed to be small, $\vartheta$ has a high probability of being in the interval [-1,+1]. We also note that for $\Delta t→ 0$ , $\text{sech}^{2}(\gamma\Delta t)→ 1\Rightarrow\vartheta≈eq\epsilon^{\eta}(\Delta t)$ . Moreover, the gaussian distribution is more representative of the thermal noise found in physical devices. We also consider the case when both $\gamma$ and $K_{s}$ are small. Under this assumption, the dynamics for $\epsilon\ll 1$ can be approximated as, $$ \frac{d\epsilon}{dt}\approx-\gamma+2K_{s}\epsilon $$ Unlike the previous case, using the SDE framework here yields a tractable and elegant solution that offers a clear and intuitive picture of the relative competition between the stochastic and deterministic components. We begin by considering the linear Itô SDE $$ \begin{split}&d\epsilon(t)=\big(-\gamma+2K_{s}\epsilon(t)\big)\,dt+\sqrt{2K_{n}}\,dW_{t},\\ \\ &\text{where we note that }\epsilon(0)=0,\;K_{s}>0.\end{split} \tag{0} $$ Integrating factor: Let $M(t)=e^{-2K_{s}t}$ . Since $M$ is deterministic, $$ \displaystyle d\!\big(M(t)\epsilon(t)\big) \displaystyle=M(t)\,d\epsilon(t)+\epsilon(t)\,dM(t) \displaystyle=-\gamma M(t)\,dt+\sqrt{2K_{n}}\,M(t)\,dW_{t}. \tag{27} $$ Integration from $0$ to $t$ yields $$ \displaystyle M(t)\epsilon(t) \displaystyle=-\gamma\int_{0}^{t}M(s)\,ds \displaystyle\quad+\sqrt{2K_{n}}\int_{0}^{t}M(s)\,dW_{s}. \tag{28} $$ Since $$ \int_{0}^{t}e^{-2K_{s}s}\,ds=\frac{1-e^{-2K_{s}t}}{2K_{s}},\quad M(t)^{-1}=e^{2K_{s}t}, $$ we obtain $$ \epsilon(t)=\frac{\gamma}{2K_{s}}\,(1-e^{2K_{s}t})+\sqrt{2K_{n}}\int_{0}^{t}e^{2K_{s}(t-s)}\,dW_{s}. \tag{29} $$ Mean: The stochastic integral in (29) has zero mean: $$ \mathbb{E}[\epsilon(t)]=\frac{\gamma}{2K_{s}}\,(1-e^{2K_{s}t}). \tag{30} $$ Variance: By Itô isometry, $$ \displaystyle\mathrm{Var}[\epsilon(t)] \displaystyle=2K_{n}\int_{0}^{t}e^{4K_{s}(t-s)}\,ds \displaystyle=\frac{K_{n}}{2K_{s}}\,\big(e^{4K_{s}t}-1\big). \tag{31} $$ The above results imply that for $K_{s}>0$ , both mean and variance diverge exponentially. Thus, no stationary distribution exists. Two limiting cases further clarify the dynamics: - Case $\gamma=0$ . In the absence of the synaptic input, the deterministic contribution reduces to $d\epsilon=2K_{s}\epsilon\,dt$ . With $\epsilon(0)=0$ , this term alone would maintain $\epsilon(t)\equiv 0$ . However, in the presence of Gaussian noise, random perturbations are continuously injected and then exponentially amplified by the unstable drift. Consequently, the growth of $\epsilon(t)$ is seeded by stochastic fluctuations and deterministically amplified over time. - Case $\gamma≠ 0$ . When $\gamma≠ 0$ , both deterministic and stochastic contributions govern the dynamics. As seen from Eq. (30) and Eq. (31), both the drift-induced component and the noise-driven fluctuations diverge as $t$ increases. The relative dominance of these two effects depends on the balance between the deterministic drive, set by $|\gamma|$ , and the stochastic forcing, quantified by $K_{n}$ . A stronger deterministic bias leads to more predictable exponential growth, whereas larger noise intensity results in greater dispersion across trajectories. APPENDIX E ROLE OF SHI DURING STOCHASTIC SAMPLING To analyze how the SHI can help the oscillator maintain the phase trajectory resulting from the stochastic sampling process, we examine Eq. (5): $$ \frac{d\epsilon}{dt}=-\gamma\cos(\epsilon)+K_{s}\sin(2\epsilon) $$ As discussed in the main text, the role of SHI is particularly critical when the initial direction of phase relaxation is counter to that expected from synaptic feedback, owing to noise. As noted earlier, the synaptic input (alone) tends to push the phase toward the fixed point opposite to the sign of $\gamma$ . In other words, we seek to analyze the conditions under which the initial flow of the dynamics is towards $\epsilon=+\frac{\pi}{2}(-\frac{\pi}{2})$ , even when though the synaptic input is $\gamma>0\>\>\left(\gamma<0\right)$ , respectively. Assuming that the initial perturbation results in a phase magnitude given by $|\epsilon_{th}|$ , the critical condition on $K_{s}$ to ensure that the phase continues to flow in the same direction can be expressed as, $$ 2K_{s}\sin(|\epsilon_{th}|)>\left|\gamma\right| \tag{32} $$ This condition ensures the RHS of Eq. (5) maintains the same sign as the phase at $\epsilon=\epsilon_{th}$ . Further, since $\sin(.)$ is monotonic in the region, $|\epsilon|∈\{0,\frac{\pi}{2}\}$ , the sign of the RHS terms will not change until dynamics reach the corresponding fixed point. We note that in the presence of noise, the above condition should be interpreted in a probabilistic sense. We also illustrate this behavior using a simple example: a negatively coupled two-oscillator system (with $K=1$ ) in which the phase of oscillator 2 is fixed at $\epsilon_{2}=-\frac{\pi}{2}$ yielding $\gamma_{1}=-1\Rightarrow|\gamma_{1}|=1$ . In this configuration, the energetically favorable state for oscillator 1 is $\epsilon_{1}=\frac{\pi}{2}$ , and synaptic feedback is expected to drive the system toward this fixed point. To explore the system’s dynamics, oscillator 1 is initialized at various discrete phase values within the range $\epsilon_{\text{th}}∈[-0.45\pi,\,0.45\pi]$ . Figures 8 (a–c) show the evolution of $\epsilon$ for different values of $K_{s}$ . In the absence of synaptic hysteresis (i.e., $K_{s}=0$ ), the phase does not preserve the direction of its initial perturbation. Instead, it eventually aligns with the direction of the synaptic input, converging to $\epsilon=\frac{\pi}{2}$ . This behavior is expected, as the inequality in Eq. (32) is never satisfied in this regime. Introducing SHI (i.e., $K_{s}>0$ ) enables the oscillator phase to continue evolving along the direction of $\epsilon_{th}$ . However, the magnitude of $K_{s}$ must exceed a certain threshold. This behavior is illustrated in Figs. 8 (b) and 8 (c). In Fig. 8 (b), when $K_{s}$ is below the threshold required for that $\epsilon_{th}=-0.05\pi$ , the phase evolution eventually aligns with the direction favored by the synaptic input. In contrast, when $K_{s}$ is sufficiently, as in Fig. 8 (c)), all initial phase values considered here evolve along the direction of their original perturbation. Thus, the SHI injection strength, $K_{s}$ , must be carefully engineered to ensure that one phase magnitudes beyond a certain threshold evolve continue to evolve in the direction of their initial perturbation. APPENDIX F TIME INCREMENT $\Delta t$ AND IMPACT OF SLEW RATE OF THE SHI SIGNAL As described in the main text, the infinitesimal time $\Delta t$ is defined as the interval required for the oscillator phase to reach the critical threshold $|\epsilon_{\text{th}}|$ , beyond which the phase trajectory cannot reverse and the sampling event is effectively complete. The value of $\Delta t$ is determined by properties of the SHI signal, most notably its slew rate. Qualitatively, this dependence can be understood as follows: a higher slew rate drives the oscillator phase more rapidly toward the fixed points $\epsilon∈\{-\tfrac{\pi}{2},\tfrac{\pi}{2}\}$ , thereby reducing $\Delta t$ . To provide a quantitative illustration, we consider the SHI signal implemented as a linear ramp with slew rate $\rho$ , such that $K_{s}=\rho t$ . For this case, the phase dynamics (Eq. (5)), under the assumption $|\epsilon|\ll 1$ , can be expressed as: $$ \begin{split}&\frac{d\epsilon}{dt}=-\gamma+(\rho t)(2\epsilon)\\ \\ \Rightarrow\>&\frac{d\epsilon}{dt}-2\rho\,\epsilon\,t=-\gamma\end{split} \tag{33} $$ where, $\sin(2\epsilon)\sim 2\epsilon$ since $|\epsilon|\ll 1$ . Using the integrating factor, $\mu(t)=e^{-\rho t^{2}}$ , Eq. (33) can be expressed as, $$ \frac{d}{dt}\left(\epsilon e^{-\rho t^{2}}\right)=-\gamma e^{-\rho t^{2}} \tag{34} $$ With the initial condition $\epsilon(0)=0$ , the resulting phase evolution is described by, $$ \epsilon(t)=-\frac{\sqrt{\pi}\gamma e^{\rho t^{2}}\,\mathrm{erf}(\sqrt{\rho}t)}{2\sqrt{\rho}} \tag{35} $$ Using Eq. 35, we compute the time, $\Delta t$ taken by the oscillator phase to reach $|\epsilon_{\text{th}}|≈\frac{\gamma}{2K_{s,\text{th}}}$ $(|\epsilon_{th}|\ll 1),\;\text{where}\;K_{s,\text{th}}=\rho(\Delta t)$ . Equating Eq. (35) to $|\epsilon_{th}|$ yields, $$ \begin{split}\frac{\gamma}{2\rho(\Delta t)}=\frac{\sqrt{\pi}\gamma e^{\rho(\Delta t)^{2}}\,\mathrm{erf}(\sqrt{\rho}\Delta t)}{2\sqrt{\rho}}\\ \\ \Rightarrow\sqrt{\rho}\,\Delta t\,\pi\,e^{\rho\Delta t^{2}}\,\mathrm{erf}(\sqrt{\rho}\Delta t)=1\end{split} \tag{36} $$ Note that the relationship between $\rho$ and $\Delta t$ does not depend on $\gamma$ . To express this relationship in terms of $\rho$ , we define $y=\sqrt{\rho}\,\Delta t\Rightarrow\rho=\frac{y^{2}}{(\Delta t)^{2}}$ where $y$ obeys the relationship, $$ \sqrt{\pi}\,y\,e^{y^{2}}\,\mathrm{erf}(y)=1 \tag{37} $$ which is Eq. (36) expressed in terms of $y$ . Equation (37) implies that $y$ is a constant, whose value can be approximated as $y≈ 0.62$ . Thus, the relationship between the ramp rate of the SHI signal and $\Delta t$ , can be approximated as, $$ \begin{split}\rho\approx\frac{0.3844}{(\Delta t)^{2}}\\ \\ \equiv\Delta t\approx\frac{0.62}{\sqrt{\rho}}\end{split} \tag{38} $$ Equation (38) shows that $\Delta t$ is inversely proportional to the square root of the slew rate, under the approximations stated above. While this relationship was derived assuming a linear ramp signal for analytical convenience, we expect the same qualitative dependence of $\Delta t$ on $\rho$ to hold for other SHI signal schemes as well. <details> <summary>FigureKLD.png Details</summary>  ### Visual Description ## Scatter Plot: KL Divergence vs. τ for 5-state Adder ### Overview The image is a scatter plot showing the relationship between **KL Divergence** (y-axis) and **τ** (x-axis) for a system labeled "5-state Adder." Three green diamond markers are plotted at specific (τ, KL Divergence) coordinates. ### Components/Axes - **X-axis (τ)**: Labeled "τ" with values ranging from 0.01 to 0.05 in increments of 0.01. - **Y-axis (KL Divergence)**: Labeled "KL Divergence" with values from 0.00 to 0.15 in increments of 0.05. - **Data Points**: Three green diamond markers at: - (τ = 0.01, KL Divergence = 0.00) - (τ = 0.03, KL Divergence = 0.15) - (τ = 0.05, KL Divergence = 0.15) - **Text Label**: "5-state Adder" is positioned in the bottom-right corner of the plot. ### Detailed Analysis - **Data Points**: - At τ = 0.01, KL Divergence is 0.00 (minimum value). - At τ = 0.03 and τ = 0.05, KL Divergence reaches its maximum value of 0.15. - **Trend**: KL Divergence increases sharply from τ = 0.01 to τ = 0.03, then remains constant at τ = 0.05. ### Key Observations - The KL Divergence is zero at τ = 0.01, indicating no divergence at this parameter value. - A significant increase in KL Divergence occurs between τ = 0.01 and τ = 0.03, suggesting a critical threshold or transition in the system. - The divergence plateaus at τ = 0.05, implying stability or saturation beyond this parameter value. ### Interpretation The plot demonstrates that the KL Divergence of the 5-state Adder is highly sensitive to τ in the range 0.01–0.03, with a sharp increase. Beyond τ = 0.03, the divergence stabilizes, indicating that the system’s behavior becomes less variable or more predictable. The "5-state Adder" label suggests this plot is part of a study on a computational or algorithmic system with five states, where τ likely represents a parameter (e.g., time constant, threshold, or input value) influencing the system’s entropy or information divergence. The absence of additional data points or error bars limits conclusions about variability or confidence intervals. </details> Figure 9: Impact of SHI rise time and $\Delta t$ . Measured KL Divergence as a function of $\tau$ , where $K_{s}(t)=K_{s,\text{max}}(1-e^{-\frac{t}{\tau}})$ . $\tau$ impacts the time scale over which the SHI signal is asserted. The practical implications of the relationship between $\Delta t$ and $\rho$ are: (a) $\Delta t$ directly impacts the effective inverse temperature ( $\beta$ ) of the oscillator-based BSN and sets the requirements on the required slew rate. (b) Increasing $\Delta t$ suppresses the effect of noise, thereby making the system “less” stochastic. To elucidate this, we consider the expression for the updated spin state given by: $$ s^{+}=\mathrm{sgn}\left[\tanh(\gamma\Delta t)-\epsilon^{\eta}(\Delta t)\,\mathrm{sech}^{2}(\gamma\Delta t)\right] $$ The noise term, $\epsilon^{\eta}(\Delta t)\,\mathrm{sech}^{2}(\gamma\Delta t)$ , indicates that the noise perturbation is scaled by $\text{sech}^{2}(\gamma\Delta t)$ which has a maximum value (=1) when $\gamma\Delta t=0$ , and diminishes (with the value approaching 0) as the value of the argument of the $\text{sech}^{2}(.)$ increases. Consequently, a large value of $\Delta t$ (for a given $\gamma$ ), effectively diminishes the impact of noise making the system less stochastic. We evaluate the impact of the time-scale over which the SHI signal is asserted (and consequently, $\Delta t$ ) on the sampling properties. Fig. 9 shows the measured KL divergence parameter for a 5-state adder (considered above) as a function of $\tau$ – the time-constant associated with the SHI signal; $K_{s}(t)=K_{s,\text{max}}(1-e^{-\frac{t}{\tau}})$ . As noted above, this also impacts $\Delta t$ . Consistent with the preceding analysis, an increase in $\tau$ leads to a corresponding rise in the KL divergence, indicating degraded sampling quality at longer $\tau$ durations. APPENDIX G PHASE CONFIGURATION OF NON-SAMPLED OSCILLATORS As noted in the main text, it is important to ensure that the phases of the non-sampled oscillators are maintained at $$ \epsilon_{j}\in\left\{-\tfrac{\pi}{2},\tfrac{\pi}{2}\right\}\;\;\equiv\;\;\phi_{j}\in\{0,\pi\},\quad\forall\,j\in\{1,2,\dots,N\}\setminus\{i\}, $$ (where, $i$ is the sampled oscillator) to ensure that the oscillator dynamics map to Gibbs sampling. To explain this requirement, we being with Eq. (11), which is rewritten below for reference: $$ \frac{d\epsilon_{i}}{dt}=-K\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\epsilon_{i}-\epsilon_{j})+K_{s}\sin(2\epsilon_{i}) \tag{39} $$ Expanding the $\sin(\epsilon_{i}-\epsilon_{j})$ term, Eq. (39) can be expressed as, $$ \begin{split}&\frac{d\epsilon_{i}}{dt}=\\ &-K\bigg(-\cos(\epsilon_{i})\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\sin(\epsilon_{j})+\sin(\epsilon_{i})\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\cos(\epsilon_{j})\bigg)\\ &+K_{s}\sin(2\epsilon_{i})\end{split} \tag{40} $$ When $\epsilon_{j}∈\{-\frac{\pi}{2},\frac{\pi}{2}\}\equiv\phi_{j}∈\{0,\pi\}$ , $$ \sin(\epsilon_{i})\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\cos(\epsilon_{j})=0 $$ reducing Eq. (40) to Eq. (12) in the main text, which in turn establishes that the OIM dynamics can perform Gibbs sampling. In contrast, when $\epsilon_{j}∉\{-\frac{\pi}{2},\frac{\pi}{2}\}\equiv\phi_{j}∉\{0,\pi\}$ , then $$ \sin(\epsilon_{i})\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{N}J_{ij}\cos(\epsilon_{j})\neq 0 $$ which introduces an additional component orthogonal to $\cos(\epsilon_{i})$ that breaks the direct Gibbs mapping (except in special cases where $\sum_{\begin{subarray}{c}j=1\\ j≠ i\end{subarray}}^{N}J_{ij}\cos(\epsilon_{j})$ vanishes or averages to zero). APPENDIX H ADDITIONAL MAXCUT RESULTS We present the MaxCut results for ten additional randomly generated 15-node graphs. Figure 10 shows the distribution of the obtained cut values using a box plot, where each cut value is normalized with respect to the corresponding Maximum Cut. Each graph is simulated 10 times. <details> <summary>FigureMC_distribution.png Details</summary>  ### Visual Description ## Line Chart: Normalized Cut Values Across Graphs ### Overview The chart displays normalized cut values for 10 graphs (Graph #1 to Graph #10) with 15 nodes each, based on 10 trials per graph. The y-axis represents the normalized cut metric (0.5–1.0), while the x-axis enumerates graph instances. Red horizontal lines at y=1.0 serve as reference points, and blue boxes with red internal lines highlight specific graph ranges. Red dots indicate outlier or trial-specific values. ### Components/Axes - **X-axis (Graph #)**: Labeled "Graph #" with integer ticks from 1 to 10. - **Y-axis (Normalized cut)**: Labeled "Normalized cut" with increments of 0.1 from 0.5 to 1.0. - **Legend**: Located in the bottom-right corner, stating: - "N = 15 nodes" (graph size) - "10 trials per graph" (data granularity) - **Key Elements**: - **Red horizontal lines**: Fixed at y=1.0 across all graphs. - **Blue boxes**: Enclose data points for Graphs #3, #4, #7, #8, and #9. - **Red dots**: Scattered below the red lines (e.g., Graphs #2, #5, #8). - **Black crosshair**: Annotated at Graph #7, y≈0.9. ### Detailed Analysis 1. **Graph #1**: No visible data points; implied baseline at y=1.0 (red line). 2. **Graph #2**: Red dot at y≈0.95, below the red line. 3. **Graph #3**: Blue box spans y≈0.95–1.0, with a red line at y=1.0. 4. **Graph #4**: Blue box spans y≈0.95–1.0, with a red line at y=1.0. 5. **Graph #5**: Red dot at y≈0.95, below the red line. 6. **Graph #6**: No visible data points; implied baseline at y=1.0. 7. **Graph #7**: Blue box spans y≈0.9–1.0, with a red line at y=1.0. Black crosshair marks y≈0.9. 8. **Graph #8**: Blue box spans y≈0.95–1.0, with a red line at y=1.0. Red dot at y≈0.92. 9. **Graph #9**: Blue box spans y≈0.95–1.0, with a red line at y=1.0. 10. **Graph #10**: No visible data points; implied baseline at y=1.0. ### Key Observations - **Consistency**: Most graphs (1, 3, 4, 6, 9, 10) align with the red line at y=1.0, suggesting optimal normalized cut values. - **Variability**: Graphs #2, #5, and #8 show red dots below y=1.0, indicating lower cut values in specific trials. - **Confidence Intervals**: Blue boxes (Graphs #3, #4, #7, #8, #9) likely represent trial-based variability, with narrower ranges for #3, #4, and #9. - **Outliers**: The black crosshair at Graph #7 (y≈0.9) highlights a significant deviation from the baseline. ### Interpretation The data suggests that normalized cut values are generally high (close to 1.0), implying effective graph partitioning for most instances. However, the presence of red dots and variable blue box widths indicates trial-to-trial inconsistency. The red line at y=1.0 acts as a performance benchmark, while the blue boxes may reflect confidence intervals or interquartile ranges from the 10 trials. The outlier at Graph #7 (y≈0.9) warrants further investigation, as it deviates notably from the trend. The uniformity of Graphs #3, #4, and #9 (narrow blue boxes) suggests stable performance across trials, whereas Graphs #7 and #8 show greater variability. This could imply differences in graph structure or trial conditions affecting cut quality. </details> Figure 10: Computing MaxCut using OIM-based p-bit engine. Box plot showing distribution of graph cuts obtained for 10 randomly generated 15-node graphs. Each graph is simulated 10 times. The SHI scheme is the same as that used in the main text. APPENDIX I DYNAMICS OF DYNAMICAL ISING MACHINE Figure 11 evaluates the graph considered in Fig. 4 using the analog dynamics of the DIM (without stochastic sampling). A bifurcation similar to that exhibited by other models such as SBM [35] is observed. <details> <summary>Figure7r.png Details</summary>  ### Visual Description ## Line Graph: φ(π) vs. Time with Divergent Paths ### Overview The image depicts a line graph showing the behavior of a function φ(π) over time. Multiple colored lines originate from a common point at (Time = 0.5, φ(π) = 0.5) and diverge into distinct paths. The graph includes a legend labeled "DIM" in red, though the colors of the lines are not explicitly explained in the legend. The x-axis (Time) ranges from 0 to 2, and the y-axis (φ(π)) ranges from 0.0 to 1.0. --- ### Components/Axes - **X-axis (Time)**: Labeled "Time" with values from 0 to 2. - **Y-axis (φ(π))**: Labeled "φ(π)" with values from 0.0 to 1.0. - **Legend**: Located in the bottom-right corner, labeled "DIM" in red. No explicit color-to-label mapping is provided for the lines. - **Lines**: Multiple colored lines (red, blue, green, purple, orange, etc.) originate from the point (0.5, 0.5) and diverge into different trajectories. --- ### Detailed Analysis 1. **Initial Convergence**: All lines start at the same point (Time = 0.5, φ(π) = 0.5), suggesting a shared initial condition or critical threshold. 2. **Divergence Paths**: - **Red Line**: Ascends sharply to φ(π) = 1.0 by Time = 1.0, then remains flat. - **Blue Line**: Descends sharply to φ(π) = 0.0 by Time = 1.0, then remains flat. - **Green Line**: Rises to φ(π) ≈ 0.8 by Time = 1.0, then plateaus. - **Purple Line**: Drops to φ(π) ≈ 0.2 by Time = 1.0, then plateaus. - **Orange Line**: Remains flat at φ(π) = 0.5 throughout the entire time range. 3. **Unlabeled Colors**: Additional lines (e.g., teal, brown, pink) exhibit intermediate trends, but their exact paths are not clearly distinguishable due to overlapping or low contrast. --- ### Key Observations - **Critical Point at Time = 0.5**: All lines originate from the same point, indicating a shared initial state or bifurcation event. - **Divergence Patterns**: Lines exhibit distinct behaviors (ascending, descending, flat) after Time = 0.5, suggesting varying responses to an underlying parameter or condition. - **Unlabeled Legend**: The legend "DIM" in red does not clarify the meaning of the colored lines, leaving their interpretation ambiguous. - **Flat Lines**: The orange line (and possibly others) remains constant, implying a steady-state or equilibrium condition. --- ### Interpretation The graph likely represents a dynamical system or model where φ(π) evolves over time under different scenarios labeled as "DIM." The convergence at Time = 0.5 suggests a critical transition or bifurcation point, after which the system branches into multiple outcomes. The flat lines (e.g., orange) may represent stable equilibria, while the ascending/descending lines indicate dynamic changes. The absence of a clear legend for the colored lines limits the ability to assign specific meanings to each path. However, the overall structure implies that "DIM" could be a parameter or model influencing the system's behavior, with different trajectories reflecting varying conditions or inputs. </details> Figure 11: Dynamical Ising Machine. Evolution of $\phi$ in the DIM model for the graph considered in Fig. 4 (K=1; $K_{s}(t)=0.4t$ ). References - [1] T. Wang, L. Wu, P. Nobel, and J. Roychowdhury. Solving combinatorial optimisation problems using oscillator based Ising machines. Natural Computing 20 (2), 287–306 (2021). - [2] N. Mohseni, P. L. McMahon, and T. Byrnes. Ising machines as hardware solvers of combinatorial optimization problems. Nature Reviews Physics 4 (6), 363–379 (2022). - [3] A. Lucas, Ising formulations of many NP problems, Frontiers in Physics 2, 5 (2014). - [4] R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli, A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba, T. Honjo, et al. Experimental investigation of performance differences between coherent Ising machines and a quantum annealer. Science Advances 5 (5), eaau0823 (2019). - [5] T. Honjo, T. Sonobe, K. Inaba, T. Inagaki, T. Ikuta, Y. Yamada, T. Kazama, K. Enbutsu, T. Umeki, R. Kasahara, et al. 100,000-spin coherent Ising machine. Science Advances 7 (40), eabh0952 (2021). - [6] A. Litvinenko, R. Khymyn, R. Ovcharov, and J. Åkerman. A 50-spin surface acoustic wave Ising machine. Communications Physics 8 (1), 1–11 (2025). - [7] A. Mallick, M. K. Bashar, D. S. Truesdell, B. H. Calhoun, and N. Shukla. Overcoming the accuracy vs. performance trade-off in oscillator ising machines. In 2021 IEEE International Electron Devices Meeting (IEDM), 40–2 (2021). - [8] W. Moy, I. Ahmed, P.-W. Chiu, J. Moy, S. S. Sapatnekar, and C. H. Kim. A 1,968-node coupled ring oscillator circuit for combinatorial optimization problem solving. Nature Electronics 5 (5), 310–317 (2022). - [9] M. K. Bashar, A. Mallick, D. S. Truesdell, B. H. Calhoun, S. Joshi, and N. Shukla. Experimental demonstration of a reconfigurable coupled oscillator platform to solve the max-cut problem. IEEE Journal on Exploratory Solid-State Computational Devices and Circuits 6 (2), 116–121 (2020). - [10] J. Vaidya, R. S. Kanthi, and N. Shukla. Creating electronic oscillator-based Ising machines without external injection locking. Scientific Reports 12 (1), 981 (2022). - [11] O. Maher, M. Jiménez, C. Delacour, N. Harnack, J. Núñez, M. J. Avedillo, B. Linares-Barranco, A. Todri-Sanial, G. Indiveri, and S. Karg. A CMOS-compatible oscillation-based VO 2 Ising machine solver. Nature Communications 15 (1), 3334 (2024). - [12] H. Cılasun, W. Moy, Z. Zeng, T. Islam, H. Lo, A. Vanasse, M. Tan, M. Anees, A. Kumar, S. S. Sapatnekar, et al. A coupled-oscillator-based Ising chip for combinatorial optimization. Nature Electronics 1–10 (2025). - [13] A. Litvinenko, R. Khymyn, V. H. González, R. Ovcharov, A. A. Awad, V. Tyberkevych, A. Slavin, and J. Åkerman. A spinwave Ising machine. Communications Physics 6 (1), 227 (2023). - [14] A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, et al. Coherent quantum annealing in a programmable 2,000 qubit Ising chain. Nature Physics 18 (11), 1324–1328 (2022). - [15] M. K. Bashar, Z. Lin, and N. Shukla. Stability of oscillator Ising machines: Not all solutions are created equal. J. Appl. Phys. 134 (14), 144901 (2023). - [16] Y. Cheng, M. K. Bashar, N. Shukla, and Z. Lin. A control theoretic analysis of oscillator Ising machines. Chaos: An Interdisciplinary Journal of Nonlinear Science 34 (7) (2024). - [17] A. Allibhoy, A. N. Montanari, F. Pasqualetti, and A. E. Motter. Global Optimization Through Heterogeneous Oscillator Ising Machines. arXiv preprint arXiv:2505.17027 (2025). - [18] K. Y. Camsari, B. M. Sutton, and S. Datta. P-bits for probabilistic spin logic. Applied Physics Reviews 6 (1) (2019). - [19] K. Y. Camsari, R. Faria, B. M. Sutton, and S. Datta. Stochastic p-bits for invertible logic. Physical Review X 7 (3), 031014 (2017). - [20] N. A. Aadit, A. Grimaldi, M. Carpentieri, L. Theogarajan, J. M. Martinis, G. Finocchio, and K. Y. Camsari. Massively parallel probabilistic computing with sparse Ising machines. Nature Electronics 5 (7), 460–468 (2022). - [21] S. Chowdhury, A. Grimaldi, N. A. Aadit, S. Niazi, M. Mohseni, S. Kanai, H. Ohno, S. Fukami, L. Theogarajan, G. Finocchio, et al. A full-stack view of probabilistic computing with p-bits: Devices, architectures, and algorithms. IEEE Journal on Exploratory Solid-State Computational Devices and Circuits 9 (1), 1–11 (2023). - [22] W. Whitehead, Z. Nelson, K. Y. Camsari, and L. Theogarajan. CMOS-compatible Ising and Potts annealing using single-photon avalanche diodes. Nature Electronics 6 (12), 1009–1019 (2023). - [23] C. Duffee, J. Athas, Y. Shao, N. D. Melendez, E. Raimondo, J. A. Katine, K. Y. Camsari, G. Finocchio, and P. K. Amiri. Integrated probabilistic computer using voltage-controlled magnetic tunnel junctions as its entropy source. arXiv preprint arXiv:2412.08017 (2024). - [24] W. A. Borders, A. Z. Pervaiz, S. Fukami, K. Y. Camsari, H. Ohno, and S. Datta. Integer factorization using stochastic magnetic tunnel junctions. Nature 573 (7774), 390–393 (2019). - [25] N. S. Singh, K. Kobayashi, Q. Cao, K. Selcuk, T. Hu, S. Niazi, N. A. Aadit, S. Kanai, H. Ohno, S. Fukami, et al. CMOS plus stochastic nanomagnets enabling heterogeneous computers for probabilistic inference and learning. Nature Communications 15 (1), 2685 (2024). - [26] J. Si, S. Yang, Y. Cen, J. Chen, Y. Huang, Z. Yao, D.-J. Kim, K. Cai, J. Yoo, X. Fong, et al. Energy-efficient superparamagnetic Ising machine and its application to traveling salesman problems. Nature Communications 15 (1), 3457 (2024). - [27] J. Jhonsa, W. Whitehead, D. McCarthy, S. Chowdhury, K. Camsari, and L. Theogarajan. A CMOS Probabilistic Computing Chip With In-situ hardware Aware Learning. arXiv preprint arXiv:2504.14070 (2025). - [28] F. Böhm, D. Alonso-Urquijo, G. Verschae, and G. Van der Sande. Noise-injected analog Ising machines enable ultrafast statistical sampling and machine learning. Nature Communications 13 (1), 5847 (2022). - [29] K. Lee, S. Chowdhury, and K. Y. Camsari. Noise-augmented chaotic Ising machines for combinatorial optimization and sampling. Communications Physics 8 (1), 35 (2025). - [30] R. Faria, K. Y. Camsari, and S. Datta. Low-barrier nanomagnets as p-bits for spin logic. IEEE Magnetics Letters 8, 1–5 (2017). - [31] R. Adler. A study of locking phenomena in oscillators. Proceedings of the IRE 34 (6), 351–357 (2006). - [32] P. Bhansali and J. Roychowdhury. Gen-Adler: The generalized Adler’s equation for injection locking analysis in oscillators. In 2009 Asia and South Pacific Design Automation Conference, 522–527 (2009). - [33] G. E. P. Box, G. M. Jenkins, G. C. Reinsel, and G. M. Ljung, Time Series Analysis: Forecasting and Control, 5th ed. (John Wiley & Sons, 2015). - [34] E. M. H. E. B. Ekanayake and N. Shukla. Different paths, same destination: Designing physics-inspired dynamical systems with engineered stability to minimize the Ising Hamiltonian. Phys. Rev. Appl. 24 (2), 024008 (2025). - [35] H. Goto, K. Tatsumura, and A. R. Dixon. Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems. Science Advances 5 (4), eaav2372 (2019). - [36] S. Niazi, S. Chowdhury, N. A. Aadit, M. Mohseni, Y. Qin, and K. Y. Camsari, Training deep Boltzmann networks with sparse Ising machines, Nature Electronics 7 (7), 610–619 (2024). - [37] M. K. Bashar and N. Shukla. Designing Ising machines with higher order spin interactions and their application in solving combinatorial optimization. Scientific Reports 13 (1), 9558 (2023). - [38] D. Kleyko, D. Nikonov, A. Khosrowshahi, B. Olshausen, C. Bybee, and F. Sommer. Efficient optimization with higher-order ising machines. Nature Communications 14 (1) (2023). - [39] C. Duffee, J. Athas, A. Grimaldi, D. Volpe, G. Finocchio, E. Wei, and P. K. Amiri. Extended-variable probabilistic computing with p-dits. arXiv preprint arXiv:2506.00269 (2025). - [40] M. Honari-Latifpour and M.-A. Miri. Optical Potts machine through networks of three-photon down-conversion oscillators. Nanophotonics 9 (13), 4199–4205 (2020). - [41] N. Berloff and J. Cummins. Vector Ising Spin Annealer for Minimizing Ising Hamiltonians. (2025), Nature Portfolio. - [42] A. Mallick, M. K. Bashar, Z. Lin, and N. Shukla. Computational models based on synchronized oscillators for solving combinatorial optimization problems. Physical Review Applied 17 (6), 064064 (2022). - [43] C. Delacour, B. Haverkort, F. Sabo, N. Azemard, and A. Todri-Sanial. Lagrange Oscillatory Neural Networks for Constraint Satisfaction and Optimization. arXiv preprint arXiv:2505.07179 (2025). - [44] M. Ercsey-Ravasz and Z. Toroczkai. Optimization hardness as transient chaos in an analog approach to constraint satisfaction. Nature Physics 7 (12), 966–970 (2011).