Probabilistic computing with NbOx metal-insulator transition-based self-oscillatory pbit

Energy-based computing is a promising approach for addressing the rising demand for solving NP-hard problems across diverse domains, including logistics, artificial intelligence, cryptography, and optimization. Probabilistic computing utilizing pbits, which can be manufactured using the semiconductor process and seamlessly integrated with conventional processing units, stands out as an efficient candidate to meet these demands. Here, we propose a novel pbit unit using an NbOx volatile memristor-based oscillator capable of generating probabilistic bits in a self-clocking manner. The noise-induced metal-insulator transition causes the probabilistic behavior, which can be effectively modeled using a multi-noise-induced stochastic process around the metal-insulator transition temperature. We demonstrate a memristive Boltzmann machine based on our proposed pbit and validate its feasibility by solving NP-hard problems. Furthermore, we propose a streamlined operation methodology that considers the autocorrelation of individual bits, enabling energy-efficient and high-performance probabilistic computing.


Area-dependent device capacitance and current characteristics
Fig. S1. a C-V measurement result of NbO x memristor with diameter of 40 nm.b Capacitances of NbO x memristors with different diameter (100 nm, 60 nm, 40 nm) before (pristine) and after elecroforming.c NDR behaviors with different diameter (100 nm, 60 nm, 40 nm) after electroforming.They exhibited almost identical NDR behaviors considering the area variation, especially at the off state.d Device conductances before (pristine) and after elecroforming at the off state read (at -0.1 V).These results strongy support the core-shell model, where the localized core dominates the NDR behaivors.Some variations in NDR behaviors were noted, likely attributable to shell characteristics, including shell area and conductivity.

Fig. S2. Probabilistic oscillation dataset of NbO x oscillator.
Five tests were performed for each external voltage (V ext ) condition.For each test, the current oscillations were measured under 30 μs width of a single V ext pulse.The data collection process for V ext cases involved collecting data once for all V ext cases initially and then repeating the process four more times for each V ext case individually.The measurement system we used to evaluate the NbO x oscillator is illustrated in Fig. S3a.Two PMU cards connected to a Keithley 4200 SCS semiconductor analyzer apply voltage to both terminals of the NbO x memristor and read the output.PMU 1 is connected to a series resistor module, and the probe makes contact with the top electrode of the NbO x memristor.Another probe contacts the bottom electrode and is directly connected to PMU 2. All components, except for the semiconductor analyzer, are located within a dark box.

Detailed description of the measurement system's effect on device characterization
Based on this setup, an equivalent NbO x oscillator circuit is presented in Fig. S3b. Circuit elements corresponding to components in the measurement system are colored to match those in Fig. S3a.We used a noise-free NbO x memristor model for the simulation, and the results under various conditions are summarized in Table S1.
For case 1, we configured the oscillator circuit identically to the one described in the main text, setting only the load resistor (R L ) to 1200 Ω and removing all other elements.When a DC voltage (V ext ) of 1.32 V and 1.33 V was applied, we observed oscillations in the output current, consistent with Fig. 2d in the main text.
In case 2, we introduced a parasitic capacitance (C p ) to reflect the electrically isolated state of the probe in the measurement system.As a result, the node voltage at the front end of the NbO x device becomes time-dependent due to C p .Under these conditions, we observed oscillations at V ext = 1.32 V but not at 1.33 V, despite the total series resistance being the same 1200 Ω as in case 1.This indicates that C p influences the oscillation conditions of the NbO x oscillator.In particular, it can be seen that the nonoscillation of V ext = 1.33 V, despite the circuit being configured with the same series resistor, is entirely an effect of C p compared to case 1.
This trend is even more severe in case 3, which incorporates the resistances of the top and bottom electrodes as well as the contact resistance of the probe, thus closely mimicking the actual experimental setup.In this case, the circuit did not oscillate even at V ext = 1.32 V, so it can be seen that the effect of voltage dividing due to the series resistances further increases the influence of the C p on oscillation.
From these comparative simulation results, in a real-world experimental setup, the presence of C p has the effect of making the system to converge rather than oscillate from a lower point than when calculated considering only the equilibrium (static) state ( ) of the system.Consequently, even though the equilibrium point is formed at a point slightly below the NDR-2 region, as shown in the red curves in Figure 1a, it can be understood that probabilistic oscillations can occur in the measurement environment due to the resulting non-oscillatory state caused by C p .At the given load line, the operating point (blue dot) is at the positive differential resistance region above the NDR-2 region, so the self-oscillator is not oscillating but stays at the on-state.
At the equilibrium state, the perturbation may change the operating point to the metastable space. 1,2nce it happens, the time-dependent temperature gradient (dT m /dt) and voltage gradient (dv m /dt) at the point act as a restoring force that drives it to an equilibrium point.Fig. S5c-d show the calculated dv m /dt and dT m /dt vectors on the T-V space, respectively.The black lines are equilibrium lines giving dv m /dt and dT m /dt = 0, which are equal to the load line and equilibrium T-V curve.In Fig. S5d, the dT m /dt was abruptly changed at the metal-insulator transition temperature (T MIT ) due to the drastic change of R th between the metallic phase (R th,M ) and the insulating phase (R th,I ).By using the two gradient vectors, the trajectory from the metastable state to the equilibrium state over time can be estimated.Fig. S5e shows the schematics of two perturbation results in a completely different trajectories; one is a larger voltage perturbation case (black arrow), and the other is a smaller voltage perturbation case (red arrow).The  0 , ( 0 ,  0 ) is the equilibrium point which is defined by the intersection of dv m /dt and dT m /dt = 0 at a given V ext .When external voltage perturbation (v pert , blue vector) is engaged to t 0 , it tends to go to  0 , ( 0 +  pert ,  0 ).However, due to the restore forces (yellow vector), the metastable state after  can be defined by  1 ( , ) (black vector).As such, the sum of perturbation and restore forces gives the metastable points trajectory toward the equilibrium point over time.This results in a spiral  n ( n ,  n ) trajectory at the T-V space.Here, when the perturbation was small, the metastable state temperature was kept higher than the T MIT and it eventually converged to the initial equilibrium state (red trajectory).Whereas, when the perturbation was large, the metastable temperature dropped below the T MIT , so the temperature gradient drastically increased negatively due to the metallic-to-insulator transition.This resulted in a cooling down of the NbO x device to the off-state (black trajectory).Then, due to the positive voltage gradient at the off-state, the  increased until it exceeded the V th , and the device returned to the on-state, generating an oscillation peak.Then, it converged to the initial equilibrium state.As a result, perturbation-induced oscillation could be possible when the perturbation was sufficient.In addition to voltage perturbation, temperature perturbation can be considered similarly.
Figs. S5f-h show the perturbation-induced oscillation generation by intentionally applying a perturbing pulse in the numerical simulation model.Before that, we investigated the boundary condition of V ext between the periodic oscillation and no oscillation as shown in Fig. S5f.The model confirmed that the V ext under 1.33 V generated the periodic oscillation, and over 1.40 V showed no oscillation at the onstate.Next, we set the device to the on-state (V ext = 1.40 V) and then inserted the following three perturbing pulses: 1.30 V for 5 ns (P1, red), 1.33 V for 5 ns (P2, green), and 1.30 V for 1 ns (P3, blue) as shown in the top panel of Fig. S5g.The lower panel shows the output currents for three perturbating cases.Despite all voltage conditions corresponding to the oscillation condition, only P1 generated an oscillation.Fig. S5h shows the   (  ,   ) trajectory of the three cases.The perturbation by P1 was sufficiently large to cause the temperature to go down below the T MIT .Whereas the perturbation by P2 and P3 was weak, making the non-equilibrium state converge to the equilibrium point via a small spiral trajectory.Thus, these simulation results show that T MIT corresponds to the threshold for the occurrence of oscillation generation for the metastable state of the NbO x oscillator.

Probabilistic oscillation result of p-osc model
Fig. S6. a Probabilistic oscillation result of the numerical oscillator model with the OU process in the same case of Fig. 3 in main text as increasing the V ext from 1.25 V to 1.49 V with a 0.02 V interval.b Oscillation probability distribution calculated from the above oscillation results.

Derivation of general solution for the simultaneous equations of the OU process of the NbO x oscillator
The two stochastic differential equations (SDE) of the Ornstein-Uhlenbeck (OU) process are given as Equations 5 and 6 in main text: Equation 5: From Equation S9, we can calculate the probability distribution.Because Wiener process is the definite integral of a white noise Gaussian process, it is characterized by having zero mean, unit variance 3,4 .Then the mean and variance of Equation S9 are calculated as below: , where ϵ has a standard normal distribution (0, 1).As a consequence, the temperature variation induced by voltage noise in infinitesimal time interval  is described as below equation.In addition, the Wiener process is defined as a Markov process with zero mean and variance of ∆, thus a change during a small period of time , (), is equal to ϵ√.Then, Equation S15 can be modified as below: As a result, the temperature variation induced by voltage noise can be included in the noise term, and the stochastic differential equation of temperature considering collective noises of temperature and voltage is described in the following equation.
In the same way of calculating the mean (Equation S12) and variance (Equation S13) from Equation S2, the long-term mean and variance of Equation S17 can be calculated as below: which are Equations 7 and 8 in main text.The correspondence between this derived general solution and the numerical model (p-osc model) is illustrated as below.

In-depth description of the MBM to solving MVC problem
As described in the manuscript, in the operation of the Boltzmann machine, the negative derivative of the Hamiltonian (ℎ()) for each pbit at the current time step (   = −ℎ(  )/  ) is used as the input value for that pbit at the next time step.However, as can be seen from the ℎ() in Equation 10of main text,    depends on the absolute size of the ℎ() which is affected by α and β (or additional parameters in other problems 5 ) or the way the developed ℎ() is defined.Therefore, when simply using    as an input value for the next step, an excessively large or small    may make the stochastic behavior too deterministic, or may not induce significant differences in others.Thus, to mitigate this problem previously, we arithmetically calculated the maximum and minimum values of    that could be obtained from the initially set ℎ() and normalized    into the device operating range  ext .
In case of minimum vertex covering (MVC) problem 5   Fig. S9a shows the minimal vertex coverings (MVCs) of the G(6, 7) graph presented in the main text.Given the relatively small scale of the problem, it is feasible to identify the correct ( = '011001' and '101001') and incorrect solutions using a brute-force algorithm.
The energy maps according to the Hamiltonian (Equation 10) is shown in Figs.S9b-c.This is implemented in the MBM, employing the Ising model approach 5 .The energy map reveals that the MBM system assigns global minimum energy states to  = '011001' and '101001', the same as the ground truth from the brute-force algorithm.There are also local minimum states, such as '010111', '011011', '011101', '101011', '101101', '111110', '111001', and '111010', which are energetically favorable compared to their immediate surroundings but hold higher energy than the global minimum states.In accordance with the Boltzmann law, global minimum states appear most frequently, followed by the local minimum states.Thus, we anticipated that Ȳ incor , which represents the average occurrence of the five most frequently appearing non-correct cases, would primarily consist of these local minimum states.Figs.S10b-d validate this assumption by displaying the MBM simulation results for every configuration  in both the presence and absence of autocorrelation.
Autocorrelation-induced non-independent probabilistic behavior hinders the MBM from effectively exploring the solution space represented by the energy map, leading to an increase in occurrence of local minimum states.Consequently, we observe this effect as a rise in Ȳ incor in Figs.S10c-d as well as Figs.4c-d in the main text.More specifically, as shown in Fig. S10c, without the inclusion of a delay period, the '011101' configuration becomes most prominent in the presence of autocorrelation.This highlights the risk that autocorrelation could cause a local minimum state to appear dominantly within a small number of iterations, thereby deviating from the ground truth.Fig. S11 shows the probabilistic oscillations of an NbO x memristor-based oscillator fabricated under different conditions compared to the device described in the main text.For this device, the dimensions are 5 x 5 μm 2 , featuring a Pt bottom electrode and a Ti top electrode.A 40 nm thick NbO x layer was deposited using the reactive sputtering method, maintaining an Ar:O 2 ratio of 13:7.

Fig. S3 .
Fig. S3.NbO x oscillator measurement system.a A schematic of measurement system with Keithley 4200 SCS semiconductor analyzer and dark box.b Equivalent circuit diagram of a

5 . 6 .
Fig. S4. a Circuit diagram of NbO x oscillator and measured applied external voltage (V ext ) to time in case of V ext = 1.46 V. b V ext vs time plot magnifying the region from 5 us to 30 us (gray-colored region in a).c Histogram of V ext .d Power spectral density (PSD) plot of V ext to frequency, showing the relationship between PSD and frequency, indicating that the closer a frequency is to zero, the closer it is to white noise.

Fig
Fig. S7. a T-V plot of pseudo-equilibrium states of numerical oscillator model at V ext = 1.47 V, R L = 1.2 kΩ (blue dot is an equilibrium point).b Temperature histogram (red) from T-V plot and normal distribution (blue) from compact stochastic oscillator model (Equations 7 and 8) under the same condition.

11.
Fig. S8.Autocorrelation function of all V ext conditions of experimental NbO x oscillator.

Fig. S10. a
Fig. S10. a Energy map of MVC problem for given graph G(6, 7).MBM simulation results for every  states at 1,000 iterations in case of b un-autocorrelated, c autocorrelated, no delay period (Δ), d autocorrelated, Δ = 2.

13.
Fig. S11.Probabilistic oscillation in Ti/NbO x /Pt volatile memristor a Threshold switching I-V curve of Ti/NbO x /Pt volatile memristor (black) and the load lines of 2.0 kΩ series resistor under various V ext .b Probabilistic current oscillations by ramping the V ext from 1.50 V to 1.62 V at 0.02 V interval.c Selfoscillator circuit consists of a memristor, series resistor, and internal capacitive volume of the NbO x layer.d Optical microscope image of memristor used in this experiment and e schematic of Ti/NbO x /Pt structure with each layer's thickness.f Oscillation probability distribution to the V ext showing sigmoidal decreases to the V ext .

Table S2 . Parameters used in NbO x oscillator model.
1 ,  1 ) =  1 ( 0 +  pert + From the simulation results of Supplementary Information section 6, we conclude that oscillation is generated due to the drastic acceleration at T MIT .Thus the temperature threshold (T MIT ) is assumed as a boundary condition for oscillation generation.By solving the simultaneous equations of Equations S1 and S2, we derived the temperature variation due to collective interference by temperature and voltage noises.This approximated calculation is based on linear dependence and assumes stationarity and no feedback.First, we derive the general solution of Equation S2, the stochastic process of voltage, as follows.
Next, we calculate the temperature distribution induced by collective interference of temperature and voltage by inserting the result of voltage variation into the voltage term in the OU process equation of temperature.