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 ## Diagram: Synaptic Input and Energy Landscape ### Overview This diagram illustrates the effect of synaptic input on a system's energy landscape, likely related to neural dynamics or oscillatory networks. It depicts a system with first and second harmonic injections, a phase space representation of its dynamics, and a 3D visualization of its energy function. The diagram is divided into four sub-panels: (a) a schematic of the input signals, (b) a phase portrait, (c) a 3D energy landscape, and (d) cross-sectional energy profiles. ### Components/Axes * **(a) Input Schematic:** * Labels: "FHI" (First Harmonic Injection), "SHI" (Second Harmonic Injection), "V_out(ø)", "Noise", "Synaptic input". * Arrows indicate signal flow. * **(b) Phase Portrait:** * Axes: ε (vertical, ranging from -0.5 to 0.5), γ (horizontal, ranging from -0.5 to 0.5). * Label: "-∇E = ε" at the top-right. * Line: ε = 0 (dashed horizontal line). * Arrows: Represent phase trajectories resulting from synaptic input. * Parameter: K_s = 0.15 * **(c) Energy Landscape:** * Axes: γ (horizontal, ranging from -0.2 to 0.5), ε (π) (vertical, ranging from -0.25 to 0.1). * Label: "Energy (a.u.)" with an arrow pointing upwards. * **(d) Energy Profiles:** * Axes: ε(π) (horizontal, ranging from -0.5 to 0.5), Energy (a.u.) (vertical, ranging from -0.3 to 0.15). * Parameter: K_s = 0.15 * Lines: Represent energy profiles for different values of γ (-0.1, -0.15, -0.2). ### Detailed Analysis or Content Details **(a) Input Schematic:** The schematic shows two input signals: First Harmonic Injection (FHI) and Second Harmonic Injection (SHI) feeding into an output signal V_out(ø). Noise and synaptic input are also indicated as influencing the system. **(b) Phase Portrait:** The phase portrait displays the system's dynamics in the ε-γ plane. The dashed line ε = 0 divides the space. The arrows indicate trajectories influenced by synaptic input. The trajectories appear to converge towards the ε = 0 line, with some exhibiting oscillatory behavior. The parameter K_s = 0.15 is noted. **(c) Energy Landscape:** The 3D plot shows the energy of the system as a function of γ and ε(π). The energy surface is complex, with a minimum around γ ≈ 0.2 and ε(π) ≈ -0.15. The color gradient indicates energy levels, with blue representing lower energy and yellow/green representing higher energy. **(d) Energy Profiles:** This panel presents cross-sectional energy profiles for different values of γ. * **γ = -0.1:** The energy profile is approximately a parabola, with a minimum around ε(π) ≈ 0. Energy values range from approximately -0.2 to 0.05. * **γ = -0.15:** The energy profile is shifted to the left, with a minimum around ε(π) ≈ -0.1. Energy values range from approximately -0.25 to 0.05. * **γ = -0.2:** The energy profile is further shifted to the left, with a minimum around ε(π) ≈ -0.2. Energy values range from approximately -0.3 to 0.05. The profiles demonstrate that the minimum energy state shifts with changes in γ. ### Key Observations * The energy landscape in (c) is not symmetric, suggesting a directional bias in the system's dynamics. * The energy profiles in (d) show a clear relationship between γ and the optimal ε(π) for minimizing energy. As γ decreases, the optimal ε(π) also decreases. * The phase portrait in (b) suggests that synaptic input can stabilize the system around the ε = 0 line. * The parameter K_s = 0.15 appears to be a constant influencing the system's behavior. ### Interpretation This diagram likely represents a model of a neural oscillator or a similar dynamical system. The FHI and SHI represent driving forces, while the synaptic input modulates the system's behavior. The energy landscape (c) provides a visualization of the system's stability and potential attractors. The minimum of the energy landscape corresponds to a stable state. The cross-sectional energy profiles (d) reveal how the stable state changes with variations in the synaptic input strength (γ). The phase portrait (b) shows how synaptic input influences the system's trajectory in the ε-γ plane. The convergence of trajectories towards the ε = 0 line suggests that synaptic input can drive the system towards a specific state. The diagram as a whole suggests that the system can adapt its dynamics based on synaptic input, potentially enabling it to perform computations or maintain stable oscillations. The parameter K_s likely represents a scaling factor or gain parameter influencing the overall system dynamics. The diagram is a theoretical exploration of the interplay between input signals, energy landscapes, and system dynamics, potentially relevant to understanding neural circuits or other oscillatory systems. </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 \n ## Charts: Probability vs. Injection Volume with Fitting Curves ### Overview The image presents two charts (labeled (a) and (b)) displaying the relationship between probability (p) and injection volume (Vinj). Both charts include fitted curves and data points, with varying parameters (Kn and βeff). A color-coded table on the right side of chart (b) relates βeff values to a color scale. ### Components/Axes * **Chart (a):** * X-axis: Vinj (Injection Volume), ranging from approximately -0.8 to 0.8. * Y-axis: Probability (p), ranging from 0.0 to 1.0. * Title: (a) * Equation: p = (1 + tanh(βeff * Vinj)) / 2 * Ks,max = 0.15 is displayed. * **Chart (b):** * X-axis: Vinj (Injection Volume), ranging from approximately -1.0 to 1.0. * Y-axis: Probability (p), ranging from 0.0 to 1.0. * Title: (b) Ks,max = 0.15, Kn = 0.15 * Equation: p = (1 + tanh(βeff * Vinj)) / 2 * **Color Table:** * Title: βeff ≈ * Values: 2, 4, 6, 8 * Color Gradient: Dark Blue -> White -> Red. * **Legend:** Located in the bottom-left of chart (b). * Kn: 0.01, βeff (from fit): 88.76 (Dark Blue) * Kn: 0.05, βeff (from fit): 12.631 (Blue) * Kn: 0.10, βeff (from fit): 6.289 (Green) * Kn: 0.15, βeff (from fit): 3.837 (Yellow) * Kn: 0.20, βeff (from fit): 2.704 (Orange) ### Detailed Analysis or Content Details **Chart (a):** * The data points (blue circles) show an S-shaped curve, increasing from approximately 0.0 at Vinj = -0.8 to approximately 1.0 at Vinj = 0.8. * A fitted curve (orange line) closely follows the data points. * The curve appears to be centered around Vinj = 0. * No specific data points are explicitly labeled, but the curve suggests p = 0.5 when Vinj is approximately 0. **Chart (b):** * The data points (blue circles) also exhibit an S-shaped curve, similar to chart (a), but with a slightly different range and shape. * A fitted curve (orange line) closely follows the data points. * The curve appears to be centered around Vinj = 0. * The legend indicates different data series corresponding to different Kn values (0.01, 0.05, 0.10, 0.15, 0.20) and their corresponding βeff values obtained from the fit. * The color of the data points corresponds to the βeff value as indicated in the color table. For example, the data points with Kn = 0.01 (dark blue) correspond to a high βeff value of 88.76. * The data points for Kn = 0.20 (orange) have a lower βeff value of 2.704. * The data points for Kn = 0.15 (yellow) have a βeff value of 3.837. * The data points for Kn = 0.10 (green) have a βeff value of 6.289. * The data points for Kn = 0.05 (blue) have a βeff value of 12.631. **Color Table:** * The color table shows a gradient from dark blue (βeff ≈ 8) to red (βeff ≈ 2). * The table provides a visual mapping between βeff values and colors. ### Key Observations * Both charts demonstrate a sigmoidal relationship between probability and injection volume. * The fitted curves accurately represent the data in both charts. * Chart (b) shows how the shape of the curve changes with different Kn values, which are correlated with βeff values. * Higher Kn values (and lower βeff values) result in a flatter S-curve. * The color table provides a convenient way to visually assess the βeff value associated with each data point. ### Interpretation The charts likely represent a model of a binary process where the probability of an event occurring (p) depends on the injection volume (Vinj). The parameter βeff controls the steepness of the sigmoid curve, indicating the sensitivity of the process to changes in Vinj. The parameter Kn appears to modulate the overall response. The data suggests that as Vinj increases, the probability of the event occurring increases rapidly, reaching a saturation point at p = 1.0. The shape of the curve is influenced by the values of Kn and βeff. The color table allows for a quick visual assessment of the βeff values associated with different data points, providing insights into the sensitivity of the process under different conditions. The relationship between Kn and βeff is inverse, with higher Kn values corresponding to lower βeff values. The fitted equation p = (1 + tanh(βeff * Vinj)) / 2 is a standard sigmoid function, commonly used to model binary outcomes. The use of this function suggests that the underlying process is likely governed by a threshold effect, where the probability of the event occurring increases sharply once Vinj exceeds a certain threshold. </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 \n ## Bar Chart: Probability Distribution Comparison - Boltzmann Law vs. Measured Oscillators ### Overview This image presents a bar chart comparing the probability distribution predicted by the Boltzmann law with the probability distribution measured from oscillators. The chart displays probability on the y-axis against a discrete variable represented by [Cᵢ B A S C₀] on the x-axis. A table is included at the top-left of the chart, mapping binary values (Cᵢ, B, A, S, C₀) to decimal values. ### Components/Axes * **X-axis:** Labeled "[Cᵢ B A S C₀]", representing a discrete variable. The axis is marked with values 0, 6, 10, 13, 18, 21, 25, and 31. * **Y-axis:** Labeled "Probability", with a scale ranging from 0 to 0.22. * **Legend:** Located in the top-right corner. * "Boltzmann law" - represented by a blue line/bars. * "Measured (oscillators)" - represented by an orange/red line/bars. * **Table:** Located in the top-left corner. Columns labeled Cᵢ, B, A, S, C₀, and Decimal. ### Detailed Analysis The chart consists of paired bars for each x-axis value, representing the Boltzmann law prediction and the measured oscillator data. * **X = 0:** Boltzmann law probability is approximately 0.002. Measured probability is approximately 0.003. * **X = 6:** Boltzmann law probability is approximately 0.12. Measured probability is approximately 0.125. * **X = 10:** Boltzmann law probability is approximately 0.12. Measured probability is approximately 0.115. * **X = 13:** Boltzmann law probability is approximately 0.02. Measured probability is approximately 0.015. * **X = 18:** Boltzmann law probability is approximately 0.12. Measured probability is approximately 0.115. * **X = 21:** Boltzmann law probability is approximately 0.02. Measured probability is approximately 0.01. * **X = 25:** Boltzmann law probability is approximately 0.02. Measured probability is approximately 0.015. * **X = 31:** Boltzmann law probability is approximately 0.005. Measured probability is approximately 0.01. The table provides the following mapping: | Cᵢ | B | A | S | C₀ | 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 | 0 | 0 | 18 | | 1 | 0 | 1 | 0 | 1 | 21 | | 1 | 1 | 0 | 0 | 1 | 25 | | 1 | 1 | 1 | 1 | 1 | 31 | ### Key Observations * The Boltzmann law and measured probabilities are generally in good agreement, particularly at x = 6 and x = 10. * There are slight discrepancies at x = 13, x = 21, and x = 25, where the measured probabilities are lower than the Boltzmann law predictions. * The probability values are highest at x = 6 and x = 10, and decrease significantly for other values. * The probability values are very low at x = 0, x = 13, x = 21, x = 25, and x = 31. ### Interpretation The chart demonstrates a comparison between the theoretical probability distribution predicted by the Boltzmann law and the experimentally measured probability distribution of oscillators. The close agreement between the two distributions suggests that the Boltzmann law is a reasonable model for the behavior of these oscillators. The slight discrepancies observed at certain values may indicate deviations from the assumptions underlying the Boltzmann law, or may be due to experimental error. The binary representation and its decimal equivalent suggest that the variable [Cᵢ B A S C₀] represents a state or configuration of the oscillators, and the probability distribution describes the likelihood of observing each state. The peaks at x=6 and x=10 indicate that these states are the most probable configurations for the oscillators. The low probabilities at the extreme values (0, 31) suggest that these configurations are rare. </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 \n ## Diagram: Oscillator Sampling Sequence & Dynamics ### Overview The image presents a series of three sub-diagrams (labeled a, b, and c) illustrating an oscillator sampling sequence and its associated dynamics. The top diagram (a) defines the sampling sequence. The middle diagram (b) shows the phase (Φ) of the oscillator over time, and the bottom diagram (c) depicts the Order-of-Interaction Metric (OIM) and the corresponding state transitions of the oscillator. The diagrams are arranged vertically, with time as the primary independent variable in diagrams (b) and (c). ### Components/Axes * **Diagram (a):** Text block defining the oscillator sampling sequence. * **Diagram (b):** * **Y-axis:** Phase (Φ) with a scale ranging from approximately -1.5 to 1.5. * **X-axis:** Time, with tick marks at intervals of approximately 10 units. * **Data Series:** A single, oscillating line representing the phase of the oscillator. * **Diagram (c):** * **Y-axis:** Order-of-Interaction Metric (OIM) with a scale ranging from approximately 25 to 40. * **X-axis:** Time, with tick marks at intervals of approximately 10 units. * **Data Series:** A step-like function representing the OIM, with several discrete jumps. * **State Transitions:** Numbered circles (1-15) representing the oscillator's states, connected by arrows indicating the sampling sequence. * **Red Arrows:** Two red arrows pointing downwards, indicating specific points in time. ### Detailed Analysis or Content Details **Diagram (a): Oscillator Sampling Sequence** The sequence is: 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 **Diagram (b): Phase (Φ)** The phase oscillates between approximately -1.2 and 1.2. The oscillation appears roughly periodic, with a period of approximately 20 time units. The amplitude of the oscillation appears to decrease slightly over time. The phase is approximately 0 at time = 0. **Diagram (c): OIM and State Transitions** * The OIM starts at approximately 35 and remains constant until approximately time = 20, where it drops to approximately 28. * The OIM then increases to approximately 38 by time = 40. * A significant drop in OIM occurs at approximately time = 45, falling to around 30. * The OIM then fluctuates, reaching a peak of approximately 38 again around time = 70. * The state transitions follow the sequence defined in diagram (a). The numbered circles represent the oscillator's state at each point in time. * The first red arrow points to a time of approximately 40, and the second red arrow points to a time of approximately 60. ### Key Observations * The OIM appears to be related to the phase of the oscillator, with drops in OIM coinciding with changes in the oscillator's state. * The sampling sequence dictates the order in which the oscillator transitions between states. * The phase oscillation appears to be dampened over time. * The red arrows highlight specific points in time where the OIM experiences a significant change. ### Interpretation The diagrams illustrate a system where an oscillator's state is sampled according to a predefined sequence. The OIM provides a measure of the complexity or interaction within the oscillator's state space. The phase diagram shows the oscillator's internal dynamics, while the OIM diagram reveals how the sampling process influences these dynamics. The drops in OIM likely correspond to transitions between different modes of oscillation or changes in the system's stability. The sampling sequence appears to drive the oscillator through a series of states, potentially exploring different regions of its state space. The decreasing amplitude of the phase oscillation suggests that the system may be losing energy or approaching a stable state. The red arrows likely indicate critical points in the sampling process where the system's behavior changes significantly. The entire system seems to be designed to explore the state space of the oscillator in a controlled manner, using the sampling sequence and OIM as tools for analysis and manipulation. </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 \n ## Scatter Plots: Energy vs. Log(p) and ACF vs. Lag ### Overview The image presents two scatter plots, labeled (a) and (b). Plot (a) shows the relationship between Energy and log(p), while plot (b) depicts the Autocorrelation Function (ACF) against Lag. Each plot contains three data series, distinguished by color and represented with associated β and R² values in a legend. ### Components/Axes **Plot (a): Energy vs. Log(p)** * **X-axis:** Energy, ranging from approximately -20 to 10. * **Y-axis:** log(p), ranging from approximately -15 to 0. * **Data Series 1:** Purple, β = 1.145, R² = 0.99 * **Data Series 2:** Orange, β = 0.662, R² = 0.977 * **Data Series 3:** Green, β = 0.362, R² = 0.99 **Plot (b): ACF vs. Lag** * **X-axis:** Lag, ranging from approximately 0 to 30. * **Y-axis:** ACF, ranging from approximately 0.0 to 1.0. * **Data Series 1:** Purple, β = 1.145 * **Data Series 2:** Orange, β = 0.662 * **Data Series 3:** Green, β = 0.362 Both plots share a common legend positioned in the top-right corner. ### Detailed Analysis or Content Details **Plot (a): Energy vs. Log(p)** * **Purple Data Series:** The data points exhibit a strong downward trend. Starting at approximately log(p) = -1.0 when Energy = -20, the series decreases to approximately log(p) = -14.5 when Energy = 10. The points are tightly clustered around a decreasing line. * **Orange Data Series:** This series also shows a downward trend, but less pronounced than the purple series. Starting at approximately log(p) = -3.5 when Energy = -20, it decreases to approximately log(p) = -13.5 when Energy = 10. The points are more scattered than the purple series. * **Green Data Series:** This series exhibits a similar downward trend to the orange series, but with a shallower slope. Starting at approximately log(p) = -2.5 when Energy = -20, it decreases to approximately log(p) = -14.0 when Energy = 10. The points are relatively well-aligned. **Plot (b): ACF vs. Lag** * **Purple Data Series:** The ACF starts at approximately 1.0 at Lag = 0 and decays rapidly to approximately 0.2 at Lag = 30. The decay is relatively slow initially, then accelerates. * **Orange Data Series:** The ACF starts at approximately 1.0 at Lag = 0 and decays to approximately 0.1 at Lag = 30. The decay is faster than the purple series. * **Green Data Series:** The ACF starts at approximately 1.0 at Lag = 0 and decays to approximately 0.05 at Lag = 30. This series exhibits the fastest decay among the three. ### Key Observations * In both plots, the purple data series consistently exhibits the strongest trend (steepest slope in (a) and fastest decay in (b)). * The R² values in plot (a) indicate a strong linear relationship between Energy and log(p) for all three series, with the purple series having the highest R² value (0.99). * The β values are different for each series, indicating different rates of change. * The ACF values in plot (b) all start at 1.0, as expected for autocorrelation at lag 0. ### Interpretation The plots likely represent the analysis of a time series or a similar sequential dataset. Plot (a) suggests a power-law relationship between Energy and probability (represented by log(p)), where higher energy levels correspond to lower probabilities. The different β values indicate that the rate of this decrease varies depending on the specific series. The high R² values suggest that this relationship is well-modeled by a linear function in log-log space. Plot (b) shows the autocorrelation function, which measures the correlation between the series and its lagged values. The rapid decay of the ACF indicates that the series is not strongly autocorrelated, meaning that values at different time points are largely independent. The different decay rates for the three series suggest that they have different levels of temporal dependence. The purple series, with the slowest decay, exhibits the strongest autocorrelation. The combination of these two plots provides insights into the statistical properties of the underlying data. The power-law relationship in plot (a) suggests a scale-free behavior, while the ACF in plot (b) reveals the temporal structure of the data. The differences between the three series suggest that they may represent different aspects of the same phenomenon or different phenomena altogether. </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 ## Chart: Bifurcation Diagram and Cut Value Over Time ### Overview The image presents two time-series plots, labeled (a) and (b). Plot (a) displays a series of fluctuating lines representing a quantity denoted as Φ (in radians, π), plotted against Time. Plot (b) shows a step-like function representing "Cut" value, also plotted against Time. A red arrow indicates a point of interest in both plots. The text "Before Bifurcation" is placed above the initial portion of plot (a). The label "DIM" is present in the bottom-right corner of the image. ### Components/Axes * **Plot (a):** * X-axis: Time (units not specified, ranging from 0 to 300) * Y-axis: Φ (π) (ranging from 0 to 1.0) * Multiple lines representing different trajectories or conditions. The lines are colored: blue, orange, yellow, red, purple, green, and several shades of red/brown. * **Plot (b):** * X-axis: Time (units not specified, ranging from 0 to 300) * Y-axis: Cut (ranging from approximately 32 to 40) * A single green line representing the Cut value. * **Annotations:** * Red arrow pointing to approximately Time = 225 in both plots. * Text: "Before Bifurcation" positioned above the initial portion of plot (a). * Text: "DIM" in the bottom-right corner. ### Detailed Analysis or Content Details **Plot (a): Φ vs. Time** The plot shows multiple lines fluctuating between approximately 0 and 1.0. Initially (before Time = 50), the lines exhibit relatively high frequency oscillations. After Time = 50, the lines settle into more stable, but still fluctuating, patterns. The lines are not all identical; there is variation in their amplitude and frequency. * **Blue Line:** Starts at approximately Φ = 0.8, drops to near 0, then fluctuates between 0.2 and 0.6. * **Orange Line:** Starts at approximately Φ = 0.7, drops to near 0, then fluctuates between 0.2 and 0.7. * **Yellow Line:** Starts at approximately Φ = 0.6, drops to near 0, then fluctuates between 0.1 and 0.5. * **Red Line:** Starts at approximately Φ = 0.5, drops to near 0, then fluctuates between 0.1 and 0.4. * **Purple Line:** Starts at approximately Φ = 0.4, drops to near 0, then fluctuates between 0.1 and 0.3. * **Green Line:** Starts at approximately Φ = 0.3, drops to near 0, then fluctuates between 0.1 and 0.2. * **Brown/Red Lines:** Several lines with similar behavior, fluctuating between 0.1 and 0.5. Around Time = 225 (indicated by the red arrow), there appears to be a change in the behavior of the lines, with some lines exhibiting larger fluctuations. **Plot (b): Cut vs. Time** The plot shows a step-like function. The Cut value starts at approximately 37.5, drops to approximately 33 at Time = 0, remains relatively constant until Time = 200, then increases to approximately 39 at Time = 200. The value remains at approximately 39 until Time = 300. The red arrow indicates the step change at Time = 200. ### Key Observations * The "Before Bifurcation" label suggests that the initial portion of plot (a) represents a state before a bifurcation event. * The red arrow in both plots highlights a potential correlation between the change in the Cut value and a change in the behavior of the Φ values. * The lines in plot (a) do not converge to a single value, indicating a complex system with multiple possible states. * The Cut value in plot (b) is a control parameter that appears to influence the behavior of the system represented by plot (a). ### Interpretation The data suggests a system undergoing a bifurcation, where a small change in a control parameter ("Cut") leads to a significant change in the system's behavior (represented by the fluctuating Φ values). The "Before Bifurcation" label indicates that the initial state is relatively stable, but as the Cut value is changed, the system enters a more complex and potentially chaotic regime. The multiple lines in plot (a) represent different possible trajectories or states of the system, and the fluctuations in Φ indicate sensitivity to initial conditions. The DIM label is unclear without further context, but it may refer to a specific dimension or parameter of the system. The correlation between the step change in Cut and the change in Φ suggests that the Cut value is a key driver of the bifurcation. The system appears to be sensitive to changes in the Cut parameter, leading to a shift in the dynamics of the system. </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 \n ## Chart: Energy vs. Strain for Different K_s Values ### Overview The image presents a line chart illustrating the relationship between Energy (in arbitrary units, a.u.) and Strain (ε, in units of π) for three different values of K_s (0.10, 0.15, and 0.20) with γ = 0. The chart depicts the energy landscape as a function of strain, showing how the energy changes with deformation for each K_s value. ### Components/Axes * **X-axis:** Strain (ε) measured in units of π. The axis ranges from approximately -0.55 to 0.55. * **Y-axis:** Energy measured in arbitrary units (a.u.). The axis ranges from approximately -0.12 to 0.11. * **Legend:** Located in the top-left corner, the legend identifies the three lines by their corresponding K_s values and colors: * K_s = 0.10 (Orange) * K_s = 0.15 (Green) * K_s = 0.20 (Purple) * **Annotation:** γ = 0 is written in the top-right corner. ### Detailed Analysis The chart displays three curves, each representing a different K_s value. * **K_s = 0.10 (Orange Line):** This line exhibits a symmetrical parabolic shape. The minimum energy is located at ε = 0. The energy at ε = 0 is approximately 0.00 a.u. The energy at ε = -0.5 and ε = 0.5 is approximately 0.08 a.u. * **K_s = 0.15 (Green Line):** This line also shows a symmetrical parabolic shape, but the curve is narrower than the orange line. The minimum energy is located at ε = 0. The energy at ε = 0 is approximately 0.00 a.u. The energy at ε = -0.5 and ε = 0.5 is approximately 0.06 a.u. * **K_s = 0.20 (Purple Line):** This line is the narrowest of the three, indicating a steeper energy landscape. The minimum energy is located at ε = 0. The energy at ε = 0 is approximately 0.00 a.u. The energy at ε = -0.5 and ε = 0.5 is approximately 0.04 a.u. All three lines intersect at the point (0, 0), where the strain is zero and the energy is zero. ### Key Observations * As K_s increases, the energy landscape becomes steeper, meaning that a larger force is required to deform the material by the same amount. * The minimum energy for all three K_s values occurs at zero strain (ε = 0). * The curves are all symmetrical around the y-axis, indicating that the energy landscape is the same for positive and negative strain. ### Interpretation The chart demonstrates the relationship between energy and strain for different values of a parameter K_s, likely representing a stiffness or spring constant. The parabolic shapes suggest a simple harmonic oscillator-like behavior, where the energy increases quadratically with strain. The higher the K_s value, the stiffer the system, and the more energy is required to achieve a given strain. The γ = 0 annotation suggests that damping is absent in this model. The data suggests that increasing K_s leads to a more resistant material, requiring more energy to deform. The symmetry of the curves indicates that the material behaves identically under tensile and compressive strain. This type of plot is common in solid-state physics or materials science to model the potential energy surface of a crystal lattice or a material under stress. </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 Chart: Oscillation Damping vs. Time with Varying Spring Constants ### Overview The image presents three line charts, arranged horizontally, illustrating the damping of oscillations over time. Each chart corresponds to a different spring constant (Ks) value: 0, 2, and 5. The y-axis represents the angular frequency (ω), and the x-axis represents time. Each chart displays multiple lines representing different initial conditions. A label indicates an initial condition of εinitial = -0.05π. ### Components/Axes * **Y-axis Label:** ω (Angular Frequency) - Scale ranges from approximately -0.5 to 0.5. * **X-axis Label:** Time - Scale ranges from 0 to 6. * **Chart (a):** Ks = 0. Displays multiple oscillating lines. * **Chart (b):** Ks = 2. Displays multiple lines that quickly dampen to around 0. * **Chart (c):** Ks = 5. Displays lines that almost immediately settle at 0. * **Annotation:** εinitial = -0.05π, with an arrow pointing to a line in chart (a). * **Sub-labels:** (a), (b), (c) denoting the different spring constant values. ### Detailed Analysis or Content Details **Chart (a) - Ks = 0:** * **Trend:** Multiple lines exhibit sustained oscillations with varying amplitudes and frequencies. * **Data Points (approximate):** * Black Line: Starts at ~0.45, oscillates around 0, with a period of ~1.5. * Red Line: Starts at ~0.45, oscillates around 0, with a period of ~1.5. * Green Line: Starts at ~-0.45, oscillates around 0, with a period of ~1.5. * Blue Line: Starts at ~-0.45, oscillates around 0, with a period of ~1.5. * Purple Line: Starts at ~0.45, oscillates around 0, with a period of ~1.5. * Cyan Line: Starts at ~-0.45, oscillates around 0, with a period of ~1.5. * The line indicated by the annotation (εinitial = -0.05π) is a cyan line, starting at approximately -0.45. **Chart (b) - Ks = 2:** * **Trend:** Lines quickly dampen towards ω = 0, with some initial oscillations. The damping is significantly faster than in chart (a). * **Data Points (approximate):** * Black Line: Starts at ~0.45, quickly decays to 0 within ~2 time units. * Red Line: Starts at ~0.45, quickly decays to 0 within ~2 time units. * Green Line: Starts at ~-0.45, quickly decays to 0 within ~2 time units. * Blue Line: Starts at ~-0.45, quickly decays to 0 within ~2 time units. * Purple Line: Starts at ~0.45, quickly decays to 0 within ~2 time units. * Cyan Line: Starts at ~-0.45, quickly decays to 0 within ~2 time units. **Chart (c) - Ks = 5:** * **Trend:** Lines almost immediately settle at ω = 0, indicating very rapid damping. * **Data Points (approximate):** * Black Line: Starts at ~0.45, almost immediately decays to 0. * Red Line: Starts at ~0.45, almost immediately decays to 0. * Green Line: Starts at ~-0.45, almost immediately decays to 0. * Blue Line: Starts at ~-0.45, almost immediately decays to 0. * Purple Line: Starts at ~0.45, almost immediately decays to 0. * Cyan Line: Starts at ~-0.45, almost immediately decays to 0. ### Key Observations * Increasing the spring constant (Ks) dramatically reduces the oscillation period and increases the damping rate. * When Ks = 0, the system exhibits sustained oscillations. * As Ks increases, the oscillations are quickly suppressed, and the system returns to equilibrium (ω = 0). * The initial condition (εinitial = -0.05π) does not appear to significantly affect the overall damping trend, only the initial direction of the oscillation. ### Interpretation The charts demonstrate the effect of spring constant (Ks) on the damping of oscillations. A higher spring constant introduces a stronger restoring force, which dissipates energy more rapidly, leading to faster damping. The system transitions from undamped oscillations (Ks = 0) to critically damped or overdamped behavior (Ks = 5). The initial condition influences the starting point of the oscillation but doesn't alter the fundamental damping behavior dictated by the spring constant. This suggests a model of a damped harmonic oscillator where the spring constant plays a crucial role in controlling the system's stability and response to disturbances. The consistent behavior across different initial conditions indicates the system's robustness to small variations in starting state. </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 \n ## Scatter Plot: KL Divergence vs. Tau for 5-state Adder ### Overview This image presents a scatter plot illustrating the relationship between KL Divergence and a parameter denoted as 'τ' (tau) for a "5-state Adder". The plot displays three data points, each represented by a green diamond marker. ### Components/Axes * **X-axis:** Labeled "τ" (tau). The scale ranges from approximately 0.00 to 0.06, with tick marks at 0.01, 0.02, 0.03, 0.04, and 0.05. * **Y-axis:** Labeled "KL Divergence". The scale ranges from approximately -0.02 to 0.16, with tick marks at 0.00, 0.05, 0.10, and 0.15. * **Title/Label:** "5-state Adder" is positioned in the bottom-right corner of the plot area. * **Data Markers:** Green diamond shapes represent individual data points. ### Detailed Analysis The plot contains three data points: 1. **Point 1:** Located at approximately τ = 0.01 and KL Divergence = 0.00. 2. **Point 2:** Located at approximately τ = 0.03 and KL Divergence = 0.16. 3. **Point 3:** Located at approximately τ = 0.05 and KL Divergence = 0.16. The data points show a non-linear relationship. KL Divergence remains at approximately 0.00 for τ = 0.01, then rapidly increases to approximately 0.16 for τ = 0.03 and remains constant at 0.16 for τ = 0.05. ### Key Observations * The KL Divergence is minimal at the lowest value of τ (0.01). * There is a significant jump in KL Divergence between τ = 0.01 and τ = 0.03. * KL Divergence plateaus at approximately 0.16 for τ values of 0.03 and 0.05. * The data is sparse, with only three points, making it difficult to determine the exact nature of the relationship. ### Interpretation The plot suggests that the KL Divergence, a measure of how one probability distribution differs from a reference probability distribution, is sensitive to changes in the parameter τ for the 5-state adder. Initially, small changes in τ do not significantly affect the KL Divergence. However, beyond a certain threshold (around τ = 0.03), further increases in τ do not lead to further increases in KL Divergence, indicating a saturation point. This could imply that the 5-state adder's behavior changes significantly around τ = 0.03. Below this value, the system operates in a regime where the probability distribution is relatively stable. Above this value, the system reaches a state where further adjustments to τ do not alter the probability distribution significantly. The high KL Divergence values suggest a substantial difference between the probability distribution at τ = 0.03 and 0.05 and the reference distribution (presumably at τ = 0.01). Further investigation with more data points would be needed to fully characterize the relationship and understand the underlying reasons for this behavior. The choice of KL Divergence as the metric suggests an interest in quantifying the information loss or difference in distributions as τ is varied. </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 vs. Graph Number ### Overview The image presents a series of box plots illustrating the distribution of "Normalized Cut" values across ten different graphs. Each graph was generated using 15 nodes and 10 trials. The box plots display the median, quartiles, and outliers for each graph's normalized cut values. ### Components/Axes * **X-axis:** "Graph #" ranging from 1 to 10. * **Y-axis:** "Normalized cut" ranging from 0.5 to 1.0. * **Box Plot Elements:** Each box plot represents the distribution of normalized cut values for a single graph. The box represents the interquartile range (IQR), the line inside the box represents the median, and the whiskers extend to the most extreme data points within 1.5 times the IQR. Points beyond the whiskers are considered outliers and are plotted individually. * **Annotation:** "N = 15 nodes\n10 trials per graph" located in the bottom-right corner. ### Detailed Analysis The box plots are arranged horizontally, one for each graph number from 1 to 10. * **Graph 1:** The median normalized cut is approximately 0.98. The IQR extends from roughly 0.96 to 1.0. There is one outlier at approximately 0.94. * **Graph 2:** The median normalized cut is approximately 0.98. The IQR extends from roughly 0.96 to 1.0. There is one outlier at approximately 0.94. * **Graph 3:** The median normalized cut is approximately 0.98. The IQR extends from roughly 0.96 to 1.0. There is one outlier at approximately 0.94. * **Graph 4:** The median normalized cut is approximately 0.98. The IQR extends from roughly 0.96 to 1.0. There is one outlier at approximately 0.94. * **Graph 5:** The median normalized cut is approximately 0.98. The IQR extends from roughly 0.96 to 1.0. There is one outlier at approximately 0.94. * **Graph 6:** The median normalized cut is approximately 0.98. The IQR extends from roughly 0.96 to 1.0. * **Graph 7:** The median normalized cut is approximately 0.94. The IQR extends from roughly 0.92 to 0.97. * **Graph 8:** The median normalized cut is approximately 0.98. The IQR extends from roughly 0.96 to 1.0. There is one outlier at approximately 0.94. * **Graph 9:** The median normalized cut is approximately 0.98. The IQR extends from roughly 0.96 to 1.0. * **Graph 10:** The median normalized cut is approximately 0.98. The IQR extends from roughly 0.96 to 1.0. ### Key Observations * Graphs 1, 2, 3, 4, 5, 8, 9, and 10 exhibit very similar distributions of normalized cut values, with medians close to 1.0 and minimal variation. * Graph 7 shows a noticeably lower median normalized cut value (approximately 0.94) compared to the other graphs. * Outliers are present in Graphs 1, 2, 3, 4, 5, and 8, all falling below 0.95. * The normalized cut values are consistently high across most graphs, suggesting a strong degree of connectivity or separation within the graphs. ### Interpretation The data suggests that the normalized cut values are generally high and consistent across most of the graphs generated. The normalized cut is a measure of the "strength of cut" between two sets of nodes in a graph. A higher normalized cut indicates a stronger separation between the sets. The consistency across graphs 1, 2, 3, 4, 5, 8, 9, and 10 suggests that the underlying graph structures are similar in terms of their connectivity patterns. Graph 7 stands out as an anomaly, exhibiting a lower median normalized cut. This could indicate that the graph structure for Graph 7 is different from the others, potentially having weaker connections between node sets or a different overall topology. The outliers observed in some graphs may represent unusual configurations or edge cases within those specific graph instances. The annotation "N = 15 nodes, 10 trials per graph" indicates that each graph consists of 15 nodes and that the normalized cut was calculated based on 10 independent trials for each graph. This suggests that the observed distributions represent the variability in normalized cut values due to random factors or variations in the graph generation process. </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: Time Series of Φ(π) ### Overview The image presents a line chart depicting the relationship between "Time" on the x-axis and "Φ(π)" on the y-axis. Multiple lines are plotted, representing different time series or conditions. The chart appears to model a dynamic process where Φ(π) changes over time, exhibiting a complex pattern of increase and decrease. A label "DIM" is present in the bottom-right corner. ### Components/Axes * **X-axis:** Labeled "Time", ranging from approximately 0 to 2. * **Y-axis:** Labeled "Φ(π)", ranging from approximately 0 to 1. * **Lines:** Multiple colored lines are plotted, each representing a different time series. The colors observed are: black, green, purple, blue, cyan, magenta, and orange. * **Label:** "DIM" is located in the bottom-right corner of the chart. ### Detailed Analysis The chart contains several lines, each representing a different trajectory of Φ(π) over time. * **Black Line:** Starts at approximately Φ(π) = 0.5, remains relatively constant until Time ≈ 0.7, then decreases to near 0 by Time ≈ 1.2, and remains near 0 until Time = 2. * **Green Line:** Starts at approximately Φ(π) = 0.5, increases rapidly to a peak around Φ(π) = 0.95 at Time ≈ 0.6, then decreases to near 0 by Time ≈ 1.2, and remains near 0 until Time = 2. * **Purple Line:** Starts at approximately Φ(π) = 0.5, decreases rapidly to near 0 by Time ≈ 0.8, and remains near 0 until Time = 2. * **Blue Line:** Starts at approximately Φ(π) = 0.5, increases to a peak around Φ(π) = 0.95 at Time ≈ 0.6, then decreases to near 0 by Time ≈ 1.2, and remains near 0 until Time = 2. * **Cyan Line:** Starts at approximately Φ(π) = 0.5, increases to a peak around Φ(π) = 0.95 at Time ≈ 0.6, then decreases to near 0 by Time ≈ 1.2, and remains near 0 until Time = 2. * **Magenta Line:** Starts at approximately Φ(π) = 0.5, increases to a peak around Φ(π) = 0.95 at Time ≈ 0.6, then decreases to near 0 by Time ≈ 1.2, and remains near 0 until Time = 2. * **Orange Line:** Starts at approximately Φ(π) = 0.5, increases to a peak around Φ(π) = 0.95 at Time ≈ 0.6, then decreases to near 0 by Time ≈ 1.2, and remains near 0 until Time = 2. Most lines exhibit a similar pattern: an initial value around 0.5, a rapid increase to near 1, a subsequent rapid decrease to near 0, and then a plateau near 0. The timing of the peak and the rate of decrease vary slightly between the lines. ### Key Observations * The lines generally follow a similar trend, suggesting a common underlying process. * There is variation in the timing of the peak and the rate of decay, indicating different conditions or parameters. * The "DIM" label's significance is unknown without further context. It could represent a specific condition, parameter setting, or data subset. * The y-axis is labeled with a function Φ(π), suggesting a mathematical or physical context. ### Interpretation The chart likely represents a dynamic system where a quantity Φ(π) evolves over time. The initial value of 0.5 suggests an initial state, and the subsequent increase and decrease indicate a transient process. The plateau near 0 at later times suggests a stable state or equilibrium. The variations between the lines could be due to different initial conditions, parameter settings, or external influences. The "DIM" label might indicate a specific dimension or aspect of the system being analyzed. Without additional context, it's difficult to determine the precise meaning of Φ(π) and the underlying physical or mathematical process. However, the chart provides valuable insights into the system's behavior and the factors that influence it. The consistent pattern across the lines suggests a robust underlying mechanism, while the variations highlight the sensitivity of the system to different conditions. </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).