Figure 1b shows an optical microscope plane-view image (top) and a cross-sectional transmission electron microscope (TEM) image of an OTS device, which has a pore-type structure with a pore size (*d*) of 300 nm. The fabrication process of the OTS device is described in detail in the Experimental Section. Figure 1c shows the characteristic current-voltage (*I*-*V*) curves of the device, which is nominally the same as Fig. 1b. In ten repeated measurements, the OTS device shows the variation in *V*th within 2.5 ~ 3.0 V while *V*H is around 1 V with a negligible variation. In addition, *R*off and *R*on of OTS devices have been read around ~ 108 and ~ 102 Ω, respectively. Figure 1d shows *V*out waveforms of a typical OTSNO with varying *V*G, where the OTSNO consists of an OTS device and a FET (n-MOSFET, LND150N3-G, Microchip Inc.). It is clearly shown that the natural frequency (*f*nat) of the OTSNO increases with *V*G as expected. To quantify *f*nat at each *V*G, the fast Fourier transform (FFT) is performed as shown in Fig. 1e and *f*nat, defined as the frequency at the primary peak, is plotted as a function of *V*G in Fig. 1f. It is observed that *f*nat is well defined and linearly proportional to *V*G being adjustable in the range of 0.5 ~ 2 MHz. These results clearly show that the oscillation in the OTSNO is composed of fundamental-frequency components and their harmonics while keeping other components negligible.

Next, we investigate the synchronization behavior of coupled OTSNOs (see Fig. 2a). Figure 2b shows the output waveforms of two OTSNOs coupled by a capacitor. It is observed that the output waveforms of two oscillators look nearly the same with a certain degree of phase difference, demonstrating that two oscillators are synchronized to each other (for a detailed study of the synchronization behavior, see section S1 in Supplementary Information). In Fig. 2c, *V*out2 is plotted as a function of *V*out1, which presents the so-called “phase portrait” of a coupled-oscillator system. It shows a butterfly-shaped attractor curve, implying the AP relationship between the two oscillators. The attractor curve appears as a band implying the chaotic nature of the OTSNO [21, 23, 27, 31], which is desirable for the application to the IM because it helps the system to find the solution of the network configuration with the global minimum energy [32].

To examine the phase stability of two coupled oscillators, we have performed a simulation study where the circuit parameters of the system are systematically varied. Especially, we have investigated two types of coupling - capacitive and resistive coupling - with varying their values, *C*C and *R*C, respectively. In the SPICE (Simulation Program with Integrated Circuit Emphasis) simulation, the OTS device is modeled as a voltage-controlled switch characterized by four parameters (*V*th, *V*H, *R*on, *R*off) in parallel connection with a parasitic capacitor. Figure 2d and 2e show the phase difference (*Δφ*) between oscillators #1 and #2 as a function of (*C*C, *R*L) and (*R*C, *R*L), respectively. Here, the parameters of both OTS devices are set to (*V*th, *V*H, *R*on, *R*off)=(3.3 V, 0.7 V, 150 Ω, 10 MΩ). *Δφ* is calculated by 2*π*(*T*2-*T*1)**f*sync after ~ 200 cycles, where *T*2-*T*1 and *f*sync are the time difference between adjacent peaks of oscillator #1 and #2 and the oscillation frequency in the synchronized state, respectively. It is found that the AP relationship between two oscillators is more stable in the capacitive coupling with *C*C for AP relationship spanning almost two orders (2 ~ 200 pF) independent of *R*L while the range of *R*C for AP relationship is relatively narrow with a strong dependence on *R*L.

For the practical application, an important aspect is the tolerance of the coupled oscillator system to the inevitable variation of the OTS device. Figure 2f and 2g show *Δφ* as a function of *ΔV*th (= *V*th2-*V*th1) and *ΔV*H (= *V*H2-*V*H1), where we have varied *V*th and *V*H of OTS #2 systematically while keeping those of OTS #1 fixed. It is shown that, in the case of the capacitive coupling, the AP relationship is robust against the *V*th-variation of up to ± 20%, which is significantly higher compared to resistive coupling where the limit is only ± 5%. The superior stability of the AP relationship in the *C*C-coupled oscillator system can be explained qualitatively by the ability of the coupling capacitor to store energy. In detail, if a voltage difference is generated between two oscillators, then it can be stored in the coupling capacitor. The stored energy is returned to the oscillators making them repel each other, stabilizing the AP relationship. In contrast, since the resistor can’t store energy in the *R*C-coupled oscillator system, the voltage difference between oscillators generates heat dissipation in *R*C. This makes the oscillators lose their energy destabilizing the AP relationship in the *R*C-coupled oscillator system.

Another intriguing finding is that the phase relationship of the coupled OTSNOs in both cases is hardly sensitive to the variation of *V*H compared to the variation of *V*th. It is related to the asymmetry of the output waveform (see Fig. 1d), which shows a much longer falling period compared to the rising period. Recalling the mechanism of the oscillating behavior of the OTSNO, it is easily found that the peak of the oscillation aligns with the ON-to-OFF transition of the OTS device whereas the valley point corresponds to the OFF-to-ON transition. Consequently, when the OTS reaches its peak, it transitions into a highly resistive state with *R*off ~ 10 MΩ. As mentioned above, since the OTS can be described by a parallel connection of a voltage-controlled switch and a parasitic capacitor, such a high *R*off results in a large RC delay, and consequently, the elongation of the falling period compared to the rising period. During the falling period, the parasitic capacitor is charged and the voltage across it runs into *V*th starting from *V*H. Since the charging time is mainly determined by *V*th, the variation in *V*H has little effect on the phase relationship of the coupled OTSNOs. The output waveforms and the phase portraits at points in various regions in Fig. 2d and 2e are presented in the Supplementary Information (Fig. S2).

Based on these results, we have built an IM using OTSNOs and coupling capacitors (*C*c=100 pF) as shown in Fig. 3a and, as a benchmark test, have tried to solve a MaxCut problem with a Mobius ladder geometry composed of 14 nodes and cubic connections (see Fig. 3b). We have also tried another geometry, which is presented in Fig. S3 in the Supplementary Information. Tlhe details of the measurement is described in the Experimental Section. A representative output waveforms of the 14 OTSNOs are presented in Fig. 3c. We have employed both the second-harmonic injection locking (SHIL) and simulated annealing (SA) techniques [16, 33–35] to improve the reliability of the solution by keeping the system from being stuck to a local minimum. In this technique, all the oscillators are driven by a DC voltage (*V*dc) and a modulational AC voltage (*V*SHsin(*ω*SH*t*)), which locks the oscillation frequency to *ω*SH (= 2*ω*sync ~ 2.75 MHz, where *ω*sync is the oscillation frequency in the synchronized state). In the right panel, which is the expansion of the waveforms in the range of 40 ~ 43 µs, it is shown that the odd-numbered OTSNOs oscillate with the nearly same phase while the even-numbered ones also show the synchronized behavior with a clearly distinguished phase.

To quantify the phase difference between oscillators, we have calculated the amplitude of the in-phase component of each *V*out(*t*) by applying the low-pass filter (the cut-off frequency of 0.1 MHz) to the product of a reference wave (cos(*ω*SH*t*)) and *V*out(*t*). Such obtained in-phase amplitude (Ampin) of the oscillators are plotted as a function of time in Fig. 3d. It is clearly observed that the oscillators are separated into two groups with the opposite signs of Ampin implying the AP relationship between those two groups of oscillators. In addition, it shows that those two groups are clearly separated 35 µs after applying the bias. Considering that *V*SH gradually increases over 35 µs in the simulated annealing scheme, it indicates that the time-to-solution (*T*sol) is expected to be at most 35 µs. From repeated implementations of the same experiment, we have observed *T*sol is in the range of 35 ~ 40 µs. It is much shorter than *T*sol ~ 4 ms reported in a previous work on the hardware implementation of a PTNO-based Ising solver [16] for solving an 8-node MaxCut problem although, in a simulation study, *T*sol=30 µs was presented for solving a 100-node MaxCut problem with random cubic connections.

The energy-to-solution (*E*sol), as another performance metric of the IM, is calculated by integrating the instantaneous power consumption (*P*(*t*) = *V*(*t*)×*I*(*t*), where \(V(t)={V_{dc}}+{V_{SH}}\sin (\omega t)\) applied to the system (refer to Fig. S3a) and \(I(t)=\sum\limits_{{i=1}}^{{14}} {{I_{i,out}}(t)}\), the sum of the output currents (\({I_{i,out}}(t)\)) at each node) over *T*sol. *P*(*t*) is plotted as a function of time in Fig. 4a, leading to an estimation of *E*sol around 2.3 µJ. It can not be compared with that of the PTNO-based IM as it was not known in the previous work. Instead, we have tried to compare the energy consumption in a single OTSNO and PTNO because of the negligible energy consumption in the coupling capacitor (see Fig. S4 in the Supplementary Information) as expected in the pure capacitive circuit. The energy consumption per cycle (*E*cycle) is estimated to be ~ 1.25 nJ/cycle (see Fig. 5b) for a 500 nm (*d*: diameter of the pore)-sized OTSNO and ~ 0.9 nJ/cycle for a 200 nm (the length of VO2 channel)-sized PTNO, respectively. In the OTS device, *E*cycle is expected to be scaled with the dimension of the device because of the shrinkage of the switching volume. Therefore, we have investigated the dependence of *E*cycle of the OTSNO on the diameter of the pore as presented in Fig. 4c (for the waveforms of *V*(*t*), *I*(*t*), and *P*(*t*) corresponding to OTSNOs with various pore sizes, see Fig. S5 in Supplementary Information). For comparison, *E*cycle of the PTNO [16, 36] is also located in the same graph, clearly showing that the OTSNO consumes less energy by about an order than the PTNO with similar device dimensions. The difference is attributed to the difference in the switching mechanisms of the PTNO and the OTSNO. In detail, the PTNO is based on the phase transition of the Mott insulator, which requires heating of a part in the channel material up to the transition temperature. In contrast, the OTSNO is believed to be based on the filling and evacuation of the trap states in the OTS, the electronic process in nature, although there is still controversy about its switching mechanism [21, 23, 31, 37]. In addition, Fig. 4c shows that *E*cycle of the OTSNO scales as *E*cycle ~ *d*1.81 with the diameter of the pore, being close to the expected relation of *E*cycle ~ *d*2 for the case of the uniform current density in the switching material. Therefore, it seems to imply that the reduction in *E*cycle is attributed to the shrinkage of the switching volume inside the cylinder. As a result, *E*cycle of the OTSNO extrapolates to ~ 1 pJ/cycle at *d* = 10 nm, which enables the development of a highly efficient large-scale IM.

Finally, we have performed a simulation study to investigate the scalability of the OTSNO-based IM. As a typical benchmark task, we have studied the Mobius ladder geometry composed of variable numbers of OTSNOs and connections (see Fig. 5a). Figures 5b ~ 5d show the representative examples of Ampin of oscillators with varying the number of nodes (*N*) from 22 to 102. The correctness of the obtained solution has been verified by comparing the energy of the solution with that of the ground state. From 20 repetitions of the simulation with varying the initial phases of each oscillator randomly, we have obtained the performance metrics, the success probability (*P*success) and *T*sol, as shown in Fig. 5e. Note that, as *N* increases, *T*sol is linearly proportional to log(*N*) with a slight decrease in *P*success. We have repeated similar simulations with varying the number of connections (*N*c) in the way of adding connections to the next nearest neighbors in order with a diagonal connection being fixed. Representative examples of Ampin of oscillators for various geometries are presented in the Supplementary Information (see Fig. S6). In Fig. 5f, it is observed that *T*sol shows a sublinear dependence on *N*c while *P*success decreases slightly with increasing *N*c. It implies that, as *N*c increases, local minima are likely to be formed in the energy landscape and the oscillators are trapped in those local minima. Due to the intrinsic chaotic nature of the OTSNO device as mentioned in the discussion of Fig. 2c, we expect that the problem of such local minima might be alleviated in the hardware implementation of the IM using OTSNOs.