1 Introduction

The recent progress in resistive switching devices [15, 20, 22, 28, 42, 46, 47] gives hope for adoption of this technology in various computing applications [44, 49] in the near future. The development of such applications, and in particular those utilizing analog properties of memristive devices [2, 3], will heavily rely on the availability of accurate predictive device models. Ideally, such models should describe complete I–V behavior, e.g., being able to predict current i(t 0) via device at time t 0 for an applied voltage bias v(t 0). Because resistive switching devices have memory, the current i(t 0) should also depend on the history of applied voltage bias before time t 0, or, equivalently, on the memory state variable vector w at time t 0. Such memory state variables represent certain physical parameters, which are changing upon switching the device, e.g., radius and/or length of switching filament [13, 29, 36, 50]. A very convenient method for capturing the complete I–V behavior is to use a set of two equations describing memristive system [11]. In particular, the change in memory state of a device is described as a function of applied electrical stimulus (e.g., voltage bias) and the current memory state of the device, i.e.,

$$\dot{\varvec{w}} = G\left( {v\left( t \right),\varvec{w}} \right).$$
(1)

and the static equation models current–voltage relation for a particular memory state, i.e.,

$$i = F\left( {v\left( t \right),\varvec{w}} \right)v\left( t \right).$$
(2)

While there has been an impressive progress in understanding and modeling switching behavior in memristive devices, the majority of reported models are not suitable for large-scale circuit simulations. For example, some models focus on specific aspects of the memristive behavior, e.g., static equation only [4], or a particular aspect of switching dynamics [13, 19, 21, 43], and hence are incomplete. Others are derived assuming very simple physical models [6, 18, 48] and therefore could not accurately predict experimental behavior—see also comprehensive reviews of such models in Refs. [12, 26]. Alternatively, some models are too computationally intensive, e.g., due to necessity of solving coupled differential equations [8, 16, 23–25, 27, 31, 32, 40] or running molecular dynamic simulations [9, 35, 38]. Several compact (SPICE) models, which are the most suitable for large-scale simulations, have been also very recently proposed for valence change [1, 17, 36, 39, 45] and electrochemical resistive switching devices [29, 50]. Unfortunately, models for at least the former type of devices are still not sufficiently accurate and require further improvement. Such models are typically derived by assuming a particular (typically simple) physical mechanism for resistive switching and electron transport and fitting experimental data to the equations corresponding to this mechanism. For example, in Ref. [36], both dynamic and static equations are fitted assuming modulation of the tunneling barrier width, which is certainly a simplification of the actual physical mechanism and as a result, the model does not predict accurately set transition. Because multiple mechanisms are possibly involved in resistive switching [49], such simplification may not always be adequate for accurate models. Additionally, the models are not likely to be general (e.g., judging by the diversity of reported models for the same material system even from the same authors), and future devices may require development of a completely new model from scratch. This is not a marginal issue because as the attempts are made to improve memristive devices, the I–V behavior may change significantly.

The main contribution of this paper is a development of a general approach for modeling memristive devices. The approach is mainly based on fitting experimental data and hence potentially applicable to broad class of devices featuring smooth non-abrupt switching transition. Using the proposed approach, we derive a model for specific titanium dioxide devices. The model is accurate and at the same time, simple enough to be suitable for large-scale simulations.

2 General modeling approach

The modeling approach is based on several assumptions, which simplify derivation of Eqs. 1 and 2. The first assumption is to use pulse stress for deriving a dynamic equation. The primary reason is that for a constant voltage pulse with sufficiently short duration ∆t, Eq. 1 can be written as

$$\Delta \varvec{w} \approx G(v,\varvec{w})\Delta t,$$
(3)

which simplifies derivation of G(v, w) by a fitting procedure.

The second assumption, which helps to decouple the derivation of static and dynamic equations, is that practical (nonvolatile) memristive devices have highly nonlinear kinetics [47, 49] so that it should be possible to measure I–V at relatively small biases without causing much disturbance to the memory state. The safe range of voltages depends on the particular type of devices and can be, e.g., determined by performing characterization of the switching kinetics [2, 36]. A related assumption is that the memory state is considered to be uniquely characterized by I–V measured at small non-disturbing biases (denoted as “read” biases in this paper).

Taking into account the described assumptions, the first step of dynamic equation modeling is the collection of large amounts of data by switching the device with fixed short-duration voltage pulses with different amplitudes and measuring the I–V at a non-disturbing bias after each pulse, e.g., similar to the pulse algorithms described in Refs. [2, 36]. To simplify the model, it is convenient to use as few state variables as possible, so that only a small number of measurements along non-disturbing I–V is required to characterize uniquely the internal state of the device. Ideally, this could be just one state variable, e.g., a measured resistance of the devices w ≡ R Vread at some non-disturbing bias v read. In this ideal case, the objective is to measure the change in state ΔR Vread for many different combinations of initial state R Vread and pulse amplitude v. The actual implementation details of the pulse generation algorithm are not important as long as it will cover all combinations of R Vread and v. The next step is to find G(v, R Vread) by fitting a surface to ΔR Vread (v, R Vread) data. If the ΔR Vread (v, R Vread) data are noisy and, e.g., there is a large spread of ΔR Vread for the same values of R Vread and v, then more state variables might be needed, e.g., corresponding to the measured resistance at different non-disturbing biases. In this case, separate fitting for each state variable should be performed.

The static equation is modeled by first obtaining I–V data from fast non-disturbing sweeps for the device in various initial states, and then fitting the data. For example, in case of single state variable, F(v, R Vread) is found by fitting a surface to i(v, R Vread) data, where v is within a range of voltages used for sweep experiment. Similar to the modeling of dynamic equation, more state variables must be introduced if data are noisy and, e.g., if the device has different I–Vs for the same R Vread. Note that it is important to use as large voltage range as possible without disturbing the state of the device (which can be ensured by checking that the currents for rising and falling directions of voltage sweep overlap), because nonlinear features in static I–V are typically prominent at high voltages.

3 Model for Pt/TiO2−x/Pt memristive devices

Let us now demonstrate the proposed modeling approach on the example of Pt/TiO2−x/Pt memristive devices, whose structure and fabrication methods are described in Ref. [3].Footnote 1 For such devices, a single state variable R 0.5, which represents resistance measured at non-disturbing read bias v read = 0.5 V, turns out to be sufficient for good accuracy. Following the proposed approach, the device is switched into different intermediate states by applying a sequence of positive and negative write voltage pulses with Δt = 10 µs and different amplitudes (Fig. 1a, b). Each write voltage pulse is followed by a read pulse to measure the new device state R 0.5 + ΔR 0.5. The process is repeated for sufficiently large number of different write voltage pulses and initial device states R 0.5 (with more than 15,000 measurements in total) to gather enough points for fitting procedure. Note that write voltage amplitudes were limited to >−2.5 and <5 V for set and reset transitions, respectively, to avoid breaking the device.

Fig. 1
figure 1

a Evolution of device resistance measured at 0.5 V as a result of (b) application of sequence of voltage pulses with 10 μs write pulse duration and 1 s time between pulses. c Same as a shown as a normalized 3D plot (in percent) and (d) fitted surface described by Eq. 4 with fitting parameters shown in the inset. To reduce the effect of random telegraph noise [14], the resistance measurement is averaged over 20,000 samples taken over 1 ms

Full size image

The resulting 3D plot for resistance change (Fig. 1c) is smooth and features an effective voltage threshold, which justifies using 0.5 V bias for non-disturbing read, and convergence to zero for both very high and very low resistances. These features are likely related to Joule heating-assisted resistive switching [49], e.g., super-linear dependence of temperature increase on applied voltage, redistribution of dissipated power from active region to series resistance upon decrease in resistance [7], and decrease in total dissipated power when resistance increases. Instead of relying on accurate physical model, we use fitting functions \(\sinh {{\left[ v \right]} \mathord{\left/ {\vphantom {{\left[ v \right]} {\left( {1 + \exp \left[ {\chi v + \zeta } \right]} \right)}}} \right. \kern-0pt} {\left( {1 + \exp \left[ {\chi v + \zeta } \right]} \right)}}\), \({{R_{0.5} } \mathord{\left/ {\vphantom {{R_{0.5} } {\left( {1 + \exp \left[ {\delta R_{0.5} + \theta } \right]} \right)}}} \right. \kern-0pt} {\left( {1 + \exp \left[ {\delta R_{0.5} + \theta } \right]} \right)}}\), and \(\exp \left[ {\lambda R_{0.5} } \right]\) to mimic these three features, respectively, i.e.,

$$\Delta R_{0.5} = \alpha \frac{\sinh [v]}{{1 + \exp \left[ {\chi v + \zeta } \right]}}\frac{{R_{0.5} }}{{1 + \exp \left[ {\delta R_{0.5} + \theta } \right]}} \exp \left[ {\lambda R_{0.5} } \right] \Delta t,$$
(4)

where α, λ, δ, θ are fitting parameters specific to the direction of switching (inset of Fig. 1d). Figure 1d shows the 3D surface based on least square error fitting of Eq. 4 to the experimental data. For such fitting, the first term grows super-exponentially with the voltage and emulates an effective switching threshold, the second term is super-linear with R 0.5 for R 0.5 < 10 kΩ for set switching and linear with R 0.5 for both set and reset switching in the remaining range, while the last term introduces an exponential decrease with respect to R 0.5 in the whole range of resistances for reset switching and for R 0.5 > 50 kΩ for set switching. Note that such choice of fitting function is ad hoc and primarily motivated by having as few fitting coefficients as possible.

To obtain the static model, the device is first switched into several intermediate states (represented by R 0.5) by applying positive and negative triangular sweep voltage stimuli (Fig. 2). The static portion of the experimental data i(v, R 0.5) is then fitted with the following function of v and R 0.5

$$\log_{10} |i| = g_{1} \tanh \left( {1.5\log_{10} \left| v \right|} \right) + \log_{10} |v| + g_{2} ,$$
(5)

where g 1 and g 2 are functions of R 0.5 (Fig. 3). Similar to a dynamic model derivation, the choice of fitting function for static equation is ad hoc and motivated by a tilted tangent hyperbolic shape of the static curves in the log–log scale (Fig. 2b, d). It is worth noting that for either small or large voltages, Eq. 5 simplifies to the following linear relation between i and v

$$i \approx \left\{ {\begin{array}{*{20}c} {10^{{g_{2} - g_{1} }} v,} & {\left| v \right| \ll 1} \\ {10^{{g_{2} + g_{1} }} v,} & {\left| v \right| \gg 1} \\ \end{array} } \right..$$
(6)

The physical explanation for linear behavior at small voltages is self-evident, while at high voltages, it is likely due to the dominant effect of series resistance in memristive devices [7].

Fig. 2
figure 2

Switching I–V characteristics and corresponding fitting of static I–Vs for a, b reset and c, d set operations. For clarity, figures are shown using a, c linear and b, d logarithmic scales

Full size image
Fig. 3
figure 3

Fitting for g 1 and g 2 functions for a, b positive and c, d negative voltages. In each panel, blue dots are experimental data, while red curve is fitting according to the specific formula

Full size image

Given G(v, w) from Eq. 4, the device response to an arbitrary time-varying voltage stimulus can be in principle calculated by solving differential Eq. 1. However, another approach is to approximate time-varying voltage stimulus with a sequence of corresponding fixed-duration voltage pulses and use Eq. 3 instead. Figure 4 shows simulation results for the full I–V sweep using approximated stimuli with the only input parameters to the simulation are initial (measured) state of the device R 0.5 and applied voltage stimulus. In particular, two cases are simulated and shown on Fig. 4. In the first case, only the dynamic equation is utilized (which predicts change in the resistance at 0.5 V), and linear current–voltage dependence i = v/R 0.5 is assumed to get the current at the specific applied voltage. As expected, in this case model is somewhat accurate for small voltages (and hence the static equation may not be needed for modeling) but significantly underestimates current at high voltages. On the other hand, the simulated switching I–V characteristics are in a good agreement with experimental data in the whole range of voltages if both static and dynamic equations are employed.

Fig. 4
figure 4

Experimental and simulated switching I–V characteristics using dynamic equation only and combined dynamic and static equations for triangular voltage sweep

Full size image

4 Discussion and summary

Let us now discuss some limitations and potential improvements for the proposed modeling approach. One reservation concerning using pulse stress (the first assumption of the modeling approach) is that transient effects with slow characteristic times are challenging to model. For example, such transient may be due to a relatively slow heating transient in the device and could be represented by a state variable corresponding to the internal temperature [5, 37]. In this case, the device response to a train of pulses greatly depends on an interpulse delay, even if Δt is very small. In the proposed modeling approach, slow transients are neglected assuming that there is sufficiently long time ≫Δt between applied write and read voltage pulses. Nevertheless, because voltage pulse stimulus is easy to implement in a hardware, this simplification is justified for practical applications. Another compelling reason to use pulse train stimulus with large interpulse delay is to eliminate the effect of secondary volatile switching, which is often present in metal oxide devices [5, 10, 30, 33, 34].

It is clear from Fig. 1 that there are not much meaningful data for the applied voltages below effective switching threshold and for the device states close to extreme on or off values. In these regions, the changes in R 0.5 are insignificant and typically too noisy to be used for reliable fitting. To address this issue, the modeling approach can be extended by applying variable duration pulse stress [2, 36]. Measured changes in the device state upon application of long-duration voltage pulses can be normalized with respect to ∆t duration and used reliably for the proposed fitting approach. Additionally, subthreshold behavior might be modeled in ad hoc fashion by introducing an additional term in dynamic equation, e.g., based on activation energy of switching [41]. For example, the model predicts that the switching rate for set transition changes by as little as a factor of 900 when applied voltage is reduced from 1 to 0.5 V, which is not appropriate for true nonvolatile memory. By multiplying right hand side of Eq. 4 by 1/(1 + exp[(v 1 − |v|)/v 2]), the ratio of switching rates at the same applied voltages is increased to >1016 when using v 1 = 0.75 and v 2 = 0.008, which is expected retention performance for the considered devices [3]. At the same time, it is easy to check that such extra term will have no effect on the simulation shown on Fig. 4.

More generally, the switching model can be approximated with an equivalent circuit shown on Fig. 5. It is natural to expect that the maximum resistance of the memristive device is limited by the leakage through the film, which is modeled with resistor in parallel to the active part of the device. Its minimum resistance is determined by a resistance (e.g., corresponding to a filament) connected in series with an active part of the device, which is described by memristive equations. The advantage of such equivalent circuit is that it naturally bounds the device resistance and provides physically plausible switching behavior of the device near extreme on and off states. Additionally, the equivalent circuit can be extended to include secondary volatile switching behavior [5, 10, 30, 33, 34] and electronic noise [14]. Finally, a practical way to include switching variations is to modify parameter ζ in the model, for example, by adding a zero-mean Gaussian random variable to mimic device-to-device and cycle-to-cycle variations.

Fig. 5
figure 5

Generalized equivalent circuit for modeling metal oxide memristive devices

Full size image

In summary, this paper outlines a general approach for deriving memristive equations for resistive switching devices. The approach is purely phenomenological and is based on fitting of the experimental data and hence can be applied to a broad class of memristive devices. The knowledge of the switching mechanisms and electron transport can be helpful for finding the best fitting functions; however, it is not a requirement, which further simplifies modeling approach. The proposed approach is tested on a particular metal oxide device by comparing simulated and experimental I–Vs for a full sweep. The model shows good accuracy and at the same time, because of explicit form of equations, is computationally inexpensive, which makes it suitable for simulation of large-scale memristive circuits.