Compact modeling of the switching dynamics and temperature dependencies in TiOx-based memristors: Part I-Behavioural Model

Memristor is a promising device as a fundamental building block for future unconventional system architectures such as neuromorphic computing, reconfigurable logic and multibit memories. Therefore, to facilitate circuit design using memristors, accurate and efficient models spanning a wide range of programming voltages and temperatures are required. In the first part of this series, we propose a behavioural model for temperature dependence of non-volatile switching dynamics of TiOx memristors. We begin by describing pulsed resistance transients of the memristors and then we use a multi-stage methodology to establish bias and temperature dependence of the model parameters. The proposed model is then shown to accurately describe the pulsed resistance transient characteristics of Pt/TiOx/Au and Pt/TiOx/Pt memristors.


I. INTRODUCTION
The demonstration of memristors [1] in the past decade has opened up new research frontiers. These devices have been shown to be promising for applications in memory [2], neuromorphic computing [3], [4], [5], [6], [7], [8], reconfigurable logic [9], [10], [11], sensing etc. To realise several potential applications, CMOS integration of memristor technology must be achieved. And for circuit designers to make use of this hybrid memristor-CMOS technology, it is important to capture the electrical characteristics of memristors with reasonable accuracy and computational efficiency. In the past decade, several models have been proposed to describe the currentvoltage (I-V) characteristics and mechanisms of change of the internal state variable(s) of memristors. These mechanisms are often related by modelers to electrochemical principles. For example, in the ion drift model the internal state variable is interpreted as the position of the boundary between oxygenrich and oxygen-deficient regions in the active layer oxide, whereas in TiO x memristors, it is the barrier height or width at the electrode oxide interface [1], [12], [13], [14], [15], [16], [17]. In general, internal state variable dictates the low field I-V relationships of memristor devices, whilst the rate of change of internal state variable is dictated by the applied programming or switching bias. Previously, it has been shown that in general for a variety of voltage controlled memristor The authors are with the Centre for Electronics Frontiers, University of Southampton, Southampton, SO171BJ, UK (e-mail: t.prodromakis@soton.ac.uk) devices, a threshold exists for the programming voltage. The threshold voltage in this context is roughly defined as a voltage below which, it is not possible to obtain a significant change in resistance in a practical time-frame. Therefore, it is possible to separate the two regimes viz. low field regime or non-switching regime and programming or switching regime for the modeling purposes. Such an approach was used in the VTEAM model [15]. More recently it was shown that it is possible to describe resistive switching in memristors independently of the low field device operation i.e. a physics agnostic approach is used to define the state variable [14]. Because the model is behavioural and stresses speed over any other concerns, the memristor is modelled as having a single state variable: its resistive state. Next, in the same work, a protocol involving pulsed resistance transient (PRT) measurement of the state variable (device resistance) at a fixed read-out voltage is proposed for extracting the parameters of the switching dynamics model.
In this work we begin from the memristive model described in [14] and propose a method to extract improved switching dynamics components (leaving the static, i.e. non-switching component the same) that also take into account the temperature dependence. The proposed components do not affect the static I-V relationship of the device. Moreover, model parameter extraction is carried out using the same form of PRT measurements proposed in [14], thus maintaining testing routine compatibility. The proposed model describes resistive switching in terms of fractional resistance change d∆R/dt as opposed to simply dR/dt from our previous work. This allows us to reduce the size of the parameter space substantially.
This paper is the first of a 2-part series: in part I we demonstrate a methodology for extracting behavioural models for the switching dynamics of TiO x memristors including temperature dependence. We also show a simple technique used to reduce the model parameter size. In part II we propose a physics-based model for the temperature dependence, where physical realism has been stressed over performance.

II. EXPERIMENTAL METHODS
In this work we have used two different types of TiO x based memristor devices. Schematic representation of both the type of memristor devices is shown in Fig. 1. Memristor I-a and I-b (jointly referred to as memristor I-* in the remainder of the paper) are Pt/TiO x /Pt with Pt as top and bottom electrode and memristor II is Pt/TiO x /Au with Pt as top electrode and Au as bottom electrode. Both the type of devices are fabricated on two separate 6-inch oxidized (200 nm dry oxide i.e. SiO 2 ) silicon wafers. Bottom electrodes are defined using photolithography followed by e-beam evaporation of adhesive 5 nm Ti layer and then 12 nm Pt (20 nm Au) for memristor I-* (memristor II). This deposition is followed by lift-off process in N-Methyl-2-pyrrolidone (NMP). Followed by subsequent photolithography, TiO x active layer of thickness 10 nm (25 nm) for memristor I-* (memristor II) is then deposited using reactive magnetron sputtering of Ti target. Lift-off process was then carried out using NMP followed by photolithography for the top electrode. The e-beam evaporation of 15 nm Pt (12 nm Pt) for memristor I-* (memristor II) and lift-off processes are then carried out.
For electrical characterization we use an Arc One from Arc Instruments [18], [19]. The electrical characterization of To achieve memristive operation in our TiO x devices, we use a pulse-based electroforming protocol [20], [21]. And then, to obtain the PRT characteristics, we apply train of pulses of fixed amplitude. The amplitude of the pulse train is swept across biases in the switching regime and at each bias step both polarities are visited. An example of such a bias sequence is 0.8, -0.8, 0.9, -0.9, 1.0, -1.0 V. In this paper, the application of train of pulses at a particular amplitude is referred to as the switching bias. Each switching bias consists of 500 pulses (for memristor I-*) and 200 pulses (for memristor II) of 100 µs pulse width. In between the consecutive programming pulses a read pulse of 200 mV is applied to measure the resistance of the device. Since the model presented in this work is intended for application in non-volatile memristors, for a short interpulse and read pulse duration, device retention is guaranteed and device resistance is unchanged. Fig. 2 shows the biasing protocol employed in the PRT measurements. Fig. 2(a) shows the typical response of a memristor device observed in PRT measurements. Fig. 2 shows the typical sequence of switching biases employed in the PRT measurements and the outset of the figure shows pulse and read scheme which constitute as an elementary unit of the pulse train. For memristor I-b and for memristor II, we have also performed PRT measurements at temperatures ranging from 313K -353 K and 300 K -360 K correspondingly. For memristor I-b, these temperature-dependent measurements are achieved using a chuck-based temperature controller whereas for memristor II a micro-chamber based temperature control system [22] is used.

III. SWITCHING DYNAMICS MODEL
In the switching regime voltage biases between top and bottom electrodes greatly exceed the read-out voltage: |V | > |V read |. If sufficiently high, this voltage may cause the resistance of the memristive device to change. In the case of pulse train stimulation, it is typical that the first switching pulse will . change the resistive state significantly with consecutive pulses adding increasingly small increments (or decrements) to the resistance value. This saturation effect can be seen in Fig. 2(a) and provides a natural bound on dR/dt and it is described by an exponential decay. The similar PRT behaviour is reported in [13], [14], however instead of applying the exponential on the absolute resistance we use the change in resistance vs. an initial reference point: where R 0 is the resistance value of the memristor before the application of switching bias. The rate of change of the "R 0referred resistance difference" is now modelled by: where, s and R p are the fitting parameters. V is the magnitude of switching bias and T is the temperature. If we assume that prior to applying the switching bias, memristor has reached equilibrium or steady state then, d∆R/dt| t=0 − = 0. Then, the solution to eq. (2) is: where, ∆R 0 is the initial discrepancy from R 0 . For the switching bias which consists of pulse trains, eq. (3) can be expanded to express the resistance at the end of the n th pulse as, where, t w is the pulse width. For the very first pulse in the switching bias initial relative resistance is zero, i.e. ∆R 0 = 0. Note that this "R 0 -referred resistance empirical model" uses two parameters where as the previous models used three parameters [14], [13]. This impacts the overall parameter space significantly as we shall see in the next section, which discusses the bias and temperature dependence of the model parameters and how to establish these relationships using a multi-stage methodology.
IV. MODEL FITTING METHODOLOGY The parameters, s and R p presented in eq. (2) have a dependence on temperature and applied switching bias. To establish the nature of s(V, T ) and R p (V, T ), we use a threestage methodology as illustrated in Fig. 3. With every stage the size of the parameter space expands. In this methodology, except for stage I, the mathematical expressions used to expand the parameter space are not assumed a priori. Rather, an appropriate mathematical expression is chosen based on the extracted parameters from the previous stage. Initially, pulsed transients of ∆R measured at a range of temperatures are fitted at each switching bias independently. With this we obtain a scatter of s and R p for various temperatures and switching biases. In stage II the functional shape of bias-dependence is established for each temperature, i.e. we look at how s and R p vary with bias, establish the functional form of this dependence and extract parameters for each temperature. In stage III, we determine the temperature-dependence of the parameters extracted during stage II. This process is repeated for each voltage bias and polarity.
Depending on the mathematical functions used at each stage, the parameter space may expand to a very large one. The number of parameters expands in a multiplicative manner at each stage, therefore it is important for computational efficiency to reduce the number of parameters at all stages of the methodology. Previous work uses a three-parameter model in eq. (2), ending up with a parameter space of 3 × 2 × 3 = 18 parameters/polarity (ppp). As we discuss in the next section, we use a two-parameter model in eq. (3), reducing the total size of the parameter space to 2×2×3 = 12 ppp and increasing computational efficiency.

V. RESULTS AND DISCUSSION
In this section we present our switching dynamics models for memristor I-* and memristor II. We begin with demonstrating our new switching dynamics model at room temperature on memristor I-a. For this, we measure PRT as in Fig. 2. Switching biases of 1.4-2.0 V of alternate polarity are applied consecutively for 20 times on memristor I-a to check model robustness. The PRT and stage I fitting are shown in Fig. 4. Note that, in Fig. 4(a) the pulsed transients of ∆R are shown with the corresponding absolute resistance values shown in Fig. 4(b). Even though the values of absolute resistance do change in consecutive measurements, ∆R does not (to within approximation). Therefore, if we fit the model to the average values of the ∆R transients, we can obtain a robust set of model parameters.
The fittings to the average of pulsed transients of ∆R for memristor I-a at stage I and at stage II are shown in Fig. 4(a). The average of ∆R is obtained from twenty consecutive PRT measurements data shown in Fig. 4(b). Stage I fitting yields the parameter values (s and R p ) at the measured switching biases. Then, in order to establish the bias dependence of s and R p we use the exponential expression for s and second order polynomials for R p in stage II as follows: Fig. 4(a) shows the fitting of pulsed transients of relative resistance at the switching biases of ±1.4 V to ±2.0 V swept with a ±0.04 V step size. Fig. 5(a-d) shows the bias dependence of s and R p and the corresponding fits using eq. (5) and (6). Fig. 4(a) also shows the simulated pulsed transients of relative resistance using stage II fitting parameters i.e. s A , s k , p 0 , p 1 and p 2 . It can be seen that the simulated pulsed ∆R transients match very well to the average of the pulsed transients of ∆R and the stage I fit. Fig. 6 presents the comparison of previous work [14] and this work. Fig. 6(a) shows the model fittings at stage I in  both the works. We can notice that both the model fittings overlap indicating that there is negligible difference. The percentage fitting error is presented in Fig. 6(b). Obviously, both the models present no difference in the fitting error. The maximum error at stage I is nearly 5% in both the models. Given that the error difference between the two models is negligible, the two-parameter model presented in this work has an obvious advantage in terms of computational efficiency over the previous three parameter model.
We now illustrate the switching dynamics model for memristor II. Fig. 7 shows the PRT characteristics measured at different temperatures from 300 K to 360 K. It can be observed  from Fig. 7(a-g) that, the absolute resistance of the device decreases with the increase in temperature. We observe such a trend in all our TiO x devices. This thermal behaviour is a result of Schottky conduction mechanism and is discussed in detail in the part II of the series. Fig. 7(h) also shows the stage I model fit to the PRT characteristics at the corresponding temperatures. Fig. 8 shows s(V ) and R p (V ) plotted at various temperatures. For memristor II, the bias dependence of s (at different temeperatures) for both positive and negative switching biases is modeled as constant. The average values are used as the stage II model parameters. R p (V ) are modeled using exponentials as, The stage III fittings to s(T ), A(T ) and k(T ) are shown in Fig. 9. Second order polynomials are used to parameterise their temperature dependence. Second order polynomials are used for its generic applicability and simplicity, however this strictly data-driven approach in using polynomials would impose limitations on extrapolating outside of fitted range. For example, extrapolating s p,n (T ) and A p,n (T ) to very high temperatures might change the sign suggesting that model would no longer be applicable beyond a certain temperature. All the coefficients of these second order polynomials constitute a full temperature-dependent switching dynamics model of memristor II. Therefore at this final stage of the modeling methodology we must simulate the PRT entirely using the stage III parameters and compare it with the PRT data. The resistance transient characteristics determined by only using the stage III parameters are plotted (with solid black line) in Fig. 7 (a-g). It can be seen that the simulated characteristics are well matched with the measured experimental data. s(V, T ) is not expanded in stage II and modeled as scalar value equaling the average of all stage I values. In stage III the thermal dependence of s(V, T ) is modeled as a second order polynomial. Therefore modeling of s(V, T ) takes in total 1x1x3=3 ppp. On the other hand R p (V, T ) is modeled as a two-parameter exponential mentioned in eq. (6) in stage II and then thermal dependence of each stage II parameter (A(T ) and κ(T )) is modeled with a three-parameter second order polynomial. Therefore the modeling of R p (V, T ) takes 1x2x3=6 ppp. The final size of the parameter space at the end of the stage III for memristor II is 9 ppp. The model summary for memristor II is presented in Table II in the appendix. These parameters can be used in verilog-A for circuit design.
We also present fitted compact model for memristor I-b. Model fittings at stage-I and at stage-III are shown in Fig.  10. We observe that the model describes the PRT behaviour of memristor I-b well. At 313 K and at 353 K, the model predictions for positive polarity switching biases seem to produce significant errors at stage III whereas stage I fitting errors are within reasonable limits (Fig. 10(f)). This increase of the fitting error can be due to outliers in the s(V, T ) and R p (V, T ) scatter obtained in stage-I.

VI. CONCLUSIONS
We have presented a methodology for extracting improved switching dynamics models for memristive devices and demonstrate how it is used to extract models for two memristive devices. The methodology can be used to extract thermal characteristics of resistive switching, which previous works did not include. We observe that although the absolute values of the resistance change in the consecutive measurements, the "R 0 -referred" change in resistance does not, thus giving rise to a robust switching dynamics model. This technique also reduces the size of the parameter space at a fundamental level i.e. in stage I of the fitting methodology, making the model intrinsically more computationally efficient. The switching dynamics is modelled purely behaviourally in this work, using a multi-stage fitting methodology to extract the complete set of model parameters. Since the bias-and temperature-dependence of the model parameters is not assumed a-priori and instead allowed to emerge through a multi-stage methodology, the behavioural modelling methodology is generally applicable and can in principle be used for different types of TiO x memristors.