A Simple, Robust, and Versatile MATLAB Formulation of the Dynamic Memdiode Model for Bipolar-Type Resistive Random Access Memory Devices

Abstract

Modeling in an emerging technology like RRAM devices is one of the pivotal concerns for its development. In the current bibliography, most of the models face difficulties in implementing or simulating unconventional scenarios, particularly when dealing with complex input signals. In addition, circuit simulators like Spice require long running times for high-resolution results because of their internal mathematical implementation. In this work, a fast, simple, robust, and versatile model for RRAM devices built in MATLAB is presented. The proposed model is a recursive and discretized version of the dynamic memdiode model (DMM) for bipolar-type resistive switching devices originally implemented in LTspice. The DMM model basically consists of two coupled equations: one for the current (non-linear current generator) and a second one for the memory state of the device (time-dependent differential equation). This work presents an easy-to-use tool for researchers to reproduce the experimental behavior of their devices and predict the outcome from non-trivial experiments. Three study cases are reported, aimed at capturing different phenomenologies: a frequency effect study, a cycle-to-cycle variability fit, and a stochastic resonance impact analysis.

1. Introduction

In recent years, the scientific community in the microelectronics field has expressed huge interest in RRAM (Resistive Random Access Memory) devices due to their unique properties, such as low power consumption, high scalability, endurance, and CMOS compatibility, among others [1]. As the switching material, these devices can be based on metal oxides like TiO₂ or HfO₂, chalcogenide glasses, lower dimensional materials like nanowires or nanodots, as well as other novel materials like metal–oxide frameworks [2,3,4]. Their application range includes logic circuits, neuromorphic systems, cryptology, and information storage [5,6,7,8,9]. Since RRAM devices are still an emerging technology, device-level characterization, and powerful modeling tools are key for enhancing its evolution to massive commercialization. To date, a lot of models have been implemented in simulators like Spice or Verilog [10]. These simulators are useful for complex architectures involving several devices but also lead to long waiting times for running the simulations and post-processing the data. In this work, a simple and effective tool is presented for reproducing the behavior of single RRAM devices in MATLAB, a common and well-known programming language. The advantages of the tool presented consist of easy implementation (since the script is openly included in the manuscript), straightforward data analysis, short simulating running times, and versatility in terms of parameters definition and input possibilities (applying non-trivial signals like noise). In the current literature, a few works using MATLAB for RRAM modeling are discussed. In [11], different existing models were studied using MATLAB simulations. In [12], a Kinetic Monte Carlo simulation in MATLAB to estimate and optimize the reliability of RRAM devices was reported. In [13], the implementation of a memristor-based neural network using MATLAB, Simulink, and LTspice for simulations and analysis was presented. In this work, the presented tool consists of a recursive version of the dynamic memdiode model (DMM) for RRAM devices, which was originally built for LTspice [14]. Its main features and the corresponding script are introduced in Section 2. In Section 3, a study about the effect of the input signal frequency is shown. Section 4 discusses a method for incorporating cycle-to-cycle variability in the model. Finally, a stochastic resonance impact study is presented in Section 5.

2. Dynamic Memdiode Model

This section introduces the DMM, which is defined by means of mainly two coupled equations, one for the current conduction and another for the memory state. The DMM is a compact behavioral model built under the idea that the current flows through a conductive filament (CF) running across the insulating layer of a metal–insulator–metal/semiconductor structure. The CF is created and destroyed as a consequence of the application of an external field, and this is the origin of the switching behavior [15,16]. For a bipolar-type resistive switching device, the DMM reads as follows:

$$I\left(V\right)=I_0\left(\lambda \right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\left\{\alpha \left(\lambda \right)\left[V-\left(R_i+R\left(\lambda \right)\right)I\right]\right\}$$

(1)

where $I_0\left(\lambda \right)={(I}_{on}-I_{off})\lambda +I_{off}$ is the current amplitude, and $\lambda$ is the memory state (l in the model script), which ranges from 0 to 1. $\alpha \left(\lambda \right)$ (a in the model script) and $R\left(\lambda \right)$ are defined in the same way as $I_0$. ‘off’ and ‘on’ terms are linked to the high and low-resistance states (HRS and LRS), respectively. $V$ is the input voltage, and $R_i$ is a constant series resistance [14]. After some considerations thoroughly discussed in [17], the discrete recursive approach for the I-V characteristic reads

$$I\left(V\right)\approx \frac{1}{\alpha R}\{W\left[\frac{\alpha R{I}_0}{2}exp\left(\alpha V_i\right)\right]-W\left[\frac{\alpha R{I}_0}{2}exp\left(-\alpha V_i\right)\right]\}$$

(2)

where W is the Lambert function (lambertw in MATLAB), and $V_i=V-I·R_i$. Next, the evolution of the memory state in the DMM is written as the following balance equation:

$$\frac{d\lambda }{dt}=\frac{1-\lambda }{{\tau }_S}-\frac{\lambda }{{\tau }_R}$$

(3)

(3) is the general equation for a redox process and is a nondimensional differential equation. ${\tau }_{S,R}$ are the characteristic switching times for the set and reset transitions, respectively. According to our convention, the set occurs at positive bias and the reset at negative bias. Again, after some approaches discussed in [17], the recursively expressed and discretized version of (3) reads

$${\lambda }_{t+1}=\left[{\lambda }_t-H\left(V_{it}\right)\right]\exp \left\{-\frac{\mathsf{\Delta }t}{{\tau }_{iS,R}}\right\}+H\left(V_{it}\right)$$

(4)

${\lambda }_t$ is the memory state, $\mathsf{\Delta }t$ is the selected timestep, and $H$ is the Heaviside step function.

$${\tau }_{iS}\left(V_{it}\right)=\mathrm{e}\mathrm{x}\mathrm{p}[-{\eta }_S\left(V_{it}-V_S\right)]$$

(5)

$${\tau }_{iR}\left(V_{it}\right)=exp[-{\eta }_R{\lambda }^{\gamma }\left(V_{it}-V_R\right)]$$

(6)

where ${\eta }_{S,R}$ and $V_{S,R}$ (vs and vr in the model script) are the transition rates (${\eta }_S>0, {\eta }_R<0$, etas and etar in the model script, respectively) and reference switching voltages ($V_S>0, V_R<0$). Physically, these dependences are linked to the movement of metal ions or oxygen vacancies. $\gamma$ or gam, in the model script, controls the reset transition rate. The model is ruled by Equations (2) and (4). Equations (5) and (6) are auxiliary functions. The complete script is reported in Table 1 for the reader’s possibility to copy and use it for its own purposes. In the model script, the first section (in black) builds the input signal. The second section (in green) is dedicated to the model parameters and initial conditions. The third section of the script (in blue) presents the auxiliary functions required by the model equations. The fourth section (in red) calculates the current I(t + 1) and the memory state l(t + 1). Finally, (in violet) the simulated currents and voltages are stored in the Im and Vm variables and plotted in the last two lines of the script. In addition, the memory state is stored in the lm variable. The same color pattern is followed in Appendix A, Appendix B and Appendix C, where the scripts for the case studies in the following Section 3, Section 4 and Section 5 are presented. In the Appendix A, Appendix B and Appendix C, the violet fragment includes the data analysis. Figure 1 shows the model results for the script described in Table 1, showing the typical memristor current–voltage loop (in black) and the memory state–voltage loop (in blue), highlighting the high/low-resistance states and the set/reset processes.

In Figure 2, a study of the impact of modifying different model parameters on the current–voltage characteristic is presented: (a) for vs, (b) for vr, (c) for aon, (d) for aoff, (e) for etas and etar, and (f) for gam. This study shows the model’s sensitivity and is a useful guide for potential users of the model.

3. First Case Study: Modeling of Frequency and Ramp Rate Effects

This section is dedicated to illustrating the model’s sensibility to variations in the input signal frequency (for sinusoidal signals) or the ramp rate (for voltage ramps). In the literature, it is reported that for high frequencies or ramp rates, both the set and reset voltages of bipolar RS devices increase in absolute value, making a wider butterfly-shaped I-V curve. A linear relationship between the logarithm of the ramp rate and the set and reset voltages has been reported [18,19]. For ramped input voltages, the magnitude for expressing the timing is the ramp rate (RR), which is directly related to the signal frequency. Figure 3 presents the evolution of the conduction characteristic for a typical experiment as reported in [19], using a voltage ramp as the input with different RRs of 50 V/s (red), 500 V/s (orange), 5000 V/s (light green), and 50,000 V/s (dark green). On one hand, in Figure 3a, the study is focused on the set evolution when the reset ramp rate is kept fixed. On the other hand, Figure 3b focuses on the reset evolution when the ramp rate for the set event is kept unaltered. Figure 3c shows the evolution of the set and reset voltages as a function of the ramp rate for these cases. Notice that the experimental curves are affected by a series transistor (current limitation during the set). Further details about the devices and the measurement setup can be found in [19]. Simulations are reported in Figure 4, and they do not include the transistor effect. Here, the focus is on the frequency phenomenology. The dependence of the set and reset voltages as a function of the ramp rate is shown in Figure 4a, including I-V curves using as the input signal a voltage ramp with identical ramp rates and colors used in Figure 3a. In Figure 4b, the dependence of the set and reset voltages as a function of the signal frequency is illustrated, including I-V curves using a sinusoidal voltage with frequencies of 1 Hz (red), 10 Hz (orange), 100 Hz (light green), and 1000 Hz (dark green) as the input signal. The insets in Figure 4a,b show the set and reset voltage evolution as a function of the logarithm of the ramp rate (Figure 4a) or the logarithm of the frequency (Figure 4b). A linear dependence is found in both cases. The required modifications in the model script to achieve these results can be found in Appendix A. Notice that some parts are omitted, identified with %Ramp Rate or %Frequency at the end of the line, since in one script, the two kinds of evaluations (RR or frequency) can be performed. Appendix A corresponds to the frequency study. For the RR study, the lines ending with %Frequency are deleted, and the lines ending with %Ramp Rate are activated. The model parameters RR and $\mathsf{\Delta }t$ (At in the model script), included in Appendix A, are modified accordingly to set the desired input signal ramp rate or frequency.

4. Second Case Study: Modeling Cycle-to-Cycle Variability

In RRAM devices, the most significant challenge relates to its inherent C2C variability linked to atomic scale changes in the CF across the dielectric layer [20,21]. The inclusion of C2C variability is of utmost importance for any comprehensive model intended to reproduce the real device behavior. In the literature, different approaches for including variability in compact modeling can be found [22,23,24]. This section reproduces the experimental results reported in [25], where uncorrelated C2C variability is included in a previous version of the memdiode model. The experimental data come from 450 cycles of I-V curves measured in HfO₂-based memristors. Additional information about the electrical characterization and fabrication process of the devices can be found in [26]. Before the incorporation of the C2C variability in the model, a study of the experimental features of the I-V curves involving the fitdistrplus package for the R language was carried out [27]. The best distributions for each model parameter were found and incorporated into the model parameters definition. Appendix B presents the MATLAB model script used for the simulations discussed in this Section. Here, two different examples of parameters, including variability, are reported:

$$Maoff=2.1; Saoff=0.13; aoff=Maoff+Saoff\ast randn;$$

(7)

$$\begin{gathered} MIoff=3.7E-5; IoffLN=log\left(MIoff\right); SIoff=0.06; \\ Ioff=exp\left(IoffLN+SIoff\ast randn\right); \end{gathered}$$

(8)

In (7) and (8) (see Appendix B), two cases are presented, one following a normal distribution (7) and one following a lognormal distribution (8). The variables defined in (7), Maoff and Saoff (which represent the mean and standard deviation, respectively), are used to define the normally distributed aoff parameter. For a lognormally distributed parameter, like in (8), the same Mioff and Sioff are defined, but two extra steps are required: IoffLN, which is the logarithm of Mioff and the application of an exponential function in order to generate the lognormal distribution of Ioff. The randn function is used to generate a different random number in every cycle. Notice that not all the model parameters include variability, and some of them are constant for the sake of simplicity.

In Figure 5, the experimental data reported in [25] in blue (Figure 5a) and the simulations obtained using the recursive DMM model in MATLAB in red (Figure 5b) are compared. The good agreement between experimental and simulated results indicates the ability of the model to reproduce variability in the I-V loops. In addition, Figure 6a–d present histograms for the HRS current, LRS current, set voltage, and reset voltage, respectively, both for the experiments and simulations. All the histograms are well reproduced except the LRS current case, in which the experimental data exhibit a double peak likely caused by different atomic configurations of the CF. In [25], several trials were performed until finding the optimum combination of model parameters for achieving the observed variability. This process is time-consuming when performed in LTspice since the obtained results must be analyzed somewhere else. Now, the parameters, model equations, and data analysis are carried out with the same MATLAB script, largely reducing the time spent, thus easing the possibility of finding better-fitting results.

5. Third Case Study: Modeling the Stochastic Resonance Impact on RRAM Devices

Very often in electronics, noise is considered an undesired factor that needs to be eliminated. However, for certain non-linear systems, noise can play a beneficial role in terms of device performance [28]. This phenomenon is referred to as stochastic resonance (SR) and is present in a plethora of research fields like biology, physics, engineering, and many more [29,30]. Since RRAM devices are non-linear components, some works in the recent literature have studied the inclusion of external noise as a beneficial aspect of their performance [31,32,33,34,35]. In this section, noise is added to the applied voltage ramp. Again, the randn function in MATLAB is used to generate a normally distributed input voltage signal. The model script used for this study is presented in Appendix C, where the main difference with the basic model script shown in Table 1 is the following line:

$$ampn=\left(\right(\mathrm{g}-1)\ast 40)\ast 1\mathrm{E}-3;\mathrm{V}\left(\mathrm{t}\right)=\mathrm{V}\left(\mathrm{t}\right)+ampn\ast \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{n};$$

(9)

In (9), ampn controls the noise standard deviation (σ) added to the input signal V by means of the randn MATLAB function; g is a variable controlling the noise σ (every 200 cycles, it increases 40 mV). The simulation process follows the experimental method reported in [36], consisting of 200 I-V cycles adding different noise values to the applied voltage ramps. This is carried out with the aim of capturing the real device phenomenology under noisy signals. Appendix C presents the script for the simulations described in this section. Figure 7a shows the simulation results without added noise. The figure of merit that is analyzed here to determine the noise impact on RRAM devices is the separation between HRS and LRS, in particular, the resistance ratio extracted at a fixed voltage ($V=-0.3 \mathrm{V}$) as illustrated in Figure 7a. Figure 7b,c show the time evolution of the applied voltage and the simulated current, respectively, for a single simulated cycle using a noise amplitude σ = 200 mV. The resistance ratio evolution as a function of σ is illustrated in Figure 8. The experimental results extracted from [36] are shown in Figure 8a and represent the mean resistance ratio of 200 I-V cycles for different σ values. The mean resistance ratio of 200 simulated I-V cycles against noise σ is presented in Figure 8b, where an evident maximum is found for σ = 240 mV. Even though the experimental and simulated curves noticeably differ (a deep investigation of all the model parameter values is required), the model results qualitatively capture the device behavior when an external noise source is added to the input signal. The experimental simulated coherence in this case study could be enhanced by modifying the model equations to consider the impact of external noise.

6. Conclusions

In this work, a versatile and useful tool for modeling RRAM devices is presented. The model is a recursive version of the dynamic memdiode model (DMM) implemented in MATLAB software. This work includes the model script for the ease and free usage of the reader. This model has proven to be frequency sensitive, to reproduce the cycle-to-cycle variability of experimental devices, and to predict the impact of an external noise source in RRAM devices, where the resistance ratio against the noise amplitude follows the typical stochastic resonance curve. Compared with the DMM version implemented in LTspice, the time needed to modify parameters, simulate, and analyze the results is largely reduced. Because of these reasons, this tool could be of interest to researchers working in the RRAM field at the device level.