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

The frequency of vanadium dioxide (VO2) oscillators is a fundamental figure of merit for the realization of neuromorphic circuits called oscillatory neural networks (ONNs) since the high frequency of oscillators ensures low-power consuming, real-time computing ONNs. In this study, we perform electrothermal 3D technology computer-aided design (TCAD) simulations of a VO2 relaxation oscillator. We find that there exists an upper limit to its operating frequency, where such a limit is not predicted from a purely circuital model of the VO2 oscillator. We investigate the intrinsic physical mechanisms that give rise to this upper limit. Our TCAD simulations show that it arises a dependence on the frequency of the points of the curve current versus voltage across the VO2 device corresponding to the insulator-to-metal transition (IMT) and metal-to-insulator transition (MIT) during oscillation, below some threshold values of Cext . 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 note that this variation of the IMT/MIT points below some threshold values of Cext is due to a combination of different factors: intermediate resistive states achievable by VO2 channel and the interplay between frequency and heat transfer rate. Although the upper limit on the frequency that we extract is linked to the specific VO2 device we simulate, our findings apply qualitatively to any VO2 oscillator. Overall, our study elucidates the link between electrical and thermal behavior in VO2 devices that sets a constraint on the upper values of the operating frequency of any VO2 oscillator.


Introduction
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 2 or NbO 2 . 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 2 exhibits a steep RS due to the change from insulating monoclinic (M 1 ) to metallic tetragonal rutile-like (R) crystal structure. The IMT M 1 −→ R is driven by temperature and occurs at a threshold temperature T IMT of about 340 K. The reverse metal-to-insulator transition (MIT) R −→ M 1 occurs during cooling down of the material, when a threshold temperature T MIT is achieved. Generally, T MIT is smaller than T IMT . Thus, the experimental curve of resistance R against temperature T shows hysteresis.
RS is also induced in VO 2 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][3][4][5]. The channel of VO 2 devices switches from an high resistance state (HRS) to a low resistance state (LRS) depending on the VO 2 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 2 . Consequently, VO 2 devices are named volatile Mott memristors [6][7][8][9].
A relaxation oscillator, which generates a non-sinusoidal periodic waveform at its output, can be obtained by connecting in series a VO 2 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 ext , where these cycles are self-piloted by the VO 2 volatile memristor being periodically in the ON/OFF (LRS/HRS) state. We define (V IMT , I IMT ) to be the IMT point of the curve current (I) versus voltage across the VO 2 device (V D ), which is the point of the curve where the IMT occurs. Similarly, (V MIT , I 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 DD and external resistance R ext , are such that the load-line current (I L = (V DD − V D )/R ext ) is larger than the I IMT at V IMT , and smaller than I MIT at V MIT . It is important to observe that VO 2 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.
Recently, oscillatory neural networks (ONNs) based on VO 2 oscillators have been under intense scrutiny [10][11][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][17][18][19][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 2 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 2 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 2 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 2 device starts to depend on frequency, during oscillation, below some threshold values of C 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 2 achieved below some threshold values of C ext that make HRS conductance of VO 2 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 2 device we simulate, our findings apply qualitatively to any VO 2 oscillator. Overall, this study elucidates the link between ). Then, the IMT transition occurs, and the material switches to a low resistance state (LRS). Effect is volatile: the material resets to the original HRS if it is cooled down, usually below different threshold temperature T MIT . In the latter case, the curveρ versus T shows an hysteresis window. If the IMT is not so abrupt, then it will occur within an interval of temperatures THRS < T < TLRS. (b) During IMT (THRS(H) < T < TLRS(H)) the temperature should always be increasing. If, instead, at a given moment, the material starts to cool, thenρ remains constant, until the branch of the cooling cycle is reached. (c) During MIT (THRS(C) < T < TLRS(C)) the temperature should always be decreasing. If, instead, at a given moment the material starts to warm, thenρ remains constant, until the branch of the heating cycle is reached. This behavior emulates the observed experimental behavior of R versus T in VO2 [1]. (d) Crossbar (CB) architecture of a VO2 volatile memristor. (e) Schematics of a VO2 oscillator circuit. I, IL and IC are the currents across the VO2 device, Rext and Cext, respectively. Representative, non-sinusoidal waveform of voltage at the output node over one period is also shown. (f) Self-oscillatory electrical behavior is associated with the compliance of the load-line with the self-oscillatory regime. electrical and thermal behavior in VO 2 devices that sets a constraint on the upper values of the frequency of a VO 2 oscillator.

TCAD simulation approach for VO 2 volatile memristors and oscillators
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 2 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 2 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 2 channel is mimicked with a good approximation, which plays a key role in modulating the behavior of VO 2 volatile memristors and oscillators [26][27][28].
We use Silvaco Victory Device [29] to solve PDEs of VO 2 device physics on the generated mesh. We consider VO 2 as a conductor, and we use a customized C-Interpreter function [26][27][28] of the PCM model [29] to empirically emulate VO 2 resistivity. Overall, VO 2 resistivity depends on local temperature T and time t after the Johnson-Mehl-Avrami model [30][31][32]: where t is time, ∆t is the time step, T is the local temperature,ρ is the steady-state resistivity and τ is a lifetime. We calibrate the steady-state resistivityρ versus T (figure 1(a)) using experimental data of resistance R exp versus temperature T [26]. Further details on the VO 2 resistivity model and the calibration procedure are contained in appendix B. Since VO 2 is described as a conductor, the current density is J = E/ρ = −∇ϕ /ρ, for E the electric field and ϕ the electric potential. Thus, to calculate self-consistently the current transport and temperature distribution inside the VO 2 volatile memristor we need to solve the Poisson equation coupled to the heat flow equation: where K = 0.06 W (cm K) −1 [33] is the thermal conductivity, and C = 3 J (cm 3 K) −1 [34] is the heat capacity per unit volume of VO 2 , respectively. The term (∇ϕ) 2 /ρ corresponds to the heat generation rate ρJ 2 due to the Joule effect. For simplicity, we consider the bottom surface of the VO 2 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 ( for T 0 the external temperature, R Th the interface thermal resistance and A the bottom surface area of the device. The term 1/(R Th A ) regulates the heat dissipation of the VO 2 volatile memristor with respect to the environment. Finally, we perform simulations of the VO 2 relaxation oscillator by combining the TCAD engine to solve the PDEs of the VO 2 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.

Relationship between RS and oscillator frequency
We perform electrothermal 3D TCAD mixed-mode simulations of a VO 2 oscillator (circuit scheme in figure 1(e)), where the VO 2 CB device is made of an 80 nm thick VO 2 layer inserted between the top and bottom Pt electrodes. The electrodes are 250 nm wide, while the VO 2 layer is a square with each side 1 µm long. We set T 0 = 303 K and R Th = 8.3 × 10 5 K W −1 as boundary conditions at the thermal contact. We choose the circuit load-line specified by V DD = 2 V, R ext = 15 kΩ in order to comply with the self-oscillatory regime. We modulate the operational oscillator frequency by varying C ext in the interval 5 nF-400 pF. The simulated voltages at the output node V out are shown in figure 2(a). The oscillations are turned-off at C ext = 400 pF. Frequency versus C 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 2 oscillator where the VO 2 volatile memristor is approximated as a two-state (LRS/HRS) variable resistor, activated by threshold voltages V IMT and V MIT , respectively (see appendix D, and [35]). It can be seen that for larger values of C ext the simulated and the theoretically calculated data are almost overlapping. However, the theoretical calculations predict frequencies corresponding to C ext values equal to or below 400 pF, while we are not able to obtain oscillations in this range of C 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 2 oscillator (appendix D, [35]).
In figure 3(a), we plot curves I versus V D = V DD − V out (symbols) as extracted from simulated oscillations for C 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 ext = 5 nF) all the curves essentially overlap, and then some mismatch becomes evident at higher-frequency values (lower C ext values). In particular, 1. the IMT/MIT points become dependent from C ext below some threshold value (C ext ⩽ 800 pF for IMT points, C ext ⩽ 1 nF for MIT points); see insets of figure 3(a), 2. the HRS branch of the I-V D = V DD − V out curve is also dependent from C ext for C ext ⩽ 1 nF.
We compare the above curves with the quasi-static I-V D (azure thick line of figure 3(a)). The curve is obtained by electrothermal 3D TCAD simulation of a voltage-controlled VO 2 volatile memristor, where the VO 2 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 D are the same in quasi-static conditions and during oscillations up to when a dependence on C ext arises of I-V D = V DD − V out curves. This dependence on C ext of IMT/MIT points implies an apparent modulation of the negative differential resistance (NDR) region with frequency. Indeed, the NDR region of a Blue-diamond curve is from electrothermal 3D TCAD mixed-mode simulated oscillations. In this case, the frequency of oscillations ranges from 41.7 kHz (Cext = 5 nF) up to 221.4 kHz (Cext = 500 pF). At 400 pF the oscillations are inhibited. Red-asterisk curve is obtained by theoretical calculation through a purely circuital model of the VO2 oscillator where the VO2 device is simplified as a two-state VO2 variable resistor whose RS is activated by threshold voltages V IMT /V MIT (appendix D). In this case, frequencies for Cext values are predicted well below 400 pF. volatile memristor is rigorously defined from its DC current-voltage characteristics [36]. In our case, since the quasi-static I-V D and the DC I-V D overlap (as demonstrated in figure C2(b) in appendix C), the NDR of the quasi-static I-V D in figure 3(a) is a rigorous reference. It should be observed that the value of current at the IMT point, I IMT , increases with frequency for C ext ⩽ 800 pF, while the value of current at the MIT point, I MIT , decreases with frequency for C ext ⩽ 1 nF. Since the load-line current is independent of C ext , then it may be expected that it will become smaller than the I IMT at V IMT , and/or larger than the I MIT at V MIT for a certain value of C 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 2 oscillator, we need to investigate the mechanism of variation of the IMT/MIT points, and of the HRS branch of I-V D = V DD − V 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 D = V DD − V out for C 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 D the associated T-V D (azure thick line in figure 3(b)). It can be seen that: 3. the temperatures T IMT , associated with IMT points, become dependent on C ext for C ext ⩽ 800 pF, while the temperatures T MIT , associated with MIT points, are constant and equal to 320 K.
In particular, the points (V IMT , T IMT ) appear (approximately) located along the part of quasi-static T-V D corresponding to HRS-to-LRS transition (dashed red line). In the case of C ext = 400 pF, T is constant, with a value slightly below the HRS-to-LRS branch of quasi-static T-V D . Thus, T does not achieve the value for RS at C 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 osc ) for an oscillation cycle, where C ext = 10 nF, 5 nF, 1 nF, 800 pF, 500 pF and 400 pF. The use of t/period 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 osc for an oscillation cycle. The T versus t/period osc (and ρ versus t/period osc ) curves show dependency on frequency for C ext ⩽ 1 nF, even though the qualitative behavior is the same. Finally, figure 4(c) reports the plot of ρ versus T during an oscillation cycle. Figure 4(a) shows that the temperature increases abruptly to the maximum T 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 ρ 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 2 is not related to any removal of electrical stimulus V 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 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 ext up to C ext = 1 nF. For C ext = 1 nF, ρ at the beginning of discharging has the same value of ρ IMT = ρ(T IMT ) of quasi-static I-V 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 D = V DD − V out differs from the one of the quasi-static I-V D . Finally, for C ext ⩽ 800 pF, the value of ρ at the beginning of discharging is smaller than ρ IMT . Thus, the HRS heating branch is never recovered, and not both the HRS part of I-V D = V DD − V out and also the IMT point start to depend on C 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 dT/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]. is not correlated to an increase in P/C th , and thus it can only be explained as being due to a reduction in the efficiency of thermal dissipation in this range of frequencies.
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 dT/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: where P is the electrical power from simulated oscillations, C th = CV CB is the thermal capacity where V CB is the volume of the CB device, A is the heat transfer surface area and 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 CB = 8 × 10 −14 cm 3 ) we calculate C th = 2.4 × 10 −13 J K −1 . We extract Ah by following the approach used in [28]. We match the quasi-static simulated T to T 0 + P/(Ah), where P is the quasi-static simulated power. We obtain that (Ah) = 1.0753 × 10 −6 W K −1 (see appendix E). It should be noted that equation (3) holds under the approximation that T is spatially uniform [37]. We have compared dT/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 D = V DD − V out . It can be seen that both P/C th and −(T − T 0 )(Ah)/C th are almost unaffected up to C 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 2 is in the LRS), the term −(T − T 0 )(Ah)/C th begins to depend on frequency for C ext ⩽ 1 nF. In particular, the absolute value of −(T − T 0 )(Ah)/C th increases with frequency for each value of V D = V DD − V out when it is ⪅ 0.2 V. However, this is not due to the corresponding increase in P/C th , given that the latter is not dependent on C ext . Indeed, this is consistent with the fact that P = V 2 D /R met with R met is constant when VO 2 is in the LRS (see table B2 in appendix B). Thus, the increase in the absolute value of −(T − T 0 )(Ah)/C 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 th /(Ah) = 2.23 × 10 −7 s, which thus becomes comparable to the period of charging (1.28 × 10 −6 s) for C ext = 500 pF. In this way, the temperature T MIT = 320 K is achieved at lower and lower power for C ext ⩽ 1 nF, which explains the observed dependence on V MIT on C 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 ext ⩽ 1 nF. This is consistent with the above observation that the HRS branch of I-V D varies with C ext for C ext ⩽ 1 nF.
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 MIT decreases with frequency (C ext ⩽ 1 nF) while the value of I IMT increases with it (C ext ⩽ 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 ext = 400 and 350 pF is achieved with a load-line specified by V DD = 10 V, R ext = 80 kΩ. The corresponding oscillation frequencies are 294 kHz (C ext = 400 pF) and 291 kHz (C ext = 350 pF). Oscillatory behavior at C ext = 320 pF is achieved with a load-line specified by V DD = 30 V, R ext = 240 kΩ, with a corresponding oscillation frequency of 283 kHz. Thus, the load-line parameters V DD and R 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 ext = 300 pF and below.

Conclusion
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 2 oscillator. Our TCAD simulations show that the IMT/MIT points of current versus voltage across the VO 2 device depend on frequency, during oscillation, below some threshold values of C 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 2 oscillator, where the VO 2 device is simplified as a two-state VO 2 variable resistor whose RS is activated by threshold voltages V IMT and V 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 2 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 2 modulated by frequency, making HRS conductance of VO 2 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 2 device we simulate, our findings apply qualitatively to any VO 2 oscillator. Overall, this study elucidates the link between electrical and thermal behavior in VO 2 devices that sets a constraint on the upper values of the frequency of a VO 2 oscillator.

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

Acknowledgments
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 2 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.

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

Appendix A. VO 2 devices as memristive one-ports
VO 2 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: and the heat transfer equation, 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, 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 2 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 2 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 = −∇ϕ = ρ(T)J and equation (2), respectively, with the functional form of ρ(T) given in appendix B.

Appendix B. TCAD resistivity model of VO 2 and calibration procedure
In our TCAD approach, VO 2 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ρ(T) is described by the set of equations listed in table B1. This set of equations apply for both (heating and cooling) branches ofρ(T), where by a heating (or cooling) branch is meant the curveρ(T) when the external temperature is increased (or decreased) monotonically. The temperature interval T HRS < T < T LRS corresponds to the one of transition from HRS to LRS (for the heating branch), or of transition from LRS to

ρ(T)
Condition   The experimental data are obtained as follows. The VO 2 resistance is measured in steady-state condition at different temperatures as provided by heating VO 2 through an external heater, without any contribution to heat generation in the VO 2 layer due to the Joule effect. The experimental device is made of a square VO 2 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ρ = R × A Ceff L , where L is the 'channel' length (that is the VO 2 thickness), and A Ceff = A C × 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 2 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ρ(T) is implemented by using different parameters for the heating/cooling branches ofρ(T). The calibrated parameters are shown in table B2.
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ρ(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ρ(T(t)) is frozen:ρ(T(t)) =ρ(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ρ(T(t)) is frozen:ρ(T(t)) =ρ(T old ).  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 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 D /max(V D ) versus time for any bias applied. It can be observed the change of V 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.

Appendix C. Simulation of quasi-static I-V
Overall, the DC current-voltage characteristics is obtained by repeating the simulations at different, constant applied bias. The DC I-V D is reported in figure C2(b) (star symbols). It can be seen that the quasi-static I-V D and the DC I-V D overlap. Thus, the NDR region determined from the quasi-static I-V D falls into the rigorous definition of the NDR region of a volatile memristor (see for instance [36]). Figure C1. I-VD obtained by sweeping an applied voltage of 0.5 V for different ramp times: (a) equal to or greater than 0.1875 ms, or (b) equal to or smaller than 37.5 µs. In both panels, the cyan line is the curve corresponding to a ramp time of 1.875 ms.
A final remark. We have derived resistance values R HRS = 4.0763 kΩ and R LRS = 722.5434 Ω from the quasi-static I-V D of VO 2 memristor. We have used these values to approximate the memristor as a two-state variable resistor in the purely circuital model of the VO 2 oscillator, which we have applied to calculate the theoretical frequencies (red asterisks) of figure 2(b) (appendix D).

Appendix D. Frequencies calculated by theoretical circuit model of the VO 2 oscillator
Theoretical frequencies shown in figure 2(b) (red asterisks) are calculated through the purely circuital model of the VO 2 oscillator where the VO 2 volatile memristor is approximated as a two-state (LRS/HRS) variable resistor, activated by threshold voltages V IMT and V MIT , respectively. For generality, we consider the two possible topologies of the VO 2 oscillator circuit, where C ext is connected across R ext ( figure D1(a)), or across the VO 2 device ( figure D1(b)). The equivalent Thevenin circuit (figure D1(c)) is similar, since the equivalent Thevenin resistance is R thev = RVO2Rext RVO2+Rext , while the equivalent Thevenin voltages are slightly different: with R VO2 = R LRS , R HRS depending on whether the VO 2 channel is in the ON, OFF state. It should be observed that C ext is charging/discharging when VO 2 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: , The time of charging/discharging can be derived from (D.1) as, and the period can be estimated as τ = τ ch + τ disch , where τ ch = t initial − t switch is the time interval of C ext charging and τ disch = t ′ initial − t ′ switch is the time interval of C ext discharging. This brings us to the following set of equations for the circuit topology of figure  We have calculated the theoretical frequencies extracted from the periods τ = β ′ C ext of the VO 2 oscillator with the circuit topology drawn in figure D1(b). Indeed, we have found that β = 3896.762 624 63s / F = β ′ , which means that the two circuit topologies yield the same period for any given C ext value. The parameters used to perform these latter calculations are listed in table D1, with the following: