Um estudo coletivo sobre modelagem e simulação de memória de acesso aleatório resistiva

Modelos RRAM para dispositivos bipolares

Modelo Chua

Modelo Linear de Deriva de Íons

Modelo de deriva de íons não linear

Exponential Ion Drift Model

Simmons Tunneling Barrier Model

Yakopcic Model

TEAM/VTEAM Model

A sinusoidal input of 1 V applied to the TEAM model using the same fitting parameters as used in Fig. 10 [76]. The values of R _ON e R _OFF are set as 50 Ω and 1 kΩ, and an ideal rectangular window function is applied in Eqs. (38) and (37). a I - V curve and b state variable. It is to be noted that the device is asymmetric, i.e., switching OFF is slower than switching ON [76]

Variation of the state variable with time is asymmetrical in nature, as shown in Fig. 8b. This means that the ON and OFF switching times are not equal. In the Eq. (36), i _em e i _desligado act as the current thresholds. Functions f _em and f _desligado are window functions which bound the internal state variable x ( t ) within [w _em , w _desligado ] Window functions are described as [76]:
$$ {f}_{\mathrm{off}}(w)=\exp \left[-\exp \left(\frac{w-{a}_{\mathrm{off}}}{w_c}\right)\right], $$ (38) $$ {f}_{\mathrm{on}}(w)=\exp \left[-\exp \left(-\frac{w-{a}_{\mathrm{on}}}{w_c}\right)\right], $$ (39)
The window functions describe the dependence of the derivative in the state variable x . They work well within the described boundaries, but the problem arises when the device goes beyond the boundaries. There are no limiting parameters here, and the window function only describes the state variable inside a particular limit. If the device goes beyond the boundaries, it can cause convergence issues with the simulator and it does not make sense for good modeling practice in case of analog devices.

I - V relationship in this model is derived from the tunneling barrier model, as discussed in the “Simmons Tunneling Barrier Model” section. Due to the non-linear nature of the tunneling current, the change in resistance varies exponentially with the state variable. So, it is assumed that any change in the tunnel barrier width changes the memristance in an exponential manner which deduces to [76]:
$$ v(t)={R}_{\mathrm{ON}}{e}^{\left(\lambda /{w}_{\mathrm{off}}-{w}_{\mathrm{on}}\right)\left(w-{w}_{\mathrm{on}}\right)}\times i(t) $$ (40)
Here, λ is a fitting parameter and R _ON the equivalent effective resistance at the bounds.

I - V relationship for this model can be seen in Fig. 8a [76]. Although there is a presence of a pinched hysteresis, the form and structure of the curve are not well-defined. The model is driven with a sinusoidal input of 1 V. The verification done for this model is different from the tunneling model [70,71,72] in terms of the platform used to simulate it. The latter model uses a SPICE macro model [72] to describe the equations, but SPICE takes up a significant amount of computation time. Modeling in Verilog-A [140,141,142,143] is much more efficient, and the TEAM model [75] utilizes this functionality to model the equations presented by them.

A slightly modified version of the TEAM model with the introduction of voltage threshold levels was reported by the same group, called Voltage Threshold Adaptive Memristor model (VTEAM) [77]. Discussed TEAM model was based on threshold currents, whereas VTEAM is based on threshold voltages. The major advantages cited for using threshold voltages is that comparison with current causes performance and reliability issues if the condition is not satisfied, i.e., a low-current threshold will automatically have a low-voltage threshold as well. This might affect the overall performance of the device. Also with a threshold voltage, there is no risk with going overboard with high power and voltage destroying the device as the values are automatically controlled.

The VTEAM follows a similar concept to the TEAM model, being based on an expression of the derivative of an internal state variable. The current is dependent on the state variable itself. The only difference is inclusion of a threshold voltage. The internal state variable (w ) is defined as [77]:
$$ \frac{dw}{dt}=\left\{\begin{array}{c}{k}_{\mathrm{off}}\times {\left(\frac{v(t)}{v_{\mathrm{off}}}-1\right)}^{\alpha_{\mathrm{off}}}\times {f}_{\mathrm{off}}(w),\kern0.5em 0<{v}_{\mathrm{off}} Similar to the TEAM model, the functions f _em and f _desligado act as window functions which bound the internal state variable w within [w _em , w _desligado ] As has been assumed in the model, current varies exponentially with the internal state variable on most occasions which is defined by [77]:
$$ i(t)=\frac{e^{-\frac{\lambda }{w_{\mathrm{off}}-{w}_{\mathrm{on}}}\times \left(w-{w}_{\mathrm{on}}\right)}}{R_{ON}}\times v(t) $$ (42)
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]. a 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]:
$$ \frac{dg}{dt}={\nu}_0\exp \left(\frac{-{E}_{a,m}}{k_bT}\right)\sinh \left(\frac{q{a}_h\gamma v}{L{k}_bT}\right) $$ (43)
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]:
$$ g\mid t+\Delta t=F\left[g|t,\frac{dg}{dt}\right]+{\delta}_g\times \overset{\sim }{X}(n)\Delta t,\kern2.25em n=\left\lfloor \frac{t}{T_{GN}}\right\rfloor $$ (44)
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]:
$$ {\delta}_g(T)=\frac{\delta_g^0}{\left\{1+\exp \left[\frac{\left({T}_{\mathrm{crit}}-T\right)}{T_{\mathrm{amb}}}\right]\right\}} $$ (45)
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\frac{\partial T}{\partial t}=v(t)i(t)-k\left(T-{T}_{\mathrm{amb}}\right) $$ (46)
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].
$$ i\left(g,v\right)={i}_0\exp \left(\frac{-g}{g_0}\right)\sinh \left(\frac{v}{v_0}\right) $$ (47)
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. I - 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].

a 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].
$$ \frac{\partial {n}_D}{\partial t}=\nabla \times \left({D}_s\nabla {n}_D-\mu v{n}_D\right)+G $$ (48)
onde 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]:
$$ {D}_s=\frac{1}{2}\times {a}^2\times {f}_e\times \exp \left(-{E}_a/{k}_{\boldsymbol{B}}T\right) $$ (49) $$ v={a}_h\times f\times \exp \left(-{E}_a/{k}_BT\right)\times \sinh \left(q{a}_hE/{k}_BT\right) $$ (50) $$ G=A\times \exp \left(-\left({E}_a-q{l}_mE\right)/{k}_BT\right) $$ (51)
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. a V _o density (n _D ) map. Calculated profiles of b n _D , c T , e 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 ²⁻ e 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]:
$$ \frac{dx}{dt}=a\times f\times \exp \left(-\frac{E_i-\gamma ZeV}{k_BT}\right) $$ (52)
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]:
$$ \frac{dx}{dt}=a\times f\times \exp \left(-\frac{E_h}{k_BT}\right)\sinh \left(\frac{a_h ZeE}{k_BT}\right) $$ (53)
The RESET process when dominated by the recombination between O ²⁻ and V_o is written as [89]:
$$ \frac{dx}{dt}=a\times f\times \exp \left(-\frac{\Delta {E}_r}{k_BT}\right) $$ (54)
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]:
$$ \frac{dx}{dt}=-a\times {f}_e\times \exp \left(-\frac{E_a-{\alpha}_a ZeE}{k_BT}\right) $$ (55)
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]:
$$ \frac{dw}{dt}=\left(\Delta w+\frac{\Delta {w}^2}{2w}\right)\times {f}_e\times \exp \left(-\frac{E_a-\gamma Zev}{k_BT}\right) $$ (56)
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]:
$$ i={i}_0\exp \left(-x/{x}_T\right)\sinh \left(v/{v}_T\right) $$ (57)
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]:
$$ {\tau}_{\mathrm{form}}={\tau}_{\mathrm{form}0}\times {e}^{\frac{E_{a\mathrm{Form}}-q\times {\alpha}_s\times {v}_{\mathrm{Cell}}}{k_b\times T}} $$ (58) $$ \frac{d{r}_{\mathrm{CFmax}}}{dx}=\frac{r_{\mathrm{work}}-{r}_{\mathrm{CFmax}}}{\tau_{\mathrm{form}}} $$ (59)
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]:
$$ {\tau}_{\mathrm{Red}}={\tau}_{\mathrm{Red}\mathrm{ox}}\times {e}^{\frac{E_a-q\times {\alpha}_s\times {V}_{\mathrm{cell}}}{k_b\times T}} $$ (60) $$ {\tau}_{Ox}={\tau}_{Redox}\times {e}^{\frac{E_a+q\times \left(1-{\alpha}_s\right)\times {V}_{\mathrm{Cell}}}{k_b\times T}} $$ (61)
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]:
$$ \frac{d{r}_{CF}}{dt}=\frac{r_{CF max}-{r}_{CF}}{\tau_{red}}-\frac{r_{CF}}{\tau_{Ox}} $$ (62)
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]:
$$ \sigma (x)\times E{(x)}^2=-k\times \frac{\partial^2T(x)}{\partial {x}^2} $$ (63) $$ T(x)={T}_{\mathrm{amb}}+\frac{v_{\mathrm{Cell}}^2}{2\times {L}_x^2\times k}\times \left(\frac{L_x^2}{4}-{x}^2\right)\times {\sigma}_{\mathrm{eq}} $$ (64) $$ T={T}_{\mathrm{amb}}+\frac{v_{\mathrm{Cell}}^2}{8\times k}\times {\sigma}_{\mathrm{eq}} $$ (65) $$ {\sigma}_{\mathrm{eq}}={\sigma}_{\mathrm{CF}}\times \frac{r_{\mathrm{CF}}^2}{r_{\mathrm{work}}^2}-{\sigma}_{Ox}\times \frac{r_{\mathrm{CF}\mathrm{max}}^2-{r}_{\mathrm{CF}}^2}{r_{\mathrm{work}}^2} $$ (66)
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].

I - 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]:
$$ {i}_{\mathrm{cell}}={i}_{\mathrm{sub}-\mathrm{oxide}}+{i}_{\mathrm{CF}}+{i}_{\mathrm{pristine}} $$ (67) $$ {i}_{\mathrm{CF}}=E\times \pi \times {\sigma}_{\mathrm{CF}}\times {r}_{\mathrm{CF}}^2 $$ (68) $$ {i}_{\mathrm{sub}-\mathrm{oxide}}=E\times \pi \times {\sigma}_{ox}\times \left({r}_{\mathrm{CF}\mathrm{max}}^2-{r}_{\mathrm{CF}}^2\right) $$ (69) $$ {i}_{\mathrm{pristine}}={S}_{\mathrm{cell}}\times {A}_e\times {E}^2\times \exp \frac{-{B}_e}{E} $$ (70) $$ {A}_e=\frac{m_e\times {q}^3}{8\pi \times h\times {m}_e^{ox}\times {\phi}_b} $$ (71)
a , 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]:
$$ \mathrm{if}\ {\phi}_b\ge q{L}_xE:{B}_e=\frac{8\pi \sqrt{2{m}_e^{ox}}}{3\times h\times q}\left[{\phi}_b^{\frac{3}{2}}-{\left({\phi}_b-q{L}_xE\right)}^{\frac{3}{2}}\right] $$$$ \mathrm{otherwise}:{B}_e=\frac{8\pi \sqrt{2{m}_e^{ox}}}{3\times h\times q}\times {\phi}_b^{3/2} $$ (72)
Here, S_cell is a section of the RRAM cell. A _e e 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]:
$$ {r}_{{\mathrm{CFmax}}_{i+1}}=\left({r}_{{\mathrm{CFmax}}_i}-{r}_{\mathrm{work}}\right)\times {e}^{\frac{-\Delta t}{\tau_{\mathrm{Form}}}}+{r}_{\mathrm{work}} $$ (73) $$ {r}_{{\mathrm{CF}}_{i+1}}=\left({r}_{{\mathrm{CF}}_i}-{r}_{{\mathrm{CF}\mathrm{max}}_i}\times \frac{\tau_{\mathrm{eq}}}{\tau_{\mathrm{Red}}}\right)\times {e}^{\frac{-\Delta t}{\tau_{\mathrm{eq}}}}+{r}_{{\mathrm{CF}\mathrm{max}}_i}\times \frac{\tau_{\mathrm{eq}}}{\tau_{\mathrm{Red}}} $$ (74) $$ \mathrm{where}\ {\tau}_{\mathrm{eq}}=\frac{\tau_{\mathrm{Red}}\times {\tau}_{\mathrm{Ox}}}{\tau_{\mathrm{Red}}+{\tau}_{\mathrm{Ox}}} $$ (75)
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]:
$$ {\sigma}_n={\sigma}_0\exp \left(-\frac{E_{ac}}{k_B{T}_n}\right) $$ (76)
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]:
$$ {P}_{g,n}=\left(1-{C}_n\right)\exp \left(-\frac{E_a-{a}_a ZeE}{k_B{T}_n}\right) $$ (77) $$ {P}_{h,n\to m}={C}_n\left(1-{C}_m\right)\exp \left(-\frac{E_h-{a}_h ZeE}{k_B{T}_n}\right) $$ (78) $$ {P}_r={C}_n\exp \left(-\frac{\Delta {E}_r}{k_B{T}_n}\right) $$ (79)
Aqui, 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 e E _h are the activation energies for oxygen species generation and hopping respectively. Similarly, a _a e 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]:
$$ {\sigma}_{mn}={f}_e\frac{e^2{c}_m\left(1-{c}_n\right)}{d_{mn}{k}_B{T}_n}\exp \left(-\frac{2{d}_{mn}}{\alpha}\right)\exp \left(-\frac{e\mid \nabla {\varphi}_{mn}\mid }{k_B{T}_n}\right) $$ (80)
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]:
$$ \frac{1}{\tau_{{\mathrm{Red}}_{\left(T,B\right)}}}={A}_{\mathrm{Red}\mathrm{ox}}\times {e}^{-\frac{E_a-q\times {\alpha}_s\times {v}_{\left(T,B\right)}}{k_B\times {T}_{{\mathrm{CF}}_{\left(T,B\right)}}}} $$ (81) $$ \frac{1}{\tau_{{\mathrm{Ox}}_{\left(T,B\right)}}}={A}_{\mathrm{Redox}}\times {e}^{-\frac{E_a-q\times \left(1-{\alpha}_s\right)\times {v}_{\left(T,B\right)}}{k_b\times {T}_{{\mathrm{CF}}_{\left(T,B\right)}}}} $$ (82)
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]:
$$ \frac{d{r}_{CF_T}}{dt}=\frac{r_{CF_{TM}}-{r}_{CF_T}}{\tau_{Red_T}}-\frac{r_{CF_T}}{\tau_{Ox_T}} $$ (83) $$ \frac{d{r}_{CF_B}}{dt}=\frac{r_{CF_T}-{r}_{CF_B}}{\tau_{Red_B}}-\frac{r_{CF_B}}{\tau_{Ox_B}} $$ (84)
CF radius is set a boundary which defines the limits for the CF growth/rupture which is given as [93]:
$$ {r}_{CF_{Tm}}\le {r}_{CF_T}\le {r}_{CF_{TM}} $$ (85) $$ {r}_{CF_B}\le {r}_{CF_{BM}} $$ (86)
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 _ox ) 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]:
$$ {i}_{\mathrm{CF}}=\frac{v_{\mathrm{app}}}{R_{\mathrm{CF}}\Big\Vert {R}_{\mathrm{ox}}} $$ (87) $$ {R}_{\mathrm{CF}}=\frac{L}{\pi \times {\sigma}_{\mathrm{CF}}\times {r}_{{\mathrm{CF}}_T}\times {r}_{{\mathrm{CF}}_B}} $$ (88) $$ {R}_{\mathrm{ox}}=\left\{\begin{array}{c}\frac{L\times \beta }{2\times {\sigma}_{\mathrm{ox}}\times \pi \times {r}_{\mathrm{CF}\mathrm{ext}}\times \left({r}_{{\mathrm{CF}}_B}-{r}_{{\mathrm{CF}}_T}\right)},\kern0.75em {r}_{{\mathrm{CF}}_T}\ne {r}_{{\mathrm{CF}}_B}\\ {}\frac{L}{\sigma_{\mathrm{ox}}\times \pi \times \left({r}_{\mathrm{CF}\mathrm{ext}}^2-{r}_{{\mathrm{CF}}_T}^2\right)},\kern4.25em {r}_{{\mathrm{CF}}_T}={r}_{{\mathrm{CF}}_B}\end{array}\right. $$ (89)
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]:
$$ {i}_{\mathrm{ox}}=\operatorname{sign}\ \left({v}_{\mathrm{app}}\right)\times {A}_{\mathrm{HRS}}\times {S}_{\mathrm{Cell}}\times {\left\{\frac{\mid {v}_{\mathrm{app}}\mid }{L}\right\}}^{\alpha_{\mathrm{HRS}}} $$ (90)
Here, A _HRS e α _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.

a 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 I - 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]:
$$ {\sigma}_{\mathrm{CF}}\times E{(x)}^2=\frac{2{c}_h}{r_{\mathrm{CF}}(x)}\left({T}_{\mathrm{CF}}(x)-{T}_{\mathrm{ox}}\right)-k\frac{\partial^2{T}_{\mathrm{CF}}(x)}{\partial {x}^2} $$ (91)
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]:
$$ {T}_{{\mathrm{CF}}_{\left(T,B\right)}}={T}_{\mathrm{amb}}+\frac{\sigma_{{\mathrm{CF}}_{\left(T,B\right)}}\times {r}_{{\mathrm{CF}}_{\left(T,B\right)}}\times {E}_{\left(T,B\right)}^2}{2\times {c}_h}\left(1-\cos h\left(\frac{\alpha_{\left(T,B\right)}L}{2}\right)\right)+\frac{d{T}_{0_{\left(T,B\right)}}}{\alpha_{\left(T,B\right)}}\sin h\left(\frac{\alpha_{\left(T,B\right)}L}{2}\right) $$ (92) $$ {\alpha}_{\left(T,B\right)}=\sqrt{\frac{2h}{k\times {r}_{{\mathrm{CF}}_{\left(T,B\right)}}}} $$ (93) $$ d{T}_{0_{\left(T,B\right)}}=\frac{\sigma_{{\mathrm{CF}}_{\left(T,B\right)}}\times {r}_{{\mathrm{CF}}_{\left(T,B\right)}}\times {\xi}_{\left(T,B\right)}^2\tan h\left(\frac{\alpha_{\left(T,B\right)}L}{2}\right)}{\sqrt{2\times k\times {c}_h\times {r}_{{\mathrm{CF}}_{\left(T,B\right)}}}} $$ (94) $$ {\sigma}_{{\mathrm{CF}}_{\left(T,B\right)}}=\frac{\sigma_{{\mathrm{CF}}_0}}{1+{\alpha}_T\left({T}_{{\mathrm{CF}}_{\left(T,B\right)}}-{T}_{\mathrm{amb}}\right)} $$ (95)
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.

a 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

Filament Dissolution Model

Bocquet Unipolar Model

Window Function Models

• 
  ➢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.
  

Joglekar Window

Biolek Window Function

Benderli-Way Window Function

Shin Window Function

Prodromakis Window Function

Boundary Condition-Based Model (BCM) Function

Model Verification

Well-Posed Memristive System Definitions

Accurate Description of Hysteresis

Proper Definition of Internal Unknown Variables 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].
  

Developing Generic RRAM Models

Improving Solution Convergence

Improving Existing Models

Novel RRAM Applications