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 ## Combined Plots: Harmonic Injection Analysis ### Overview The image presents a series of plots analyzing harmonic injection, including a circuit diagram, phase trajectory, energy landscape, and energy profiles. The plots explore the relationship between parameters like gamma (γ), epsilon (ε), and energy. ### Components/Axes **Plot (a): Circuit Diagram and Phase Trajectory** * **Circuit Diagram:** Shows a block labeled "SHI" (Second Harmonic Injection) connected to a sinusoidal wave generator labeled "FHI" (First Harmonic Injection), with an output Vout(Ø). * **Phase Trajectory:** A circular plot representing the trajectory of Ø,ε. The circle is marked with π on the left and 0 on the right. A blue dot indicates ε = 0. Arrows indicate the direction of the trajectory. The plot also shows the effect of noise and synaptic input, with arrows indicating -γ and +γ. **Plot (b): Heatmap of -∇E = ε̇** * **X-axis:** γ, ranging from -0.5 to 0.5. * **Y-axis:** ε(π), ranging from -0.5 to 0.5. * **Color Scale:** Represents -∇E = ε̇, ranging from approximately -0.5 (dark blue) to 0.5 (yellow). * **Annotations:** Includes a dashed line indicating ε̇ = 0, a yellow arrow indicating the phase trajectory resulting from synaptic input, and the label Ks = 0.15. Vertical dotted lines are present at approximately γ = -0.15 and γ = 0.15. **Plot (c): 3D Energy Landscape** * **X-axis:** ε(π), ranging from approximately -0.2 to 0.5. * **Y-axis:** γ, ranging from approximately -0.2 to 0.2. * **Z-axis:** Energy (a.u.), ranging from -0.25 to 0.1. * **Color Scale:** Represents energy, ranging from approximately -0.25 (blue) to 0.1 (yellow). **Plot (d): Energy Profiles** * **X-axis:** ε(π), ranging from -0.5 to 0.5. * **Y-axis:** Energy (a.u.), ranging from -0.30 to 0.15. * **Subplots:** Three subplots representing different values of γ. * Top subplot: Ks = 0.15, with γ = -0.1 (pink), γ = -0.15 (purple), and γ = -0.2 (green). * Middle subplot: γ = 0 (orange). * Bottom subplot: γ = 0.1 (pink), γ = 0.15 (purple), and γ = 0.2 (green). ### Detailed Analysis or ### Content Details **Plot (a):** * The phase trajectory shows a circular path, indicating the cyclical nature of the system. The synaptic input influences the trajectory, shifting it along the γ axis. **Plot (b):** * The heatmap shows the relationship between γ, ε, and -∇E = ε̇. The region around ε̇ = 0 (dashed line) represents stable states. The phase trajectory arrow indicates the system's movement in response to synaptic input. * The color gradient indicates the rate of change of ε (ε̇) with respect to the energy gradient (-∇E). Yellow regions indicate high positive values, while blue regions indicate high negative values. **Plot (c):** * The 3D energy landscape shows a potential well, with the minimum energy point representing a stable state. The shape of the well is influenced by both γ and ε. * The energy surface has a clear minimum, suggesting a stable equilibrium point for the system. **Plot (d):** * The energy profiles show how energy varies with ε for different values of γ. The curves shift and change shape depending on γ. * For Ks = 0.15, as γ decreases from -0.1 to -0.2, the energy minimum shifts downward. * For γ = 0, the energy profile is a symmetrical curve. * For positive γ values (0.1, 0.15, 0.2), the energy profiles show a similar trend to the negative γ values, with the energy minimum shifting downward as γ increases. ### Key Observations * The system exhibits a stable state characterized by a minimum energy point. * Synaptic input (γ) influences the phase trajectory and energy landscape. * The energy profiles show how the energy landscape changes with different values of γ. ### Interpretation The plots provide a comprehensive analysis of a system influenced by harmonic injection and synaptic input. The heatmap shows the dynamics of the system, with the phase trajectory indicating how the system evolves over time. The energy landscape and profiles provide insights into the stability and behavior of the system under different conditions. The data suggests that the system has a stable equilibrium point that can be influenced by external factors like synaptic input. The relationship between γ, ε, and energy is crucial in understanding the system's behavior and stability. </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. V_inj Plots ### Overview The image presents two plots, (a) and (b), each showing the probability 'p' as a function of V_inj (injection voltage). Both plots display sigmoidal curves, with 'p' increasing as V_inj increases. The plots explore the impact of different parameters on this relationship. Plot (a) shows the effect of varying K_n while keeping K_s,max constant, and plot (b) shows the effect of varying β_eff. Each plot includes a table summarizing the parameter values used for each curve. ### Components/Axes * **Plot Titles:** (a) and (b) * **Y-axis Title:** Probability, p * Scale: 0.0 to 1.0, with tick marks at 0.0, 0.5, and 1.0 * **X-axis Title:** V_inj * Plot (a) Scale: -0.8 to 0.8, with tick marks at -0.8, -0.4, 0.0, 0.4, and 0.8 * Plot (b) Scale: -1.0 to 1.0, with tick marks at -1.0, -0.5, 0.0, 0.5, and 1.0 * **Equations:** Both plots include the equation: p = (1 + tanh(β_eff \* V_inj)) / 2 * **Plot (a) Table:** | 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 | * **Plot (b) Table:** | β_eff ≈ | |---------| | 8 (Purple) | | 6 (Orange) | | 4 (Green) | | 2 (Dark Blue) | * **Constants:** * Plot (a): K_s,max = 0.15 * Plot (b): K_s,max = 0.15, K_n = 0.15 ### Detailed Analysis or ### Content Details **Plot (a):** * **Dark Blue Line (K_n = 0.01):** This line has the steepest slope and reaches a probability of 1.0 at approximately V_inj = 0.2. It starts at a probability of approximately 0.0 for V_inj < -0.4. * **Orange Line (K_n = 0.05):** This line has a shallower slope than the dark blue line. It reaches a probability of 1.0 at approximately V_inj = 0.4. * **Green Line (K_n = 0.10):** This line has a shallower slope than the orange line. It reaches a probability of 1.0 at approximately V_inj = 0.6. * **Purple Line (K_n = 0.15):** This line has a shallower slope than the green line. It reaches a probability of 1.0 at approximately V_inj = 0.8. * **Olive Line (K_n = 0.20):** This line has the shallowest slope. It reaches a probability of 1.0 at approximately V_inj = 1.0. **Plot (b):** * **Dark Blue Line (β_eff ≈ 2):** This line has the shallowest slope. It reaches a probability of 1.0 at approximately V_inj = 1.0. * **Green Line (β_eff ≈ 4):** This line has a steeper slope than the dark blue line. It reaches a probability of 1.0 at approximately V_inj = 0.6. * **Orange Line (β_eff ≈ 6):** This line has a steeper slope than the green line. It reaches a probability of 1.0 at approximately V_inj = 0.4. * **Purple Line (β_eff ≈ 8):** This line has the steepest slope. It reaches a probability of 1.0 at approximately V_inj = 0.2. ### Key Observations * In plot (a), as K_n increases, the slope of the curve decreases, indicating a less sensitive response of probability to changes in V_inj. * In plot (b), as β_eff increases, the slope of the curve increases, indicating a more sensitive response of probability to changes in V_inj. * The equation p = (1 + tanh(β_eff \* V_inj)) / 2 is used to fit the data in both plots. ### Interpretation The plots demonstrate the relationship between injection voltage (V_inj) and probability (p), influenced by parameters K_n and β_eff. Plot (a) shows that increasing K_n, while keeping K_s,max constant, reduces the sensitivity of the probability to changes in V_inj. This suggests that a higher K_n value requires a larger change in V_inj to achieve the same change in probability. Plot (b) shows that increasing β_eff increases the sensitivity of the probability to changes in V_inj. This indicates that a higher β_eff value results in a larger change in probability for the same change in V_inj. The equation used to fit the data suggests a sigmoidal relationship between V_inj and probability, which is consistent with the observed curves. The plots provide insights into how different parameters can modulate the response of a system to an input voltage. </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: Probability Distribution Comparison ### Overview The image presents a bar chart comparing the probability distribution of a system as predicted by the Boltzmann law and as measured from oscillators. The x-axis represents different states of the system, labeled using a combination of binary inputs (Ci, B, A, S, Co) and their corresponding decimal values. The y-axis represents the probability of each state. A table is included that maps the binary inputs to their decimal equivalents. ### Components/Axes * **Y-axis:** Probability, ranging from 0 to 0.2, with tick marks at 0, 0.05, 0.1, 0.15, and 0.2. * **X-axis:** \[Ci B A S Co], representing the state of the system. The x-axis is labeled with the decimal equivalents of the binary states: 0, 6, 10, 13, 18, 21, 25, 31. * **Legend:** Located on the right side of the chart. * Blue line: Boltzmann law * Orange line: Measured (oscillators) * **Table:** Located at the top-left of the chart. It maps the binary inputs (Ci, B, A, S, Co) to their decimal equivalents. The table's header row is "Ci", "B", "A", "S", "Co", and "Decimal". ### Detailed Analysis or ### Content Details **Table Data:** | Ci | B | A | S | Co | Decimal | |---|---|---|---|---|---| | 0 | 0 | 0 | 0 | 0 | 0 | | 0 | 0 | 1 | 1 | 0 | 6 | | 0 | 1 | 0 | 1 | 0 | 10 | | 0 | 1 | 1 | 0 | 1 | 13 | | 1 | 0 | 0 | 1 | 0 | 18 | | 1 | 0 | 1 | 0 | 1 | 21 | | 1 | 1 | 0 | 0 | 1 | 25 | | 1 | 1 | 1 | 1 | 1 | 31 | **Bar Chart Data:** * **State 0:** * Boltzmann law (blue): Probability ~0.12 * Measured (orange): Probability ~0.12 * **State 6:** * Boltzmann law (blue): Probability ~0.008 * Measured (orange): Probability ~0.008 * **State 10:** * Boltzmann law (blue): Probability ~0.12 * Measured (orange): Probability ~0.12 * **State 13:** * Boltzmann law (blue): Probability ~0.008 * Measured (orange): Probability ~0.008 * **State 18:** * Boltzmann law (blue): Probability ~0.008 * Measured (orange): Probability ~0.008 * **State 21:** * Boltzmann law (blue): Probability ~0.12 * Measured (orange): Probability ~0.12 * **State 25:** * Boltzmann law (blue): Probability ~0.008 * Measured (orange): Probability ~0.008 * **State 31:** * Boltzmann law (blue): Probability ~0.12 * Measured (orange): Probability ~0.12 **Trend Verification:** The Boltzmann law (blue) and Measured (orange) bars closely follow each other for each state. The probabilities are high for states 0, 10, 21, and 31, and low for states 6, 13, 18, and 25. ### Key Observations * The Boltzmann law and the measured probabilities are very similar for all states. * The states 0, 10, 21, and 31 have significantly higher probabilities compared to the other states. * The states 6, 13, 18, and 25 have very low probabilities. ### Interpretation The data suggests that the Boltzmann law accurately predicts the probability distribution of the system as measured from the oscillators. The high probabilities for states 0, 10, 21, and 31 indicate that these states are more likely to occur in the system. The low probabilities for states 6, 13, 18, and 25 suggest that these states are less likely to occur. The close agreement between the Boltzmann law and the measured probabilities validates the theoretical model and indicates that the system is behaving as expected. </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 and Cut Analysis ### Overview The image presents three subplots (a, b, and c) illustrating the behavior of oscillators and a related cut analysis over time. Subplot (a) shows the oscillator sampling sequence when FHI° > 0 and Ks = 0. Subplot (b) shows the phase (phi) of the oscillators when FHI° = 0 and Ks > 0. Subplot (c) shows the cut value over time, along with a network graph representing oscillator connections. ### Components/Axes **Top Text:** * "Oscillator sampling sequence:" * Sampling 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 **Subplot (a):** * Label: (a) * Text: FHI° > 0, Ks = 0 * Vertical axis: No explicit label, but represents oscillator state (on/off). * Horizontal axis: Time (implicit, aligned with other subplots). **Subplot (b):** * Label: (b) * Text: FHI° = 0, Ks > 0 * Vertical axis label: φ (π) * Scale: -1 to 1 * Horizontal axis: Time (implicit, aligned with other subplots). **Subplot (c):** * Label: (c) * Text: OIM (in red) * Vertical axis label: Cut * Scale: 0 to 40 * Horizontal axis label: Time * Scale: 0 to 80 **Network Graph (within subplot c):** * Nodes: Labeled 1 to 15 (representing oscillators). * Edges: Represent connections between oscillators. ### Detailed Analysis **Subplot (a): Oscillator Sampling Sequence (FHI° > 0, Ks = 0)** * This subplot shows a series of vertical bars, each representing the state of an oscillator over time. * The bars are either at a high level (on) or a low level (off). * Each oscillator has a unique color. * The oscillators switch on and off in a repeating sequence. * The sequence appears to be periodic. **Subplot (b): Oscillator Phase (FHI° = 0, Ks > 0)** * This subplot shows the phase (φ) of each oscillator over time. * The vertical axis is labeled "φ (π)" and ranges from -1 to 1. * Each oscillator is represented by a line of a specific color, matching the colors in subplot (a). * The lines show the phase evolution of each oscillator. * Some oscillators exhibit stable phases, while others oscillate or fluctuate. * There are several downward-pointing red arrows indicating a drop in phase for certain oscillators. **Individual Oscillator Phase Trends (Subplot b):** * **Black:** Stays relatively constant around 0. * **Purple:** Starts around -1, fluctuates, and then remains around -1. * **Olive:** Starts around 1, drops to 0, then fluctuates around 0. * **Teal:** Starts around 1, drops to 0, then fluctuates around 0. * **Orange:** Starts around 1, drops to 0, then fluctuates around 0. * **Magenta:** Stays relatively constant around 0. * **Brown:** Starts around 1, drops to 0, then fluctuates around 0. * **Dark Blue:** Stays relatively constant around 0. **Subplot (c): Cut Value and Network Graph** * The green line shows the "Cut" value over time. * The Cut value generally increases in steps over time. * There is a downward-pointing red arrow indicating a step increase in the Cut value. * The network graph shows connections between oscillators. * Each node in the graph is labeled with a number from 1 to 15, corresponding to the oscillators. * The edges (blue lines) represent connections between the oscillators. **Network Graph Node Positions (Approximate):** * 1: (65, 5) * 2: (75, 25) * 3: (70, 10) * 4: (70, 35) * 5: (80, 10) * 6: (65, 5) * 7: (55, 30) * 8: (55, 25) * 9: (85, 35) * 10: (75, 20) * 11: (50, 20) * 12: (75, 40) * 13: (50, 5) * 14: (65, 25) * 15: (65, 15) ### Key Observations * Subplot (a) shows a clear, repeating pattern of oscillator activation. * Subplot (b) shows varying phase behaviors among the oscillators, with some stabilizing and others fluctuating. * Subplot (c) shows a stepwise increase in the Cut value over time, suggesting a change in the network configuration. * The network graph provides a visual representation of the connections between oscillators. ### Interpretation The image presents a multi-faceted analysis of oscillator behavior and network connectivity. The oscillator sampling sequence in (a) provides a baseline activation pattern. The phase dynamics in (b) reveal how the oscillators interact and synchronize (or fail to synchronize) under specific conditions (FHI° = 0, Ks > 0). The Cut value in (c), combined with the network graph, suggests an optimization process where connections are modified to achieve a certain objective, possibly minimizing the "Cut" across the network. The red arrows highlight specific events where the phase or Cut value changes significantly, indicating key moments in the system's evolution. The relationship between the oscillator phases and the cut value is not immediately clear, but the data suggests that changes in oscillator synchronization may influence the overall network configuration and cut value. </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 Plot: Energy vs. Log Probability and Lag vs. Autocorrelation Function ### Overview The image presents two scatter plots side-by-side. Plot (a) shows the relationship between Energy and the logarithm of probability (log(p)), while plot (b) illustrates the relationship between lag and the Autocorrelation Function (ACF). Both plots display three data series, each corresponding to a different value of β (1.145, 0.662, and 0.362). The plots include legends indicating the β values and their corresponding R-squared values for plot (a). ### Components/Axes **Plot (a): Energy vs. log(p)** * **X-axis:** Energy, ranging from approximately -20 to 10. Axis markers are present at -20, -10, 0, and 10. * **Y-axis:** log(p), ranging from approximately -15 to 0. Axis markers are present at -15, -10, -5, and 0. * **Data Series:** * **Purple:** β = 1.145, R² = 0.99 * **Orange:** β = 0.662, R² = 0.977 * **Green:** β = 0.362, R² = 0.99 * **Legend:** Located in the top-right corner of plot (a). It displays the β values and their corresponding R² values in a table format. * **Plot Label:** (a) is located in the bottom-left corner of the plot. **Plot (b): Lag vs. ACF** * **X-axis:** lag, ranging from 0 to 30. Axis markers are present at 0, 10, 20, and 30. * **Y-axis:** ACF, ranging from 0.0 to 1.0. Axis markers are present at 0.0, 0.5, and 1.0. * **Data Series:** * **Purple:** β = 1.145 * **Orange:** β = 0.662 * **Green:** β = 0.362 * **Legend:** Located in the top-right corner of plot (b). It displays the β values. * **Plot Label:** (b) is located in the bottom-left corner of the plot. ### Detailed Analysis **Plot (a): Energy vs. log(p)** * **Purple (β = 1.145, R² = 0.99):** The data points are approximately: (-18, -3), (-13, -7), (-8, -10), (-3, -14). The line slopes downward. * **Orange (β = 0.662, R² = 0.977):** The data points are approximately: (-18, -4), (-13, -6), (-3, -11), (2, -16). The line slopes downward. * **Green (β = 0.362, R² = 0.99):** The data points are approximately: (-18, -5), (-13, -7), (-3, -11), (7, -14). The line slopes downward. **Plot (b): Lag vs. ACF** * **Purple (β = 1.145):** The data points are approximately: (0, 1), (5, 0.9), (10, 0.7), (15, 0.5), (20, 0.4), (25, 0.3), (30, 0.25). The line slopes downward. * **Orange (β = 0.662):** The data points are approximately: (0, 1), (5, 0.75), (10, 0.5), (15, 0.3), (20, 0.15), (25, 0.05), (30, 0.01). The line slopes downward. * **Green (β = 0.362):** The data points are approximately: (0, 1), (5, 0.7), (10, 0.4), (15, 0.2), (20, 0.08), (25, 0.02), (30, 0.001). The line slopes downward. ### Key Observations * In plot (a), as Energy increases, log(p) decreases for all values of β. * In plot (b), as lag increases, ACF decreases for all values of β. * The R-squared values for plot (a) are all very high (close to 1), indicating a strong correlation between Energy and log(p). * The rate of decrease in ACF with increasing lag varies depending on the value of β. Higher β values result in a slower decrease in ACF. ### Interpretation The plots suggest that there is a strong negative correlation between energy and the log of probability, as indicated by the high R-squared values. The autocorrelation function decreases with increasing lag, which is expected. The different β values influence the rate at which the autocorrelation decays. A higher β value corresponds to a slower decay of the autocorrelation function, suggesting a longer memory or persistence in the system being modeled. The data suggests that the system with β = 1.145 retains its autocorrelation for a longer duration compared to systems with lower β values. </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 ## Time Series Charts: Phase and Cut Parameter Evolution ### Overview The image presents two time series charts, (a) and (b), stacked vertically. Chart (a) displays the evolution of phase (φ) in units of π over time, with multiple colored lines representing different data series. Chart (b) shows the evolution of a "Cut" parameter over time, represented by a single green line. A vertical light blue shaded region highlights a period "Before Bifurcation". ### Components/Axes **Chart (a):** * **Y-axis:** φ(π), ranging from 0.0 to 1.0 in increments of 0.5. * **X-axis:** Time, implicitly ranging from 0 to 300 (based on chart (b)). * **Data Series:** Multiple colored lines (red, purple, green, orange, dark red, cyan, blue, pink) each representing a different phase evolution. * **Title:** (a) in the top-right corner. **Chart (b):** * **Y-axis:** Cut, ranging from 35 to 40 in increments of 5. * **X-axis:** Time, ranging from 0 to 300 in increments of 50. * **Data Series:** A single green line representing the "Cut" parameter. * **Title:** (b) in the top-right corner. * **Additional Text:** "DIM" in red text, located in the bottom-right. * **Annotation:** "Before Bifurcation" in blue text, above the plot, with a left-pointing arrow indicating the region. ### Detailed Analysis **Chart (a): Phase Evolution** * **Initial Phase (Time = 0):** The phase values for all series start near 0.5π. * **Before Bifurcation (Time < ~20):** Within the light blue shaded region, the phase values rapidly diverge. Some series increase towards 1.0π, while others decrease towards 0.0π. * **After Bifurcation (Time > ~20):** The phase values stabilize around either 0.0π or 1.0π, with some series exhibiting switching behavior between these two states. * **Specific Series Trends:** * **Red Line:** Starts at approximately 0.5π, increases rapidly to approximately 1.0π, and remains relatively stable with minor fluctuations. * **Purple Line:** Starts at approximately 0.5π, decreases rapidly to approximately 0.0π, and remains relatively stable with minor fluctuations. * **Green Line:** Starts at approximately 0.5π, increases rapidly to approximately 1.0π, and exhibits switching behavior, dropping to approximately 0.5π around time = 100, then returning to 1.0π. * **Orange Line:** Starts at approximately 0.5π, decreases rapidly to approximately 0.0π, and remains relatively stable with minor fluctuations. * **Dark Red Line:** Starts at approximately 0.5π, decreases rapidly to approximately 0.0π, and remains relatively stable with minor fluctuations. * **Cyan Line:** Starts at approximately 0.5π, decreases rapidly to approximately 0.0π, and exhibits switching behavior, increasing to approximately 0.5π around time = 75, then returning to 0.0π. * **Blue Line:** Starts at approximately 0.5π, increases rapidly to approximately 1.0π, and exhibits switching behavior, dropping to approximately 0.5π around time = 225, then returning to 1.0π. * **Pink Line:** Starts at approximately 0.5π, decreases rapidly to approximately 0.0π, and remains relatively stable with minor fluctuations. **Chart (b): Cut Parameter Evolution** * **Initial Value (Time = 0):** The "Cut" parameter starts at approximately 37. * **Step Changes:** The parameter exhibits step-like increases at approximately Time = 10 and Time = 50. * **Final Value (Time > ~225):** The parameter reaches a final value of approximately 39 and remains stable. ### Key Observations * **Bifurcation Point:** The light blue shaded region highlights a critical point where the system undergoes a bifurcation, leading to distinct phase states. * **Phase Locking:** After the bifurcation, the phase values tend to lock into either 0.0π or 1.0π. * **Switching Behavior:** Some phase series exhibit switching between the two locked states. * **Correlation:** The "Cut" parameter increases in steps, potentially influencing the switching behavior observed in the phase series. ### Interpretation The charts illustrate the dynamics of a system undergoing a bifurcation. Before the bifurcation, the phase values are similar. After the bifurcation, the system splits into two distinct phase states (0.0π and 1.0π), indicating a symmetry breaking. The "Cut" parameter, which increases in steps, may represent an external control parameter that influences the system's behavior and potentially triggers the switching between phase states. The "DIM" label is likely an identifier for the specific experimental or simulation setup. The data suggests that the system's behavior is sensitive to the "Cut" parameter, which drives the system towards different stable states after the initial bifurcation. </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 Chart: Energy vs. ε(π) for different Ks values ### Overview The image is a line chart displaying the relationship between Energy (in arbitrary units) and ε(π) for three different values of Ks (0.10, 0.15, and 0.20). The chart shows how energy changes as a function of ε(π) for each Ks value, with γ held constant at 0. ### Components/Axes * **X-axis:** ε(π), ranging from -0.5 to 0.5. Markers are present at -0.5, 0.0, and 0.5. * **Y-axis:** Energy (a.u.), ranging from -0.1 to 0.1. Markers are present at -0.1, 0.0, and 0.1. * **Legend (Top-Left):** * Orange line: Ks = 0.10 * Green line: Ks = 0.15 * Purple line: Ks = 0.20 * **Title (Top-Right):** γ = 0 ### Detailed Analysis * **Orange Line (Ks = 0.10):** * Trend: Starts at approximately -0.04 at ε(π) = -0.5, rises to a peak of approximately 0.06 at ε(π) = 0.0, and then decreases back to approximately -0.04 at ε(π) = 0.5. * **Green Line (Ks = 0.15):** * Trend: Starts at approximately -0.07 at ε(π) = -0.5, rises to a peak of approximately 0.08 at ε(π) = 0.0, and then decreases back to approximately -0.07 at ε(π) = 0.5. * **Purple Line (Ks = 0.20):** * Trend: Starts at approximately -0.1 at ε(π) = -0.5, rises to a peak of approximately 0.1 at ε(π) = 0.0, and then decreases back to approximately -0.1 at ε(π) = 0.5. ### Key Observations * All three lines exhibit a similar U-shaped curve, with a minimum at ε(π) = -0.5 and 0.5, and a maximum at ε(π) = 0.0. * As the value of Ks increases, the amplitude of the curve also increases. The purple line (Ks = 0.20) has the largest amplitude, while the orange line (Ks = 0.10) has the smallest. * The curves are symmetrical around ε(π) = 0.0. ### Interpretation The chart illustrates the relationship between energy and ε(π) for different values of Ks, while γ is held constant. The data suggests that as Ks increases, the energy range also increases. The symmetrical U-shape of the curves indicates a specific type of relationship between energy and ε(π), possibly related to a physical system with a periodic potential or a similar characteristic. The fact that γ = 0 might indicate a specific condition or simplification in the model being represented. </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 Charts: ε vs. Time for Varying Ks ### Overview The image presents three line charts arranged horizontally, each displaying the relationship between ε (epsilon) and Time. The charts differ in their Ks values (Ks = 0, Ks = 2, and Ks = 5), which are indicated on each chart. Each chart contains multiple lines, each likely representing a different initial condition or parameter. The initial value of epsilon is indicated as εinitial = -0.05π on the second chart. ### Components/Axes * **X-axis (Time):** All three charts share the same x-axis, labeled "Time," ranging from 0 to 6. * **Y-axis (ε):** All three charts share the same y-axis, labeled "ε," ranging from -0.5 to 0.5. A dashed horizontal line is present at ε = 0.0. * **Chart Titles:** Each chart is labeled with its corresponding Ks value: * Chart (a): Ks = 0 * Chart (b): Ks = 2 * Chart (c): Ks = 5 * **Initial Condition:** Chart (b) includes the annotation "εinitial = -0.05π" with an arrow pointing to the region where the lines originate. * **Line Colors:** The charts contain multiple lines of different colors, including magenta, teal, blue, red, black, purple, orange, brown, and green. There is no explicit legend, so the meaning of each color is not directly specified. ### Detailed Analysis **Chart (a): Ks = 0** * Trend: All lines start at different values below ε = 0 and increase over time, converging towards a value near ε = 0.5. * Specific Values: * The magenta line starts at approximately ε = -0.45 and gradually increases to approximately ε = 0.45. * The teal line starts at approximately ε = -0.35 and increases to approximately ε = 0.45. * The blue line starts at approximately ε = -0.25 and increases to approximately ε = 0.45. * The red line starts at approximately ε = -0.15 and increases to approximately ε = 0.45. * The black line starts at approximately ε = -0.05 and increases to approximately ε = 0.5. **Chart (b): Ks = 2** * Trend: All lines start at different values below ε = 0 and quickly converge to either a value near ε = 0.5 or a value near ε = -0.5. * Specific Values: * The orange line starts at approximately ε = -0.05π and decreases to approximately ε = -0.5. * The purple line starts at approximately ε = -0.15 and decreases to approximately ε = -0.5. * The black line starts at approximately ε = -0.05 and increases to approximately ε = 0.5. * The blue line starts at approximately ε = -0.25 and increases to approximately ε = 0.5. * The red line starts at approximately ε = -0.15 and increases to approximately ε = 0.5. * The green line starts at approximately ε = -0.35 and increases to approximately ε = 0.5. * The magenta line starts at approximately ε = -0.45 and decreases to approximately ε = -0.5. **Chart (c): Ks = 5** * Trend: All lines quickly converge to either a value near ε = 0.5 or a value near ε = -0.5. * Specific Values: * Most lines start near ε = -0.5 and remain there. * A few lines start near ε = 0.5 and remain there. ### Key Observations * As Ks increases, the convergence of ε to either 0.5 or -0.5 becomes more rapid. * When Ks = 0, ε gradually increases towards 0.5. * When Ks = 2, ε either increases to 0.5 or decreases to -0.5. * When Ks = 5, ε remains close to its initial value of either 0.5 or -0.5. ### Interpretation The charts illustrate the effect of the parameter Ks on the evolution of ε over time. The data suggests that Ks influences the stability and convergence of the system. When Ks is low (Ks = 0), the system gradually evolves towards a stable state. As Ks increases (Ks = 2), the system exhibits bistability, converging to one of two stable states (ε = 0.5 or ε = -0.5). When Ks is high (Ks = 5), the system becomes highly stable, remaining close to its initial state. The initial condition εinitial = -0.05π in Chart (b) highlights the starting point for some of the lines, indicating that the initial value of ε influences the final state when Ks = 2. </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 ## Chart: KL Divergence vs. τ for 5-state Adder ### Overview The image is a scatter plot showing the relationship between KL Divergence and τ for a 5-state Adder. The plot contains three data points, all represented by green diamonds. The KL Divergence appears to increase sharply from τ = 0.01 to τ = 0.03, then plateaus or slightly decreases at τ = 0.05. ### Components/Axes * **X-axis (Horizontal):** τ, with tick marks at 0.01, 0.02, 0.03, 0.04, and 0.05. * **Y-axis (Vertical):** KL Divergence, with tick marks at 0.00, 0.05, 0.10, and 0.15. * **Data Points:** Three green diamond markers. * **Text Label:** "5-state Adder" is located near the bottom-right of the plot. ### Detailed Analysis * **Data Point 1:** Located at approximately (τ = 0.01, KL Divergence = 0.00). The marker is a green diamond. * **Data Point 2:** Located at approximately (τ = 0.03, KL Divergence = 0.165). The marker is a green diamond. * **Data Point 3:** Located at approximately (τ = 0.05, KL Divergence = 0.16). The marker is a green diamond. ### Key Observations * The KL Divergence increases significantly between τ = 0.01 and τ = 0.03. * The KL Divergence plateaus or slightly decreases between τ = 0.03 and τ = 0.05. * There are only three data points in the plot. ### Interpretation The plot suggests that for the 5-state Adder, the KL Divergence is highly sensitive to changes in τ at lower values (around 0.01 to 0.03). As τ increases beyond 0.03, the KL Divergence stabilizes, indicating a potential saturation point or diminishing returns in terms of divergence. The limited number of data points makes it difficult to draw definitive conclusions about the overall trend, but it highlights a region of interest for further investigation around τ = 0.01 to 0.05. </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 ## Box Plot: Normalized Cut for Different Graphs ### Overview The image is a box plot visualizing the normalized cut values for 10 different graphs. Each graph was tested with 10 trials, and each graph has 15 nodes. The y-axis represents the normalized cut, ranging from 0.5 to 1. The x-axis represents the graph number, from 1 to 10. The box plots show the distribution of normalized cut values for each graph. ### Components/Axes * **Y-axis:** "Normalized cut", ranging from 0.5 to 1.0, with tick marks at 0.5, 0.6, 0.7, 0.8, 0.9, and 1. * **X-axis:** "Graph #", ranging from 1 to 10, with tick marks at each integer value. * **Box Plots:** Each box plot (in blue) represents the distribution of normalized cut values for a specific graph. The red line within each box represents the median. The top and bottom edges of the box represent the upper and lower quartiles. The red lines extending from the box represent the range of the data. Red dots represent outliers. * **Text:** "N = 15 nodes" and "10 trials per graph" are located in the bottom-right corner of the chart. ### Detailed Analysis Here's a breakdown of the box plots for each graph: * **Graph 1:** The box plot is a horizontal red line at approximately 0.99. * **Graph 2:** The box plot is a horizontal red line at approximately 0.99, with one outlier at approximately 0.95. * **Graph 3:** The box plot ranges from approximately 0.94 to 1.0, with a median around 0.98. * **Graph 4:** The box plot ranges from approximately 0.95 to 1.0, with a median around 0.98. * **Graph 5:** The box plot is a horizontal red line at approximately 0.99, with one outlier at approximately 0.95. * **Graph 6:** The box plot is a horizontal red line at approximately 0.99. * **Graph 7:** The box plot ranges from approximately 0.91 to 1.0, with a median around 0.98. * **Graph 8:** The box plot ranges from approximately 0.96 to 1.0, with a median around 0.98, and one outlier at approximately 0.92. * **Graph 9:** The box plot ranges from approximately 0.96 to 1.0, with a median around 0.98. * **Graph 10:** The box plot ranges from approximately 0.97 to 1.0, with a median around 0.98. ### Key Observations * The normalized cut values are generally high, clustered around 1.0 for most graphs. * Graphs 3, 4, 7, 8, 9, and 10 show more variability in their normalized cut values compared to the others. * Graphs 2, 5, and 8 have outliers, indicating some trials resulted in significantly lower normalized cut values. ### Interpretation The data suggests that for most of the tested graphs, the normalized cut is consistently high, indicating a strong separation of nodes within the graph. The variability observed in some graphs suggests that the normalized cut is more sensitive to the specific structure of those graphs. The outliers indicate that in some trials, the graph partitioning was less effective, resulting in a lower normalized cut. The fact that each graph has 15 nodes and was tested with 10 trials provides context for understanding the distribution of normalized cut values. </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 Chart: Phase Transition Over Time ### Overview The image is a line chart depicting phase transitions over time. Several lines of different colors show how a variable, phi (φ), changes from an initial value of approximately 0.5 to either 0 or 1 over a time period from 0 to 2. The chart includes the label "DIM" in the bottom-right corner. ### Components/Axes * **X-axis:** * Label: "Time" * Scale: 0 to 2, with markers at 0, 1, and 2. * **Y-axis:** * Label: "φ (π)" * Scale: 0.0 to 1.0, with markers at 0.0, 0.5, and 1.0. * **Data Series:** There are approximately 9 distinct colored lines, each representing a different condition or parameter setting. The colors include dark blue, green, pink, blue, olive, orange, red, teal, and purple. There is no explicit legend. * **Text Label:** "DIM" in the bottom-right corner, colored red. ### Detailed Analysis * **General Trend:** The lines start at φ = 0.5. Some lines transition downwards to φ = 0, while others transition upwards to φ = 1. The transitions occur within the time frame of approximately 0.2 to 1.0. * **Dark Blue Line:** This line starts at φ = 0.5 and decreases to φ = 0. The transition begins around Time = 0.2 and stabilizes at φ = 0 around Time = 0.8. * **Green Line:** This line starts at φ = 0.5 and decreases to φ = 0. The transition begins around Time = 0.3 and stabilizes at φ = 0 around Time = 0.9. * **Pink Line:** This line starts at φ = 0.5 and decreases to φ = 0. The transition begins around Time = 0.35 and stabilizes at φ = 0 around Time = 0.9. * **Blue Line:** This line starts at φ = 0.5 and decreases to φ = 0. The transition begins around Time = 0.4 and stabilizes at φ = 0 around Time = 0.95. * **Olive Line:** This line starts at φ = 0.5 and decreases to φ = 0. The transition begins around Time = 0.45 and stabilizes at φ = 0 around Time = 1.0. * **Orange Line:** This line starts at φ = 0.5 and increases to φ = 1. The transition begins around Time = 0.4 and stabilizes at φ = 1 around Time = 0.9. * **Red Line:** This line starts at φ = 0.5 and increases to φ = 1. The transition begins around Time = 0.35 and stabilizes at φ = 1 around Time = 0.8. * **Teal Line:** This line starts at φ = 0.5 and increases to φ = 1. The transition begins around Time = 0.3 and stabilizes at φ = 1 around Time = 0.75. * **Purple Line:** This line starts at φ = 0.5 and increases to φ = 1. The transition begins around Time = 0.45 and stabilizes at φ = 1 around Time = 0.9. ### Key Observations * The lines exhibit a clear bifurcation, with some transitioning to 0 and others to 1. * The transitions occur within a relatively narrow time window. * The "DIM" label is present but its significance is not immediately clear without additional context. ### Interpretation The chart likely represents a system undergoing a phase transition, where the initial state (φ = 0.5) is unstable, and the system evolves towards one of two stable states (φ = 0 or φ = 1). The different colored lines could represent different initial conditions or parameter settings that influence the direction and speed of the transition. The "DIM" label might refer to the specific system or model being studied. The data suggests that the system is sensitive to initial conditions or parameter variations, leading to the observed bifurcation. </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).