How fast can vanadium dioxide neuron-mimicking devices oscillate? Physical mechanisms limiting the frequency of vanadium dioxide oscillators

Memristors are emerging electronic devices with, usually, a simple and compact architecture made of an active layer and two electrodes. Their channel conductance can be programmed into two or more conductance states by triggering with an external stimulus. This resistive switching (RS) makes these devices ideal candidates to fabricate the elemental bricks of neuromorphic circuits, neurons or synapses, respectively, depending on whether RS is volatile or not. Recently, among the many possible physical mechanisms that give rise to RS, the insulator-to-metal transition (IMT) has attracted huge interest. IMT occurs in transition-metal oxides (TMOs), such as VO₂ or NbO₂. TMOs should behave like metals since they possess partially filled (3d, 4d or 4f) orbitals that are expected to give rise to a conductive band of electronic states. However, these materials are insulators since the incomplete band is indeed split into an empty and a filled band. Most importantly, the physical properties of TMOs may vary abruptly under the influence of external excitation. In this respect, VO₂ exhibits a steep RS due to the change from insulating monoclinic (M₁) to metallic tetragonal rutile-like (R) crystal structure. The IMT M₁ $\longrightarrow$ R is driven by temperature and occurs at a threshold temperature $T_{\mathrm{IMT}}$ of about 340 K. The reverse metal-to-insulator transition (MIT) R $\longrightarrow$ M₁ occurs during cooling down of the material, when a threshold temperature $T_{\mathrm{MIT}}$ is achieved. Generally, $T_{\mathrm{MIT}}$ is smaller than $T_{\mathrm{IMT}}$. Thus, the experimental curve of resistance R against temperature T shows hysteresis.

RS is also induced in VO₂ when the latter is the channel of a two-terminal device, by applying voltage or injecting current across the channel. In this case, experimental evidence supports the idea that the driving mechanism is self-heating [2–5]. The channel of VO₂ devices switches from an high resistance state (HRS) to a low resistance state (LRS) depending on the VO₂ temperature, which is modulated by the Joule effect. The RS is volatile since the resetting to the original HRS occurs spontaneously upon a decrease in electrical power, which induces cooling down of VO₂. Consequently, VO₂ devices are named volatile Mott memristors [6–9].

A relaxation oscillator, which generates a non-sinusoidal periodic waveform at its output, can be obtained by connecting in series a VO₂ volatile memristor with an RC parallel circuit [7, 10] (figure 1(e)). The self-oscillatory behavior at the output node of the oscillator circuit arises from cyclic charging and discharging of the external capacitor $C_{\mathrm{ext}}$, where these cycles are self-piloted by the VO₂ volatile memristor being periodically in the ON/OFF (LRS/HRS) state. We define ($V_{\mathrm{IMT}}$, $I_{\mathrm{IMT}}$) to be the IMT point of the curve current (I) versus voltage across the VO₂ device (V_D ), which is the point of the curve where the IMT occurs. Similarly, ($V_{\mathrm{MIT}}$, $I_{\mathrm{MIT}}$) is the MIT point of the I–V_D curve where the MIT occurs. Then, the condition for the self-oscillatory regime is that the fixed oscillator circuit parameters, applied bias $V_{\mathrm{DD}}$ and external resistance $R_{\mathrm{ext}}$, are such that the load-line current ($I_\mathrm{L} = (V_{\mathrm{DD}} - V_{\mathrm{D}}) / R_{\mathrm{ext}}$) is larger than the $I_{\mathrm{IMT}}$ at $V_{\mathrm{IMT}}$, and smaller than $I_{\mathrm{MIT}}$ at $V_{\mathrm{MIT}}$. It is important to observe that VO₂ oscillators are quite compact, which is an advantage for the scaling up of circuits with oscillators, considering that CMOS-based oscillators are realized with up to tens of transistors.

Zoom In Zoom Out Reset image size

Recently, oscillatory neural networks (ONNs) based on VO₂ oscillators have been under intense scrutiny [10–12] as neuromorphic computing engines. Neuromorphic computing is an emerging discipline that aims to mimic highly efficient, brain-inspired computing strategies. There are different approaches to imitating brain activity. The mainstream implementation of neuromorphic circuits is through artificial neural networks [13]. In these systems, artificial neurons accumulate weighted input values, and the output is generated from the inputs through the application of an activation function. In spiking neural networks [14] the information is transmitted at the end of each cycle only at a given threshold. This closely imitates natural neurons, which fire when their membrane potential reaches a threshold, sending a signal to neighboring neurons. A third approach is represented by ONNs [15], which are systems of coupled oscillators that have been demonstrated to be powerful computing engines to emulate some brain functionalities [16–20].

ONNs encode information into the stable phase difference between each oscillator of the network and a reference oscillator, achieved if/when the system settles into a synchronized state [21, 22]. The energy that the ONN takes to compute is proportional to the ONN settling time [23], which is the number of cycles to synchronize multiplied by the oscillation period. The increase in oscillator frequency should bring about an overall reduction in energy consumption for ONN operativity. Thus, it is of interest for ONN applications to explore the range of frequencies at which oscillators can operate. Specifically, given the electrothermal mechanism driving the RS of VO₂ volatile memristors, it can be expected that not only the electrical time constant, but also the thermal time constant plays a role in determining the oscillation period.

In this study, we perform electrothermal 3D technology computer-aided design (TCAD) mixed-mode simulations of a VO₂ oscillator. We find that there exists an upper limit to its operating frequency, where this limit is not predicted from a purely circuital model of the VO₂ oscillator. We investigate the intrinsic physical mechanisms that give rise to this upper limit. Our TCAD simulations show that the IMT/MIT points of current versus voltage across the VO₂ device starts to depend on frequency, during oscillation, below some threshold values of $C_{\mathrm{ext}}$. This implies that the condition for the self-oscillatory regime may be satisfied by a given load-line in the low-frequency range but no longer at higher frequencies, with consequent suppression of oscillations. We find that this variation of the IMT/MIT points is due to a combination of different factors: 1) intermediate resistive states of VO₂ achieved below some threshold values of $C_{\mathrm{ext}}$ that make HRS conductance of VO₂ channel variable with frequency and 2) interplay between frequency and heat transfer rate. Although the upper limit on the frequency that we extract is linked to the specific VO₂ device we simulate, our findings apply qualitatively to any VO₂ oscillator. Overall, this study elucidates the link between electrical and thermal behavior in VO₂ devices that sets a constraint on the upper values of the frequency of a VO₂ oscillator.

We use finite element analysis (FEA) to solve the partial differential equations (PDEs) that model the device physics. In particular, we use the commercial Silvaco TCAD software suite to perform FEA calculations. This facilitates the implementation of multi-physics simulations, with many easily accessible built-in models of physical phenomena that can be integrated into a single simulation flow. We simulate VO₂ volatile memristors with so-called crossbar (CB) geometry (figure 1(d)). These devices have been recently investigated since they show enhanced reliability and repeatability [11]. In FEA, the simulation domain (device architecture) has to be split into a discrete number of elements, whose ensemble is called a mesh, where the PDEs are numerically solved. In the present case, the VO₂ device configuration is 3D and cannot be straightforwardly reduced to a 2D or 1D geometry. In addition, it is well known that 1D or 2D approximation of heat transfer analysis for a 3D problem may introduce errors (see for instance Kamala and Goldak [24]). Thus, we build a 3D mesh using Silvaco Victory Mesh [25] to reproduce the CB geometry. In this way, the coupled electrothermal behavior inside the VO₂ channel is mimicked with a good approximation, which plays a key role in modulating the behavior of VO₂ volatile memristors and oscillators [26–28].

We use Silvaco Victory Device [29] to solve PDEs of VO₂ device physics on the generated mesh. We consider VO₂ as a conductor, and we use a customized C-Interpreter function [26–28] of the PCM model [29] to empirically emulate VO₂ resistivity. Overall, VO₂ resistivity depends on local temperature T and time t after the Johnson–Mehl–Avrami model [30–32]:

$$\begin{equation} \rho (t + \Delta t, T) = \rho (t, T) + [\hat{\rho} (T) - \rho (t, T)] \times (1 - \mathrm e^{- \Delta t / \tau}) ,\end{equation} \tag{ 1 }$$

where t is time, $\Delta t$ is the time step, T is the local temperature, $\hat{\rho}$ is the steady-state resistivity and τ is a lifetime. We calibrate the steady-state resistivity $\hat{\rho}$ versus T (figure 1(a)) using experimental data of resistance $R_{\mathrm{exp}}$ versus temperature T [26]. Further details on the VO₂ resistivity model and the calibration procedure are contained in appendix B.

Since VO₂ is described as a conductor, the current density is $J = E / \rho = - \nabla \phi / \rho$, for E the electric field and φ the electric potential. Thus, to calculate self-consistently the current transport and temperature distribution inside the VO₂ volatile memristor we need to solve the Poisson equation coupled to the heat flow equation:

$$\begin{equation} C \frac{\partial T}{\partial t} = \nabla (K\, \nabla T) + \frac{(\nabla \phi)^{2}}{\rho} ,\end{equation} \tag{ 2 }$$

where K = 0.06 W (cm K)⁻¹ [33] is the thermal conductivity, and C = 3 J (cm³ K)⁻¹ [34] is the heat capacity per unit volume of VO₂, respectively. The term $(\nabla \phi)^{2} / \rho$ corresponds to the heat generation rate $\rho J^{2}$ due to the Joule effect. For simplicity, we consider the bottom surface of the VO₂ device as the only thermal contact (boundary). Therein, we specify the thermal boundary condition for the solution of (2) by setting that the projection of the energy flux, onto the unit normal, external to the boundary surface is $(T - T_{0}) / (R_{\mathrm{Th}} \mathscr{A}$), for T₀ the external temperature, $R_{\mathrm{Th}}$ the interface thermal resistance and $\mathscr{A}$ the bottom surface area of the device. The term $1 / (R_{\mathrm{Th}} \mathscr{A}$) regulates the heat dissipation of the VO₂ volatile memristor with respect to the environment. Finally, we perform simulations of the VO₂ relaxation oscillator by combining the TCAD engine to solve the PDEs of the VO₂ device with a circuit solver (the so-called TCAD mixed mode), in such a way that the device physics is self-consistently solved at each iteration of the solver of the circuit equations.

We perform electrothermal 3D TCAD mixed-mode simulations of a VO₂ oscillator (circuit scheme in figure 1(e)), where the VO₂ CB device is made of an 80 nm thick VO₂ layer inserted between the top and bottom Pt electrodes. The electrodes are 250 nm wide, while the VO₂ layer is a square with each side 1 µm long. We set $T_{0} = 303$ K and $R_\mathrm{Th} = 8.3 \times 10^{5}$ K W⁻¹ as boundary conditions at the thermal contact. We choose the circuit load-line specified by V$_\mathrm{DD} =$ 2 V, R$_\mathrm{ext} =$ 15 kΩ in order to comply with the self-oscillatory regime. We modulate the operational oscillator frequency by varying $C_{\mathrm{ext}}$ in the interval 5 nF–400 pF. The simulated voltages at the output node $V_{\mathrm{out}}$ are shown in figure 2(a). The oscillations are turned-off at $C_{\mathrm{ext}} =$ 400 pF. Frequency versus $C_{\mathrm{ext}}$ of TCAD simulated oscillations are plotted (blue diamonds) in figure 2(b) and reported in table D3 in appendix D. Theoretical frequencies (red asterisks) that are also shown in the same figure, have been calculated through the purely circuital model of the VO₂ oscillator where the VO₂ volatile memristor is approximated as a two-state (LRS/HRS) variable resistor, activated by threshold voltages $V_\textrm{IMT}$ and $V_\textrm{MIT}$, respectively (see appendix D, and [35]). It can be seen that for larger values of $C_{\mathrm{ext}}$ the simulated and the theoretically calculated data are almost overlapping. However, the theoretical calculations predict frequencies corresponding to $C_{\mathrm{ext}}$ values equal to or below 400 pF, while we are not able to obtain oscillations in this range of $C_{\mathrm{ext}}$ values. Thus, an upper limit to the achievable frequencies is not explained within the framework of a theoretical, purely circuital model of the VO₂ oscillator (appendix D, [35]).

Zoom In Zoom Out Reset image size

In figure 3(a), we plot curves I versus $V_{\mathrm{D}} = V_{\mathrm{DD}} - V_{\mathrm{out}}$ (symbols) as extracted from simulated oscillations for $C_{\mathrm{ext}} =$ 10 nF, 5 nF, 1 nF, 800 pF, 500 pF and 400 pF. The circuit load-line is also drawn (dashed black line). It can be seen that up to a certain frequency ($C_{\mathrm{ext}} =$ 5 nF) all the curves essentially overlap, and then some mismatch becomes evident at higher-frequency values (lower $C_{\mathrm{ext}}$ values). In particular,

• 1.  
  the IMT/MIT points become dependent from $C_{\mathrm{ext}}$ below some threshold value ($C_{\mathrm{ext}} \leqslant$ 800 pF for IMT points, $C_{\mathrm{ext}} \leqslant$ 1 nF for MIT points); see insets of figure 3(a),
• 2.  
  the HRS branch of the I–$V_{\mathrm{D}} = V_{\mathrm{DD}} - V_{\mathrm{out}}$ curve is also dependent from $C_{\mathrm{ext}}$ for $C_{\mathrm{ext}} \leqslant$ 1 nF.

Zoom In Zoom Out Reset image size

We compare the above curves with the quasi-static I–$V_{\mathrm{D}}$ (azure thick line of figure 3(a)). The curve is obtained by electrothermal 3D TCAD simulation of a voltage-controlled VO₂ volatile memristor, where the VO₂ device is connected in series with an external resistance of 1 kΩ, and the voltage sweeping rate is such to guarantee quasi-static conditions (see appendix C). The positions of the IMT/MIT points as well as the HRS branches of I–$V_{\mathrm{D}}$ are the same in quasi-static conditions and during oscillations up to when a dependence on $C_{\mathrm{ext}}$ arises of I–$V_{\mathrm{D}} = V_{\mathrm{DD}} - V_{\mathrm{out}}$ curves. This dependence on $C_{\mathrm{ext}}$ of IMT/MIT points implies an apparent modulation of the negative differential resistance (NDR) region with frequency. Indeed, the NDR region of a volatile memristor is rigorously defined from its DC current–voltage characteristics [36]. In our case, since the quasi-static I-$V_{\mathrm{D}}$ and the DC I-$V_{\mathrm{D}}$ overlap (as demonstrated in figure C2(b) in appendix C), the NDR of the quasi-static I-$V_{\mathrm{D}}$ in figure 3(a) is a rigorous reference. It should be observed that the value of current at the IMT point, $I_{\mathrm{IMT}}$, increases with frequency for $C_{\mathrm{ext}} \leqslant$ 800 pF, while the value of current at the MIT point, $I_{\mathrm{MIT}}$, decreases with frequency for $C_{\mathrm{ext}} \leqslant$ 1 nF. Since the load-line current is independent of $C_{\mathrm{ext}}$, then it may be expected that it will become smaller than the $I_{\mathrm{IMT}}$ at $V_{\mathrm{IMT}}$, and/or larger than the $I_{\mathrm{MIT}}$ at $V_{\mathrm{MIT}}$ for a certain value of $C_{\mathrm{ext}}$. This will make the load-line no longer compliant with the conditions for a self-oscillatory regime, which explains our observation of inhibition of self-oscillations at certain frequencies.

To understand better the upper limit to the frequency of a VO₂ oscillator, we need to investigate the mechanism of variation of the IMT/MIT points, and of the HRS branch of I–$V_{\mathrm{D}} = V_{\mathrm{DD}} - V_{\mathrm{out}}$, above threshold frequencies. From the electrothermal 3D TCAD mixed-mode simulated oscillations, we extract the local temperature as probed in the geometrical center of the device, T. We consider it as a good descriptor of the temperature of the 'channel', given that the temperature is approximately uniform therein [26, 28]. We plot in figure 3(b) the curves (symbols) of T vs $V_{\mathrm{D}} = V_{\mathrm{DD}} - V_{\mathrm{out}}$ for $C_{\mathrm{ext}} =$ 10 nF, 5 nF, 1 nF, 800 pF, 500 pF and 400 pF. We also extract from the electrothermal 3D TCAD simulations of quasi-static I–$V_{\mathrm{D}}$ the associated T–$V_{\mathrm{D}}$ (azure thick line in figure 3(b)). It can be seen that:

• 3.  
  the temperatures $T_{\mathrm{IMT}}$, associated with IMT points, become dependent on $C_{\mathrm{ext}}$ for $C_{\mathrm{ext}} \leqslant$ 800 pF, while the temperatures $T_{\mathrm{MIT}}$, associated with MIT points, are constant and equal to 320 K.

In particular, the points ($V_{\mathrm{IMT}}$, $T_{\mathrm{IMT}}$) appear (approximately) located along the part of quasi-static T–V$_{\mathrm{D}}$ corresponding to HRS-to-LRS transition (dashed red line). In the case of $C_{\mathrm{ext}} =$ 400 pF, T is constant, with a value slightly below the HRS-to-LRS branch of quasi-static T–V$_{\mathrm{D}}$. Thus, T does not achieve the value for RS at $C_{\mathrm{ext}} =$ 400 pF, which is consistent with the suppression of oscillations in this case.

In figure 4(a), it is shown the plot of T against time normalized to the oscillation period ($t / period_{\mathrm{osc}}$) for an oscillation cycle, where $C_{\mathrm{ext}} =$ 10 nF, 5 nF, 1 nF, 800 pF, 500 pF and 400 pF. The use of $t / period_{\mathrm{osc}}$ allows to compare the temperature profiles during oscillations at different frequencies. In figure 4(b), the corresponding local resistivity as probed in the geometrical center of the device (ρ) is plotted against $t / period_{\mathrm{osc}}$ for an oscillation cycle. The T versus $t / period_{\mathrm{osc}}$ (and ρ versus $t / period_{\mathrm{osc}}$) curves show dependency on frequency for $C_{\mathrm{ext}} \leqslant$ 1 nF, even though the qualitative behavior is the same. Finally, figure 4(c) reports the plot of ρ versus T during an oscillation cycle.

Zoom In Zoom Out Reset image size

Figure 4(a) shows that the temperature increases abruptly to the maximum $T_{\mathrm{peak}}$ value very shortly after the IMT point, which is consistent with a sharp increase in the current due to IMT. This value is well above the temperature needed to complete the transition to the metallic state, thus ρ becomes equal to $\rho_{\mathrm{LRS}}$ (which, for simplicity, we have assumed to be constant, see table B2 in appendix B). Then the external capacitor starts to charge and the temperature to decrease. It is important to observe that the cooling in VO₂ is not related to any removal of electrical stimulus $V_{\mathrm{DD}}$ that in contrast is kept constant. Instead, it is due to the change (decrease) in electrical power (the heat source) after IMT switching.

The resetting to the HRS state occurs at the MIT point, and starts the discharging of the external capacitor. Overall, figure 4(c) shows that ρ at the beginning of discharging has not recovered the frequency-independent values of the HRS heating branch (figure B1(b) in appendix B). Instead, it assumes some constant, slightly frequency-dependent values (red crosses in figures 4(b) and (c)). However, up to $C_{\mathrm{ext}} =$ 1 nF, it always recovers the HRS heating branch during the heating cycle (gray circles in figures 4(b) and (c)), before the HRS-to-LRS transition. This explains why IMT points do not depend on $C_{\mathrm{ext}}$ up to $C_{\mathrm{ext}} =$ 1 nF. For $C_{\mathrm{ext}} =$ 1 nF, ρ at the beginning of discharging has the same value of $\rho_{\mathrm{IMT}} = \rho (T_{\mathrm{IMT}})$ of quasi-static I–$V_{\mathrm{D}}$. This value is retained during the entire heating cycle and the HRS heating branch is intercepted just at the IMT point. This explains why the IMT points are the same, although the HRS part of I–$V_{\mathrm{D}} = V_{\mathrm{DD}} - V_{\mathrm{out}}$ differs from the one of the quasi-static I–$V_{\mathrm{D}}$. Finally, for $C_{\mathrm{ext}} \leqslant$ 800 pF, the value of ρ at the beginning of discharging is smaller than $\rho_{\mathrm{IMT}}$. Thus, the HRS heating branch is never recovered, and not both the HRS part of I–$V_{\mathrm{D}} = V_{\mathrm{DD}} - V_{\mathrm{out}}$ and also the IMT point start to depend on $C_{\mathrm{ext}}$.

The above-described behavior is clarified by observing the ρ versus T during oscillations in figure 4(c). The ρ values at the beginning of discharging (red crosses in figures 4(b) and (c)) correspond to the ones where the temperature trend switches from decreasing to increasing (sign flipping of $\mathrm dT/\mathrm dt$; red crosses in figures 4(a) and (c)). Since this occurs before the transition to the insulating state is completed, the TCAD ρ value is frozen at an intermediate resistivity value until when/if the HRS heating branch of ρ (figure B1(b) in appendix B) is intercepted (figure 1(c)). In this respect, our TCAD model follows the experimental finding of intermediate resistance states attainable in resistance versus temperature curves when the sweep in temperature is reversed before the full transition is achieved [1].

The above discussion demonstrates the key role played by the change of T versus t with frequency. Thus, we need to look at the variation of $\mathrm dT/\mathrm dt$ with frequency. In principle, the temporal derivative of T can be calculated directly from T versus t of simulated oscillations. However, we can also approximate it from the heat diffusion equation [37] as the sum of two components:

$$\begin{equation} C_\mathrm {th} \frac{\mathrm d T}{\mathrm d t} = P - (T - T_{0})(A \overline{h}) ,\end{equation} \tag{ 3 }$$

where P is the electrical power from simulated oscillations, $C_\mathrm {th} = C V_{\mathrm{CB}}$ is the thermal capacity where $V_{\mathrm{CB}}$ is the volume of the CB device, A is the heat transfer surface area and $\overline{h}$ is the average over A of the heat transfer coefficient [37]. Through equation (3), the t-behavior of T can be simply correlated with the electrical power and the heat transfer rate during the oscillation cycle. For the present device geometry ($V_{\mathrm{CB}} = 8 \times 10^{-14}$ cm³) we calculate $C_{th} = 2.4 \times 10^{-13}$ J  K⁻¹. We extract $A \overline{h}$ by following the approach used in [28]. We match the quasi-static simulated T to $T_{0} + P/(A \overline{h})$, where P is the quasi-static simulated power. We obtain that $(A \overline{h}) = 1.0753 \times 10^{-6}$ W K⁻¹ (see appendix E). It should be noted that equation (3) holds under the approximation that T is spatially uniform [37]. We have compared $\mathrm dT/\mathrm dt$ calculated from equation (3) with the temporal derivative determined directly from T versus t as obtained from simulated oscillations (figure F1 in appendix F). The match between the two curves is good overall, even though there are some discrepancies. This is expected, given that the 3D temperature distribution is not perfectly uniform.

Figure 5 shows the two components of equation (3) as plotted against $V_{\mathrm{D}} = V_{\mathrm{DD}} - V_{\mathrm{out}}$. It can be seen that both $P / C_\mathrm {th}$ and $-(T - T_{0})(A \overline{h})/C_\mathrm {th}$ are almost unaffected up to $C_{\mathrm{ext}} =$ 5 nF, and then start to depend on frequency. In particular, we can distinguish between the charging or discharging parts of an oscillation cycle. During the charging part of the oscillation cycle (when VO₂ is in the LRS), the term $-(T - T_{0})(A \overline{h})/C_\mathrm {th}$ begins to depend on frequency for $C_{\mathrm{ext}} \leqslant$ 1 nF. In particular, the absolute value of $-(T - T_{0})(A \overline{h})/C_\mathrm {th}$ increases with frequency for each value of $V_{\mathrm{D}} = V_{\mathrm{DD}} - V_{\mathrm{out}}$ when it is $\lessapprox$ 0.2 V. However, this is not due to the corresponding increase in $P / C_\mathrm {th}$, given that the latter is not dependent on $C_{\mathrm{ext}}$. Indeed, this is consistent with the fact that $P = V_{\mathrm{D}}^{2} / R_{\mathrm{met}}$ with $R_{\mathrm{met}}$ is constant when VO₂ is in the LRS (see table B2 in appendix B). Thus, the increase in the absolute value of $-(T - T_{0})(A \overline{h})/C_\mathrm {th}$ means a decrease in thermal dissipation with consequent thermal build-up in this range of frequencies. This conclusion agrees with the thermal time constant of the oscillator estimated to be $C_\mathrm {th} / (A \overline{h}) = 2.23 \times 10^{-7}$ s, which thus becomes comparable to the period of charging ($1.28 \times 10^{-6}$ s) for $C_{\mathrm{ext}} =$ 500 pF. In this way, the temperature $T_{\mathrm{MIT}} = 320$ K is achieved at lower and lower power for $C_{\mathrm{ext}} \leqslant$ 1 nF, which explains the observed dependence on $V_{\mathrm{MIT}}$ on $C_{\mathrm{ext}}$ in this range. Regarding the discharging part of the oscillation cycle, it can be seen that the term proportional to P begins to depend on frequency for $C_{\mathrm{ext}} \leqslant$ 1 nF. This is consistent with the above observation that the HRS branch of I–$V_{\mathrm{D}}$ varies with $C_{\mathrm{ext}}$ for $C_{\mathrm{ext}} \leqslant$ 1 nF.

Zoom In Zoom Out Reset image size

Finally, it is worth observing that, in principle, the modulation of IMT/MIT points from frequencies above a threshold value can be compensated by changing the load-line. Since the value of $I_{\mathrm{MIT}}$ decreases with frequency ($C_{\mathrm{ext}} \leqslant$ 1 nF) while the value of $I_{\mathrm{IMT}}$ increases with it ($C_{\mathrm{ext}} \leqslant$ 800 pF), as a consequence, 'flatter' and 'flatter' load-lines are necessary to continue to comply with the self-oscillatory regime. This is shown in figure 6. Oscillatory behavior at $C_{\mathrm{ext}} =$ 400 and 350 pF is achieved with a load-line specified by V$_{\mathrm{DD}} =$ 10 V, R$_{\mathrm{ext}} =$ 80 kΩ. The corresponding oscillation frequencies are 294 kHz ($C_{\mathrm{ext}} =$ 400 pF) and 291 kHz ($C_{\mathrm{ext}} =$ 350 pF). Oscillatory behavior at $C_{\mathrm{ext}} =$ 320 pF is achieved with a load-line specified by V$_{\mathrm{DD}} =$ 30 V, R$_{\mathrm{ext}} =$ 240 kΩ, with a corresponding oscillation frequency of 283 kHz. Thus, the load-line parameters $V_{\mathrm{DD}}$ and $R_{\mathrm{ext}}$ have to be increased more and more to maintain compliance with the self-oscillatory regime, with a not negligible cost in terms of power consumption. Finally, it becomes practically impossible to find a load-line yielding a self-oscillatory regime for $C_{\mathrm{ext}} =$ 300 pF and below.

Zoom In Zoom Out Reset image size

In this study, we use electrothermal 3D TCAD simulations to investigate the factors that may set an upper limit, if any, over the maximum frequency achievable by a given VO₂ oscillator. Our TCAD simulations show that the IMT/MIT points of current versus voltage across the VO₂ device depend on frequency, during oscillation, below some threshold values of $C_{\mathrm{ext}}$. This implies that the condition for the self-oscillatory regime may be satisfied by a given load-line in the low-frequency range, but no longer at higher frequencies, with consequent suppression of oscillations. It is important to observe that a theoretical, purely circuital model of the VO₂ oscillator, where the VO₂ device is simplified as a two-state VO₂ variable resistor whose RS is activated by threshold voltages $V_\textrm{IMT}$ and $V_\textrm{MIT}$, respectively (appendix D, [35]), is not able to predict this behavior, since such a model does not consider any thermal effects. Our more accurate TCAD simulations, where device simulation is integrated into the solution of the circuit equation, consider the electrothermal mechanism of the VO₂ oscillator, allowing us to observe the evolution at the IMT/MIT points. We explain this behavior with frequency as a combination of effects: 1) the achievement of intermediate resistive states of VO₂ modulated by frequency, making HRS conductance of VO₂ channel dependent on frequency, and 2) the interplay between frequency and heat transfer rate. Although he upper limit on the frequency that we extract is linked to the specific VO₂ device we simulate, our findings apply qualitatively to any VO₂ oscillator. Overall, this study elucidates the link between electrical and thermal behavior in VO₂ devices that sets a constraint on the upper values of the frequency of a VO₂ oscillator.

Present address of S C: Silvaco France, 55 Rue Blaise Pascal, 38330 Montbonnot-Saint-Martin, France

S C and A T S conceived the idea. S C designed and performed the TCAD simulations. S C analyzed and validated the simulated data. S C, G B and A T S critically interpreted the data. A P and A N implemented the customized version of the PCM model used to simulate the VO₂ material, and gave useful suggestions to implement TCAD simulations. S K supplied the experimental data used for the calibration of the TCAD model, and provided useful discussions to interpret them. All the authors discussed the results. S C drafted the manuscript and A T S critically revised it. The manuscript has been read and approved by all authors. A T S had the management, coordination responsibility and supervision for the planning and execution of the research activity leading to this publication. A T S acquired the financial support for the project leading to this publication.

This study was supported by the European Union's Horizon 2020 research and innovation programme, EU H2020 NEURONN (www.neuronn.eu) project under Grant No. 871501.

The data that support the findings of this study are openly available at the following URL/DOI: https://osf.io/sw7va/.

The input files to implement TCAD simulations with Silvaco Victory Mesh and Victory Device tools are provided as supporting material.

VO₂ is a thermistor material with a negative temperature coefficient, since it is a thermally sensitive resistor exhibiting a decrease in electrical resistance with an increase in its temperature. Noticeably, thermistor devices are among 'some examples of physical devices which should be modeled as memristive one-ports' [38]. The equations describing the system are as follows:

$$\begin{equation} V = R(T)I, \end{equation} \tag{ A.1 }$$

and the heat transfer equation,

$$\begin{equation} \frac{\mathrm dT}{\mathrm dt} = -\frac{\delta}{C} (T - T_{0}) - \frac{P(t)}{C} ,\end{equation} \tag{ A.2 }$$

provided that C is the heat capacitance, δ is the dissipation constant and P(T) is the electrical power. Equation (A.2) can indeed be rewritten as,

$$\begin{equation} \frac{\mathrm dT}{\mathrm dt} = -\frac{\delta}{C} (T - T_{0}) - \frac{R(T)I^{2}}{C} \equiv f(T,I) ,\end{equation} \tag{ A.3 }$$

and the combination of (A.1) and (A.3) implies that a thermistor device is a first-order time-invariant current-controlled memristive one-port.

It is important to observe that the above description is an ideal 0D model, where indeed the whole channel is considered to be a homogeneous resistor dependent on the temperature T. This description would apply, rigorously, if the VO₂ channel was self-heated homogeneously, which corresponds to the current flowing through the channel in a homogeneous way. This is not the case for a CB VO₂ device, where just the part between overlapping top and bottom contacts is involved in device conduction, making it a true 3D physical system. Consequently, (A.1) and (A.3) appear in this study in their microscopic, local version as $E = -\nabla \phi = \rho (T) J$ and equation (2), respectively, with the functional form of $\rho (T)$ given in appendix B.

In our TCAD approach, VO₂ is considered to be a conductor, whose local resistivity depends on local temperature T and time t according to equation (1). In particular, the steady-state resistivity $\hat{\rho}(T)$ is described by the set of equations listed in table B1. This set of equations apply for both (heating and cooling) branches of $\hat{\rho}(T)$, where by a heating (or cooling) branch is meant the curve $\hat{\rho}(T)$ when the external temperature is increased (or decreased) monotonically. The temperature interval $T_{\mathrm{HRS}}\lt T\lt T_{\mathrm{LRS}}$ corresponds to the one of transition from HRS to LRS (for the heating branch), or of transition from LRS to HRS (for the cooling branch). The parameters $m_{\mathrm{HRS}}$, $T_{\mathrm{HRS}}$, $\hat{\rho}(T_{\mathrm{HRS}})$, $m_{\mathrm{LRS}}$, $T_{\mathrm{LRS}}$, $\hat{\rho}(T_{\mathrm{LRS}})$, are calibrated against available experimental measurements of resistance R versus temperature T (figure B1(a)) of a real CB device. The experimental data are obtained as follows. The VO₂ resistance is measured in steady-state condition at different temperatures as provided by heating VO₂ through an external heater, without any contribution to heat generation in the VO₂ layer due to the Joule effect.

Zoom In Zoom Out Reset image size

Table B1. Analytical expressions of heating/cooling branches of TCAD $\hat{\rho}(T)$.

$\hat{\rho} (T)$	Condition
$\hat{\rho}(T_{\mathrm{HRS}}) + m_{\mathrm{HRS}} \times (T - T_{\mathrm{HRS}})$	if $T \lt T_{\mathrm{HRS}}$
$\hat{\rho}(T_{\mathrm{LRS}}) + m_{\mathrm{LRS}} \times (T - T_{\mathrm{LRS}})$	if $T \gt T_{\mathrm{LRS}}$
$\hat{\rho}(T_{\mathrm{HRS}}) + \frac{(T - T_{\mathrm{HRS}}) \times (\hat{\rho}(T_{\mathrm{LRS}}) - \hat{\rho}(T_{\mathrm{HRS}}))}{T_{\mathrm{LRS}} - T_{\mathrm{HRS}}}$	Otherwise

The experimental device is made of a square VO₂ layer with each side 5 µm and a thickness of 80 nm, while the width of the contacts is 250 nm. To perform the calibration, first we have to extract the resistivity from the experimental resistance. To do this, we use the formula $\hat{\rho} = R \times \frac{A_{C\textrm {eff}}}{L}$, where L is the 'channel' length (that is the VO₂ thickness), and $A_{C\textrm{eff}} = A_{C} \times a$, for A_C the cross-sectional area of overlapping electrodes, and a is a parameter. We determine a by performing a voltage-controlled 3D electrothermal simulation of the VO₂ device in series with an external resistor of 1 kΩ, and by fitting the HRS parts of 1) simulated I versus V_D and 2) experimental I–V_D . If the current flow is coincident with the region of overlapping between the two contacts, then it would be a = 1, while we find a = 2.7. Finally, the hysteresis of $\hat{\rho}(T)$ is implemented by using different parameters for the heating/cooling branches of $\hat{\rho}(T)$. The calibrated parameters are shown in table B2.

Table B2. Calibrated parameters of TCAD $\hat{\rho}(T)$.

	$m_{\mathrm{HRS}}$	$T_{\mathrm{HRS}}$	$\hat{\rho}(T_{\mathrm{HRS}})$	$m_{\mathrm{LRS}}$	$T_{\mathrm{LRS}}$	$\hat{\rho}(T_{\mathrm{LRS}})$
Branch	(${\mu\Omega}$cm K)−1	(K)	(${\mu\Omega}$cm)	(${\mu\Omega}$cm K)−1	(K)	(${\mu\Omega}$cm)
Heating	$-9.81 \times 10^{3}$	327.6	$5.543 \times 10^{5}$	0	340.9	$1.33 \times 10^{5}$
Cooling	$-9.81 \times 10^{3}$	304.6	$5.543 \times 10^{5}$	0	320	$1.33 \times 10^{5}$

The plots of figures 1(b) and (c) schematize the capability of the TCAD model to reproduce the experimentally observed intermediate resistive states in R versus T curves obtained when the sweep of temperature is inverted before the full HRS-to-LRS transition (or the reverse) is completed [1]. The authors of [1] have found that if, at a given intermediate resistance value during the HRS-to-LRS transition, the trend of temperature values is inverted (so that it is no longer increasing), then the resistance remains approximately constant at such an intermediate value. The same applies if, at a given intermediate resistance value during the LRS-to-HRS transition, the trend of temperature values is inverted. In our TCAD model, this experimental behavior is implemented as follows: at each simulation time t step, the $\hat{\rho}(T(t))$ is calculated at each mesh point for the local temperature T(t) at this point, and:

• I.  
  during the heating cycle, if the current T(t) is no longer higher than the temperature at previous time step T_old , then $\hat{\rho}(T(t))$ is frozen: $\hat{\rho}(T(t)) = \hat{\rho}(T_{old})$.
• II.  
  during the cooling cycle, if the current T(t) is no longer lower than the temperature at previous time step T_old , then $\hat{\rho}(T(t))$ is frozen: $\hat{\rho}(T(t)) = \hat{\rho}(T_{old})$.

The quasi-static I-$V_{\mathrm{D}}$ (azure line of figure 3(a) and cyan lines of figure C1) has been simulated by sweeping an applied voltage of 0.5 V for a ramp time of 1.875 ms onwards (and sweeping backwards to 0 V in the same time interval), where the VO₂ device is connected in series with an external resistance of 1 kΩ. Indeed, this setup corresponds to a quasi-static condition, as it is demonstrated by the fact that overlapping curves are obtained by applying the same voltage for ramptimes of 0.1875, 18.75, 187.5 ms and 1.875 s (black lines in figure C1(a)). The quasi-static condition does not apply any longer for shorter ramp times, as shown when the ramp time is 37.5 µs (orange triangles in figure C1(b)) or less. We have also simulated the current–voltage characteristic in DC mode. We have determined each point (I, $V_{\mathrm{D}}$) of the DC curve by applying a bias that is constant with time. Whenever the bias is applied, the current starts to flow across the channel, and the Joule effect produces an increase in local temperature which, in turn, modifies the local resistivity. Thus, the system needs some time to achieve equilibrium. Figure C2(a) shows $V_{\mathrm{D}}$/max($V_{\mathrm{D}}$) versus time for any bias applied. It can be observed the change of $V_{\mathrm{D}}$ with time, until an equilibrium state is achieved. In the inset, the equilibrium time is plotted, defined as the time necessary to achieve equilibrium, versus the applied bias. Overall, the DC current–voltage characteristics is obtained by repeating the simulations at different, constant applied bias. The DC I-$V_{\mathrm{D}}$ is reported in figure C2(b) (star symbols). It can be seen that the quasi-static I-$V_{\mathrm{D}}$ and the DC I-$V_{\mathrm{D}}$ overlap. Thus, the NDR region determined from the quasi-static I-$V_{\mathrm{D}}$ falls into the rigorous definition of the NDR region of a volatile memristor (see for instance [36]).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A final remark. We have derived resistance values $R_\mathrm {HRS} = 4.0763$ kΩ and $R_\mathrm {LRS} = 722.5434\,\Omega$ from the quasi-static I–$V_{\mathrm{D}}$ of VO₂ memristor. We have used these values to approximate the memristor as a two-state variable resistor in the purely circuital model of the VO₂ oscillator, which we have applied to calculate the theoretical frequencies (red asterisks) of figure 2(b) (appendix D).

Theoretical frequencies shown in figure 2(b) (red asterisks) are calculated through the purely circuital model of the VO₂ oscillator where the VO₂ volatile memristor is approximated as a two-state (LRS/HRS) variable resistor, activated by threshold voltages $V_\mathrm {IMT}$ and $V_\mathrm {MIT}$, respectively. For generality, we consider the two possible topologies of the VO₂ oscillator circuit, where $C_{\mathrm{ext}}$ is connected across $R_{\mathrm{ext}}$ (figure D1(a)), or across the VO₂ device (figure D1(b)). The equivalent Thevenin circuit (figure D1(c)) is similar, since the equivalent Thevenin resistance is $R_{\mathrm{thev}} = \frac{R_{\mathrm{VO2}} R_{\mathrm{ext}}}{R_{\mathrm{VO2}} + R_{\mathrm{ext}}}$, while the equivalent Thevenin voltages are slightly different:

• $V_{\mathrm{thev}} = \frac{R_{\mathrm{ext}} V_{\mathrm{DD}}}{R_{\mathrm{VO2}} + R_{\mathrm{ext}}}$ (circuit topology of figure D1(a))
• $V^\prime_{\mathrm{thev}} = \frac{R_{\mathrm{VO2}} V_{\mathrm{DD}}}{R_{\mathrm{VO2}} + R_{\mathrm{ext}}}$ (circuit topology of figure D1(b))

Zoom In Zoom Out Reset image size

with $R_{\mathrm{VO2}} = R_{\mathrm{LRS}}$, $R_{\mathrm{HRS}}$ depending on whether the VO₂ channel is in the ON, OFF state.

It should be observed that $C_{\mathrm{ext}}$ is charging/discharging when VO₂ is in the ON/OFF state, respectively, if circuit topology is the one shown in figure D1(a). The reverse applies, instead, if the circuit topology is the one shown in figure D1(b) (see, for instance [39]). Overall, the voltage transient at the output node is as follows:

$$\begin{equation} V_{\mathrm{out}}(t) = V_{\mathrm{final}} - (V_{\mathrm{final}} - V_{\mathrm{initial}}) \times \exp\bigg(- \frac{t}{R_{\mathrm{thev}}C_{\mathrm{ext}}}\bigg) ,\end{equation} \tag{ D.1 }$$

with $V_{\mathrm{final}} = V_{\mathrm{thev}}$. $V_{\mathrm{initial}}$ is equal to $\left\{ \begin{array}{l} V_{\mathrm{L}} = \min\limits_{t \in \tau}{V_{\mathrm{out}}} \text{if} \ \text{C}_{\mathrm{ext}}\ \text{is charging}\\ V_{\mathrm{H}} = \max\limits_{t \in \tau}{V_{\mathrm{out}}} \text{if} \ \text{C}_{\mathrm{ext}}\ \text{is discharging}\\ \end{array} \right.$ for τ oscillation period, where $\left\{ \begin{array}{l} V_{\mathrm{L}} = V_{\mathrm{DD}} - V_{\mathrm{IMT}} \text{and } V_{\mathrm{H}} = V_{\mathrm{DD}} - V_{\mathrm{MIT}} \text{for circuit topology of figure D1(a)}\\ V_{\mathrm{L}} = V_{\mathrm{MIT}} \text{and } V_{\mathrm{H}} = V_{\mathrm{IMT}} \text{for circuit topology of figure D1(b)}\\ \end{array} \right.$.

The time of charging/discharging can be derived from (D.1) as,

$$\begin{equation} t = \bigg[\ln\bigg( \frac{V_{\mathrm{initial}} - V_{\mathrm{final}}}{V_{\mathrm{out}}(t) - V_{\mathrm{final}}}\bigg)\bigg]R_{\mathrm{thev}}C_{\mathrm{ext}} ,\end{equation} \tag{ D.2 }$$

and the period can be estimated as $\tau = \tau _{\mathrm{ch}} + \tau _{\mathrm{disch}}$, where $\tau _{\mathrm{ch}} = t_{\mathrm{initial}} - t_{\mathrm{switch}}$ is the time interval of $C_{\mathrm{ext}}$ charging and $\tau _{\mathrm{disch}} = t^\prime_{\mathrm{initial}} - t^\prime_{\mathrm{switch}}$ is the time interval of $C_{\mathrm{ext}}$ discharging.

This brings us to the following set of equations for the circuit topology of figure D1(a):

$$\begin{eqnarray} \kern-22.9pt \tau_{\mathrm{ch}} & = \bigg[\ln\bigg( \frac{V_{\mathrm{L}} - V_{\mathrm{thev}}(LRS)}{V_{\mathrm{H}} - V_{\mathrm{thev}}(LRS)}\bigg)\bigg]R_{\mathrm{thev}}(LRS)C_{\mathrm{ext}},\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray} & = \bigg[\ln\bigg( \frac{V_{\mathrm{DD}} - V_{\mathrm{IMT}} - \frac{R_{\mathrm{ext}} V_{\mathrm{DD}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}}{V_{\mathrm{DD}} - V_{\mathrm{MIT}} - \frac{R_{\mathrm{ext}} V_{\mathrm{DD}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}}\bigg)\bigg]\frac{R_{\mathrm{LRS}} R_{\mathrm{ext}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}C_{\mathrm{ext}}. \end{eqnarray} \tag{ D.4 }$$

$$\begin{eqnarray} \kern-31.2pt \tau_{\mathrm{disch}} & = \bigg[\ln\bigg( \frac{V_{\mathrm{H}} - V_{\mathrm{thev}}(HRS)}{V_{\mathrm{L}} - V_{\mathrm{thev}}(HRS)}\bigg)\bigg]R_{\mathrm{thev}}(HRS)C_{\mathrm{ext}} ,\end{eqnarray} \tag{ D.5 }$$

$$\begin{eqnarray} & = \bigg[\ln\bigg( \frac{V_{\mathrm{DD}} - V_{\mathrm{MIT}} - \frac{R_{\mathrm{ext}} V_{\mathrm{DD}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}}{V_{\mathrm{DD}} - V_{\mathrm{IMT}} - \frac{R_{\mathrm{ext}} V_{\mathrm{DD}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}}\bigg)\bigg]\frac{R_{\mathrm{HRS}} R_{\mathrm{ext}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}C_{\mathrm{ext}}. \end{eqnarray} \tag{ D.6 }$$

Thus, $\tau = \beta C_{\mathrm{ext}}$ for $\beta = \ln\bigg( \frac{V_{\mathrm{DD}} - V_{\mathrm{IMT}} - \frac{R_{\mathrm{ext}} V_{\mathrm{DD}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}}{V_{\mathrm{DD}} - V_{\mathrm{MIT}} - \frac{R_{\mathrm{ext}} V_{\mathrm{DD}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}}\bigg)\frac{R_{\mathrm{LRS}}R_{\mathrm{ext}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}} + \ln\bigg( \frac{V_{\mathrm{DD}} - V_{\mathrm{MIT}} - \frac{R_{\mathrm{ext}} V_{\mathrm{DD}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}}{V_{\mathrm{DD}} - V_{\mathrm{IMT}} - \frac{R_{\mathrm{ext}} V_{\mathrm{DD}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}}\bigg)\frac{R_{\mathrm{HRS}} R_{\mathrm{ext}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}$.

For the circuit topology of figure D1(b):

$$\begin{eqnarray} \tau_{\mathrm{ch}} & = \bigg[\ln\bigg(\frac{V_{\mathrm{L}} - V_{\mathrm{thev}}(HRS)}{V_{\mathrm{H}} - V_{\mathrm{thev}}(HRS)}\bigg)\bigg]R_{\mathrm{thev}}(HRS)C_{\mathrm{ext}} ,\end{eqnarray} \tag{ D.7 }$$

$$\begin{eqnarray} & = \bigg[\ln\bigg( \frac{V_{\mathrm{MIT}} - \frac{R_{\mathrm{HRS}} V_{\mathrm{DD}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}}{V_{\mathrm{IMT}} - \frac{R_{\mathrm{HRS}} V_{\mathrm{DD}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}}\bigg)\bigg]\frac{R_{\mathrm{HRS}} R_{\mathrm{ext}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}C_{\mathrm{ext}}. \end{eqnarray} \tag{ D.8 }$$

$$\begin{eqnarray} \kern-13.8pt \tau_{\mathrm{disch}} & = \bigg[\ln\bigg( \frac{V_{\mathrm{H}} - V_{\mathrm{thev}}(LRS)}{V_{\mathrm{L}} - V_{\mathrm{thev}}(LRS)}\bigg)\bigg]R_{\mathrm{thev}}(LRS)C_{\mathrm{ext}} ,\end{eqnarray} \tag{ D.9 }$$

$$\begin{eqnarray} & = \bigg[\ln\bigg( \frac{V_{\mathrm{IMT}} - \frac{R_{\mathrm{LRS}} V_{\mathrm{DD}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}}{V_{\mathrm{MIT}} - \frac{R_{\mathrm{LRS}} V_{\mathrm{DD}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}}\bigg)\bigg]\frac{R_{\mathrm{LRS}} R_{\mathrm{ext}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}C_{\mathrm{ext}}. \end{eqnarray} \tag{ D.10 }$$

Thus, $\tau = \beta^{\prime} C_{\mathrm{ext}}$ for $\beta^{\prime} = \ln\bigg( \frac{V_{\mathrm{MIT}} - \frac{R_{\mathrm{HRS}} V_{\mathrm{DD}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}}{V_{\mathrm{IMT}} - \frac{R_{\mathrm{HRS}} V_{\mathrm{DD}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}}}\bigg)\frac{R_{\mathrm{HRS}}R_{\mathrm{ext}}}{R_{\mathrm{HRS}} + R_{\mathrm{ext}}} + \ln\bigg( \frac{V_{\mathrm{IMT}} - \frac{R_{\mathrm{LRS}} V_{\mathrm{DD}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}}{V_{\mathrm{MIT}} - \frac{R_{\mathrm{LRS}} V_{\mathrm{DD}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}}\bigg)\frac{R_{\mathrm{LRS}} R_{\mathrm{ext}}}{R_{\mathrm{LRS}} + R_{\mathrm{ext}}}$.

The theoretical frequencies plotted in figure 2(b) have been extracted from the periods $\tau = \beta C_{\mathrm{ext}}$ of VO₂ oscillator with the circuit topology drawn in figure D1(a). They are listed in table D2. The parameters used to perform calculations are listed in table D1, with the following: $V_{\mathrm{thev}}(HRS) =$ 1.572 632 678 8559 V, $V_{\mathrm{thev}}(LRS) =$ 1.908 088 235 2941 V, $V_{\mathrm{L}} = V_{\mathrm{DD}} - V_{\mathrm{IMT}} =$ 1.712 915 98 V, $V_{\mathrm{H}} = V_{\mathrm{DD}} - V_{\mathrm{MIT}} =$ 1.881 303 54 V. We have calculated the theoretical frequencies extracted from the periods $\tau = \beta^{\prime} C_{\mathrm{ext}}$ of the VO₂ oscillator with the circuit topology drawn in figure D1(b). Indeed, we have found that $\beta = 3896.762\,624\,63 \text{s / F} = \beta^{\prime}$, which means that the two circuit topologies yield the same period for any given $C_{\mathrm{ext}}$ value. The parameters used to perform these latter calculations are listed in table D1, with the following: $V^{^{\prime}}_{\mathrm{thev}}(HRS) =$ 0.427 367 321 144 V, $V^\prime_{\mathrm{thev}}(LRS) =$ 0.091 911 764 7059 V, $V_{\mathrm{L}} = V_{\mathrm{MIT}}$, $V_{\mathrm{H}} = V_{\mathrm{IMT}}$.

Table D1. Parameters used to calculate theoretical frequencies.

$R_{\mathrm{HRS}}$	$R_{\mathrm{LRS}}$	$V_{\mathrm{IMT}}$	$V_{\mathrm{MIT}}$	$R_{\mathrm{thev}}(HRS)$	$R_{\mathrm{thev}}(LRS)$
4.0763 kΩ	722.5434 Ω	$0.287\,084\,02$ V	$0.118\,696\,46$ V	$3.205\,254\,908\,581$ kΩ	$689.338\,235\,2941$ Ω

Table D2. Frequencies calculated by theoretical circuit model.

$C_{\mathrm{ext}}$ (pF)	104	$5 \times 10^{3}$	103	800	600	500	400	300
Frequency (kHz)	25.7	51.3	256.6	320.8	427.7	513.2	641.6	855.4
$C_{\mathrm{ext}}$ (pF)	200	100	10
Frequency (kHz)	1.283 $\times 10^{3}$	2.566 $\times 10^{3}$	2.566 $\times 10^{4}$

Table D3. TCAD Simulated Frequencies.

$C_{\mathrm{ext}}$ (pF)	104	$5 \times 10^{3}$	103	800	600	500
Frequency (kHz)	21.6	41.678	157.9	182.935	210.2	221.4

To determine $A \overline{h}$, we follow [28]. We match the temperature T, as probed in the center of the active region of the device, with $T_{0} + P/(A \overline{h})$, for P the device power. A is the heat transfer surface area and $\overline{h}$ the heat transfer coefficient averaged over A [37]. Both T and P are extracted from TCAD 3D electrothermal simulation of quasi-static I versus $V_{\mathrm{D}}$. The use of $T = T_{0} + P/(A \overline{h})$ is equivalent to modeling the heating effects in the VO₂ device with the theory for thermistors [40]. We achieve an excellent agreement (figure E1) by $A \overline{h} = 1.0753 \times 10^{-6}$ W K⁻¹.

Zoom In Zoom Out Reset image size