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

Stefania Carapezzi, Andrew Plews, Gabriele Boschetto, Ahmed Nejim, Siegfried Karg, Aida Todri-Sanial

► To cite this version:

Stefania Carapezzi, Andrew Plews, Gabriele Boschetto, Ahmed Nejim, Siegfried Karg, et al.. How fast can vanadium dioxide neuron-mimicking devices oscillate? Physical mechanisms limiting the frequency of vanadium dioxide oscillators. Neuromorphic Computing and Engineering, 2023, 3 (3), pp.034010. 10.1088/2634-4386/acf2bf. lirmm-04230790

HAL Id: lirmm-04230790

https://hal-lirmm.ccsd.cnrs.fr/lirmm-04230790

Submitted on 6 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED
2 May 2023

REVISED
1 August 2023

ACCEPTED FOR PUBLICATION
22 August 2023

PUBLISHED
11 September 2023

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

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

S Carapezzi¹,*, A Plews², G Boschetto¹, A Nejim², S Karg³ and A Todri-Sanial¹,⁴,*

¹ Microelectronics Department, LIRMM, Univ. of Montpellier, CNRS, 161 Rue Ada, 34095 Montpellier, France

² Silvaco Europe Ltd, St Ives PE27 5JL, United Kingdom

³ Department of Science and Technology, IBM Research Europe–Zurich, Ruschlikon 8803, Switzerland

⁴ Department of Electrical Engineering, Eindhoven University of Technology, Eindhoven 5612 AP, Netherlands

* Authors to whom any correspondence should be addressed.

E-mail: stefania.carapezzi@silvaco.com and a.todri.sanial@tue.nl

Keywords: relaxation oscillator, vanadium dioxide, frequency, device simulation, oscillatory neural network, neuromorphic circuits

Supplementary material for this article is available online

## Abstract

The frequency of vanadium dioxide (VO₂) 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 VO₂ 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 VO₂ 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 VO₂ device corresponding to the insulator-to-metal transition (IMT) and metal-to-insulator transition (MIT) during oscillation, below some threshold values of Cₐₓ. 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 Cₐₓ is due to a combination of different factors: intermediate resistive states achievable by VO₂ 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 VO₂ device we simulate, our findings apply qualitatively to any VO₂ oscillator. Overall, our study elucidates the link between electrical and thermal behavior in VO₂ devices that sets a constraint on the upper values of the operating frequency of any VO₂ oscillator.

## 1. 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₂ 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₁ → 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₁ 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₂ 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 a 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ₐₓ, where these cycles are self-piloted by the VO₂ 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₂ device (V<sub>D</sub>), which is the point of the curve where the IMT occurs. Similarly, (VᵢMIT, IᵢMIT) is the MIT point of the I-V<sub>D</sub> curve where the MIT occurs. Then, the condition for the self-oscillatory regime is that the fixed oscillator circuit parameters, applied bias V<sub>DD</sub> and external resistance R<sub>ext</sub>, are such that the load-line current (I<sub>L</sub> = (V<sub>DD</sub> − V<sub>D</sub>)/R<sub>ext</sub>) is larger than the I<sub>IMT</sub> at V<sub>IMT</sub>, and smaller than I<sub>MIT</sub> at V<sub>MIT</sub>. 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.

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<sub>ext</sub>. 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<sub>ext</sub> 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
<figure><img src="image_4.png" /><figcaption>Figure 1. (a)–(c) TCAD model of steady-state resistivity $\hat{\rho}$ for temperature-triggered RS of a VO$_2$ volatile memristor. (a) Material is in a high resistance state (HRS) below threshold temperature $T_{IMT}$ (about 340 K). 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 $\hat{\rho}$ versus $T$ shows an hysteresis window. If the IMT is not so abrupt, then it will occur within an interval of temperatures $T_{HRS} < T < T_{LRS}$. (b) During IMT ($T_{HRS}(H) < T < T_{LRS}(H)$) the temperature should always be increasing. If, instead, at a given moment, the material starts to cool, then $\hat{\rho}$ remains constant, until the branch of the cooling cycle is reached. (c) During IMT ($T_{HRS}(C) < T < T_{LRS}(C)$) the temperature should always be decreasing. If, instead, at a given moment the material starts to warm, then $\hat{\rho}$ remains constant, until the branch of the heating cycle is reached. This behavior emulates the observed experimental behavior of R versus T in VO$_2$ [1]. (d) Crossbar (CB) architecture of a VO$_2$ volatile memristor. (e) Schematics of a VO$_2$ oscillator circuit. $I$, $I_L$ and $I_C$ are the currents across the VO$_2$ device, $R_{ext}$ and $C_{ext}$, 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.</figcaption></figure>

electrical and thermal behavior in VO$_2$ devices that sets a constraint on the upper values of the frequency of a VO$_2$ oscillator.

## 2. 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–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–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–32]:

$$ \rho(t + \Delta t, T) = \rho(t, T) + [\hat{\rho}(T) - \rho(t, T)] \times (1 - e^{-\Delta t/\tau}), \quad (1) $$
where *t* is time, Δ*t* is the time step, *T* is the local temperature, $\hat{\rho}$ is the steady-state resistivity and $\tau$ is a lifetime. We calibrate the steady-state resistivity $\hat{\rho}$ versus *T* (figure 1(a)) using experimental data of resistance $R_{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 $\phi$ 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:

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

where $K = 0.06 \text{ W (cm K)}^{-1}$ [33] is the thermal conductivity, and $C = 3 \text{ J (cm}^3 \text{ K)}^{-1}$ [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_{Th}\mathcal{A})$, for $T_0$ the external temperature, $R_{Th}$ the interface thermal resistance and $\mathcal{A}$ the bottom surface area of the device. The term $1/(R_{Th}\mathcal{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.

### 3. Relationship between RS and oscillator frequency

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 \text{ K}$ and $R_{Th} = 8.3 \times 10^5 \text{ K W}^{-1}$ as boundary conditions at the thermal contact. We choose the circuit load-line specified by $V_{DD} = 2 \text{ V}$, $R_{ext} = 15 \text{ k}\Omega$ 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 \text{ 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₂ oscillator where the VO₂ 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₂ 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 \text{ 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 \text{ 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} \le 800 \text{ pF}$ for IMT points, $C_{ext} \le 1 \text{ 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} \le 1 \text{ 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₂ 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_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
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Electrothermal 3D TCAD mixed-mode simulations of voltages at the output node $V_{out}$ for $C_{ext}$ in the interval 5 nF–400 pF. Intrinsic device capacitance is accurately considered by simulation of device electrostatics. Parameters of circuit load-line: $V_{DD}$ = 2 V, $R_{ext}$ = 15 kΩ. (b) Frequency versus $C_{ext}$. Blue-diamond curve is from electrothermal 3D TCAD mixed-mode simulated oscillations. In this case, the frequency of oscillations ranges from 41.7 kHz ($C_{ext}$ = 5 nF) up to 221.4 kHz ($C_{ext}$ = 500 pF). At 400 pF the oscillations are inhibited. Red-asterisk curve is obtained by theoretical calculation through a purely circuitual 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_{IMT}/V_{MIT}$ (appendix D). In this case, frequencies for $C_{ext}$ values are predicted well below 400 pF.</figcaption></figure>

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} \le 800$ pF, while the value of current at the MIT point, $I_{MIT}$, decreases with frequency for $C_{ext} \le 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₂ 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} \le 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.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Curves of *I* versus *V*<sub>D</sub> = *V*<sub>DD</sub> − *V*<sub>out</sub> (symbols) as extracted from simulated oscillations (10 nF–400 pF). For comparison, quasi-static *I* versus *V*<sub>D</sub> (azure thick line) is also shown. Details about simulation of quasi-static *I* versus *V*<sub>D</sub> are reported in appendix C. Two insets above the panel show a magnification of the IMT (on the left; green dots) and IMT points (on the right; red dots). (b) Corresponding quasi-static *T* versus *V*<sub>D</sub> (azure thick line) and *T* versus *V*<sub>D</sub> = *V*<sub>DD</sub> − *V*<sub>out</sub> from simulated oscillations (symbols). In both panels (a) and (b) the HRS/LRS branches are highlighted.</figcaption></figure>

In figure 4(a), it is shown the plot of *T* against time normalized to the oscillation period (*t/period*<sub>osc</sub>) for an oscillation cycle, where *C*<sub>ext</sub> = 10 nF, 5 nF, 1 nF, 800 pF, 500 pF and 400 pF. The use of *t/period*<sub>osc</sub> 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*<sub>osc</sub> for an oscillation cycle. The *T* versus *t/period*<sub>osc</sub> (and *ρ* versus *t/period*<sub>osc</sub>) curves show dependency on frequency for *C*<sub>ext</sub> ≲ 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_{\text{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 $\rho$ becomes equal to $\rho_{\text{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_{\text{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 $\rho$ 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_{\text{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_{\text{ext}}$ up to
<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) *T* versus *t/period*<sub>osc</sub> (b) *ρ* vs *t/period*<sub>osc</sub> and (c) *ρ* versus *T* for an oscillation cycle. *C*<sub>ext</sub> = 10 nF–400 pF. Red crosses represent the points at the beginning of *C*<sub>ext</sub> discharging, after MIT. From panels (a) and (c) it can be seen that these points are associated with *dT/dt* sign flipping during LRS-to-HRS transition. From panels (b) and (c) it can be seen that the sign flipping occurs before the HRS heating branch is recovered. This makes accessible intermediate HRS values of resistivities, until the HRS heating branch is intercepted (gray circles) and the quasi-static HRS resistivity is recovered. If the HRS heating branch is not intercepted before HRS-to-LRS transition, the intermediate values of resistivities determine the IMT points, making them frequency-dependent.</figcaption></figure>

$C_{\text{ext}} = 1 \text{ nF. For } C_{\text{ext}} = 1 \text{ nF, } \rho \text{ at the beginning of discharging has the same value of } \rho_{\text{IMT}} = \rho(T_{\text{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_{\text{DD}} - V_{\text{out}}$ differs from the one of the quasi-static $I-V_D$. Finally, for $C_{\text{ext}} \le 800 \text{ pF}$, the value of $\rho$ at the beginning of discharging is smaller than $\rho_{\text{IMT}}$. Thus, the HRS heating branch is never recovered, and not both the HRS part of $I-V_D = V_{\text{DD}} - V_{\text{out}}$ and also the IMT point start to depend on $C_{\text{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].
<figure><img src="image_8.png" /><figcaption>Figure 5. The two components of equation (3), P/C<sub>th</sub> and -(T - T<sub>0</sub>)(A<sup>h</sup>)/C<sub>th</sub>, plotted against V<sub>D</sub> = V<sub>DD</sub> - V<sub>out</sub>. Red ellipses highlight the LRS part of the oscillation cycle where the increase in absolute value of -(T - T<sub>0</sub>)(A<sup>h</sup>) is not correlated to an increase in P/C<sub>th</sub>, and thus it can only be explained as being due to a reduction in the efficiency of thermal dissipation in this range of frequencies.</figcaption></figure>

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:

$$C_{\text{th}} \frac{dT}{dt} = P - (T - T_0)(A\bar{h}), \quad (3)$$

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 $\bar{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 \times 10^{-14}$ cm<sup>3</sup>) we calculate $C_{th} = 2.4 \times 10^{-13}$ J K<sup>-1</sup>. We extract $A\bar{h}$ by following the approach used in [28]. We match the quasi-static simulated *T* to $T_0 + P/(A\bar{h})$, where *P* is the quasi-static simulated power. We obtain that $(A\bar{h}) = 1.0753 \times 10^{-6}$ W K<sup>-1</sup> (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)(A\bar{h})/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)(A\bar{h})/C_{th}$ begins to depend on frequency for $C_{ext} \le 1$ nF. In particular, the absolute value of $-(T - T_0)(A\bar{h})/C_{th}$ increases with frequency for each value of $V_D = V_{DD} - V_{out}$ when it is $\lesssim 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_D^2/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)(A\bar{h})/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}/(A\bar{h}) = 2.23 \times 10^{-7}$ s, which thus becomes comparable to the period of charging ($1.28 \times 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} \le 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} \le 1$ nF. This is consistent with the above observation that the HRS branch of $I-V_D$ varies with $C_{ext}$ for $C_{ext} \le 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} \le 1$ nF) while the value of $I_{IMT}$ increases with it ($C_{ext} \le 800$ pF), as a consequence, 'flatter'
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Curves of *I* versus *V*<sub>D</sub> = *V*<sub>DD</sub> − *V*<sub>out</sub> (symbols) as extracted from simulated oscillations. Oscillatory behavior at *C*<sub>ext</sub> = 500 pF corresponds to load-line specified by *V*<sub>DD</sub> = 2 V, *R*<sub>ext</sub> = 15 kΩ (black dashed line). Oscillatory behavior at *C*<sub>ext</sub> = 400 and 350 pF is achieved with a load-line specified by *V*<sub>DD</sub> = 10 V, *R*<sub>ext</sub> = 80 kΩ (red line). Oscillatory behavior at *C*<sub>ext</sub> = 320 pF is achieved with a load-line specified by *V*<sub>DD</sub> = 30 V, *R*<sub>ext</sub> = 240 kΩ (blue line). For comparison, quasi-static *I* versus *V*<sub>D</sub> (thickazure line) is also shown.</figcaption></figure>

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<sub>ext</sub> = 400 and 350 pF is achieved with a load-line specified by V<sub>DD</sub> = 10 V, R<sub>ext</sub> = 80 kΩ. The corresponding oscillation frequencies are 294 kHz (C<sub>ext</sub> = 400 pF) and 291 kHz (C<sub>ext</sub> = 350 pF). Oscillatory behavior at C<sub>ext</sub> = 320 pF is achieved with a load-line specified by V<sub>DD</sub> = 30 V, R<sub>ext</sub> = 240 kΩ, with a corresponding oscillation frequency of 283 kHz. Thus, the load-line parameters V<sub>DD</sub> and R<sub>ext</sub> 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<sub>ext</sub> = 300 pF and below.

**4. 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₂ 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ₐ₆₉. 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 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₂ 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 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.

**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₂ 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₂ devices as memristive one-ports

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:

$$V = R(T)I, \qquad (A.1)$$

and the heat transfer equation,

$$\frac{dT}{dt} = -\frac{\delta}{C}(T - T_0) - \frac{P(t)}{C}, \qquad (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,

$$\frac{dT}{dt} = -\frac{\delta}{C}(T - T_0) - \frac{R(T)I^2}{C} \equiv f(T, I), \qquad (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 = −∇φ = ρ(T)J and equation (2), respectively, with the functional form of ρ(T) given in appendix B.

Appendix B. TCAD resistivity model of VO₂ and calibration procedure

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 ρ̂(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
<table><caption>**Table B1.** Analytical expressions of heating/cooling branches of TCAD $\hat{\rho}(T)$.</caption><thead><tr><th>ρ̂(T)</th><th>Condition</th></tr></thead><tbody><tr><td>ρ̂(T<sub>HRS</sub>) + m<sub>HRS</sub> × (T - T<sub>HRS</sub>)</td><td>if T < T<sub>HRS</sub></td></tr><tr><td>ρ̂(T<sub>LRS</sub>) + m<sub>LRS</sub> × (T - T<sub>LRS</sub>)</td><td>if T > T<sub>LRS</sub></td></tr><tr><td>ρ̂(T<sub>HRS</sub>) + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">(T - T<sub>HRS</sub>) × (ρ̂(T<sub>LRS</sub>) - ρ̂(T<sub>HRS</sub>))</div><div style="padding: 0 5px;">T<sub>LRS</sub> - T<sub>HRS</sub></div></div></td><td>Otherwise</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure B1.** (a) Experimental resistance *R* versus temperature *T* data used for calibration of the resistivity model. (b) Calibrated TCAD $\hat{\rho}(T)$. Red solid line is the heating branch of the curve, corresponding to the variation of $\hat{\rho}(T)$ when *T* is monotonically increasing with time. Dotted blue line is the cooling branch of the curve, giving the variation of $\hat{\rho}(T)$ when *T* is monotonically decreasing with time. HRS/LRS parts of heating/cooling branches are highlighted.</figcaption></figure>

<table><caption>**Table B2.** Calibrated parameters of TCAD $\hat{\rho}(T)$.</caption><thead><tr><th>Branch</th><th>m<sub>HRS</sub><br>(μΩcm K)<sup>-1</sup></th><th>T<sub>HRS</sub><br>(K)</th><th>ρ̂(T<sub>HRS</sub>)<br>(μΩcm)</th><th>m<sub>LRS</sub><br>(μΩcm K)<sup>-1</sup></th><th>T<sub>LRS</sub><br>(K)</th><th>ρ̂(T<sub>LRS</sub>)<br>(μΩcm)</th></tr></thead><tbody><tr><td>Heating</td><td>-9.81 × 10<sup>3</sup></td><td>327.6</td><td>5.543 × 10<sup>5</sup></td><td>0</td><td>340.9</td><td>1.33 × 10<sup>5</sup></td></tr><tr><td>Cooling</td><td>-9.81 × 10<sup>3</sup></td><td>304.6</td><td>5.543 × 10<sup>5</sup></td><td>0</td><td>320</td><td>1.33 × 10<sup>5</sup></td></tr></tbody></table>

HRS (for the cooling branch). The parameters $m_{HRS}$, $T_{HRS}$, $\hat{\rho}(T_{HRS})$, $m_{LRS}$, $T_{LRS}$, $\hat{\rho}(T_{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.

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_{Ceff}}{L}$, where *L* is the 'channel' length (that is the VO₂ thickness), and $A_{Ceff} = 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.

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})$.

## Appendix C. Simulation of quasi-static I-V<sub>D</sub> and DC I-V<sub>D</sub>

The quasi-static I-V<sub>D</sub> (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<sub>2</sub> 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<sub>D</sub>) 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<sub>D</sub>/max(V<sub>D</sub>) versus time for any bias applied. It can be observed the change of V<sub>D</sub> 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<sub>D</sub> is reported in figure C2(b) (star symbols). It can be seen that the quasi-static I-V<sub>D</sub> and the DC I-V<sub>D</sub> overlap. Thus, the NDR region determined from the quasi-static I-V<sub>D</sub> falls into the rigorous definition of the NDR region of a volatile memristor (see for instance [36]).

<figure><img src="image_11.png" /><figcaption>**Figure C1.** I-V<sub>D</sub> 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.</figcaption></figure>

A final remark. We have derived resistance values $R_{HRS} = 4.0763$ kΩ and $R_{LRS} = 722.5434$ Ω from the quasi-static I-V<sub>D</sub> of VO<sub>2</sub> memristor. We have used these values to approximate the memristor as a two-state variable resistor in the purely circuitual model of the VO<sub>2</sub> oscillator, which we have applied to calculate the theoretical frequencies (red asterisks) of figure 2(b) (appendix D).
<figure><img src="image_12.png" /></figure>

**Figure C2.** (a) $V_D$/max($V_D$) versus time for any bias applied. Change of $V_D$ with time is observed, until an equilibrium state is achieved. In the inset, the equilibrium time versus the applied bias is plotted, where the equilibrium time is the earliest time at which the $V_D$ values stabilize to fixed values. (b) DC $I-V_D$ (star symbols). Quasi-static $I-V_D$, obtained by sweeping an applied voltage of 0.5 V along a ramp time of 1.875 ms onwards, and sweeping backwards to 0 V in the same time interval, is plotted for comparison. It can be seen that the quasi-static $I-V_D$ and the DC $I-V_D$ overlap.

<figure><img src="image_13.png" /></figure>

**Figure D1.** (a) and (b) are the two possible oscillator circuit topologies. In both cases, the blue part of the $V_{out}$ waveform corresponds to cooling of the VO₂ channel (VO₂ is ON), the red part corresponds to heating of the VO₂ channel (VO₂ is OFF). (c) Equivalent Thevenin circuit.

## Appendix D. Frequencies calculated by theoretical circuit model of the VO₂ oscillator

Theoretical frequencies shown in figure 2(b) (red asterisks) are calculated through the purely circuitial 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_{IMT}$ and $V_{MIT}$, respectively. For generality, we consider the two possible topologies of the VO₂ oscillator circuit, where $C_{ext}$ is connected across $R_{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_{thev} = \frac{R_{VO2}R_{ext}}{R_{VO2}+R_{ext}}$, while the equivalent Thevenin voltages are slightly different:

* $V_{thev} = \frac{R_{ext}V_{DD}}{R_{VO2}+R_{ext}}$ (circuit topology of figure D1(a))

* $V'_{thev} = \frac{R_{VO2}V_{DD}}{R_{VO2}+R_{ext}}$ (circuit topology of figure D1(b))

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

It should be observed that $C_{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:
$$V_{\text{out}}(t) = V_{\text{final}} - (V_{\text{final}} - V_{\text{initial}}) \times \exp\left(-\frac{t}{R_{\text{they}}C_{\text{ext}}}\right), \quad (\text{D.1})$$

$$ \text{with } V_{\text{final}} = V_{\text{they}}, V_{\text{initial}} \text{ is equal to } \begin{cases} V_L = \min_{t \in \tau} V_{\text{out}} & \text{if } C_{\text{ext}} \text{ is charging} \\ V_H = \max_{t \in \tau} V_{\text{out}} & \text{if } C_{\text{ext}} \text{ is discharging} \end{cases} \quad \text{for } \tau \text{ oscillation period, where} \\ \begin{cases} V_L = V_{\text{DD}} - V_{\text{IMT}} \text{ and } V_H = V_{\text{DD}} - V_{\text{MIT}} \text{ for circuit topology of figure D1(a)} \\ V_L = V_{\text{MIT}} \text{ and } V_H = V_{\text{IMT}} \text{ for circuit topology of figure D1(b)} \end{cases} $$

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

$$t = \left[ \ln \left( \frac{V_{\text{initial}} - V_{\text{final}}}{V_{\text{out}}(t) - V_{\text{final}}} \right) \right] R_{\text{they}} C_{\text{ext}}, \quad (\text{D.2})$$

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

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

$$\tau_{\text{ch}} = \left[ \ln \left( \frac{V_L - V_{\text{they}}(\text{LRS})}{V_H - V_{\text{they}}(\text{LRS})} \right) \right] R_{\text{they}}(\text{LRS}) C_{\text{ext}}, \quad (\text{D.3})$$

$$= \left[ \ln \left( \frac{V_{\text{DD}} - V_{\text{IMT}} - \frac{R_{\text{ext}}V_{\text{DD}}}{R_{\text{LRS}} + R_{\text{ext}}}}{V_{\text{DD}} - V_{\text{MIT}} - \frac{R_{\text{ext}}V_{\text{DD}}}{R_{\text{LRS}} + R_{\text{ext}}}} \right) \right] \frac{R_{\text{LRS}} R_{\text{ext}}}{R_{\text{LRS}} + R_{\text{ext}}} C_{\text{ext}}. \quad (\text{D.4})$$

$$\tau_{\text{disch}} = \left[ \ln \left( \frac{V_H - V_{\text{they}}(\text{HRS})}{V_L - V_{\text{they}}(\text{HRS})} \right) \right] R_{\text{they}}(\text{HRS}) C_{\text{ext}}, \quad (\text{D.5})$$

$$= \left[ \ln \left( \frac{V_{\text{DD}} - V_{\text{MIT}} - \frac{R_{\text{ext}}V_{\text{DD}}}{R_{\text{HRS}} + R_{\text{ext}}}}{V_{\text{DD}} - V_{\text{IMT}} - \frac{R_{\text{ext}}V_{\text{DD}}}{R_{\text{HRS}} + R_{\text{ext}}}} \right) \right] \frac{R_{\text{HRS}} R_{\text{ext}}}{R_{\text{HRS}} + R_{\text{ext}}} C_{\text{ext}}. \quad (\text{D.6})$$

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

For the circuit topology of figure D1(b):

$$\tau_{\text{ch}} = \left[ \ln \left( \frac{V_L - V_{\text{they}}(\text{HRS})}{V_H - V_{\text{they}}(\text{HRS})} \right) \right] R_{\text{they}}(\text{HRS}) C_{\text{ext}}, \quad (\text{D.7})$$

$$= \left[ \ln \left( \frac{V_{\text{MIT}} - \frac{R_{\text{HRS}}V_{\text{DD}}}{R_{\text{HRS}} + R_{\text{ext}}}}{V_{\text{IMT}} - \frac{R_{\text{HRS}}V_{\text{DD}}}{R_{\text{HRS}} + R_{\text{ext}}}} \right) \right] \frac{R_{\text{HRS}} R_{\text{ext}}}{R_{\text{HRS}} + R_{\text{ext}}} C_{\text{ext}}. \quad (\text{D.8})$$

$$\tau_{\text{disch}} = \left[ \ln \left( \frac{V_H - V_{\text{they}}(\text{LRS})}{V_L - V_{\text{they}}(\text{LRS})} \right) \right] R_{\text{they}}(\text{LRS}) C_{\text{ext}}, \quad (\text{D.9})$$

$$= \left[ \ln \left( \frac{V_{\text{IMT}} - \frac{R_{\text{LRS}}V_{\text{DD}}}{R_{\text{LRS}} + R_{\text{ext}}}}{V_{\text{IMT}} - \frac{R_{\text{LRS}}V_{\text{DD}}}{R_{\text{LRS}} + R_{\text{ext}}}} \right) \right] \frac{R_{\text{LRS}} R_{\text{ext}}}{R_{\text{LRS}} + R_{\text{ext}}} C_{\text{ext}}. \quad (\text{D.10})$$

$$\text{Thus, } \tau = \beta' C_{\text{ext}} \text{ for } \beta' = \ln \left( \frac{V_{\text{IMT}} - \frac{R_{\text{HRS}}V_{\text{DD}}}{R_{\text{HRS}} + R_{\text{ext}}}}{V_{\text{IMT}} - \frac{R_{\text{HRS}}V_{\text{DD}}}{R_{\text{HRS}} + R_{\text{ext}}}} \right) \frac{R_{\text{HRS}} R_{\text{ext}}}{R_{\text{HRS}} + R_{\text{ext}}} + \ln \left( \frac{V_{\text{IMT}} - \frac{R_{\text{LRS}}V_{\text{DD}}}{R_{\text{LRS}} + R_{\text{ext}}}}{V_{\text{IMT}} - \frac{R_{\text{LRS}}V_{\text{DD}}}{R_{\text{LRS}} + R_{\text{ext}}}} \right) \frac{R_{\text{LRS}} R_{\text{ext}}}{R_{\text{LRS}} + R_{\text{ext}}}. $$

The theoretical frequencies plotted in figure 2(b) have been extracted from the periods $\tau = \beta C_{\text{ext}}$ of VO$_2$ 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_{\text{they}}(\text{HRS}) = 1.572\,632\,678\,8559$ V, $V_{\text{they}}(\text{LRS}) = 1.908\,088\,235\,2941$ V, $V_L = V_{\text{DD}} - V_{\text{IMT}} = 1.712\,915\,98$ V, $V_H = V_{\text{DD}} - V_{\text{MIT}} = 1.881\,303\,54$ V. We have calculated the theoretical frequencies extracted from the periods $\tau = \beta' C_{\text{ext}}$ of the VO$_2$ oscillator with the circuit topology drawn in figure D1(b). Indeed, we have found that $\beta = 3896.762\,624\,63s / F = \beta'$, which means that the two circuit topologies yield the same period for any given $C_{\text{ext}}$ value. The parameters used to perform these latter calculations are listed in table D1, with the following: $V'_{\text{they}}(\text{HRS}) = 0.427\,367\,321\,144$ V, $V'_{\text{they}}(\text{LRS}) = 0.091\,911\,764\,7059$ V, $V_L = V_{\text{MIT}}$, $V_H = V_{\text{IMT}}$.
<table><caption>**Table D1.** Parameters used to calculate theoretical frequencies.</caption><thead><tr><th>R<sub>HRS</sub></th><th>R<sub>LRS</sub></th><th>V<sub>IMT</sub></th><th>V<sub>MIT</sub></th><th>R<sub>they</sub>(HRS)</th><th>R<sub>they</sub>(LRS)</th></tr></thead><tbody><tr><td>4.0763 kΩ</td><td>722.5434 Ω</td><td>0.287 084 02 V</td><td>0.118 696 46 V</td><td>3.205 254 908 581 kΩ</td><td>689.338 235 2941 Ω</td></tr></tbody></table>

<table><caption>**Table D2.** Frequencies calculated by theoretical circuit model.</caption><tbody><tr><td>C<sub>ext</sub> (pF)</td><td>10<sup>4</sup></td><td>5 × 10<sup>3</sup></td><td>10<sup>3</sup></td><td>800</td><td>600</td><td>500</td><td>400</td><td>300</td></tr><tr><td>Frequency (kHz)</td><td>25.7</td><td>51.3</td><td>256.6</td><td>320.8</td><td>427.7</td><td>513.2</td><td>641.6</td><td>855.4</td></tr><tr><td>C<sub>ext</sub> (pF)</td><td>200</td><td>100</td><td>10</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Frequency (kHz)</td><td>1.283 × 10<sup>3</sup></td><td>2.566 × 10<sup>3</sup></td><td>2.566 × 10<sup>4</sup></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>**Table D3.** TCAD Simulated Frequencies.</caption><tbody><tr><td>C<sub>ext</sub> (pF)</td><td>10<sup>4</sup></td><td>5 × 10<sup>3</sup></td><td>10<sup>3</sup></td><td>800</td><td>600</td><td>500</td></tr><tr><td>Frequency (kHz)</td><td>21.6</td><td>41.678</td><td>157.9</td><td>182.935</td><td>210.2</td><td>221.4</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Figure E1.** TCAD-simulated local temperatures probed in the geometrical center of the device (solid lines) for $T_0$ equal to 293 and 313 K, and temperatures calculated from $T = T_0 + P/(A\bar{h})$ (symbols).</figcaption></figure>

## Appendix E. Determination of $A\bar{h}$

To determine $A\bar{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\bar{h})$, for $P$ the device power. $A$ is the heat transfer surface area and $\bar{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_D$. The use of $T = T_0 + P/(A\bar{h})$ is equivalent to modeling the heating effects in the $VO_2$ device with the theory for thermistors [40]. We achieve an excellent agreement (figure E1) by $A\bar{h} = 1.0753 \times 10^{-6} \text{ W K}^{-1}$.
## Appendix F. Simulated and calculated $dT/dt$

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

**Figure F1.** Temporal derivative of *T* determined by curve *T* versus *t* from simulated oscillations, and calculated from equation (3) of the manuscript, for *C*<sub>ext</sub> equal to (a) 10 nF, (b) 5 nF, (c) 1 nF, (d) 800 pF and (e) 500 pF. In each case, the match between the two curves is good overall, especially at lower frequencies. Discrepancies start to appear at higher frequencies (*C*<sub>ext</sub> below or equal 1 nF), especially in the LRS part of the oscillations.

## ORCID iDs

S Carapezzi https://orcid.org/0000-0002-9271-1189

A Plews https://orcid.org/0000-0002-4979-6754

G Boschetto https://orcid.org/0000-0001-8830-3572

A Nejim https://orcid.org/0000-0001-7656-5874

S Karg https://orcid.org/0000-0001-9873-2717

A Todri-Sanial https://orcid.org/0000-0001-8573-2910

## References

[1] Rana A, Li C, Koster G and Hilgenkamp H 2020 Resistive switching studies in VO₂ thin films *Sci. Rep.* 10 3293

[2] Zimmers A, Aigouy L, Mortier M, Sharoni A, Wang S, West K G, Ramirez J G and Schuller I K 2013 Role of thermal heating on the voltage induced insulator-metal transition in VO₂ *Phys. Rev. Lett.* 110 056601
[3] Kumar S, Pickett M D, Strachan J P, Gibson G, Nishi Y and Williams R S 2013 Local temperature redistribution and structural transition during joule-heating-driven conductance switching in VO₂ *Adv. Mater.* **25** 6128–32

[4] Li D, Sharma A A, Gala D K, Shukla N, Paik H, Datta S, Schlom D G, Bain J A and Skowronski M 2016 Joule heating-induced metal-insulator transition in epitaxial VO₂/TiO₂ devices *ACS Appl. Mater. Interfaces* **8** 12908–14

[5] Bortnikov S G, Aliev V S, Badmaeva I A and Mzhelskiy I V 2018 VO₂ film temperature dynamics at low-frequency current self-oscillations *J. Appl. Phys.* **123** 075701

[6] Driscoll T, Kim H-T, Chae B-G, Di Ventra M and Basov D N 2009 Phase-transition driven memristive system *Appl. Phys. Lett.* **95** 043503

[7] Driscoll T, Quinn J, Di Ventra M, Basov D N, Seo G, Lee Y-W, Kim H-T and Smith D R 2012 Current oscillations in vanadium dioxide: evidence for electrically triggered percolation avalanches *Phys. Rev. B* **86** 094203

[8] Wang R, Yang J-Q, Mao J-Y, Wang Z-P, Wu S, Zhou M, Chen T, Zhou Y and Han S-T 2020 Recent advances of volatile memristors: devices, mechanisms and applications *Adv. Intell. Syst.* **2** 2000055

[9] Zhou G et al 2022 Volatile and nonvolatile memristive devices for neuromorphic computing *Adv. Electron. Mater.* **8** 2101127

[10] Shukla N et al 2014 Synchronized charge oscillations in correlated electron systems *Sci. Rep.* **4** 4964

[11] Corti E, Cornejo Jimenez J A, Niang K M, Robertson J, Moselund K E, Gotsmann B, Ionescu A M and Karg S 2021 Coupled VO₂ oscillators circuit as analog first layer filter in convolutional neural networks *Front. Neurosci.* **15** 628254

[12] Todri-Sanial A et al 2022 How frequency injection locking can train oscillatory neural networks to compute in phase *IEEE Trans. Neural Netw. Learn. Syst.* **33** 1996–2009

[13] Jain A K, Mao J and Mohiuddin K M 1996 Artificial neural networks: a tutorial *Computer* **29** 31–44

[14] Paugam-Moisy H and Bohte S 2012 *Computing With Spiking Neuron Networks* (Springer) pp 335–76

[15] Hoppensteadt F C and Izhikevich E M 1999 Oscillatory neurocomputers with dynamic connectivity *Phys. Rev. Lett.* **82** 2983–6

[16] Hoppensteadt F C and Izhikevich E M 2012 *Weakly Connected Neural Networks* (Applied Mathematical Sciences) (Springer)

[17] Pikovsky A, Kurths J and Rosenblum M 2001 *Synchronization. A universal Concept in Nonlinear Sciences* (Cambridge University Press)

[18] Dörfler F and Bullo F 2014 Synchronization in complex networks of phase oscillators: a survey *Automatica* **50** 1539–64

[19] Engel A K, Fries P and Singer W 2001 Dynamic predictions: oscillations and synchrony in top-down processing *Nat. Rev. Neurosci.* **2** 704–16

[20] Shamsi J, Avedillo M J, Linares-Barranco B and Serrano-Gotarredona T 2021 Hardware implementation of differential oscillatory neural networks using VO₂-based oscillators and Memristor-Bridge circuits *Front. Neurosci.* **15** 674567

[21] Hoppensteadt F C and Izhikevich E M 2000 Pattern recognition via synchronization in phase-locked loop neural networks *IEEE Trans. Neural Netw.* **11** 734–8

[22] Corti E, Khanna A, Niang K, Robertson J, Moselund K E, Gotsmann B, Datta S and Karg S 2020 Time-delay encoded image recognition in a network of resistively coupled VO₂ on Si oscillators *IEEE Electron Device Lett.* **41** 629–32

[23] Delacour C, Carapezzi S, Abernot M and Todri-Sanial A 2023 Energy-performance assessment of oscillatory neural networks based on VO₂ devices for future edge AI computing *IEEE Trans. Neural Netw. Learn. Syst.* **1**–14

[24] Kamala V and Goldak J A 1993 Error due to two dimensional approximation in heat transfer analysis of welds *Weld. J.* (available at: www.osti.gov/biblio/6100681)

[25] 2021 Victory Mesh User Manual version 1.8.2.R Silvaco Inc.

[26] Carapezzi S, Boschetto G, Delacour C, Corti E, Plews A, Nejim A, Karg S and Todri-Sanial A 2021 Advanced design methods from materials and devices to circuits for brain-inspired oscillatory neural networks for edge computing *IEEE J. Emerg. Sel. Topics Power Electron.* **11** 586–96

[27] Carapezzi S, Boschetto G, Karg S and Todri-Sanial A 2022 Electro-thermal simulations of beyond-CMOS vanadium dioxide devices and oscillators *MRS Commun.* **12** 427–33

[28] Carapezzi S, Delacour C, Plews A, Nejim A, Karg S and Todri-Sanial A 2022 Role of ambient temperature in modulation of behavior of vanadium dioxide volatile memristors and oscillators for neuromorphic applications *Sci. Rep.* **12** 19377

[29] 2021 Victory Device User Manual version 1.19.1.C Silvaco Inc.

[30] Avrami M 1939 Kinetics of phase change I. General theory *J. Chem. Phys.* **7** 1103–12

[31] Avrami M 1940 Kinetics of phase change II. Transformation-time relations for random distribution of nuclei *J. Chem. Phys.* **8** 212–24

[32] Avrami M 1941 Granulation, phase change and microstructure kinetics of phase change. III *J. Chem. Phys.* **9** 177–84

[33] Lee S et al 2017 Anomalously low electronic thermal conductivity in metallic vanadium dioxide *Science* **355** 371–4

[34] Hamaoui G et al 2019 Thermophysical characterisation of VO₂ thin films hysteresis and its application in thermal rectification *Sci. Rep.* **9** 8728

[35] Maffezzoni P, Daniel L, Shukla N, Datta S and Raychowdhury A 2015 Modeling and simulation of vanadium dioxide relaxation oscillators *IEEE Trans. Circuits Syst.* **I** 62 2207–15

[36] Ascoli A, Demirkol A S, Tetzlaff R, Slesazeck S, Mikolajick T and Chua L O 2021 On local activity and edge of chaos in a NaMLab memristor *Front. Neurosci.* **15** 651452

[37] Lienhard J H IV and Lienhard J H V 2011 *A Heat Transfer Textbook* 4th edn (Phlogiston Press)

[38] Chua L O and Kang S-M 1976 Memristive devices and systems *Proc. IEEE* **64** 209–23

[39] Carapezzi S, Delacour C and Todri-Sanial A 2022 Simulation toolchain for neuromorphic oscillatory neural networks based on beyond-CMOS vanadium dioxide devices 2022 *IEEE Int. Conf. on Flexible and Printable Sensors and Systems (FLEPS)* pp 1–4

[40] Pandey R K 2019 *Fundamentals of Electroceramics. Materials, Devices and Applications* (Wiley)