nanomaterials Review On the Thermal Models for Resistive Random Access Memory Circuit Simulation Juan B. Roldán 1, * , Gerardo González-Cordero 1 , Rodrigo Picos 2 , Enrique Miranda 3 , Félix Palumbo 4 , Francisco Jiménez-Molinos 1 , Enrique Moreno 5 , David Maldonado 1 , Santiago B. Baldomá 6 , Mohamad Moner Al Chawa 7 , Carol de Benito 2 , Stavros G. Stavrinides 8 , Jordi Suñé 3 and Leon O. Chua 9 1 Departamento de Electrónica y Tecnología de Computadores, Facultad de Ciencias, Universidad de Granada, Avd. Fuentenueva s/n, 18071 Granada, Spain; g2cordr@gmail.com (G.G.-C.); jmolinos@ugr.es (F.J.-M.); dmaldonado@ugr.es (D.M.) 2 Industrial Engineering and Construction Department, University of Balearic Islands, 07122 Palma, Spain; rodrigo.picos@uib.es (R.P.); carol.debenito@uib.es (C.d.B.) 3 Department Enginyeria Electrònica, Universitat Autònoma de Barcelona, Edifici Q., 08193 Bellaterra, Spain; enrique.miranda@uab.cat (E.M.); jordi.sune@uab.cat (J.S.) 4 Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Godoy Cruz 2290, Buenos Aires C1425FQB, Argentina; felix.palumbo@conicet.gov.ar 5 UJM-St-Etienne, CNRS, Laboratoire Hubert Curien UMR 5516, Institute of Optics Graduate School, University Lyon, F-42023 St-Etienne, France; enrique.manuel.moreno.perez@univ-st-etienne.fr 6 Unidad de Investigación y Desarrollo de las Ingenierías (UIDI), Facultad Regional Buenos Aires,  Universidad Tecnológica Nacional, Medrano 951, Buenos Aires C1179AAQ, Argentina; sanboyeras@gmail.com  7 Institute of Circuits and Systems, Technische Universität Dresden, 01062 Dresden, Germany; mohamad_moner.al_chawa@tu-dresden.de Citation: Roldán, J.B.; 8 School of Science and Technology, Thermi University Campus, International Hellenic University, González-Cordero, G.; Picos, R.; 57001 Thessaloniki, Greece; s.stavrinides@ihu.edu.gr Miranda, E.; Palumbo, F.; 9 Electrical Engineering and Computer Science Department, University of California, Jiménez-Molinos, F.; Moreno, E.; Berkeley, CA 94720-1770, USA; chua@berkeley.edu Maldonado, D.; Baldomá, S.B.; Moner * Correspondence: jroldan@ugr.es Al Chawa, M.; et al. On the Thermal Models for Resistive Random Access Abstract: Resistive Random Access Memories (RRAMs) are based on resistive switching (RS) oper- Memory Circuit Simulation. ation and exhibit a set of technological features that make them ideal candidates for applications Nanomaterials 2021, 11, 1261. https:// related to non-volatile memories, neuromorphic computing and hardware cryptography. For the full doi.org/10.3390/nano11051261 industrial development of these devices different simulation tools and compact models are needed in order to allow computer-aided design, both at the device and circuit levels. Most of the different Academic Editor: Antonio Di RRAM models presented so far in the literature deal with temperature effects since the physical Bartolomeo mechanisms behind RS are thermally activated; therefore, an exhaustive description of these effects is essential. As far as we know, no revision papers on thermal models have been published yet; and Received: 11 April 2021 Accepted: 1 May 2021 that is why we deal with this issue here. Using the heat equation as the starting point, we describe Published: 11 May 2021 the details of its numerical solution for a conventional RRAM structure and, later on, present models of different complexity to integrate thermal effects in complete compact models that account for the Publisher’s Note: MDPI stays neutral kinetics of the chemical reactions behind resistive switching and the current calculation. In particular, with regard to jurisdictional claims in we have accounted for different conductive filament geometries, operation regimes, filament lateral published maps and institutional affil- heat losses, the use of several temperatures to characterize each conductive filament, among other iations. issues. A 3D numerical solution of the heat equation within a complete RRAM simulator was also taken into account. A general memristor model is also formulated accounting for temperature as one of the state variables to describe electron device operation. In addition, to widen the view from different perspectives, we deal with a thermal model contextualized within the quantum point Copyright: © 2021 by the authors. contact formalism. In this manner, the temperature can be accounted for the description of quantum Licensee MDPI, Basel, Switzerland. effects in the RRAM charge transport mechanisms. Finally, the thermometry of conducting filaments This article is an open access article and the corresponding models considering different dielectric materials are tackled in depth. distributed under the terms and conditions of the Creative Commons Keywords: resistive memories; thermal model; heat equation; thermal conductivity; circuit simula- Attribution (CC BY) license (https:// tion; compact modeling; resistive switching; nanodevices creativecommons.org/licenses/by/ 4.0/). Nanomaterials 2021, 11, 1261. https://doi.org/10.3390/nano11051261 https://www.mdpi.com/journal/nanomaterials Nanomaterials 2021, 11, 1261 2 of 46 1. Introduction Resistive memories (also known as resistive random access memories or RRAMs) base their operation on resistive switching mechanisms to modulate their conductance in a non-volatile manner [1–3]. Their promising potential is subject to scrutiny both in academia and in industry. These devices could be alternatives to flash devices in applications related to their non-volatility, such as storage class memories [4]. RRAMs show short read/write times, good enough endurance and retention behavior, low power operation, CMOS technology compatibility and the possibility to be built on 3D stacks [1,2]. Apart from non- volatile memory circuits, where certain applications are currently on the market, RRAMs show great potential for cryptographic circuits due to their inherent stochasticity, which can be employed for the design of random number generators and physical unclonable functions [5–9]. Nevertheless, currently, the hottest application for these devices is linked to neuromorphic circuits [10–15]. RRAMs can mimic biological synapses within a fully compatible CMOS technology context to facilitate the fabrication of hardware neural networks. This approach allows the use of a lower number of electronic components (by means of RS device crossbars) and low power consumption [1,10,11,14,16] than the purely CMOS alternative to neuromorphic circuits firstly introduced by Mead in the late eighties [17]. The number of papers on this subject including resistive switching devices is growing by leaps and bounds. Nevertheless, different issues should be improved for these devices in order to overcome the difficulties connected to massive industrial use; among them, the following should be counted: great temperature sensitivity, manufacturing processes, variability and the lack of well-established electronic design automation (EDA) tools, including in the latter issue strategies for parameter extraction. In relation to these essential aspects, we will deal here with the device thermal description and its modeling from a circuit simulation perspective. The development of RRAM simulation and modeling infrastructure is key to step forward into industrial mature applications. In this respect, in the memory realm, it is important to highlight that DRAM and NAND Flash technologies will continue to hold their ground and even advance despite the slowing of Moore’s Law in the short-medium term. Although a great number of papers have been published so far on modeling and simulation issues [18–41], there are many open questions that need to be addressed to improve RRAM position in the EDA context. One of the pressing modeling questions is connected with the device inherent stochasticity that produces cycle-to-cycle variability [19,42–51]. This type of variability has to be managed to achieve RRAM technological maturity; however, even if a variability reduction is achieved, taking into consideration its nature in the modeling is a must [52]. Another modeling battlefield lies in what is linked to thermal effects. It is known that most RS mechanisms are thermally activated [19,23,32,53–57]. The temperature evolution produced by Joule heating in device-relevant regions, where charge conduction is concentrated, determines the operation in many different types of RRAMs [18,19,38,53,58–60]. Consequently, an accurate thermal description is essential to implement both good RRAM physical simulators and compact models. It is also important to highlight that in cross-bar architectures, an optimum topology for hardware neural network implementation among other applications, the thermal evolution of one device might be influenced by neighbor cells. This effect, known as thermal crosstalk [60–63], has to be considered at the integration stage of chips based on RS devices. Modeling of RRAM, as well as its physical simulation in general, shows notorious differences with respect to other electron devices. In this respect, in devices such as MOSFETs [64,65] (even multigate FETs [66]) or diodes [67,68], once a reasonable grid is established, the main differential equations are discretized and, in general, but for very scarce cases, convergence is searched for to obtain a solution. This is the main scheme to follow in drift-diffusion, hydrodynamic and Monte Carlo simulation approaches [69]. By contrast, in RRAMs the modeling paradigm is different due to the particularities of their operation. For the initial forming process (a pristine dielectric is assumed) and further set processes within RS cycles (in the common filamentary conduction regime, the case Nanomaterials 2021, 11, 1261 3 of 46 we will deal with henceforth), a conductive filament (CF) through the dielectric is formed which shorts the device electrodes and greatly reduce the device resistance [2,70]. Non- linear physical mechanisms (mostly thermally activated, following an Arrhenius’ equation relationship) come into play in a positive feedback loop that leads certain magnitudes, such as the temperature, to shoot up. This process has an obvious reflection at the experimental level, the current has to be limited to avoid the device hard breakdown and its consequent destruction. All this leads to the caveat that we have to deal with numerical divergence from the simulation viewpoint. The current is limited in the context of forming or set processes to avoid the device rupture; in this manner, an uncontrolled temperature rise and the corresponding CF quick growth is avoided because after a hard breakdown process RS operation does not hold. At the simulation level, a current limit is also needed. The reset process is also linked to RRAM simulation divergence. As the reset goes on, the CF shrinks and the same current is strangled in a CF narrow section. Hence, Joule heating produces a temperature increase at this CF narrowing. The thermally activated mechanisms in the CF narrowing surroundings rise exponentially and the CF rupture takes place. If the rupture is thermally controlled, as in many unipolar devices, it is based upon a thermal run-away process that leads to the CF destruction. In summary, for a correct RRAM simulation description (both at the device and circuit levels, in the latter case it would be a compact modeling approach), we have to face numerical divergence in addition to nonlinearity in some of the physical magnitudes at hand in one way or other. This means, in general, as the reader can imagine, a numerical nightmare. For instance, when you set up a limit for the compliance current, there might be an important difference between the iteration prior to achieving the current limit and the following step (once the compliance current has been exceeded) in terms of the CF size and temperature in the hottest spot. If you miss to stop the iterations just when the compliance current is reached, and for some reason you iterate once more within the positive feedback loop the device is going through, you can end up wondering what is the dielectric melting temperature. There are different RRAM models in the literature, and most of them consider thermal effects in one way or another. Some assume a simplified version of the steady-state heat equation (HE); others account for a non-steady-state model where the heat capacitance is included. In general, the main differences are found in the boundary conditions formulation to describe the device physics in the context of the HE or an equivalent energy balance equation. In any case, all these models can be formulated using the standard description proposed by [71], where they fit naturally. By using this formulation, it is clear that all these devices can be considered as extended memristors, this being the most general class of memristor devices according to usual classification [72]. Notice that, from this point of view, RRAM devices can be used inside the memristor framework as circuital elements for purposes further than memory applications [73–75]. Different flavors of these thermal models have been reported; however, as far as we know, they have not been brought together under one roof yet as has been done, for example, with electrical models [41]. We do so here. In Section 2, we comment on the HE in the RRAM operation context, in particular, we explain 3D physical simulation and compact modeling focusing on thermal effects. Section 3 is devoted to the RRAM thermal description through a general memristor model; the RRAM quantum point contact modeling including thermal effects is unfolded in Section 4, and finally, the thermometry of conducting filaments and corresponding modeling is developed in Section 5. 2. Mathematical Description of RRAM Thermal Effects 2.1. Heat Equation We start by taking the 3D heat equation in the device into consideration, Equation (1): . ∂T (r, t) ∇·(k th (r)∇ T (r, t)) + e generated (r) = ρ(r)c(r) (1) ∂t Nanomaterials 2021, 11, 1261 4 of 46 where kth —stands for the thermal conductivity (W/(m K)). This parameter depends on the temperature and geometry (assuming different material layers, which is the case for the usual RRAM architecture), c—stands for the specific heat or specific heat capacity (J/(kg K)). It is assumed that we are considering the specific heat capacity at constant pressure, which is why it is also denoted as cp , ρ—stands for the material density (kg/m3 ) . and e generated —stands for the power density generated (rate of heat generation by means of Joule heating, per unit volume inside the system we are considering). It can be calculated as σ(r)E2 (r), where σ(r) is the local electrical conductivity and E(r) is the local electric field [18,38]; i.e., the product of the field and the current density. The RRAM thermal description requires the solution of this equation in the whole device active region. Nevertheless, in the compact modeling approach (CMA) some simplifying assumptions are made. Among others, we consider this equation in the region close to the conductive filament, where charge conduction takes place after a successful set or forming process (when the CF is created the device is said to be in the low resistance state, LRS, while if the CF is ruptured, after a successful reset process, the device resistance is much higher, the device enters in the high resistance state, HRS). If all the different device material layers are included (dielectric, possibly a multilayer stack, electrodes, etc.), the thermal conductivity, the density and specific heat of the different materials have to be considered [29,76,77]. A step forward consists in using the 1D HE version (a simplifying assumption that works well in many cases). If the x coordinate is assumed to be parallel to the CF lon- gitudinal axis, from one of the dielectric-electrode interfaces to the other, and the CF is considered to be narrow enough to consider the same temperature in the transverse sections perpendicular to the x-axis; then, the HE could be written as follows:   ∂ ∂T ( x, t) . ∂T ( x, t) k (x) + e generated ( x ) = ρ( x )c( x ) (2) ∂x th ∂x ∂t This equation is most of the times solved within the device conductive filament, whose geometry is assumed to be a cylinder or a truncated-cone. For the HE particularization in the RRAM CMA some further assumptions are employed for the sake of simplicity: (a) Constant thermal conductivity, i.e., kth (x,T) = kth . Neither geometric nor temperature dependencies are considered. In most cases, the CF thermal conductivity is the one considered. (b) A single temperature in the whole conductive filament [38,78] is taken into account (this means a strong simplifying approach). Some models for circuit simulation can account for two different temperatures [79], this is a good strategy since the key (higher) temperature at the CF narrowing, where the CF is ruptured, is decoupled from the main CF bulk temperature; this latter temperature does not increase in the same manner. See Figure 4c in [80], where the CF temperature along the dielectric is plotted for different voltages. It is clear that the temperature is much higher in the CF narrowing while it shows a different behavior for the main CF body. The model with two different CF temperatures is more complex, hence, this issue has to be taken into account when dealing with circuits including hundreds or thousands of components. Different RRAM cell schematics are shown in Figure 1. Assuming filamentary con- duction, the CF evolution has to be calculated to describe RS operation and determine the device current. Several CF types (shapes) from the CMA employed in the literature are represented in Figure 1. Nanomaterials 2021, 11, 1261 5 of 47 Nanomaterials 2021, 11, 1261 device current. Several CF types (shapes) from the CMA employed in the literature are 5 of 46 represented in Figure 1. Figure 1. Different conductive filament shapes employed in RRAM compact models for circuit simulation. We will use Figure 1. Different conductive filament shapes employed in RRAM compact models for circuit simulation. We will use this CF shapes in the thermal models described below. (a) CF that occupies all the modeling domain, (b) cylindrical CF, this CF shapes in the thermal models described below. (a) CF that occupies all the modeling domain, (b) cylindrical CF, (c) (c) truncated-cone truncated-cone shaped shaped CF. CF. 2.2. A Numerical Approach for the Heat Equation 2.2. A Numerical Approach for the Heat Equation Equation (1) is essential to all types of RRAM physical simulators. It is usually auto- Equation consistently (1) is with solved essential other todifferential all types ofequations RRAM physical (Poisson simulators. equation, Itkinetic is usually auto- equations consistently solved with other differential equations (Poisson equation, for the chemical reactions that control the CF evolution, etc.) [21,23,34,53,70,77,81]. Since kinetic equations for the the most chemical common reactions device that controlisthe structure notCF evolution, based etc.) [21,23,34,53,70,77,81]. on curved surfaces or volumes, Since even the most common device structure is not based on in the more scaled cases, a finite difference approach could be a reasonable curved surfaces or volumes, even choice in that the more scaled cases, a finite difference approach could be simplifies the grid and the differential equations discretization. Boundary conditions a reasonable choice that sim- plifies are keythe togrid and the describe differential correctly equations the physics of discretization. the device operation.BoundaryInconditions are key many simplified to describe correctly the physics of the device operation. In models, the electrodes are supposed to be perfect heat sinks and, therefore, Dirichlet’s many simplified models, the electrodes boundary are supposed conditions aretoestablished be perfect heat sinks at the and, therefore, Dirichlet’s electrode-dielectric interface boundary (in most casescon- ditions are established room temperature at the at is assumed electrode-dielectric this interface). However, interfaceto (in most cases correctly room describe heattempera- transfer ture is assumed between at this interface). the conductive filament, However, to correctly the dielectric describe layer and heat transfer the electrodes, between some partstheof conductive filament, the dielectric layer and the electrodes, the latter should be included in the simulation domain (SD). The SD lateral interfaces some parts of the latter should be included in the (perpendicular to simulation domain (SD).interface) the dielectric-electrode The SD lateral interfaces are usually (perpendicular described by Neumann to the dielectric-electrode boundary conditions. interface) Sometimes are usually describedderivative the temperature by Neumann boundary in the conditions. normal direction of Sometimes these interfaces the temperature is assumed to derivative be null. This in the meansnormal direction adiabatic of theseaccordingly conditions interfaces iswith as- sumed Fourier’s to belawnull. This conduction. for heat means adiabatic conditions accordingly with Fourier’s law for heat conduction. We have included here some results to illustrate the HE role in a RRAM simulator. Our We have included simulator calculates here thesome RRAM results to illustrate current in the LRS. the The HE role in a resistance internal RRAM simulator. calcula- Our tion simulator is performed calculates assuming the RRAM fully formedcurrentmetallic-like in the LRS. The internal filaments conductive resistance ofcalculation different shapes. is performedIn addition assuming to the current fully formed calculation, we solve metallic-like the 3D HE conductive [29]. Inofthe filaments SD we different have included shapes. In additionthe dielectric stackcalculation, to the current and part of we thesolve electrodes; the 3Dprecisely, HE [29]. In10the nmSD of the Si-n+ we have (bottom electrode) included the dielectric andstack Ni (top and electrode). Consequently, part of the electrodes; the SD precisely, 10thermal nm of the boundary Si-n+ (bot- is shifted from the dielectric/electrode interfaces to the electrodes, tom electrode) and Ni (top electrode). Consequently, the SD thermal boundary is shifted as commented above. A 40 nmthe from axis) × 40 nm (Y axis) (X dielectric/electrode × 30 nm interfaces to(Z the axis, vertical axis electrodes, running from as commented the SiAlayer above. 40 nm to theaxis) (X Ni layer) × 40 nm SD (Yis considered, axis) × 30 nm see(ZFigure 2. axis, vertical axis running from the Si layer to the Ni layer)Dirichlet’s boundary SD is considered, seeconditions Figure 2. were supposed at the outer electrode layer surfaces. The constant temperature located outside the device was room temperature (a reasonable assumption accounting for the high electrode thermal conductivity). For the SD lateral faces we employed perfectly matched layers (PML) [82], an improved implementation of Neumann boundary conditions since PML are particularly appropriate for differential equation solution to deal with open boundary problems, our case here [83]. Joule heating takes place in the CFs. Some devices (schematics shown in Figure 2) were simulated. In our simulations we included two CFs to account for a different device configuration with respect to conventional studies with just one CF. Nanomaterials2021, Nanomaterials 2021,11, 11,1261 1261 66 of of4647 Figure2.2. Simulated Figure Simulated RRAMs. RRAMs.Four Fourdifferent differentLRS LRS situations situationsareare considered, considered, in in all all cases a double cases a double conductivefilament conductive filamentwaswasemployed employed assuming assuming a distance, a distance, d, between d, between them them (the(the bottom bottom andand top top figures figures are the same in each case from different perspectives). The bottom electrode is assumed to+ are the same in each case from different perspectives). The bottom electrode is assumed to be Si-n be Si-n+ and the top electrode is made of Ni, the dielectric consists of HfO2. The conductive fila- and the top electrode is made of Ni, the dielectric consists of HfO2 . The conductive filament shapes ment shapes employed are shown for the different cases under study, they are assumed to be me- employed tallic-like,are shownby formed forNithe different atom cases clusters under study, [24,55,80]. The they are assumed physical to beare parameters metallic-like, formed by the same employed Ni atom clusters [24,55,80]. in the simulation in [29]. The physical parameters are the same employed in the simulation in [29]. In Figure 3,boundary Dirichlet’s symmetricconditions temperatureweredistributions supposed ataretheobserved in (a) and outer electrode (b), layer corre- surfaces. sponding to similar CFs, and therefore, to alike Joule heating effects. See the effects of the The constant temperature located outside the device was room temperature (a reasonable low thermal conductivity in the dielectric. According to Fourier’s law for heat conduction, assumption accounting for the high electrode thermal conductivity). For the SD lateral the heat flux, q, is equal to the product of thermal conductivity, kth , and the negative faces we employed perfectly matched layers (PML) [82], an improved implementation of local temperature gradient (q = −k th ∇ T). Since the HfO2 thermal conductivity is around Neumann boundary conditions since PML are particularly appropriate for differential 1 W/(K m), the temperature drops off rapidly around the CF. See also the effects of the equation solution to deal with open boundary problems, our case here [83]. Joule heating dielectric-electrode boundary at z = 10 nm and z = 30 nm, the temperature reduction takes place in the CFs. Some devices (schematics shown in Figure 2) were simulated. In is different from what is seen in the CF perpendicular direction along the x-axis. The our simulations we included two CFs to account for a different device configuration with maximum temperature is obtained in the narrowest section for the symmetrical truncated- respect to conventional studies with just one CF. cone shaped CF, as seen in (c) and (d). At this point, the physical mechanisms behind RS In Figure 3, symmetric temperature distributions are observed in (a) and (b), corre- are more active and, consequently, they trigger the CF rupture at this location. See the sponding to similar CFs, and therefore, to alike Joule heating effects. See the effects of the thermal connection between the CFs for the different shapes and distances in between; low in thermal this conductivity respect, in the in tree-branch dielectric. shaped According filaments, to Fourier’sof the destruction lawthefor heat conduction, branches and the Nanomaterials 2021, 11, 1261 the heat flux, q, is equal to the product of thermal conductivity, k th, and the negative thermal and current redistribution in the remaining intertwined branches makes the7 reset local of 47 temperature gradient a complicated process. (𝒒 = −𝑘 ∇𝑇). Since the HfO 2 thermal conductivity is around 1 W/(K m), the temperature drops off rapidly around the CF. See also the effects of the die- lectric-electrode boundary at z = 10 nm and z = 30 nm, the temperature reduction is dif- ferent from what is seen in the CF perpendicular direction along the x-axis. The maximum temperature is obtained in the narrowest section for the symmetrical truncated-cone shaped CF, as seen in (c) and (d). At this point, the physical mechanisms behind RS are more active and, consequently, they trigger the CF rupture at this location. See the thermal connection between the CFs for the different shapes and distances in between; in this re- spect, in tree-branch shaped filaments, the destruction of the branches and the thermal and current redistribution in the remaining intertwined branches makes the reset a com- plicated process. Figure 3. Figure 3. Temperature Temperature plots plots for for some some ofof the the devices devices simulated simulated as as explained explained above above (see (see the the insets), insets), the the cross-section cuts cross-section cuts corresponds to y = 20 nm in our SD. (a) RRAM with two cylindrical CFs (diameter = 3 nm) separated 1 nm apart for corresponds to y = 20 nm in our SD. (a) RRAM with two cylindrical CFs (diameter = 3 nm) separated 1 nm apart for a bias of 0.6 a bias V. of 0.6 V. (b) RRAM with two cylindrical CFs (diameter = 3 nm) separated 5.5 nm apart for a bias of 0.7 V. (c) RRAM with one cylindrical CF (diameter = 6 nm) and a symmetrical truncated-cone shaped CF (low diameter = 3 nm, high diameter = 6 nm) separated 6 nm apart for a bias of 0.5 V. (d) RRAM with one cylindrical CF (diameter = 6 nm) and a symmetrical truncated-cone shaped CF (low diameter = 3 nm, high diameter = 6 nm) separated 6.5 nm apart for a bias of 0.8 V. For the Nanomaterials 2021, 11, 1261 7 of 46 Figure 3. Temperature plots for some of the devices simulated as explained above (see the insets), the cross-section cuts corresponds to y = 20 nm in our SD. (a) RRAM with two cylindrical CFs (diameter = 3 nm) separated 1 nm apart for a bias (b) RRAM with two cylindrical CFs (diameter = 3 nm) separated 5.5 nm apart for a bias of 0.7 V. (c) RRAM with one cylindrical of 0.6 V. (b) RRAM with two cylindrical CFs (diameter = 3 nm) separated 5.5 nm apart for a bias of 0.7 V. (c) RRAM with CF one(diameter = 6CF cylindrical nm) and a symmetrical (diameter = 6 nm) andtruncated-cone a symmetricalshaped CF (low diameter truncated-cone shaped CF = 3(low nm,diameter high diameter = 3 nm, = 6high nm) diameter separated= 66nm apart for a bias of 0.5 V. (d) RRAM with one cylindrical CF (diameter = 6 nm) and a symmetrical truncated-cone nm) separated 6 nm apart for a bias of 0.5 V. (d) RRAM with one cylindrical CF (diameter = 6 nm) and a symmetrical shaped CF (low diametershaped truncated-cone = 3 nm,CFhigh diameter (low = 6 =nm) diameter separated 3 nm, 6.5 nm apart high diameter for aseparated = 6 nm) bias of 0.86.5 V. For nm the sake apart forofa visibility, bias of 0.8 some of the V. For the 3D plots sake are rotatedsome of visibility, withofrespect the 3Dtoplots the 2D areCF scheme. rotated with respect to the 2D CF scheme. In In Figure Figure 44 simulations simulations ofof other otherRRAM RRAMconfigurations configurationsareareshown. shown. SeeSee that that small small changes in the device voltage can lead to important temperature variations changes in the device voltage can lead to important temperature variations in the CF bot-in the CF bottleneck where RS mechanisms are thermally enhanced. In (c) and (d), the thermal tleneck where RS mechanisms are thermally enhanced. In (c) and (d), the thermal cross- cross- talk talk between between CFsCFs is is observed observed when when thethe distance distance between between CFs CFsisisaround around11nm.nm.ForFormuch much higher inter CF distances, a lower thermal connection results even for higher currents. higher inter CF distances, a lower thermal connection results even for higher currents. Figure Figure 4. Temperature plots 4. Temperature plots for for some some of of the devices devices simulated simulated as as explained explained above above (see (see the the insets), insets), the the cross-section cross-sectioncuts cuts correspondstotoyy==2020nm corresponds nminin our our SD. SD. (a)(a) RRAMRRAM with with oneone cylindrical cylindrical CF (diameter CF (diameter = 3 nm) = 3 nm) and and a symmetrical a symmetrical truncated- truncated-cone cone shaped shaped CF (lowCFdiameter (low diameter = 3 nm,=high 3 nm, high diameter diameter = 6 nm) =separated 6 nm) separated 2 nm apart2 nm for apart a bias for a bias of 0.8 V. (b)ofRRAM 0.8 V. (b) withRRAM with one one cylindrical cylindrical CF (diameter = 3 nm) and a symmetrical truncated-cone shaped CF (low diameter = CF (diameter = 3 nm) and a symmetrical truncated-cone shaped CF (low diameter = 3 nm, high diameter = 6 nm) separated 3 nm, high diameter =6 nm) separated 6 nm apart for a bias of 0.5 V. (c) RRAM with one cylindrical CF (diameter = 3 nm) and a symmetrical 6 nm apart for a bias of 0.5 V. (c) RRAM with one cylindrical CF (diameter = 3 nm) and a symmetrical truncated-cone shaped truncated-cone shaped CF (low diameter = 5 nm, high diameter = 7 nm) separated 6 nm apart for a bias of 0.75 V. (d) CF (low diameter = 5 nm, high diameter = 7 nm) separated 6 nm apart for a bias of 0.75 V. (d) RRAM with one cylindrical CF RRAM with one cylindrical CF (diameter = 3 nm) and a symmetrical truncated-cone shaped CF (low diameter = 5 nm, (diameter = 3 nm) high diameter = 7and nm)a separated symmetrical truncated-cone 1 nm shaped apart for a bias CFV. of 0.6 (low Fordiameter the sake=of5 visibility, nm, high diameter some of the = 7 nm) separated 3D plots 1 nm are rotated apart for a bias of 0.6 V. For the with respect to the 2D CF scheme. sake of visibility, some of the 3D plots are rotated with respect to the 2D CF scheme. 2.3. Explicit Heat Equation Solutions 2.3.1. RRAM with a Cylindrical Filament (Steady-State Operation, No Heat Transfer Term) In the case under consideration, we assume a cylindrical CF with constant electrical and thermal conductivities. The boundary conditions are established at the extremes of the filament (x = 0 and x = tox , respectively, with temperatures T(x = 0) = T(x = tox ) = T0 , where tox stands for the dielectric thickness), see Figure 5. From Equation (2), we obtain [38]: ∂2 T ( x  σE2 = −k th (3) ∂x2 where E is the constant electric field in the CF (E = VRRAM /tox ). Figure 5 shows the different elements taken into account to solve the HE. The solution for the maximum temperature in the middle of the CF, with the boundary conditions sketched in Figure 5, is: 2.3.1. RRAM with a Cylindrical Filament (Steady-State Operation, No Heat Transfer Term) In the case under consideration, we assume a cylindrical CF with constant electrical and thermal conductivities. The boundary conditions are established at the extremes of Nanomaterials 2021, 11, 1261 the filament (x = 0 and x = tox, respectively, with temperatures T(x = 0) = T(x = tox) = T80,of 46 where tox stands for the dielectric thickness), see Figure 5. Figure 5. Schematic Figure Schematicshowing thethe showing different elements different considered elements to solve considered to the simplified solve HE in a cy- the simplified HE in a lindrical and cylindrical homogenous and homogenousconductive filament. conductive filament. From Equation (2), we obtain [38]: σ t2 E2 Tmax = 𝜕T0𝑇(𝑥)+ ox (4) 𝜎𝐸 = −𝑘 8 k th (3) 𝜕𝑥 The single temperature that will represent the whole CF thermal state is chosen to be Twhere max , asEshown below (see is the constant Appendix electric field inD.3 thein CF[84]), (E = VRRAM/tox). Figure 5 shows the different elements taken into account to solve the HE. The solution for the maximum temperature 2 σ VRRAM in the middle of the CF, with the boundary conditions sketched in Figure 5, is: Tmax = T0 + (5) 8 k th 𝜎𝑡 𝐸 𝑇 = 𝑇 + (4) Henceforth, this will be our thermal model 8𝑘 1, TM1. This assumption is justified by the fact that this maximum value controls the physical mechanisms that lead to the CF The single temperature that will represent the whole CF thermal state is chosen to be narrowing in a reset process and finally to the CF rupture. Nevertheless, among other Tmax, as shown below (see Appendix D.3 in [84]), issues in this model, the calculation of the ohmic resistance as the CF heats up is not correctly solved since the main CF body 𝜎 remains 𝑉 at a temperature much lower than the hottest spot where the maximum 𝑇 = 𝑇 + temperature is achieved. (5) The Verilog-A implementation 8𝑘 is shown in Table 1 (TM1). Henceforth, this will be our thermal model 1, TM1. This assumption is justified by 2.3.2. RRAM the fact withmaximum that this a Cylindricalvalue Filament controlsIncluding the physical a Heat Transfer that mechanisms Termlead to the CF (Steady-State narrowing in Operation) a reset process and finally to the CF rupture. Nevertheless, among other Nanomaterials 2021, 11, 1261 issuesWeinhave this model, added athe calculation term to account of the for ohmic the heat resistance transferas the the from CF heats CF to up theissurrounding not 9 ofcor- 47 rectly solved insulator since the (3), to Equation main CFFigure see body remains 6. The heatat a losses temperature much lower are included than the by means hot-heat of the test spotcoefficient transfer where the (h) maximum temperature [31,38,79,85], is achieved. see Equation (6): The Verilog-A implementation is shown in Table 1 (TM1). 2 𝜕 𝑇(𝑥)∂2 T ( x ) 𝑇(𝑥) −T𝑇( x ) − T0 𝜎 𝐸 =σ E 2.3.2. RRAM with a Cylindrical Filament−𝑘 = + 22 ℎ + 2ahHeat Transfer Term (Steady-State −k th Including (6) (6) 𝜕𝑥 ∂x 𝑟 rCF Operation) where whereWe rrCF stands stands CFhave for forthe added atheconductive termconductive to accountfilament for theradius. filament radius. heat transfer from the CF to the surrounding insulator to Equation (3), see Figure 6. The heat losses are included by means of the heat transfer coefficient (ℎ) [31,38,79,85], see Equation (6): Figure Figure6.6.Sketch Sketchofofa acylindrical cylindricalfilament with filament thethe with different terms different in the terms inHE theshown in Equation HE shown (6), (6), in Equation including the heat transfer term. including the heat transfer term. Equation (6) can be solved and the result is given below [27], 𝜎 𝐸 𝑟 (𝑒 − 1) 𝑇 = 𝑇 + (7) 2 ℎ (𝑒 + 1) where: Nanomaterials 2021, 11, 1261 9 of 46 Table 1. Verilog-A implementation of some of the thermal models described. Thermal Model Verilog-A Code for the Temperatures Calculation TM1 T = T0 + sigma*V**2/(8.0*kth); E = V/tox; TM2 alpha = tox/2.0*sqrt(2*h/(kth*rcf)); T = T0 + (sigma*E**2*rcf*(exp(alpha)−1)**2)/(2*h*(exp(2*alpha)+1)); LCF = tox−g; rg = sqrt(rt*rb); eta = rt/rb; TM3 E = V/LCF; alpha = LCF*sqrt(2*h/(kth*rg)); T = T0 + rg*sigma*E**2/(eta*h)*(0.5 − (exp(alpha/2.0)/(exp(alpha) + 1)); analog function real fdt0 ; real sigmat,rcf,E,alpha; input sigmat,rcf,E,alpha; begin fdt0 = sigmat*rcf*E**2*tanh(alpha*tox/2.0)/(sqrt(2*kth*h*rcf)); end endfunction analog function real fT; real sigmat,rcf,eta,alpha,dt0 ; input sigmat,rcf,eta,alpha,dt0 ; begin fT = T0 + sigmat*rcf*E**2/(2*h)*(1-cosh(alpha*tox/2.0)) + dt0 /alpha*sinh(alpha*tox/2.0); end endfunction TM4 analog function real falpha; real rcf; input rcf; begin falpha = sqrt(2*h/(kth*rcf)); end endfunction analog function real fsigmat; real T; input T; begin fsigmat = sigma/(1 + alphat * (T − T0 )); end endfunction Equation (6) can be solved and the result is given below [27], σ E2 rCF (eα − 1)2 Tmax = T0 + (7) 2 h (e2α + 1) where: s tox 2h α= (8) 2 k th rCF In this case, as before, the electric field is assumed constant, and the CF temperature is considered to be Tmax , as calculated in Equation (7). We will consider this the thermal model 2, TM2 (see the Verilog-A code in Table 1). Nanomaterials 2021, 11, 1261 10 of 46 2.3.3. RRAM with a Truncated-Cone Shaped Filament Including Heat Transfer Coefficient (Steady-State Operation) The boundary conditions and the CF geometry for this model are shown in Figure 7. Notice that the CF radius rCF ( x) of the truncated cone depends on the x position. Conse- quently, the heat equation can be expressed now as follows: 2h ∂2 T ( x ) ials 2021, 11, 1261 σ ( x ) E ( x )2 − ( T ( x ) − T0 ) = −k th 10 of 47 (9) rCF ( x ) ∂x2 a) b) Figure 7. (a) Energy dissipation terms included in the heat equation and geometrical domain for the CF thermal description, Figure 7. (a) Energy dissipation terms included in the heat equation and geometrical domain for (b) cylindrical CF equivalent employed to simplify the HE solution and obtain a compact analytical expression for the the CF thermal description, (b) cylindrical CF equivalent employed to simplify the HE solution CF temperature. Note that the conductive filament length (LCF ) can be lower than the oxide layer (tox ), due to the gap (g) and obtain a compact analytical expression for the CF temperature. Note that the conductive fila- between ment the top electrode length (LCF)and can the conductive be lower filament than the oxide tip (see layer (tox[27,78,86–88]). ), due to the gap (g) between the top elec- trode and the conductive filament tip (see [27,78,86–88]). This equation cannot be solved analytically because of the variable CF radius, rCF ( x). In order to obtain an approximated analytical solution, we considered a transformation Using this parameter: to simplify. A truncated-cone shaped CF with constant conductivity was found to be √ approximately equivalent in terms of this calculation to a cylinder with radius (r g = r T r B ) and variable conductivity 𝛼 = 𝐿 σ2ℎ CF ( x) [27], see Figure 7. For a fixed applied voltage, we include (11) the maximum electric field. 𝑟 value is affected by the ratio between the two truncated- 𝑘 This cone radii (η = r T /r B ). Under this assumption, the following simplified HE was obtained: The value of Tmax can be obtained as follows 2 (the one we assume for the whole CF, σ VRRAM this will be the thermal model 3, TM3, see the 2h code in Table 1), ∂2 T ( x ) Verilog-A − ( T ( x ) − T0 ) = −k th (10) L2CF η rg ∂x2 Using this 𝑟 𝜎𝐸 1 𝑒 𝑇 = 𝑇parameter: + − (12) 𝜂ℎ 2 𝑒 +1 s 2h α = LCF (11) k th r g The value of Tmax 2.3.4. RRAM with a Truncated-Cone can be Shaped obtained Filament as follows Including (the Heat one weCoeffi- Transfer assume for the whole CF, cient (Steady-State Operation) and Two Temperature Values to Represent the CFTable this will be the thermal model 3, TM3, see the Verilog-A code in Ther-1), mal Behavior " # r g σ E2 1 α e2 The use of two different CF temperatures Tmax = Ta0 more allows + accurate−thermal description (12) ηh 2 eα + 1 of the CF narrowing, where the temperature rises due to the thermal run-away process that leads to the RS operation, and a CF bulk temperature that accounts for the tempera- 2.3.4. RRAM with a Truncated-Cone Shaped Filament Including Heat Transfer Coefficient ture in the CF wider part. ThisOperation) (Steady-State CF regionand could TwobeTemperature left almost untouched in the sequence Values to Represent the CF Thermal Behavior of set and reset processes. We will consider this version the thermal model 4, TM4 (see in The use of two different CF temperatures allows a more accurate thermal description Table 1 a Verilog-A implementation of the first thermal models proposed here). In the of the CF narrowing, where the temperature rises due to the thermal run-away process that approach we followed, the electric field is calculated by means of two components asso- leads to the RS operation, and a CF bulk temperature that accounts for the temperature in ciated with the CF top and bottom parts [79]. In this respect, the electric field can be cal- ( , ) culated as: 𝐸( , ) = , where 𝑉( , ) is obtained considering the voltage divider / formed by the resistances associated with the CF top and bottom portions (see the Appen- dix in [79]). In the approximation strategy developed here, the maximum temperature (see the Nanomaterials 2021, 11, 1261 11 of 46 the CF wider part. This CF region could be left almost untouched in the sequence of set and reset processes. We will consider this version the thermal model 4, TM4 (see in Table 1 a Verilog-A implementation of the first thermal models proposed here). In the approach we followed, the electric field is calculated by means of two components associated with Nanomaterials 2021, 11, 1261 11 of 47 the CF top and bottom parts [79]. In this respect, the electric field can be calculated as: V ( T,B) E(T,B) = L/2 , where V(T,B) is obtained considering the voltage divider formed by the resistances associated with the CF top and bottom portions (see the Appendix in [79]). ) 𝑟 the 𝜎( , In 𝐸 ) ( , ) ( , approximation 𝛼( , ) 𝑡 𝑑𝑇 ( here, 𝛼( , ) 𝑡 1 − 𝑐𝑜𝑠ℎstrategy developed 𝑠𝑖𝑛ℎthe maximum temperature (see (13) the , ) 𝑇( , ) = 𝑇 + + previous2ℎthermal models) is obtained 2 𝛼( , cylinders for two ) of 2different radii (associated to the top and bottom CF temperatures), see Figure 8. These values are assumed to be the where, parameters 𝛼( , ) and 𝑑𝑇 ( , ) are given in the equations below [79]: temperatures at the main CF volumes linked to the cylinders: for the thicker CF section (Figure 8b), TT , and for the narrow CF (Figure 8c), TB . The analytical expression for the 2ℎ temperature calculation is given 𝛼( below: , ) = (14) 𝑘 𝑟 ( , ) σ(T,B) rCF(T,B) E(2T,B)   α(T,B) tox  dT 0(T,B)  α(T,B) tox  T(T,B) = T0 + 1 − cosh 𝛼( , ) 𝑡+ sinh (13) 2h 𝜎( , ) 𝑟 ( , ) 𝐸( , ) 𝑡𝑎𝑛ℎ2 2 2 α(T,B) 𝑑𝑇 ( , ) = (15) where, parameters α(T,B) and dT0(T,B) are2given 𝑘 ℎ 𝑟in the ( , ) equations below [79]: s Equation (16) includes the CF conductance temperature 2h dependence, assuming a me- tallic behavior: α(T,B) = (14) k th rCF ( T,B) 𝜎 = 𝜎( ( T,B) tox α  , ) 2 (16) 1 +) 𝛼CF(T,B σ(T,B r 𝑇() , ()T,B E −)𝑇tanh 2 dT0(T,B) = q (15) where 𝜎 stands for the CF conductivity at 2room k th htemperature rCF(T,B) (𝑇 ) and 𝛼 is the conduc- tivity temperature coefficient. The Verilog-A code is shown in Table 1 (TM4). a) b) c) Figure 8. RRAM cell CF scheme for the thermal model based on two different CF temperatures. (a) Original filament Figure 8. RRAM cell CF scheme for the thermal model based on two different CF temperatures. (a) Original filament with with the corresponding boundary conditions. Cylindrical CFs (shown in dashed lines) employed to compute the (b) top the corresponding boundary conditions. Cylindrical CFs (shown in dashed lines) employed to compute the (b) top tem- temperature TT and perature TT and (c) the (c) the TB. TB . bottom bottom Equation (16) includes the CF conductance temperature dependence, assuming a Table 1. Verilog-A implementation of some of the thermal models described. metallic behavior: σ Thermal Model Verilog-A Code for) the σ(T,B  0 Calculation = Temperatures  (16) TM1 1 + α T T(T,B) − T0 T = T0 + sigma*V**2/(8.0*kth); E = V/tox; TM2 where σ0 stands for the CFalpha conductivity at room temperature (T0 ) and α T is the conductivity = tox/2.0*sqrt(2*h/(kth*rcf)); temperature coefficient. The Verilog-A code is shown in Table 1 (TM4). T = T0 + (sigma*E**2*rcf*(exp(alpha)−1)**2)/(2*h*(exp(2*alpha)+1)); LCF = tox−g; rg = sqrt(rt*rb); eta = rt/rb; TM3 E = V/LCF; alpha = LCF*sqrt(2*h/(kth*rg)); T = T0 + rg*sigma*E**2/(eta*h)*(0.5 − (exp(alpha/2.0)/(exp(alpha) + 1)); Nanomaterials 2021, 11, 1261 12 of 46 2.4. Energy Balance in the Device If we apply the first law of thermodynamics in terms of an energy balance in the device active region, we would obtain Equation (17) [78]: . ∂T (t) e generated − κ ( T (t) − T0 ) = Cth (17) ∂t where T0 stands for the temperature of the dielectric that surrounds the CF (usually assumed at room temperature) and κ is the inverse of the thermal resistance. In this case we do not account for a temperature distribution along the CF as in Equation (2). Our perspective here accounts for the whole CF, and even its surroundings, which is represented by a single temperature. In addition, we do not follow, as above, a scheme based on the HE solution and the association of the maximum temperature with the CF. The power generated can be calculated as VRRAM (t)IRRAM (t), although we will employ V(t)I(t) for short. We can study this differential equation accounting for different operation regimes. 2.4.1. Steady-State This regime works well under the conventional ramped voltage stress, RVS. With long ramps to switch the device, Equation (17) can be written as follows: . e generated = κ ( T (t) − T0 ) (18) The power generated in the conductive filament (CF) (Joule heating effects calculated as current x voltage) equals the power dissipated toward the electrodes and the dielectric. Under the consideration of a single temperature to characterize the device active region, assuming that the electrodes and the dielectric are perfect heat sinks at a fixed temperature, Nanomaterials 2021, 11, 1261 T0 , Equation (18) can be written as follows: 13 of 47 T − T0 Rth = (19) VI 𝑇−𝑇 𝑅 = (19) where Rth stands for the effective 𝑉𝐼 thermal resistance (it depends on the device physical features and is associated with heat conduction [89]). Using this simple model (thermal where Rth stands for the effective thermal resistance (it depends on the device physical model 5, TM5), the device temperature can be estimated from Equation (19) [88,90–94]. features and is associated with heat conduction [89]). Using this simple model (thermal model At circuit 5, TM5), thesimulation level, Equation device temperature (19) can can be estimated frombe implemented Equation with the equivalent (19) [88,90–94]. electrical sub-circuit At circuit simulationshown in Figure level, Equation (19)9.canThe heat dissipated be implemented with power is represented by a the equivalent dependent electrical current sub-circuit source shown whose in Figure value 9. The heat isdissipated described power = V I, where as Gis1 represented by athe voltage is dependent determined current source by the twowhose inputvalue is described sub-circuit as 𝐺 = V terminals 𝑉𝐼, + where and Vthe− , voltage and theis current de- is sensed termined by the two input sub-circuit terminals 𝑉 and 𝑉 , and by a null voltage source Vsense connected in series between the input terminalsthe current is sensed by I + and a null − voltage source 𝑉 connected in series between the input terminals 𝐼 and 𝐼 . I . The thermal resistance is represented by an electrical resistance, Rth , and the room The thermal resistance is represented by an electrical resistance, 𝑅 , and the room tem- temperature by a constant perature by a constant voltage source, , with a T voltage𝑇source, 0 , with value thata equals value the thatroom equals the room temperature temperature (T ), in K. (𝑇 ),0 in 𝐾. The output sub-circuit terminal (𝑇) provides a voltage that represents therepresents The output sub-circuit terminal (T) provides a voltage that device the device temperature temperature 𝐾). K). T (in T (in .param T0= 300 ; K V+ T V- Rth I+ G1 Vsense VT0 0 {T0} I- value={(V(V+)-V(V-))*I(Vsense)} ; G1 value Figure Figure9. 9. Equivalent electric Equivalent circuit electric for thefor circuit RSthe device RS thermal model based device thermal on abased model thermal onresistance a thermal resistance Rth . Rth. 2.4.2. Non-Steady-State Approach Models based only on thermal resistances do not provide capacitive effects in terms of transient operation (mostly related to pulsed voltage stress, PVS, that can be employed to tune the device resistance in a multilevel operation regime, as needed in neuromorphic Nanomaterials 2021, 11, 1261 13 of 46 2.4.2. Non-Steady-State Approach Models based only on thermal resistances do not provide capacitive effects in terms of transient operation (mostly related to pulsed voltage stress, PVS, that can be employed to tune the device resistance in a multilevel operation regime, as needed in neuromorphic circuits). In order to include thermal inertia in the model, a capacitor (in the electrically equivalent thermal circuit) is added, Cth , the thermal or heat capacitance. This approach is used in different RRAMs compact models (thermal model 6, TM6), [78,95,96]. We can reformulate Equation (17) to the following expression: T − T0 dT VI = + Cth (20) Rth dt From Equation (20), the temperature can be obtained as: dT T = T0 + V I Rth − τth (21) dt where τth (thermal time constant) is defined as follows: τth = Cth Rth (22) Equation (20) can be solved analytically if we assume a constant voltage (see Equation (23)).   − t T (t) = T0 + V I Rth 1 − e τth (23) If we apply a time-dependent voltage, V(t), Equation (23) does not work. In this case, we have to solve numerically Equation (20). Assuming a new function, X(t) = T(t) − T0 , we obtain, X (t) dX (t) V (t) I (t) = + Cth (24) Rth dt Discretizing under equally distant temporal points ti+1 = ti + ∆t: Xi X − Xi Vi Ii = + Cth i+1 (25) Rth ∆t which gives us: Xi + 1 − Xi Rth Vi Ii = Xi + τth (26) ∆t and consequently, ∆t Xi + 1 = ( R V I − Xi ) + Xi (27) τth th i i Finally, the device temperature could be calculated as Ti+1 = Xi+1 + T0 . From a circuit simulation point of view, Equation (20) can be implemented with the equivalent electrical sub-circuit shown in Figure 10. The values of the different electric elements and the role of the pins is the same as in Figure 9. However, a capacitor has been added to account for the thermal capacitance. The values κ = 2.8–25 µJ/K s (Rth = 4 × 104 –3.5 × 105 K/W) and Cth = 0.04–1.1 pJ/K were employed in [78]. In [95], the following values were given: Cth = 0.318 fJ/K and τth = 0.23 ns, from them the thermal resistance can be extracted Rth = τth /Cth = 7.23 × 105 K/W. The ther- mal resistance in [91] was taken Rth = 5 × 105 K/W. An estimation of 33 ps for the thermal time constant is reported in [18], which has to be compared to the electric pulse-width to assess the importance of thermal transient effects in conventional RRAM operation. The heat capacitance can be calculated as Cth = Cp tox A [18], where A is the CF effective area, the effective CF length that approximately corresponds to the dielectric thickness (tox ) and and consequently, ∆𝑡 𝑋 = (𝑅 𝑉 𝐼 − 𝑋 ) + 𝑋 (27) 𝜏 Nanomaterials 2021, 11, 1261 Finally, the device temperature could be calculated as Ti+1 = Xi+1 + T0. From a circuit 14 of 46 simulation point of view, Equation (20) can be implemented with the equivalent electrical sub-circuit shown in Figure 10. The values of the different electric elements and the role ofCthe pins is the same as in Figure 9. However, a capacitor has been added to account for p is the volumetric heat capacity (calculated as ρ × c, Equation (2)) [89]. As detailed in [18], the thermal capacitance. the thermal resistance can be calculated as follows: .param T0= 300 ; K T V+ V- Rth Cth I+ G1 Vsense VT1 Nanomaterials 2021, 11, 1261 15 of 47 0 {T0} I- to assess the importance of thermal transient effects in conventional RRAM operation. The heat capacitance canValue={(V(V+)-V(V-))*I(Vsense)} be calculated as Cth = Cp tox A [18], where ; G1 value A is the CF effective area, the effective CF length that approximately corresponds to the dielectric thickness (tox) and Figure Figure 10.10. Equivalent electric Equivalent circuit electric of a RRAM circuit thermalthermal of a RRAM model based modelonbased a thermal on aresistance thermaland resistance and C p is the volumetric heat capacity (calculated as ρ × c, Equation (2)) [89]. As detailed in capacitance to to implement Equation (20). (20). capacitance implement Equation [18], the thermal resistance can be calculated as follows: The values κ = 2.8–25 µJ/K s (Rth = 4 × 104–3.5 × 105 K/W) and Cth = 0.04–1.1 pJ/K were 𝑡 employed in [78]. In [95], the following ≈R were 𝑅 values tox given: Cth = 0.318 fJ/K and τth = 0.23 (28) ns, from them the thermal resistance can be𝑘th ≈𝐴 k ARth = τth/Cth = 7.23 × 105 K/W. The extracted (28) th thermal resistance in [91] was taken Rth = 5 × 105 K/W. An estimation of 33 ps for the ther- Assuming Assuming mal time isasreported constant as thereference the reference material material in [18], which Hf, Hf,towe has we the behave have compared tothe followingfollowing values the electric kthvalues kth−=1 K23 = 23 Wm pulse-width −1 , and −1CKp −1=, 1.92 Wm and CJKp−=11.92 cm−JK3 . −1cm−3. Therefore, Therefore, the thermal the thermal time constant time constant can be can be estimated estimated as: as: 𝐶 𝑡 2 =Cth = C p tox 𝑅 𝐶 Rth (29) (29) 𝑘 k th which which is is 33 33 ps ps for for the the case case considered considered (t (tox = 20 nm). ox = 20 nm). A A quick estimation to assess the role ofofthe quick estimation to assess the role thethermal thermalresistance resistance (R(R th = 4 × 1044 K/W) can th = 4 × 10 K/W) can be be performed performed ifif we we assume assume an an ideal ideal device device in in the the LRS under aa steady-state LRS under steady-state regimeregime withwith IIRRAM RRAM = 1 mA and VRRAM = 1 V; using Equation (19), we would have T = T0 +×RIth× × = 1 mA and V RRAM = 1 V; using Equation (19), we would have T = T 0 + Rth V I= × T0V+ 40 = TK, or40T K, 0 +T500 = Tor = TK, 5 0 +if500 Rth =K,5 if× R 10th K/W. 5 0+ = 5 × 10 K/W. The The role role of the heat capacitance can be easily seen in Figure 11. For a pulsed voltage signal signal ofof amplitude amplitude(peak (peaktotopeak)= peak)=0.025 0.025V,V,offset offset= =0.0125 VV 0.0125 and and f =f 0.5 GHz = 0.5 GHz applied applied in an ideal in an idealdevice devicewhosewhose resistance resistance is assumed is assumed to to bebeRRRAM RRRAM Ω, the = 1 =Ω,1 the temperature temperature cancan be obtained be obtained fromfromEquation Equation (27) (27) for Rfor th =R2th×= 1025 × K/W105and K/W Cthand = 0.1CfJ/K th = (τ = 20 ps), 0.1th fJ/K (τ th C=th20 = 0.25 ps), fJ/K Cth =(τ0.25 th = 50 ps),(τCthth== 50 fJ/K 0.5ps), fJ/KC(τ th th==0.5100fJ/K ps).(τ th = 100 ps). Figure 11. (a) Voltage applied to the the device device versus time, (b) Temperature Temperature versus versus time time obtained obtained for for different different values values of of C Cthth (assuming (assuming RRthth ==22××1010 5 K/W). 5 K/W). As can be seen in the previous figure, the device thermal response should be faster than the electric pulses employed for an operation free of delays (“inertia”) due to thermal effects. In fact, experimental techniques have been developed to estimate the device tem- perature with the application of ultra-short electrical pulses to RRAMs [97,98]. Nanomaterials 2021, 11, 1261 16 of 47 Table 2. Model parameters employed for the simulations performed with the Stanford model, in Nanomaterials 2021, 11, 1261 particular for Figures 12 and 13. 15 of 46 Stanford-PKU Model Parameters Device Parameters Unit Resistive Switching As can be seen in the previous figure, the device thermal response should be faster SET RESET than the electric pulses employed for an operation free of delays (“inertia”) due to thermal Vo V 0.4 effects. In fact, experimental techniques have been developed to estimate the device I0 mA 0.2 temperature with the application of ultra-short electrical pulses to RRAMs [97,98]. g0 nm 0.35 We have made use of some of the thermal models reported above. They cannot exist ν0 m/s 106 on their own since the conductive filament kinetics need to be considered to describe the α - 1 device RS operation and, in doing so, calculate the current. The use of the thermal models β - 3 TM1-TM4 has been presented previously in [27,56,79]. We show here some results of the γ0 - 10 Stanford model [78,88,93,95] along with the thermal models TM5 and TM6. For the device description, taking into consideration that no experimental data fitting is considered (since Several I-V curves were plotted in Figure 12 considering different thermal re- it was performed in the references where the models were introduced), we assumed a set sistances. As can be seen, the main dissimilarities are found in the set and reset regions. of model parameters close to the one suggested in [88], see Table 2. For the Stanford model, the highest temperatures are found in the set process, see that the Joule heating Table 2. Model is higher inemployed parameters the set process and in the performed for the simulations interval ofwith positive voltages the Stanford above model, in V SET. This is coherent with experimental findings that show that the Joule heating role is particular for Figures 12 and 13. essential to describe the SET kinetics [81,99]. On average, these maximum temperatures achieved are in line with the estimations Stanford-PKU performed above for the thermal resistances con- Model Parameters sidered. It is important Device Parameters to highlight Unit the importance of the thermal Resistiveresistance Switching in the device design. See that higher Rth values lead to lower set and reset voltages to produce RS, and SET RESET hence, lower power consumption. From this viewpoint, higher thermal resistances (i.e., V devices showingo V 0.4 a character thermodynamically more adiabatic with respect to the sur- I roundings) 0might be more interesting.mA Since the heat flux from the0.2 CF to the metallic elec- trodes (operating g0 these as heat sinks, nm a reasonable assumption, allows 0.35 to calculate the heat flux with the material 3D thermal conductivity) could be in the same order of magnitude ν0 m/s 106 for different RRAMs, the dielectric could be the key (save the role of the contact thermal - resistance).αIn this respect, we have to call the reader’s attention to1the fact that the usual RRAM dielectrics β - are grown in nanometric layers; consequently, the 3 real thermal conduc- tivity due toγ0 phonon quantization- is far from the values corresponding 10 to the correspond- ing 3D dielectric materials, as it was shown in [100–102]. Figure 12. Simulations Simulations performed performed making making use use of of the the Stanford Stanford model model including including the theTM5 TM5with withdifferent different thermal thermalresistances. resistances. (a) Current versus voltage applied to the device, (b) voltage signal applied to the device, (c) temperature versus time. (a) Current versus voltage applied to the device, (b) voltage signal applied to the device, (c) temperature versus time. Several I-V curves were plotted in Figure 12 considering different thermal resistances. In the case of TM6 along with the Stanford model (using Table 2), we have plotted As can be seen, the main dissimilarities are found in the set and reset regions. For the the results of a simulation using a pulsed voltage signal in Figure 13. A voltage signal such Stanford model, the highest temperatures are found in the set process, see that the Joule as the one in Figure 12b is close to DC, therefore the use of TM6 instead of TM5 is irrele- heating is higher in the set process and in the interval of positive voltages above VSET . This vant. However, if a fast voltage signal is employed (Figure 13a), a high enough thermal is coherent with experimental findings that show that the Joule heating role is essential to capacitance makes a difference in terms of the device temperature transient (see Figure describe the SET kinetics [81,99]. On average, these maximum temperatures achieved are 13b). in line with the estimations performed above for the thermal resistances considered. It is important to highlight the importance of the thermal resistance in the device design. See that higher Rth values lead to lower set and reset voltages to produce RS, and hence, lower power consumption. From this viewpoint, higher thermal resistances (i.e., devices showing a character thermodynamically more adiabatic with respect to the surroundings) might be more interesting. Since the heat flux from the CF to the metallic electrodes (operating these as heat sinks, a reasonable assumption, allows to calculate the heat flux with the Nanomaterials 2021, 11, 1261 16 of 46 material 3D thermal conductivity) could be in the same order of magnitude for different RRAMs, the dielectric could be the key (save the role of the contact thermal resistance). In this respect, we have to call the reader’s attention to the fact that the usual RRAM dielectrics are grown in nanometric layers; consequently, the real thermal conductivity Nanomaterials 2021, 11, 1261 17 of 47 due to phonon quantization is far from the values corresponding to the corresponding 3D dielectric materials, as it was shown in [100–102]. Figure Figure 13.13. Simulations Simulations performed performed with with the the Stanford Stanford model model (including (including TM6) TM6)making makinguse useof ofdifferent differentthermal thermalcapacitances, capacitances, C , assuming a common value of the thermal resistance, R = 4 × 105 5K/W. (a) Voltage applied to the device versus time, Cth , assuming a common value of the thermal resistance, Rth = 4 × 10 K/W. (a) Voltage applied to the device versus time, th th (b) (b) Temperature versus time. Temperature versus time. The In the thermal case oftime TM6constants along with (Equation (22)) model the Stanford produced (usingobvious Tabledelays 2), we for have a voltage plotted signal such of the results as athe one in Figure simulation using13. a These pulseddelaysvoltage could seriously signal in Figure affect 13. the device opera- A voltage signal tion suchunder as thepulsed one in input Figuresignals 12b is since close RS is most to DC, of the the therefore times uselinked of TM6 to instead thermalofeffects TM5 [1,24,53,81,99]. The thermal is irrelevant. However, if atime fast constants corresponding voltage signal is employed to Figure (Figure13b areathe 13a), following: high enough 0.4, 0.2 and thermal 0.1 ns. These capacitance makes values are largeincompared a difference terms of the to the value device (33 ps) given temperature by Ielmini transient (see [18] (we13b). Figure do not go into details of device materials at this point). An estimation of the Cth The thermal associated to a CF time constants (Equation in a conventional RRAM, (22)) produced as described obvious in [18], could delays help for a voltage to shed light signal on this such issue.as Letthe usone assume in Figure 13. These a cylindrical CF delays radius of could 5 nm seriously affect layer in a dielectric the device of 20 nm op- eration thick, if under the heat pulsed capacity input of Hfsignals since RS (C is considered is pmost = 1.92ofJKthe −1cm times −3) thelinked to thermal CF thermal capaci-ef- fects [1,24,53,81,99]. tance would be (Cth =The Cp tthermal ox A) 0.003 time constants fJ/K. Therefore, corresponding a value τth = 1.2 to Figure 13b are the ps is expected fol- for the lowing: same 0.4, 0.2resistance thermal and 0.1 ns. These values employed are large in Figure comparedthis 13. Although to the value time thermal (33 ps) given by constant is Ielmini [18] (we do not go into details of device materials short, different authors [78,95] have used thermal capacitances values that lead to devices at this point). An estimation of the with Cth associated higher to a CF in thermal constants. a conventional A thermal device model RRAM, as described described in [18],(20) by Equation couldandhelpthe to shed light on this issue. Let us assume a cylindrical CF thermal capacitance of an average CF produces so low thermal time constants that no tran- radius of 5 nm in a dielectric layerterm of 20innm thick, if(20)thewould heat capacity −1 −3 ) the CF sient Equation be worthofbeing Hf istaken considered (Cp = 1.92 into account. FromJKthe cm experimental thermal capacitance viewpoint, no delayswould linkedbe to(C th = Cp tinertia thermal ox A) 0.003 wouldfJ/K. be Therefore, a value τ th =memory seen for conventional 1.2 ps is expected for the same thermal resistance employed in Figure pulsed signals. Nevertheless, current transients on longer time scales than the previously 13. Although this thermal time constant calculated is short, τth, linked different to some extentauthors [78,95] to thermal havehave effects usedbeen thermal reportedcapacitances previously values [78]. that lead to devices with higher thermal constants. A thermal In this respect, the thermal model, i.e., Equation (20), might not be enough to accurately device model described by Equation describe RRAM (20)thermal and theresponse. thermal capacitance of an average CF produces so low thermal timeTheconstants that no transient coupling of the CF temperature term in to Equation a heat sink (20)atwould room be worth being temperature with taken into just two account. From parameters the experimental (thermal resistance and viewpoint, no delays capacitance) couldlinked not betodescribed thermal inertiaby thewouldthermal be seen for conventional memory pulsed signals. Nevertheless, features of a nanometric CF. We could use an intermediate temperature corresponding to current transients on longer time an scales than average region thesurrounding previously calculated the CF that , linked τ thcould to some help to build extent twotodifferent thermalthermal effects have cir- been reported previously [78]. In this respect, the thermal cuits, one accounting for the coupling between the CF and this intermediate surrounding model, i.e., Equation (20), might not be enough to accurately describe RRAM thermal response. region (it could include the dielectric and electrode zones closer to the CF), and a second The coupling of the CF temperature to a heat sink at room temperature with just two one coupling between this intermediate region and the outside world. This is in line with parameters (thermal resistance and capacitance) could not be described by the thermal previous thermal descriptions (although introduced in another context) employed for features of a nanometric CF. We could use an intermediate temperature corresponding magnetic current sensors [103] and AlGaAs/GaAs HBTs [104]. For this purpose, we pre- to an average region surrounding the CF that could help to build two different thermal sent the following model. We will show below that the model based on Equation (20) can be used if a thermal capacitance that accounts for the whole device is taken into consider- ation. This means using a parameter much higher than the exclusively associated to the CF. In doing so, as always in compact modeling, the simplicity and accuracy trade-off has to be considered. Nanomaterials 2021, 11, 1261 17 of 46 circuits, one accounting for the coupling between the CF and this intermediate surrounding region (it could include the dielectric and electrode zones closer to the CF), and a second one coupling between this intermediate region and the outside world. This is in line with previous thermal descriptions (although introduced in another context) employed for magnetic current sensors [103] and AlGaAs/GaAs HBTs [104]. For this purpose, we present the following model. We will show below that the model based on Equation (20) can be Nanomaterials 2021, 11, 1261 used if a thermal capacitance that accounts for the whole device is taken into consideration. 18 of 47 This means using a parameter much higher than the exclusively associated to the CF. In doing so, as always in compact modeling, the simplicity and accuracy trade-off has to be considered. 2.4.3. Non-Steady-State Approach with Two Different Temperatures Associated to the Device (Second-Order Memristor) 2.4.3. Non-Steady-State Approach with Two Different Temperatures Associated to the This(Second-Order Device thermal model, as suggested above, employs two different temperatures to de- Memristor) scribe the device from an energy balance perspective: the internal device temperature (ap- This thermal model, as suggested above, employs two different temperatures to proximately the CF temperature, T) that affects the RS mechanisms linked to the CF crea- describe the device from an energy balance perspective: the internal device temperature tion and destruction and a second temperature, associated to the CF surrounding regions (approximately the CF temperature, T) that affects the RS mechanisms linked to the CF (an effective creation andtemperature, destruction Tand S). The latter influences the internal device temperature T but a second temperature, associated to the CF surrounding itregions shows (ana different effective temperature, TSThe time evolution. device ). The latterintermediate influences the surrounding region internal device (at tem- temperature perature T S) is characterized by an outer boundary assumed to be at room temperature T but it shows a different time evolution. The device intermediate surrounding region (at (T 0 = 300 K). This outer boundary is considered to be far away from the RS active region. temperature TS ) is characterized by an outer boundary assumed to be at room temperature The intermediate (T0 = surrounding 300 K). This outer boundary region can include is considered to bedifferent far awaymaterial from thelayers; RS activetherefore, ef- region. The fective thermal constants are employed to account for the heat flux intermediate surrounding region can include different material layers; therefore, effective between this region and the exterior thermal constants zone. areBesides, employed theto coupling accountbetween the inner for the heat (CF volume) flux between and theand this region inter- the mediate CF surrounding region could be modeled by an effective exterior zone. Besides, the coupling between the inner (CF volume) and the intermediate thermal resistance and thermal capacitance: CF surrounding Rth1 could region and Cth1 be. Under modeled thisbyapproach, an effectivethe thermal device can be described resistance by the and thermal following capacitance:two R equations th1 and C(we th1 . assume Under this this procedure approach, to thebe the device thermal can bemodel 7, described TM7;by see the the circuital implementation in Figure 14). following two equations (we assume this procedure to be the thermal model 7, TM7; see the circuital implementation in Figure 14). 𝑑(𝑇 − 𝑇 ) 1 𝐶 = 𝑉(𝑡)𝐼(𝑡) − (𝑇 − 𝑇 ) (30) 𝑑𝑡 d( T − TS ) 𝑅 1 Cth1 = V (t) I (t) − ( T − Ts ) (30) dt Rth1 𝑑𝑇 1 𝐶 = dT 𝑉(𝑡)𝐼(𝑡) s − (𝑇 1− 𝑇 ) (31) 𝑑𝑡 th2 C = V ( ) 𝑅( ) − t I t ( Ts − T0 ) (31) dt Rth2 . where whereV(t)I(t) V(t)I(t)allows allowsthe thedetermination of e𝑒generated . .The determination of Thepower powerdissipated dissipatedbybyJoule Jouleheating heat- ing affects the intermediate surrounding region temperature (modeled with affects the intermediate surrounding region temperature (modeled with parameters Rth2 parameters Rand th2 and C)th2that Cth2 ) that accountsfor accounts forthe thecoupling couplingbetween betweenthis this region region and the thermalized thermalizeddevice device exterior exteriorregion. region.The Theapproach approachdescribed describedhere hereisisininline linewith withthe thedescription descriptionof ofaasecond second order ordermemristor memristor[100].[100]. .param T0= 300 ; K T Rth1 Cth1 V+ Ts V- G1 Rth2 Cth2 I+ Vsense 0 VT1 {T0} I- value={(V(V+)-V(V-))*I(Vsense)} ; G1 value Figure14. Figure Equivalentelectric 14.Equivalent electriccircuit circuitofofa aRRAM RRAM thermal thermal model model based based onon a double a double thermal thermal circuit circuit described describedby byEquations Equations(30) (30)and and(31). (31). The simulation was performed with this thermal model integrated in the Stanford compact model that was employed previously, see Figure 15. Nanomaterials 2021, 11, 1261 18 of 46 Nanomaterials 2021, 11, 1261 19 of 47 Nanomaterials 2021, 11, 1261 19 of 47 The simulation was performed with this thermal model integrated in the Stanford compact model that was employed previously, see Figure 15. Figure 15. (a) Three-dimensional view of the Stanford scheme to model the different device areas (TE: top electrode, Oxide Layer, CF: conductive filament and BE: bottom electrode), (b) model parameters, g: gap between the TE and the filament Figure 15. Figure (a)Three-dimensional 15. (a) Three-dimensional viewview of of the theStanford Stanford scheme scheme toto model model the thedifferent different device device areas areas (TE: (TE: top top electrode, electrode, Oxide Oxide tip and tox: dielectric thickness, (c) subcircuit representation for the implemented model. The connection between blocks Layer, Layer, CF: conductive filament and BE: bottom electrode), (b) model parameters, g: gap gap between between the the TE TE and and the the filament filament represents the states variables used: g, which depends on kinetic block and it is linked to the two temperatures (T, TS). tip tip and and tox ox:: dielectric dielectricthickness, thickness, (c) (c) subcircuit subcircuit representation representation for for the the implemented implemented model. model. The The connection connection between between blocks blocks represents the states variables used: represents the states variables used:The g, which g, which depends depends on kinetic on kinetic block and it is linked to the two temperatures (T, T ). (T, TSSparameters results obtained for theblock andparameters set of it is linked toofthe thetwo temperatures standard Stanford ). [88]The and the thermal model shown setofin Figure 14; Figure 15standard are givenStanford in Figure 16. The results obtained obtained forthe for theset ofparameters parameters of of thethe standard parameters Stanford parameters [88] [88] and and the thermal the thermal model model shownshown in Figure in Figure 14; Figure 14; Figure 15 are15given are given in Figure in Figure 16. 16. 2 Cth2=10fJ/K Cth2=10fJ/K Cth2=1fJ/K 2 Cth2=1fJ/K 2 Cth2=10fJ/K Cth2=10fJ/K Current (mA) Cth2=1fJ/K 2 Cth2=1fJ/K Voltage (V) 1 1 Current (mA) Voltage (V) 1 1 0 0 0 0 -1 -1 -1 a) -1 -2 b) a) 2 0 4 6 8 10 12 14 16 18 20 22 b) 20 4 6 8 10 12 14 16 18 20 22 Time (ns) -2 0 2 4 6 8 Time (ns) 10 12 14 16 18 20 22 0 2 4 6 8 10 12 14 16 18 20 22 Time (ns) Time (ns) Figure 16. RRAM simulation making use of the Stanford model including a double RC thermal model Rth1 = 40 kK/W, Rth2 = 40 kK/W, Cth1 = 0.003 fJ/K and Cth2 with values 1 fJ/K and 10 fJ/K. (a) Applied voltage pulses for consecutive set Figure and reset,16. (b)RRAM temporalsimulation current making evolution,use(c)oftemporal the Stanford modelof evolution including a doubletemperature device filament RC thermal(T) model and Rth1 = the 40 kK/W, Rth2 intermediate = 40 kK/W, region surrounding Cth1 = 0.003 (TS ) fJ/K with and Cth2 C=th21with fJ/K values and 1 (d) fJ/K C and = 10 10 fJ/K. fJ/K. (a) Applied voltage pulses for consecutive set and reset, th2 model including a double RC thermal model Rth1 = 40 kK/W, Rth2 Figure 16. RRAM simulation making use of the Stanford (b) temporal current evolution, (c) temporal evolution of device filament temperature (T) and the intermediate surround- = 40 kK/W, Cth1 = 0.003 fJ/K and Cth2 with values 1 fJ/K and 10 fJ/K. (a) Applied voltage pulses for consecutive set and reset, ing region (TS) with Cth2 = 1 fJ/KDepending and (d) Cth2 =on 10CfJ/K. th2 ,of different CF and intermediate surrounding region temperatures (b) temporal current evolution, (c) temporal evolution device filament temperature (T) and the intermediate surround- transient responses ing region (TS) with Cth2 = 1 fJ/K and (d) Cth2 = 10 fJ/K. are obtained, producing the corresponding effects on the device current. This isDepending on Cth2a, consecutive noticeable when different CFseries and intermediate surrounding of set and reset pulses are region temperatures applied, as shown transient responses are obtained, producing the corresponding effects on the Depending on Cth2, different CF and intermediate surrounding region temperatures device cur- rent. This is noticeable when a consecutive series of set and reset pulses transient responses are obtained, producing the corresponding effects on the deviceare applied, cur-as rent. This is noticeable when a consecutive series of set and reset pulses are applied, as Nanomaterials 2021, 11, 1261 19 of 46 in Figure 16a, in which a sequence of set pulses (1.5 V with 1 ns on time and 0.1 ns rise and fall times) and reset pulses (−1.5 V with 1 ns on time and 0.1 ns rise and fall times). In the first configuration, with the thermal capacities Cth1 = 0.003 fJ/K, Cth2 = 1 fJ/K and thermal resistances Rth1 = Rth2 = 40 kK/W, the current evolution is shown in Figure 16b. Figure 16c shows the CF (T) and intermediate surrounding region (TS ) evolution. In this first configuration, the corresponding transient shows low thermal inertia; after the pulse application, both temperatures reach room temperature (T = TS = T0 ). The devices show a slight increase in the maximum temperatures obtained in the set (TSET ) and reset (TRESET ) processes (Figure 16c). In the second configuration, thermal capacities Cth1 = 0.003 fJ/K, Cth2 = 10 fJ/K and thermal resistances Rth1 = Rth2 = 40 kK/W, the current evolution is plotted in Figure 16b. Figure 16d shows the CF (T) and intermediate surrounding region (TS ) evolution. In this second configuration, the thermal inertia is higher in the second thermal circuit, after each set/reset pulse, the temperatures cannot go back to room temperature (T, TS 6= T0 ). As a result, each new cycle starts from a higher temperature than in the previous cycle; therefore, the maximum set (TSET ) and reset (TRESET ) temperatures show a growing trend (see Figure 16d). This temperature increase over the cycles implies that the device CF gap decreases in each new cycle, then the current increases (Figure 16b). This effect suggests the consideration of the temperature, in addition to the CF gap, as a state variable in line with the approach presented in [100] for second-order memristor. The temperature increase reported above could be employed with a series of pulses to tune the device con- ductivity in set cycles within a neuromorphic circuit context [100,105,106]. It is noteworthy that a third-order approach has been introduced in modeling devices for neuromorphic engineering [107]. 2.4.4. SPICE-Based Circuital Models with Two or More CF Temperatures In Section 2.3.4, a compact model based on two different CF representative tempera- tures was reviewed. In that model, some simplifications were performed in order to obtain analytical expressions for calculating these two temperatures although the flexibility linked to the consideration of different temperatures for the CF narrowing and CF main body added high value to the modeling process. In this section, an approach based on the electri- cal equivalent representation of this thermal framework is proposed in order to calculate numerically the temperature. It is noteworthy that most of the thermal models based on an electrical equivalent circuit (Sections 2.4.2 and 2.4.3) assume constant thermal resistances and capacitances, which can be used as fitting parameters. However, heat conduction through the CF and lateral heat dissipation depends on the filament size. This effect is in- cluded in physically-based simulators, such as those based on a kMC approach or on finite differences/elements, as far as the heat equation is solved using thermal parameters which are updated at simulation time according to the current filament size. Thermal models based on the heat equation analytical solutions in simplified CF geometries (Section 2.3) incorporate also this effect because the temperature analytical expression depends on the actual geometry, and it is evaluated at each simulation time step. In this section, we present two thermal models based on an equivalent electrical representation (SPICE based), but including the effects of the CF evolution on the thermal properties, which also evolve as the simulation proceeds. The first one is a steady-state model, while the latter includes thermal capacitances. Furthermore, both models account for longitudinal heat conduction and lateral heat dissipation by means of several thermal resistances (or RC networks in the non-steady approach). Figure 17 shows a general schema of the overall model. As explained before, the temperature of the CF main body is, in general, lower than that of the small portion of the filament where the filament evolves faster. Therefore, the CF is modeled by two different subcircuits (Figure 17), each of them characterized by different state variables (CF radius and temperature). Nanomaterials 2021, 11, 1261 21 of 47 Nanomaterials 2021, 11, 1261 20 of 46 Figure 17. Figure 17. Schema Schemaofofthethecircuital circuital compact compact model. model. A truncated-cone A truncated-cone shaped shaped conductive conductive filament filament is represented is represented by con- by connected nected cylinders cylinders for modeling for modeling purposespurposes (a). The behaviour (a). The behaviour of eachofportion of each portion of the(cylinder) the filament filament (cylinder) is modeledisby modeled by the the subcircuit subcircuit inside the blue rectangle (b), which has electrical connections (EC) and thermal connections (TC). Each cylinder inside the blue rectangle (b), which has electrical connections (EC) and thermal connections (TC). Each cylinder (subcircuit) (subcircuit) is characterized by different state variables (radius, r_cf, and temperature, Temp). The cylinder subcircuit is characterized by different state variables (radius, r_cf, and temperature, Temp). The cylinder subcircuit consists of several consists of several more subcircuits: a kinetic block for calculating the transient CF evolution; an electrical block for current more subcircuits: calculation; a kineticthe and, finally, block for calculating thermal subcircuit,the transient which CF evolution; includes an electrical the equivalent block circuit for thefor currentmodel. thermal calculation; As canand, be finally, the thermal subcircuit, which includes the equivalent circuit for the thermal model. As can be seen, the subcircuits (thermal, kinetic and electrical blocks) are connected all together because they are interdependent. Ifseen, the subcircuits (thermal, necessary,kinetic and electrical a last subcircuit blocks) is added are connected in series (a) in orderalltotogether account because they are interdependent. for the conduction If necessary, through a constriction a last by means of subcircuit the quantum is added pointin series model contact (a) in order to account (see Section for the conduction through a constriction by means of the quantum point 4) [108]. contact model (see Section 4) [108]. Now, we focus on the description of the thermal subcircuit (Figure 18). The longitu- dinalNow, we focus onisthe heat conduction description modeled by Rof Thl1the andthermal subcircuit RThl2. For the sake(Figure 18). Theitlongitudi- of generality, has been nal splitheat conduction off into is modeledinbyorder two contributions RThl1toand RThl2 make . Forthe easier theconnection sake of generality, with otherit has been thermal split off intoand sub-circuits twotocontributions in order to build more complex make models. thermal easier theTheir connection with values are otherbythermal given [89]: sub-circuits and to build more complex thermal models. Their values are given by [89]: 1 𝐿 𝑅 = 𝑅 = (32) R Thl1i = R Thl2i 2 𝑘= 1𝜋 𝑟 Li (32) 2 k th π ri 2 where kth is the CF thermal conductivity, Li is the length of the portion of the filament where modeled kth by is the the CF thermaland sub-circuit conductivity, ri, its radius Li (index is the length i refersoftothe theportion cylinderoforthe filament sub-circuit modeled by the sub-circuit and number 1 or 2 in Figure 17a). On the r i , its radius (index i refers to the cylinder or sub-circuit other hand, RThn accounts for the lateral heat dissipa- number tion and 1 orcalculated it is 2 in Figure 17a). On following thisthe other hand, expression [89]: RThn accounts for the lateral heat dissipation and it is calculated following this expression [89]: 1 𝑅 = (33) 2 ℎ=𝐿 𝜋 𝑟 1 R Thni (33) 2 h Li π ri Note that in this model (TM8), the thermal resistances are not directly the fitting pa- Note that since rameters in thisthey model are(TM8), the thermal calculated according resistances are notfilament to the actual directly size. the fitting parameters Therefore, as far since they are calculated according to the actual filament as the complete device model is able to reproduce the geometrical evolution of size. Therefore, as farthe as fila- the complete device model is able to reproduce the geometrical evolution ment (kinetic block in Figure 17), the thermal resistances are not fixed, but their values of the filament (kinetic evolve asblock the in Figure 17), simulation the The runs. thermal resistances are implementation ofnot suchfixed, but theirresistances dependent values evolve as can be the simulation runs. The implementation of such dependent resistances can be made by made by means of behavioural sources. means of behavioural sources. Although the model shown in Figure 17 uses two cylinders (and the corresponding Although the model shown in Figure 17 uses two cylinders (and the corresponding two two sub-circuits), the CF could be represented by more cylinders in order to get a more sub-circuits), the CF could be represented by more cylinders in order to get a more detailed detailed CF description [31]. Indeed, in the limit with an infinite number of differential CF description [31]. Indeed, in the limit with an infinite number of differential length length cylinders, the circuit simulation would be equivalent to the resolution of the differ- cylinders, the circuit simulation would be equivalent to the resolution of the differential ential Equation (9). In fact, with a reduced number of filaments, the results are very similar Equation (9). In fact, with a reduced number of filaments, the results are very similar to to those obtained with a finite differences simulator [31,109]. Furthermore, complex fila- those obtained with a finite differences simulator [31,109]. Furthermore, complex filaments ments such as those with several branches forming a tree structure or interlaced between such as those with several branches forming a tree structure or interlaced between them ials 2021, 11, 1261 22 of 47 Nanomaterials 2021, 11, 1261 21 of 46 them could be simulated using more blocks and interconnecting them following a given structure (Figures 5 and could be6simulated in [31]). using more blocks and interconnecting them following a given structure On the contrary, if only (Figures 5 andone6 block is used, this thermal model (TM8) would be equiv- in [31]). alent to TM5, although Onthe thethermal resistance contrary, evolution if only one block isisused, allowed this in TM8. model (TM8) would be equiva- thermal lent to TM5, although the thermal resistance evolution is allowed in TM8. Figure 18. Circuital equivalent for thermal model TM8. The circuit inputs are the dissipated power Figure 18. Circuital equivalent for thermal model TM8. The circuit inputs are the dissipated power (pw) and the CF radius (r_cf), while the output is the temperature (T_CF). The actual values of the (pw) and the CF radius (r_cf), while the output is the temperature (T_CF). The actual values of the thermal resistances depend on the filament radius (it is assumed to be a cylinder) and, therefore, thermal resistances depend on the filament radius (it is assumed to be a cylinder) and, therefore, they are updated asthey theare CF updated as the evolves. The CF evolves. subcircuit has The beensubcircuit prepared has been prepared for being connectedfortobeing other connected to other thermal subcircuits (through TC1, TC2 and TCox) in order to obtain a more complex thermal complex thermal model thermal subcircuits (through TC1, TC2 and TCox) in order to obtain a more model of the wholeofdevice the whole (withdevice several(with several temperatures temperatures along the along the filament or filament different or different temperatures for the temperatures for the surrounding insulator orinsulator surrounding bulk insulator, or bulkFigure 17). If insulator, only 17). Figure one block If onlyisone used, all the block re- all the resistances are is used, sistances are in parallel and the model is equivalent to TM5, although the thermal resistance in parallel and the model is equivalent to TM5, although the thermal resistance value value keeps the keeps the dependency on the filament size in TM8 and it changes during the simulation. dependency on the filament size in TM8 and it changes during the simulation. Figure 19 shows the results Figure provided 19 shows by TM8provided the results when is coupled to kinetic by TM8 when and conduc- is coupled to kinetic and conduc- tive blocks fitted tive blocks fitted to simulate reset transitions in unipolar Ni/HfO2 /Si-n+ resistive switching to simulate reset transitions in unipolar Ni/HfO 2 /Si-n + resistive switching devices [80,109,110]. The[80,109,110]. devices QPC block has Thealso QPC been added. block Thebeen has also two added. cylinders-TM8 The twomodel cylinders-TM8 model results have beenresults comparedhavewith beenthose provided compared by those with a finiteprovided differences bysimulator that was simulator that a finite differences was used todata used to fit the experimental fit the experimental [80]. data dissipation The lateral heat [80]. The lateral h, has beenparameter, h, has heat dissipation parameter, changed in orderbeen changed to check in order on its influence to check the i-vits influence curve. As canonbe the i-v curve. seen, As can more heat be seen, more heat dissi- pation requires adissipation requires higher voltage a higher to reach voltage totriggered the thermally reach thereset thermally triggered transition, reset transition, a a well- well-known known effect in RRAMs. effect in RRAMs. Thermal inertia can also be considered following this approach if thermal capacitances are added (Figure 20, thermal model TM9) [37]. In this new model, thermal capacitances could also be made dependent on the filament size instead of being fixed fitting parame- ters [37,111]. As previously mentioned, several blocks can be connected for modeling the CF. On the contrary, if only one segment is used, the model is equivalent to model TM6 (although TM9 lets the thermal components evolve at simulation time). Figure 21 shows the simulation of the transient response of a Ni/20 nm-HfO2 /Si-n+ resistive switching device [110] when a 3 V reset pulse is applied (for 100 ns) [37]. Several values of the thermal capacities have been considered. The simulation context here is different from the one shown in Section 2.4.3 since all the modeling components are linked to the device conductive filaments. Note also that only values higher than 0.2 fJ/K influence the device response. Fixed and variable thermal capacities have been used in order to analyze the role of size-dependent thermal capacities, which evolve at simulation time. As expected, variable thermal capacitors, whose value is reduced during a reset process, produce lower thermal inertia. Nanomaterials2021, Nanomaterials 2021,11, 11,1261 1261 23 of 46 22 of 47 Figure 19. Simulation of a reset transition in unipolar Ni/20 nm-HfO2/Si-n+ resistive switching de- vices [109]. The i-v curve provided by a finite differences simulator used to fit the experimental data [80] is compared with the i-v curve calculated by means of the two cylinders model with TM8. Note that with only two subcircuits (Figure 17a) and taking variable electrical and thermal resistances into account, both types of simulators provide very close results, as far as the circuital model includes variable electric and thermal resistances. Two values of the lateral heat dissipation parameter, h, have been used for the sake of comparison. Figure 19.be of of aa reset transition transitionininunipolar Ni/20 nm-HfO 2/Si-n + resistive + resistive switching de- Thermal inertia can also Figure 19. considered Simulation following reset this approach unipolar Ni/20 if thermal nm-HfO 2 /Si-n capaci- switching vices [109]. devices [109].TheThei-vi-v curve curveprovided provided byby a finite a differences finite simulator differences simulatorused usedto fit to the fit experimental the experimental tances are added (Figure 20,[80] data[80] thermal model TM9) [37].calculated In this bynew model, thermal capac- data isiscompared comparedwith withthe thei-v i-vcurve curvecalculated meansof by means ofthe thetwo twocylinders cylindersmodel modelwith withTM8. itances could also be made dependent TM8. Note on the that with only filament(Figure two subcircuits size instead of being 17a) and taking variablefixed fitting electrical and thermal Note that with only two subcircuits (Figure 17a) and taking variable electrical and thermal resistances resistances parameters [37,111]. Asinto previously into account, both types of simulators provide very close results,foras far as the circuital account, bothmentioned, types of simulators several provideblocks canresults, very close be connected as far as the circuital mod- model includes model includes variable electric and thermal resistances. Two values of the lateral heat dissipation eling the CF. On the contrary, variable if electric only and one thermal segment resistances. is Two used, parameter, h, have been used for the sake of comparison.values the of the model lateral is heat equivalent dissipation to h, have parameter, model TM6 (although TM9 lets the thermal components evolve at simulation time). been used for the sake of comparison. Thermal inertia can also be considered following this approach if thermal capaci- tances are added (Figure 20, thermal model TM9) [37]. In this new model, thermal capac- itances could also be made dependent on the filament size instead of being fixed fitting parameters [37,111]. As previously mentioned, several blocks can be connected for mod- eling the CF. On the contrary, if only one segment is used, the model is equivalent to model TM6 (although TM9 lets the thermal components evolve at simulation time). Figure 20. Circuital equivalent for thermal model TM9. It is similar to TM8 (Figure 18), but thermal Figure 20. Circuital equivalent for thermal model TM9. It is similar to TM8 (Figure 18), but ther- inertia has been added by means of capacitances. As in TM8, the actual values of the thermal mal inertia has been added by means resistances andofcapacitances capacitances. As on depend in TM8, the actual the filament radiusvalues of the thermal (it is assumed to be a cylinder) and, resistances and capacitances depend on the filament radius (it is assumed therefore, they are updated as the CF evolves [37]. to be a cylinder) and, Figure 20. Circuital equivalent for thermal model TM9. It is similar to TM8 (Figure 18), but ther- therefore, they are updated as the CF evolves [37]. mal inertia has been added by means of capacitances. As in TM8, the actual values of the thermal resistances and capacitances depend on the filament radius (it is assumed to be a cylinder) and, therefore, they are updated as the CF evolves [37]. to the device conductive filaments. Note also that only values higher than 0.2 fJ/K influ- ence the device response. Fixed and variable thermal capacities have been used in order to analyze the role of size-dependent thermal capacities, which evolve at simulation time. As expected, variable thermal capacitors, whose value is reduced during a reset process, Nanomaterials 2021, 11, 1261 produce lower thermal inertia. 23 of 46 Figure 21. Simulated current in a Ni/20 nm-HfO2 /Si-n+ resistive switching device [37,110] when a 3 V reset pulse is applied (for 100 ns). Different values of the thermal capacities have been assumed. Figure 21. Simulated current in a Ni/20 nm-HfO2/Si-n+ resistive switching device [37,110] when a 3 Fixedpulse V reset thermal capacitances is applied (solid (for 100 lines) andvalues ns). Different size-dependent thermal of the thermal capacitances capacities (symbols) have been have assumed. been used [37,111]. Fixed thermal capacitances (solid lines) and size-dependent thermal capacitances (symbols) have been used [37,111]. Finally, it was previously seen that TM7 deals with two temperatures (at the filament andFinally, at the surrounding it was previouslyinsulator). Note seen that TM7that although deals TM9 with two (following the temperatures twofilament (at the cylinders schema shown in Figure 17) also consider two temperatures, they are both and at the surrounding insulator). Note that although TM9 (following the two cylinders linked to the conductive filament. TM7 can be obtained as a particular case of TM9 considering schema shown in Figure 17) also consider two temperatures, they are both linked to the only one cylinder, but connecting another RC network between the node TCox conductive filament. TM7 can be obtained as a particular case of TM9 considering only and a voltage source for representing the bulk insulator temperature. It is important to note that in the one cylinder, but connecting another RC network between the node TCox and a voltage TM9 case, variable (at simulation time) thermal components are allowed. source for representing the bulk insulator temperature. It is important to note that in the TM9 case, variable 3. General (at simulation Memristor Modeling time) thermal with Framework components Thermal are allowed. Effects Emphasis In this section, a different approach to model RRAM thermal features is proposed. 3. General Memristor Modeling Framework with Thermal Effects Emphasis For this purpose, we make use of the general memristor modeling workspace that was In this section, introduced a different by Corinto and Chuaapproach to This in [71]. model RRAM thermal alternative features perspective is proposed. complements the Fordevelopments presented in the previous section. In particular, the authors [71] that this purpose, we make use of the general memristor modeling workspace was developed introduced a unified by Corinto and theoretical Chua inand framework [71].also This alternative discussed theperspective advantagescomplements of using the the flux- developments presented in the previous section. In particular, the authors charge (ϕ-Q) domain to study memristor elements. Within this framework, memristors, [71] developed a in unified theoretical the taxonomy framework proposed and(the in [72] alsoideal, discussed the advantages the generic, of usingmemristor), and the extended the flux- charge (𝜑-Q) domain to study memristor elements. Within this framework, are described as the result of different approximations in the equations. This extended memristors, incategorization the taxonomy proposed emerged as in a[72] (the ideal, necessity the generic, in order and theoretically to include the extendedthememristor), description areof described as the result of different approximations in the equations. This extended pinched, hysteretic behaviours demonstrated by various elements, not only in circuit theory cate- gorization emerged and electronics but as alsoa innecessity nature. in order to include theoretically the description of pinched, hysteretic Among behaviours the different demonstrated categories presentedby various above, theelements, not onlyclass most general in circuit the-to is linked ory and electronics extended but also memristors, whichin refers nature. to memristors that have extra state variables (in addition to ϕ and Q). For the specific case of charge-controlled memristors, they are described by the following equations: v = M( Q, i, X) i (34) dX = gQ ( Q, i, X) (35) dt dQ =i (36) dt The memristance M of an extended memristor is implied in Equation (34), apparently bearing the feature of nonlinearity; v is the voltage across the memristor, i is the current flowing through it, and Q is the charge, i.e., the current first momentum. The vector X stands for a set of extra state variables, including all the necessary physical magnitudes Nanomaterials 2021, 11, 1261 24 of 46 according to the implemented memristive system; indicatively, they could be the device internal temperature, the conducting filament radius, or any other non-electrical variable influencing the memristor state that ultimately affects the device charge conduction. Ap- parently, the dynamics of the state variables X are governed by gQ (Q,i,X) and Equation (35). It is noted that the importance of the class of extended memristors comes from the fact that all real-world memristor devices known until now are indeed extended memristors. Notice that from Equation (34), we can define the memristance as follows: v dϕ/dt dϕ M ( Q, i, X ) = = = (37) i dQ/dt dQ where ϕ is the voltage first momentum, usually called flux in analogy with the charge. See, however, that a more useful way to obtain this relation is through derivation by assuming that the charge and the flux are related by a function f : ϕ = f ( Q, i, X ) (38) Then, by deriving with respect to time, and using the chain rule, we obtain: dϕ ∂ f ( Q, i, X ) dQ ∂ f ( Q, i, X ) di ∂ f ( Q, i, X ) dX ∂ f ( Q, i, X ) ∂ϕ =v= + + = i= i (39) dt ∂Q dt ∂i dt ∂X dt ∂Q ∂Q The latter part of Equation (39) is obtained under the assumption [71] described below, ∂ f ( Q, i, X ) di ∂ f ( Q, i, X ) dX + =0 (40) ∂i dt ∂X dt The system memory capability under no excitation is determined by a special case of Equation (35), often referred to as the power-off plot (POP) equation; in this case we have i = 0 or equivalently Q = constant. If a single state variable is considered, it is clear that if the POP equation is nil under these conditions, the system presents then a long-term memory since the state variable will not change with time; while if it is different to zero, the system is capable of exhibiting only short-term memory. The consideration of the temperature, T, as the state variable is a peculiar case since there is an influence from external sources (i.e., the ambient –room– temperature, T0 ); this implying a possible energy input in some cases if T0 is not constant. In addition, as it has been discussed in Section 2, it is difficult to determine a single value for the device internal temperature (some models, as shown above, include two different temperatures to better describe the device operation). We can write the equations by separating the temperature as follows: v = M ( Q, i, X, T, T0 ) i (41) dX = gQ ( Q, i, X, T, T0 ) (42) dt dT = gT ( Q, i, X, T, T0 ) (43) dt These equations include both temperatures: the internal device temperature, T, and the external T0 . Obviously, there is no equation governing T0 dynamics since it can be considered as an external signal. Temperature, T, may be a position-dependent temperature T(x,y,z), as already presented in Section 2.1. In this case, Equation (43) would correspond to the heat equation. As an additional note, it is important to highlight that a device will not present long-term memory characteristics associated solely with temperature, since the device will tend to reach thermal equilibrium with the external medium. In absence of any external electrical input, this would mean that the POP equation related to the evolution of the internal temperature is not zero in the general case. This does not preclude, however, that the system may have other internal variables that do present a long-term memory capability. As an example, we can think if the case of a phase change memory (PCM), Nanomaterials 2021, 11, 1261 25 of 46 where a phase change is activated by temperature, and it remains even after the device has cooled back to room temperature. A similar situation occurs when a RRAM conductive filament is ruptured because of an enhanced diffusion process favored by a temperature rise [24,80]. Example of Application As an easy example, we can consider the case of the thermistor, which is one of the first elements identified as an extended memristor, i.e., a device whose resistance is dependent on some internal state variable that presents memory effects [112,113]. In this respect, and for the modeling developments henceforth, they display parallelism with RRAMs. The model has been known since Steinhart and Hart published a function which fitted the variation of thermistor-resistance according to temperature [114] as a Taylor’s expansion of the device conductance in terms of the temperature logarithm. This function, along with the equivalent Ohm’s law in Equation (44), has proved to be suitable in a wide variety of thermistors, for ranges from a few degrees to a few hundred degrees, and it has been widely used to model this kind of devices when used as temperature sensors. The most usual way to describe it is by means of a simplification shown in Equation (45). In addition, a key thermistor characteristic is linked to self-heating, which can be described by Equation (46), by neglecting radiative heat dissipation: v = M ( T, T0 ) i (44)    1 1 M ( T, T0 ) = R0 exp B − (45) T T0 dT R δ = 0 i2 − ( T − T0 ) (46) dt C C In the above equations, T0 is the room temperature, and T the device internal tem- perature. The rest of the symbols are parameters of the thermistor model and can be considered as constants for all practical purposes. At this point, it is noteworthy to point out that these equations bear exactly the same form as Equations (41) and (43) and, thus, identify the thermistor as a memristor. That is, the thermistor is a device whose resistance depends on its electrical history, and it has an internal state variable that governs the overall behaviour (the device internal temperature). Thus, the device can be classified as an extended memristor. In addition, if we look at the POP equation (Equation (46) with i = 0), we see that the temperature derivative is different from zero for any situation other than the device thermal equilibrium with the surrounding environment. Thus, the device does not possess a feature linked to long-term memory. In this respect, it is also illustrative to point out that there are other memristive systems [115] that are also extended memristors due to self-heating, even if the thermal specific mechanisms are different from those of thermistors. As an example to illustrate the memristive behavior of this device, we have simulated it when driven by a triangular current waveform, using specific values extracted from a datasheet for the thermistor constants: δ = 4 × 10−3 W K−1 , C = 60 × 10−3 J K−1 , B = 3950 K, T0 = 298 K, R0 = 10 kΩ. Additionally, and in accordance with a typical thermistor datasheet, we have set a maximum current of 4.5 mA, and we have also used 5 different ramp slopes, as plotted in Figure 22. Figure 23 plots the evolution of the memresistance as well as the current input versus time. The input signal shown in Figure 22 has been employed as well as the Equations (44)–(46). In addition, we have also plotted these two magnitudes used in the Y axis previously against each other in Figure 24, showing the effects of the different slopes in the device memristance. datasheet for the thermistor constants: δ = 4 × 10 W K , C = 60 × 10 J K , B = 3950 K, T0 = 298As K, an R0 example = 10 kΩ. to illustrate theand Additionally, memristive behavior in accordance ofathis with device, typical we havedatasheet, thermistor simulated it when driven by a triangular current waveform, using specific values extracted we have set a maximum current of 4.5 mA, and we have also used 5 different ramp slopes, from a datasheet for the thermistor as plotted in Figure 22. constants: δ = 4 × 10 −3 W K−1, C = 60 × 10−3 J K−1, B = 3950 K, T0 = 298 K, R0 = 10 kΩ. Additionally, and in accordance with a typical thermistor datasheet, Nanomaterials 2021, 11, 1261 we have set a maximum current of 4.5 mA, and we have also used 5 different ramp slopes, 26 of 46 as plotted in Figure 22. Figure 22. Input current waveform versus time for five different ramps. Colours are coherently employed in the following plots. Figure Figure Figure 22. 23 plots 22.Input Input thewaveform current currentevolution waveform ofversus thetime versus memresistance for five time as well fordifferent asramps. ramps. five different the current Colours input areversus are coherently Colours coherently time. The employed input in signal the followingshown plots. in Figure 22 has been employed as well as the Equations employed in the following plots. (44)–(46). Figure 23 plots the evolution of the memresistance as well as the current input versus time. The input signal shown in Figure 22 has been employed as well as the Equations (44)–(46). Nanomaterials 2021, 11, 1261 28 of 47 In addition, we have also plotted these two magnitudes used in the Y axis previously against Figure Figure23. each 23. other in Figure Memristance Memristanceversus 24, showing versustime timefor forfive the effects fivedifferent different inputof input the different voltage voltage signals slopes signals under under in thevoltage ramped ramped devicestress. voltage memristance. stress. Figure 23. Memristance versus time for five different input voltage signals under ramped voltage stress. Figure24. Figure Memristanceversus 24.Memristance versus input input current, current, forfor five five different different slopes. slopes. Colours Colours are are coherent coherent withwith the the results results shown shown in in thethe previous previous figures. figures. Then, the typical i-v memristor characteristics showing the famous loop are drawn in Figure 25. See that the shapes are in line with what is expected for memristors [71]. Nanomaterials 2021, 11, 1261 Figure 24. Memristance versus input current, for five different slopes. Colours are coherent with 27 of 46 the results shown in the previous figures. Then, the typical i-v memristor characteristics showing the famous loop are drawn in See that Figure 25. See that the shapes are in line with what is expected for memristors [71]. i-v thermistor Figure 25. i-v thermistor characteristics for five different input voltages. Figure 25 Figure 25 shows showsthe thetwo twofingerprints fingerprintsofofa amemristor memristor [116]: [116]: (1)(1) a pinched a pinched loop; loop; andand (2) (2) an area that varies with frequency, tending to a line at high frequencies. an area that varies with frequency, tending to a line at high frequencies. In fact, the behav- In fact, the behaviour iour at the at the highest highest frequencyfrequency (cyan (cyan line) isline) nearlyis nearly that ofthatan of an ohmic ohmic resistor, resistor, with with no no loop loop area, while at lower frequencies the behaviour is different. area, while at lower frequencies the behaviour is different. It has to be noted that (see It has to be noted that (see [115]) [115]) the situation the situation is more is more interesting interesting thanthan simply simply anan areaincrease area increaseatatlower lowerfrequencies. frequencies. As highlighted in [115], the thermistor control variable is the internal As highlighted in [115], the thermistor control variable is the internal temperature and temperature andit it tends to be in thermal equilibrium with the surroundings, which tends to be in thermal equilibrium with the surroundings, which causes the loop area to causes the loop area to reduce reduce at aatvery a very lowlow frequency. frequency. As As we we useuse it initthe in the corresponding corresponding equations, equations, we can we can see see that for very low input signal slopes, the device always reaches thermal equilibrium, that for very low input signal slopes, the device always reaches thermal equilibrium, since since the dT/dt is nearly zero, and its behaviour tends to be close to a nonlinear resistor, as shown in Figure 24; Figure 25 (red lines), with a null area enclosed in the loop. At this point, the concept of dynamic route map (DRM) [36] comes up since it is quite useful to represent the device time evolution (again, in this facet, a full parallelism with RRAM is observed [36]). In fact, the DRM is a concept arising from non-linear dynamical systems, and represents the trajectories a system follows in the phase space of the state variable versus its derivative for different control parameter values. If we plot it as a 3D diagram, we find that all the trajectories, for a given constant T0 , are bound to fall on the same surface, as seen in Figure 26. Using this representation may provide very interesting insights into the device dynamics. Considering our thermistor, if a trajectory goes from a state with positive derivative to another characterized by a negative derivative with increasing temperature, then it will reach a stable equilibrium point at the temperature where the derivative nullifies. An equilibrium point is a state where the device tends to remain at even if it drifts from it in its operation; this idea resembles a similar concept related to the DC quiescent point in circuit theory, or memory in memristors. In the opposite case, when the trajectory goes from negative to positive, an equilibrium point might seem to come up at the zero-crossing temperature, but it is unstable, which means that the slightest change will force the system to come out of it. where the derivative nullifies. An equilibrium point is a state where the device tends to remain at even if it drifts from it in its operation; this idea resembles a similar concept related to the DC quiescent point in circuit theory, or memory in memristors. In the op- posite case, when the trajectory goes from negative to positive, an equilibrium point might Nanomaterials 2021, 11, 1261 seem to come up at the zero-crossing temperature, but it is unstable, which means that28 the of 46 slightest change will force the system to come out of it. Figure 26. Memristor Dynamic Route Map (surface), showing as lines the five trajectories corre- Figure 26. Memristor Dynamic Route Map (surface), showing as lines the five trajectories corre- spondingtotothe sponding theprevious previousfigures, figures, using using the the same same colour colour code. code. ItIt can can bebe seen seen that that all all these these trajectories trajecto- fallfall ries on on thethe surface, which surface, defines which univocally defines the the univocally device behaviour. device behaviour. Notice thatthat Notice this this surface is, inis,fact, surface ina fact, family of surfaces a family that depend of surfaces on theon that depend room the temperature T0 . room temperature T 0. 4. RRAM Quantum Point Contact Modeling, Thermal Effects 4. RRAM Quantum Point Contact Modeling, Thermal Effects So far, we have studied thermal effects from a classical physics viewpoint. This So far, we have studied thermal effects from a classical physics viewpoint. This is the is the approach commonly used in most compact models. However, the scale of the approach commonly used in most compact models. However, the scale of the material material layers that form part of a RRAM makes feasible for certain devices and for layers that form part of a RRAM makes feasible for certain devices and for particular op- particular operation conditions the observation of quantum effects. In particular, when eration conditions the observation of quantum effects. In particular, when the conducting the conducting filament (CF) cross-sectional area in a RRAM device becomes very narrow, filament (CF) cross-sectional area in a RRAM device becomes very narrow, only a few only a few atoms wide, the continuum description is expected to break down, so that atoms wide, the continuum description is expected to break down, so that the direct ap- the direct application of the heat equation becomes questionable [117,118]. We enter plication of the heat equation becomes questionable [117,118]. We enter into the regime of into the regime of heat transport and dissipation at the nanoscale [119]. Before starting heat transport and dissipation at the nanoscale [119]. Before starting with a discussion with a discussion about these topics, it is important to discriminate between the role that about these topics, it is important to discriminate between the role that temperature plays temperature plays on the ion/vacancy movement across the insulating films (essentially on the ion/vacancy governed by Kramers’movement across and theory [120]), the insulating films (essentially the temperature dependencegoverned by Kra- of the mechanism mers’ theory [120]), and the temperature dependence of the mechanism adopted adopted for the electronic transport in the CF itself (Schottky, Poole-Frenkel, tunneling, for the electronic transport variable range in the hopping, CF charge space itself (Schottky, Poole-Frenkel, limited conduction, quantumtunneling, variable point contact, etc.)range [121]. hopping, space charge limited conduction, quantum point contact, etc.) [121]. Although both descriptions are intrinsically connected and they are at the heart of the Although both descriptions complexity of theare intrinsically RRAM connected behavior, they are and they often are atseparately treated the heart of forthe complexity simplicity, of or one of them completely dropped. To deepen into these intertwined approaches, a detailed RRAM dynamic simulation at the microscopic level is imperative. Ion/vacancy diffusion process, which defines the filamentary structure, depends on temperature, and the size of the structure determines the magnitude of the electron current that governs the power dissipation (Joule heating effects), which in turn affects the local temperature that drives the diffusion process. In this section we will exclusively focus the attention on the modeling of the electron transport at the nanoscale and the influence of temperature and power dissipation on it. The role of ions/vacancies will be indirectly addressed (more on this issue is explained in Section 5). Temperature evolution in the structures under study (e.g., Figure 1; Figure 2) is extensively covered in Section 2. First, a fixed ion/vacancy arrangement will be considered for simplicity, so that two particular extreme cases, LRS and HRS, will be examined. Second, a phenomenological model for the transition from HRS to LRS which takes into account the power dissipated at the CF bottleneck will be presented. As already discussed previously, while in LRS, the CF is completely formed establishing a metallic-like connection in between the electrodes; in HRS, the filament presents a gap along its structure as a consequence of the absence of conduction states [122]. Thereby, a common theory for both cases able to demonstrate consistency with the experimental observations is required. Mesoscopic physics, a subdiscipline of condensed-matter physics, has resulted in a suitable framework for semi-empirically describing the temperature dependence of the Nanomaterials 2021, 11, 1261 29 of 46 electron flow in narrow constrictions and, in particular, in RRAMs. Mesoscopic physics describes the conducting properties of systems whose size lays in between the macro- scopic (bulk material) and the microscopic (atoms and molecules) worlds. Since we are talking about atomic dimensions, quantum mechanics is at the foundations of meso- scopic physics [123]. Within this framework, the quantum confinement associated with filamentary conduction is described in terms of the electrochemical potential, potential wells, potential barriers, and bands (valence, conduction, gap). We also talk about a semi- empirical approach because the local temperature is frequently unknown and the external temperature is considered as a control parameter. This objection can be partly overcome if models including two different temperatures (for the cold and hot part of the CF, and for the CF and its surrounding region, see Sections 2.3.4 and 2.4.3) are employed. The idea of using a mesoscopic approach is the consequence of the observation of experimental RRAM conductance values around integer and non-integer multiples of the quantum conductance unit G0 = 2 e2 /h, where e is the electron charge and h the Planck’s constant [124]. In terms of resistance, this unit is R0 = 1/G0 = 12.9 KΩ. The experimental conductance values for many cycles measured at a fixed bias, or the conductance measured at consecutive steps in one cycle at different or constant biases are often displayed using histogram plots with the x-axes normalized to G0 . In many cases, these histograms reveal a peak structure which is interpreted as an indicator of the number of channels available for conduction or as the occurrence of preferred atomic configurations for the CF [125]. Although the detection of peaks in the device histograms is recognized as the signature of quantum point-contact conduction, it should be taken into account that measurements can be seriously affected by a number of factors such as the existence of multiple conduction paths, series resistance, roughness and scattering caused by the granularity of matter, in general, non-adiabatic (non-smooth) potential profiles. Caution should also be exercised with the use of the term conductance quantization: only for simple s-electron metals, the transmission probability for the conductance channels is expected to open close to integer values [126]. For this reason, observations in the field of RRAM should be more appropriately considered to be in the quantum (rather than in the quantized) regime of conductance [125]. Many experimental results on resistive switching materials have been interpreted in terms of conduction through atom-sized filamentary structures [127]. This is the case of a wide variety of binary and ternary oxides such as SiOx [128–131], HfO2 [132–139], Ta2 O5 [140–142], NiO [143,144], ZnO [145,146], a-Si:H [147], TiO2 [148], V2 O5 [149], YOx [150], and BiVO4 [151]. Nonlinear effects in HfO2 were also reported by Degraeve et al. [152] and in CeOx /SiO2 -based structures by Miranda et al. [153]. From the point of view of theory, it is worth mentioning that the CF formation in monoclinic- and amorphous-HfO2 was investigated from first principles by Cartoixa et al. [154] and by Zhong et al. [155]. The filamentary paths are built from oxygen vacancies and using a Green’s function formalism coupled to a density functional theory code, the conductance of filaments of different lengths was calculated. According to the obtained results, even the thinnest CFs can sustain conductive channels exhibiting signs of quantum conduction. Very often, LRS is associated with conductance values G ≥ G0 and with a linear I-V curve (not to be confused with Ohmic behavior). In this case, the device conductance can reach values from 10 to 100 times G0 which indicates the large number of atoms partici- pating in the filament formation. On the other hand, HRS is associated with conductance values G < G0 and with a non-linear I-V curve (mainly with exponential behavior). This state is characterized by a gap or potential barrier which acts as a blocking element for the electron flow. As the starting point for the inclusion of the thermal effects in RRAMs, the Buttiker-Landauer approach for quantum point contacts is considered [156]. Importantly, the analysis does not discriminate between CBRAMs and OxRAMs, so they are treated on equal grounds. curve (not to be confused with Ohmic behavior). In this case, the device conductance can reach values from 10 to 100 times G0 which indicates the large number of atoms partici- pating in the filament formation. On the other hand, HRS is associated with conductance values G < G0 and with a non-linear I-V curve (mainly with exponential behavior). This Nanomaterials 2021, 11, 1261 state is characterized by a gap or potential barrier which acts as a blocking element for the 30 of 46 electron flow. As the starting point for the inclusion of the thermal effects in RRAMs, the Buttiker-Landauer approach for quantum point contacts is considered [156]. Importantly, the analysis Accordingdoesto notthe discriminate finite-bias between Landauer’s CBRAMs formulaand[157], OxRAMs, so they the I-V are treatedofona characteristic equal grounds. mesoscopic conductor can be expressed as: According to the finite-bias Landauer’s formula [157], the I-V characteristic of a 2e be expressed as: Z mesoscopic conductor I= can D ( E)[ f ( E − βeV ) − f ( E + (1 − β)eV )]dE (47) h 2𝑒 𝐼 = 𝐷(𝐸) 𝑓(𝐸 − 𝛽𝑒𝑉) − 𝑓(𝐸 + (1 − 𝛽)𝑒𝑉) 𝑑𝐸 (47) where E is the energy, D ℎthe tunneling probability, f the Fermi-Dirac (FD) distribution function, and 0 ≤ β ≤ 1 the fraction of the applied voltage that drops at the source side where E is the energy, D the tunneling probability, f 1the Fermi-Dirac (FD) distribution of the constriction. For a symmetrical structure β = 2 . Assuming an inverse parabolic function, and 0 ≤ β ≤ 1 the fraction of the applied voltage that drops at the source side of potential barrier for the constriction bottleneck, D is given by [158]: the constriction. For a symmetrical structure β = ½. Assuming an inverse parabolic poten- tial barrier for the constriction bottleneck, D is given by [158]: D ( E) = {1 + exp[−ν( E − ϕ)]}−1 (48) 𝐷(𝐸) = 1 + 𝑒𝑥𝑝 −𝜈(𝐸 − 𝜑) (48) where ν is a coefficient related to the curvature of the potential barrier and ϕ the height of where ν is a coefficient the potential barrier that related to thethe represents curvature of theeffect confinement potential barrier27). (see Figure andFor φ the T =height of 0 K, (47) the potential barrier that represents and (48) yield [159] (see Figure 27): the confinement effect (see Figure 27). For T = 0 K, (47) and (48) yield [159] (see Figure  27): 1 + exp[ν( ϕ − βeV )]   1 I (V ) = G0 V + 1ln 1 + 𝑒𝑥𝑝 𝜈(𝜑 − 𝛽𝑒𝑉) (49) 𝐼(𝑉) = 𝐺 𝑉 +eν 𝑙𝑛1 + exp[ν( ϕ + (1 − β)eV )] (49) 𝑒𝜈 1 + 𝑒𝑥𝑝 𝜈(𝜑 + (1 − 𝛽)𝑒𝑉) Figure Figure 27. 27. Schematic Schematicof ofthe theenergy energystructure structureof ofthe theconducting conductingfilament. filament.InInLRS LRS(high (highcurrent), current),the the CF CF isiscompletely completelyformed formedandandthe theconfinement confinementpotential potentialbarrier is low. barrier In HRS is low. (low In HRS current), (low the the current), filament filament isis broken broken and and the the confinement confinement potential potential barrier barrierisishigh. high.The Thegreen greenarrows arrowswidth widthindicates indicates the electron current magnitude. the electron current magnitude. Figure Figure 28 28 shows shows some some typical modeling results using Equation (49). Any additional potential potential drop along along the the confinement confinementstructure structurecan canbebeaccounted accountedforfor using using thethe transfor- transforma- mation V → tion V →V-I RS in (49), where RS is a series resistance. Equation (49) can be modified so so V-I R S in (49), where RS is a series resistance. Equation (49) can be modified as to include many parallel conducting channels [138]. For LRS, we can consider that there is no blocking element along the CF so that assuming ϕ→−∞ (D→1) in (49), we obtain: I (V ) = G0 V (50) which is the celebrated Landauer formula for a monomode ballistic conductor [160]. Fol- lowing [134], the temperature dependence can be introduced into (50), assuming RS (T) = RS0 ·[1 + αT (T − T0 )], where RS0 = RS (T0 ), αT is a temperature coefficient, and T0 the room temperature. In this case, the I-V characteristic still follows a linear relationship but with a lower slope given by: G0 I (V ) = V (51) 1 + G0 RS ( T ) Nanomaterials 2021, 11, 1261 33 of 47 Nanomaterials 2021, 11, 1261 31 of 46 10-4 a) T b) T LRS LRS 10-5 T Current [A] T HRS 10-6 HRS αΤ=3E-2K-1 200K 200K ϕ0=1eV 250K 250K 10-7 ν=3eV-1 300K 300K θ=1E-3eVK-1 350K 350K 0.0 0.5 1.0 1.5 2.0 0.01 0.1 1 Voltage [V] Voltage [V] Figure 28. Effects of the temperature on the HRS and LRS I-V characteristics. (a) log-linear axis and Figure (b) 28. log-log axis. Effects of the temperature on the HRS and LRS I-V characteristics. (a) log-linear axis and (b) log-log axis. If αT is a positive coefficient, as expected for a metallic-like conductor, the current Todecreases concludeasthis thesection, temperature increases it is worth (see Figure mentioning that 28). This behavior according is in agreement to the standard theory with the experimental of mesoscopic devices, observations heat largely [134]. Notice dissipates at that here the emphasis the electrodes is putand (reservoirs) on the connection thermal of the ballistic region with the rest of the device (internal or external) and electrical conductances are proportional through the Wiedemann-Franz law [119]. and in particular the contacts. Nevertheless, R This idea is to heuristically explain why a CF of atomic dimensions is able to reach a sta- factor with S can also be viewed as the momentum relaxation tionaryalong the state current filamentary structure. If values with conductance we move of to thethe opposite order of G0limit, for that . Notice HRS,theand we consider current specifically the case E << ϕ, (48) reads: density flowing through a nanoscale CF can be extraordinarily high. The question can be summarized as, where is power dissipated in a RRAM system exhibiting quantum prop- D ( E) ≈ exp[ν( E − ϕ)] (52) erties? This is a fundamental question in mesoscopic physics [123]. Let us consider here the progressive so that (47)increase can be of the current integrated flowinto taking as aaccount function theoftemperature-dependent time when the device issmearing sub- of jected the to aFDconstant voltage stress after electroforming. distributions at the contacts. The result is [161]: This process corresponds to the transition HRS→LRS which arises because of the CF widening. Following [162], we can write first the following phenomenological2e exp(−νϕequation ) for the current evolution: I (V ) ≈ {exp[νβeV ] − exp[−ν(1 − β)eV ]} (53) hν sinc(πνkT 𝑑𝐼 ) = 𝜂𝑃 (54) which provides the exponential behavior 𝑑𝑡 observed for HRS. Now, notice that the tempera- where ture η isappears explicitlyand a temperature- in (53) through the sinc function. material-dependent couplingIncoefficient this case, the andcurrent increases PC is the powerwith the temperature dissipated as expected at the constriction from the bottleneck. Foravailability the simplest ofcase more ofenergetic a constantelectrons applied at the bias V,injecting Equationelectrode. (54) expressesHowever, it can that the be shown current levelsthat off for a set in the of typical long fittingpower run because parameters, the smearing dissipation of the FD first increases andfunctions is not enough then progressively to account transfers fromfor thethe observed temperature constriction to the effects electrodes. in HRS.according Second, In order to tocircumvent Landauer’sthis problem, formula, the atransmission new parameter is introduced probability D (av- into the erage) model can be Equation expressed(53) as aso that a larger function of thevariation of the current current flowing throughcan thebe achieved. structure as: Following experimental observations for the soft breakdown conduction mode in SiO2 [161], the confinement potential barrier𝐷height = 𝐺 ϕ𝐺can = be(𝐺parameterized 𝑉) 𝐼 as ϕ(T) = ϕ0 -θ(T − T(55) 0 ), where ϕ = ϕ(T and the 0power ) and θ > 0 is a linear temperature coefficient. This correction 0dissipated at the constriction can be calculated from the voltage drop V term arises from C the thermal movement of occurring at the constriction using: ions/vacancies in the CF around their equilibrium positions. In this case, as the temperature increases, the tunneling current increases because of the reduction 𝐼 of the effective barrier 𝑃 = height 𝑉 𝐼 (see = 𝑉Figure1 − 𝐷 28). 𝐼 =The𝑉𝐼 temperature 1− effect on the barrier (56)profile was recently investigated in detail in [139] using inverse𝐺modeling 𝑉 in combination with the WKB approximation for the tunneling probability. To conclude this section, it is worth mentioning that according to the standard theory of mesoscopic devices, heat largely dissipates at the electrodes (reservoirs) and thermal Nanomaterials 2021, 11, 1261 32 of 46 and electrical conductances are proportional through the Wiedemann-Franz law [119]. This idea is to heuristically explain why a CF of atomic dimensions is able to reach a stationary current state with conductance values of the order of G0 . Notice that the current density flowing through a nanoscale CF can be extraordinarily high. The question can be summarized as, where is power dissipated in a RRAM system exhibiting quantum properties? This is a fundamental question in mesoscopic physics [123]. Let us consider here the progressive increase of the current flow as a function of time when the device is subjected to a constant voltage stress after electroforming. This process corresponds to the transition HRS→LRS which arises because of the CF widening. Following [162], we can write first the following phenomenological equation for the current evolution: dI = ηPC (54) dt where η is a temperature- and material-dependent coupling coefficient and PC is the power dissipated at the constriction bottleneck. For the simplest case of a constant applied bias V, Equation (54) expresses that the current levels off in the long run because power dissipation first increases and then progressively transfers from the constriction to the electrodes. Second, according to Landauer’s formula, the transmission probability D (average) can be expressed as a function of the current flowing through the structure as: e = G −1 G = ( G0 V )−1 I D (55) 0 and the power dissipated at the constriction can be calculated from the voltage drop VC occurring at the constriction using:     I PC = VC I = V 1 − D I = V I 1 − e (56) G0 V Then, from expression (54), an explicit differential equation for the current evolution is obtained:   dI I = ηV I 1 − (57) dt G0 V Equation (57) is nothing but the logistic equation with effective transition rate ηV and carrying capacity G0 V. The solution to Equation (57) for a constant bias reads: G0 I0 Vexp(ηVt) I (t) = (58) I0 [exp(ηVt) − 1] + G0 V which complies with I(t = 0) = I0 and I(t = ∞) = G0 V, the initial and stationary conditions, respectively. Equation (58) expresses that, when a mesoscopic channel with conductance G0 is formed, the power fundamentally dissipates at the electrodes and not at the constriction’s bottleneck (see Figure 29). Power is indeed dissipated at the constriction during the CF formation as discussed in the next section. Of course, this is a simplistic view of a much more complex process. respectively. Equation (58) expresses that, when a mesoscopic channel with conduct G0 is formed, the power fundamentally dissipates at the electrodes and not at the striction’s bottleneck (see Figure 29). Power is indeed dissipated at the constriction du the CF formation as discussed in the next section. Of course, this is a simplistic view Nanomaterials 2021, 11, 1261 much more complex process. 33 of 46 0.14 a) η=10(Vs)-1 I0=1E-7A b) V=1V I0=1E-7A 0.08 0.12 P P Dissipated power [mW] 0.10 0.06 0.08 η 0.04 10(Vs)-1 0.06 25(Vs)-1 50(Vs)-1 0.04 V PC PC 0.75V 0.02 0.02 1.00V 1.25V 0.00 0.00 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 Time [s] Time [s] Figure 29. Evolution of the power dissipated in the structure P (solid lines), and at the constriction PC (dashed Figure 29. Evolution lines). of the power(a) Corresponds dissipated in to thedifferent structureapplied voltages P (solid V, and lines), and (b) constriction at the correspondsPtoC (dashed differentlines). (a) Corresponds transition rates to different applied η. voltages V, and (b) corresponds to different transition rates η. Although most of the models Although most addressed in this of the models manuscript addressed are devoted in this to devices manuscript that to de are devoted show single filamentary conduction, some of them could also be applied to area-dependent that show single filamentary conduction, some of them could also be applied to area devices and even to devices pendent devicesthat anddo notto even require deviceselectroforming [1,2,70]. that do not require In particular,[1,2,70]. electroforming the In pa Landauer’s approach can be extended to area-dependent devices by assuming multifila- ular, the Landauer’s approach can be extended to area-dependent devices by assum mentary conduction. This is theconduction. multifilamentary case when the ThisLandauer formula is the case when the(49)Landauer includes aformula prefactor(49) includ N dealing with the number of identical filaments assumed [138,159]. In addition, a model- prefactor N dealing with the number of identical filaments assumed [138,159]. In addi ing procedure in line with the general memristor framework presented in Section 3 could a modeling procedure in line with the general memristor framework presented in Se also be possible [71,72]. 3 could also be possible [71,72]. 5. Thermometry of Conducting Filaments In addition to the developments described above in relation to the HE solutions in different modeling approaches and levels of complexity, it is interesting to understand the dielectric breakdown (BD) phenomenon in the context of RRAM operation. In this respect, we shed light in this section on the structural damage that occurs during the BD current transient. In ultra-thin dielectric layers (<5 nm), there is a wide consensus around considering that the intrinsic BD is related to the generation of defects in the dielectric film [163–166]. When the density of bulk defects is high enough, the BD event is triggered by the local connection of the electrodes through a defect related conduction path. Once a defect percolation path is formed, the main feature is a progressive increase of the current across the dielectric. This phenomenon, often referred to as progressive breakdown (PBD), is an universal process that lies under a wide variety of dielectric materials, ranging from traditional oxides, such as SiO2 , SiOX NY [167] to innovative 2D dielectrics, such as h-BN [168], passing through high-k materials such as Al2 O3 and HfO2 [169]. The physical structure of the filament formed during the PBD regime has been deeply studied. In the case of poly-Si/SiON/Si MOS devices during PBD, it was demonstrated that the filament is, at least in part, made of Si atoms, through the mechanism of dielectric breakdown induced epitaxy (DBIE) [170,171]. The filament sizes were directly observed by scanning transmission electron microscopy (STEM) after the BD of MIM structures with either Ti/HfO2 /TiN or with Hf/HfO2 /TiN devices, in which the top electrode (Ti Nanomaterials 2021, 11, 1261 34 of 46 or Hf) acted as cathode. Clear evidence of the formation of a metallic filament made of, respectively, Ti or Hf was reported by using electron energy loss spectroscopy (EELS) imaging [172,173]. In the case of 2D h-BN (CVD) dielectric layers, there is strong evidence that the CFs are formed by metal ions that penetrate from the electrodes into the h-BN stack under the action of the electric field [174,175]. Different experiments have probed that the SET event in RRAM devices [176,177] and the dielectric BD of gate oxides [169,178] have some common aspects. In addition to the clear dependence of the CF characteristics on the maximum current flowing through the device [18,179], TEM imaging of Si-based MOS capacitors prior to and post dielectric BD [163,170,180] and HfO2 -based RRAM cells after forming and cycling [172,181] show comparable microstructural changes in the oxide, suggesting the diffusion of the anodic atomic species into the oxide layer in both cases. Thus, these two phenomena share not only similar electrical characteristics, but also generate comparable microstructural changes, suggesting a common underlying physical mechanism. In such scenario, we propose to model the results for the SET event in RRAM devices similarly to the gate-oxide BD in MOSFETs. The PBD effect has been captured by the model proposed by Palumbo et al. in [169,182], and later expanded by Lombardo in [183], clarifying the primary role played by the carrier energy loss through the PBD spot. Such model accounts for the physics behind the progressive evolution of the current, where the BD process is closely linked to the energy transfer from the CF itself to its surrounding atomic lattice. According to this idea, the high temperature associated with the localized current flow (being the BD spot area of 1–50 nm2 , the current density can reach a few MA/cm2 [184–186]) would contribute to the generation and enlargement of the BD filament connecting the electrodes of the stack, enabling the promotion of the electro-migration of the fastest available atomic species. Since this technique unambiguously relate the transition rate (dITr /dt) to the heat dissipation properties during the atomic diffusion of the cathode or anode atoms into the gate dielectric in the region of the percolation path, it is possible estimate the CF temperature. Considering the model reported in [169], we can express the current transition rate (marked as TR) as: dITr q V f1 TR = = DISET (59) dt k B T t2ox where T is the temperature of the CF, tox is the dielectric thickness, kB is the Boltzmann constant, D is the diffusion constant of the atomic species responsible for the generation of CF, ISET is the current level at the onset of the transition, and f 1 = ne λe σe , with ne being the electron density, λe the electron mean free path and σe the cross-section for the electron- atom collision (responsible for the momentum transfer). V is the applied voltage across the BD spot which has been assumed to be equal to the overall externally applied bias between the metal contacts of the stack. f 1 value is around the unity since the defect concentration in the CF is most likely very high [172]. According to Equation (59), dITr /dt is proportional to D × ISET . This means that the BD growth rate rises either by increasing the dominant diffusivity D of the fastest atomic species or by increasing the charge carrier flux. The I-V characteristics under the PBD regime can be explained by some well-known physical models, for example invoking trap-assisted tunneling (TAT) current [187], co- tunneling [188,189], and the quantum point contact model [159]. Although the transport properties of stacks with different materials are fitted by considering different transport models, the underlying concept is similar in all cases, the electrons passing through the PBD spots experience a very large energy loss. To simplify the approach, the BD spot ISET -V curve is usually modeled by assuming a simple analytical dependence as described in [190]. It is important to mention that independent experiments have probed that power dissipation taking place inside the dielectric layer is a reasonable assumption. According to Takagi [191], electrons tunneling through defects responsible for stress induced leakage current (SILC) in thin oxynitrides loose a fraction of energy, and as shown by Blochl and Stathis [192], this is caused by defect relaxation. It is reasonable to assume that a Nanomaterials 2021, 11, 1261 35 of 46 similar effect takes place for electron transport through the BD spot, since there is a clear dependence of the CF characteristic on the maximum current flowing through the device, both for the BD of gate oxides [179] and the SET event in RRAM devices [18,163,193], as well as for layered dielectrics, such as h-BN [168,190]. A simplification of spherical symmetry around the CF can be assumed to model the temperature in the BD spot. Within this approach, the temperature can be described by Equation (60), where the dissipated power at the CF is proportional to ISET × V, kth is the thermal conductivity of the dielectric, T0 is the room temperature and f 2 is the fraction of the energy lost at the constriction: f 2 V IBD T= + T0 (60) 2 π tox k th The fitting parameters are f 1 , f 2 , D0 , and Eact ; (taking into consideration that D = D0 *exp(-Eact /kB T)) and they describe the main features of the progressive BD effect on different stacks [190]. In this section we considered a RRAM device based on a MIM stack with a 10 nm thick atomic layer deposited HfO2 film sandwiched between Ti and TiN electrodes [193]. During the HRS to LRS transition the current (ITr ) increases gradually with time evidencing the progressive nature of the SET event (see Figure 30a,b). It is a noisy and progressive process well in agreement with the literature [176–178,194] whose duration shows a strong voltage dependence and dispersion. The time evolution of the HRS to LRS transition is Nanomaterials 2021, 11, 1261 37 of 47 quantified by the slope dITr /dt, as defined in [169,185,195]. TR values were experimentally evaluated through measurements such as those shown in Figure 30a,b (approximately 100 measurements for each voltage value). Figure 30. Figure Currenttransient 30. Current transientof ofthe theDUTs DUTs under underconstant constantvoltage voltagestress stress(CVS). (CVS).(a) (a)represents representsthe thelowest lowest–450 –450mV– mV– voltage voltage whereas (b) whereas (b) the the highest highest –650–650 mV–. mV–. Ball Ball marker marker 1 points out the initial current (IInit Init ), ), whereas whereas 22 the the onset onset of of the the progressive progressive increase increase of of current current (I (IOn On ) )and and33the thefinal finaljump jumptotothe thecompliance compliancelevel level(I(IEnd). ). End Thefitting The fittingresults resultsforforTRTR as as a function a function of the of the applied applied voltage, voltage, obtained obtained with with the the pro- proposed model in Equations (59) and (60), are illustrated in Figure 31. posed model in Equations (59) and (60), are illustrated in Figure 31. The proposed fit ac- The proposed fit accounts for both the t counts for both the tox ox reduction (t reduction (tox ox considered is equal to considered is equal to tgapgap t ~2 nm due to the ~2 nm due to the forming forming step)and and the the increase increase in diffusivity −6 cm22 /sec as other species are step) diffusivity (D (D00 is is in in the theorder orderof of~10 ~10−6 cm /sec as other species are consideredto considered tocomplete completethe theCF, CF, i.e., i.e.,oxygen oxygenvacancies). vacancies). The Therest restof ofthe theparameters parametersinvolved involved remainas remain aspreviously previouslymentioned mentioned in in the the literature literature (E(Eact eV,f2f~0.1 ~0.3–0.7eV, act~0.3–0.7 2 ~0.1 f 1 ~1) andf1~1) and [193]. [193]. The fitting results for TR as a function of the applied voltage, obtained with the pro- posed model in Equations (59) and (60), are illustrated in Figure 31. The proposed fit ac- counts for both the tox reduction (tox considered is equal to tgap~2 nm due to the forming step) and the increase in diffusivity (D0 is in the order of ~10−6 cm2/sec as other species are considered to complete the CF, i.e., oxygen vacancies). The rest of the parameters involved Nanomaterials 2021, 11, 1261 36 of 46 remain as previously mentioned in the literature (Eact~0.3–0.7 eV, f2~0.1 and f1~1) [193]. Figure Figure TRTR 31.31. data data (squaremarkers (square markers–ball –ballmarkers markersrepresent representthe the mean mean value–) fitted fitted assuming assumingoxygen oxy- vacancies gen vacancies andand = =toxtox(cyan tgaptgap (cyandashed dashedline line– curve NN° – curve ◦ 1).1). Additionally, Additionally,TRTR assuming assuming literature values literature values for toxfor andtoxDand D is plotted is plotted –curves N◦ 2, 3N° –curves 2, 34–. and and TR4–. TR presents presents a strong a strong dependence dependence with ap- with applied voltage, plied voltage,almost increasing increasing almost one order of one order offor magnitude magnitude every 50 mVfor every step. 50 mV step. Nanomaterials 2021, 11, 1261 38 of 47 ToToimplement implementthis thismodel modelforforthe theSET SETevent, event,the theeffect effectofofthe thetoxtoxreduction reductionafter afterthe the forming formingstep stepwas wasconsidered. considered.First, First, aa forming forming operation (a (a controlled controlleddielectric dielectricBD) BD)creates cre- thethe ates CFCF through through the thefresh freshoxide oxidelayer.layer.Then, Then, the the switching mechanism mechanismisisdriven drivenby bythe the creation creation of a gap (RESET) and and the restoration restoration ofof the CF (SET). (SET). In the case under study, it has been demonstrated using the statistics of of the the SET SET switching switching time time (t(tSet Set) (i.e., the time to complete the HRS to LRS transition) that tgap gap≈2≈nm2 nm is is a reasonable a reasonable value value [193]. [193]. Figure Figure 32 presents the temperature calculations as a function of voltage provided Equation (60). f22 represents the fractionfraction of of energy energylostlostby bythe thecarriers, carriers,which whichrangesrangesfrom from0 0toto1.1.This Thisparameter parameter also depends also depends ononthethe temperature, temperature, mainly mainlybecause because of phonon-electron of phonon-elec- scattering tron [183]. scattering Therefore, [183]. f2 isf2 aisfunction Therefore, a function ofofvoltage voltageand andtemperature temperaturewhose whose behavior is found by a best fitting procedure. The influence of f2 on the temperature is shown in Figure 32 for different f2 values. values. Temperature estimation according Figure 32. Temperature according to to Equation Equation (60) (60) as as function function of of voltage voltage and and the the energy en- ergy loss loss in theinCF. theThe CF.red Theshaded red shaded zone zone indicates indicates the voltages the voltages employed employed forCVS for the the (0.45 CVS to (0.45 0.65toV). 0.65 V). The best fitting diffusivity required to reproduce the TR vs. voltage is observed in Figure 33. The data have been calculated assuming that D0 is in the order of 3 × 10−6 cm2/s and Eact = 0.52 eV as indicated in [193]. In HfO2-based RRAM devices, the SET event is explained as the completion of the gap due to the migration of O2-ions through a field-assisted and thermally activated effect, which creates the oxygen vacancies that fill the gap along the Nanomaterials 2021, 11, 1261 37 of 46 The best fitting diffusivity required to reproduce the TR vs. voltage is observed in Figure 33. The data have been calculated assuming that D0 is in the order of 3 × 10−6 cm2 /s and Eact = 0.52 eV as indicated in [193]. In HfO2 -based RRAM devices, the SET event is explained as the completion of the gap due to the migration of O2 -ions through a field- assisted and thermally activated effect, which creates the oxygen vacancies that fill the gap along the CF [38,44,53,195,196]. This is quite a relevant point to notice, as the diffusivity of Nanomaterials 2021, 11, 1261 39 of 47 oxygen vacancies (OVs), in a HfO2 layer of thickness like tgap spread over a range similar to the fitted diffusivity for the TR [197] (see Figure 33). diffusivity for Figure 33. Reported values of diffusivity for different different atomic atomic species species vs. vs. reciprocal reciprocalof oftemperature. temperature. diffusivity D The diffusivity D required requiredforforTR TRfitting fittingisisshown showntotobebeininthe thesame samerange rangeasas thethe OVsOVs diffusivity. diffusivity. OVs OVs diffusivity diffusivity [197] [197] is ~10is 4 times ~104 times higherhigher than forthan thefor the metallic metallic speciesspecies as Cu:SiO as Cu:SiO [198],2 [198], Ag:SiC,Ag:SiC, Ag:Si and 2 Ag:Si and Ag:Al 2 O 3 [199,200]. VB diffusivity in h-BN data Ag:Al2 O3 [199,200]. VB diffusivity in h-BN data corresponds to [201]. corresponds to [201]. The diffusivity required to fit the the experimental experimental data data ofof catastrophic catastrophic BD BD in in both both SiO SiO22 high-k (Al and high-k (Al22OO33 or HfO HfO22) stacks stacks with with metal metalelectrodes electrodesare areofofthetheorder orderofof − 13 1010 cm2 /s −13cm 2/s at at 1000 K, with activation energies energies ranging ranging fromfrom 0.3 0.3 to to 0.7 eV [169], where where suchsuch values values areare inin a range range compatible compatiblewith withthe thediffusivity diffusivityofofmetals metalsinindielectrics dielectrics (see (seein in Figure Figure33.33.thethe case of case Cu of Cu diffusion diffusion into SiOSiO into 2 layers 2 layers [198]). [198]). It should be mentioned that the state state transition transition for for aa non-volatile non-volatile regime regime in in metal/h- metal/h- BN/metal BN/metal stacks stacks (i.e., (i.e.,non-reversible non-reversible event) event) hashas been beenstudied studiedwith with a similar approach a similar where approach the best fitting value where the best fitting value E act is also much larger Eact is also much larger (E = 1.3 eV) [184,189]. Such discrepancies act (Eact = 1.3 eV) [184,189]. Such discrep- may ancieslay mayon the lay fact on thethatfactthethatparticular species involved the particular in the electromigration species involved in the electromigrationand/or diffusion process may change, depending on the severity of and/or diffusion process may change, depending on the severity of the SET event (volatile the SET event (volatile and non- volatile), as it occurs and non-volatile), asbetween it occursthe BD and between theSETBDevents and SET in HfO events 2 stacks. in HfO While in the 2 stacks. WhileBD in event the in BDHfO event2 the in diffusing HfO2 the ion speciesion diffusing arespecies considered to be the metallic are considered to be the ions from the metallic ionselectrodes from the −13 2 (D 0 = 1 × 10 electrodes (D0 =cm 10−13 E 1 × /sec, cm = 0.3–0.7 act2/sec, Eact =eV [159,169,189]), 0.3–0.7 the migration eV [159,169,189]), of oxygen of the migration vacancies oxygen from the TMO layer are responsible of the SET transition event (D = 1 × 10 −6 cm2 /s, vacancies from the TMO layer are responsible of the SET transition event (D0 = 1 × 10−6 0 E cm = 0.52 act2/s, Eact =eV) 0.52[188]. eV) [188]. To To make clear the make clear the impact impact of of both both thethe ttoxox reduction reduction and and thethe diffusivity diffusivity increase increase (D ), (D00), alternative fitting values were used to plot curves N ◦ 2, N◦ 3 and N◦ 4 in Figure 31. Curve alternative fitting values were used to plot curves N° 2, N° 3 and N° 4 in Figure 31. Curve N ◦ coincides with the TR for gate-oxide BD, as it considers diffusivity of metals in oxide N°44 coincides with the TR for gate-oxide BD, as it considers diffusivity of metals in oxide layers layers andand no no ttox reduction. The comparison of curves N◦ 1 and N◦ 4 evidence that TR is ox reduction. The comparison of curves N° 1 and N° 4 evidence that TR is significantly significantly higherhigher than than the the TR TR expected expected for for the the voltage voltage rangerange considered. considered. It It is is important important to point out that TR calculated considering only a t ox reduction to point out that TR calculated considering only a tox reduction or an increase in diffusivity or an increase in diffusivity (D0 ) cannot meet the experimental data. This can be interpreted as that the two factors (D0) cannot meet the experimental data. This can be interpreted as that the two factors determine the dependence with the TR voltage, since none of them can separately adjust the results independently. 6. Conclusions A comprehensive and exhaustive revision of RRAM thermal models is presented. Nanomaterials 2021, 11, 1261 38 of 46 determine the dependence with the TR voltage, since none of them can separately adjust the results independently. 6. Conclusions A comprehensive and exhaustive revision of RRAM thermal models is presented. Different approaches have been considered and described, including different conductive filament geometries, operation regimes, filament lateral heat losses, several temperatures to characterize each conductive filament, etc. A 3D numerical solution of the heat equation within a RRAM simulator was used. In addition, analytical models have been developed using equations describing the relevant physics behind the heat equation accounting for steady-state and non-steady-state device operation. A general memristor modeling framework was formulated considering the temperature as a state variable. Moreover, the quantum perspective was included in a mesoscopic context; this is a must due to the nanometric dimensions of the devices under study. In this respect, thermal effects were considered in the formulation of the quantum point contact model. Finally, conductive filament thermometry in usual RRAM technology was studied in detail. Since the physics underlying resistive switching is mainly based on thermally activated physical mechanisms, an accurate description of thermal effects is essential. As far as we know, there are not similar manuscripts that put together these types of effects under one roof. Author Contributions: J.B.R., E.M., R.P., F.P., F.J.-M., S.G.S., J.S., G.G.-C., D.M., E.M. (Enrique Moreno) performed software development and simulations. D.M., G.G.-C., E.M. (Enrique Moreno), E.M. (Enrique Miranda), R.P., M.M.A.C., S.B.B., F.J.-M., J.B.R. wrote the original draft. J.B.R., E.M. (Enrique Miranda), R.P., F.P., F.J.-M., S.G.S., C.d.B., J.S., L.O.C. reviewed and edited the text. All the authors contributed to the English and reference checking. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Spanish Ministry of Science, Innovation and Universities, grant number TEC2017-84321-C4-3-R, TEC2017-84321-C4-4-R and by the Consejería de Economía y Conocimiento de la Junta de Andalucía and European Regional Development Fund (ERDF) under projects A-TIC-117-UGR18. It was also funded by MINCyT of Argentina (Contracts PICT2013/1210, PICT2016/0579 and PME2015-0196), CONICET (Project PIP- 11220130100077CO) and UTN.BA (Projects PIDUTN EIUTIBA4395TC3, CCUTIBA4764TC, MATUNBA4936 and CCUTNBA5182). Data Availability Statement: The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. Acknowledgments: We would like to thank F. Campabadal and M. B. González from the IMB-CNM (CSIC) in Barcelona for fabricating and measuring the devices which were employed to fit some of the compact models whose thermal blocks are presented here. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. References 1. Lanza, M.; Wong, H.-S.P.; Pop, E.; Ielmini, D.; Strukov, D.; Regan, B.C.; Larcher, L.; Villena, M.A.; Yang, J.J.; Goux, L.; et al. Recommended Methods to Study Resistive Switching Devices. Adv. Electron. Mater. 2019, 5, 1800143. [CrossRef] 2. Pan, F.; Gao, S.; Chen, C.; Song, C.; Zeng, F. Recent progress in resistive random access memories: Materials, switching mechanisms, and performance. Mater. Sci. Eng. R Rep. 2014, 83, 1–59. [CrossRef] 3. Villena, M.A.; Roldán, J.B.; Jiménez-Molinos, F.; Miranda, E.; Suñé, J.; Lanza, M. SIM2 RRAM: A physical model for RRAM devices simulation. J. Comput. Electron. 2017, 16, 1095–1120. [CrossRef] 4. IRDS. The International Roadmap for Devices and Systems: More Moore IEEE; IRDS: New York, NY, USA, 2020. 5. Carboni, R.; Ielmini, D. Stochastic memory devices for security and computing. Adv. Electron. Mater. 2019, 5, 1900198. [CrossRef] 6. Puglisi, F.M.; Larcher, L.; Padovani, A.; Pavan, P. A Complete Statistical Investigation of RTN in HfO2 -Based RRAM in High Resistive State. IEEE Trans. Electron Devices 2015, 62, 2606–2613. [CrossRef] 7. Wei, Z.; Katoh, Y.; Ogasahara, S.; Yoshimoto, Y.; Kawai, K.; Ikeda, Y.; Eriguchi, K.; Ohmori, K.; Yoneda, S. True random number generator using current difference based on a fractional stochastic model in 40-nm embedded ReRAM. In Proceedings of the 2016 IEEE International Electron Devices Meeting (IEDM), IEEE, San Francisco, CA, USA, 3–7 December 2016; pp. 4.8.1–4.8.4. Nanomaterials 2021, 11, 1261 39 of 46 8. Puglisi, F.M.; Zagni, N.; Larcher, L.; Pavan, P. random telegraph noise in resistive random access memories: Compact modeling and advanced circuit design. IEEE Trans. Electron Devices 2018, 65, 2964–2972. [CrossRef] 9. Lanza, M.; Wen, C.; Li, X.; Zanotti, T.; Puglisi, F.M.; Shi, Y.; Saiz, F.; Antidormi, A.; Roche, S.; Zheng, W.X.; et al. Advanced data encryption using two-dimensional materials. Adv. Mater. 2021, in press. 10. Yao, P.; Wu, H.; Gao, B.; Tang, J.; Zhang, Q.; Zhang, W.; Yang, J.J.; Qian, H. Fully hardware-implemented memristor convolutional neural network. Nat. Cell Biol. 2020, 577, 641–646. [CrossRef] [PubMed] 11. Merolla, P.A.; Arthur, J.V.; Alvarez-Icaza, R.; Cassidy, A.S.; Sawada, J.; Akopyan, F.; Jackson, B.L.; Imam, N.; Guo, C.; Nakamura, Y.; et al. A million spiking-neuron integrated circuit with a scalable communication network and interface. Sci- ence 2014, 345, 668–673. [CrossRef] 12. Yu, S.; Gao, B.; Fang, Z.; Yu, H.; Kang, J.; Wong, H.-S.P. A neuromorphic visual system using RRAM synaptic devices with Sub-pJ energy and tolerance to variability: Experimental characterization and large-scale modeling. In Proceedings of the 2012 International Electron Devices Meeting, San Francisco, CA, USA, 10–13 December 2012; pp. 10.4.1–10.4.4. 13. Zidan, M.A.; Strachan, J.P.; Lu, W.D. The future of electronics based on memristive systems. Nat. Electron. 2018, 1, 22–29. [CrossRef] 14. Prezioso, M.; Merrikh-Bayat, F.; Hoskins, B.D.; Adam, G.C.; Likharev, K.K.; Strukov, D.B. Training and operation of an integrated neuromorphic network based on metal-oxide memristors. Nature 2015, 521, 61–64. [CrossRef] [PubMed] 15. Romero-Zaliz, R.; Pérez, E.; Jiménez-Molinos, F.; Wenger, C.; Roldán, J. Study of Quantized Hardware Deep Neural Networks Based on Resistive Switching Devices, Conventional versus Convolutional Approaches. Electronics 2021, 10, 346. [CrossRef] 16. Quesada, E.P.-B.; Romero-Zaliz, R.; Pérez, E.; Mahadevaiah, M.K.; Reuben, J.; Schubert, M.; Jiménez-Molinos, F.; Roldán, J.; Wenger, C. Toward Reliable Compact Modeling of Multilevel 1T-1R RRAM Devices for Neuromorphic Systems. Electronics 2021, 10, 645. [CrossRef] 17. Mead, C.; Ismail, M. Analog VLSI Implementation of Neural Systems; Springer: Berlin/Heidelberg, Germany, 1989. 18. Ielmini, D. Resistive switching memories based on metal oxides: Mechanisms, reliability and scaling. Semicond. Sci. Technol. 2016, 31, 063002. [CrossRef] 19. Ielmini, D.; Milo, V. Physics-based modeling approaches of resistive switching devices for memory and in-memory computing applications. J. Comput. Electron. 2017, 16, 1121–1143. [CrossRef] 20. Huang, P.; Gao, B.; Chen, B.; Zhang, F.; Liu, L.; Du, G. Stochastic simulation of forming, SET and RESET process for transition metal oxide-based resistive switching memory. Proc. SISPAD 2011, 2012, 312–315. 21. Aldana, S.; Roldán, J.B.; García-Fernández, P.; Suñe, J.; Romero-Zaliz, R.; Jiménez-Molinos, F.; Long, S.; Gómez-Campos, F.; Liu, M. An in-depth description of bipolar resistive switching in Cu/HfOx /Pt devices, a 3D Kinetic Monte Carlo simulation approach. J. Appl. Phys. 2018, 123, 154501. [CrossRef] 22. Garcia-Redondo, F.; Gowers, R.P.; Crespo-Yepes, A.; Lopez-Vallejo, M.; Jiang, L. SPICE Compact modeling of bipolar/unipolar memristor switching governed by electrical thresholds. IEEE Trans. Circuits Syst. I Regul. Pap. 2016, 63, 1255–1264. [CrossRef] 23. Dirkmann, S.; Kaiser, J.; Wenger, C.; Mussenbrock, T. Filament growth and resistive switching in hafnium oxide memristive devices. ACS Appl. Mater. Interfaces 2018, 10, 14857–14868. [CrossRef] 24. Aldana, S.; García-Fernández, P.; Rodríguez-Fernández, A.; Romero-Zaliz, R.; González, M.B.; Jiménez-Molinos, F.; Campabadal, F.; Gómez-Campos, F.; Roldán, J.B. A 3D kinetic monte carlo simulationstudy of resistive switching processes in Ni/HfO2 /Si-n+-based RRAMs. J. Phys. D 2017, 50, 335103. [CrossRef] 25. Jagath, A.L.; Nandha Kumar, T.; Almurib, H.A.F. Modeling of Current Conduction during RESET Phase of Pt/Ta2 O5 /TaOx /Pt Bipolar Resistive RAM Devices. In Proceedings of the 2018 IEEE 7th Non-Volatile Memory Systems and Applications Symposium (NVMSA), Hakodate, Japan, 28–31 August 2018; pp. 55–60. 26. Fang, X.; Yang, X.; Wu, J.; Yi, X. A Compact SPICE model of unipolar memristive devices. IEEE Trans. Nanotechnol. 2013, 12, 843–850. [CrossRef] 27. González-Cordero, G.; González, M.B.; García, H.; Campabadal, F.; Dueñas, S.; Castán, H.; Jiménez-Molinos, F.; Roldán, J.B. A physically based model for resistive memories including a detailed temperature and variability description. Microelectron. Eng. 2017, 178, 26–29. [CrossRef] 28. Karpov, V.; Niraula, D.; Karpov, I. Thermodynamic analysis of conductive filaments. Appl. Phys. Lett. 2016, 109, 093501. [CrossRef] 29. Maestro-Izquierdo, M.; Gonzalez, M.B.; JimenezMolinos, F.; Moreno, E.; Roldan, J.B.; Campabadal, F. Unipolar resistive switching behavior in Al2 O3 /HfO2 multilayer dielectric stacks: Fabrication, characterization and simulation. Nanotechnology 2020, 31, 135202. [CrossRef] [PubMed] 30. Larentis, S.; Nardi, F.; Balatti, S.; Gilmer, D.C.; Ielmini, D. Resistive Switching by Voltage-Driven Ion Migration in Bipolar RRAM—Part II: Modeling. IEEE Trans. Electron Devices 2012, 59, 2468–2475. [CrossRef] 31. Jimenez-Molinos, F.; Villena, M.A.; Roldan, J.B.; Roldan, A.M. A SPICE compact model for unipolar RRAM reset process analysis. IEEE Trans. Electron Devices 2015, 62, 955–962. [CrossRef] 32. Menzel, S.; Kaupmann, P.; Waser, R. Understanding filamentary growth in electrochemical metallization memory cells using kinetic Monte Carlo simulations. Nanoscale 2015, 7, 12673–12681. [CrossRef] [PubMed] 33. Picos, R.; Roldán, J.B.; al Chawa, M.M.; García-Fernández, P.; García-Moreno, F.J.Y.E. Semiempirical modeling of reset transitions in unipolar resistive-switching based memristors. Radioeng. J. 2015, 24, 420–424. [CrossRef] Nanomaterials 2021, 11, 1261 40 of 46 34. Vandelli, L.; Padovani, A.; Larcher, L.; Bersuker, G. Microscopic modeling of electrical stress-induced breakdown in poly- crystalline hafnium oxide dielectrics. IEEE Trans. Electron Devices 2013, 60, 1754–1762. [CrossRef] 35. Blasco, J.; Ghenzi, N.; Suñé, J.; Levy, P.; Miranda, E. Equivalent circuit modeling of the bistable conduction characteristics in electroformed thin dielectric films. Microelectron. Reliab. 2015, 55, 1–14. [CrossRef] 36. Maldonado, D.; Gonzalez, M.B.; Campabadal, F.; JimenezMolinos, F.; Al Chawa, M.M.; Stavrinides, S.G.; Roldan, J.B.; Tetzlaff, R.; Picos, R.; Chua, L.O. Experimental evaluation of the dynamic route map in the reset transition of memristive ReRAMs. Chaos Solitons Fractals 2020, 139, 110288. [CrossRef] 37. Jiménez-Molinos, F.; González-Cordero, G.; Cartujo-Cassinello, P.; Roldán, J.B. SPICE modeling of RRAM thermal reset transition for circuit simulation purposes. In Proceedings of the Spanish Conference on Electron Devices, Barcelona, Spain, 8–10 February 2017. 38. Bocquet, M.; Deleruyelle, D.; Aziza, H.; Muller, C.; Portal, J.-M.; Cabout, T.; Jalaguier, E. Robust compact model for bipolar oxide-based resistive switching memories. IEEE Trans. Electron. Devices 2014, 61, 674–681. [CrossRef] 39. al Chawa, M.M.; Picos, R.; Tetzlaff, R. A Simple Memristor Model for Neuromorphic ReRAM Devices. In Proceedings of the 2020 IEEE International Symposium on Circuits and Systems (ISCAS), Seville, Spain, 10–21 October 2020. 40. al Chawa, M.M.; Picos, R. A simple quasi-static compact model of bipolar ReRAM memristive devices. IEEE Trans. Circuits Syst. II 2020, 67, 390–394. [CrossRef] 41. Panda, D.; Sahu, P.P.; Tseng, T.Y. A collective study on modeling and simulation of resistive random access memory. Nanoscale Res. Lett. 2018, 13, 1–48. [CrossRef] [PubMed] 42. Reuben, J.; Biglari, M.; Fey, D. Incorporating Variability of Resistive RAM in Circuit Simulations Using the Stanford–PKU Model. IEEE Trans. Nanotechnol. 2020, 19, 508–518. [CrossRef] 43. Mikhaylov, A.; Guseinov, D.; Belov, A.; Korolev, D.; Shishmakova, V.; Koryazhkina, M.; Filatov, D.; Gorshkov, O.; Maldonado, D.; Alonso, F.; et al. Stochastic resonance in a metal-oxide memristive device. Chaos Solitons Fractals 2021, 144, 110723. [CrossRef] 44. Aldana, S.; Pérez, E.; JimenezMolinos, F.; Wenger, C.; Roldán, J.B. Kinetic Monte Carlo analysis of data retention in Al:HfO2 -based resistive random access memories. Semicond. Sci. Technol. 2020, 35, 115012. [CrossRef] 45. Roldán, J.B.; Alonso, F.J.; Aguilera, A.M.; Maldonado, D.; Lanza, M. Time series statistical analysis: A powerful tool to evaluate the variability of resistive switching memories. J. Appl. Phys. 2019, 125, 174504. [CrossRef] 46. Miranda, E.; Mehonic, A.; Ng, W.H.; Kenyon, A.J. Simulation of cycle-to-cycle instabilities in SiOx -based ReRAM devices using a self-correlated process with long-term variation. IEEE EDL 2019, 40, 28–31. [CrossRef] 47. Kvatinsky, S.; Ramadan, M.; Friedman, E.G.; Kolodny, A. VTEAM: A general model for voltage-controlled memristors. IEEE Trans. Circuits Syst. II 2015, 62, 786–790. [CrossRef] 48. Picos, R.; Roldan, J.B.; Al Chawa, M.M.; JimenezMolinos, F.; Garcia-Moreno, E. A physically based circuit model to account for variability in memristors with resistive switching operation. In Proceedings of the 2016 Conference on Design of Circuits and Integrated Systems (DCIS), Granada, Spain, 23–25 November 2016; pp. 1–6. 49. al Chawa, M.M.; de Benito, C.; Picos, R. A simple piecewise model of reset/set transitions in bipolar ReRAM memristive devices. IEEE Trans. Circuits Syst. I 2018, 65, 3469–3480. [CrossRef] 50. al Chawa, M.M.; Tetzlaff, R.; Picos, R. A Simple Monte Carlo Model for the Cycle-to-Cycle Reset Transition Variation of ReRAM Memristive Devices. In Proceedings of the 9th International Conference on Modern Circuits and Systems Technologies (MOCAST), Bremen, Germany, 7–9 September 2020. 51. Alonso, F.J.; Maldonado, D.; Aguilera, A.M.; Roldan, J.B. Memristor variability and stochastic physical properties modeling from a multivariate time series approach. Chaos Solitons Fractals 2021, 143, 110461. [CrossRef] 52. Pérez, E.; Maldonado, D.; Acal, C.; Ruiz-Castro, J.E.; Alonso, F.J.; Aguilera, A.M.; Jiménez-Molinos, F.; Wenger, C.; Roldán, J.B. Analysis of the statistics of device-to-device and cycle-to-cycle variability in TiN/Ti/Al:HfO2 /TiN RRAMs. Microelectron. Eng. 2019, 214, 104–109. [CrossRef] 53. Aldana, S.; García-Fernández, P.; Romero-Zaliz, R.; González, M.B.; Jiménez-Molinos, F.; Gómez-Campos, F.; Campabadal, F.; Roldán, J.B. Resistive switching in HfO2 based valence change memories, a comprehensive 3D kinetic Monte Carlo approach. J. Phys. D 2020, 53, 225106. [CrossRef] 54. Guy, J.; Molas, G.; Blaise, P.; Bernard, M.; Roule, A.; Le Carval, G.; Delaye, V.; Toffoli, A.; Ghibaudo, G.; Clermidy, F.; et al. Investigation of forming, SET, and data retention of conductive-bridge random-access memory for stack optimization. IEEE Trans. Electron Devices 2015, 62, 3482–3489. [CrossRef] 55. Villena, M.A.; Roldán, J.B.; González, M.B.; González-Rodelas, P.; Jiménez-Molinos, F.; Campabadal, F.; Barrera, D. A new parameter to characterize the charge transport regime in Ni/HfO2 /Si-n+ -based RRAMs. Solid State Electron. 2016, 118, 56–60. [CrossRef] 56. González-Cordero, G.; Roldán, J.B.; Jiménez-Molinos, F. SPICE simulation of RRAM circuits. A compact modeling perspective. In Proceedings of the 2017 Spanish Conference on Electron Devices, Barcelona, Spain, 8–10 February 2017; pp. 26–29. 57. Huang, P.; Zhu, D.; Chen, S.; Zhou, Z.; Chen, Z.; Gao, B.; Liu, L.; Liu, X.; Kang, J. Compact model of HfOX -based electronic synaptic devices for neuromorphic computing. IEEE Trans. Electron. Devices 2017, 64, 614–621. [CrossRef] 58. Kwon, S.; Jang, S.; Choi, J.-W.; Choi, S.; Jang, S.-J.; Kim, T.-W.; Wang, G. Controllable switching filaments prepared via tunable and well-defined single truncated conical nanopore structures for fast and scalable SiOx memory. Nanoletters 2017, 17, 7462–7470. [CrossRef] Nanomaterials 2021, 11, 1261 41 of 46 59. Villena, M.; Gonzalez, M.B.; Roldán, J.; Campabadal, F.; Jiménez-Molinos, F.; Gómez-Campos, F.; Suñé, J. An in-depth study of thermal effects in reset transitions in HfO2 based RRAMs. Solid-State Electron. 2015, 111, 47–51. [CrossRef] 60. Lohn, A.J.; Mickel, P.R.; Marinella, M.J. Analytical estimations for thermal crosstalk, retention, and scaling limits in filamentary resistive memory. J. Appl. Phys. 2014, 115, 234507. [CrossRef] 61. Sun, P.; Lu, N.; Li, L.; Li, Y.; Wang, H.; Lv, H.; Liu, Q.; Long, S.; Liu, S.; Liu, M. Thermal crosstalk in 3-dimensional RRAM crossbar array. Sci. Rep. 2015, 5, 13504. [CrossRef] [PubMed] 62. Deshmukh, S.; Islam, R.; Chen, C.; Yalon, E.; Saraswat, K.C.; Pop, E. Thermal modeling of metal oxides for highly scaled nanoscale RRAM. In Proceedings of the 2015 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD), Washington, DC, USA, 9–11 September 2015; Volume 2015, pp. 281–284. 63. Wang, D.-W.; Chen, W.; Zhao, W.-S.; Zhu, G.-D.; Kang, K.; Gao, P.; Schutt-Aine, J.E.; Yin, W.-Y. Fully Coupled Electrothermal Simulation of Large RRAM Arrays in the “Thermal-House”. IEEE Access 2018, 7, 3897–3908. [CrossRef] 64. Rodríguez, N.; Roldán, F.G.y.J.B. Modeling of inversion layer centroid and polysilicon depletion effects on ultrathin-gate-oxide MOSFET behaviour: The influence of crystallographic orientation. IEEE Trans. Electron. Devices 2007, 54, 723–732. [CrossRef] 65. González, B.; Roldán, J.; Iniguez, B.; Lazaro, A.; Cerdeira, A. DC self-heating effects modelling in SOI and bulk FinFETs. Microelectron. J. 2015, 46, 320–326. [CrossRef] 66. Roldán, J.B.; Gámiz, F.; JimenezMolinos, F.; Sampedro, C.; Godoy, A.; Rodríguez, N. An analytic I-V model for surrounding-gate MOSFET including quantum and velocity overshoot effects. IEEE Trans. Electron. Devices 2010, 57, 2925–2933. [CrossRef] 67. Blanco-Filgueira, B.; Roldán, P.L.Y.J.B. Analytical modeling of size effects on the lateral photoresponse of CMOS photodiodes. Solid State Electron. 2012, 73, 15–20. [CrossRef] 68. Blanco-Filgueira, B.; Roldán, P.L.y.J.B. A closed-form and explicit analytical model for crosstalk in CMOS photodiodes. IEEE Trans. Electron Devices 2013, 60, 3459–3464. [CrossRef] 69. Gámiz, F.; Godoy, A.; Donetti, L.; Sampedro, C.; Roldán, J.B.; Ruiz, F.; Tienda, I.; Jiménez-Molinos, N.R.Y.F. Monte Carlo simulation of nanoelectronic devices. J. Comput. Electron. 2009, 8, 174–191. [CrossRef] 70. Ielmini, D.; Waser, R. Resistive Switching: From Fundamentals of Nanoionic Redox Processes to Memristive Device Applications; Wiley-VCH: Hoboken, NJ, USA, 2015. 71. Corinto, F.; Civalleri, P.P.; Chua, L.O. A theoretical approach to memristor devices. IEEE J. Emerg. Sel. Top. Circuits Syst. 2015, 5, 123–132. [CrossRef] 72. Chua, L.O. Everything you wish to know about memristors but are afraid to ask. Radioengineering 2015, 24, 319–368. [CrossRef] 73. James, A.P. A hybrid memristor–CMOS chip for AI. Nat. Electron. 2019, 2, 268–269. [CrossRef] 74. Volos, C.K.; Kyprianidis, I.M.; Stavrinides, S.G.; Stouboulos, I.N.; Anagnostopoulos, A.N. Memristors: A new approach in nonlinear circuits design. In Proceedings of the 14th WSEAS International Conference on Communication, Cape Town, South Africa, 23–27 May 2010; pp. 25–30. 75. Li, Y.; Wang, Z.; Midya, R.; Xia, Q.; Yang, J.J. Review of memristor devices in neuromorphic computing: Materials sciences and device challenges. J. Phys. D Appl. Phys. 2018, 51, 503002. [CrossRef] 76. Padovani, A.; Larcher, L.; Pirrotta, O.; Vandelli, L.; Bersuker, G. Microscopic Modeling of HfOx RRAM operations: From forming to switching. IEEE Trans. Electron. Devices 2015, 62, 1998–2006. [CrossRef] 77. Cazorla, M.; Aldana, S.; Maestro, M.; González, M.B.; Campabadal, F.; Moreno, E.; Jiménez-Molinos, F.; Roldán, J.B. A thermal study of multilayer RRAMs based on HfO2 and Al2 O3 oxides. J. Vac. Sci. Technol. B 2019, 37, 012204. [CrossRef] 78. Guan, X.; Yu, S.; Wong, H.-S.P. A SPICE compact model of metal oxide resistive switching memory with variations. IEEE Electron Device Lett. 2012, 33, 1405–1407. [CrossRef] 79. González-Cordero, G.; Roldan, J.B.; Jiménez-Molinos, F.; Suñé, J.; Liu, S.L.y.M. A new model for bipolar RRAMs based on truncated cone conductive filaments, a Verilog-A approach. Semicond. Sci. Technol. 2016, 31, 115013. [CrossRef] 80. Villena, M.A.; González, M.B.; Jiménez-Molinos, F.; Campabadal, F.; Roldán, J.B.; Suñé, J.; Romera, E.; Miranda, E. Simulation of thermal reset transitions in resistive switching memories including quantum effects. J. Appl. Phys. 2014, 115, 214504. [CrossRef] 81. Von Witzleben, M.; Fleck, K.; Funck, C.; Baumkötter, B.; Zuric, M.; Idt, A.; Breuer, T.; Waser, R.; Böttger, U.; Menzel, S. Investigation of the impact of high temperatures on the switching kinetics of redox-based resistive switching cells using a high-speed nanoheater. Adv. Electron. Mater. 2017, 3, 1700294. [CrossRef] 82. Lantos, N.; Nataf, F. Perfectly matched layers for the heat and advection–diffusion equations. J. Comput. Phys. 2010, 229, 9042–9052. [CrossRef] 83. Moreno, E.; Hemmat, Z.; Roldan, J.B.; Pantoja, M.F.; Bretones, A.R.; Garcia, S.G.; Faez, R. Implementation of open boundary problems in photo-conductive antennas by using convolutional perfectly matched layers. IEEE Trans. Antennas Propag. 2016, 64, 4919–4922. [CrossRef] 84. González-Cordero, G. Compact Modeling of Memristors Based on Resistive Switching Devices. Ph.D. Thesis, Universidad de Granada, Granada, Spain, 2020. 85. Villena, M.A.; Jiménez-Molinos, F.; Roldan, J.B.; Suñe, J.; Long, S.; Lian, X.; Gamiz, F.; Liu, M. An in-depth simulation study of thermal reset transitions in resistive switching memories. J. Appl. Phys. 2013, 114, 144505. [CrossRef] 86. Guan, X.; Yu, S.; Wong, H.S.P. On the Variability of HfOx RRAM: From Numerical Simulation to Compact Modeling Technical. In Proceedings of the 2012 NSTI Nanotechnology Conference and Expo, NSTI-Nanotech, Santa Clara, CA, USA, 18–21 June 2012; Volume 2, pp. 815–820. Nanomaterials 2021, 11, 1261 42 of 46 87. Jiang, Z.; Yu, S.; Wu, Y.; Engel, J.H.; Guan, X.; Wong, H.-S.P. Verilog-A compact model for oxide-based resistive random access memory (RRAM). In Proceedings of the 2014 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD), Yokohama, Japan, 9–11 September 2014; pp. 41–44. 88. Jiang, Z.; Wu, Y.; Yu, S.; Yang, L.; Song, K.; Karim, Z.; Wong, H.-S.P. A compact model for metal–oxide resistive random access memory with experiment verification. IEEE Trans. Electron Devices 2016, 63, 1884–1892. [CrossRef] 89. Moran, M.J.; Shapiro, H.N.; Munson, B.R.; Dewitt, D.P.; Wiley, J.; Hepburn, K.; Fleming, L. Introduction to Thermal Systems Engineering: And Heat Transfer; John Wiley & Sons Inc.: New York, NY, USA, 2003; Volume 169, ISBN 0471204900. 90. Sheridan, P.; Kim, K.-H.; Gaba, S.; Chang, T.; Chen, L.; Lu, W. Device and SPICE modeling of RRAM devices. Nanoscale 2011, 3, 3833–3840. [CrossRef] [PubMed] 91. Huang, P.; Liu, X.Y.; Chen, B.; Li, H.T.; Wang, Y.J.; Deng, Y.X.; Wei, K.L.; Zeng, L.; Gao, B.; Du, G.; et al. A physics-based compact model of metal-oxide-based RRAM DC and AC operations. IEEE Trans. Electron Devices 2013, 60, 4090–4097. [CrossRef] 92. Li, H.; Huang, P.; Gao, B.; Chen, B.; Liu, X.; Kang, J. A SPICE Model of Resistive Random Access Memory for Large-Scale Memory Array Simulation. IEEE Electron Device Lett. 2013, 35, 211–213. [CrossRef] 93. Li, H.; Jiang, Z.; Huang, P.; Wu, Y.; Chen, H.Y.; Gao, B.; Liu, X.Y.; Kang, J.F.; Wong, H.S. Variation-Aware, Reliability-Emphasized Design and Optimization of RRAM using SPICE Model. In Proceedings of the Design, Automation & Test in Europe Conference & Exhibition, Grenoble, France, 9–13 March 2015; pp. 1425–1430. 94. Chen, A. A review of emerging non-volatile memory (NVM) technologies and applications. Solid State Electron. 2016, 125, 25–38. [CrossRef] 95. Chen, P.-Y.; Yu, S. Compact modeling of RRAM devices and its applications in 1T1R and 1S1R array design. IEEE Trans. Electron. Devices 2015, 62, 4022–4028. [CrossRef] 96. Chiang, M.-H.; Hsu, K.-H.; Ding, W.-W.; Yang, B.-R. A predictive compact model of bipolar RRAM cells for circuit simulations. IEEE Trans. Electron. Devices 2015, 62, 2176–2183. [CrossRef] 97. Kwon, J.; Sharma, A.A.; Chen, C.M.; Fantini, A.; Jurczak, M.; Herzing, A.A.; Bain, J.A.; Picard, Y.N.; Skowronski, M. Transient thermometry and high-resolution transmission electron microscopy analysis of filamentary resistive switches. ACS Appl. Mater. Interfaces 2016, 8, 20176. [CrossRef] 98. Sharma, A.A.; Noman, M.; Skowronski, M.; Bain, J.A. Technology, Systems and Applications (VLSI-TSA). In Proceedings of the 2014 International Symposium on VLSI Technology, Systems and Applications, Hsinchu, Taiwan, 28–30 April 2014; p. 1. 99. Nishi, Y.; Menzel, S.; Fleck, K.; Boettger, U.; Waser, R. Origin of the SET kinetics of the resistive switching in tantalum oxide thin films. IEEE Electron. Device Lett. 2013, 35, 259–2061. [CrossRef] 100. Kim, S.; Du, C.; Sheridan, P.; Ma, W.; Choi, S.; Lu, W.D. Experimental demonstration of a second-order memristor and its ability to biorealistically implement synaptic plasticity. Nano Lett. 2015, 15, 2203–2211. [CrossRef] 101. Panzer, M.A.; Shandalov, M.; Rowlette, J.A.; Oshima, Y.; Chen, Y.W.; McIntyre, P.C.; Goodson, K.E. Thermal Properties of Ultrathin Hafnium Oxide Gate Dielectric Films. IEEE Electron Device Lett. 2009, 30, 1269–1271. [CrossRef] 102. Scott, E.A.; Gaskins, J.T.; King, S.W.; Hopkins, P.E. Thermal conductivity and thermal boundary resistance of atomic layer deposited high-k dielectric aluminum oxide, hafnium oxide, and titanium oxide thin films on silicon. APL Mater. 2018, 6, 058302. [CrossRef] 103. Roldán, A.M.; Roldán, J.B.; Reig, C.; Cubells-Beltrán, M.-D.; Ramírez, D.; Cardoso, S.; Freitas, P.P. A DC behavioral electrical model for quasi-linear spin-valve devices including thermal effects for circuit simulation. Microelectron. J. 2011, 42, 365–370. [CrossRef] 104. Busani, M.; Menozzi, R.; Borgarino, M.; Fantini, F. Dynamic thermal characterization and modeling of packaged AlGaAs/GaAs HBTs. IEEE Trans. Compon. Packag. Technol. 2000, 23, 352–359. [CrossRef] 105. Pedro, M.; Martin-Martinez, J.; Gonzalez, M.; Rodriguez, R.; Campabadal, F.; Nafria, M.; Aymerich, X. Tuning the conductivity of resistive switching devices for electronic synapses. Microelectron. Eng. 2017, 178, 89–92. [CrossRef] 106. González-Cordero, G.; Pedro, M.; Martin-Martinez, J.; González, M.; Jiménez-Molinos, F.; Campabadal, F.; Nafría, N.; Roldán, J. Analysis of resistive switching processes in TiN/Ti/HfO2 /W devices to mimic electronic synapses in neuromorphic circuits. Solid-State Electron. 2019, 157, 25–33. [CrossRef] 107. Kumar, S.; Williams, R.S.; Wang, Z. Third-order nanocircuit elements for neuromorphic engineering. Nat. Cell Biol. 2020, 585, 518–523. [CrossRef] 108. González-Cordero, G.; Jiménez-Molinos, F.; Roldán, J.B.; González, M.B.; Campabadal, F. Transient SPICE Simulation of Ni/HfO2 /Si-n+ Resistive Memories. In Proceedings of the Design of Circuits and Integrated Systems Conference, DCIS, Granada, Spain, 23–25 November 2016. 109. González-Cordero, G.; Jiménez-Molinos, F.; Villena, M.A.; Roldán, J.B. SPICE Simulation of Thermal Reset Transitions in Ni/HfO2 /Si-n+ RRAMs Including Quantum Effects. In Proceedings of the 19th Workshop on Dielectrics in Microelectronics, WoDiM, Catania, Italy, 27–30 June 2016. 110. González, M.B.; Rafí, J.M.; Beldarrain, O.; Zabala, M.; Campabadal, F. Analysis of the switching variability in Ni/HfO2 -based RRAM devices. IEEE Trans. Device Mater. Reliab. 2014, 14, 769–771. 111. Wang, W.; Laudato, M.; Ambrosi, E.; Bricalli, A.; Covi, E.; Lin, Y.-H.; Ielmini, D. Volatile Resistive Switching Memory Based on Ag Ion Drift/Diffusion—Part II: Compact Modeling. IEEE Trans. Electron Devices 2019, 66, 3802–3808. [CrossRef] 112. Chua Leon, O.; Kang, S.M. Memristive devices and systems. Proc. IEEE 1976, 64, 209–223. [CrossRef] Nanomaterials 2021, 11, 1261 43 of 46 113. Ginoux, J.M.; Muthuswamy, B.; Meucci, R.; Euzzor, S.; Di Garbo, A.; Ganesan, K. A physical memristor based Muthuswamy- Chua-Ginoux system. Sci. Rep. 2020, 10, 1–10. [CrossRef] 114. Steinhart, J.S.; Hart, S.R. Calibration curves for thermistors. Deep Sea Res. Oceanogr. Abstr. 1968, 15, 497–503. [CrossRef] 115. Theodorakakos, A.; Stavrinides, S.G.; Hatzikraniotis, E.; Picos, R. A non-ideal memristor device. In Proceedings of the 2015 International Conference on Memristive Systems (MEMRISYS), Paphos, Cyprus, 8–10 November 2015; pp. 1–2. 116. Biolek, D.; Biolek, Z.; Biolkova, V.; Kolka, Z. Some fingerprints of ideal memristors. In Proceedings of the 2013 IEEE International Symposium on Circuits and Systems (ISCAS), Beijing, China, 19–23 May 2013; pp. 201–204. 117. Wagner, G.; Jones, R.; Templeton, J.; Parks, M. An atomistic-to-continuum coupling method for heat transfer in solids. Comput. Methods Appl. Mech. Eng. 2008, 197, 3351–3365. [CrossRef] 118. Xu, Z. Heat transport in low-dimensional materials: A review and perspective. Theor. Appl. Mech. Lett. 2016, 6, 113–121. [CrossRef] 119. Mosso, N.; Drechsler, U.; Menges, F.; Nirmalraj, P.; Karg, S.; Riel, H.; Gotsmann, B. Heat transport through atomic contacts. Nat. Nanotech. 2017, 12, 430–433. [CrossRef] [PubMed] 120. Hanggi, P.; Talkner, P.; Borkovec, H. Reaction-rate theory: Fifty years after Kramers. Rev. Mod. Phys. 1990, 62, 251. [CrossRef] 121. Chiu, F. A review on conduction mechsnisms in dielectric films. Adv. Mater. Sci. Eng. 2014, 2014, 578168. [CrossRef] 122. Waser, R.; Dittmann, R.; Saikov, G.; Szot, K. Redox-based resistive switching memories nanoionic mechanisms, prospects, and challenges. Adv. Mater. 2009, 21, 2632–2663. [CrossRef] 123. Datta, S. Electronic Transport in Mesoscopic Systems; Cambridge University Press: Cambridge, UK, 1995. 124. Kouwenhoven, L.P.; Van Wees, B.J.; Harmans, C.J.P.M.; Williamson, J.G.; Van Houten, H.; Beenakker, C.W.J.; Foxon, C.T.; Harris, J.J. Nonlinear conductance of quantum point contacts. Phys. Rev. B 1989, 39, 8040–8043. [CrossRef] [PubMed] 125. van Ruitenbeek, J.; Masis, M.M.; Miranda, E. Quantum point contact conduction. In Resistice Switching: From Fundamentals of Nanoinic Redox Processes to Memristive Device Applications; Ielmini, D., Waser, R., Eds.; John Wiley & Sons Inc.: New York, NY, USA, 2016; pp. 197–224. 126. Agrait, N.; Yeyati, A.L.; van Ruitenbeek, J.M. Quantum properties of atomic-sized conductors. Phys. Rep. 2003, 377, 81. [CrossRef] 127. Li, Y.; Long, S.; Liu, Y.; Hu, C.; Teng, J.; Liu, Q.; Lv, H.; Suñé, J.; Liu, M. Conductance quantization in resistive random access memory. Nanoscale Res. Lett. 2015, 10, 420. [CrossRef] [PubMed] 128. Suñé, J.; Miranda, E.; Nafría, M.; Aymerich, X. Point contact conduction at the oxide breakdown of MOS devices. In Proceedings of the IEEE International Electron Device Meeting (IEDM), San Francisco, CA, USA, 6–9 December 1998; p. 191. 129. Suñé, J.; Miranda, E.; Nafría, M.; Aymerich, X. Modeling the breakdown spots in silicon dioxide films as point contacts. Appl. Phys. Lett. 1999, 75, 959–961. [CrossRef] 130. Mehonic, A.; Vrajitoarea, A.; Cueff, S.; Hudziak, S.; Howe, H.; Labbe, C.; Rizk, R.; Pepper, M.; Kenyon, A.J. Quantum conductance in silicon oxide resistive memory devices. Sci. Rep. 2013, 3, 2708. [CrossRef] [PubMed] 131. Nandakumar, S.R.; Minvielle, M.; Nagar, S.; Dubourdieu, C.; Rajendran, B. A 250 mV Cu/SiO2 /W Memristor with Half-Integer Quantum Conductance States. Nano Lett. 2016, 16, 1602–1608. [CrossRef] [PubMed] 132. Miranda, E.; Walczyk, C.; Wenger, C.; Schroeder, T. Model for the resistive switching effect in HfO2 MIM structures based on the transmission properties of narrow constrictions. IEEE Electron. Device Lett. 2010, 31, 609–611. [CrossRef] 133. Degraeve, R.; Roussel, P.; Goux, L.; Wouters, D.; Kittl, J.; Altimime, L.; Jurczak, M.; Groeseneken, G. Generic learning of TDDB applied to RRAM for improved understanding of conduction and switching mechanism through multiple filaments. In Proceedings of the 2010 International Electron Devices Meeting, San Francisco, CA, USA, 6–8 December 2010. 134. Walczyk, C.; Walczyk, D.; Schroeder, T.; Bertaud, T.; Sowinska, M.; Lukosius, M.; Fraschke, M.; Wolansky, D.; Tillack, B.; Miranda, E.; et al. Impact of temperature on the resistive switching behavior of embedded HfO2 -Based RRAM devices. IEEE Trans. Electron. Dev. 2011, 58, 3124–3131. [CrossRef] 135. Long, S.; Lian, X.; Cagli, C.; Cartoixà, X.; Rurali, R.; Miranda, E.; Jiménez, D.; Perniola, L.; Liu, M.; Suñé, J. Quantum-size effects in hafnium-oxide resistive switching. Appl. Phys. Lett. 2013, 102, 183505. [CrossRef] 136. Prócel, L.M.; Trojman, L.; Moreno, J.; Crupi, F.; Maccaronio, V.; Degraeve, R.; Goux, L.; Simoen, E. Experimental evidence of the quantum point contact theory in the conduction mechanism of bipolar HfO2 -based resistive random access memories. J. Appl. Phys. 2013, 114, 074509. [CrossRef] 137. Rahavan, N. Performance and reliability trade-offs for high-K RRAM. Microelectron. Reliab. 2014, 54, 2253–2257. [CrossRef] 138. Roldán, J.B.; Miranda, E.; González-Cordero, G.; García-Fernández, P.; Romero-Zaliz, R.; González-Rodelas, P.; Aguilera, A.M.; González, M.B.; Jiménez-Molinos, F. Multivariate analysis and extraction of parameters in resistive RAMs using the Quantum Point Contact model. J. Appl. Phys. 2018, 123, 014501. [CrossRef] 139. Calixto, D.; Maldonado, E.; Miranda, J.B.; Roldán, M. Modeling of the temperature effects in filamentary-type resistive switching memories using quantum point-contact theory. J. Phys. D Appl. Phys. 2020, 53, 295106. [CrossRef] 140. Tsuruoka, T.; Hasegawa, T.; Terabe, K.; Aono, M. Conductance quantization and synaptic behavior in a Ta2 O5 -based atomic switch. Nanotechnology 2012, 23, 435705. [CrossRef] [PubMed] 141. Chen, C.; Gao, S.; Zeng, F.; Wang, G.; Li, S.; Song, C.; Pan, F. Conductance quantization in oxygen-anion-migration-based resistive switching memory devices. Appl. Phys. Lett. 2013, 103, 043510. [CrossRef] 142. Yi, W.; Savelev, S.; Medeiros-Ribeiro, G.; Miao, F.; Zhang, M.; Yang, J.; Bratkovsky, A.; Williams, R.S. Quantized conductance coincides with state instability and excess noise in tantalum oxide memristors. Nat. Commun. 2016, 7, 11142. [CrossRef] Nanomaterials 2021, 11, 1261 44 of 46 143. Ye, J.Y.; Li, Y.Q.; Gao, J.; Peng, H.Y.; Wu, S.X.; Wu, T. Nanoscale resistive switching and filamentary conduction in NiO thin films. Appl. Phys. Lett. 2010, 97, 132108. [CrossRef] 144. Nishi, Y.; Sasakura, H.; Kimoto, T. Appearance of quantum point contact in Pt/NiO/Pt resistive switching cells. J. Mater. Res. 2017, 32, 2631–2637. [CrossRef] 145. Zhu, X.-J.; Shang, J.; Li, R.-W. Resistive switching effects in oxide sandwiched structures. Front. Mater. Sci. 2012, 6, 183–206. [CrossRef] 146. Zhu, X.; Su, W.; Liu, Y.; Hu, B.; Pan, L.; Lu, W.; Zhang, J.; Li, R. Observation of conductance quantization in oxide-based resistive switching memory. Adv. Mater. 2012, 24, 3941–3946. [CrossRef] 147. Hajto, J.; Rose, M.J.; Snell, A.J.; Osborne, I.S.; Owen, A.E.; Lecomber, P.G. Quantised electron effects in metal/a-Si:H/metal thin film structures. J. Non-Cryst. Solids 1991, 137, 499–502. [CrossRef] 148. Samardzic, N.; Mionic, M.; Dakic, B.M.; Hofmann, H.; Dautovic, S.; Stojanovic, G. Analysis of Quantized Electrical Characteristics of Microscale TiO2 Ink-Jet Printed Memristor. IEEE Trans. Electron Devices 2015, 62, 1898–1904. [CrossRef] 149. Yun, E.-J.; Becker, M.F.; Walser, R.M. Room temperature conductance quantization in Vkamorphous-V2 O5 kV thin film structures. Appl. Phys. Lett. 1993, 63, 2493–2495. [CrossRef] 150. Petzold, S.; Piros, E.; Eilhardt, R.; Zintler, A.; Vogel, T.; Kaiser, N.; Radetinac, A.; Komissinskiy, P.; Jalaguier, E.; Nolot, E.; et al. Tailoring the Switching Dynamics in Yttrium Oxide-Based RRAM Devices by Oxygen Engineering: From Digital to Multi-Level Quantization toward Analog Switching. Adv. Electron. Mater. 2020, 6, 2000439. [CrossRef] 151. Zhao, M.; Yan, X.; Ren, L.; Zhao, M.; Guo, F.; Zhuang, J.; Du, Y.; Hao, W. The role of oxygen vacancies in the high cycling endurance and quantum conductance in BiVO4 -based resistive switching memory. InfoMat 2020, 2, 960–967. [CrossRef] 152. Degraeve, R.; Fantini, A.; Clima, S.; Govoreanu, B.; Goux, L.; Chen, Y.Y.; Wouters, D.; Roussel, P.; Kar, G.; Pourtois, G.; et al. Dynamic ‘hour glass’ model for SET and RESET in HfO2 RRAM. In Proceedings of the Symposium on VLSI Technology, Honolulu, HI, USA, 12–14 June 2012. 153. Miranda, E.; Kano, S.; Dou, C.; Kakushima, K.; Suñé, J.; Iwai, H. Nonlinear conductance quantization effects in CeOx /SiO2 -based resistive switching devices. App. Phys. Lett. 2012, 101, 012910. [CrossRef] 154. Cartoixa, X.; Rurali, R.; Suñé, J. Transport properties of oxygen vacancy filaments in metal/crystalline or amorphous HfO2 /metal structures. Phys. Rev. B 2012, 86, 165445. [CrossRef] 155. Zhong, X.; Rungger, I.; Zapol, P.; Heinonen, O. Oxygen modulated quantum conductance for ultra-thin HfO2 -based memristive switching devices. Phys. Rev. B 2016, 94, 165160. [CrossRef] 156. Büttiker, M. Quantized transmission of a saddle-point constriction. Phys. Rev. B 1990, 41, 7906. [CrossRef] 157. Hu, P. One-dimensional quantum electron system under a finite voltage. Phys. Rev. B 1987, 35, 4078. [CrossRef] 158. Senz, V.; Heinzel, T.; Ihn, T.; Lindermann, S.; Held, R.; Ensslin, K.; Wegscheider, W.; Bichler, M. Analysis of the temperature- dependent quantum point contact conductance in view of the metal-insulator transition in two dimensions. J. Phys. Cond. Mat. 2001, 13, 3831. [CrossRef] 159. Miranda, E.; Suñé, J. Electron transport through broken down ultra-thin SiO2 layers in MOS devices. Microelectron. Reliab. 2004, 44, 1–23. [CrossRef] 160. Landauer, R. Spatial variation of currents and fields due to localized scatterers in metallic conduction. IBM J. Res. Dev. 1957, l, 223–231. [CrossRef] 161. Avellán, A.; Miranda, E.; Schroeder, D.; Krautschneider, W. Model for the voltage and temperature dependence of the soft- breakdown current in ultrathin gate oxides. J. Appl. Phys. 2005, 97, 14104. [CrossRef] 162. Miranda, E. The role of power dissipation on the progressive breakdwon dynamics of ultra-thin gate oxides. In Proceedings of the IEEE Proc International Reliability Physics Simposium, Phoenix, AZ, USA, 15–19 April 2007; p. 572. 163. Lombardo, S.; Stathis, J.H.; Linder, B.P.; Pey, K.L.; Palumbo, F.; Tung, C.H. Dielectric breakdown mechanisms in gate oxides. J. Appl. Phys. 2005, 98, 121301. [CrossRef] 164. Stathis, J.H. Percolation models for gate oxide breakdown. J. Appl. Phys. 1999, 86, 5757–5766. [CrossRef] 165. Stathis, J.H. Reliability limits for the gate insulator in CMOS technology. IBM J. Res. Dev. 2002, 46, 265–286. [CrossRef] 166. Dumin, D.J. Oxide Reliability: A Summary of Silicon Oxide Wearout, Breakdown, and Reliability; World Scientific: Singapore, 2002. 167. Linder, B.; Lombardo, S.; Stathis, J.; Vayshenker, A.; Frank, D. Voltage dependence of hard breakdown growth and the reliability implication in thin dielectrics. IEEE Electron Device Lett. 2002, 23, 661–663. [CrossRef] 168. Palumbo, F.; Liang, X.; Yuan, B.; Shi, Y.; Hui, F.; Villena, M.A.; Lanza, M. Bimodal Dielectric Breakdown in Electronic Devices Using Chemical Vapor Deposited Hexagonal Boron Nitride as Dielectric. Adv. Electron. Mater. 2018, 4, 1700506. [CrossRef] 169. Palumbo, F.; Lombardo, S.; Eizenberg, M. Physical mechanism of progressive breakdown in gate oxides. J. Appl. Phys. 2014, 115, 224101. [CrossRef] 170. Tung, C.H.; Pey, K.L.; Tang, L.J.; Radhakrishnan, M.K.; Lin, W.H.; Palumbo, F.; Lombardo, S. Percolation path and dielectric- breakdown-induced-epitaxy evolution during ultrathin gate dielectric breakdown transient. Appl. Phys. Lett. 2003, 83, 2223–2225. [CrossRef] 171. Pey, K.L.; Ranjan, R.; Tung, C.H.; Tang, L.J.; Lin, W.H.; Radhakrishnan, M.K. Gate dielectric degradation mechanism associated with DBIE evolution. In Proceedings of the IEEE International Reliability Physics Symposium, Phoenix, AZ, USA, 25–29 April 2004; Volume 2004, pp. 117–121. Nanomaterials 2021, 11, 1261 45 of 46 172. Privitera, S.; Bersuker, G.; Butcher, B.; Kalantarian, A.; Lombardo, S.; Bongiorno, C.; Geer, R.; Gilmer, D.; Kirsch, P. Microscopy study of the conductive filament in HfO2 resistive switching memory devices. Microelectron. Eng. 2013, 109, 75–78. [CrossRef] 173. Privitera, S.; Bersuker, G.; Lombardo, S.; Bongiorno, C.; Gilmer, D. Conductive filament structure in HfO2 resistive switching memory devices. Solid-State Electron. 2015, 111, 161–165. [CrossRef] 174. Yang, Y.; Gao, P.; Li, L.; Pan, X.; Tappertzhofen, S.; Choi, S.; Waser, R.; Valov, I.; Lu, W.D. Electrochemical dynamics of nanoscale metallic inclusions in dielectrics. Nat. Commun. 2014, 5, 4232. [CrossRef] [PubMed] 175. Pan, C.; Ji, Y.; Xiao, N.; Hui, F.; Tang, K.; Guo, Y.; Xie, X.; Puglisi, F.M.; Larcher, L.; Miranda, E.; et al. Coexistence of Grain- Boundaries-Assisted Bipolar and Threshold Resistive Switching in Multilayer Hexagonal Boron Nitride. Adv. Funct. Mater. 2017, 27. [CrossRef] 176. Rodriguez-Fernandez, A.; Cagli, C.; Perniola, L.; Suñé, J.; Miranda, E. Identification of the generation/rupture mechanism of filamentary conductive paths in ReRAM devices using oxide failure analysis. Microelectron. Reliab. 2017, 76–77, 178–183. [CrossRef] 177. Nishi, Y.; Fleck, K.; Böttger, U.; Waser, R.; Menzel, S. Effect of RESET Voltage on Distribution of SET Switching Time of Bipolar Resistive Switching in a Tantalum Oxide Thin Film. IEEE Trans. Electron Devices 2015, 62, 1561–1567. [CrossRef] 178. Palumbo, F.; Shekhter, P.; Weinfeld, K.C.; Eizenberg, M. Characteristics of the dynamics of breakdown filaments in Al2 O3 /InGaAs stacks. Appl. Phys. Lett. 2015, 107, 122901. [CrossRef] 179. Palumbo, F.; Miranda, E.; Ghibaudo, G.; Jousseaume, V. Formation and Characterization of Filamentary Current Paths in HfO2 -Based Resistive Switching Structures. IEEE Electron Device Lett. 2012, 33, 1057–1059. [CrossRef] 180. Palumbo, F.; Condorelli, G.; Lombardo, S.; Pey, K.; Tung, C.; Tang, L. Structure of the oxide damage under progressive breakdown. Microelectron. Reliab. 2005, 45, 845–848. [CrossRef] 181. Du, H.; Jia, C.-L.; Koehl, A.; Barthel, J.; Dittmann, R.; Waser, R.; Mayer, J. Nanosized conducting filaments formed by atomic-scale defects in redox-based resistive switching memories. Chem. Mater. 2017, 29, 3164–3173. [CrossRef] 182. Palumbo, F.; Eizenberg, M.; Lombardo, S. General features of progressive breakdown in gate oxides: A compact model. In Proceedings of the IEEE International Reliability Physics Symposium, Monterey, CA, USA, 19–23 April 2015; Volume 2015, pp. 5A11–5A16. 183. Lombardo, S.; Wu, E.Y.; Stathis, J.H. Electron energy dissipation model of gate dielectric progressive breakdown in n- and p-channel field effect transistors. J. Appl. Phys. 2017, 122, 085701. [CrossRef] 184. Lombardo, S.; Stathis, J.H.; Linder, B.P. breakdown transients in ultrathin gate oxides: Transition in the degradation rate. Phys. Rev. Lett. 2003, 90, 167601. [CrossRef] [PubMed] 185. Pagano, R.; Lombardo, S.; Palumbo, F.; Kirsch, P.; Krishnan, S.; Young, C.; Choi, R.; Bersuker, G.; Stathis, J. A novel approach to characterization of progressive breakdown in high-k/metal gate stacks. Microelectron. Reliab. 2008, 48, 1759–1764. [CrossRef] 186. Palumbo, F.; Lombardo, S.; Stathis, J.; Narayanan, V.; Mcfeely, F.; Yurkas, J. Degradation of ultra-thin oxides with tungsten gates under high voltage: Wear-out and breakdown transient. In Proceedings of the 2004 IEEE International Reliability Physics Symposium, Phoenix, AZ, USA, 25–29 April 2004; pp. 122–125. [CrossRef] 187. Larcher, L. Statistical simulation of leakage currents in mos and flash memory devices with a new multiphonon trap-assisted tunneling model. IEEE Trans. Electron Devices 2003, 50, 1246–1253. [CrossRef] 188. Nigam, T.; Martin, S.; Abusch-Magder, D. Temperature Dependence and Conduction Mechanism after Analog Soft Breakdown. In Proceedings of the 41st Annual Symposium 2003 IEEE International Reliability Physics, Dallas, TX, USA, 30 March–4 April 2003; pp. 417–423. 189. Condorelli, G.; Lombardo, S.A.; Palumbo, F.; Pey, K.-L.; Tung, C.H.; Tang, L.-J. Structure and conductance of the breakdown spot during the early stages of progressive breakdown. IEEE Trans. Device Mater. Reliab. 2006, 6, 534–541. [CrossRef] 190. Palumbo, F.; Wen, C.; Lombardo, S.; Pazos, S.; Aguirre, F.; Eizenberg, M.; Hui, F.; Lanza, M. A Review on Dielectric Breakdown in Thin Dielectrics: Silicon Dioxide, High-k, and Layered Dielectrics. Adv. Funct. Mater. 2020, 30, 1900657. [CrossRef] 191. Takagi, S.; Yasuda, N.; Toriumi, A. Experimental evidence of inelastic tunneling in stress-induced leakage current. IEEE Trans. Electron Devices 1999, 46, 335–341. [CrossRef] 192. Blöchl, P.E.; Stathis, J.H. Hydrogen electrochemistry and stress-induced leakage current in Silica. Phys. Rev. Lett. 1999, 83, 372–375. [CrossRef] 193. Aguirre, F.L.; RodriguezFernandez, A.; Pazos, S.M.; Sune, J.; Miranda, E.; Palumbo, F. Study on the connection between the set transient in RRAMs and the progressive breakdown of thin oxides. IEEE Trans. Electron. Devices 2019, 66, 1–7. [CrossRef] 194. Ielmini, D. Modeling the universal set/reset characteristics of bipolar RRAM by field-and temperature-driven filament growth. IEEE Trans. Electron Devices 2011, 58, 4309–4317. [CrossRef] 195. Pazos, S.; Aguirre, F.L.; Miranda, E.; Lombardo, S.; Palumbo, F. Comparative study of the breakdown transients of thin Al2 O3 and HfO2 films in MIM structures and their connection with the thermal properties of materials. J. Appl. Phys. 2017, 121, 094102. [CrossRef] 196. Rodríguez-Fernández, A.; Cagli, C.; Sune, J.; Miranda, E. Switching Voltage and Time Statistics of Filamentary Conductive Paths in HfO2 -Based ReRAM Devices. IEEE Electron Device Lett. 2018, 39, 656–659. [CrossRef] 197. Zafar, S.; Jagannathan, H.; Edge, L.F.; Gupta, D. Measurement of oxygen diffusion in nanometer scale HfO2 gate dielectric films. Appl. Phys. Lett. 2011, 98, 152903. [CrossRef] Nanomaterials 2021, 11, 1261 46 of 46 198. Shacham-Diamand, Y.; Dedhia, A.; Hoffstetter, D.; Oldham, W.G. Copper Transport in Thermal SiO2 . J. Electrochem. Soc. 1993, 140, 2427–2432. [CrossRef] 199. Nason, T.C.; Yang, G.; Park, K.; Lu, T. Study of silver diffusion into Si(111) and SiO2 at moderate temperatures. J. Appl. Phys. 1991, 70, 1392–1396. [CrossRef] 200. Kim, B.G.; Yeo, S.; Lee, Y.W.; Cho, M.S. Comparison of diffusion coefficients and activation energies for Ag diffusion in silicon carbide. Nucl. Eng. Technol. 2015, 47, 608–616. [CrossRef] 201. Zobelli, A.; Ewels, C.P.; Gloter, A.; Seifert, G. Vacancy migration in hexagonal boron nitride. Phys. Rev. B 2007, 75, 094104. [CrossRef]