Event-Driven Stochastic Compact Model for Resistive Switching Devices

A stochastic compact model for resistive switching (RS) devices is presented. The motivation is twofold: first, introducing variability in a natural way, and second, accounting for the discrete jumps of conductance observed during set and reset transitions. The model is based on an event generation rate, and it is an “on-the-fly” procedure because events are randomly generated as the simulation proceeds in time. For the generation of events, we assume a mixed nonhomogeneous Poisson process (NHPP). Before considering RS, we deal with the generation of successive breakdown (BD) events in metal-insulator-semiconductor structures. This confirms the validity of the approach by comparing it with experimental data in which discrete events are evident. To deal with RS, we transform a previous compact model into a stochastic model. Comparison with experiments in TiN/Ti/HfO2/W devices shows the validity of the approach. Current-voltage loops and potentiation-depression transients in pulsed experiments are captured with a single set of parameters. Moreover, the model is an adequate framework to deal with both cycle-to-cycle and device-to-device variability.

Although there are several types of memristive devices, we focus on Resistive Random Access Memories (RRAMs) which involve the transport of oxygen/metal ions forming a conducting filament (CF) formed across a metal-insulatormetal (MIM) structure [2].These devices have a large number of applications including nonvolatile embedded memories in microprocessors [3], brain-inspired computing for the implementation of deep neural networks and spiking neural networks (SNNs) [4], secured communication and chip authentication using true random number generators and physically unclonable functions [5], in-memory computing [1], and field-programable gate arrays [6], among others.The present status and future prospects of industrial applications of RRAM have been recently reviewed by Wang et al. [7].
In RRAMs, the CF is initially generated (electroforming) by inducing a soft breakdown (BD) event in the insulating layer which can be partially destroyed (reset) and reformed (set) by applying voltage (or current) signals.In this way, the device can change from a low resistance state (LRS) to a high resistance state (HRS), and vice versa.Controlling the applied voltage also allows programming intermediate resistance states.The set and reset transitions are intrinsically stochastic because they are related to the random addition or extraction of defects to and from the CF.This suggests that modeling of CF-based memristors might require the consideration of random discrete changes of conductance.Experimental observations such as conductance jumps of the order of G 0 = 2q 2 / h [8] and large amplitude random telegraph noise [9] further confirm the convenience of this stochastic approach.
In the following, we present a method for the random generation of conductance jump events.As a first check, we consider the generation of multiple BD events in MIM structures.Then, we develop a behavioral model for the stochastic simulation of set and reset transitions in resistive switching (RS) devices.

II. ON-THE-FLY METHOD FOR RANDOM EVENT GENERATION
Assume that a device subjected to voltage stress shows a series of events (conductance jumps) occurring at random times.If the average number of events at time t is n 0 (t), the defect generation rate is λ (t) = dn 0 (t)/dt.
It can be shown that the cumulative distribution of the time, t, to the first successive event after t is where F t is the probability for an event to occur between t and t + t.If during the simulation we consider a sufficiently small time interval t so that λ (t) remains practically constant during t, then F t ( t) ≈ 1 − exp(−λ (t) t).Inverting this equation and generating a random number, u, uniformly distributed between 0 and 1, we can generate a random time for the first event after time t as t u = −Ln(u)/λ (t).If t u < t, an event is generated at time t.Otherwise, the event is rejected.This is an "on-the-fly" method because events are generated as the simulation proceeds in time [10].

III. STOCHASTIC MODEL FOR MULTIPLE SUCCESSIVE BD EVENTS IN OXIDE LAYERS
If a MIM structure is subjected to a constant (high) voltage stress, several successive BD events can be induced, as shown in Fig. 1 [11], [12].For modeling these results with the "on-the-fly" method of Section II, we need a model for the generation rate, λ (t).First, we consider a nonhomogeneous Poisson process (NHPP) based on the Weibull distribution, which is a standard in the modeling of the BD statistics [13].
In the Weibull model, the evolution of the average number of events is given by n 0 (t) = (t/τ ) β , where β is the shape factor and τ , the characteristic scale time.Thus, the generation rate is λ (t) ≡ dn 0 /dt = (β/τ )(t/τ ) β−1 .
In addition, in a NHPP, the cumulative distribution for the K th event is given by the following equation [13]: To check the validity of the procedure, we have simulated successive BD distributions with the "on-the-fly" method and compared the results with the analytical results of (2).Fig. 1(a) shows that the agreement is perfect, thus confirming the procedure.

A. Mixture Poisson Model
To compare the "on-the-fly" method with experimental data, we have analyzed the evolution of the number of BD events in Al 2 O 3 /HfO 2 /Al 2 O 3 /HfO 2 /Al 2 O 3 nanolaminate oxide structures (∼2 nm per layer) which show multiple and well-defined BD events in each sample.The dielectric films were grown by atomic-layer deposition (ALD) on a p + −Si (100) wafer and the top electrode is Al.Details of the fabrication process can be found in [14].
The evolution of the conductance when 40 devices are subjected to a −7 V constant voltage stress is shown in Fig. 1(b).In each sample, the BD events are identified by conductance jumps of ∼1.7G 0 .Consistently, in the "on-thefly" simulation, we will assume that each BD event causes a jump of 1.7G 0 .(a) Successive BD distributions: "on-the-fly" method (dots) and analytical result for the Weibull distribution.(b) Multiple BD events for different samples subjected to a constant voltage stress: experiment (yellow) versus "on-the-fly" method results (black).(c) Successive BD distributions (TDCM): experiment (dots) versus "on-the-fly" method (lines).(d) Evolution of the number of events under the application of a ramped voltage for different voltage acceleration factors by using the on-the-fly method.
The Weibull model of the previous section is not able to provide a good fitting of the experimental results.In particular, the dispersion of the n(t) trajectories are much smaller than experimentally found.For this reason, a generalization of the approach is required.
In previous work, we established that the statistical distribution of successive BD events in these investigated samples is consistent with the Time-Dependent Clustering Model (TDCM) [11].This is a particular type of a Mixture Poisson Model (MPM) [15] which is required when, in addition to the intrinsic randomness of time events in each sample, there are sample-to-sample variations that introduce an extra uncertainty.In this case, we need to assume that there is a distribution of the average density of defects in the different samples, with probability density f (n).The case of the TDCM is obtained by using the gamma distribution, f (n) = γ (n, α, n 0 /α), with α being a shape parameter and n 0 = (t/τ ) β , the average number of events at time t.In this case, the probability of having K events is calculated by averaging the Poisson result with the gamma distribution.The cumulative distribution of the K th event is given by the following equation [11]: To implement the "on-the-fly" method we need to consider both the gamma and Weibull distributions.First, for each sample, we select a random value of the density of defects.This is equivalent to choosing a random value of the characteristic time τ .This time can be calculated by inverting the cumulative gamma distribution (n, α, n 0 /α) and generating a random Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
number, u ′ , so that Then, the "on-the-fly" method proceeds exactly as previously described but with a different time scale τ (u ′ ) for each sample.This is equivalent to a weighted average of Weibull distributions.
Following this generalized "on-the-fly" method, we have simulated the event-time trajectories, n(t), which nicely compares with experiments as shown in Fig. 1(b).As shown in Fig. 1(c), the successive BD statistics are also in good agreement with the experimental data and with (3).
Up to this point, we considered the case of constant-voltage stresses.However, this can be generalized to arbitrary voltage signals, V (t).Although several voltage dependence models have been proposed in the literature, we consider τ (V ) = τ 0 exp(−γ V ), i.e., the so-called E-model [11], where γ is the acceleration parameter and τ 0 a constant.The generalization of the "on-the-fly" method is straightforward since the generation rate at time t is calculated according to τ = τ (V (t)).As an example, Fig. 1(d) shows the evolution of conductance for a ramped signal, V (t) = Rt, where R is the ramp rate and different values of the acceleration factor, γ .

IV. EVENT-DRIVEN MODEL OF RS A. Continuous RS Model
To implement the event-driven model, we depart from a continuous compact model which has been demonstrated to reproduce the device's electrical properties [16] and to simulate large realistic neuromorphic crossbar circuits for pattern recognition [17].
This continuous model is based on two equations.One for the memory (the internal state variable) and one for the current.For the memory equation, a first-order balance equation was proposed where is the memory variable (0 < < 1), and τ S y τ R are the set and reset characteristic times, respectively.Since the set transition resembles the BD process [18], [19], an exponential dependence for τ S was assumed where γ S is the acceleration factor, τ S0 the time scale prefactor, and R S the series resistance.On the other hand, the reset transition is controlled both by ion drift and out diffusion to the CF surroundings.Thus, in general, one has to consider both acceleration by the applied voltage and by the local temperature in the CF.This temperature is related to the dissipated power, P = I (V − I R S ).Thus, for the sake of generality, and assuming an Arrhenius temperature dependence, the characteristic reset time, τ R , can be described by the following equation: where τ R0 is the reset scale prefactor, γ R the reset voltage acceleration factor, E a the activation energy, and R TH the thermal resistance.Notice that either E a = 0 or R TH = 0 corresponds to the case of voltage acceleration exclusively while γ R = 0 describes a pure thermal process.This equation also incorporates the dependence on the ambient temperature, T but, since have not been able to check its validity experimentally, we have kept the Ahrrenius model as a first-order approximation.
The thermal resistance has been described in the literature in terms of two parallel paths for heat evacuation [20].The longitudinal thermal resistance, R L , corresponding to heat transport along the channel and the transversal resistance, R T , to heat transport toward the surrounding material.The latter is considered to be independent of the filament size to the first order, while R L is inversely proportional to the filament area, represented here by n ch that will be defined later on.Thus, R L = K L /n ch , where K L is a constant.The total thermal resistance is given by the parallel combination of R L and R T , so that The other equation of the model is the current equation This current equation is based on the Landauer approach to electron transport through narrow constrictions [21].In this approach, it is assumed that the CF is narrow enough so as to behave as a quasi-1-D nanowire.Under these conditions, conduction is described by n ch quantized channels, and the I(V) characteristic corresponds to the LRS, and it is rougly linear.Each channel contributes with ∼ G 0 to the total conductance.When the CF is broken, a potential barrier (often identified as a gap in the literature) remains, and transport is assumed to be related to tunneling through this barrier.If the voltage drops symmetrically at the two ends of the constriction, the I -V curve follows a hyperbolic sine dependence with the applied voltage.For a more detailed description of our quantum-mechanical approach leading to the heuristic current equation presented here, refer to our previous publication [22].

B. Stochastic RS Model
A natural approach to the description of RS in terms of the random generation of events is to consider small conductance jumps corresponding to the creation or destruction of one channel.For simplicity, we assume that each event increases or decreases the conductance by the same amount ( G = G 0 ) during set and reset, respectively.However, this can be generalized by considering an arbitrary value G (a model parameter) and even by introducing random changes of G.However, for the sake of simplicity, we stick to the assumption G = G 0 .The creation/destruction of single channels will occur at random times during the application of the external voltage.In this picture, the internal state of the device is n ch , which is considered as the new memory parameter.Thus, for the event-driven description of the set and reset transitions, we will first derive the memory equation by relating n ch to as: n ch = (n max −n min ) +n min , where n max y n min correspond to the maximum ( = 1) and minimum ( = 0) number of channels.These parameters play the same role as I max and I min in the continuous approach.In this work, we assume n min = 0 (full reset) and n max is left as a free parameter representing the area of the CF created during electroforming.Actually, n min has been introduced for generality (to establish a one-to-one relation with the continuous model) but allowing n min ̸ = 0 is not required because if the applied voltage is large enough reset will always be complete.That is the reason to keep n min = 0 in all our calculations.From ( 5), it follows: The two terms of the RHS represent the set and reset event generation rates which fully dominate for positive and negative bias, respectively.This is due to the strong opposite exponential dependence of τ S and τ R on voltage.On the other hand, since n min ≤ n ch ≤ n max , these generation rates are positive and negative for set and reset, so that n ch increases and decreases, respectively.
As far as the current is concerned, we have considered where I B and η are scale and shape parameters.The first term corresponds to the conduction though the n ch channels, and the second to the background tunneling regime, i.e., when the filament has a gap.Although the functional form of this second term has been successfully used to model the current in the HRS in different types of devices, other mechanisms such as trap-assisted tunneling, Poole-Frenkel emission, or hopping have been invoked in the literature.To deal with these cases, the model can be easily adapted by changing the functional dependence on the voltage of the second term.Finally, notice that n ch couples the current and memory equations, as does in the continuous approach.
For the generation of events, we follow the "on-the-fly" method.If the number of events is n(t), the event generation rate is λ (t) = dn(t)/dt.During the set transition, n ch = n min + n(t) so that λ (t) = dn ch /dt and during reset n ch = n max −n(t), so that λ (t) = −dn ch /dt.Thus, the event generation rates can be obtained from the RHS of (8).

C. Comparison With Experiment
TiN/Ti/HfO 2 /W structures (5 × 5 µm 2 ) in a cross-point configuration have been stressed under positive and negative voltage ramps as well as with trains of pulses.
The devices were fabricated on silicon wafers with a thermally grown 200 nm thick SiO 2 layer.A 200 nm-thick W back electrode was deposited by magnetron sputtering and patterned by photolithography and dry etching.Subsequently, a 10 nm thick HfO 2 layer was deposited by atomic layer deposition at 225 o C using TDMAH and H 2 O as precursors and N 2 as a carrier and purge gas.Next, the top electrode, consisting of a metal stack of 200 nm-thick TiN on a 10 nmthick Ti, was deposited by magnetron sputtering and patterned by photolithography and dry etching [23].

TABLE I MODEL PARAMETERS EXTRACTED FROM EXPERIMENTAL DATA
First, we consider the RS characteristics obtained by applying dc double sweeps for the set and reset processes with ramp rates of 0.125 V/s and 0.21 V/s, and voltage limits of + 1.1 V and −1.4 V, respectively.In Ti/HfO 2 -based memristors, both transitions are usually found to be gradual (Fig. 2), though in some particular cases, the set transition shows a rather abrupt first current jump.
We have found that it is possible to fit the experimental results of the reset transition considering only temperature acceleration (temperature is related to power and power to voltage).Thus, in this case, we have considered γ R = 0 to reduce the number of fitting parameters.
In Fig. 2, we compare the model and experiment for the RS switching loop.The extracted parameters are those of Table I.
The model parameters have been allowed to freely vary to provide the best fit of the experimental set-reset loop except for the activation energy, which has been kept fixed (E act = 1.2 eV) which is well in the range of values reported in the literature [20], [24].The final thermal resistance R T H depends on the number of channels and ranges from 20 to 70 K/mW.This is perfectly compatible with what is reported in the literature [20], [25].On the other hand, it is difficult to judge whether the set voltage acceleration factor γ S = 21V −1 is reasonable because the values reported in the literature show large dispersion [19], [26].The parameters that model the background current are similar for the set and the reset transitions indicating a rather symmetric conduction for positive and negative polarities in the HRS.
To further compare with experiments, we have considered the effect of applying a current compliance limit I comp during the set transition.I comp limits the size of the reformed filament and hence the number of channels.Without changing any of the parameters of Table I, we have simulated the effect of I comp to compare with experiments as shown in Fig. 3.Although the agreement is not quantitative, the general trends are nicely captured.In particular, the reduction of the reset peak and the displacement of this transition toward lower voltages.
Maybe the most important problems that limit the practical application of RS devices are Cycle-to-Cycle (C2C) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and Device-to-Device (D2D) variability.Fig. 4(a) shows the experimental variability of ten I(V) loops.Our goal here is to show how the model deals with both types of variability, rather than obtaining a quantitative agreement with experiments.Being based on the random generation of current jumps, our model intrinsically provides a certain degree of C2C variability.To reveal this variability, we simulate different I(V) loops keeping the same model parameters (gray lines).The random variables used in the "on-the-fly" procedure change from cycle to cycle and this gives rise to the variability shown in Fig. 4(b).The final results are qualitatively compatible with the experimental data shown in Fig. 4(a).Similarly, the model can also deal with D2D variability.To do this, we assume (as we did in Section III-A for multiple BDs) that each device has different characteristic times τ S and τ R which are randomly generated according to (4).Considering β = 1 is required to be consistent with the generation rates of (8).As an example, we have arbitrarily chosen α = 2 (the shape factor of the gamma distribution).As shown in Fig. 4(c), the variability of the I(V) loops increases with respect to considering only C2C variations.Unfortunately, no statistically significant D2D variability experimental results are available.
Up to this point, we have considered the simulation of devices subjected to ramped voltages.Now, we focus on the device behavior under the application of trains of pulses.
The experiments consist of applying 1000 pulses of 10 µs and different voltage amplitudes.In between each pair of pulses, a low-voltage pulse (−0.1 V) is used to measure the conductance in the linear regime.The separation between pulses was also 10 µs.Starting from the HRS, positive voltage pulses cause the progressive set of the device (potentiation), while departing from the LRS, negative voltage pulses progressively induce the reset transition (depression).Both the experimental and model transients are shown in Fig. 5 with excellent agreement.It is remarkable that the simulation of these transients has been performed with the same set of parameters extracted from the fitting of I(V) loops, i.e., those of Table I.These results demonstrate that the event-driven model is able to deal with any type of voltage signal and that the results are consistent for different experiments.

V. CONCLUSION
A stochastic compact model for the set and reset transitions in RS devices consistent with quantum observations has been presented for the first time.The model is a discrete version of a previous continuous behavioral model and is based on the random generation of conductance jumps (events) during the application of voltage signals of arbitrary shape.The implementation of the model follows an "on-the-fly" method in which the events are successively and randomly generated as the simulation time evolves.This method appears to be adequate to implement the model SPICE or Verilog-A and to consider discrete conductance levels as those revealing quantization in terms of G 0 .However, the actual implementation of circuit simulation still requires careful examination.First, we have applied the model to the generation of successive BD events, showing that it reproduces the experimental random event-time trajectories and the statistical distribution of successive BDs.This preliminary study sets up the validity of the approach.However, our main focus is RS devices.After the formulation of our model for these devices, we have reproduced the experimental results for both ramped voltage signals and trains of pulses with the same set of parameters.Being a stochastic approach, it has been shown that the model can deal with C2C and D2D variability, a key element for the realistic simulation of memristive circuits.

Fig. 1 .
Fig. 1.(a) Successive BD distributions: "on-the-fly" method (dots) and analytical result for the Weibull distribution.(b) Multiple BD events for different samples subjected to a constant voltage stress: experiment (yellow) versus "on-the-fly" method results (black).(c) Successive BD distributions (TDCM): experiment (dots) versus "on-the-fly" method (lines).(d) Evolution of the number of events under the application of a ramped voltage for different voltage acceleration factors by using the on-the-fly method.
with ζ a shape parameter and I 0 being related to as I 0 ( ) = (I max − I min ) + I min , where I max and I min correspond to the maximum ( = 1) and minimum ( = 0) currents.Notice that couples the two model equations.

Fig. 2 .
Fig. 2. Comparison of model (black lines) and experimental (red lines) set/reset current-voltage loops in (a) linear and (b) logarithmic current scales.

Fig. 3 .
Fig. 3. Effects of compliance current in the current-voltage loops.Experimental (a) and model (b) results.

Fig. 4 .
Fig. 4. Variability of I(V) loops.Ten different cycles are considered in the three cases.Red curves correspond to experiments and gray lines represent model calculations.(a) Experimental results.(b) Model results obtained from the generation rates of (8).(c) Model results considering an MPM based on the gamma distribution (see Section II) with α = 2.The latter shows a much larger dispersion of curves since it includes device-to-device variability.