Un estudio colectivo sobre modelado y simulación de memoria resistiva de acceso aleatorio

The comparative analysis of the VTEAM model with the Yakopcic model [73, 74], BCM model [99] (discussed further in this article), and the TEAM model are presented in Fig. 9 [77]. It represents the flexibility that the model possesses, as it can be tuned to fit all the three models. It shows good agreement with all the three models illustrated, respectively, in Fig. 9a–c [77]. Fundamentally, the TEAM/VTEAM models are quite generalized physics-based models. This means that with the help of fitting parameters, they can be comparable with the multitude of other models, and fit to a variety of experimental characterization data from memristive systems.

The VTEAM model is compared with previously proposed memristor models [77]. un Yakopcic model [73]. b BCM model [99]. c TEAM model [76]

Stanford/ASU Model

A physics-based model which has become very popular is the one developed by Guan et al. and Chen et al. of Stanford University and ASU, known as Stanford/ASU model [78,79,80]. This model is exclusively developed for RRAM devices, rather than a generalized one for memristive systems which was fitted for those particular devices. It included the effect of critical phenomenon of switching such as Joule heating and temperature change, which had been neglected before. The developed model was applied in the I - V switching characteristics of HfO₂ RRAM [144]. Along with it, Verilog-A [79] and SPICE [81] implementations of the model are also presented.

This model is based on the growth of conductive filament. The CF growth leaves a gap with the top electrode which is called as the filament gap. This growth of the filament gap is considered as internal state variable in this case. So, the rate of filament growth and the filament gap govern the dynamics of the model. The filament growth is explained due to the movement of oxygen ions and vacancy regeneration and recombination [145]. Considering the gap value g (nominally in the range of 0–3 nm) to be the state variable, the rate of change of g is defined as [78]:

The parameter E _a is the activation energy for vacancy generation and oxygen vacancy migration in the SET and RESET processes, respectively. v is the applied voltage across the device, ν ₀ the velocity containing the attempt-to-escape frequency, L the switching material thickness and a _h , the hopping site distance.

A significant feature of this model is the inclusion of variations in the model caused due to the stochastic property of the ion process and the spatial variation in the gap size among multiple filaments. To account for these variations in the model, a noise signal is added to the gap distance as [78]:

The variation in the gap size δ _g is defined as a function of the ions’ kinetic energy and invariably on the temperature in the filament and is given as [78]:

Here, T _crit is defined as a threshold temperature beyond which there is a significant change in the gap size. This can be understood as the point where the device undergoes a physical transformation such as transitioning into a SET or RESET state. In this case, threshold is considered in terms of temperature, rather than voltage or current, whatever employed in the previous models [75,76,77]. So, the equation basically depicts the resistance fluctuation that occurs when the CF temperature is increased beyond the room temperature.

Now that temperature can be considered a critical driving force in the model, a modified form of the steady-state Fourier heat flow equation is implemented in this model. Rather than considering heat flow throughout the filament, the vicinity of the tip of the filament is considered. There is a dynamic inner domain temperature T which significantly changes with change in the cell characteristics, and an outer domain remains at an ambient room temperature T _amb , related as [78]:

c _p is the effective heat capacitance of the inner domain, and k the effective thermal conductivity are both fitted based on the type of oxide and electrodes used in the RRAM system. RESET transition from LRS to HRS generally has higher temperature associated with it across the device, while the SET transition has a considerably lower temperature. The current inside the device is modeled using a generalized conduction mechanism where the tunneling distance and field strength have an exponential relationship. This is true in case of tunneling current conduction mechanisms such as Poole-Frenkel, Fowler-Nordhiem, trap-assisted, or direct tunneling [9, 16, 46, 49, 51, 55]; these are the mechanisms most commonly associated with RRAM systems [51, 55, 61, 66]. The current conduction is defined as [78].

The advantage with a generalized current equation is that for a particular device if some other mechanism is fitting better, it can be incorporated easily by adding the required parameters and adjusting their values accordingly. yo - V response of the model compared with experimental data is shown in Fig. 10. The experimental response is shown in Fig. 10a while the simulated curve is shown in Fig. 10b. Simulated transient response shows the capabilities of the model in taking variations into account. Developed model was verified using Ngspice [146] as a macrocircuit. Ngspice is an open source SPICE simulator which is quite efficient and convenient for doing DC and AC analysis. This model can be implemented in MLC memory circuits and also to verify the efficiency of programming strategies and error correction codes [78].

un Experimental and b simulated transient responses of a HfO _x RRAM device to the −  2.3 V 50 ns input pulses. The experimental result is reported elsewhere [144] and included here in a for convenience. c In a larger time range, the simulated transient response for the same device including the gap size and temperature is shown. Current compliance set at 200 μ A in simulation [78]

A major feature of this model is implemented in the neuromorphic systems and RRAM synaptic device design [147]. This model has been tested against a HfO _x /TiO _x multi-stack RRAM system [148] which is implemented in a neuromorphic system. This gives the model great flexibility and wide applications as there are only a few models that are actually applicable for neuromorphic systems. Also, the model defined for these systems has been deemed tolerant to training error caused by device variation [149]. The gradual resistance modulation which is critical to the learning process in a synaptic device can be quantified in the model [150] which marks a significant development in using RRAM synaptic stacks in neuromorphic computing systems.

Physical Electro-Thermal Model

This model is an extensive physical model which describes the bipolar operation in RRAM devices using equations closely resembling the physical mechanisms. This model was reported by Kim et al. [87], and it was verified with a tantalum pentoxide (Ta₂ O₅ )-based bi-layered RRAM structure [15, 151, 152]. It makes use of the finite element solving method employed in the previous model to solve the differential equations. The major value addition by this model over the model proposed by Larentis et al. [86] was the proper description provided for the SET state in the bipolar RRAM device. The previous model was inadequate in accommodating the complete transition and explaining it properly but this model makes up for that. Also, it improved upon a physical electro-thermal model reported by Menzel et al. [153] which attempts at calculating the CF temperature precisely.

It also uses the electro-thermal physics phenomenon approach for modeling which we have seen in the previous model [86]. The major advantage with models based on this concept is their ease of use owing to the simple fundamental equations and the flexibility to employ a proper finite element method (FEM) solver to simulate the system very accurately. But a major disadvantage is that the model becomes very difficult to implement in circuit solvers based on SPICE and providing an equivalent implementation in Verilog. This is because of the lack of support in SPICE and Verilog for properly defining partial differential equations which make up for the vastness of the model. Normal ordinary differential equations and the ones which are in analytical form can be solved in circuit solvers but partial differential equations (PDE) cannot be solved.

Electro-thermal models are equally important as compared to the other physics-based models discussed before because temperature is an important factor governing the set and reset processes. Ion and vacancy migration plays a dominant role for switching mechanism [16, 46], although the governing factors are behind this process and the exact type of ions is still up for debate. So, the fact that temperature is a governing factor in this process makes these models attention worthy. Also, experiments [85, 154] in this regard suggest that there is significant change in the temperature in the CF during the switching process. Some of the previous models discussed above have neglected this effect by considering conducting filament-oxide interface to be at room temperature or by taking constant conducting filament temperature [39, 86, 88, 89, 144].

The major difference between this model and the previously discussed electro-thermal model is in the expressions used to describe the drift-diffusion process. CF is described as a doped region where the oxygen vacancies act as dopants, and the CF runs from the top electrode to the bottom electrode. This is an assumption that many models take that the CF runs from one end of the electrode to the other when the state variable is considered as the length of CF. A few models discussed previously [78, 80] have used the filament gap to the top electrode as state variable. So, the assumptions generally vary from system to system and are dependent on what mechanism is employed to describe the device.

Another assumption taken to describe the drift-diffusion of vacancy migration is that the same equation used can describe both the oxygen ions and vacancies. This is generally the case to simplify the model and reduce the complexity of the equations. The rigid point ion model by Mott and Gurney [155] is employed here to describe the process given as [87].

donde D _s describes the diffusion process, v gives the drift velocity of the vacancies, and G is the generation rate of vacancy or the CF growth rate which actually describes the SET process. The G term is a specialized parameter added to better describe the complete switching process [156, 157]. The parameters are defined as [87]:

Here, l _m is the mesh size. So, using the Eqs. (48)–(50), the oxygen vacancy transport given in Eq. (47) can be defined which contains all the factors of drift-diffusion as well as the vacancy regeneration. These equations govern the CF growth and rupture which defines the physical transformation of the device during the SET and RESET transition of the device. So, it basically acts as a dynamic internal state variable which controls the switching rate of the device.

The simulation results for the reset transition is shown in Fig. 11 [87]. Concentration of the oxygen ions is shown at different voltages in Fig. 11a [87] which invariably governs the switching in the device. The point C (3.0 V) is the point where the reset transition occurs, so the concentration of ions is also the highest at the interfaces for that voltage point as evident in Fig. 11b [87]. On similar lines, the temperature and flux are on the higher side which can be seen in Fig. 11c, d, respectively [87].

Simulation results for the reset transition of the device. un V _o density (n _D ) map. Calculated profiles of b n _D , c T y d y for states A (1.0 V), B (1.7 V), and C (3.0 V). The position of z  = 15 nm indicates the Ta₂ O₅ /TaO _x interface in the structure schematic. The shaded area shows the depleted gap, defined for n _D  < 5 × 10 ²¹ cm ^-3 [87]

Equations (95) and (98)mentioned further are also used in the model to describe the current conduction and the temperature change due to Joule heating in the device. The equations are simultaneously solved in COMSOL to generate the required simulated profiles. The obtained simulated profiles are compared and verified against a TaO_x bi-layered RRAM system [87]. In addition to the DC I - V characteristics the model was also used to generate time-dependent reset characteristics by investigating its response to square pulses.

Huang’s Physical Model

A very comprehensive physical model of RRAM devices is developed by Huang et al. [88, 89]. Its major feature is its consideration of the multitude of factors affecting the CF dynamics in the RRAM device. This model is comprehensive in the sense that it considers both the width of CF as well as filament gap to the electrode as factors affecting the state variable dynamics. The model was validated in a TiO₂ based device and also applied in a 2 × 2 RRAM array cell [88].

Covering bipolar devices primarily, it also accounts for the temperature distribution in the device with multiple heating sources. SET/RESET process is considered to be caused due to generation/recombination process of the oxygen ions (O ²⁻ ) and oxygen vacancies (V _o ). Top electrode (TE) is the active electrode and acts as an oxygen reservoir for the release or absorption of oxygen ions [88]. The CF evolution during the SET process is modeled based on the width of the CF. Growth of the CF is thought to start from the tip of the active electrode. With an increase of voltage the CF enlarges along the radius resulting in a final width of the CF as w . So, the value of w is critical to determine the LRS resistance in the SET process. Huang et al. [88] assumed that the CF grows in a symmetrical cylindrical shape which is simplifying at best. While the cylinder has been the most popular to describe the shape of the CF, it might not be the most accurate.

Rupture of the CF during the reset process is considered to start from the TE first. CF disconnects from the starting point and then dissolves internally with increase in the voltage. Distance between the tip of the CF and the active electrode layer is defined as the filament gap distance (x ). The value of x determines the resistance of HRS during the RESET process. x and dx/dt are thus critical in defining the RESET process. A very important feature of the model is that there are two parameters defining the state of the system, in place of one parameter. The parameter w acts as the state variable for the SET process and x for the RESET process. So dx/dt and dw/dt define the dynamics of the device during the SET/RESET transition. Analytical model for a RRAM cell presented by Huang et al. [88] is developed by modeling the parameters x, w and their evolving speeds.

This model also presents one of the most detailed descriptions for the processes involved behind the RESET process. The rate of the CF shortening is affected by three processes, (a) O ²⁻ release by the electrode, (b) O ²⁻ hopping in the oxide layer, and (c) recombination between O ²⁻ y V _o . Slowest process among the three dominates the CF reduction process which is defined by the parameter x . Speed of the processes is affected by the specific device characteristics and the oxide used.

CF reduction rate during first reset process, i.e., O ²⁻ release by the electrode can be given as [89]:

In case of the O ²⁻ hopping in the oxide layer, the CF with a being the distance between two V_o , reduction rate is described by [89]:

The RESET process when dominated by the recombination between O ²⁻ and V_o is written as [89]:

The value of x is fixed to x ₀ after the RESET process. This invariably will act as the boundary condition for the model. But the problem here is the value and the role of x ₀ is not clearly defined here. This will possibly create ambiguities while defining the states of the device or switching between two states. In the first step of the SET process which is dominated by recombination of oxygen vacancies and where a thin CF is initially grown is described by [89]:

Here, Z and α_a are fitting parameters. In the second step, the CF grows along the radial direction of the CF is defined as [89]:

Current flowing through the device has been taken in the model due to the hopping conduction and metallic conduction. The current in CF region can be calculated using the basic structures of Ohm’s law and Arrhenius law [158]. But the current in the gap region as a result of hopping conduction is given a little different. It is modeled as a correlation of the hopping current with the voltage and gap distance is given by [147]:

Temperature effects in the model are considered from the Filament Dissolution model [82, 83] discussed further in the “Filament Dissolution Model” section. Validation of the model is performed in HfO_x /TiO_x system [88, 89]. Transient results obtained from simulating the model are compared against the data from the device, which shows a good match as demonstrated by Huang et al. [88]. The model is also validated against devices fabricated by other groups [144, 159] and the parameters are adjusted accordingly. A pretty accurate match between the simulation and the experimental results suggests a good level of flexibility with the model. The model also demonstrates that the switching speed of the device is highly dependent on the input voltage sweep rate.

Although the model is very comprehensive and takes into account a variety of detailed processes affecting the RRAM operation; it has some critical shortcomings. A major one is the non-compatibility with SPICE or Verilog-A. Implementations in any of the circuit simulators based on these platforms has not been demonstrated which raises a question on its readiness for simulations. Also, boundary conditions and non-linear effects have not been applied in the model which leaves it open to unphysical solutions. There has been no attempt to fit a window function with the model to account for this effect. These shortcomings make the model difficult for application for simulations, but its physics give a lot of insights into the functioning of RRAM devices.

Bocquet Bipolar Model

A very interesting and unique model from Bocquet et al. [90, 92] which utilizes a physics based modeling approach to describe bipolar oxide based resistive switching memories. This was a model developed exclusively for the RRAM devices. Although a point of speculation still exists, it has been more or less accepted that the bipolar resistive switching mechanism is governed by the valence change mechanism which occurs in specific transition metal oxides and the field-assisted motion of oxygen ions O ²⁻ [160].

This is also one of the few models that can describe electroforming process. This process basically initiates the CF growth for the first time when the device is in a pristine state. It requires significantly higher voltage as compared to the set or reset voltage because the CF formation requires an electric breakdown of the oxide and this requires higher voltage and energy. However, forming free RRAM devices have been reported [85] by adjusting the oxygen stoichiometry of the active layer. Removal of the forming process will reduce the voltage requirement of the device and make it more energy efficient.

Bocquet bipolar model uses some concepts from the Bocquet unipolar model [90] and modifies it significantly according to the bipolar switching characteristics. Major features of the model are its intrinsic simplicity in the model equations, full compatibility with SPICE based electric simulators and inclusion of voltage and time dependencies of the device. Internal state variable here is the radius of the CF which governs the switching rate. Radius of the CF varies with growth/rupture mechanism of the CF which is explained in the model with the help of local electrochemical redox processes [82, 83, 105, 161] which are dependent on the applied bias polarity. A single master equation in which both the SET and RESET processes are accounted for simultaneously is controlled by the CF radius which thus gives the switching rate of the device.

Electroforming stage is modeled using electroforming rate which describes the process of conversion of the pristine oxide into a switchable sub-oxide layer. CF radius (r _CF ) varies from a minimum value of 0 to a maximum value of r _CFmax . The electroforming stage is modeled as [92]:

Some of the simplifying assumptions in the model are regarding the current conduction in the LRS and HRS. During the LRS, the conduction is assumed to be Ohmic, i.e., it follows Ohm’s law. In the HRS region, the current is dominated by a leakage current in the sub-oxide region which is basically due to trap-assisted conduction, but for simplicity sake, Ohmic conduction is considered here. The SET/RESET operation in the model is described by the electrochemical redox reaction derived from the Butler-Volmer equation [162] given as [92]:

Here, τ_Red and τ_Ox are the reduction and oxidation reaction rates, respectively. τ_Redox is the effective reaction rate considering both the reduction and oxidation reactions. Above two equations are coupled together in a master equation which define the switching rate given as [92]:

This is quite a comprehensive model in the sense that it includes the temperature effects as well. Temperature plays a significant role in the redox reaction rates [163, 164] and thus the local temperature in the filament is a very important parameter in this regard. The basic heat equation is used in this model and modified it accordingly given as [92]:

On the face of it, the equations seem pretty complex to evaluate. But in reality, they are analytical in nature which makes them easily solvable in a numeric solver and can be implemented in an electric simulator. This is a major advantage of this model. Almost all of the models which employ the concept of temperature change in the filament follow the basic principles of the filament dissolution model [82, 83] discussed further in the “Filament Dissolution Model” section. During set operation, the temperature rises due to the increase in the CF radius, while it falls due to a decrease in the CF radius during the reset operation. This creates a positive feedback loop between the two processes leading to a self-accelerated reaction. This forms the basis of the filament dissolution model and all models incorporating the temperature effects in the device converge on this phenomenon [82, 83, 86,87,88,89, 92].

yo - V characteristics of NiO based RRAM along with simulated curve using Bocquet model is presented in Fig. 12 [91]. Figure 12a represents the set and reset transitions of the device while Fig. 12b highlights the forming process. The current conduction in the Bocquet bipolar model is treated a little differently from what we have seen from previous models [87, 88, 90]. It considers the current as a combination of contributions from three different sources. The first one is the current from the conductive area (i _CF ), the second is the conduction through the switchable sub-oxide (i _sub-oxide ) and then the conduction through the pristine device (i _pristine ). The total current is described as [92]:

un , b Experimental I (V) characteristics for Electroforming , Set , and Reset processes measured on a large number of memory elements to understand the device-to-device variability. The experimental device-to-device variability is accounted for in Monte Carlo simulations with a ± 5% standard deviation on parameters α and L _x [21]

The parameter B_e is metal-oxide barrier height (ϕ _b ), dependent which is given as [92]:

Here, S_cell is a section of the RRAM cell. A _e y B _e are additional parameters defined to make the equations concise. To implement the model in an electrical simulator, discrete solutions are required which are well provided in this model. This makes the model suitable for proper simulation involving electrical circuits, thus widening its use case scenario. In this model, the equations are implemented in an Eldo circuit simulator [165, 166]. Memory effect of the device was replicated in the form that for each cell of the RRAM instance during transient simulation, the previous state of the filament as well as the applied voltage are given as the present state of the device [92]. New state gets solved as a function of these new inputs and the time step considers in the transient simulation. The discrete solutions are given as [92]:

The model has been verified against electrical characterization from an HfO₂ based system [167]. To better judge the model on circuit level, a 2T/1R bipolar OxRAM [168] cell was simulated using Eldo, as shown in Fig. 13a [92]. Simulation of this type helps check the stability of the model when applied to a system environment. Current variation, voltage, and CF radius (r _CF ) are shown with respect to time. Voltage follows a triangular wave form, which is the input sweep. Current in the device transitions from high to low and vice-versa depend on the voltage levels. The sudden drop in the current levels, as shown in Fig. 13b [92], indicates the device transition. CF radius follows a similar path as the current which is expected behavior of the internal state variable.

Simulating the electrical characteristics of the considered 2T-1R OxRAM structure [174]

Another important feature of OxRAM that has been highlighted in the model is the soft-reset [168]. It mainly induces the dependence between resistance in HRS and the stop voltage during the preceding reset operation. This phenomenon is basically due to the incomplete destruction of the CF during the reset process. So, the CF radius and temperature decrease during this process, leading to a decrease in the reaction rate. This means a self-limited reaction rate thus getting the name soft-Reset . This model can account for the device to device variability very efficiently [169, 170]. The standard deviation obtained for the important parameters such as the length of the oxide (L_x ) is well within the accepted range, thus accounting for the variations when the materials change in different devices [167, 168].

A shortcoming in the model which can be highlighted is the lack of a voltage or current threshold. Also, it works on the simplifying assumption that the CF radius grows from one end of the top electrode to the other end of the bottom electrode. This makes the model immune to significant fluctuations if the growth of the CF is not complete, thus leaving a filament gap. There is no provision to account for the effect of the filament gap if it occurs.

Berco-Tseng Model

The proposed model and simulation approach [171,172,173,174,175] by Berco-Tseng for RRAM devices is based on describing the CF growth process. The Gibbs free energy criteria [174, 175] is used as an indicator to model the growth dynamics of the CF. Simulation approach for the forming, set and reset process in the model is based on the Metropolis Monte Carlo algorithm [174]. This approach importantly does not rely on time evolution of the CF, thus increasing the efficiency of comparison of the relative retention properties of MIM structures.

The model is quite comprehensive in terms of describing the underlying physical parameters which affect the CF kinetics in the resistive switching layer. It also introduces the concept of “hot-spots” [172,173,174] which are random localized initial clustering of oxygen vacancies which facilitate the formation of the CF. The major parameter governing the Gibbs free energy is the enthalpy of formation of an oxygen vacancy [174] is used to define the CF growth dynamics in the switching layer and integrate it into the Monte Carlo simulator. As a result, all the CF processes, namely forming, set and reset can be effectively simulated.

A monoclinic HfO₂ switching layer is primarily used to implement the model. But other material stacks, such as ZrO₂ , Cu-HfO₂ are also studied and compared with by us. Electrical conduction in the RRAM device is modeled on the Arrhenius relation, given as [174]:

Typical boundary condition, such as V _dd at top electrode and ground at bottom electrode is applied to the device. For modeling the CF accurately, it is divided into a grid structure to discretize it, which is in line with the finite element analysis (FEA) method. The various parameters defining each grid site are its spatial coordinates (x , y ), local potential φ , temperature T , N _o , N _ov , trap occupancy c , electrical conductivity σ and thermal conductivity k _t h . The various processes associated with the evolution of CF within the oxide layer involves generation, recombination and hopping of oxygen (O ), oxygen vacancies (OV ) and electrons.

As a result, these processes are defined in terms of probabilities in the MC simulator. The probabilities are defines according to the minimum energy criteria as discussed earlier. This approach smoothens out the iterative steps in the simulator rather than the abrupt + 1 or − 1 levels in the discrete approach. The generation (P _g ), hopping (P _h ), and recombination (P _r ) probabilities for the oxygen species are given as [174]:

Here, E is the local electric field, C _n represents the ratio of N _ov (density of oxygen vacancies) in the low state to the maximal one at site n (n ^th grid site). E _a y E _h are the activation energies for oxygen species generation and hopping respectively. Similarly, a _a and a _h are the field lowering factor for O generation and hopping.

The mechanism of the current conduction is considered to be trap-assisted tunneling [171, 173, 174] through the switching layer. The Mott variable hopping model [174, 175] is taken into account here to model the tunneling effect. Mott DC conduction considered at any two grid points m, n are given as [174]:

Here, d _mn is the distance between m and n , α is the typical attenuation length of the electron wave function in the trap and c being the trap occupancy.

The model is simulated and the results are compared with many different material stacks of RRAM devices from other groups as well. The Fig. 14 shows the simulated heat map of the N _o (density of oxygen) concentration during the CF rupture process. The simulation result corroborates with the experimental results, showing wide applicability potential for the model.

Resistive switching layer and Ti/HfO₂ interface oxygen spatial density (No) plot during the CF rupture using negative bias:a initial and b final [174]

Gonzalez-Cordero et al. Bipolar Model

It is a compact physical model proposed by Gonzalez-Cordero et al. [93] describing the working of bipolar RRAM systems. The model is unique because it considers the CF as a truncated cone, which is a significant departure from previous models considering the CF shape generally as a cylinder. Also, the model is validated by implementing it in Verilog-A which gives us a closer look into the description and simulation of RRAM devices on the circuit level using Verilog-A. The proposed model builds upon the concepts introduced in the previous Bocquet bipolar model [91, 92] and modifies it accordingly to suit the new CF shape proposed.

One of the important aspects about the model is the consideration of a truncated cone shaped CF [176,177,178,179]. Majority of the models we have encountered till now consider the CF as a symmetrical cylinder which is more of a simplifying assumption [91, 92]. This is because it has been shown that the CF can grow [39, 51, 55] from one end of either electrode to the other depending on the active electrode. So, it is quite possible that the CF in this case might not be a perfect cylinder. So, a truncated cone is equipped to account for any variability and fluctuations arising due to the shape of the CF. Shapes other than simplified geometrical shapes are not considered in the models because of algebraic complexities. In previous models [91, 92], we noticed that device to device and cycle to cycle variability’s have a significant effect on the application of particular models to devices. So, by taking a truncated cone as the CF shape provides this model more flexibility than the others.

Another significant feature is the role of temperature in the CF and the reset process. Majority of the models which describe the CF rupture due to the self-accelerated dissolution, consider that the process takes place at the CF narrowing point and temperature increases at that particular point [82,83,84, 91, 92]. This point is generally in the middle of the cylindrical CF due to its symmetry. So, when we look at it from a physical standing point, the temperature at each of the points in all the RRAMs stacked together in a circuit has to be evaluated. Realizing this from the circuital standing point and simulating thousands of devices in the circuits is a very time-consuming process and slows down the simulation. This problem can be circumvented by considering two temperatures in the CF instead of the general single temperature; this approach also keeps the simplicity of the model intact. Two temperatures represent the main CF body that is not destroyed during the reset operation and the CF narrowing. This has been implemented in this model by considering the two temperatures as the wide region and the narrow region of the truncated cone, respectively.

This model extends the previously discussed Bocquet bipolar model [91, 92] in the “Bocquet Bipolar Model” section. In the previous model, the equations were defined keeping a cylindrical CF in mind, so the equations here have been modified to account for the change in the CF shape. The truncated cone CF is described by two different radii. CF is considered to grow from the top electrode to the bottom; the interface radius with the top electrode (TE) is r _CFT which is always greater than the radius of the interface with the bottom electrode (BE) r _CFB . This adheres to the structure of the truncated cone. An assumption is made here that during CF rupture, height of the cone is not affected; this makes the model open to fluctuations if there is any filament gap produced due to premature growth of CF. Although a forming process is considered for the device, it is not included in the model making the model not suitable for application to devices where forming is a significant factor. A possible explanation for leaving out the forming process is to avoid adding more complexity to the model because the forming parameters have to be included in the set/reset equations as well.

Similar to the previous model, the set/reset processes are described by an electrochemical redox reaction and diffusion processes [82,83,84, 105, 161, 162] which control the growth and rupture of the CF, respectively. The reduction and oxidation reaction rates are given as [93]:

Velocity of the CF radius increase/decrease which is basically the switching rate is controlled by two master differential equations which include the top (r _CF,T ) and bottom (r _CF,B ) radius of the CF. They are given as [93]:

CF radius is set a boundary which defines the limits for the CF growth/rupture which is given as [93]:

Here r_CF,TM /r_CF,BM is the maximum top/bottom radius that can be achieved and r_CF,Tm the minimum value of the top radius. This equation indicates CF geometry following a truncated cone structure and the top radius is greater than the bottom radius at all times. The model employs a numeric solving method similar to the one used in the previous model to find the discrete solutions for the master differential equations. But the solution for this model is not found, which means it is difficult to validate the reliability of the equations.

A very interesting point here is that a separate local diffusion process is not added in this model to describe the reset process in addition to the oxidation/reduction process. Many of the previously discussed models [91, 92] have a separate equation for the diffusion. But in this model diffusion has been integrated into the Eqs. (80) and (81) for redox reactions by considering different activation energies for the reduction and oxidation rates. This has been deliberately done considering the fact that the equation used in previous models to describe the diffusion velocity is similar in structure to the redox equations. As a result, the activation energies for both are combined together to consider a single activation energy which includes diffusion as well.

The current in the RRAM device is considered to be comprised of two components, the CF current (i _CF ) and the oxide current (i _buey ). In CF current, contribution of the CF resistance (R_CF ) and the oxide surrounding the CF included in the CF radius are given as [93]:

The other component in the device current, i.e., the oxide current represents the current throughout the oxide accounting for the whole area of the device except the one occupied the CF. The oxide current is described as [93]:

Here, A _HRS and α _HRS are fitting constants. The experimental and simulated results from the model are illustrated in Fig. 15 [93]. The simulated I - V curves were tuned using device parameters to fit the reported experimental data [180]. The major device parameter which was tuned the CF radius. Simulated results are almost a perfect fit to the experiment results which highlights the flexibility of the model with the proper tuning of parameters.

un RRAM current versus voltage for different cycles; the CF features are tuned to fit different curves such as r _CFTM  = 45 nm, r _CFTm  = 15 nm, and r _CFBM  = 7 nm for cycle 3, r _CFTM  = 45 nm, r _CFTm  = 25 nm, and r _CFBM  = 7 nm for cycle 31, and r _CFTM  = 40 nm, r _CFTm  = 35 nm, and r _CFBM = 7.8 nm for cycle 60; b yo - V curves for the devices under consideration along with different modeled data (the CFs’ dimensions considered to fit the whole set of curves have been generated randomly). The temperature at the most representative points of the I–V curve, highlighted in a cycle 3:T _{SET (T,B)}  = [579,656] K, T _{LRS+ (T,B)}  = [349,368] K, T _LRS-(T,B)  = [466,525] K, T _{RESET (T,B)}  = [565,621] K; cycle 31:T _{SET (T,B)}  = [573,669] K, T _{LRS+ (T,B)}  = [349,368] K, T _{LRS- (T,B)}  = [460,525] K, T _{RESET (T,B)}  = [566,622] K; cycle 60:T _{SET (T,B)}  = [564,685] K, T _{LRS+ (T,B)}  = [347,367] K, T _{LRS- (T,B)}  = [460,525] K, T _{RESET (T,B)}  = [545,641] K [93]

When considering electrochemical redox reactions, temperature is a very important factor. As can be seen from the equations describing the set/reset process, temperature parameter is a defining factor. Modeling using this concept, temperature in the CF is described through a one-dimensional approach as given as [93]:

However, this approach is not the best one for truncated cones. It leads to increased complexity and improper calculation of temperature. So, Gonzelez-Cordero et al. [93] proposed a different approach where they have developed simplified analytical equations which are suitable for simulation. In the previous Bocquet model [92] they have assumed a cylindrical uniform geometry where the calculated unique temperature is uniform throughout the CF. But to consider a more detailed physical description in this model; they have considered two temperatures in the CF. One temperature is at the hottest CF point where the reset process ruptures the CF and the other temperature is considering the main CF volume. This is a more reasonable model for a truncated cone structure as the two radii grow independently of each other depending on the oxidation/reduction processes [171, 181].

To mask this concept into a simplified analytic expression, a simplifying assumption is made. Uniform cylinders are considered to represent the truncated cone CF; one corresponding to the larger radius of the cone, i.e., r _CF,T and the other to the smaller radius, i.e., r _CF,B . So, the maximum temperature calculated for each cylinder considering the main CF volume (T _CF,T ) and at the narrowest part (T_CF,B ) is given as [93]:

Here, α _T is the conductivity temperature coefficient. The model is simulated and compared with the results from the previous Bocquet bipolar model [93] on which it is based. The results compare the findings from the model considering a cylindrical CF to the one considering the truncated cone CF. There is some evidence [93] presenting a better fit with the experimental data for this particular model as compared to the previous models where cylindrical CFs are considered and also results pertaining to the cases where multiple CFs are also presented; this shows the model’s flexibility in accommodating devices where multiple CFs is existent.

Detailed simulation results presenting the variation of current, voltage and CF radius are shown in Fig. 16 [93]. The variation of the applied voltage and the device current with time is shown in Fig. 16a. The complete simulated I-V curve can be seen in Fig. 16b which follows the form of a hysteresis loop, suggesting resistive switching behavior. Variations of the CF radius with time and applied voltage are presented in Fig. 16c, d, respectively.

un Applied triangular voltage signal (blue) and corresponding device current (red) versus time. b The RRAM current against the applied voltage. c Top and bottom CF radii versus time for the devices under consideration. d Top and bottom CF radii corresponding to figure 7(c) versus applied voltage [93]

Primary aim of the model was to be simplistic enough to be implemented in electric circuit simulators. Analytical equations are properly laid out to be used in SPICE simulations and it has been represented through a 1T/1R circuit. The model is also represented through a Verilog-A representation [93] which shows its applicability in digital circuits as well.

RRAM Models Based on Unipolar Devices

Random Circuit Breaker Network Model

In 2008, Noh’s et al. [182] proposed a random circuit breaker model to explain the switching in unipolar resistive switching devices. This model evolved to clear the considerable debate regarding the switching mechanism in unipolar devices, at the early stage of the study mechanism of resistive switching memory. Some reported that switching is the result of homogeneous/non-homogeneous transition of current distribution, while some other says due to the formation and rupture of conducting filaments. A new percolation model was reported by Noh’s group [182] in this regard which was based on a network of circuit breakers with two switchable metastable states. The device used for the study is a polycrystalline TiO₂ RRAM device. It shows wide distributions of SET and RESET voltage with uniform resistance change at the particular transition voltage. Conductive atomic force microscopy (C-AFM) tip was used as a top electrode, and external voltage was applied through it for the resistive switching operations.

The C-AFM images along with the switching curve are shown in Fig. 17a, b. The current map shows the state of the system in the LRS and HRS. In LRS, many conductive spots can be observed in the current map which can be considered to be conducting filaments. Corresponding to the HRS, the conducting spots are vanishing, which could be translated as the rupture of the filaments. This behavior was qualitatively described by a percolation model comprising of random circuit breakers (RCB), termed as RCB network model, and is shown in Fig. 17c, d. In this model, the RRAM system is considered as a combination of a number of circuit breaker, which can have either two resistance states, high (OFF state) or low (ON state). When the ON state circuit breaker receives the RESET voltage, it switches to the OFF state. Conversely, when the OFF state circuit breaker receives the SET voltage, it switches to the ON state.

Conducting AFM (C-AFM) topographs of TiO₂ films. un Schematic of C-AFM measurements (b ) Typical I - V curves of TiO₂ films using C-AFM tip as a top electrode. c At V _tip = 8 V, mapping of the current flow through the surface just after the forming operationshows locally distributed conducting regions. d TiO ₂ surface in the HRS shows locally distributed conducting regions disappear after reset operation with V _tip  = 1 V [182]

This model basically laid the foundation for future percolation models used to describe RRAM switching behavior, since it is dealt with the stochastic reversible dynamic processes. Most percolation models either investigate static cluster topology problems or dynamic percolation problems. A combination of reversible and dynamic processes is quite interesting. This also enabled future model developers to account for stochastic switching in the physical equations describing unipolar RRAM devices.

Filament Dissolution Model

Filament dissolution model was proposed by Russo et al. [82,83,84] exclusively for unipolar RRAM devices, although later revisions by the same group [85, 86] made this model suitable for bipolar devices as well. This model is based on the fundamental concept of Joule heating and filament temperature change. The model primarily focuses on the RESET transition of the devices, i.e., the transition from LRS to HRS. This is because of the high resistance associated with the RESET transition and this is where the major physical transformation in the device takes place. The model is based on the concept of conductive filament ruptured or dissolute under the effect of significant temperature change [84]. This temperature change in the filament is caused due to Joule heating. The proposed filament dissolution model has been deemed as self-accelerated due to the process of the rupture of filament accelerates by itself under suitable conditions.

Major advantage of this model is that it makes use of simple well-known coupled partial differential equations which describe the various effects in the device. The model is applied on a NiO based unipolar system [82, 83] where the oxide layer is sandwiched between two Pt electrodes and the filament is considered to grow from one end of the electrode to the other. Temperature profile in the oxide layer across its geometry is considered as parabolic; meaning that the temperature in the filament is minimum at the electrodes and maximum at the middle.

The mechanism behind the filament dissolution can be explained by the basic concept of Joule heating and dissolution which acts as an activator for the CF rupture. With the application of bias across the top electrode of the device, heat is produced in the filament due to the current flowing. The temperature steadily rises with an increase in the bias and when the bias reaches a significant level called reset voltage, the temperature rises above a value called the critical temperature. At this point, the dissolution of the filament is activated and the filament gets ruptured at a very fast rate leading to the device reaching a HRS.

Filament dissolution model uses coupled partial differential equations to describe the current in the device, the temperature changes due to Joule heating and the dissolution velocity. The current conduction in the device is described by the Poisson’s equation [83]:

$$ \nabla \times \left(\frac{1}{\rho}\nabla v\right)=0 $$ (96)

Here, ρ is the resistivity of the oxide and v the electric potential developed in the device due to the application of an external bias voltage v _term . The voltage bias is applied at one of the electrodes, while the other electrode is connected to ground which act as boundary conditions for the device. NiO is the switching oxide, and the CF is formed as a sub-layer comprising primarily of metal ions and oxygen vacancies. The CF is considered to have a diameter of φ _d .

As a result of the potential across the electrodes, heat is produced in the device due to Joule heating which leads to an increase in the temperature inside the CF. The CF geometry is divided into a number of mesh grids, to compare the ionic motions during filament formation and dissolution. The temperature is calculated at each part of the mesh grid to describe the thermal dissolution or rupture of the CF. This effect is described by the Fourier steady state heat equation given as [83]:

$$ -\nabla \times \left(k\nabla T\right)=\rho {J}^2 $$ (97)

Where k is the thermal conductivity of the oxide layer, T is the device temperature and J is current density. Thermal conductivity and electrical resistivity values are dependent on the position they are applied in. So, ρ  = ρ _CF inside the CF while it is equal to ρ _OX in the oxide layer. The same analogy applies for the values of thermal conductivity as well. The temperature is considered to be equal to room temperature T ₀ at the electrodes, i.e., they act as heat sinks.

The device temperature T increases up to the point of reset, where it reaches critical temperature, after which the CF dissolution takes place. As a result of the rupture of CF, the current conduction is interrupted which transitions the device into HRS or reset state. This dissolution factor is modeled as [83]:

$$ {\nu}_{\mathrm{DIS}}={\nu}_{\mathrm{DIS}-\mathrm{F}}{e}^{-\frac{E_a}{k_BT}} $$ (98)

Where, E _a is the activation energy, k _B is the Boltzmann constant, v _DIS-F is a fitting parameter and v _DIS is the velocity of the CF boundary toward the symmetry axis.

The temperature dependent resistivity in the CF can be described as [83]:

$$ {\rho}_{\mathrm{CF}}(T)={\rho}_{\mathrm{CF}-\mathrm{RT}}\left[1+c\left(T-{T}_0\right)\right] $$ (99)

where c is the experimentally calculated temperature coefficient of resistivity and ρ _CF-RT the standard CF resistivity at room temperature.

The coupled equations defined in the model are self-continuously solved using the numerical solver COMSOL [183, 184] Multiphysics. This software is well suited to handle these types of simulations owing to its Multiphysics capabilities. The model here uses mechanisms from electrostatics as well as heat transfer, so to simultaneously being able to handle multiple physical phenomenon’s is a big advantage on the part of the software. The obtained simulation results are shown in Fig. 18 [185]. The results show the variation in CF temperature with the change in voltage for three samples [185]. It follows the expected pattern of the model where the CF is dissolved after reaching the reset state.

Calculated temperature maps extracted from the electro thermal simulations for applied voltages a 0.4 V, b 0.45 V, c 0.575 V, and d 0.6 V for sample P (500 °C annealed); e 0.3 V, f 0.375 V, g 0.45 V, and h 0.475 V for sample Q (400 °C annealed); i 0.3 V, j 0.4 V, k 0.55 V, and l 0.575 V for sample R (300 °C annealed); m 0.15 V, n 0.25 V, o 0.375 V, and p 0.4 V for sample S (200 °C annealed) across the CF geometry of length t _buey = 20 nm. Different bias points are considered for the devices so that the reset transition for each device can be visualized easily [185]

Filament dissolution model discussed here has been modified and presented for bipolar RRAM devices by Larentis et al. [85, 86]. It is based on the same temperature and field accelerated ion migration. The set and reset processes in the device are defined by the mechanisms of drift migration induced by local electric field, ionic/electronic conduction and Joule heating. This is a point of departure from the model for the unipolar devices where the switching mechanisms for the set state were not properly defined and understood.

As a result, an equation explaining the ion migration given in Eq. (48) is introduced in the model. It includes both drift and diffusion components which rely on ion hopping. Rate of drift and diffusion is generally governed by the external applied bias and the amount of barrier lowering caused by it which is critical for the process of ion hopping. Barrier lowering generally takes place in the direction of applied electric field F , and the ion hopping depends exponentially on this energy barrier value. This causes significant migration of ions in the direction of the electric field F . This is a very simplistic and acceptable physical approach to explain the filament formation due to the charge carriers such as metal ions or oxygen vacancies. This also explains the dependence of the filament formation on the polarity of the external bias in bipolar devices. The ion-migration factor can be defined as [86]:

$$ {J}_D={J}_{\mathrm{diff}}+{J}_{\mathrm{drift}}=-{D}_s\nabla {n}_D+\mu E{n}_D $$ (100)

Where, n _D is the doping density, D _s is ion diffusivity, E is applied electric field and μ is the ion mobility. The temperature activated ion diffusivity, based on the Arrhenius law and ion mobility [85] is given as [86]:

$$ {D}_s={D}_0{e}^{-\frac{E_a}{kT}} $$ (101) $$ \mu =\frac{q{D}_s}{kT} $$ (102)

Using the above equations, the continuity equation for the drift-diffusion is given as [86]:

$$ \frac{\partial {n}_D}{\partial t}=\nabla \times \left({D}_s\nabla {n}_D-\mu E{n}_D\right) $$ (103)

Another difference from the model applied for unipolar devices is the use of conductivity values rather than resistivity. Also, resistivity had a linear variance with temperature in the previous model [82, 83] while in this case conductivity is modeled to be exponentially dependent on temperature given as [86]:

$$ \sigma ={\sigma}_0{e}^{-\frac{E_{ac}}{k_BT}} $$ (104)

The current conduction defined in Eq. (95) and Joule heating described in Eq. (96) are used in the model for the bipolar devices. This suggests an assumption that the temperature profile for both types of devices follows a similar pattern. Along with it, the current conduction mechanism is also assumed to be similar. This in a sense might be an over-simplified assumption because many of the models described for bipolar resistive switching devices have been unable to be used for unipolar RRAM devices [90,91,92, 186] due to a marked mismatch in the conduction and switching mechanisms.

The differential Eqs. (95)–(98) are simulated in COMSOL Multiphysics considering a 3-D cylindrical symmetrical geometry. The oxide material system considered here for reference is HfO₂ and the obtained simulation results are verified with the experimental results as shown in Fig. 19 [86]. Device shows bipolar switching characteristics owing to the extra terms added in the model equations. There is somewhat of a good match between the experimental data and the calculated results. This suggests that this model could be potentially used to model bipolar devices. It is quite comprehensive in its definition of device parameters and also agrees with experimental data.

un Measured and b calculated I - V curves, for reset (V > 0) and set (V < 0) transitions, obtained by applying triangular voltage sweeps. For preparing a reset state with variable R , the reset sweep applied to an initial set state with R  = 400 Ω, is interrupted at V _stop . Then, the set sweep is applied, showing that V _set increases with V _stop , hence with R [86]

Bocquet Unipolar Model

This model was developed by Bocquet et al. [90] for describing both the set and reset processes in unipolar RRAM devices. It is basically a modified extension of the model proposed by Russo et al. [82, 83] in the sense that it can model both the transitions of the RRAM device while the former only considers reset transition. For set process, a local electrochemical reduction of the oxide is considered to be responsible for formation of conductive filaments. However, the reset mechanism follows the tried and tested formula for unipolar devices which considers thermally assisted destruction of the formed metallic filaments by Joule heating as the primary mechanism. Also, it has to be mentioned that the model proposes equations which are analytical in nature and can be conveniently solved in an electric circuit solver.

The CF growth and rupture in the device is described by a local redox reaction and a thermally assisted diffusion respectively [82, 83, 105, 161] for a NiO based unipolar system. The reaction velocities for the reduction and oxidation processes during CF growth are based on the Butler-Volmer equation [162]. They are given by [90]:

$$ {\nu}_{\mathrm{red}}={k}_0{e}^{-\frac{\Delta r{G}_0+2\left(1-{\alpha}_s\right)F\left(E-{E}_{eq}\right)}{P\times {T}_{\mathrm{ox}}}}\left(1-{C}_{Ni}\right) $$ (105) $$ {\nu}_{\mathrm{ox}}={k}_0{e}^{-\frac{\Delta r{G}_0-2{\alpha}_sF\left(E-{E}_{\mathrm{eq}}\right)}{P\times {T}_{\mathrm{CF}}(x)}}\times {C}_{Ni} $$ (106)

where ΔrG₀ is the reaction free energy at equilibrium, E _eq is the equilibrium constant, α _s is the asymmetry factor, P is the ideal gas constant, k ₀ kinetics constant of chemical reaction, C _Ni is the dimensionless concentration of metallic species, and F is Faraday constant.

The above redox reaction describing Eqs. (104) and (105) can account for the CF growth, i.e., the SET process. But to explain the CF dissolution or the RESET process, it uses the filament dissolution model by Russo et al. [82, 83] as described further in the “Filament Dissolution Model” section. As discussed, the filament dissolution model uses the Joule heating mechanism to account for the reset transition in the device. The local diffusion velocity (ν_diff ) which is the governing equation for the filament dissolution model is given as [90]:

$$ {\nu}_{\mathrm{diff}}={k}_{\mathrm{diff}}\times {e}^{-\frac{E_q}{k_b\times {T}_{\mathrm{CF}}(x)}}\times {C}_{Ni} $$ (107)

As the dissolution velocity is exponentially dependent on temperature, it gets activated only when there is significant amount of temperature. Temperature value is high only when there is a comparable amount of voltage, i.e., the reset voltage has to cross a critical value during CF dissolution. This acts as means of a self-activation voltage threshold where the voltage controls the CF dissolution.

The set of reduction, oxidation and dissolution equations are coupled together in a master equation which controls the switching rate of the device given as [90]:

$$ \frac{d{C}_{Ni}}{dt}={\nu}_{\mathrm{red}}-{\nu}_{\mathrm{ox}}-{\nu}_{\mathrm{diff}} $$ (108)

The coupled equations in this model have to be solved simultaneously and continuously due to the fact that the model relies on self-consistent kinetic equations accounting for both CF growth and destruction mechanisms. This is a key feature which has to be implemented when using simulation tools to attain numerical accuracy.

Current conduction and temperature change in the device are described using simple current and heat flow equations. Current in the device on application of a voltage v _cell across the electrodes is given by [90]:

$$ {i}_0=\frac{v_{\mathrm{cell}}}{\int_0^{t_{\mathrm{ox}}}R(x) dx} $$ (109) $$ R(x)=\frac{1}{r_{CF}^2(x)\pi \left({\sigma}_{CF}(x)-{\sigma}_{\mathrm{ox}}\right)+{r}_{{\mathrm{CF}}_{\mathrm{max}}}^2\times \pi {\sigma}_{\mathrm{ox}}} $$ (110)

Current flowing in the device gives rise to a temperature due to Joule heating, and this effect is modeled as [90]:

$$ {\sigma}_{\mathrm{CF}}(x)\times E{(x)}^2=-k\frac{\partial^2{T}_{\mathrm{CF}}(x)}{\partial {x}^2}+{c}_h\frac{T_{\mathrm{CF}}(x)-{T}_{\mathrm{ox}}}{t_{\mathrm{ox}}} $$ (111)

Numerical values obtained from simulation in general profoundly interlinked between the set and reset transitions. This is because from a practical stand point the CF profile obtained after the set operation is used as the initial state to simulate the subsequent reset operación. Also, the reset current and LRS resistance depends significantly on the maximum current reached during the previous set operation [187,188,189]. This is basically due to the minimization of the CF radius which subsequently increases the resistance of the device [190].

The Bocquet unipolar model is compared against a NiO system similar to the one used in the previous Filament dissolution model [82, 83]. It is to be noted that the model is applicable for a unipolar device only. But the comparison with the NiO system is limited to a single system using a numerical solver. This is a major shortcoming in this particular model regarding the non-availability of exact experimental characteristics data comparison from other sources to calibrate the model. It means that the fitting parameters have not been tested for a variety of characterization data or other models as well. So, it is difficult to judge the accuracy and viability of the model even though it uses some interesting concepts to explain the switching process in unipolar devices.

Window Function Models

Window functions are introduced in the “Linear Ion Drift Model” section [3]. These functions are generally required to limit the values that the internal state variable can reach. The dynamics of the state variable governs the switching property of the device. So, the state variable has to be set bounds within which it can grow so that the device always remains in the permissible state and does not go out of bounds. For example, if the growth/rupture of CF is being modeled, the CF physically can only grow from one electrode end to the other. If the model growth crossed that limit, it suggests a mismatch between the physical phenomenon and the model. As a result, certain window functions [94,95,96,97,98,99] which acts as limiting functions is introduced into the model to set bounds for the device.

This is further required as near to the boundaries non-linear dopant drift effects take charge and heavily suppress the speed of the ions. As a result, it leads to a non-uniform and capricious rate of change of the state variable. To properly account for all the effects, an effective window function should have the following characteristics [98]:

• ➢Consider and take into account the boundary conditions at the top and bottom electrodes of the device.

• ➢Be capable of imposing non-linear drift over the entire active core of the device.

• ➢Provide linkage between the linear and non-linear dopant drift models.

• ➢Be scalable, meaning a range of fmax ( x ) can be obtained such that 0 ≤ fmax ( x ) ≤ 1.

• ➢Utilize an in-built control parameter for adjusting the model.

Window functions are generally multiplied with the particular I - V equation or state variable equation and are designed in a way so as to bind the complete equation within a certain boundary. There has been various window function based models proposed. In this section, we will review the most important and popular among them. We have also reported a comparative analysis among all the window functions in Table 2.

Joglekar Window

One of the very first window functions proposed by Yogesh Joglekar and Stephen Wolf [94] is based on the linear ion drift model [3]. It was developed when memristors were still in its early stages of development after the breakthrough by the HP team proposed linear model. Window functions are aim to generalize the behavior of the model around the device boundary.

Linear ion drift model [3] in general accounted for majority of the characteristics of the memristor but where it fell short was at the physical device boundary. Behavior of the devices at their physical boundaries was much more non-linear and this was left unaccounted in linear drift model [3]. Joglekar et al. [94] window function aimed to counter this limitation and generalize the non-linear behavior. The ions mobility is significantly higher in bulk of the memristor, but when it comes closed to boundary, the speed is highly suppressed. Joglekar et al. [94] proposed a modified equation for the linear model including the window function which accounted for this effect, is given as [94]:

$$ \frac{dw}{dt}=\eta \frac{\mu_D{R}_{\mathrm{ON}}}{D}i(t)F\left(\frac{w}{D}\right) $$ (112)

This modified equation reflects the speed suppression at the edges, i.e., w~ 0 or w~D . Here, F ( x ) is a window function which satisfies the conditions F (0) = F (1) = 0 so that there is no drift at the boundaries. The function is also symmetric about x = ½ and increases monotonically over the interval 0 ≤ x  ≤ 1/2, 0 ≤ F ( x ) ≤ 1 = F ( x  = 1/2). The window function F ( x ) parameterized by a positive integer p, is defined as [94]:

$$ F(x)=1-{\left(2x-1\right)}^{2p} $$ (113)

Simulated results of the window function are shown in Fig. 20 [94]. The results obtained from the window function implementation, it was understood that with a significantly high value of p , the function F ( x ) provides excellent generalization of the linear drift model ³ without any of its constraints. So, the window function works best at a high value of p, where it models the linear model accurately and also accounts for the non-linear effects at the boundary. This can be particularly seen in the variation of the function. With an increase in p , F ( x ) stays constant over an increasing interval around x  = 1/2 and when r  → ∞, F ( x ) = 1 for all x except 0 and 1 which are the boundary conditions. So, Eq. (112) models the general memristive systems perfectly at these particular conditions of p .

Theoretical i -v curves for a memristor with (realistic) dopant drift modeled by window functions F _p ( x ) = (1 − (2x − 1)) ^2p considering p  = 1 (red solid) and p  = 10 (green dashed). An external voltage v ( t ) = 2v₀ sin(ω ₀ t /2) is applied. The memristor parameters are w ₀ /D  = 0.5 and R _OFF / R _ON  = 50. The memristive behavior at p  = 10 has been enhanced. The slope of the i -v curve at a small time period is the same, R ₀ ⁻¹ , in both cases whereas the slope on return sweep depends on the window function. For a large p value, the return-sweep slope is R _ON ⁻¹  = 1> > R ₀ ⁻¹ , corresponding to a fully doped memristor [94]

This characteristic also acts as a significant limitation for the model. On the one hand, at low values of p , the window function does not perform as per expectation. But at significantly high values of p where the non-linear effects are taken into account, the difference between the linear drift and non-linear drift in the model disappears. This means that there is no proper way to account for both the linear and non-linear drift effects at the same time, while implementing this window function. Also, the mobility at the boundaries was suppressed down to zero, which made the function to be stuck at the 0 value at the terminal states. Shortcomings of this function and the improvements made over it by the Biolek window function [95] are discussed further.

Window function implemented here should be understood as extensions to the physics based models. General limitation with the physics based models is that the models do not account for the effects at the physical boundaries of the device. Memristive devices have been found to behave differently in the bulk and the boundary of the device. So, these window functions can help to overcome these limitations in the device by setting a boundary for the model and properly accounting the boundary effects.

Joglekar et al. [94] also investigated memristive device implementation in standard fundamental circuits along with the other basic elements R, L and C. Combining these four basic elements (memristor, inductor, resistor, and capacitor), the functioning of standard circuits such as MC, MLC, etc. were studied. They came to a similar conclusion as the HP team [3] that the primary property of a memristor is the memory of the charge that has passed through it. The memristor was dimensionally characterized as a magnetic flux D ² /μ_D donde D is the memristor size and μ is the mobility.

Biolek Window Function

Biolek window function [95] developed by Zdenek and Dalibor Biolek in 2009 was modeled on the proposed memristor equations by Strukov et al. [3] Its primary aim was to provide a marked improvement over the previous Joglekar model [94] and also provided a SPICE implementation of the Linear drift model [3]. They proposed changes to the way the window functions are defined so that a closer approximation between the model and the real circuit element can be achieved. They also reported the SPICE implementation of the linear drift model [3], which opened up the model to a wide range of circuit applications at that time.

Strukov et al. [3] first pin-pointed the pertinent problems in the Joglekar function while implementing in SPICE and then proposed improvements over it. First major problem with Joglekar function is its way of setting up of the terminal state R_ON and R_OFF . State equation and the window function defined, respectively, in Eqs. (113) and (114) bound the value of the variable to 0 at the boundary and it is forced to hold that value. This state cannot even be changed by an external stimulus. This happens due to the HP memristor remembering the x -coordinate of the boundary between two layers and not the amount of electric charge passed through it. As a result, when a new set or reset transition is to be started from a terminal value, the device has to start from 0 and not the actual value it had in its previous state.

Second problem of the window function is noticed when the model is implemented as an actual circuit element. The circuit component exactly remembers the entire charge which is passing through it. So, in case of the Joglekar window function, to transpose the memristor from a state x ₀ to x ₁ , a certain amount of charge q se requiere. Now the same amount of charge but in the opposite polarity, i.e., –q is required to bring the memristor back from x ₁ to x ₀ . Thus, when a memristor is being driven by a constant current with a time interval, say t , the same time t is also required for restoring the device back to its original state.

This occurs regardless of the fact that the device could be in its terminal state all the while when the current flows. This leads to significant operating delays as documented by the SPICE simulation presented by Biolek et al. [95]. Also, when the current direction is reversed, the boundaries start to move in an opposite direction regardless of the past state, thus the state is lost along another curve.

Window function being added to a particular model should be able to enhance its accountability for most of the device characteristics and get rid of the arising discrepancies. Window function proposed by Biolek et al. [95] works in this regard, defined by a positive integer p , memristor current i , and is given as [95]:

$$ f(x)=1-{\left(x- stp\left(-i\right)\ \right)}^{2p} $$ (114) $$ stp(i)=\left\{\begin{array}{c}1\kern0.5em \mathrm{for}\ i\ge 0\\ {}0\kern0.5em \mathrm{for}\ i<0\end{array}\right. $$ (115)

If the function increases the width of the doped layer, or x  → 1, the current is positive. The function value is 0 at the boundaries. When increasing the value of p , the function yields a flat window with steep troughs to zeros at x  = 0 and x  = 1.

Figure 21 [95] shows the simulated results of the window function implementation. Tuning parameter p can be used to fine tune the function accordingly to fit to different models. The obtained results satisfy all the conclusions obtained by the linear drift model [3]. The only critical limitation is the inability to account for the hard-switching effects governed by non-linear ionic drift. This means that a symmetrical hysteresis loop obtained in the models cannot be achieved. This is an inherent disadvantage the Biolek [95] window function has in its definition. So, in hindsight none of the functions were truly applicable to the proper extent in the linear drift model [3]. It was a limitation of the time which was further accentuated by the lack of detailed information about the non-linear ionic drift.

A proposed new window function demonstrated in the model (taking p  = 2) [95]

Benderli-Way Window Function

Another window function based on the basic HP model was proposed by Benderli and Wey [96] in 2009. The end result they set out to get was similar to the Biolek function [95]. The developers wanted to develop a SPICE compatible macro model based on the HP memristor which would be suitable for applications in circuit simulations. They proposed a clipping circuit which will bind it within the constraints of the length of the device (D ).

The proposed clipping circuit was comprised of four comparators whose job was to ensure that the state variable function w ( t ) did not go beyond its limits. The comparators clipped w ( t ) at its top and bottom boundaries. It basically acts as a switch in which if the comparators detect a certain value in the device they activate a switch and set the device at a particular voltage. So, when w ( t ) reached the upper boundary of the device, the device is connected to a voltage source of value D , which effectively clips w ( t ) at D . This operation occurs when the voltage bias is positive.

At negative voltage bias when w ( t ) approaches the lower boundary, the circuit is connected to ground, thus clipping w ( t ) at 0. This clip is enforced until the voltage polarity is changed suggesting the correct operation is being performed. Non-linear effects at the boundaries are modeled by a proposed window function which takes into account the non-linear dopant drift. The window function is defined as [96]:

$$ f(x)=\frac{w(t)\left(D-w(t)\right)}{D^2} $$ (116)

Also by increasing the capacitance near the boundaries, the non-linear effects could be accounted for it in the circuit. The major shortcomings of this function are its simplifying approximations and the lack of a clear description of how the linear and non-linear drift can be modeled in the circuit. Although it manages to obtain a hysteresis relation for the device, it suffers from similar limitations as the Joglekar model [94]. Additionally, lack of clear information regarding the non-linear effects was an equal deterrent to the application of the function in circuits.

Shin Window Function

Shin, Kim, and Kang [97] in 2010 tried to circumvent the issue of window functions by proposing a constitutive relationship derived from the basics of the memristors developed by Chua [1]. This is different from the previously reported window functions in the sense that they tried to model the memristors perfectly by relating charge and flux together. This was the fundamental essence of the Chua model and is a stark contrast to the linear ion drift mechanics proposed by Strukov et al. [3].

Utilizing the relationships between flux (φ ) and charge (q ) in a current-controlled and voltage-controlled memristor, the window function developed keeps the model checked within the bounds. Chua model for a current controlled memristor is defined as [97]:

$$ \frac{d\varphi}{d t}=\frac{d}{d t}\left[f(q)\right]=\left\{\frac{d}{d q}f(q)\right\}\times \frac{d q}{d t} $$ (117) $$ {v}_M=\left\{\frac{d}{dq}f(q)\right\}\times {i}_M\equiv {R}_M(q)\times {i}_M $$ (118)

Here R _M (q ) is the memristance defined by a derivative of the charge-flux relationship with respect to the charge. Thus, R _M (q ) = df (q )/dq defines it as a current controlled memristor.

Similar to the above equations, a voltage controlled memristor is defined in terms of the charge-flux relationship q  = g (φ ) as [97]:

$$ \frac{d q}{d t}=\frac{d}{d t}\left[g\left(\varphi \right)\right]=\left\{\frac{d}{d\varphi}g\left(\varphi \right)\right\}\times \frac{d\varphi}{d t} $$ (119) $$ {i}_M=\left\{\frac{d}{d\varphi}g\left(\varphi \right)\right\}\times {v}_M\equiv {G}_M\left(\varphi \right)\times {v}_M $$ (120)

So, on similar terms, G_M (φ ) is a voltage-controlled memductance whose values can be calculated by measuring slopes of charge-flux relationship g (φ ).

The above written Eqs. (118) and (119) can be used in compact models for circuit implementations. But they are inadequate when it is needed to define it within a bounded resistance range. This is where window functions are written to modify the circuit parameters so that the model operates within the resistance range of R _MIN y R _MAX . Thus, in mathematical terms it means that R _M ∈ [R _MIN, R _MAX ]. Memristor needs to be confined within the available range of resistance so as to adhere to design requirements. When the device reaches one of its boundary values, it has to stay in that state after any excess charge or flux is applied to the device. This has to be ensured so that the device does not violate its boundary conditions under hard switching conditions.

This condition is obtained in this case by ignoring the excess input current or voltage. So, in the model it is masked by a window function H so that the value of q or φ always stays within the available range. The masked input current and input voltage in the case of a current controlled and voltage controlled memristor are respectively given by [97]:

$$ H\left({i}_M\right)=\left\{\begin{array}{c}{i}_M,\kern5.5em \mathrm{if}\ {R}_M\in \left({R}_{\operatorname{MIN},\kern0.75em }{R}_{\mathrm{MAX}}\right)\\ {}0,\kern0.5em \mathrm{else}\ \mathrm{if}\ {i}_M\ \mathrm{does}\ \mathrm{not}\ \mathrm{pass}\ \mathrm{zero}\end{array}\right. $$ (121) $$ H\left({v}_M\right)=\left\{\begin{array}{c}{v}_M,\kern5.5em if\ {R}_M\in \left({R}_{\operatorname{MIN},\kern0.75em }{R}_{\mathrm{MAX}}\right)\\ {}0,\kern0.5em \mathrm{else}\ \mathrm{if}\ {v}_M\ \mathrm{does}\ \mathrm{not}\ \mathrm{pass}\ \mathrm{zero}\end{array}\right. $$ (122)

Above equations disregard any excess input current or voltage in the model space. The boundary state of the device is held until the polarity of the input source is reversed. When the polarity is reversed, it indicates the start of a new transition; the function will force the memristor to move back into memristive region. The operation is similar to the clipping circuit proposed in the Benderli-Way function [96] discussed previously, but here, there is no requirement of a complicated comparator circuit. The major purpose of using this approach to model the devices was to remove the usage of a special window function as done previously but still be able to adhere to the boundary conditions implicitly.

But the developers were also wary of the fact that some devices may not be modeled properly by their proposed constitutive relationship, which have a non-constant dynamic state variable. This was made in order to account for the non-linear drift effects discussed and experimentally shown in the HP memristor device [3].Hence, in order to take into account, the non-linear drift effects as well, they proposed a special window function. The window function aimed to overcome the backing problem faced by the Joglekar [94] function discussed previously. It has been described as [97].

$$ W(q)=\delta +\left\{1-{\left(2z-1\right)}^{2p}\right\} $$ (123)

where δ is a non-zero positive constant (δ < < 1), z is a normalized memristive charge, and z defined as z  = (q  − q _MIN )/(q _MAX  − q _MIN ). The non-zero value of δ ensures that the state of the model returns to the normal memristive region when the input polarity is reversed. This removes the backing problem suffered by previous window functions.

Prodromakis Window Function

Prodromakis window function was proposed in 2011 by Themis Prodromakis et al. [98] of Imperial College, London, which aimed for a simple and efficient function modeled the memristor device characteristics [191] effectively. Some of the limitations and constraints of the previous models were alleviated which made the function easy and accurate to use.

The window function is considered to be parabolic in nature. It also employs a control parameter (p ) in the exponent which provides the model with the required scalability and flexibility. It also makes the window function f ( x ) scale upwards or downwards which helps create a family of distinct curves. The function f ( x ) is given in terms of p ∈ R ⁺ as [98]:

$$ f(x)=1-{\left[{\left(x-0.5\right)}^2+0.75\right]}^p $$ (124)

Control parameter is critical in this function as it helps to remove many of the constraints and limitations of the previous functions. The function can scale upwards due to the control parameter, which suggests that f _máx ( x ) can take any value between 0 and 1 inclusive. Also a very large value of p provides a linkage with the linear dopant drift effects. A serious limitation with the previous Joglekar [94] and Biolek [95] models was that the control parameter was allowed to take only integer values. But here, p could have real values as well which added more flexibility to the model.

Results obtained from the function suggest it returning a zero value at the active bi-layer edges. The drift of the dopants is also suppressed near the metal interfaces. This accounts for the non-linear drift effects at the boundaries. Other major problem of zero value stuck at the terminal state is tackled by implementing a feedback path as suggested in Eq. (123). The function can also be adjusted for any peculiar cases such as when f _máx  > 1 with the help of a second control parameter j given as [98]:

$$ f(x)=j\left(1-{\left[{\left(x-0.5\right)}^2+0.75\right]}^p\right) $$ (125)

Hysteresis loop obtained using the window function is asymmetrical which has been explained by Prodromakis et al. [98] as a result of the different switching rates of the ON and OFF rates which is quite reasonable. The hysteresis also suggests there is no terminal state problem as highlighted in the Joglekar function. Width of the doped region does not go higher than D , and the memristance is correctly limited. In the case of reverse polarity, it does not get stuck at a zero value and does not take any error states as highlighted by the results.

Results presented by Prodromakis et al., shown in Fig. 22 [98], suggest that the function can be adjusted and scaled effectively using the two control parameters in the function. The p parameter supports lateral scaling while the j parameter supports vertical scaling. Value of p has been increased from p  = 1 to p  = 80 and the results are shown in Fig. 22, which showcase the scaling features of the window function. This window function is one of the most efficient and accurate among the others which have been reported so far.

Plot of the proposed window function f ( x ) = 1 − [(x  − 0.5)2 + 0.75] ^p against the normalized width of the doped region x (p is considered as a variable) [98]

Boundary Condition-Based Model (BCM) Function

This window function developed in 2012 by Fernando Corinto and Alan Ascoli [99] was aimed at improving the models proposed by Joglekar [94] and Biolek [95]. They identified possible limitations with the previous functions with respect to their exhibition of single-valuedness and multi-valuedness, respectively. Also tuning the range and the boundary conditions were not possible with the Joglekar [94] and Biolek [95] functions. This was handled by the BCM window function by deriving novel methods to propose closed-form solutions for memristor devices. Along with that they also added tuning parameters to increase the flexibility of the boundary conditions used in the models.

Design of the Joglekar function limits it to a single value of memductance-flux characteristics at all input values. Similarly, the input dependent Biolek function limits it to only multi-values of the function under sign varying input. But the BCM function allows for both single-valued and multi-valued memductance-flux characteristics under a single sign varying input. The function assumes a linear dopant drift effect, which simplifies the analytic integration as well as makes it suitable for closed form solutions under any initial condition state. But this invariably neglects the non-linear boundary effects in the device. So, on the one hand, the BCM function proposes a very simplified expression for defining the boundaries of the device but it misses out on accounting for the non-linear effects due to the simplifying assumptions.

BCM model uses tuning parameters in the window function equations in an attempt to account for the non-linear effects [3, 46, 172, 173]. But they are not accurate to the same level as implicit definition of those effects. Another assumption which further simplifies the model to allow for closed form solutions is that the ionic drift rate remains constant. The BCM model is based on a window function having unitary values for all x ( t ) ∈ (0,1), and also exhibiting the vertical transitions as [99]:

$$ 1\longrightarrow 0\left\{\begin{array}{c}\mathrm{if}\exists t\mid x(t)=1\kern0.5em \mathrm{for}\ \eta v(t)\ge -{v}_{th,1}\\ {}\mathrm{if}\exists t\mid x(t)=0\kern0.5em \mathrm{for}\ \eta v(t)\le {v}_{th,0}\end{array}\right. $$ (126)

and

$$ 0\longrightarrow 1\left\{\begin{array}{c}\mathrm{if}\exists t\mid x(t)=1\kern0.5em \mathrm{for}\ \eta v(t)<-{v}_{th,1}\\ {}\mathrm{if}\exists t\mid x(t)=0\kern0.5em \mathrm{for}\ \eta v(t)>{v}_{th,0}\end{array}\right. $$ (127)

Here, η is a linear control parameter and ∃ is a quantifier denoting “there exists” which signifies that for x ( t ) there exists exactly one solution. Values of the non-negative parameters v _th,1 y v _th,0 determine the occurrence of such transitions. The conditions in the first Eq. (125) are established when x ( t ) obeys the boundary conditions x  = 1 and x  = 0. But the conditions in the second Eq. (126) is established when the function x ( t ) no longer obeys these boundaries.

These dynamics are encapsulated in the mathematical description defined by the tunable boundary conditions C _n (n  = 1, 2, 3) as [99]:

$$ {C}_1=\left\{\begin{array}{c}x(t)\in \left(0,1\right)\ \mathrm{or}\ \left(x(t)=0\ \mathrm{and}\ \eta v(t)>{v}_{th,0}\right)\\ {}\mathrm{or}\ \left(x(t)=1\ \mathrm{and}\ \eta v(t)<-{v}_{th,1}\right)\end{array}\right\} $$ (128) $$ {C}_2=\left\{x(t)=0\ \mathrm{and}\ \eta v(t)\le {v}_{th,0}\right\} $$ (129) $$ {C}_2=\left\{x(t)=1\ \mathrm{and}\ \eta v(t)\ge {v}_{th,1}\right\} $$ (130)

So, in case of η = + 1, v _th,0 is the threshold voltage. This is the minimum value of the input needs to cross, after it enters the positive region. Similarly, v _th,1 represents the negative region, for η = − 1. This means that the conditions in Eqs. (127) and (128) have to be met first, before the conditions in Eq. (129) holds. The window function F is thus defined as [99]:

$$ F\left(x,\eta v\right)=\left\{\begin{array}{c}1;\mathrm{if}\ (127)\ \mathrm{holds}\\ {}0;\mathrm{if}\ (128)\ \mathrm{or}\ (129)\ \mathrm{holds}\end{array}\right. $$ (131)

The window function qualitatively works similar to other functions. At the boundary conditions, the vertical transition from 0 to 1 or vice versa occurs depending on the polarity of the input stimulus. Thus, the input used here is sign varying in nature.

Model Verification

Several ways are there to verify the working of the presented models, in this work. Some of the implementations and verification have been included with the description of the model. Models which have been described quantitatively using mathematical equations can be verified by solving the equations in a simulator. Generally, I - V characteristics of the simulated model are compared with the corresponding experimental data from the device. This gives a fair idea on the reliability of the model. Physical models described by mathematical equations can be solved by a multitude of solvers such as MATLAB, Mathematica, COMSOL, etc. Compact models which have been translated to work in the circuit level are generally simulated using the SPICE framework. There are a variety of SPICE based simulators in the market such as HSPICE, Ngspice, etc. which could be utilized. Corresponding output characteristics can be matched with the experimental results.

Physically verifying the switching mechanism in a model is trickier. It generally involves in-situ observation of the switching process [192] which requires a lot of precision and high-end equipment. However, it is very solid evidence regarding the viability of the switching mechanism presented in the model. There have been some novel methods reported to observe the growth of the conductive filaments during the switching process. Conductive atomic force microscopy (C-AFM) [182] has been used to visualize the formation and rupture of conducting filaments. yo - V switching curve is shown in Fig. 17, which is clearly shows HRS/LRS states and the corresponding state of the filaments. Electrostatic force microscopy (EFM) [193] can also be used to visualize the migration and accumulation of oxygen ions by calculating the electrostatic force between the probe and the sample. It is one kind of in-situ TEM, where the focus is primarily on the charge of the carrier. Formation of conductive channels can be observed by high resolution energy-dispersive X-ray spectroscopy (EDS) which can provide accurate detailing of the composition in the filaments. These two methods have been proven to be quite effective in verifying the physical switching mechanism and the visualization of conductive filaments in RRAM devices.

Well-Posed Memristive System Definitions

An excellent work published in 2015 by Wang and Chowdhury [100] of UC Berkeley set a new paradigm for memristor and RRAM modeling. It was a push in the right direction for the whole memristor modeling community. The major features of the work were the significant improvements made on the pre-existing models. Some tweaks were proposed in the fundamental understanding of memristor models which contribute to eradicating some of the long-standing issues which has plagued the memristor models. Also, they demonstrated implementation of the models in SPICE, Verilog-A as well as their own prototyping platform based on MATLAB called as MAPP [191].

The root cause of many of the issues affecting memristor models were improper mathematical implementation. As a result, it limited their application in a variety of simulation and design scenarios. Simulation of the models is a very critical factor in determining the applicability of the model; however, the existing models were unable to be applied in a variety of simulation studies such as DC, transient etc. This work aimed at modifying the models into a form where simulating them would be a simple task.

The common ill-posed or erroneous definitions that many of the previous models suffer are not being properly defined at all biases, outputs not being unique or continuity problems. These basic problems are to be avoided in the models for wide application. All the various problems that the authors have encountered and the improvements they have presented are discussed.

A very valid point highlighted by them is that a well-written model for a particular circuit should be able to replicate its characteristics or be valid beyond its actual boundaries as well. Even if getting outputs beyond the applicable range might not be physically possible, but in simulation environments like SPICE it is imperative that the circuits work at all level of biases. This will lead to efficient simulation and produce smooth varying outputs at all biases. Many of the models we have discussed previously have sought to ignore the operation of the model beyond the available range, leading to their incompatibility for use in circuit simulators. This requirement can also be understood by the underlying algorithms of circuit simulators such as the Newton-Ralphson (NR) algorithm [194] which is commonly used for solving non-linear equations.

The NR algorithm [194] works on the principle of applying a sequence of biases to devices so that convergence can be reached for a valid solution of the circuit. So even if the bias is physically possible or not, for the NR algorithm to find a solution the model must evaluate at all bias. This invariably means that if the NR finds a solution to a well-designed model, the input bias will be physically reasonable. In the cases where the converged solution is physically not possible, it provides insights into the problematic areas of the models and is critical to troubleshooting. Therefore, in order for the NR algorithm to work correctly and find a proper solution, models should be designed to be evaluated at all biases.

Another fundamental problem is the divide-by-zero error. Many of the models have terms such as 1/(x − a ) which cause these errors. It leads to the solutions getting unbounded and causing discrepancies. Along with that, some expressions use square root (√) with negative arguments which give rise to complex arguments. Non-real numbers are not valid arguments for models and can cause non-convergence to a valid solution in simulators.

Almost all of the models we have discussed earlier do not account for a very important aspect of device model simulations. Any mathematically viable input must produce a mathematically viable output, and the most basic among this is the DC analysis. It is commonly the crux and starting point of any analysis and a proper model should produce an accurate DC solution. A proper well-designed model should work consistently with all kind of analysis. But almost all of the models suffer from significant DC response problems. So, this is another area that needs to be improved. Wang y col. [100] also addresses the problems of the models generally faced when being defined in circuit level languages such as Verilog-A.

The crux of the improvement to the models Wang et al. [100] have proposed revolve around the correct way of modeling hysteresis itself, i.e., using internal unknowns and implicit equations. This is because the dynamics of the filament in a RRAM closely follows a hysteresis characteristic. Also, the improvements make the models simulation ready with all the major analyses like DC, AC, and transient providing acceptable results. Various techniques are also proposed to aid convergence in electric simulators including a proposed new limiting function which replaces the functionality of window functions and overcomes all their limitations.

Accurate Description of Hysteresis

Memristive systems like RRAM devices have been proven to follow the features of hysteresis [1, 2] very closely. So, tweaking the models should start from the very basic. Modeling and defining the hysteresis characteristics accurately are critical to the proper functioning of RRAM models. Devices with hysteresis do not have a simple algebraic mapping between the voltage v ( t ) and current i ( t ). A state variable s ( t ) which defines the state of the device is required for the mapping given as [100]:

$$ i(t)={f}_1\left(v(t),s(t)\right) $$ (132)

A differential equation is used to govern the dynamics of the state variable, i.e., the rate of change of the state described as [100]:

$$ \frac{d}{dt}s(t)={f}_2\left(v(t),s(t)\right) $$ (133)

The above equations serve as the model template for modeling hysteresis. Value of the s ( t ) at a particular time instant t is governed by the history or the state of v ( t ). Thus, the device is considered to be having an internal memory of the input voltage. Functions f₁ and f₂ are chosen accordingly to define the characteristics of the device. Choice of these functions could be termed as critical in defining the dynamic of the device.

There has been a very common thought that hysteresis shows up in transient analysis only [1,2,3]. But Wang et al. [100] have demonstrated hysteresis based on the DC solutions of the models. This has been achieved by the proper selection of the governing functions f ₁ and f ₂ . As discussed earlier, obtaining a DC response is imperative for the proper analysis and simulation of the devices. So, obtaining a hysteresis in a DC solution makes this concept highly noteworthy. The governing functions are defined as [100]:

$$ {f}_1\left(v(t),s(t)\right)=\frac{v(t)}{R}.\left(\tanh \left(s(t)\right)+1\right) $$ (134) $$ {f}_2\left(v(t),s(t)\right)=\frac{1}{\tau }.\left(v(t)-{s}^3(t)+s(t)\right) $$ (135)

In the function f₁ , hyperbolic tangent function tanh is chosen because of its monotonically increasing properties with a range of (− 1,1). The dynamics of s ( t ) is governed by the choice of f ₂ . When the value of f ₂ is 0, the corresponding (v , s ) values are the DC solution of the circuit. As shown in the contour plot in the Fig. 23a, [100] the curve f ₂  = 0 folds back in the middle and crosses the v  = 0 axis three times. Thus, it has three stable states or three possible values of s in the DC solution which forms the foundation of the DC analysis of hysteresis. The operation of the variable s as a hysteresis curve is shown in Fig. 23b [100]. As s modulates the current, the I - V relationship will result in a hysteresis as well. Multiple stability of the state variable and abrupt changes in the DC solutions leads to the formation of a hysteresis in DC analysis. This sets a very strong foundation for accurate and efficient modeling of RRAM devices.

f ₂ function in (6) plotted in contour form and predicted s -v hysteresis curve depending on the sign of f ₂ [100]

Proper Definition of Internal Unknown Variables in Verilog-A

It has been discussed several times during review of the various models, implementation of the models accurately in SPICE [116, 117] and Verilog-A [142] is critical for their acceptability. This is because SPICE is the most commonly used circuit simulation platform, and Verilog-A is the most widely used hardware description language. So, simulating in these platforms is as close as it gets to the real physical devices. A major shortcoming of the previous models is the way that internal unknown variables were handled in Verilog-A.

In a memristor model, the state variable is an internal unknown because its value gets changed with different states. Verilog-A does not have a straightforward way of handling and defining these unknowns. As a result, it can get very confusing while dealing with a constantly changing value. Wang y col. [100] proposed to declare the state variable s as a voltage or potential rather than any “real value.” Some very critical points are mentioned below while handling internal unknowns in Verilog-A.

• Different Verilog compilers handle variables declared using “real” differently. Then, this will lead to very inconsistent results.

• Differential equations should not be defined by using the in-built idt() function. Because this function has very inconsistent support in the compilers and causes many limitations [140, 142].

• Time integration to obtain analytical solutions should not be coded inside the model. The process is pretty simple but has serious pitfalls as given below.

• ➢This method makes use of “abstime” function. To define the starting point of the integration it also uses “initial_step.” These have been termed as bad modeling practices in analog simulation [140, 143].

• ➢The internal unknowns are defined as a memory state in this method, which can create problems for periodic steady state (PSS) analysis.

• ➢This method bypasses many of the simulators built in facilities such as the convergence aids, time step control etc.

• ➢It can cause serious convergence issues for stiff systems due to its dependence on the Forward Euler (FE) algorithm [195].

These problems are generally a combination of bad modeling practices and the incapability of Verilog-A to handle internal unknowns efficiently. As a result, declaring s as a voltage has been demonstrated as an effective way of getting around the problem.

Developing Generic RRAM Models

Taking the previously discussed hysteresis equations as a template, Wang et al. [100] presented a generic way of developing compact models for RRAM devices. To demonstrate the RRAM model, ASU/Stanford model [78, 80] is considered.

The filament gap is used as state variable. Current across the device is considered as itb , the voltage as vtb and the unknown state variable as gap . The equations are defined based on the previously stated hysteresis template given as [100]:

$$ itb(t)={f}_1\left( vtb(t), gap(t)\right) $$ (136) $$ \frac{d}{dt} gap(t)={f}_2\left( vtb(t), gap(t)\right) $$ (137)

The above equations define the physical contexts of the RRAM device. Now, choosing the proper functions for f ₁ and f ₂ is critical to capture the physical properties. In most of the models, we have encountered till now the application of consistent equations for f ₁ with some changes in the internal unknowns. This fact is corroborated by Wang and Roychowdhury’s work [100]. They have considered the f ₁ function as [100]:

$$ {f}_1\left( vtb, gap\right)={i}_0\times \exp \left(-\frac{gap}{g_0}\right)\times \sin h\left(\frac{vtb}{V_0}\right) $$ (138)

The f ₂ function is considered according to the ASU/Stanford model [78]:

$$ {f}_2\left( vtb, gap\right)={v}_0\times \exp \left(-\frac{E_a}{V_T}\right)\times \sinh \left(\frac{vtb\times \gamma \times {a}_0}{t_{ox}\times {V}_T}\right) $$ (139) $$ \gamma ={\gamma}_0-\beta \times {gap}^3 $$ (140)

The γ here is the local field enhancement factor [196] which contributes to the abrupt SET (filament growth) and gradual RESET (filament dissolution). A common property among most of the RRAM models is the fact that the sign of f ₂ is same as that of –sinh (vtb ). This means in terms of the gap that it starts to decrease whenever vtb is positive and vice versa. But this growth or dissolution cannot be indefinite for numerical simulation to work. For the simulations to work in reality, they have to be bounded which has been discussed in depth in the next section.

Various methods have been proposed to account the boundary effects of the devices. It will come up short with the methods having some serious limitations. Some of the models have implemented direct “if-then-else” statements in the Verilog code [80, 81]. But the problem is that the use of “if-else” statements removes the model from the differential equations framework which is not acceptable. It also introduces hard discontinuities in the model whereas we need smooth continuous curves.

A very popular way of modeling boundary conditions is the use of window functions which we have discussed in the “Window Function Models” section. The window functions as discussed work on the principle of setting f ₂  = 0 when gap = maxGap and minGap . We have seen that improvements made with the window functions make them suitable for transient simulations and they produce smooth and continuous results. But the real problem with these functions is actually deep-rooted. The problem can be understood by analyzing the sign and zero-crossings of the f ₂ curve as shown in Fig. 24 reported by Wang and Chowdhury [100].

Different choices of f ₂ considered in the RRAM model [100]

As is seen in Fig. 24b, the f ₂  = 0 curve has three lines, V  = 0, maxGap , and minGap . Beyond the values of maxGap and minGap , there will not be any stable DC solutions so they will not show up in transient analysis. Therefore, when sweeping between those values, the transient solution will work fine with the window function multiplied to f ₂ . But with the other analyses there are several problems. In DC operating point analysis, unphysical solutions can show up owing to the fact that all the lines consisting of the f ₂  = 0 curves are valid. So, when the voltage is zero, any arbitrary value can be the value of gap and it no longer follows the boundaries. Hence, DC analysis is a major limitation for the window functions.

This has been tackled very efficiently by Wang et al. [100] by keeping the DC solutions in a single continuous curve and then trying to bind the value of gap . As illustrated in Fig. 24c, curves A and C contain the stable solutions and B has the unstable solutions. So, when the sweeping starts after zero, the value of gap will switch between maxGap and minGap . It is mathematically represented by introducing a couple of clipping functions F _clipmin and F _clipmax represented as [100]:

$$ {f}_2^{\ast}\left(\mathrm{vtb},\mathrm{gap}\right)={f}_2\left(\mathrm{vtb},\mathrm{gap}\right)+{F}_{\mathrm{clipmin}}\left(\mathrm{vtb},\mathrm{gap}\right)+{F}_{\mathrm{clipmax}}\left(\mathrm{vtb},\mathrm{gap}\right) $$ (141) $$ {F}_{\mathrm{clip}\mathrm{min}}\left(\mathrm{vtb},\mathrm{gap}\right)=\left(\mathrm{safeexp}\left({K}_{\mathrm{clip}}\times \left(\mathrm{minGap}-\mathrm{gap}\right),\mathrm{maxslope}\right)-{f}_2\left(\mathrm{vtb},\mathrm{gap}\right)\right)\times {F}_{w1}\left(\mathrm{gap}\right) $$ (142) $$ {F}_{\mathrm{clip}\mathrm{max}}\left(\mathrm{vtb},\mathrm{gap}\right)=\left(-\mathrm{safeexp}\left({K}_{\mathrm{clip}}\times \left(\mathrm{gap}-\mathrm{maxGap}\right),\mathrm{maxslope}\right)-{f}_2\left(\mathrm{vtb},\mathrm{gap}\right)\right)\times {F}_{w2}\left(\mathrm{gap}\right) $$ (143)

The functions F _w1 and F _w2 are smoother versions of the step functions. This serves the purpose of adhering to the property required for a clipping function to work properly while maintaining a smooth continuous curve. These are described as:

$$ {F}_{w1}\left(\mathrm{gap}\right)=\mathrm{smoothstep}\left(\mathrm{minGap}-\mathrm{gap},\mathrm{smoothing}\right) $$ (144) $$ {F}_{w2}\left(\mathrm{gap}\right)=\mathrm{smoothstep}\left(\mathrm{gap}-\mathrm{maxGap},\mathrm{smoothing}\right) $$ (145)

The functions safeexp () and smoothstep () are smoother versions of the normal variants of the exp () and step () function. They have been developed by Wang et al. [100] in their MAPP [191] platform and is available to use within the platform. The clipping functions here closely mimic the actual physical effects occurring in the device. It can be termed as a huge force which keeps the state variable gap within its acceptable physical limits.

The templates provided by Wang and Chowdhury [100] for RRAM modeling is capable of widespread applicability and can be used as an ideal platform to develop other models. It consists of accurate modeling of hysteresis, includes proper handling of internal unknowns in Verilog-A and does not need to use incompatible functions like “idt()” and “initial_step” in the differential equation framework. They also circumvent the various limitations of the window functions by the use of mathematically accurate clipping functions. The model templates support a variety of analyses such as DC, AC, transient, PSS, and homotopy [197] in Verilog-A, SPICE, and MAPP. This is a very exhaustive list of advantages which should be used for development of future models by the RRAM community.

Improving Solution Convergence

Obtaining solution convergence is one of the most important features, every RRAM model ought to possess. The convergence of the solution points to the fact that it is valid and acceptable. This has been a problem for many of the compact models describing RRAM devices. Several techniques have been proposed [100] to aid convergence of solutions in these models. The use of limiting functions compatible with SPICE is very important so that it limits the solutions whenever they cross the acceptable range. A very simple way to make sure the solutions converge is to properly scale the unknowns and variables. Proper scaling makes sure that any results obtained are defined relatively accurate to the inputs.

But the major feature that aids convergence is the use of proper limiting functions in SPICE. Generally during DC operating point analysis, the Newton Raphson (NR) [194] algorithm iterations take very large values while guessing the DC solutions. This is because of the fast-growing sinh functions used in the models. This in turn leads to large values of f ₁ and f ₂ . Limiting functions are the best technique to circumvent this and prevent the NR iterations from taking large input values but keeping the sinh functions intact. Presently, SPICE includes pnjlim , fetlim , and limvds as limiting functions. But these are not enough to obtain numerical accuracy for RRAM models. A new limiting function dubbed as sinhlim has been proposed for this purpose. The function is based on the original pnjlim and is given as [100]:

$$ {x}_{\mathrm{lim}}=\mathrm{sinhlim}\left({x}_{\mathrm{new}},{x}_{\mathrm{old}},k\right)=\frac{1}{k}\times \ln \left({y}_{\mathrm{lim}}+\sqrt{1+{y}_{\mathrm{lim}}^2}\right) $$ (146)

The major feature here is that the limiting function does not use very large values of x _new , instead it increases the value in iterations. This function is fully compatible with SPICE and can be implemented in any SPICE compatible circuit simulator. This marks a new addition to the number of limiting functions available for use in circuit simulators apart from the ones developed decades ago.

Improving Existing Models

On the basis of the accurate generic model templates discussed in the previous section, Wang and Chowdhury [100] have proposed some improvements in some of the existing popular models. We discuss them here and present it in a concise form so that it becomes easy to understand the changes and implement it forward in future models.

The improvements have been proposed on the linear ion drift [3], non-linear ion drift [46], Yakopcic [73, 74], TEAM/VTEAM [75,76,77], and ASU [78] models. Many of the models have I-V relationships in common with other models. Therefore, the most common and important among those are considered and termed as f ₁ functions as discussed earlier. The state variable equations of those models are termed as f ₂ functions are corresponding improvements are proposed.

Both the f ₁ and f ₂ functions are generally non-linear and asymmetric. And the reported models use very discontinuous and fast-growing terms like exponential, sinh , and power (pow ) functions which results in difficulties during the convergence of the solutions. So, these can be overcome by using “smooth” and “safe” functions as proposed by Wang et al. [100]. The smooth functions are used in place of discontinuous functions. Major design criteria in the smooth functions used is a common smoothing factor which combines common elementary functions to approximate the original non-smooth functions. Smoothing factor controls the trade-off between better approximation and more smoothness. The safe functions are versions of the fast-growing functions which limit the maximum slope the functions can attain, and then linearize it to keep the slopes constant beyond it. For some functions like sqrt , log , etc., the “safe” versions clip the inputs using smoothclip so that non-valid outputs can be avoided.

In the particular f₁ and f₂ functions used in the models, the if-else statements are replaced by smoothswitch which removes the discontinuity of the former. The exp and sinh functions are correspondingly replaced by safeexp and safesinh . The authors have demonstrated the definition of the functions in MAPP and Verilog-A which makes it easy for future model developers to integrate it into their system. A very common problem with the f ₂ functions, i.e., the state variable dynamics is the uncertainty over the range of the internal unknown. The previous models counter this by either introducing window functions that bound it within a range or do not account for the effects at all. This has been very efficiently handled by introducing self-modeled clipping functions which define the acceptable range of the internal unknown.

Another major problem which is countered is the poor DC hysteresis response of the models. There have been some attempts in the Yakopcic and TEAM/VTEAM models to model this effect by setting the value of f ₂  = 0 within a certain voltage range. It has been discussed that DC hysteresis occurs due to the DC solution curve folding backwards. With the approach used in the other two models, when f ₂  = 0 with the voltage close to 0, there are infinite number of solutions for the state variables within that voltage range. This is a problematic approach in those unstable DC solutions are also included here which makes the equation system ill-conditioned. This will also cause problems during DC operating point analysis and homotopy analysis. The DC solutions of the models will vary from simulator to simulator because of the manner the equations are designed. This has been very efficiently handled and circumvented by Wang and Chowdhury [100] as summarized in Tables 3 and 4.

Novel RRAM Applications

There have been several new breakthroughs with RRAM architectures and applications. Among them noteworthy in case of architectures is the use of materials such as graphene, amorphous carbon films, transition metal dichalcogenides (TMDs), black phosphorous in a RRAM device. Neuromorphic computing is a novel application scheme for RRAM devices which utilize the memory retention property to use them as synapse devices.

RRAM devices based on graphene and related materials [198] have showcased performance similar to conventional metal oxide devices. These devices are different due to their unique lattice structure and belong to a completely different family of materials. So, investigation on modeling of such devices is highly necessary. Whether conventional modeling techniques such as the ones presented in this work can be used for these devices depends on the physical phenomena governing them and the corresponding I - V charcateristics. The hypotheses presented to explain the switching in graphene oxide (GO) based devices [198] are consistent with standard bipolar RRAM switching mechanism. The absorption and creation of the conductive filaments are thought to be a result of diffusion of metallic ions from the electrode to the switching layer or transport of oxygen related carriers in the switching media. RRAM devices based on amorphous carbon [198] as the switching media are thought to operate under a similar mechanism. Owing to the similarity in the nature of the physical transport mechanisms, existing physical models can be used to explain the novel GO and amorphous carbon based devices.

Neuromorphic computing [199] is a novel architecture scheme which employs RRAM as synapse devices as its fundamental building block. It is believed that these RRAM based neuromorphic systems can replicate how our brain functions, harnessing the ability of memristors to remember their state. This enables the system to be trained for specific applications just like the human brain. With RRAM forming the crux of these systems, it is critical that device characteristics for the RRAM devices are well modeled. But modeling of RRAM based synapse devices is challenging owing to the fact that the RRAM devices used need well-defined analog behavior, which is the precondition for brain like functions. Device performance under AC stress and cycle to cycle variability are factors which affect the potentiation and depression of the synapse device. Standard models reported in this work can predict the digital behavior of the RRAM, but one may implement them for the analog behavior. Though a few models are reported [200, 201] to describe the switching mechanism in analog RRAM, but it has been difficult to mathematically define them and translate it into a compact form. However, significant research is ongoing to quantitatively describe these characteristics and translate it into a compact form to be used on the circuit level.

Conclusions

In summary, the important features of all widely accepted RRAM models have been discussed. This work fulfills the requirement of the modeling community for a unified discussion on the various RRAM models. Many of the recent models, such as Stanford/ ASU model, Gonzelez-Cordero et al. model, Prodromakis model, have provided apt explanations for RRAM processes based on the early models. Implementations of different window functions like Joglekar, Biolek, and Prodromakis have been presented and compared. Various unexplained phenomena occurring in the devices are numerically validated in the models. No one model can be deemed as the perfect one, owing to the variety of materials, fabrication processes and device operations exist in the RRAM devices. Each model has been tuned accordingly to fit the device used. Researchers are still some time away from developing a generic RRAM model owing to these factors and also due to the deficiencies in the modeling techniques. Accurate and well-defined modeling techniques have been discussed in the “Well-Posed Memristive System Definitions” section, which should act as a competent template for future model development. Combined with the detailed analysis provided for past RRAM models, this review work can potentially act as a focal point for RRAM model developers.

Perovskita híbrida de vaporización secuencial para células solares de heterounión plana Síntesis en un solo recipiente de nanoplacas de Cu2ZnSnSe4 y su actividad fotocatalítica impulsada por la luz visible