micromachines Article SPICE Implementation of the Dynamic Memdiode Model for Bipolar Resistive Switching Devices Fernando Leonel Aguirre * , Jordi Suñé and Enrique Miranda * Departament d’Enginyeria Electrònica, Universitat Autònoma de Barcelona (UAB), 08193 Barcelona, Spain; jordi.sune@uab.cat * Correspondence: aguirref@ieee.org (F.L.A.); enrique.miranda@uab.cat (E.M.) Abstract: This paper reports the fundamentals and the SPICE implementation of the Dynamic Memdiode Model (DMM) for the conduction characteristics of bipolar-type resistive switching (RS) devices. Following Prof. Chua’s memristive devices theory, the memdiode model comprises two equations, one for the electron transport based on a heuristic extension of the quantum point-contact model for filamentary conduction in thin dielectrics and a second equation for the internal memory state related to the reversible displacement of atomic species within the oxide film. The DMM represents a breakthrough with respect to the previous Quasi-static Memdiode Model (QMM) since it describes the memory state of the device as a balance equation incorporating both the snapback and snapforward effects, features of utmost importance for the accurate and realistic simulation of the RS phenomenon. The DMM allows simple setting of the initial memory condition as well as decoupled modeling of the set and reset transitions. The model equations are implemented in the LTSpice simulator using an equivalent circuital approach with behavioral components and sources. The practical details of the model implementation and its modes of use are also discussed.  Keywords: memristor; resistive switching; memory; memdiode  Citation: Aguirre, F.L.; Suñé, J.; Miranda, E. SPICE Implementation of the Dynamic Memdiode Model for Bipolar Resistive Switching Devices. 1. Introduction Micromachines 2022, 13, 330. https:// Selecting a suitable generic model for an electron device is far from being simple doi.org/10.3390/mi13020330 and straightforward. The model should be able to cover not only the basic and common Academic Editors: Andrey Sokolov features of the system under study but must also have the necessary capacity of adaptation and Haider Abbas for incorporating its distinctive behavior. For circuit simulation-oriented models, this adaptation capability must be achieved by means of a few numbers of simple and robust Received: 21 January 2022 equations driven by a reduced set of parameters, if possible, with physical origin, if not, with Accepted: 14 February 2022 some degree of electrical meaning. This is the signature of a compact behavioral approach, Published: 19 February 2022 in which the central objective is not to outperform the microscopic level models (for instance Publisher’s Note: MDPI stays neutral the kinetic Monte-Carlo and Finite Element models such as those reported in [1–4]) in terms with regard to jurisdictional claims in of accuracy and fidelity to the device physics, but matching observations and simulations published maps and institutional affil- as close as possible. Clearly, accurate representation of the electron transport across the iations. investigated device under arbitrary input signals is two-fold: first, it encourages the design and assessment of more complex circuits and systems, and second, it allows to identify and establish links among the elementary modeling pieces that lead to the variety of observed behaviors (conduction characteristics). Since the first practical description of a Copyright: © 2022 by the authors. memristive device by HP in 2008 [5], a number of compact models for the current-voltage Licensee MDPI, Basel, Switzerland. (I-V) curves of resistive switching (RS) devices has been proposed [6]. Because of the This article is an open access article hysteretic nature of the phenomenon (see Figure 1a), the electron transport model involves distributed under the terms and complementary information about the previous history of the device. Here is where Prof. conditions of the Creative Commons Chua’s theory of memristive devices comes into play [7]. According to this theory, a Attribution (CC BY) license (https:// physical or electrical variable expressed as a first order differential equation controls the creativecommons.org/licenses/by/ 4.0/). current flow through the structure. Based on this framework, a plethora of memristor Micromachines 2022, 13, 330. https://doi.org/10.3390/mi13020330 https://www.mdpi.com/journal/micromachines Micromachines 2022, 13, 330 2 of 18 models have been proposed, such as the general phenomenological models (Yakopcic [8], TEAM [9], VTEAM [10], Eshraghian [11], etc.). Despite having a proven capability to successfully fit experimental data, these models rely on various internal equations or artificial window functions (commonly used for modelling the SET/RESET transitions) that can seriously affect the model’s convergence [12,13]. The reader can find further details Micromachines 2022, 11, x FOR PEER for each model in the corresponding paper or in review papers [11,13–18]. REVIEW 3 of 20 b) c) d) a) Figure 1. (a) Hysteretic behavior of the filamentary-type I-V characteristic. Filament stages: (A) Figure 1. (a) Hysteretic behavior of the filamentary-type I-V characteristic. Filament stages: (A) bar- barrier-controlled, high resistance state (HRS), (B) completion, (C) expansion, (D,E) low resistance rier-controlled, high resistance state (HRS), (B) completion, (C) expansion, (D,E) low resistance state state (LRS), (F) partial dissolution, (G) rupture and (H) HRS (b) Schematic of the filamentary structure: (LRS), (F) partial dissolution, (G) rupture and (H) HRS (b) Schematic of the filamentary structure: fixed (red) and dynamic (green, blue) sections. (c) Schematic of the tunneling barrier/gap. λ is the fixed (red) and dynamic (green, blue) sections. (c) Schematic of the tunneling barrier/gap. λ is the barrier barrier height. height. (d) (d) Representation Representation of of the the QPC QPC model model using using an an equivalent equivalent circuit circuit approach. approach. Many of these models almost exclusively focus on the popular quasi-static, pinched 2. Dynamic Memdiode Model (DMM) I-V loop, ignoring the associated time-related dependencies. However, the latter has special In this relevance Section, when the two real considering equations that definescenarios case application the DMMsuch are as presented and discussed. those described in [19], They are: (i) the current-voltage (I-V) relationship, and (ii) where the programming and reading of the device is made in terms of voltage the memory state equation pulses(λ- of t). The very varying basic idea behind frequency/duty cycle. the DMM In this is that regard, wetheaimcurrent to report flows the through a kindand fundamentals of fila- the mentary SPICE structure embedded implementation in the oxide of a revised Dynamic layerMemdiode in which some Model of (DMM) its atomic forconstituents bipolar RS can reversibly devices capablemove in and out according of incorporating to the forces the time-related exerted by dependencies, asthe well external field. This as the guidelines changing for configuration its usage. Since thisalters new in turn the version overall transmission incorporates a dynamicproperties of the structure, balance equation for the leading the memory system state and atohigher a hysteretic level ofbehavior. details Since thisof in terms is modelling at the end aaccuracy, behavioral weapproach, consider the model that can be appropriately it is a breakthrough modified with respect soprevious to the as to cover non-filamentary-type models proposed by ourconduction group: the as well. Quasi-static Memdiode Model (QMM) [20–22] and a first version of the DMM [19]. The first one relies on the double-diode circuit controlled by the Krasnosel’skii-Pokrovskii hysteresis 2.1. Current-Voltage operator Characteristic [23]. The main advantage of that simplified approach is the elimination of the time According integrationtostepthe quantumofpoint-contact in favor the use of the so-called (QPC) modelhysteron [25–28], orthememory mapflows current that (i.e., the low-voltage conductance exhibited by the device as a function through a single nanosized filamentary structure (see Figure 1b) can be calculated using of the applied signal). However, evenLandauer the finite-bias though a time module formula for acan be added conductor monomode to the base [29]: model [24], this is in essence a quasi-static approach in the sense that the memory state of the device does not change unless a hard threshold) condition 2𝑒 +∞ 𝑒𝑉𝐶 be used for𝑒𝑉 𝐶 𝐼(𝑉𝐶 = ∫ is 𝑇(𝐸) met. The QMM [𝑓 (𝐸 − can ) − 𝑓 (𝐸 + arbitrary )] 𝑑𝐸 input signals, (1) but the memory state is unable ℎ −∞ 2 to evolve in the subthreshold regime.2 Seeking to overcome this limitation, we introduced the original DMM in [19]. where 𝑉𝐶 = 𝑉 − 𝐼𝑅𝐼 is the potential drop across the constriction, V the applied voltage, 𝑅𝐼 In this work, we presented a revised version of the DMM. Compared to the previous the internal/external series resistance (permanent section of the filament/wire resistance), version, the current one incorporates modifications capable of representing a variety of E the energy, f the Fermi-Dirac function, and T the transmission coefficient for the con- details that could not be modeled by its predecessor. These are: (i) the simulation of the finement barrier. e and h are the electron charge and the Planck constant, respectively. electroforming event, which implies a first cycle different from its subsequent stationary Equation (1) can be extended to the case of asymmetric potential drops at the two ends of cycles, (ii) the modelling of the snapback (SB) effect, experimentally observed in memristors the constriction using a coefficient different from ½ for the energy window. Assuming an inverted parabolic potential barrier for the constriction’s bottleneck (scatterer), T(E) is ex- pressed as: 1 𝑇(𝐸) = (2) Micromachines 2022, 13, 330 3 of 18 during the SET event, (iii) the modelling of the snapforward (SF) event (low current tail) during the RESET event, and (iv) the de-coupling of the SET and RESET equations, which simplifies the fitting. The paper is organized as follows: Section 2 discusses the physical background of the model for both the transport and memory equations, as well as for the resulting switching dynamics. The dependences of the model with the fitting parameters are evaluated in Section 3, focusing on the memory equation, the switching dynamics, and the particular case of Complementary Resistive Switching (CRS) devices. Then, Section 4 presents an overview of the DMM capability for reproducing different sets of experimental data, from the quasi-static I-V loops as well as from the potentiation and depression dynamics for neuromorphic applications. Finally, Section 5 reports the conclusions of this work. 2. Dynamic Memdiode Model (DMM) In this Section, the two equations that define the DMM are presented and discussed. They are: (i) the current-voltage (I-V) relationship, and (ii) the memory state equation (λ-t). The very basic idea behind the DMM is that the current flows through a kind of filamentary structure embedded in the oxide layer in which some of its atomic constituents can reversibly move in and out according to the forces exerted by the external field. This changing configuration alters in turn the overall transmission properties of the structure, leading the system to a hysteretic behavior. Since this is at the end a behavioral approach, the model can be appropriately modified so as to cover non-filamentary-type conduction as well. 2.1. Current-Voltage Characteristic According to the quantum point-contact (QPC) model [25–28], the current that flows through a single nanosized filamentary structure (see Figure 1b) can be calculated using the finite-bias Landauer formula for a monomode conductor [29]: 2e +∞      eV eV Z I (VC ) = T ( E) f E − C − f E + C dE (1) h −∞ 2 2 where VC = V − IR I is the potential drop across the constriction, V the applied voltage, R I the internal/external series resistance (permanent section of the filament/wire resis- tance), E the energy, f the Fermi-Dirac function, and T the transmission coefficient for the confinement barrier. e and h are the electron charge and the Planck constant, respectively. Equation (1) can be extended to the case of asymmetric potential drops at the two ends of the constriction using a coefficient different from 12 for the energy window. Assuming an inverted parabolic potential barrier for the constriction’s bottleneck (scatterer), T(E) is expressed as: 1 T ( E) = (2) 1 + e θ ( E− ϕ) − where θ is the barrier shape factor and ϕ the barrier height (see Figure 1c). In the pure tunneling regime ( ϕ  E) and zero-temperature limit, Equations (1) and (2) yield [30]:   4e −θ ϕ θeVC I (VC ) ≈ e sinh (3) θh 2 Notice that, if the barrier width collapses ( θ → 0 ) because of the completion of the filament, Equation (3) results in the standard Landauer formula I = G0 VC , where G0 = 2e2 /h is the quantum conductance unit. The effect of θ on the I-V curve is illustrated in Figure 2a. As θ increases, the current decreases, and the I-V changes its dependence from linear to linear- exponential. Combining Equation (3) in the linear regime (eVc  2/θ) with Equation (2), we obtain: 1 1 − T ( E = 0)   1 θϕ VC ∼ e I= I (4) G0 G0 T ( E = 0) is the quantum conductance unit. The effect of 𝜃 on the I-V curve is illustrated in Figure 2a. As 𝜃 increases, the current decreases, and the I-V changes its dependence from linear to linear-exponential. Combining Equation (3) in the linear regime (𝑒𝑉𝑐 ≪ 2⁄𝜃 ) with Equa- tion (2), we obtain: Micromachines 2022, 13, 330 1 𝜃𝜑 1 1 − 𝑇(𝐸 = 0) 4 of 18 𝑉𝐶 ~ 𝑒 𝐼= ( )𝐼 (4) 𝐺0 𝐺0 𝑇(𝐸 = 0) Since in a mesoscopic system, the constriction resistance 𝑅 can be regarded as the sum of Since in a mesoscopic system, the constriction resistance RC𝐶 can be regarded as the sum of the contact contact resistance 𝑅𝐿 and resistance R the scatterer resistance 𝑅𝐵 as [29]: the L and the scatterer resistance R B as [29]: 1 1 =1 1 +1 1 (11−  − 𝑇 𝑅 = T)=𝑅 +𝑅 (5) RC = 𝐶 𝐺 =𝑇 𝐺+ 𝐺 = R𝐿L + R𝐵B G0 T 0 G0 0 G0 0 T𝑇 we can associate Equation (4) with 𝑅 . However, it is clear from Equation (5) that not only we can associate Equation (4) with 𝐵R B . However, it is clear from Equation (5) that not the confinement barrier contributes to the constriction resistance but also the way the con- only the confinement barrier contributes to the constriction resistance but also the way the striction is attached constriction to the is attached charge to the reservoirs charge (or(or reservoirs thermalizing thermalizingregion) region) through𝑅𝐿R.LThis through is isa . This well-known result [29] and is the consequence of the funneling effect of the electron a well-known result [29] and is the consequence of the funneling effect of the electron wave wave function caused function caused byby the mismatch in the mismatch the number in the number ofof available available energy energy states states when when passing passing from the reservoir to the constriction and vice versa. from the reservoir to the constriction and vice versa. 10-3 10-2 a) G0=2e2/h b) R=1K-1 I0min=1E-7A 10 -4 I=G0V 10 -3 =3eV I0max=1E-3A 10-5 10-4 10-6 10-5 Current [A] 10-7 10-6 10-8 10-7   -9 8 eV-1 -8 0 10 10  6 eV-1  1E-3 4 eV-1 1E-2 10-10 10-9 1E-1 2 eV-1 1E-3 eV-1 1E0 10-11 10-10 0.01 0.1 1 2 0.01 0.1 1 2 Voltage [V] Voltage [V] Figure 2. (a) Single filament conduction characteristic and effect of the parameter θ. θ. (b) Multiple filament conduction filament conduction characteristic characteristic and and effect effect of of the the parameter parameter λ. λ. Since, in general, for a wider constriction formed by a bunch of conducting channels neither the number N of elemental filamentary structures involved involved is known nor their specific barrier parameters 𝜃θ𝑖i and 𝜑 can be accessed individually [31], we extrapolate ϕ𝑖i can extrapolate Equation (3) to that case using the heuristic approximation: approximation: 𝐼(𝑉 I (V ) =I0𝐼sinh C )𝐶 = 0 𝑠𝑖𝑛ℎ[𝛼(𝑉 [α(VC 𝐶−−R𝑅S𝑆I𝐼)] )] (6) (6) which has the same functional asymptotes as the original Equation (3) for large 𝐼0 values which has the same functional asymptotes as the original Equation (3) for large I0 values (𝑅𝑆 𝐼 ≈ 𝑉𝐶 ) and low applied voltages (𝑅𝑆 𝐼 ≪ 𝑉𝐶 ) (see Figure 2b). While 𝑅𝑆 in Equation (6) ( RS I ≈ VC ) and low applied voltages ( RS I  VC ) (see Figure 2b). While RS in Equation (6) accounts for the contact resistance, the hyperbolic sine function expresses the barrier re- accounts for the contact resistance, the hyperbolic sine function expresses the barrier sistance as illustrated in Figure 1d. 𝑅𝑃𝑃 in the same figure deals with the device resistance resistance as illustrated in Figure 1d. R PP in the same figure deals with the device resistance before the forming event. A central difference between 𝑅𝑆 and 𝑅𝐼 , is that 𝑅𝑆 will be allowed before the forming event. A central difference between RS and R I , is that RS will be allowed to change (if necessary) according to the memory state of the device (movement of ions/vacancies). Notice also that Equation (6) does not correspond strictly to N times the current flowing through a single filament, otherwise a parallel shift of Equation (3) towards higher current values would be obtained [32]. In addition, Equation (6) complies with the pinched condition I (V = 0) = 0 and, because of its heuristic nature, overcomes the physical limitation on the voltage drop eVC /2 ≤ ϕ imposed by Equation (3). Physically, the moving species in Figure 1b represent the hopping of ions/vacancies induced by the external applied field. As schematically illustrated in Figure 1c, the opening (RESET) or closing (SET) of the atom chain raises or lowers the top of the confinement barrier for the electron flow [33,34]. According to this picture, ballistic transport would not be required along the whole filament structure but just at the narrowest section of the constriction. For a complete discussion about conductance quantization effects in RS devices, see Refs. [35–43]. Micromachines 2022, 13, 330 5 of 18 From the electrical viewpoint, Equation (6) can be envisaged as two opposite biased diodes in parallel with a single series resistance [6] (see Figure 1d). Since these are not real diodes, inverse saturation currents are disregarded. This was the approach followed in our previous works and justifies the origin of the name memdiode [23], i.e., a diode with memory. The main point is that Equation (6), as shown in Figure 2b, self-rectifies the I-V curve as the device switches from HRS (linear-exponential) to LRS (linear), which is in agreement with many experimental observations. Equation (6) is implemented in LTSpice XVII from Linear Technologies using two resistors (RS and R I ) in series with a behavioral voltage-controlled hyperbolic sine current generator with the amplitude factor given by the expression:   I0 (λ) = Ion − Io f f λ + Io f f (7) where 0 ≤ λ ≤ 1 is the memory state variable, and Io f f , Ion are calibration parameters (OFF and ON currents). λ = 0 and λ = 1 correspond to HRS and LRS, respectively. The linear relationship between I0 and λ is a key feature of the model and likely reflects the connection between the memory state and the density of conducting sites in the filamentary structure [44,45]. The effect of λ on the I-V curve is illustrated in Figure 2b. As λ increases, the I-V curve becomes more linear as expected for a fully formed conducting channel. Notice that, at low biases, the HRS I-V is in the linear tunneling regime (which must not be confused with ohmic conduction). Importantly, Io f f and Ion control the barrier resistance and do not refer directly to the minimum and maximum currents that flow through the device. This will be ultimately determined by the whole system’s dynamics (hysteresis effect and series resistances). For the sake of completeness, α and RS in Equation (6) receive a similar treatment in the LTSpice script as that given to I0 (λ). Both parameters can be swept from a minimum (OFF) to a maximum (ON) if required. If not, α and RS remain fixed. In the following Section, the memory state equation and its circuital implementation are discussed. 2.2. Memory State Equation As reported in [46], a very convenient differential equation for the memory state variable λ that complies with a number of experimental observations in memristive struc- tures is: dλ 1−λ λ = − (8) dt τS (λ, VC ) τR (λ, VC ) where τS,R are characteristic times associated with the SET (V > 0) and RESET (V < 0) transitions, i.e., with the ionic/defect movement within the dielectric film in one or the opposite direction. Equation (8) can be simply regarded as the normalized version of a birth-death process for a two-state system with transition rates τS−1 and τR−1 : τS−1 X Y (9) τR−1 in which there are n1 particles in the state X and n2 particles in the state Y, with n1 + n2 = N the total number of particles. This comes to represent for example the REDOX process in VCMs [47]. Notice that τS,R in Equation (8) are expressed as a function of VC and λ. In our case, where SET and RESET only occur for biases with opposite signs, Equation (8) can be treated as two separate differential equations, one for V > 0 and one for V < 0. This is not mandatory but simplifies the model calibration since the SET and RESET processes are completely disentangled. Under this consideration, Equation (8) can be represented by the equivalent circuit schematically depicted in Figure 3. λ corresponds to the voltage drop across the capacitor C = 1 F. Notice that the memory state behavior during the RS cycle is nothing but the alternate action of two current sources that charge and discharge a capacitor. The green and red arrows in Figure 3 indicate the position of the switches as a function of the sign of VC . In practice, the switches are modeled by a behavioral if can be treated as two separate differential equations, one for V > 0 and one for V < 0. This is not mandatory but simplifies the model calibration since the SET and RESET processes are completely disentangled. Under this consideration, Equation (8) can be represented by the equivalent circuit schematically depicted in Figure 3. λ corresponds to the voltage drop across the capacitor C = 1 F. Notice that the memory state behavior during the RS Micromachines 2022, 13, 330 cycle is nothing but the alternate action of two current sources that charge and discharge 6 of 18 a capacitor. The green and red arrows in Figure 3 indicate the position of the switches as a function of the sign of VC. In practice, the switches are modeled by a behavioral if state- ment in the statement LTSpice in the LTSpicescript (Algorithm script (Algorithm1). 1). TheTheinitial condition initial forfor condition thethe memory memory state stateis introduced is introducedthrough throughthe theinitial initialvoltage voltagedrop dropacross across the capacitor as V(t = the capacitor = 0) 0) == λλ00. In In the the LTSpice LTSpice script, script, λλisisrepresented represented by by the thecapital capital letter letterH H(for (forhysteron). hysteron). In In fact, fact, λλ== V(H) V(H) is is the the voltage voltageatatthe theHHnode. node.ForForthethesake sakeofofcompleteness, completeness, and andtoto help thethe help reader reader in the in theuseuseof the DMM, of the DMM, Supplementary Supplementary Figures FiguresS1–S3 1–3 illustrates illustrates three three simulation simulation exercises exercises inin LTSpice LTSpice obtained with the script reported in Algorithm obtained with the script reported in Algorithm 1. 1. Figure 3.3. Equivalent Figure Equivalent circuit circuit model model for for the the memory memory state state equation equation using using two two behavioral behavioral current current sourcesand sources andone onecapacitor. capacitor.The switchingτ𝜏𝑆 and Theswitching andτ𝜏𝑅 can canbe be regarded regardedas asvariable variableresistances. resistances.The Thestate state S R of the switches depends on the sign of the applied voltage. The memory state λ is the voltage of the switches depends on the sign of the applied voltage. The memory state λ is the voltage across across the capacitor. λ0 is the initial memory state. the capacitor. λ0 is the initial memory state. Algorithm Algorithm 1: 1: Memdiode Memdiodescript scriptforfor LTSpice LTSpice XVII. XVII. + and+− and are−the aredevice the device terminals. terminals. H is theH is the memory memory statestate output. output. The colors The colors indicate indicate the different the different sections: parameter sections: parameter values, values, memory equation, I-V characteristic, and auxiliary functions. memory equation, I-V characteristic, and auxiliary functions. 11 .subckt memdiode+ +−−H .subckt memdiode H 22 *created *created bybyE.Miranda, E.Miranda, F. F. Aguirre andand Aguirre J.Suñé, revised J.Suñé January , revised 2022 2022 January 3 .params 3 .params 4 + H0 = 0 ri = 50 RPP = 1E10 4 + H0 = 0 ri = 50 RPP = 1E10 5 + etas = 50 vs = 1.4 56 ++etas etar == 50 100vsvr==1.4 −0.4 67 ++etar ion = 1E-2 vr = 100 aon= =−0.4 2 ron = 10 78 ++ion ioff==1E-2 aon == 22 roff 1E-7 aoff ron==10 10 9 + vt = 0.4 isb = 2E-4 gam 8 + ioff = 1E-7 aoff = 2 roff = 10 = 1; isb = 1/gam = 0 no SB/SF 10 *Memory Equation 9 + vt = 0.4 isb = 2E-4 gam = 1; isb = 1/gam = 0 no SB/SF 11 BI 0 H I = if(V(+,-)> = 0, (1-V(H))/TS(V(C,-)),-V(H)/TR(V(C,-))) 10 *Memory Equation 12 CH H 0 1 ic = {H0} 11 13 BI 0 H I = if(V(+,-)> = 0, (1-V(H))/TS(V(C,-)),-V(H)/TR(V(C,-))) *I-V 12 14 RI +HC 0{ri} CH 1 ic = {H0} 13 15 *I-V RS C B R = K(ron,roff) 16 BF B - I = K(ion,ioff)*sinh(K(aon,aoff)*V(B,-)) 17 RB + - {RPP} 18 *Auxiliary functions 19 .func K(on,off) = off+(on-off)*limit(0,1,V(H)) 20 .func TS(x) = exp(-etas*(x- if(I(BF)>isb,vt,vs))) 21 .func TR(x) = exp(etar* if(gam = = 0,1,pow(limit(0,1,V(H)),gam))*(x-vr)) 22 .ends Both in the SET and RESET regions, the corresponding characteristic switching times can depend implicitly or explicitly on λ. This property is used to include the so-called snapback (positive bias) and snapforward (negative bias) effects in the RS I-V loop. These effects are typically present in VCM-based structures [48]. In this work, we introduce the memory state λ in the characteristic times as: τS (λ, VC ) = e−ηS (VC −VS (λ)) (10) Micromachines 2022, 13, 330 7 of 18 and γ (V −V ) τR (λ, VC ) = eηR λ C S (11) where ηS,R and VS,R are the transition rates (ηS , ηR > 0) and the reference switching voltages (VS > 0, VR < 0), respectively. γ ≥ 0 is referred to as the SF coefficient. The exponential dependences of Equations (10) and (11) on VC are a consequence of the ions/vacancies dy- namics associated with the hopping mechanism [49,50]. Deviations from these exponential laws in the low and high voltage regions have also been reported but are disregarded in this work [51]. If for any reason, the SB and SF effects do not need to be considered, taking VS (λ) = VS a constant reference SET voltage and λγ = 1 in Equations (10) and (11), respectively, the switching dynamics becomes exclusively voltage-controlled, as originally assumed in [46]. This behavior is typical of ECM cells in which abrupt RESET transitions are observed [47]. Under these latter conditions, Equation (8) has an analytic solution both for the constant and ramped voltage input signal cases. In any other case, because of the mathematical complexity involved, Equation (8) must be numerically solved with the help of a differential equation solver (in our case the circuit simulator itself). In the following Section, the practical implementation and the consequences of the above-mentioned effects on the I-V curve are discussed. 3. Simulation Results and Discussion In order to test the ability of the proposed model to deal with realistic simulations, a number of evaluation criteria must be adopted and assessed. In this work, we basically consider Linn’s criteria [52] to which we add some very specific features not included in the referred work. These criteria are: (i) capability of the compact model to reproduce the RS I-V characteristics including the SB and SF effects, (ii) realistic switching dynamics for the SET and RESET transitions including the ability of the model to deal with arbitrary input signals (continuous and discontinuous), and (iii) multi-device connectivity in the form of Complementary Resistive Switching (CRS). These three major issues are discussed in detail next. Before entering the discussion, it is worth mentioning that because of the complexity of the numerical problem involved, caution should be exercised with the selection of the model parameter values, as it happens with any other model. Although the DMM is robust enough, a certain combination of parameters could lead to fatal errors of convergence or to extremely long simulation times. Sometimes the numerical problems disappear by simply changing the maximum simulation timestep (shorter or longer) or numerical integration method used (trapezoidal, modified trap, gear). Depending on the required accuracy, simplifying the model equations by eliminating unnecessary details (SB, SF, resistances, limiting functions, etc.) is also a good strategy to follow. The roles played by parameters I0 , α, and RS in the HRS and LRS I-V curves are not analyzed here since they were discussed elsewhere in connection with the QMM [24]. Notice that the QMM can always be used as the starting point of any simulation exercise. The main difference between both models resides in how the SET and RESET transitions are modeled. While QMM requires a threshold voltage/current to induce the memory state evolution, DMM does not. The I-V expression in Equation (6) is common to both models. 3.1. Memory State Equation To begin with, Figure 4 illustrates typical I-V and λ-V loops obtained using Equations (6) and (8). The LTSpice script and model parameters for this particular exercise can be found in Algorithm 1. While in Figure 4a, the sinusoidal input voltage, the memory state, and the current flowing through the structure are plotted as a function of time, Figure 4b illustrates the current evolution and the memory map (hysteron) of the device as a function of the applied voltage [53]. More in detail, the SB effect is recognized by the sudden current increase in the SET region (red line in Figure 4c) caused by the reduction of the constriction resistance that occurs when the tunneling gap or confinement barrier vanishes (the CF is completely formed). During this phase the current also grows as a function (6) and (8). The LTSpice script and model parameters for this particular exercise can be found in Algorithm 1. While in Figure 4a, the sinusoidal input voltage, the memory state, and the current flowing through the structure are plotted as a function of time, Figure 4b illustrates the current evolution and the memory map (hysteron) of the device as a func- tion of the applied voltage [53]. More in detail, the SB effect is recognized by the sudden Micromachines 2022, 13, 330 current increase in the SET region (red line in Figure 4c) caused by the reduction of8the of 18 constriction resistance that occurs when the tunneling gap or confinement barrier van- ishes (the CF is completely formed). During this phase the current also grows as a function ofofVVC but C butfollowing following approximately approximately thetheload line load of of line thethe circuit (slope~1/𝑅 circuit (slope~1/R𝐼 ), and nextnext I ), and at anat almost an almost constant voltage called the transition voltage VT (blue line in Figure 4c)This constant voltage called the transition voltage V T (blue line in Figure 4c) [48]. [48]. second phase phase This second corresponds to the accumulation corresponds to the accumulationof ions/defects in the constriction of ions/defects in the constriction(or alter- (or natively to its to alternatively lateral expansion) its lateral withwith expansion) the the consequent consequent progressive progressive resistance resistance reduction. reduction. This Thisbehavior behaviorhas hasbeen beenreported reportedmany manytimestimesininthetheliterature literature[54] [54]but buthas hasreceived receivedscarcescarce attention attentionininthe thecompact compactsimulation simulationfield. field.VV isisthe TT theminimum minimumvoltagevoltagerequired requiredtotoactivate activate the theions/vacancies ions/vacancies movement movementand anditsitsvalue valueseems seemstotobebenot notonly onlya characteristic a characteristicparameter parameter ofofeach eachmaterial material butbut also a function function ofofmeasurement measurementvariables variablessuchsuchasas thethe voltage voltage ramp ramp rate or signal rate frequency or signal frequency [55].[55]. 1.5 10 V0=1.2 V a) 12 1.0 b) 0.010 Model I-V S1/RI LRS f=1 Hz I-VC 1.0 Snapback I-VC Voltage [V], emory state  8 0.8 5 0.5 0.005 4 Current [mA] I-V Current [A] 0.6 0.0 0 0 HRS VR ISB 0.000 -4 0.4 HRS VT VS -0.5 -8 0.2 -5 -0.005 -1.0 memory state input voltage -12 0.0 LRS current c) S1/RI -1.5 -10 -0.010 0.0 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 Time [s] Voltage [V] Voltage [V] Figure Figure4.4.(a) (a)Simulation Simulationexample examplefor for the the current andmemory current and memorystate stateasasa afunction function ofof time time of of a sinus- a sinusoidal oidal signal. (b) Evolution of the memory state (hysteron) and current as a function of the appliedapplied signal. (b) Evolution of the memory state (hysteron) and current as a function of the voltage. voltage. (c) Original I-V curve (red line) and its snapback correction I-VC (blue line). VT is the transi- (c) Original I-V curve (red line) and its snapback correction I-VC (blue line). VT is the transition tion voltage, VR the reset voltage, VS the set voltage, RI the internal series resistance, and ISB the voltage, VR the reset voltage, VS the set voltage, RI the internal series resistance, and ISB the snapback snapback triggering current. triggering current. The SB effect is incorporated into the model equations by modifying the SET refer- The SB effect is incorporated into the model equations by modifying the SET reference ence voltage VS (>VT) in (10) according to the rule: voltage VS (>VT ) in (10) according to the rule: ( VT I ≥ ISB VS ( I ) = (12) VS I < ISB where I is the current flowing through the device. Equation (12) is a switching rule based on the current value, but other rules related to the applied voltage, dissipated power, or memory state, are also admissible [56]. ISB is a threshold current for the SB effect. Equation (12) is written as an if statement for the SET voltage in the LTSpice script and expresses a collapse of the nominal SET voltage VS to a lower value VT after reaching the threshold condition ISB . This event generates a sudden current increase compatible with the voltage drop along the load line of the circuit. It is worth mentioning that the SB effect is not always observable in practice since its detection depends on a number of factors linked to the specific features of the device under test and to the measurement conditions (current magnitude, current compliance, etc.). When combined with other parameters (Ioff and RPP ), VS can also be used to represent the forming step (see the Supplementary Material). This may require code edition for a specific conduction mechanism (Schottky, Fowler-Nordheim, etc.) in the fresh device [57]. For the opposite polarity (V < 0), after the SF event (current decrease following the circuit load line with slope~1/R I ), the main difference appears at the low current region, once the filament is almost dissolved. In this case, since λ approaches zero as the current Micromachines 2022, 13, 330 9 of 18 drops, the factor λγ gains weight in Equation (11), reducing the RESET characteristic time. The result is remarkable since the current deviates from the load line generating a lobe. In other words, as the current decreases, larger voltages are required to deplete the constriction from conducting atomic species up to the point in which the initial gap or tunneling barrier is completely restored. The referred protuberance is clearly visible in many VCM-type devices but is rarely observed in ECM-type structures, which exhibit more abrupt transitions [58,59]. Although VR is considered an independent model parameter in the LTSpice script, in general VR = −VT is found, which is consistent with a field-induced activation of the SET/RESET processes in bipolar RS devices. Figure 5 illustrates the effects of some of the model parameters on the I-V curve. The analysis is carried out on the second I-V loop, i.e., once the transient effects associated with the initial loop or forming process plays no role. As shown in Figure 5a, R I mainly affects the slope of the LRS I-V curve and the apparent RESET voltage. The small shift in the SET voltage is a consequence of the modifications that occurred in the RESET region of the first loop (not shown). Remarkable changes are also observed in the I-VC curve (see below). Figure 5b illustrates the effects of the threshold current ISB . As this parameter increases, the completion of the filament takes place at a higher voltage, thus reducing the observable effects associated with VT . No change is detected in the RESET transition since the LRS I-V remains unaltered. Figure 5c illustrates the effects of the SF parameter γ. The main effect on the RESET transition is the change of the triggering point of the current lobe. Since the HRS current in the RESET region is also affected by this change, γ also alters the triggering point in the SET region. It is also important to mention that the observation of the SB effect in the simulated curve strongly depends on the memory state initial condition (λ0 ). As shown in Figure 6a, depending on the HRS current magnitude, the SB triggering point can differ (because the same current value is reached at different voltages). Once the device Micromachines 2022, 11, x FOR PEER reaches REVIEW LRS, the RESET process becomes independent of the initial condition. Subsequent 10 of 20 loops do not carry on information about the initial state of the device. 10 15 10 a) b) 1E-4 A c) 1.2 2E-4 A ISB RI 8 3E-4 A 8  1.1 1.0 1 4E-4 A 0.9 25  10 50  Current [mA] Current [mA] Current [mA] 6 6 75  RI 4 4  5 ISB 2 2 0 0 0 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 Voltage [V] Voltage [V] Voltage [V] Figure 5. Figure 5. Effect Effect of of the the DMM DMM parameters parameterson onthe theI-V I-Vcharacteristic: characteristic:(a)(a) internal/external series internal/external re- series sistance R , (b) snapback triggering current I , and (c) snapforward coefficient resistance RI , (b) snapback triggering current ISB , and (c) snapforward coefficient λ. I SB λ. AA remarkable remarkable property property of of the theI-V I-VCC curve, curve, which which results results from from thethe SB SB transformation transformation (V-R· (V-R ·I) I) of of the original I-V the original I-V curve (see Figure 6b,c), 6b,c), isis that in in addition addition to to the the current current increase increase at at aa constant voltageVVTT occurring constantvoltage occurring in in the theSET SETregion, region,the minorI-V theminor I-VCCloops loopsalso peakatat-V alsopeak -VTT in in the the RESET RESET region region (see (see Figure Figure 6d) for V 6d) for −V VRR ==−V . This T. TThis hashas been been experimentally experimentally verified verified in in [60][60] and and indicates indicates thatthat the constriction the constriction voltage voltage or alternatively or alternatively theand the field field notand the not the current current magnitude magnitude are responsible are responsible for triggering for triggering the RESET the RESET process. process. As another As another example example of the of the switching dynamics that can be achieved with the switching dynamics that can be achieved with the DMM, Figure 6e shows the case DMM, Figure 6e shows the case in in which the SB effect is not considered. Notice the hard threshold which the SB effect is not considered. Notice the hard threshold voltage for the SET con- voltage for the SET condition. dition. TheThe I-VI-V curve curve (red (red line) line) shows shows thethe intermediate intermediate current current states states (minorloops) (minor loops)forfora adamped dampedsinusoidal sinusoidalinput inputvoltage. voltage.These Thesestates statesare aregenerated generated by by the the hysteron hysteron (λ)(λ) shown shown in in green. green. Micromachines 2022,13, Micromachines2022, 11,330 x FOR PEER REVIEW 1110ofof2018 Figure 6. Effect of the initial memory state λ0 on the I-V characteristic. (a) Effect of the initial memory Figure 6. Effect of the initial memory state λ0 on the I-V characteristic. (a) Effect of the initial memory state λ on the I-V characteristic. (b) Application of a damped sinusoidal voltage on the memdiode state λ on the I-V (c) I-t characteristic, characteristic. (b)minor I-V curve with Application of a damped loops, and sinusoidal (d) I-VC curve. voltage VT is the on thevoltage. transition memdiode (e) I-t characteristic, (c) I-V curve with minor loops, and (d) I-V Simulation of intermediate states in the I-V characteristic without curve. V C considering T is the transition SB effect voltage. (e) Simulation of intermediate states in the I-V characteristic without considering SB effect. 3.2. Switching Dynamics 3.2. Switching Dynamics As the second criterion for assessing the model behavior, the DMM switching dy- As the second criterion for assessing the model behavior, the DMM switching dy- namics is discussed next. It is worth emphasizing that the switching properties of the namics is discussed next. It is worth emphasizing that the switching properties of the memory Equation (8) for constant and ramped voltage signals were reported in [46]. memory Equation (8) for constant and ramped voltage signals were reported in [46]. Briefly, Briefly, Equation (8) complies with the expected characteristic switching times (SET and Equation (8) complies with the expected characteristic switching times (SET and RESET) RESET) for a constant bias condition (Figure 7a): for a constant bias condition (Figure 7a): 𝜏𝑆,𝑅 (𝑉) = 𝜏0𝑆,𝑅 𝑒 ∓𝑉⁄𝑉0𝑆,𝑅 (13) τS,R (V ) = τ0S,R e∓V/V0S,R (13) and with the switching voltage value as a function of the applied signal ramp rate (RR) (Figure and with7a): the switching voltage value as a function of the applied signal ramp rate (RR) (Figure 7a):  𝜏𝑆,𝑅  𝑉𝑆,𝑅 = 𝑉𝑆,𝑅0 𝑙𝑛(𝑅𝑅) + 𝑉𝑆,𝑅0 𝑙𝑛 ( τS,R00) (14) VS,R = VS,R0 ln( RR) + VS,R0 ln 𝑉𝑆,𝑅0 (14) VS,R0 where 𝜏𝑆,𝑅0 and 𝑉𝑆,𝑅0 are the fitting constants. As shown in Figure 7a, as the constant ap- where τS,R 0 and VS,R0 are the fitting constants. As shown in Figure 7a, as the constant plied voltage increases, the maximum reachable current not only increases but also the applied voltage increases, the maximum reachable current not only increases but also SET switching time reduces in an exponential manner (see the inset in Figure 7a). This has the SET switching time reduces in an exponential manner (see the inset in Figure 7a). been experimentally demonstrated to occur in many material systems [61]. A similar This has been experimentally demonstrated to occur in many material systems [61]. A similar behavior is obtained for the RESET transition of the device. Although no close- Micromachines 2022, 11, x FOR PEER REVIEW 12 of 20 Micromachines Micromachines 2022, 2022, 11, 13, x330 FOR PEER REVIEW 12 11 of of 20 18 behavior is obtained for the RESET transition of the device. Although no close-form ex- behavior pression is obtained for form expression for afor a sinusoidal the RESET input sinusoidal transition is available, input of thethe the phenomenology is available, device. Although is similarno phenomenology isclose-form to that expressed similar ex- to that pression by expressed for Equation a sinusoidal by(14) input for a voltage Equation is available, (14) forramp. Figure a voltage the ramp. phenomenology 7b illustrates Figure 7b the is similar effects of illustrates thethe to that signal effects expressed offrequency the signal by on Equation frequency on(14) the I-V curve. thefor I-Vathe As voltage ramp. frequency curve. As Figure 7b the increases, the frequency illustrates SET and increases, theRESET the effects SET and of RESET the signal voltages shift frequency to higher voltages shift on the I-V values. to higherThis curve. effectAs values. is the This frequency consistent effect withincreases, withthe Equation is consistent (14) SET andand Equation has RESET (14)been voltages andexperimentally has shift to higher observed been experimentally values. using observedThis ramp effect rates using is consistent varying ramp ratesorders with varying Equation of orders magnitude (14) [62].and of magnitude has been Physically, [62]. experimentally the reason Physically, thebehind reasonobserved this behindbe- using havior ramp is therates this behavior varying voltage-time is the orders combined voltage-time of magnitude in[62]. actionaction combined Physically, theincharacteristic theswitching the characteristic reason behind this for be- timestimes switching the for havior is ionic/vacancythe the ionic/vacancy voltage-time hopping hopping combined represented represented action in by Equations the by Equations characteristic (9) and switching (10).(10). (9) and The The current times for increase current the increaseob- ionic/vacancy served observedin inthethehopping RESET RESET represented transition transition region by(VEquations region (9) and (V< <0)0)ofofFigure Figure 7b (10). 7bhas Thebeen hasalso also current been increase observed, observed, as ob- as it is served in the RESET a consequence of the transition incrementregion (V < 0) of of the current Figure lobe 7b haspoint triggering also been point [62]. observed, as it is [62]. a consequence of the increment of the current lobe triggering point [62]. Figure 7. (a) Figure 7. Effect of (a) Effect of aa constant constant voltage voltage input input signal signal on on the the I-t I-t characteristic. characteristic. The inset shows The inset shows the the time-to-switch Figure 7. (a) as a Effect function of a of the constant applied voltage voltage. input (b) signal Effect on the of I-t a sinusoidal signal characteristic. The with insetfrequency shows time-to-switch as a function of the applied voltage. (b) Effect of a sinusoidal signal with frequency theff on the I-t characteristic. time-to-switch as a function of the applied voltage. (b) Effect of a sinusoidal signal with frequency f on the I-t characteristic. on the I-t characteristic. Concerning Concerning the the switching switching dynamics dynamics for for discontinuous discontinuous signals, signals, Figure Figure 88 illustrates illustrates thethe Concerning effects sequencethe switching of equal dynamics amplitude for discontinuous voltage pulses (Vsignals, effects of a sequence of equal amplitude voltage pulses (Vapplied = 0.1, 0.3, 0.4, 0.5 of a applied = Figure 0.1, 0.3, 8 illustrates 0.4, 0.5 V) the and V) and effects period of period (T a== sequence (T 11 ss with of equal with duty duty cycleamplitude cycle == 0.5 0.5 s) onvoltage s) on the the currentpulses current (Vapplied = 0.1, magnitude. magnitude. As 0.3, 0.4,in As shown shown in0.5 V) and Figure Figure 8a, 8a, period for (T for aa device= 1 s with with an device with duty an initialcycle initial memory= 0.5 memory state s) on state λ the current λ00 ==00 (HRS), magnitude. (HRS), thethe current As shown current increases increases as in Figure as function function of8a, of for a device voltage voltage andwith and time.anThis time. initial This memory state corresponds corresponds to to the theλ0 so-called = 0 (HRS),potentiation so-called the current increases potentiation effect effect in as function of neuromorphic in neuromorphic voltage devices and devices [63].time. [63]. In This corresponds In addition, addition, for for higher to the so-called higher voltages, voltages, as as shown shown potentiation in in Figure Figure 8b, effect theincurrent 8b, the neuromorphic current not not only only devices [63]. Inincreases progressively addition,but for also higher voltages, switches to asLRS shown after in Figure the reaching 8b, the currentcondition threshold not only progressively dictated by theincreases SB effectbut also switches (pulse-induced to LRS after switching). Forreaching For negativethe negative threshold voltages voltages (see condition (see Figure Figure 8c), 8c), dictated by the SB effect (pulse-induced switching). For negative the current behaves in a similar fashion. First, the current decreases monotonically but as voltages (see Figure 8c), the sooncurrent as the behaves in a similar RESET condition fashion. is met, First, the the current currentandecreases exhibits monotonically abrupt reduction. In thisbut as latter latter soon case, as thethe RESET device condition memory stateisstarts met, the λcurrent at λ exhibits an abrupt reduction. In this latter 00 ==11(LRS). (LRS). 5 case, the device 10 memory state starts at λ0 = 1 (LRS). 7 0.5V a) 10 b) V 0.88 V applied 7 c) -0.63 5 0.4V 0.87 V 6 V -0.64 0.5V Vapplied 0.3V a) 8 b) V 0.88 V 0.86 V applied c) -0.63 -0.65 applied 4 0.4V V 0.87 V V 6 -0.64 0.1V 0.85 V applied -0.66 0.3V applied 8 0.86 V 5 -0.65 4 0.1V 0.85 V -0.66 [mA][mA] 5 [mA][mA] [mA][mA] 3 6 4 Current 6 Current 3 Current 4 3 Current Current 2 4 Current 3 2 4 2 1 2 2 1 1 2 1 0 0 0 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 0 0 0 0 20 40 Time [s] 60 80 100 Time [s] Time [s] 60 0 20 40 60 80 100 0 20 40 80 100 Time [s] Time [s] Time [s] Figure 8. (a) Effect of a pulsed signal on the SET I-t characteristic, (b) similar to (a) but with a higher voltage, Figure 8.and Figure 8. (a) (c) similar (a)Effect Effect of to (a) but ofaapulsed pulsed for on signal signal aon RESET the I-tI-t theSET SET curve. I-tcharacteristic, characteristic, (b) (b) similar similar to to (a) (a) but but with with aa higher higher voltage, and (c) similar to (a) but for a RESET I-t curve. voltage, and (c) similar to (a) but for a RESET I-t curve. Micromachines2022, Micromachines 2022,13, 11,330 x FOR PEER REVIEW 13of 12 of18 20 3.3. CRS 3.3. CRS Devices Devices Complementary Resistive Complementary Resistive Switching consists in the the anti-serial anti-serial combination combination of of two two memristive devices memristive devices [64,65]. [64,65]. This is the third criterion selected for evaluating the DMM. This is This is an an emblematic emblematic problemproblem to to demonstrate demonstrate the connection capacity of the model de- vices. CRSs are vices. are intended intended to tobebeused usedasasselector selectordevices devicesinincrossbar crossbararrays arrays with with thetheaimaimof reducing of reducing thethe crosstalk crosstalk effect [55,66]. effect [55,66].Different Different behaviors behaviors areareexperimentally experimentally observed observed de- pending onon depending thethe voltage voltage and andcurrent current window window investigated investigated including includingprogressive progressive andandab- rupt transitions abrupt transitions [67–71]. [67–71].Most Mostof the RS models of the RS models published published to date are unable to date are unable to cope with to cope with all these all these behaviors behaviors withinwithin a single a single framework. framework. As As shown shown in in Figure Figure 9a, 9a, the the top top device device (DMM1) (DMM1) isis initially initially inin HRS HRS (λ (λ00 == 0) 0) and and thethe bottom bottom device device (DMM2) (DMM2) in in LRS LRS (λ (λ00 == 1). 1). Figure 9 illustrates illustrates thethe combined combined actionaction of of both both memdiodes memdiodesin inthe thestationary stationaryloop. loop.The Thecurrent current behavior behavior is characterized is characterized by the by theappearance appear- of ancetwoof bumps two bumps (transmission (transmission windows) windows) at opposite at opposite voltages. voltages. The Thehigh highcurrent currentstatestate isis reached when both devices are in LRS. Remarkably, different reached when both devices are in LRS. Remarkably, different behaviors can be achieved behaviors can be achieved depending depending on on the the specific specific features features selected selected in in the the simulation simulation model. model. Figure Figure 9a 9a illustrates illustrates three three cases cases of of particular particular interest. interest. In In general, general, the the inclusion inclusion of of the the SB SB effect effect yields yields abrupt abrupt HRS/LRS transitions, while HRS/LRS transitions, while thethe absence absence of of the the SBSB effect effect leads leads toto progressive progressive transitions. transitions. 𝜆𝛾 ≠6=11ininEquation The γ The inclusion inclusion of of the the SF SF effect effect with with λ Equation(11) (11)results resultsininthe theappearance appearanceofofa alobe lobecurrent currentinin the LRS/HRS transition. In order to understand the LRS/HRS transition. In order to understand the complexity of the the complexity ofana- the analyzed problem, Figure 9b illustrates the potential drop distribution lyzed problem, Figure 9b illustrates the potential drop distribution across each device as across each device as a function a function ofof time time forforthe thecircuit circuitshown shownininthe theinset insetofofFigure Figure 9a 9a when when aa sinusoidal sinusoidal signalsignal is applied. The figure illustrates the cases with SB and is applied. The figure illustrates the cases with SB and without SF effects. without SF effects. Figure 9. Figure 9. (a) (a) I-V I-V characteristic characteristic for for aa CRS CRS structure structure w/o w/o SB SB and and SF. (b) Voltage SF. (b) Voltage drop drop distribution distribution as as aa function of time. function of time. 4. Experimental 4. Experimental Validation Validation The The model modeldiscussed discussedininthetheprevious previoussub-sections sub-sections was wasputput under under testtest by fitting experi- by fitting ex- mental datadata perimental extracted fromfrom extracted different published different works, published bothboth works, for Valence for Valence Change Memories Change Mem- (VCM) and Electrochemical ories (VCM) Memories and Electrochemical (ECM).(ECM). Memories In particular, Figure 10a,b In particular, Figureshows 10a,b the results shows the obtained for two for results obtained different VCM-type two different RRAM structures VCM-type with HfO RRAM structures 2 [72] with HfO and TaOand 2 [72] X [73] TaO asX the [73]dielectric layer. Figure as the dielectric 10c presents layer. Figure the fitting 10c presents results the fitting for afor results commercially a commercially available avail- ECM device able ECM [74] comprising device [74] comprisinga Wadopped Ge2Ge W dopped Se32Se active layer. 3 active layer.InInallallthe thecases theI-Vs casesthe I-Vs were were measured at room temperature and under voltage sweeps. The experimental data measured at room temperature and under voltage sweeps. The experimental data were werefitted fittedusing usingthe theSPICE SPICEmodel modeldescribed describedininAlgorithm Algorithm11 basedbased onon Equations Equations (6) (6) to to (11) (11) and and applying applying driving driving signals signals as as described described in in the the corresponding corresponding references. references. The The fitting fitting Micromachines 2022, 13, 330 13 of 18 Micromachines 2022, 11, x FOR PEER REVIEW 14 of 20 parameters are listed in Table 1 as a reference. It should be mentioned that the DMM does not onlyare parameters provide a simple listed in SPICE-compatible Table 1 as implementation a reference. It should be mentioned that forthe theDMM resistive does mem- orynot only provide devices but alsoa asimple SPICE-compatible versatile one, as it canimplementation accurately fit the for the I-V resistive memory loops experimentally devices but measured also a versatile in different RRAM one,devices, as it can while accurately fit the I-Vfor accounting loops experimentally their meas- particular features (as theured in different current RRAM limitation in devices, while Figure 10b accounting and for theireffect the snapback particular features in Figure (as the 10c). cur-that in Note rent limitation in Figure 10b and the snapback effect in Figure 10c). Note that in Figure Figure 10c, the x-axis accounts for the voltage effectively applied to the memristive device, 10c, the x-axis accounts for the voltage effectively applied to the memristive device, (V-IR) (V-IR) that results after subtracting the voltage drop in the series resistance (RS , which in that results after subtracting the voltage drop in the series resistance (RS, which in the the experiments is reported to be 46.25 kΩ [75]). experiments is reported to be 46.25 kΩ [75]). Figure 10. Experimental I-V loops of different memristor structures reported in the literature fitted Figure theExperimental with10. DMM model: (a)I-V loops2/Pt Ta/HfO of different memristor [72], (b) TaO X [73], andstructures (c) W doppedreported Ge2Se3 in the literature [74,75]. As refer- fitted ence, with thethe HRSmodel: DMM and LRS(a) curves are indicated Ta/HfO 2 /Pt [72],in (a). (b) Note TaO X that in [73], (b) and a current (c) W compliance dopped Ge 2 of3 1[74,75]. Se mA As was imposed reference, to prevent the HRS and LRSpermanent curves aredielectric breakdown, indicated in (a). which can be Note that in also (b) arepresented by the current compliance of DMM and SPICE. (d) Experimental and simulation results for the reset characteristics of SiO X from 1 mA was imposed to prevent permanent dielectric breakdown, which can be also represented UCL (data from [76]) using the QMM model. The inset shows the input signal. CRS: Experimental by and the DMM andI-V simulated SPICE. (d) Experimental characteristics. The arrows andindicate simulation results for the direction of the the reset characteristics applied bias: (e) of SiOX from UCL (data from [76]) using the QMM model. The inset shows the input signal. CRS: Experimental and simulated I-V characteristics. The arrows indicate the direction of the applied bias: (e) Pt/Ta2 O5 /Ta [69] and (f) Pt/Ta2 O4.7 /TaO1.67 /Pt [77]. (g) Pulse-enabled LTP and LTD data from [73] was also used to test the suitability of the DMM to replicate the potentiation and depression Micromachines 2022, 13, 330 14 of 18 behavior of memristor with multiple intermediate states. (h) LTP dependence on the pulses ampli- tude [78] can also be captured with the DMM. (i) The versatility of the DMM also allow fitting the current measured during the LTP of a Ag/ZnO/Pt nanowire memristive device [79]. Table 1. List of parameters used to fit the data presented in Figure 10. Imin Imax αmin αmax RSmin RSmax ηSET ηRESET VSET VRESET ISB Gam VT Work Material [A] [A] [a.u.] [a.u.] [Ω] [Ω] [V] [V] [A] [V] [72] Ta/HfO2 / 80 µ 1.1 m 2 2.75 100 150 8 10 600 m −575 m 300 µ 0 350 m Pt [73] TaOX 75 µ 1.5 m 2.4 4 120 120 40 7 375 m −130 m 1 0.05 350 m [74] W-Ge2 Se3 500 n 50 µ 4.3 1.75 10 10 50 250 200 m −20 m 700 n 0.35 50 m [76] SiOX 1µ 60 µ 3 3 1k 1 20 20 395 m −395 m 1 1 350 m Pt/Ta2 O5 / 3µ 0.9 m 3 1.75 160 160 50 50 2.4 −1.35 60 µ 0.3 0 [69] Ta 2µ 0.9 m 4 3 160 160 50 50 1.15 −1.05 40 µ 0.3 0 Pt/Ta2 O4.7 / 24.5 µ 200 µ 2 2 10 10 15 50 900 m −670 m 30 µ 2 600 m [77] TaO1.67 /Pt 17 µ 140 µ 2 2 100 100 100 50 750 m −820 m 50 µ 3 650 m [79] Ag/ZnO/ 450 p 3.5 n 2 2 200 200 2.4 10 1 −1 1 0 1 Pt A critical requirement for analog memristors is to account for intermediate states between the HRS and LRS regimes, which allow for a fine -tuning of the synaptic conduc- tance in neuromorphic circuits and the storage of more than one bit per device in memory applications. This is experimentally presented in Figure 10d for a RRAM stack comprising a SiOx dielectric layer [76], which is accurately modelled by the DMM model. In order to fully represent the intermediate states from LRS to HRS altogether with the major I-V loops, seven successive ramped voltage pulses with increasing amplitude were considered, as shown in the inset of Figure 10d. As previously discussed, the model also allows to represent CRS devices, which describes the combination of two bipolar memristive regimes with opposite polarities in a single device (as proposed by Linn et al. [64]). Figure 10e shows the experimental and simulation results for the I-V characteristic observed for the Pt/Ta2 O5 /Ta structures considered in [69]. Note the sharp opening of the ON-state window caused by the inclusion of the SB effect. For the particular case of symmetric systems (that is, those presenting two similar metal/oxide memristive structures), the representation of the device resistance as a function of the applied voltage results in the so-called “table with legs” shape [80]. In these cases, the two interfaces of the active switching layer behave in a complementary way [81]: when one switches from LRS to HRS, the other switches inversely. Such a shape can also be described by using memdiode devices, as shown in Figure 10f, for the case of the Pt/Ta2 O4.7 /TaO1.67 /Pt stack studied in [77]. Moreover, we have tested the capability of the DMM to replicate the Long-Term Potentiation (LTP) and Depression (LTD) of memristive synapses used in neuromorphic hardware. This kind of evolutionary behavior is required to achieve gradual conductance changes upon pulse applications. Figure 10g shows the LTP and LTD measurements reported in [73] for TiN/(25 nm)TaOX /Pt-based devices (the experimentally measured I-V loop is plotted in Figure 10b altogether with the corresponding fit with the DMM) by the application of 300 identical pulses of 1 V (LTP) followed by 300 identical pulses of –1.1 V (LTD). All pulses have the same width (100 ns, tON ) and they are applied every 20 msec. During the time in between pulses, a low voltage (0.1 V) pulse is applied to read the conductance (memory) state of the device. These measurement conditions were replicated in the SPICE simulations shown in Figure 10b. The simulated LTP and LTD trends are superimposed to the measurements shown in Figure 10g, showing a good fit of the experimental trends. There is also a voltage acceleration of the LTP/LTD trends, as reported for instance for the TiN/(3-nm)HfO2 /Pt-based devices measured in [78], which can be reproduced by the DMM, as shown in Figure 10h. Finally, the versatility and capability of the DMM to reproduce the evolutionary behavior of memristive devices is shown (see Micromachines 2022, 13, 330 15 of 18 Figure 10i) by also fitting the gradual current increase during LTP of Ag/ZnO/Pt-based nanowire memristors [79] by considering the same stimuli (pulses with an amplitude of 2.5 V and a duration of roughly 2 msec.). Thereby, the DMM is suitable to model the response of RRAM devices with a large number of incrementally accessible conductance states. For the sake of completeness, the fitting parameters for all the cases covered in Figure 10 are summarized in Table 1. 5. Conclusions A compact behavioral model for the I-V characteristic of bipolar resistive switching devices was presented. The Dynamical Memdiode Model relies on the combined action of two equations, one for the electron transport based on an extension of the quantum- point contact model and a second one for the internal memory effect that represents the ion/vacancy displacements. It was shown how the snapback and snapforward effects play a fundamental role in the SET and RESET processes, respectively. The model equations were implemented in the LTSpice simulator but can be easily translated to any other specific simulation language. The model can deal with arbitrary input signals, continuous or discontinuous, derivable or not. Fine tuning of the model equations could be required for specific situations. In summary, the proposed model is simple, robust, and accurate, as required for a fast and reliable simulation involving resistive switching devices. Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/mi13020330/s1, Supplementary Figure S1 “Simulation of the hysteron (top panel) and I-V characteristic (bottom panel) using a sinusoidal input “, Supplementary Figure S2 “Simulation of the forming I-V characteristic (first loop) and the stationary I-V curve (second loop) using a sinusoidal input “ and Supplementary Figure S3 “Simulation of the I-V characteristic with current compliance”. Author Contributions: Conceptualization, F.L.A., E.M. and J.S.; methodology, F.L.A. and E.M.; software, F.L.A. and E.M.; validation, F.L.A., E.M. and J.S.; formal analysis, E.M.; investigation, E.M.; resources, E.M. and J.S.; data curation, F.L.A. and E.M.; writing—original draft preparation, F.L.A., E.M. and J.S.; writing—review and editing, F.L.A., E.M. and J.S.; visualization, E.M.; supervision, E.M.; project administration, E.M.; funding acquisition, E.M. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by project 20FUN06 MEMQuD with funds from the EMPIR programme co-financed by the Participating States and from the European Union’s Horizon 2020 research and innovation programme. Data Availability Statement: The fitting data is available from the respective reference. The DMM model fitting each dataset is available in the following Bitbucket repository: https://bitbucket.org/ fernando_aguirre/uab_dmm/src/master/ (accessed on 19 January 2021). Acknowledgments: G. Patterson and A. Rodriguez are greatly acknowledged for their contributions to the development of the ideas reported in this work. Conflicts of Interest: The authors declare no conflict of interest. References 1. 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] [PubMed] 2. Banerjee, W.; Kim, S.H.; Lee, S.; Lee, S.; Lee, D.; Hwang, H. Deep Insight into Steep-Slope Threshold Switching with Record Selectivity (>4 × 1010) Controlled by Metal-Ion Movement through Vacancy-Induced-Percolation Path: Quantum-Level Control of Hybrid-Filament. Adv. Funct. Mater. 2021, 31, 1–9. [CrossRef] 3. Niraula, D.; Karpov, V. Comprehensive numerical modeling of filamentary RRAM devices including voltage ramp-rate and cycle-to-cycle variations. J. Appl. Phys. 2018, 124, 174502. [CrossRef] 4. Lin, J.; Liu, H.; Wang, S.; Zhang, S. Modeling and Simulation of Hafnium Oxide RRAM Based on Oxygen Vacancy Conduction. Crystals 2021, 11, 1462. [CrossRef] 5. Strukov, D.B.; Snider, G.S.; Stewart, D.R.; Williams, R.S. The missing memristor found. Nature 2008, 453, 80–83. [CrossRef] [PubMed] Micromachines 2022, 13, 330 16 of 18 6. 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] 7. Chua, L. Resistance switching memories are memristors. Appl. Phys. A 2011, 102, 765–783. [CrossRef] 8. Yakopcic, C.; Taha, T.M.; Subramanyam, G.; Pino, R.E. Generalized memristive device SPICE model and its application in circuit design. IEEE Trans. Comput. Des. Integr. Circuits Syst. 2013, 32, 1201–1214. [CrossRef] 9. Kvatinsky, S.; Friedman, E.G.; Kolodny, A.; Weiser, U.C. TEAM: Threshold adaptive memristor model. IEEE Trans. Circuits Syst. I Regul. Pap. 2013, 60, 211–221. [CrossRef] 10. Kvatinsky, S.; Ramadan, M.; Friedman, E.G.; Kolodny, A. VTEAM: A General Model for Voltage-Controlled Memristors. IEEE Trans. Circuits Syst. II Express Briefs 2015, 62, 786–790. [CrossRef] 11. Eshraghian, K.; Kavehei, O.; Cho, K.R.; Chappell, J.M.; Iqbal, A.; Al-Sarawi, S.F.; Abbott, D. Memristive device fundamentals and modeling: Applications to circuits and systems simulation. Proc. IEEE 2012, 100, 1991–2007. [CrossRef] 12. Biolek, D.; Biolek, Z.; Biolkova, V.; Kolka, Z. Modeling of TiO2 memristor: From analytic to numerical analyses. Semicond. Sci. Technol. 2014, 29, 2–7. [CrossRef] 13. Biolek, D.; Biolek, Z.; Biolkova, V.; Kolka, Z. Reliable modeling of ideal generic memristors via state-space transformation. Radioengineering 2015, 24, 393–407. [CrossRef] 14. Yang, J.J.; Strukov, D.B.; Stewart, D.R. Memristive devices for computing. Nat. Nanotechnol. 2012, 8, 13–24. [CrossRef] 15. 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] 16. 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] 17. Hajri, B.; Aziza, H.; Mansour, M.M.; Chehab, A. RRAM Device Models: A Comparative Analysis with Experimental Validation. IEEE Access 2019, 7, 168963–168980. [CrossRef] 18. Ielmini, D. Brain-inspired computing with resistive switching memory (RRAM): Devices, synapses and neural networks. Microelectron. Eng. 2018, 190, 44–53. [CrossRef] 19. Aguirre, F.L.; Pazos, S.M.; Palumbo, F.; Suñé, J.; Miranda, E. SPICE Simulation of RRAM-Based Crosspoint Arrays Using the Dynamic Memdiode Model. Front. Phys. 2021, 9, 548. [CrossRef] 20. Aguirre, F.L.; Gomez, N.M.; Pazos, S.M.; Palumbo, F.; Suñé, J.; Miranda, E. Minimization of the Line Resistance Impact on Memdiode-Based Simulations of Multilayer Perceptron Arrays Applied to Pattern Recognition. J. Low Power Electron. Appl. 2021, 11, 9. [CrossRef] 21. Aguirre, F.L.; Pazos, S.M.; Palumbo, F.; Suñé, J.; Miranda, E. Application of the Quasi-Static Memdiode Model in Cross-Point Arrays for Large Dataset Pattern Recognition. IEEE Access 2020, 8, 202174–202193. [CrossRef] 22. Aguirre, F.L.; Pazos, S.M.; Palumbo, F.; Antoni, M.; Suñé, J.; Miranda, E.A. Assessment and Improvement of the Pattern Recognition Performance of Memdiode-Based Cross-Point Arrays with Randomly Distributed Stuck-at-Faults. Electronics 2021, 10, 2427. [CrossRef] 23. Miranda, E. Compact Model for the Major and Minor Hysteretic I-V Loops in Nonlinear Memristive Devices. IEEE Trans. Nanotechnol. 2015, 14, 787–789. [CrossRef] 24. Patterson, G.A.; Sune, J.; Miranda, E. Voltage-Driven Hysteresis Model for Resistive Switching: SPICE Modeling and Circuit Applications. IEEE Trans. Comput. Des. Integr. Circuits Syst. 2017, 36, 2044–2051. [CrossRef] 25. Sune, J.; Miranda, E.; Nafria, M.; Aymerich, X. Point contact conduction at the oxide breakdown of MOS devices. In Proceedings of the Technical Digest—International Electron Devices Meeting; IEEE: New York, NY, USA, 1998; pp. 191–194. 26. Miranda, E.; Suñé, J. Analytic modeling of leakage current through multiple breakdown paths in SiO2 films. In Proceedings of the IEEE International Reliability Physics Symposium Proceedings; Institute of Electrical and Electronics Engineers Inc.: New York, NY, USA, 2001; Volume 2001, pp. 367–379. 27. Miranda, E.; Suñé, J. Electron transport through broken down ultra-thin SiO2 layers in MOS devices. Microelectron. Reliab. 2004, 44, 1–23. [CrossRef] 28. 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] 29. Datta, S. Electronic Transport in Mesoscopic Systems. In Cambridge Studies in Semiconductor Physics and Microelectronic Engineering, 1st ed.; Cambridge University Press: Cambridge, UK, 1997; ISBN 978-0521599436. 30. Miranda, E.; Jimenez, D.; Sune, J. The Quantum Point-Contact Memristor. IEEE Electron. Device Lett. 2012, 33, 1474–1476. [CrossRef] 31. Miranda, E.; Mehonic, A.; Suñé, J.; Kenyon, A.J. Multi-channel conduction in redox-based resistive switch modelled using quantum point contact theory. Appl. Phys. Lett. 2013, 103, 222904. [CrossRef] 32. Miranda, E.; Mehonic, A.; Blasco, J.; Sune, J.; Kenyon, A.J. Multiple diode-like conduction in resistive switching SiOx-based MIM devices. IEEE Trans. Nanotechnol. 2015, 14, 15–17. [CrossRef] 33. Zhong, X.; Rungger, I.; Zapol, P.; Heinonen, O. Oxygen-modulated quantum conductance for ultrathin HfO2 -based memristive switching devices. Phys. Rev. B 2016, 94, 165160. [CrossRef] 34. Guo, Y.; Robertson, J. Materials selection for oxide-based resistive random access memories. Appl. Phys. Lett. 2014, 105, 223516. [CrossRef] Micromachines 2022, 13, 330 17 of 18 35. Walczyk, C.; Walczyk, D.; Schroeder, T.; Bertaud, T.; Sowińska, 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. Devices 2011, 58, 3124–3131. [CrossRef] 36. Mehonic, A.; Vrajitoarea, A.; Cueff, S.; Hudziak, S.; Howe, H.; Labbé, C.; Rizk, R.; Pepper, M.; Kenyon, A.J. Quantum Conductance in Silicon Oxide Resistive Memory Devices. Sci. Rep. 2013, 3, 1–8. [CrossRef] [PubMed] 37. 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, 1–30. [CrossRef] 38. Blonkowski, S.; Cabout, T. Bipolar resistive switching from liquid helium to room temperature. J. Phys. D Appl. Phys. 2015, 48, 345101. [CrossRef] 39. Yi, W.; Savel’ev, S.E.; Medeiros-Ribeiro, G.; Miao, F.; Zhang, M.-X.; Yang, J.J.; Bratkovsky, A.M.; Williams, R.S. Quantized conductance coincides with state instability and excess noise in tantalum oxide memristors. Nat. Commun. 2016, 7, 11142. [CrossRef] 40. van Ruitenbeek, J.; Masis, M.M.; Miranda, E. Quantum Point Contact Conduction. Resist. Switch. 2016, 197–224. [CrossRef] 41. 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] 42. 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] 43. Zhao, J.; Zhou, Z.; Zhang, Y.; Wang, J.; Zhang, L.; Li, X.; Zhao, M.; Wang, H.; Pei, Y.; Zhao, Q.; et al. An electronic synapse memristor device with conductance linearity using quantized conduction for neuroinspired computing. J. Mater. Chem. C 2019, 7, 1298–1306. [CrossRef] 44. Karpov, V.; Niraula, D.; Karpov, I.; Kotlyar, R. A thermodynamic theory of filamentary resistive switching. Phys. Rev. Appl. 2017, 8, 024028. [CrossRef] 45. Dai, Y.; Zhao, Y.; Wang, J.; Xu, J.; Yang, F. First principle simulations on the effects of oxygen vacancy in HfO2-based RRAM. AIP Adv. 2015, 5, 017133. [CrossRef] 46. Miranda, E.; Sune, J. Memristive State Equation for Bipolar Resistive Switching Devices Based on a Dynamic Balance Model and Its Equivalent Circuit Representation. IEEE Trans. Nanotechnol. 2020, 19, 837–840. [CrossRef] 47. Waser, R.; Dittmann, R.; Staikov, C.; Szot, K. Redox-based resistive switching memories nanoionic mechanisms, prospects, and challenges. Adv. Mater. 2009, 21, 2632–2663. [CrossRef] 48. Wouters, D.J.; Menzel, S.; Rupp, J.A.J.; Hennen, T.; Waser, R. On the universality of the I–V switching characteristics in non-volatile and volatile resistive switching oxides. Faraday Discuss. 2019, 213, 183–196. [CrossRef] 49. Rodriguez-Fernandez, 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] 50. 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] 51. 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] 52. Linn, E.; Siemon, A.; Waser, R.; Menzel, S. Applicability of well-established memristive models for simulations of resistive switching devices. IEEE Trans. Circuits Syst. I Regul. Pap. 2014, 61, 2402–2410. [CrossRef] 53. Castán, H.; Dueñas, S.; García, H.; Ossorio, O.G.; Domínguez, L.A.; Sahelices, B.; Miranda, E.; González, M.B.; Campabadal, F. Analysis and control of the intermediate memory states of RRAM devices by means of admittance parameters. J. Appl. Phys. 2018, 124, 152101. [CrossRef] 54. Karpov, V.G.; Niraula, D.; Karpov, I.V.; Kotlyar, R. Thermodynamics of phase transitions and bipolar filamentary switching in resistive random-access memory. Phys. Rev. Appl. 2017, 8, 024028. [CrossRef] 55. Wouters, D.J.; Zhang, L.; Fantini, A.; Degraeve, R.; Goux, L.; Chen, Y.Y.; Govoreanu, B.; Kar, G.S.; Groeseneken, G.V.; Jurczak, M. Analysis of complementary RRAM switching. IEEE Electron. Device Lett. 2012, 33, 1186–1188. [CrossRef] 56. Mickel, P.R.; Hughart, D.; Lohn, A.J.; Gao, X.; Mamaluy, D.; Marinella, M.J. Power signatures of electric field and thermal switching regimes in memristive SET transitions. J. Phys. D Appl. Phys. 2016, 49, 245103. [CrossRef] 57. Lim, E.W.; Ismail, R. Conduction Mechanism of Valence Change Resistive Switching Memory: A Survey. Electron 2015, 4, 586–613. [CrossRef] 58. Celano, U.; Fantini, A.; Degraeve, R.; Jurczak, M.; Goux, L.; Vandervorst, W. Scalability of valence change memory: From devices to tip-induced filaments. AIP Adv. 2016, 6, 085009. [CrossRef] 59. Lv, H.B.; Yin, M.; Fu, X.F.; Song, Y.L.; Tang, L.; Zhou, P.; Zhao, C.H.; Tang, T.A.; Chen, B.A.; Lin, Y.Y. Resistive memory switching of Cuχ films for a nonvolatile memory application. IEEE Electron. Device Lett. 2008, 29, 309–311. [CrossRef] 60. Hardtdegen, A.; La Torre, C.; Cuppers, F.; Menzel, S.; Waser, R.; Hoffmann-Eifert, S. Improved switching stability and the effect of an internal series resistor in HfO2/TiOx Bilayer ReRAM Cells. IEEE Trans. Electron. Devices 2018, 65, 3229–3236. [CrossRef] 61. Cagli, C.; Nardi, F.; Ielmini, D. Modeling of set/reset operations in NiO-based resistive-switching memory devices. IEEE Trans. Electron. Devices 2009, 56, 1712–1720. [CrossRef] Micromachines 2022, 13, 330 18 of 18 62. Rodriguez-Fernandez, A.; Cagli, C.; Perniola, L.; Suñé, J.; Miranda, E. Effect of the voltage ramp rate on the set and reset voltages of ReRAM devices. Microelectron. Eng. 2017, 178, 61–65. [CrossRef] 63. Miranda, E.; Suñé, J. Memristors for Neuromorphic Circuits and Artificial Intelligence Applications. Materials 2020, 13, 938. [CrossRef] 64. Linn, E.; Rosezin, R.; Kügeler, C.; Waser, R. Complementary resistive switches for passive nanocrossbar memories. Nat. Mater. 2010, 9, 403–406. [CrossRef] [PubMed] 65. Burr, G.W.; Shenoy, R.S.; Virwani, K.; Narayanan, P.; Padilla, A.; Kurdi, B.; Hwang, H. Access devices for 3D crosspoint memorya). J. Vac. Sci. Technol. B Nanotechnol. Microelectron. Mater. Process. Meas. Phenom. 2014, 32, 040802. [CrossRef] 66. Aluguri, R.; Tseng, T.Y. Overview of selector devices for 3-D stackable cross point RRAM arrays. IEEE J. Electron. Devices Soc. 2016, 4, 294–306. [CrossRef] 67. Ambrogio, S.; Balatti, S.; Gilmer, D.C.; Ielmini, D. Analytical modeling of oxide-based bipolar resistive memories and complemen- tary resistive switches. IEEE Trans. Electron. Devices 2014, 61, 2378–2386. [CrossRef] 68. Yang, Y.; Mathew, J.; Shafik, R.A.; Pradhan, D.K. Verilog-a based effective complementary resistive switch model for simulations and analysis. IEEE Embed. Syst. Lett. 2014, 6, 12–15. [CrossRef] 69. Miranda, E.A.; Frohlich, K. Compact Modeling of Complementary Resistive Switching Devices Using Memdiodes. IEEE Trans. Electron. Devices 2019, 66, 2831–2836. [CrossRef] 70. La Torre, C.; Zurhelle, A.F.; Breuer, T.; Waser, R.; Menzel, S. Compact modeling of complementary switching in oxide-based ReRAM devices. IEEE Trans. Electron. Devices 2019, 66, 1265–1275. [CrossRef] 71. Petzold, S.; Miranda, E.; Sharath, S.U.; Muñoz-Gorriz, J.; Vogel, T.; Piros, E.; Kaiser, N.; Eilhardt, R.; Zintler, A.; Molina-Luna, L.; et al. Analysis and simulation of the multiple resistive switching modes occurring in HfO x -based resistive random access memories using memdiodes. J. Appl. Phys. 2019, 125, 234503. [CrossRef] 72. Jiang, H.; Han, L.; Lin, P.; Wang, Z.; Jang, M.H.; Wu, Q.; Barnell, M.; Yang, J.J.; Xin, H.L.; Xia, Q. Sub-10 nm Ta Channel Responsible for Superior Performance of a HfO2 Memristor. Sci. Rep. 2016, 6, 1–8. [CrossRef] 73. Wang, Z.; Yin, M.; Zhang, T.; Cai, Y.; Wang, Y.; Yang, Y.; Huang, R. Engineering incremental resistive switching in TaO: X based memristors for brain-inspired computing. Nanoscale 2016, 8, 14015–14022. [CrossRef] 74. Knowm SDC Memristors. Available online: https://knowm.org/downloads/Knowm_Memristors.pdf (accessed on 18 Jan- uary 2022). 75. Ostrovskii, V.; Fedoseev, P.; Bobrova, Y.; Butusov, D. Structural and Parametric Identification of Knowm Memristors. Nanomaterials 2021, 12, 63. [CrossRef] [PubMed] 76. Mehonic, A.; Shluger, A.L.; Gao, D.; Valov, I.; Miranda, E.; Ielmini, D.; Bricalli, A.; Ambrosi, E.; Li, C.; Yang, J.J.; et al. Silicon Oxide (SiOx): A Promising Material for Resistance Switching? Adv. Mater. 2018, 30, 1–21. [CrossRef] [PubMed] 77. Ferreyra, C.; Sánchez, M.J.; Aguirre, M.; Acha, C.; Bengió, S.; Lecourt, J.; Lüders, U.; Rubi, D. Selective activation of memristive interfaces in TaOx-based devices by controlling oxygen vacancies dynamics at the nanoscale. Nanotechnology 2020, 31, 155204. [CrossRef] [PubMed] 78. Matveyev, Y.; Kirtaev, R.; Fetisova, A.; Zakharchenko, S.; Negrov, D.; Zenkevich, A. Crossbar Nanoscale HfO2-Based Electronic Synapses. Nanoscale Res. Lett. 2016, 11, 147. [CrossRef] [PubMed] 79. Miranda, E.; Milano, G.; Ricciardi, C. Modeling of Short-Term Synaptic Plasticity Effects in ZnO Nanowire-Based Memristors Using a Potentiation-Depression Rate Balance Equation. IEEE Trans. Nanotechnol. 2020, 19, 609–612. [CrossRef] 80. Chen, X.; Wu, N.J.; Strozier, J.; Ignatiev, A. Direct resistance profile for an electrical pulse induced resistance change device. Appl. Phys. Lett. 2005, 87, 233506. [CrossRef] 81. Nardi, F.; Balatti, S.; Larentis, S.; Ielmini, D. Complementary switching in metal oxides: Toward diode-less crossbar RRAMs. In 2011 International Electron Devices Meeting; IEEE: New York, NY, USA, 2011; pp. 709–712. [CrossRef]