# Role of ambient temperature in modulation of behavior of vanadium dioxide volatile memristors and oscillators for neuromorphic applications

Stefania Carapezzi, Corentin Delacour, Andrew Plews, Ahmed Nejim, Siegfried Karg, Aida Todri-Sanial

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

## ► To cite this version:

Stefania Carapezzi, Corentin Delacour, Andrew Plews, Ahmed Nejim, Siegfried Karg, et al. Role of ambient temperature in modulation of behavior of vanadium dioxide volatile memristors and oscillators for neuromorphic applications. Scientific Reports, 2022, 12 (1), pp.19377. 10.1038/s41598-022-23629-4. http://hal-03936310

HAL Id: hal-03936310

https://hal.science/hal-03936310

Submitted on 7 Jun 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

OPEN

# Role of ambient temperature in modulation of behavior of vanadium dioxide volatile memristors and oscillators for neuromorphic applications

Stefania Carapezzi<sup>1,✉</sup>, Corentin Delacour<sup>1</sup>, Andrew Plews<sup>3</sup>, Ahmed Nejim<sup>3</sup>, Siegfried Karg<sup>4</sup> & Aida Todri-Sanial<sup>1,2✉</sup>

Volatile memristors are versatile devices whose operating mechanism is based on an abrupt and volatile change of resistivity. This switching between high and low resistance states is at the base of cutting edge technological implementations such as neural/synaptic devices or random number generators. A detailed understanding of this operating mechanisms is essential prerequisite to exploit the full potentiality of volatile memristors. In this respect, multi-physics device simulations provide a powerful tool to single out material properties and device features that are the keys to achieve desired behaviors. In this paper, we perform 3D electrothermal simulations of volatile memristors based on vanadium dioxide (VO₂) to accurately investigate the interplay among Joule effect, heat dissipation and the external temperature T₀ over their resistive switching mechanism. In particular, we extract from our simulations a simplified model for the effect of T₀ over the negative differential resistance (NDR) region of such devices. The NDR of VO₂ devices is pivotal for building VO₂ oscillators, which have been recently shown to be essential elements of oscillatory neural networks (ONNs). ONNs are innovative neuromorphic circuits that harness oscillators' phases to compute. Our simulations quantify the impact of T₀ over figures of merit of VO₂ oscillator, such as frequency, voltage amplitude and average power per cycle. Our findings shed light over the interlinked thermal and electrical behavior of VO₂ volatile memristors and oscillators, and provide a roadmap for the development of ONN technology.

The Internet of Things (IoTs)—one of the leading technologies for the twenty first century—presents an unprecedented opportunity to provide ubiquitous computing, sensing and learning capabilities at the edge. Many smart applications in key fields such as healthcare, transport, autonomous cars, or smart cities, to name a few, are fast adopting IoTs with embedded artificial intelligence (AI) on them to perform a diverse set of applications. Presently, most of the deployed AI on edge devices are pre-trained or send data to the cloud or data centers¹. But with the continuous increase in edge devices and the need for real-time data processing, such a procedure is becoming unsuitable due to data latency, privacy, bandwidth saturation, and power consumption². Alternatively, an efficient way is to train AI algorithms on IoT devices by handling data locally where they are amassed. However, executing AI algorithms on current computing architecture based on von Neumann one suffers from large power consumption, mainly due to energy loss of data transfer between memory and processor. This clashes with the limited power and computing resources available on the edge. Therefore, novel concepts are required to go beyond von Neumann computing era to enable AI learning and inference at the edge.

Both academia and industry are exploring beyond von Neumann architectures, aiming to compute by merging memory and processing of information. While many alternatives are being scrutinized, analog computing based on the principle of "let physics compute" presents a quite interesting solution³. At its core, it transforms

¹Microelectronics Department, LIRMM, University of Montpellier, CNRS, Montpellier 34095, France. ²Department of Electrical Engineering, Eindhoven University of Technology, Eindhoven 5612, AP, The Netherlands. ³Silvaco Europe Ltd., St Ives PE27 5JL, UK. ⁴Department of Science and Technology, IBM Research Europe - Zurich, Ruschlikon 8803, Switzerland. ✉email: stefania.carapezzi@gmail.com; aida.todri@lirmm.fr
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Example of possible ONN application for image recognition tasks. (a) Input image. (b) Initialization of ONN circuit. Each pixel of input image is mapped as oscillator's phase $\Delta\Phi_{in,j}$ by delaying of $\Delta t_{in,j}$ the biasing $V_{DD,j}$ compared to the reference (bottom-most) oscillator. The computation begins at coupling of all the oscillators (c), where the simplest coupling elements are resistances $R_{C,i,j}$ (green blocks). (d) The oscillators synchronize in a stable in-phase ($\Delta\Phi_{out,j} = 0^\circ$) or out-of-phase ($\Delta\Phi_{out,j} \neq 0^\circ$) state. The phase of each oscillator is mapped back as pixel (gray scale) color to yield the (e) output image.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Temperature-triggered resistive switching of VO₂ volatile memristor. (b) 3D mesh of VO₂ CB device. (c) TCAD simulated $I-V_D$ characteristics. The associated NDR region has been highlighted. (d) Scheme of VO₂ oscillator circuit. Self-oscillatory electrical behavior is associated to a suitable choice of the load-line. (e) ONN.</figcaption></figure>

the dynamics of complex systems such as coupled oscillatory neural networks (ONNs)²,³ into an energy-efficient computing tool. ONNs take inspiration from the human brain, which is able to perform sophisticated tasks with a reduced energy budget. Specifically, ONNs aim to emulate biological neural networks, which are made of elemental units outputting periodic signals (the neurons) and interconnected into systems (the networks)⁴⁻⁸. ONN circuits mimic neural nets as ensembles of coupled oscillators (Fig. 1), where the information is encoded in the phase difference among the oscillators (phase-based computing).

To compute with ONNs, oscillators must be initialized with different delay times compared to a reference oscillator (Fig. 1b). This establishes an initial phase difference between each oscillator and the reference one. Then, at a given time, all oscillators are connected to each other through synaptic coupling elements (Fig. 1c), where the simplest coupling elements are resistors or capacitors. These couplings mediate the ONN interaction and allow to achieve a collective phase synchronization, if any, which is the result of ONN computing (Fig. 1d)⁹.

An example of application on ONN is pattern recognition, given that the phase state of each oscillator can be mapped to a white ($\Delta\Phi_j = 0^\circ$), a black ($\Delta\Phi_j = 180^\circ$) or a gray (all other values of $\Delta\Phi_j$) pixel (Fig. 1a,e). Overall, ONNs have attracted immense interest because they have a simple structure yet fast computation and energy efficient computing paradigm¹⁰⁻¹³.

Oscillators are the basic building blocks of ONNs. In this paper, we report on the realization of beyond CMOS oscillators based on vanadium dioxide (VO₂) volatile memristors. Volatile memristors are currently widely investigated for applications as neural/synaptic devices, or as random generators for hardware security¹⁴,¹⁵. Their operation mode relies on a mechanism of volatile resistive switching. VO₂ is a Transition Metal Oxide that changes from a high-resistive monoclinic (M₁) crystal structure to a low-resistive tetragonal rutile-like (R) structure. Such abrupt resistive switching is related to temperature overcoming a threshold (340 K) close to room temperature. This process takes the name of insulator-to-metal transition (IMT), it is reversible (metal-to-insulator transition, MIT) and often hysteretic (Fig. 2a). The overall variation of VO₂ resistivity across IMT/MIT can be up to several orders of magnitude.
The resistive switching of VO₂ has been experimentally observed also in channels of two-terminal VO₂ devices due to current injection or voltage biasing. This gives rise to the so called volatile Mott memristors¹⁴⁻¹⁸ because, upon removal of electrical stimulus, these devices reset back spontaneously to the high resistive state (HRS) from the low resistive state (LRS). Experimental evidence¹⁹⁻²² points towards self-heating as driving this transition by increasing the device temperature.

A negative differential resistance (NDR) region in the current-voltage (I-V) characteristics (Fig. 2c) of VO₂ volatile memristors is the hallmark of such resistive switching. This provides the possibility to build VO₂ oscillators by simply inserting a VO₂ volatile memristor into an RC circuit (Fig. 2d)¹⁷,²²,²³, without need for inductive components. The oscillatory regime is achieved provided that the circuit load-line $I_L = (V_{DD} - V)/R_{ext}$ crosses the NDR region of current I vs. voltage across the VO₂ ($V_D$) curve, where $V_{DD}$ is the applied bias and $R_{ext}$ is the external resistance. Thus, the periodic oscillation is due to cyclic charging and discharging of external capacitor $C_{ext}$ related to VO₂ device being in LRS and HRS (circuit topology of Fig. 2d).

Due to the role played by temperature to regulate the LRS/HRS states of VO₂ volatile memristors, it is expected that external temperature $T_0$, as provided by an external heater or the ambient, modulates the behavior of both VO₂ device and oscillator. Consequently, the external temperature will also impact ONN dynamics, affecting ONNs as analog computing architecture. For instance, the process of finding the ground-state of Ising machines can be implemented by gradually lowering the effective temperature (annealing)²⁴,²⁵, where such ground-state corresponds to the solution of optimization problems. Recently, an Ising Hamiltonian solver based on VO₂ oscillators has been demonstrated²⁶.

In this work, we assess the role of the external temperature over the electrical behavior of VO₂ volatile memristors and oscillators. To this aim, we use technology Computer-Aided Design (TCAD) technique to run multi-physics simulations of VO₂ devices and oscillators. Physics-based multi-scale simulations are essential tools for investigating beyond CMOS devices, where the control of materials growth and micro/nanofabrication procedures is still to be fully mastered. Given the variability of nominally identical devices, the abstraction of trends in device behavior through experiments is quite challenging. On the contrary, physics-based simulations can precisely consider material variability and/or device geometry, which captures correlations and provides insights into device behavior. We perform 3D TCAD electrothermal simulations of VO₂ volatile memristors to accurately take into account Joule effect, heat dissipation and $T_0$. We show the dependence of I vs. $V_D$ curves from $T_0$, and in particular, of the points of occurrence of resistive switching and its resetting. We are able to provide a simplified model for the evolution of such points with $T_0$. This yields the effect of modulation with $T_0$ of associated NDR region. We assess the impact of such modulation over VO₂ oscillator, considering that the load-line of oscillatory circuit does not depend on $T_0$. It is important to observe that the sensitivity of VO₂ oscillator to $T_0$ has a number of consequences. It can bring to cross-talk effects 1 self-induced, due to the local thermal build-up of VO₂ oscillator, or 2) among nearby VO₂ oscillators during ONN operativity. However, it also *de facto* embeds into VO₂ neuron-mimicking devices a (thermal) sensory function. Thus, our work provides essential TCAD tools on one hand to ascertain the thermal cross-talk effects in VO₂ oscillators, and on the other to engineer the response characteristics of VO₂ sensory neuron devices. Finally, we determine how figures of merit of VO₂ oscillator, such as frequency, voltage amplitude and average power per cycle, vary according to $T_0$. Our findings shed light on the interlinked thermal and electrical behavior of VO₂ oscillators that are helpful for providing guidelines for the implementation of ONN technology.

## TCAD electrothermal simulations of VO₂ crossbar device

We perform 3D TCAD electrothermal simulations of VO₂ crossbar (CB) devices fabricated on a silicon platform. The experimental VO₂ device is realized on a silicon substrate covered by 1 µm-thick SiO₂ layer²⁷. The bottom Pt electrode is deposited after etching a trench in SiO₂ layer. Then, the 80-nm thick VO₂ layer is grown and subsequently shaped into a square of 5 µm-side. Finally, the top electrode is laid down to complete the CB geometry (Fig. 2b). Both electrodes have a width of 250 nm.

**Modeling of the VO₂ resistive switching.** We use Silvaco Victory Mesh²⁸ for 3D mesh rendering of the actual geometry of VO₂ device (Fig. 2b). This allows to emulate accurately 3D effects of thermal dissipation and the real geometry of current flow. Experimental VO₂ is a policrystalline, granular film due to the growth on SiO₂ substrate, with average size of grains of ≈ 50 nm²⁷. However, in simulations we treat VO₂ as homogeneous material, because the active device region (about where the two contacts overlap) has a cross-section area of 250 nm by 250 nm, which amounts to a few grains. Thus, to consider the device “channel” as uniform is a good approximation that helps to reduce the computational time and complexity. We perform 3D TCAD electrothermal simulations of VO₂ volatile memristor by Silvaco Victory Device²⁹. We model the VO₂ as a conductor material whose resistivity vs. temperature is described by a user-defined version of the Phase-Change Material (PCM) model³⁰,³¹ (Fig. 2a). In static simulations, the resistivity $\hat{\rho}$ is linear with the temperature T, with different slopes $m_{LT}$, $m_{HT}$ depending on VO₂ being insulating ($T < T_{LT}$) or metallic ($T > T_{HT}$), respectively. $\hat{\rho}(T)$ is also linear across the IMT/MIT region ($T_{LT} < T < T_{HT}$). Because experimental curves of resistance vs. temperature show hysteresis (Fig. 2a), two different sets of parameters are used for the heating/cooling cycle. In time-dependent simulations, there will be a time lag for a sufficiently fast temperature change for the phase change to occur. Therefore, the time-dependent resistivity is modeled as:

$$ \rho(t + \Delta t) = \rho(t) + (\hat{\rho}(T) - \rho(t)) \times \left[ 1 - \exp\left(-\frac{\Delta t}{\tau}\right) \right] \quad (1) $$

where $\rho(t)$ is the resistivity at previous time step t, $\rho(t + \Delta t)$ is the resistivity at new time step and $\tau$ is a lifetime. The model also takes into account incomplete heating/cooling cycles, when they are reversed before achieving
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Simulated I-V<sub>D</sub> at external temperature of T<sub>0</sub> = 293 K (solid dark cyan line) and T<sub>0</sub> = 313 K (pink dotted line). The IMT/MIT points for T<sub>0</sub> ranging from 293 K (the rightmost ones) to 313 K (the leftmost ones) are shown as black (IMT) and yellow (MIT) starred lines, respectively. The black line is the load line valid for the temperature interval [293–307] K. The inset shows the MIT points for T<sub>0</sub> = 307 K (on the right) and T<sub>0</sub> = 308 K (on the left). It is evident that the MIT point for T<sub>0</sub> = 307 K is above the load-line, while the MIT point for T<sub>0</sub> = 308 K is below it.</figcaption></figure>

the full transition to the metal (or insulator) state. For example, if a cooling cycle starts at temperature T<sub>B</sub> such that T<sub>LT,H</sub> < T<sub>B</sub> < T<sub>HT,H</sub>, before the heating cycle is completed, then the resistivity remains constant until the temperature T<sub>A</sub> is reached such that ρ(T<sub>A</sub>) belongs to the cooling branch (Fig. 2a). This reproduces recent experimental findings<sup>32</sup> on the behavior of VO<sub>2</sub> resistivity vs. temperature.

We calibrate the static PCM resistivity ρ̂<sub>PCM</sub> against experimental data of resistance R<sub>exp</sub> vs temperature T<sup>31</sup>. Given the CB geometry, the cross-section of the current flow is not limited just to the region of the overlap between the top and bottom electrodes<sup>31</sup>, thus we cannot derive directly ρ̂ from R<sub>exp</sub> × A<sub>C</sub>/L, with L the VO<sub>2</sub> thickness and A<sub>C</sub> the cross-section area of overlapping contacts. However, most of the current is indeed flowing therein<sup>31</sup>, thus, we simply assume that

$$ \rho_{\text{PCM}} = R_{\text{exp}} \times A_{\text{Ceff}}/L \quad (2) $$

for A<sub>Ceff</sub> = A<sub>C</sub> × a, a being a parameter. Because VO<sub>2</sub> is modeled as a conductor, the conduction equation is solved as J = E/ρ = -∇φ/ρ, where J is the current density, E is the electric field and φ is the electric potential. The Poisson equation and the heat flow equation

$$ C \frac{\partial T_L}{\partial t} = \nabla(K \nabla T_L) + \frac{(\nabla \phi)^2}{\rho} \quad (3) $$

are solved coupled self-consistently using the Finite Element Method (FEM). In Eq. (3) T<sub>L</sub> is the local lattice temperature, K is the thermal conductivity and C is the heat capacity per unit volume. The Joule effect is accounted through the heat generation rate ∂Q/∂t = ρJ<sup>2</sup> = (∇φ)<sup>2</sup>/ρ, where Q is the heat contributed by Joule effect. The VO<sub>2</sub> temperature is obtained by TCAD simulations through solution of Eq. (3). We use experimental current I vs. voltage across the VO<sub>2</sub> (V<sub>D</sub>) characteristics as a benchmark to validate our approach<sup>31</sup>. The experimental I - V<sub>D</sub> data have been collected at the external temperature of T<sub>0</sub> = 303 K as measured at the oxide surface close to the CB device (Fig. SI0). We emulate the experimental condition by setting a temperature of T<sub>0</sub> = 303 K at the bottom surface of the VO<sub>2</sub> layer that we consider to be the thermal contact of TCAD simulations. We use a value of 0.06 W/(cm K)<sup>33</sup> for the thermal conductivity of VO<sub>2</sub>, and of 3 J/(cm<sup>3</sup> K)<sup>34</sup> for its heat capacity per unit volume.

We obtain a good agreement with the experimental data<sup>31</sup> (Fig. SI1) provided that 1) we set a = 2.7 in Eq. (2). Additionally, 2) we need to calibrate the interface thermal resistance of thermal contact set at the bottom of the device to achieve the same values of voltage V<sub>IMT</sub> and current I<sub>IMT</sub> at which the resistive switching occurs in simulated and experimental I - V<sub>D</sub> curves. The interface thermal resistance of thermal contact accounts for heat dissipation through VO<sub>2</sub> / substrate interface, and it depends from material properties of VO<sub>2</sub> and SiO<sub>2</sub>, and also from interface roughness, for example. Such interface thermal resistance is set as thermal conductance of thermal contact, which together with T<sub>0</sub> represent the thermal boundary condition. For simplicity sake, we consider the VO<sub>2</sub> / substrate interface as the main source of heat dissipation given that it is reasonable to assume that the heat dissipation across VO<sub>2</sub> / air interface, or VO<sub>2</sub> / top contact interface, can be, in comparison, neglected.

It is important to highlight that we do not fit VO<sub>2</sub> material parameters to match I - V<sub>D</sub> curves, but we use physical meaningful values. This demonstrates that, overall, we are able to reproduce by simulations the salient features of the electrical behavior of VO<sub>2</sub> CB devices because we capture the physical mechanisms operating inside the device.

## Evolution of NDR region with external temperature

We perform 3D TCAD electrothermal simulations of VO<sub>2</sub> crossbar (CB) devices at different external temperatures T<sub>0</sub>. The simulated I - V<sub>D</sub> curves for T<sub>0</sub> = 293 K (dark cyan solid line) and T<sub>0</sub> = 313 K (pink dotted line) are plotted in Fig. 3, showing the modulation effect of T<sub>0</sub>. In particular, the values of voltage and current at which the resistive switching (V<sub>IMT</sub>, I<sub>IMT</sub>) and its reset (V<sub>MIT</sub>, I<sub>MIT</sub>) occur decrease with T<sub>0</sub>. This can be easily interpreted: by increasing the background temperature, less and less heat has to be generated by Joule effect to achieve the
<figure><img src="image_7.png" /><figcaption>Figure 4. (a) Local temperatures probed in the geometrical center of the device as extracted by TCAD simulations (dark cyan solid lines) for $T_0$ equal to 293 and 313 K. We achieve a very good agreement with temperatures calculated after Eq. (4) for $\alpha = 5.15 \times 10^5$ K/W (pink dashed lines). (b) $\Delta T_{IMT}$ (red triangles) and $\Delta T_{MIT}$ (blue triangles) as determined by TCAD simulations for different $T_0$. The dotted lines are the correspondent fits by 2nd order polynomials.</figcaption></figure>

temperature for resistive switching inside the volatile memristor. We denote ($V_{IMT}$, $I_{IMT}$) as IMT point and ($V_{MIT}$, $I_{MIT}$) as MIT point. The black (IMT) and yellow (MIT) starred lines of Fig. 3 are the TCAD simulated IMT/MIT points for $T_0$ ranging from 293 K to 313 K. The IMT and MIT points define the edges of the NDR region, which then, varies with $T_0$. The NDR region defines the operative regime of the oscillator circuit (Fig. 2d): it yields electrical oscillations only if the circuit load-line $I_L = (V_{DD} - V)/R_{ext}$ crosses it.

We are able to model the dependence of IMT and MIT points with $T_0$ as obtained by 3D TCAD electrothermal simulations. Thus, we can predict the NDR region at any given $T_0$. In the following, we show the model for IMT points, the same applies to MIT points. We approximate the device as a homogeneous resistor, whose resistance is determined by $R = \hat{\rho}_{PCM}(T_{ch}) \times L/A_{Ceff}$. $T_{ch}$ is not the local temperature inside the device but a “global” temperature such that R is the device resistance at the given time. We assume $T_{ch}$ as the uniform temperature in the region of overlapping contacts if the current flux is confined therein. We correlate $T_{ch}$ to the device power as:

$$T_{ch} = T_0 + \alpha \times I \times V \quad (4)$$

where the electric power is converted into heat energy that dissipates to the ambient obeying Newton's law of cooling. This corresponds to describing the heating effects in VO₂ with the same theory used for thermistors³⁵. The validity of such approach has been experimentally demonstrated in the case of planar devices¹⁹. $\alpha$ is $(A \times \bar{h})^{-1}$, for A the heat transfer surface area, $\bar{h}$ the heat transfer coefficient averaged over A³⁶, and it is related to how the VO₂ channel exchanges heat with the ambient, thus it depends on material properties and device geometry¹⁹. We consider that TCAD simulated temperature probed in the geometrical center of the device is a good approximation of $T_{ch}$, because most of the current flow is concentrated in the region of overlapping contacts³¹. Thus, we match the TCAD temperature to Eq. (4). We achieve a very good agreement by setting $\alpha = 5.15 \times 10^5$ K/W. Figure 4a shows the TCAD (solid lines) and the calculated temperatures (dashed lines) for $T_0$ equal to 293 and 313 K. Now, $T_{IMT} = T_0 + \alpha \times I_{IMT} \times V_{IMT} = T_0 + \Delta T_{IMT}$ is the temperature of resistive switching. It is noteworthy to observe that $T_{IMT}$ depends on $V_{IMT}$ and $I_{IMT}$. In turn, $V_{IMT}$ depends from $T_{IMT}$, as well as $I_{IMT}$ depends from $T_{IMT}$. This means that there is not a closed-form expression of the relationship between $T_{IMT}$ and $V_{IMT}$, $I_{IMT}$. We use TCAD simulation to derive $T_{IMT}$ at different $T_0$. $\Delta T_{IMT}$ vs $T_0$ is plotted in Fig. 4b (red triangles). We got an excellent match by fitting these points with a 2nd-order polynomial (red dotted line). We infer $R(T_{IMT})$ from PCM model as $R(T_{IMT}) = \rho(T_{IMT}) \times L/A_{Ceff}$. Finally, $I_{IMT} = V_{IMT}/R(T_{IMT})$. However, from Eq. (4), $I_{IMT} = \Delta T_{IMT}/(\alpha \times V_{IMT})$. By combining the two above equations together, we determine

$$V_{IMT} = \sqrt{\frac{\Delta T_{IMT} \times R(T_{IMT})}{\alpha}} \quad (5)$$

and we get $I_{IMT}$ from $V_{IMT}$ and $R(T_{IMT})$. In Fig. 5a,b are plotted $I_{IMT}$ and $V_{IMT}$ against $T_0$, respectively, as extracted from TCAD simulations (red symbols) and as determined by our model (dashed red lines). We apply the same procedure also to $I_{MIT}$ and $V_{MIT}$ vs. $T_0$. The excellent agreement between simulated and calculated points in all the cases demonstrates that the behavior of IMT/MIT points with $T_0$ is correlated to the dependence of $T_{IMT}$, $T_{MIT}$
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) *I*<sub>IMT</sub> (red pentagons), *I*<sub>MIT</sub> (blue stars) and (b) *V*<sub>IMT</sub> (red circles), *V*<sub>MIT</sub> (blue diamonds) as derived from TCAD simulations. The dashed lines are the currents and voltages as calculated from the proposed model. In (a), the black solid line is drawn parallel to *T*<sub>0</sub> axis. It intersects the *I*<sub>IMT</sub> line at *T*<sub>0</sub> = 293 K, and the *I*<sub>MIT</sub> line at about *T*<sub>0</sub> = 312 K.</figcaption></figure>

from *T*₀, expressed by Eq. (4). Thus, at least for the resistive switching and its resetting, the local in-homogeneity
in resistivity as induced by the local temperature distribution can be ignored. The device can be reduced to a
homogeneous resistor, whose resistance is determined by the temperature of the center of the device.

Behavior of VO₂ oscillator with external temperature

VO₂ oscillator (Fig. 2d) is expected to depend as well as on T₀. For instance, figures of merit of the oscillator, such as the amplitude and period of the electrical signal at the output node, depend on V IMT and V MIT, and then from T₀. This can be easily understood in the following way. In first approximation, the triggering thermal mechanism induced by Joule effect can be ignored, and VO₂ volatile memristor is just described as a two-state (LRS and HRS) variable resistor, activated by threshold voltages V IMT and V MIT, respectively. Either LRS/HRS, the oscillator circuit equation is expressed by Kirchhoff’s law: C ext dV dt = I L – I 30,37 . At V IMT , the VO₂ is switched on (LRS) and the external capacitor C ext charges, with the voltage at the output node V out varying from V L = V DD – V IMT to V H = V DD – V MIT . At V MIT , the VO₂ is switched off (HRS) and C ext discharges, with V out varying from V H to V L. The cycles of capacitor charging and discharging are repeated periodically and yield the voltage oscillation. Then, the amplitude of voltage oscillations is V H – V L = V IMT – V MIT . Moreover, the times for capacitor charging τ ch and discharging τ dis both depend on V IMT and V MIT , as well as the oscillation period τ = τ ch + τ dis 37,38 . We investigate by 3D TCAD electrothermal simulation the behavior of VO₂ oscillator with T₀. It is worthwhile to remember that I) IMT/MIT points vary with T₀, while II) the circuit parameters, that is the load-line I L = (V DD – V)/R ext , do not depend from T₀. This produces a number of consequences.

* There is a finite range of external temperatures for which a fixed load-line crosses the associated NDR regions. That is, there exists a finite range of *T*₀ within which the VO₂ oscillator can oscillate.

For instance, in Fig. 3 is plotted the load-line (black solid line) corresponding to circuit parameters V<sub>DD</sub> = 3 V and R<sub>ext</sub> = 15 kΩ. It 1) passes above all the IMT points for T<sub>0</sub> > 293 K, and 2) is above MIT point at T<sub>0</sub> = 308 K and below MIT point at T<sub>0</sub> = 307 K (inset of Fig. 3). This means the load-line crosses the NDR regions in the interval of external temperatures [293–307] K. Figure 6 shows the simulated V<sub>out</sub> for the above load-line (C<sub>ext</sub> = 150 nF) at T<sub>0</sub> = 293 K (dark cyan solid line), T<sub>0</sub> = 303 K (dark cyan dashed dotted line), T<sub>0</sub> = 307 K (dark cyan dotted line) and T<sub>0</sub> = 308 K (dark cyan dashed line). As expected, the oscillatory behavior is present for all the T<sub>0</sub> except 308 K. It also follows that there will be intervals of external temperatures for which it is impossible to find a common load-line such that it crosses all the correspondent NDRs. For instance, in Fig. 5a, the (black solid) line is drawn parallel to T<sub>0</sub> axis, which intercepts the I<sub>IMT</sub> line at T<sub>0</sub> = 293 K and the I<sub>MIT</sub> line at about T<sub>0</sub> = 312 K. This means that the NDR regions for T<sub>0</sub> > 312 K cannot overlap with NDR region at T<sub>0</sub> = 293 K. Thus, it does not exist a load-line which crosses the NDR regions from T<sub>0</sub> = 293 K to above 312 K.

* There are temperatures for which the NDR region does not exist.

For example, in Fig. 5a the *I*<sub>IMT</sub> line crosses the *I*<sub>MIT</sub> line at about *T*₀ = 315 K. This means that for *T*₀ > 315 K the NDR region does not exist because *I*<sub>MIT</sub> < *I*<sub>IMT</sub>.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** TCAD simulated voltage oscillations for $V_{DD} = 3$ V, $R_{ext} = 15$ kΩ and $C_{ext} = 1.5 \times 10^{-7}$ F, at external temperature of $T_0 = 293$ K (solid line), $T_0 = 303$ K (dash-dotted line), $T_0 = 307$ K (dotted line) and $T_0 = 308$ K (dashed line).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a) Frequency of voltage oscillations vs the external temperature for $C_{ext} = 150$ nF (dark cyan triangles) and $C_{ext} = 40$ nF (pink stars). The curves obtained by the model of capacitor charging and discharging³⁷ (pink dashed and dark cyan solid lines) are also plotted. (b) TCAD simulated periods of charging $\tau_{ch}$ (downward triangles) and discharging $\tau_{dis}$ (upward triangles) for $C_{ext} = 150$ nF. The dotted and dashed lines are the curves for $\tau_{ch}$, $\tau_{dis}$ by applying the circuit equations³⁷.</figcaption></figure>

We perform TCAD electrothermal simulations of VO₂ oscillator at different external temperatures $T_0$. The oscillator circuit is simulated by considering the external resistance $R_{ext}$ and capacitance $C_{ext}$ as lumped element boundaries of VO₂ memristor. Any intrinsic device capacitance is accurately accounted into TCAD simulation of device electrostatics. We extract from simulated $V_{out}$ the frequency against $T_0$ in interval [293–307] K, for the load-line correspondent to $V_{DD} = 3$ V, $R_{ext} = 15$ kΩ. The results for $C_{ext} = 150$ nF (dark cyan lines) and $C_{ext} = 40$ nF (pink lines) are plotted in Fig. 7a. It can be appreciated the change of oscillator frequencies induced by variation of external temperature in the range [293–307] K. In particular, the frequency increases with $T_0$, it achieves a maximum, and then it decreases again at higher $T_0$. This behavior applies to all the probed values of the external capacitance. The decrease of frequency at $T_0 = 293$ K and $T_0 = 307$ K is easily interpreted with the fact that the load-line gets close to IMT point at $T_0 = 293$ K and to MIT point at $T_0 = 307$ K. Then, at these two temperatures, $I_L - I$ is very small. This, in turn, slows down the variation of voltage with time, as determined by Kirchhoff's law. This is also coherent with the fact that at $T_0 < 293$ K and $T_0 > 307$ K the load-line does not anymore cross the associated NDR regions. Thus, there is no more oscillatory behavior, which means zero frequency. We plot in Fig. 7b the periods of capacitor charging $\tau_{ch}$ (downward triangles) and discharging $\tau_{dis}$ (upward triangles) vs $T_0$ ([293–307] K). Their sum yields the oscillator period $\tau$. Figure SI2 of Supporting Information (SI) shows the correspondent curves for $C_{ext} = 40$ nF. Overall, TCAD simulated $\tau_{ch}$ and $\tau_{dis}$ have opposite trends with $T_0$ (belonging to [293–307] K): of increasing for $\tau_{ch}$ and of decreasing for $\tau_{dis}$. The percentage of variation between maximum and minimum values compared to the minimum one for TCAD simulated $\tau_{ch}$ varies between 113% ($C_{ext} = 150$ nF) and 87% ($C_{ext} = 40$ nF), while for TCAD simulated $\tau_{dis}$ it is comprised between 46% ($C_{ext} = 150$ nF) and 35% ($C_{ext} = 40$ nF). The percentage of variation between maximum and minimum values compared to the minimum one for TCAD simulated frequency is between 22% ($C_{ext} = 150$ nF) and 29% ($C_{ext} = 40$ nF). The
decreased variability of frequency with $T_0$ is explained with the partially compensating effects of the opposed trends of $\tau_{ch}$ and $\tau_{dis}$.

We complete the survey of effects induced by $T_0$ over the VO₂ oscillator. We derive the variation of voltage amplitude (Fig. SI6a of SI) and average power per cycle (Fig. SI6b of SI) against the external temperature (range [293–307] K) for the fixed load-line $V_{DD} = 3$ V, $R_{ext} = 15$ kΩ. $C_{ext}$ ranges from 150 nF to 40 nF. The voltage amplitude (Fig. SI6a) decreases with the external temperature, due to the trend of $V_{IMT}$ and $V_{MIT}$ with $T_0$ (Fig. 5b). As expected, there is no dependence from $C_{ext}$. Finally, we determine the average power per cycle <P>. We find that it lies in the range of values from $5.1 \times 10^{-4}$ W to $5.7 \times 10^{-4}$ W within the temperature interval [293 – 307] K and for all the external capacitance values.

## Conclusions

In this paper, we investigate the effect of ambient temperature over VO₂ volatile memristors and oscillators, using TCAD multi-physics simulations. We perform 3D TCAD electrothermal simulations of VO₂ volatile memristors to accurately take into account Joule effect, heat dissipation and $T_0$. Since our simulations link in a physical way VO₂ material properties (for example: resistivity vs temperature relationship) and oscillator circuit parameters (load-line), we are able to extract a number of consequences. The most important is to predict and explain the variation of oscillator frequency with the external temperature $T_0$. In particular, the variation of frequency is due to 1) the dynamics of variation of NDR region induced by $T_0$, and, on the contrary, 2) the stability against $T_0$ of load-line. It is important to observe that this behavior has been recently experimentally demonstrated for VO₂ memristors in self-oscillatory regime³⁹. This modulation effect of $T_0$ over behavior of VO₂ oscillators has multiple repercussions. It can bring to cross-talk effects (1) self-induced, due to the local thermal build-up of VO₂ oscillator, or (2) among nearby VO₂ oscillators during ONN operativity. In this last case, the amount of cross-talk will depend on a number of factors, among them the distance among oscillators, and thus it may ultimately set a limit to integration of VO₂ oscillators. Finally, the sensitivity of VO₂ oscillator to $T_0$ also *de facto* embeds into VO₂ neuron-mimicking device a (thermal) sensory function. Thus, our work provides essential TCAD tools on one hand to ascertain the thermal cross-talk effects in VO₂ oscillators, and on the other to engineer the response characteristics of VO₂ sensory neuron devices.

## Data availability

The datasets used and/or analyzed during the current study available from the corresponding author on reasonable request.

Received: 16 August 2022; Accepted: 2 November 2022

Published online: 12 November 2022

## References

1. Chin, J., Callaghan, V. & Allouch, S. B. The Internet-of-Things: Reflections on the past, present and future from a user-centered and smart environment perspective. *J. Ambient Intell. Smart Environ.* **11**, 45–69. https://doi.org/10.3233/AIS-180506 (2019).

2. Merenda, M., Porcaro, C. & Iero, D. Edge machine learning for AI-enabled IoT devices: A review. *Sensors*. https://www.mdpi.com/1424-8220/20/9/2533 (2020).

3. Raychowdhury, A. et al. Computing with networks of oscillatory dynamical systems. *Proc. IEEE* **107**, 73–89 (2019).

4. Hoppensteadt, F. & Izhikevich, E. *Weakly Connected Neural Networks. Applied Mathematical Sciences*. (Springer, 2012).

5. Pikovsky, A., Kurths, J. & Rosenblum, M. *Synchronization. A Universal Concept in Nonlinear Sciences* (Cambridge University Press, 2001).

6. Dörfler, F. & Bullo, F. Synchronization in complex networks of phase oscillators: A survey. *Automatica* **50**, 1539–1564. https://www.sciencedirect.com/science/article/pii/S0005109814001423 (2014).

7. Engel, A. K., Fries, P. & Singer, W. Dynamic predictions: Oscillations and synchrony in top-down processing. *Nat. Rev. Neurosci.* **2**, 704–716. https://doi.org/10.1038/35094565 (2001).

8. Shamsi, J., Avedillo, M. J., Linares-Barranco, B. & Serrano-Gotarredona, T. Hardware implementation of differential oscillatory neural networks using VO₂-based oscillators and memristor-bridge circuits. *Front. Neurosci*. https://www.frontiersin.org/articles/10.3389/fnins.2021.67456 (2021).

9. Csaba, G. & Porod, W. Coupled oscillators for computing: A review and perspective. *Appl. Phys. Rev.* **7**, 011302. https://doi.org/10.1063/1.5120412 (2020).

10. Eslahi, H., Hamilton, T. J. & Khandelwal, S. Energy-efficient ferroelectric field-effect transistor-based oscillators for neuromorphic system design. *IEEE J. Explor. Solid-State Comput. Devices Circuits* **6**, 122–129 (2020).

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

12. Delacour, C. et al. Oscillatory Neural Networks for Edge AI Computing. In *2021 IEEE Computer Society Annual Symposium on VLSI (ISVLSI)* 326–331 (2021).

13. Delacour, C., Carapezzi, S., Abernot, M. & Todri-Sanial, A. Energy-performance assessment of oscillatory neural networks based on VO₂ devices for future edge AI computing. https://hal.archives-ouvertes.fr/lirmm-03591176/

14. Wang, R. et al. Recent advances of volatile memristors: Devices, mechanisms, and applications. *Adv. Intell. Syst.* **2**, 2000055. https://doi.org/10.1002/aisy.20200055 (2020).

15. Zhou, G. et al. Volatile and nonvolatile memristive devices for neuromorphic computing. *Adv. Electron. Mater.* **8**, 2101127. https://doi.org/10.1002/aelm.202101127 (2022).

16. Driscoll, T., Kim, H.-T., Chae, B.-G., Di Ventra, M. & Basov, D. N. Phase-transition driven memristive system. *Appl. Phys. Lett*. https://doi.org/10.1063/1.3187531 (2009).

17. Driscoll, T. et al. Current oscillations in vanadium dioxide: Evidence for electrically triggered percolation avalanches. *Phys. Rev. B* **86**, 094203. https://doi.org/10.1103/PhysRevB.86.094203 (2012).

18. del Valle, J. et al. Electrically induced multiple metal-insulator transitions in oxide nanodevices. *Phys. Rev. Appl.* **8**, 054041. https://doi.org/10.1103/PhysRevApplied.8.054041 (2017).

19. Zimmers, A. et al. Role of thermal heating on the voltage induced insulator-metal transition in VO₂. *Phys. Rev. Lett.* **110**, 056601. https://doi.org/10.1103/PhysRevLett.110.056601 (2013).
20. Kumar, S. et al. Local temperature redistribution and structural transition during joule-heating-driven conductance switching in VO₂. *Adv. Mater.* **25**, 6128–6132. https://doi.org/10.1002/adma.201302046 (2013).

21. Li, D. et al. Joule Heating-Induced Metal-Insulator Transition in Epitaxial vo₂/tio₂ Devices. *ACS Appl. Mater. Interfaces* **8**, 12908–12914. https://doi.org/10.1021/acsami.6b03501 (2016).

22. Bortnikov, S. G., Aliev, V. S., Badmaeva, I. A. & Mzhelskiy, I. V. VO₂ film temperature dynamics at low-frequency current self-oscillations. *J. Appl. Phys.* **123**, 075701 (2018). https://doi.org/10.1063/1.5010971.

23. Shukla, N. et al. Synchronized charge oscillations in correlated electron systems. *Sci. Rep.* **4**, 4964. https://doi.org/10.1038/srep04964 (2014).

24. Mohseni, N., McMahon, P. L. & Byrnes, T. Ising machines as hardware solvers of combinatorial optimization problems (2022). *arXiv:2204.00276*.

25. Delacour, C. et al. VO₂-based Oscillatory Ising Machine: The Role of External Temperature on Performance. https://hal-lirmm.ccsd.cnrs.fr/lirmm-03725704.

26. Dutta, S. et al. An Ising Hamiltonian solver based on coupled stochastic phase-transition nano-oscillators. *Nat. Electron.* **4**, 502–512. https://doi.org/10.1038/s41928-021-00616-7 (2021).

27. Corti, E. et al. Coupled VO–2 Oscillators Circuit as Analog First Layer Filter in Convolutional Neural Networks. *Front. Neurosci*.https://www.frontiersin.org/articles/10.3389/fnins.2021.628254 (2021).

28. Victory mesh user manual. Version 1.8.2.R (Silvaco Inc, (2021).

29. Victory device user manual. Version 1.19.1.C (Silvaco Inc, 2021).

30. Carapezzi, S. et al. Multi-Scale Modeling and Simulation Flow for Oscillatory Neural Networks for Edge Computing. In *2021 19th IEEE International New Circuits and Systems Conference (NEWCAS)* 1–5 (2021).

31. Carapezzi, S. et al. Advanced design methods from materials and devices to circuits for brain-inspired oscillatory neural networks for edge computing. *IEEE J. Emerg. Sel. Top. Circuits Syst.* **11**, 586–596 (2021).

32. Rana, A., Li, C., Koster, G. & Hilgenkamp, H. Resistive switching studies in VO₂ thin films. *Sci. Rep.* **10**, 3293. https://doi.org/10.1038/s41598-020-60373-z (2020).

33. Lee, S. et al. Anomalously low electronic thermal conductivity in metallic vanadium dioxide. *Science* **355**, 371–374 (2017). https://www.science.org/doi/abs/10.1126/science.aag0410.

34. Hamaoui, G. et al. Thermophysical characterisation of VO₂ thin films hysteresis and its application in thermal rectification. *Sci. Rep.* **9**, 8728. https://doi.org/10.1038/s41598-019-45436-0 (2019).

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

36. Lienhard, J. H. I. & Lienhard, J. H. V. *A Heat Transfer Textbook* 4th edn. (Phlogiston Press, 2011).

37. Maffezzoni, P., Daniel, L., Shukla, N., Datta, S. & Raychowdhury, A. Modeling and simulation of vanadium dioxide relaxation oscillators. *IEEE Trans. Circuits Syst.* **62**, 2207–2215 (2015).

38. Carapezzi, S., Boschetto, G., Karg, S. & Todri-Sanial, A. Electro-thermal simulations of beyond-CMOS vanadium dioxide devices and oscillators. *MRS Commun*.https://doi.org/10.1557/s43579-022-00196-3 (2022).

39. Han, C. Y. et al. Characterization and modelling of flexible VO₂ Mott Memristor for the artificial spiking warm receptor. *Adv. Mater. Interfaces* **9**, 2200394. https://doi.org/10.1002/admi.202200394 (2022).

## Acknowledgements

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

## Author contributions

S.C. conceived the idea and designed the simulations. S.C. performed all the TCAD simulations and analyzed the data. S.C., C.D. 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.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-022-23629-4.

**Correspondence** and requests for materials should be addressed to S.C. or A.T.-S.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.