Time series statistical analysis: A powerful tool to evaluate the variability of resistive switching memories

1 Time series statistical analysis: a powerful tool to evaluate the variability of resistive switching memories J.B. Roldán1*, F. J. Alonso2, A. M. Aguilera2, D. Maldonado1, M. Lanza3 1Departamento de Electrónica y Tecnología de Computadores. Universidad de Granada. Facultad de Ciencias. Avd. Fuentenueva s/n, 18071 GRANADA, Spain. *Corresponding author Email: [email protected] 2Departamento de Estadística e Investigación Operativa. Universidad de Granada. Facultad de Ciencias. Avd. Fuentenueva s/n, 18071 GRANADA, Spain. 3Institute of Functional Nano and Soft Materials (FUNSOM), Collaborative Innovation Center of Suzhou Nanoscience & Technology, Soochow University, 199 Ren-Ai Road, Suzhou, 215123, China. ABSTRACT Time series statistical analyses (TSSA) have been employed to evaluate the variability of resistive switching memories, and to model the set and reset voltages for modeling purposes. The conventional procedures behind time series theory have been used to obtain autocorrelation and partial autocorrelation functions and determine the simplest analytical models to forecast the set and reset voltages in long series of resistive switching processes. To do so, and for the sake of generality in our study, a wide range of devices have been fabricated and measured. Different oxides and electrodes have been employed, including bilayer dielectrics in devices such as: 2 Ni/HfO2/Si-n+, Cu/HfO2/Si-n+ and Au/Ti/TiO2/SiOx/Si-n+. The TSSA models obtained allowed to forecast the reset and set voltages in a series if previous values were known. The study of autocorrelation data between different cycles in the series allows estimating the inertia between cycles in long resistive switching series. Overall, TSSA seems to be a very promising method to evaluate the intrinsic variability of resistive switching memories. Index Terms—Resistive switching memory, RRAM, Conductive filaments, Variability, Time series modelling, Autocovariance, Stationary time series. 1 - INTRODUCTION Resistive random access memories (RRAM) have shown an outstanding potential for information storage, especially for internet of things (IoT) and related applications, due to their excellent scalability, low power operation, fast switching speed, easy fabrication and good compatibility with the well-established complementary metal-oxide-semiconductor (CMOS) technology [1, 2, 3, 4, 5, 6, 7]. The most promising RRAM unit cells consist of matrixes of two-terminal metal/insulator/metal (MIM) nanocells, in which the electrical resistance of the insulating film can be switched between a high resistive state (HRS) and a low resistive state (LRS) depending on the electrical impulses applied between the two metallic electrodes. However, after more than one decade of intense research, RRAM devices still have not been mass produced by the industry, nor entered in the market of non-volatile memories (NVM), mainly due to variability problems [1, 2, 5, 8, 9, 10]. In RRAM devices the resistive switching (RS) 3 is related to the formation and rupture of defect-rich conductive filaments (CFs) within the dielectric layer, which leads to a HRS-to-LRS (set) and a LRS-to-HRS (reset) transition (respectively). Set and reset transitions are related to diffusion, redox and nucleation of different chemical species within the MIM nanocells [1, 2, 5, 7, 10], which take place with a very high degree of randomness from one cycle to another, leading to an intrinsic high variability. Consequently, the electrical characteristics measured in a RRAM device reflects the stochasticity of these physical processes, and produces the so-called cycle-to-cycle variability. In the past few years, the variability of RRAM devices has been statistically analysed in most cases using the Weibull distribution (WD) [1, 11, 12]. The WD comes out in the field of reliability physics [13] and its use makes sense for RRAMs under filamentary conduction since it is a weakest-link type distribution, i.e. the failure of the whole is dominated by the degradation rate for the weakest element. However, this method does not describe all the inherent statistical particularities of RRAM devices: and although it has been previously employed to deepen on the experimental data characteristic obtained for different technologies [1, 11, 12], it does not entirely capture the essence of the RS process. In fact, classical reliability analyses with the Weibull distribution assume that times to failure, in our case set voltage (VSET) and/or reset voltage (VRESET), are independent within a RS series. This assumption may not be valid in the case of stochastic processes associated with RRAMs because successive observations could be highly dependent (in fact, a CF is formed making use of broken parts of previous ones). It is worth highlighting that other mathematical approaches have been proposed to tackle different facets of the statistical study of variability in RRAMs and, in a more general scope, thin dielectrics. Among them, the use of a clustering statistical approach complementing the use of the WD [14] can be counted; in 4 line with this, convolution-based modelling is also interesting being noticed [15]. Strategies making use of completely different distribution functions have also been reported; for example, the employment of phase-type distribution functions for certain devices led to interesting results [16]. Markov models have also been employed in the analysis of these devices [17, 18]. Kinetic Monte Carlo (KMC) and related simulations can also be considered as statistical tools to analyse RRAM variability, as they allow modifying very specific physical parameters of the devices (i.e. concentration of atomic vacancies/dopants, insulator thickness fluctuations) within a reasonable range and analyse the deviation of the electrical characteristics. A broad number of contributions have been presented in this respect [10, 19-25]. The only drawback of KMC simulations is the longer computational time, which obligates researchers to make assumptions that simplify the calculations, leading in some cases to a loss of accuracy. It is clear that variability is still an unresolved problem in RRAM devices both from the technological and the modelling viewpoint, and developing new analytical methods to shed light into this problem is highly necessary. Times series statistical analyses (TSSA) are powerful numerical methods that have been successfully applied for decades in the fields of economics and sociology, and more recently they have been also sporadically used in the field of engineering and reliability of electronic devices [26-27]. TSSA may be useful to analyse the variability of RRAM devices because: i) the data (VSET and VRESET) are collected in a continuous manner over the time (cycle-tocycle) for a long RS series [28-29]; ii) TSSA is appropriate for physical processes that exhibit any kind of inertia in some of their particular features [30, 31]. In this respect, in RS cycling, the CF is formed (set process) making use of the remnants of the CFs ruptured in the previous cycle (reset process). Therefore, from a statistical point of view it is relevant to analyse any numerical relations between neighbouring cycles and assess 5 the system “memory” in a long RS series (to put it in formal words, the parameters that characterize consecutive cycles in a RS series can be correlated and, therefore, the term autocorrelation comes up naturally). And iii) under certain mathematical conditions (i.e. time series stationarity, an assumption that our data distributions fulfil) a comprehensive analysis can be performed through a time series analysis approach [30, 31]. However, despite this strong parallelism, to the best of our knowledge, TSSA have never been employed to evaluate the variability of RRAM devices. In this work we present the first variability study of three different types of RRAM devices using the TSSA, and observe that in all cases essential RS parameters, such as VSET and VRESET, can be reasonably forecasted making use of mathematical models and the information of these parameters in previous cycles. In addition, in long RS series the dependence of forthcoming cycles on previous ones can be correctly studied by means of correlation and autocorrelation analyses, characterizing in this manner the inertia of RS operation in RRAMs for different technologies. The manuscript is organized as follows: in Section II the new model is described in depth, in Section III the fabricated devices and measurement process details are given, in Section IV the new statistical analysis is explained, and in Section V the main results and discussion are presented. Finally, the main conclusions are drawn in Section VI. 6 2 – MODEL DESCRIPTION Back in 1927, G.U. Yule introduced modern TSSA formulating a model for a pendulum dynamic movement time dependency [32]. After rearrangements in the corresponding equation describing the pendulum movement he came out with a second order autoregressive time series model where the pendulum displacement (zt) from the equilibrium position was regressed on the two previous observations (zt-1) and (zt-2) — the physics governing the pendulum is linked to a second order differential equation—. In our study, we model the values of VSET and VRESET of different RRAM devices over long series of RS cycles by considering the values of previous cycles using TSSA [30, 31, 33]. One of the main difficulties is to find the order of the model, i.e. how many VSET or VRESET values from previous cycles we need to forecast the current cycle (for a general model previous cycles are usually considered, see Equation 1). In addition, it is also necessary to find the weights (Φ1…Φp) of the autoregressive model we are seeking (see Equation 1). 𝑉𝑟𝑒𝑠𝑒𝑡𝑡= Φ1𝑉𝑟𝑒𝑠𝑒𝑡𝑡−1 + Φ2𝑉𝑟𝑒𝑠𝑒𝑡𝑡−2 + ⋯ + Φ𝑝𝑉𝑟𝑒𝑠𝑒𝑡𝑡−𝑝 + 𝜀𝑡 (Eq. 1) The order of the model (p) depends on the physics governing RS, but we will not assume any knowledge of it and we will only make use of the information within the experimental data because the underlying physics and the technological details of the fabrication process are “hidden” in the RS data collected. Therefore, the models obtained with this approach are empirical and the TSSA output will consist of the order of the model (p) and the weights set (Φ1, ..., Φp). Sometimes a model such as the one 7 described in Equation 1 works better for the centred variables; e.g., 𝑉𝑅𝐸𝑆𝐸𝑇 − 𝜇, where µ stands for VRESET mean in the RS series (this formulation is equivalent to include a constant term Φ0 in the model). The term εt, as usually employed in time series studies, stands for a residual that accounts for the model error (the difference between the measured value and the modelled value). This term is obtained in TSSA theory by generating random numbers with a normal distribution whose variance corresponds to the one calculated from the measured data. Nevertheless, for the sake of simplicity, this term will not be included in the models developed here, following the conventional notation in the engineering context. In this approach, the algebraic equations are employed in the form usually seen in compact modelling, i.e., the current value of reset or set voltages are given as a function of variables already known, such as reset and set voltages of previous cycles. The existence of a difference between the measured and modelled values is assumed. Equation 1 shows what is called an autoregressive (AR) model [30]; however, not always such an easy model can be obtained. Occasionally more complex models are needed; if this is the case, an autoregressive moving average (ARMA) model that includes AR and moving average (MA) parts are considered [30]. MA models are a linear combination of past residuals [30-32]. The general expression of an ARMA model is described in the Supplementary Information (see Equation A1). We have employed TSSA to study and model the data obtained from RRAMs made of different materials. Three types of devices were considered, two of them including HfO2 as the dielectric, sandwiched by different electrode materials (Ni and Cu), and another one based on a bilayer (TiO2/SiOX) insulating stack. 8 3 – DEVICE FABRICATION AND MEASUREMENT The HfO2-based RRAM devices consisted of Ni/HfO2/Si-n+ and Cu/HfO2/Si-n+ (20nm thick dielectric layers were used) stacks [28]. The active area of the MIM cells was 5 µm×5 µm, and the devices were patterned by photolithography. A HP-4155B semiconductor parameter analyser was used in the measurement process, which consisted of collecting long sequences (series) of current vs. voltage (I-V) curves by applying ramped voltage stress (RVS). The Si-n+ substrate (bottom) electrode was grounded and a negative voltage was applied to the Ni (0.3V/s ramps) or Cu (0.5V/s ramps) (top) electrode, although for simplicity we have assumed the absolute value of the applied voltage henceforth [28]. The RS phenomenon observed for both type of devices was unipolar, as displayed by the I-V curves (see Figures 1a, 1b) [5]. A RS series of 2800 cycles was obtained for the Ni/HfO2/Si-n+ RRAMs, and a series of 280 cycles for the Cu/HfO2/Si-n+ RRAMs. The values of VSET and VRESET were extracted from the I-V curves as reported in [1, 28, 34], and plotted in Figures 1b-1f. The variability of VSET and VRESET from one cycle to another can be clearly observed. 9 Figure 1. Typical I-V curves observed in (a) Ni/HfO2/Si-n+ and (d) Cu/HfO2/Si-n+ RRAMs. Experimental values of VSET (b, e) and VRESET (c, f) versus cycle number for a series of continuous RS cycles under RVS for RRAMs based on Ni/HfO2/Si-n+ and Cu/HfO2/Si-n+ stacks. The third type of RRAM devices was fabricated using an Au/Ti/TiO2/SiOX/Si-n+ structure. A 2nm TiO2 film was grown by atomic layer deposition (ALD), on Si-n+ wafer with a ~1.5 nm thick native SiOX layer. The device area was 5 µm x 5 µm and the 16 See that the sample PACF for Cu/HfO2/Si-n+ RRAMs has only two values (the first and second lagged cycles) above the threshold bound (Equation A9), so that an AR (2) model could be considered. However, taking into consideration that the second component is close to the threshold bound an AR (1) can be reasonable. After the AR (1) is adjusted all the validation considerations are satisfied. The VRESET model for these devices is given in Equation 5: 𝑉𝑅𝐸𝑆𝐸𝑇𝑡(𝑉) = 0.8577 + 0.5711 𝑉𝑅𝐸𝑆𝐸𝑇𝑡−1. (Eq. 5 ) For the Ni/HfO2/Si-n+ devices, the VRESET model works similarly to what was determined for the set voltage modeling (ACF and PACF decrease). Then an ARMA (1,1) model (Equation 5) holds for the VRESET time series for Ni/HfO2/Si-n+ devices: 𝑉𝑅𝐸𝑆𝐸𝑇𝑡(𝑉) = 0.2377 + 0.8573𝑉𝑅𝐸𝑆𝐸𝑇𝑡−1 + 0.6523𝜀𝑡−1. (Eq. 6) The residuals of the VRESET models depicted in Equations 5 and 6 have white noise behavior again, and confirm the appropriateness of the modeling procedure. For the Au/Ti/TiO2/SiOx/Si-n+ devices the ACF and PACF are given in Figure 6. In this case, the influence of other lagged cycles is more important than in previous cases; therefore, the reset voltage can be given by a linear combination of reset voltages obtained in the previous cycles. The model is a factorized AR (6) type. 17 Figure 6. (a) ACF and (b) PACF versus cycle lag for the reset voltage of the Au/Ti/TiO2/SiOx/Si-n+ devices under study and the RS series described in section II. The ACF and PACF minimum threshold bounds are 0.2 respectively, shown in blue lines. The analytical expression for the TSSA description of the set voltage in Au/Ti/TiO2/SiOx/Si-n+ devices is given by Equation 7. 𝑉𝑅𝐸𝑆𝐸𝑇𝑡(𝑉)= 0.3228 + 0.3198 𝑉𝑅𝐸𝑆𝐸𝑇𝑡−1(𝑉)+ 0.2197 𝑉𝑅𝐸𝑆𝐸𝑇 𝑡−5(𝑉) − 0.0703 𝑉𝑅𝐸𝑆𝐸𝑇𝑡−6(𝑉) (Eq. 7) The residuals of the VRESET models depicted in Equation 7 have white noise behavior again and consequently the modeling procedure is correct. Note that the coefficient of the (t-6) component is the product of the components (t-5) and (t-1) because of the characteristics of the parameter calculation procedure. The three technologies under study here show stationarity in the set and reset voltages series. Nevertheless, for other technologies a drift in the mean and variance shows up. In these cases, stationarity does not hold, therefore stationarity is not a general rule. If we are faced with a nonstationary data series, the methodology described in the Supplementary Information would not be appropriate and no models can be extracted. In this situations there can be other options since the TSSA theory proposes 18 changes of variables that lead the newly derived series to fulfill the stationary requirements that are needed prior to the modeling process. Autoregressive integrated moving average (ARIMA) approaches can be employed instead of the AR or ARMA modeling schemes explained above, Ref [30, 31]. 5 - RESULTS AND DISCUSSION In order to test the accuracy of the models previously developed we superposed in the same graphic the measured VSET and VRESET with the modeled ones for the devices under consideration in this manuscript. The modeling is a forecast of the actual value considering previous measured values, as it is conventionally done in TSSA. These results for VSET (Figure 7) and VRESET (Figure 8) are plotted taking into consideration Ni/HfO2/Si-n+ and Cu/HfO2/Si-n+ devices. As can be seen, the VSET mean general trend is described reasonably well by Equations 2 and 3. The main dependencies have been correctly analyzed and incorporated with our procedure; hence, within the time series context, the model works well. We have validated this point by studying the residuals correlation, and we did not obtain any significant correlation between the residuals of the current cycle and those of the lagged ones (this constitutes the validation step, as explained in the Supplementary Information). Consequently, no more dependencies have to be incorporated to the models, since no information is statistically “hidden”. Although the scales and the modeling strategies are different for VSET, i.e. AR (1) for Cu/HfO2/Si-n+ devices and ARMA (1,1) for Ni/HfO2/Si-n+ devices, the accuracy is similarly reasonable, as Figure 7 shows. The model error (εt), as highlighted before, presents a white noise structure. We have also checked that the model accuracy is maintained if the number interval is changed. 19 Figure 7. VSET versus cycle number for the RS series under consideration. The measured values are shown in blue lines and the modelled ones in red. (a) Ni/HfO2/Sin+ RRAMs, cycles 1000-1200, (b) Ni/HfO2/Si-n+ RRAMs, cycles 1800-2000 (c) Cu/HfO2/Si-n+ RRAMs, cycles 50-150, (d) Cu/HfO2/Si-n+ RRAMs, cycles 50-250. The results for VRESET are in line with those of VSET. In this case, an AR (1) model is used for the Cu/HfO2/Si-n+ devices and ARMA(1,1) for Ni/HfO2/Si-n+ devices. Again, the model reproduces accurately the VRESET mean evolution for all the cycle number intervals considered, as displayed in Figure 8. 1000 1050 1100 1150 1200 2 3 4 1800 1850 1900 1950 2000 2 3 4 50 100 150 1 2 3 4 5 6 7 8 50 100 150 200 250 1 2 3 4 5 6 7 8 (a) Ni VSET (V) Cycle number Measured Modeled (b) Ni VSET (V) Cycle number Measured Modeled Measured Modeled (c) Cu VSET (V) Cycle number VSET (V) Cycle number Measured Modeled Cu (d) 20 Figure 8. VRESET versus cycle number for the RS series under consideration. The measured values are shown in blue lines and the modeled ones in red. (a) Ni/HfO2/Si-n+ RRAMs, cycles 1000-1200, (b) Ni/HfO2/Si-n+ RRAMs, cycles 1500-1700 (c) Cu/HfO2/Si-n+ RRAMs, cycles 50-150, (d) Cu/HfO2/Si-n+ RRAMs, cycles 50-250. For the Au/Ti/TiO2/SiOx/Si-n+ devices the comparison between measured and modelled data is shown in Figure 9. The prediction is also reasonable although more values of set (Equation 4) and reset (Equation 7) voltages of previous cycles are considered. 1000 1050 1100 1150 1200 1 2 3 1500 1550 1600 1650 1700 1 2 3 50 100 150 0 1 2 3 4 50 100 150 200 250 0 1 2 3 4 Ni VRESET (V) Cycle number Measured Modeled (a) Ni VRESET (V) Cycle number Measured Modeled (b) Cu VRESET (V) Cycle number Measured Modeled (c) Cu VRESET (V) Cycle number Measured Modeled (d) 21 Figure 9. (a) Set (b) Reset voltage versus cycle number for the RS series under consideration. The measured values are shown in blue lines and the modeled ones in red. Au/Ti/TiO2/SiOx/Si-n+ devices are considered here. It can be observed in the figures above that in certain cases the fit is smaller than current values. This effect is seen when values much different than the mean show up. However, we would like to highlight that the modeling methodology we are presenting deals well with the prediction of the set and reset voltage mean. The current voltage values are also predicted well in most cases by using the statistical information of previous cycles. Therefore, taking into account that current models in the literature do not have this information, and that prediction of variability is interesting in devices (RRAMs) that show inherent stochasticity, we believe that, although it is not a modeling final solution, this technique is a step forward that can be worthwhile to characterize the device physics and help with variability modeling. It is important to highlight that Figure 3 and Figure 5 (in addition to the data needed for model building) provide information about the RS processes of the HfO2based devices analysed. In Figure 3a, we can see that the correlation between cycles is higher for the Cu/HfO2/Si-n+ devices for the first lag cycles with respect to the set voltage. In Figure 5a, a similar trend can be observed for the reset voltage series. So, in general, the influence of previous cycles in Cu/HfO2/Si-n+ devices is higher; i.e., the CF remnants from resets processes influence more the following set cycles, and that is why 020 40 60 80 100 0.5 1.0 1.5 020 40 60 80 100 -1.0 -0.8 -0.6 -0.4 -0.2 (a) VSET (V) Cycle number Measured Modeled (b) VRESET (V) Cycle number Measured Modeled 22 a higher correlation comes out for the set voltage values when several lag cycles are considered. The qualitative explanation for this behaviour could be in the nature of RS in these devices. It is known that devices with Cu and Ag electrodes are employed in conductive-bridge RAMs [2, 36, 37] because of the capacity of Cu and Ag cations to diffuse in the dielectric and form, after a reduction process, a metallic-like conductive filament. In these devices, the reduction and oxidation potential are described by a thermally activated process whose activation energy depends on the number of atoms surrounding the one which is taken into consideration [1, 23, 24]. In this manner, the formation of the percolation path introduces “inertia” as the reduced atom clusters grow denser in the dielectric since they tend to maintain their shape hindering the oxidation processes of their atoms. This behaviour could explain the higher RS “inertia” shown by Cu/HfO2/Si-n+ devices since they maintain better the CF form and size that determine VSET and VRESET values. These effects would be reflected in the ACF plot with higher autocorrelation values for cycles not distant away in the series. The lower autocorrelation (Figures 3a and 5a) for the Ni/HfO2/Si-n+ devices reflects a lower RS “inertia”, as highlighted above. That could also be linked to a mixture of RS phenomena, since oxygen vacancies could be also involved in their resistive switching operation, as suggested in Ref. [38]. In this respect, the effects linked to activation energy lowering for oxidation of clustered metal atoms from the electrode would be mitigated, producing less correlation between the set and reset values of consecutive cycles. In both cases, Ni/HfO2/Si-n+ and Cu/HfO2/Si-n+ devices, the most significant dependency (correlation) is with the previous value (first lag). In line with these latter issues, the results of Au/Ti/TiO2/SiOx/Si-n+ devices can be analysed. In Figures 4 and 6, 23 the cycles correlated in the ACF and PACF (mostly for the reset voltage) is higher than for the technologies studied previously. In particular, for the reset voltage, the correlation of cycles 9, 11, 13 is high with respect to the previous cases. This fact shows again an important inertia in the RS features of this technology in comparison to the devices based on HfO2. The nature of CFs in TiO2 for the filamentary current component and the presence of a volume current component (this component is linked to the ion distribution that is spread out in the dielectric and that can affect several consecutive cycles in a RS series, as it was shown in devices of this kind [29]), is key to explain the higher correlation highlighted above. It is important to highlight that the autocorrelation and partial autocorrelation functions employed to in the analysis presented above (Figures 3-6) can be used separately from the TSSA modeling. They reflect the correlation of set and reset voltages between the different cycles, this means extracting the dependencies of the current cycle on the previous ones. This information is useful even if no TSSA modeling is performed since it shows that the values studied as independent data, are, in fact, dependent, and this fact has implications in the study of cycle-to-cycle variability. It is important to highlight that ACF and PACF analysis could be easily performed on RRAM measured data to assess the correlation between the characteristics of successive RS cycles. From the modeling viewpoint, TSSA models could be implemented in circuit simulators with Verilog-A compilers. TSSA models could be embedded in previous models to account for the RRAM stochastic behavior and for the correlation of certain parameters such as set and reset voltages in long RS series, as explained above. We have done so making use of the RRAM Standford model [40]. In our case, a log file of previous set and reset values for the RRAMs has been employed in a model built upon the TSSA 24 approach making use of regressed values. The file was written with the current values of set and reset voltages every time a set or reset was performed. Previously, in order to simulate the RVS case we described experimentally above, we implemented a modification of the local enhancement factor (γ) [41] to be able to obtain the corresponding set and reset voltages, since these latter parameters are not model parameters. Finally, we would like to comment on the fact that the operation regime of these devices would be characterized by pulses if they are used as storage-class memory in different chips or spikes of a variety of shapes if employed to mimic synapses in neuromorphic circuits. The device conductance depends on the pulse number for each signal amplitude, because of the different thermal inertia that is produced in each operation regime [39]. In this respect, we have chosen RVS measurements to show a particular application of the TSSA in a well-known characterization approach. This approach could be considered as the DC facet of a classical compact model, further developments could be needed to deal with transient events in a more general model where thermal and capacitive effects would be needed. VI.-CONCLUSIONS Times series statistical analyses (TSSA) have been used to study long series of resistive switching processes. The experimental data analyzed here were measured in resistive random access memories with Ni/HfO2/Si-n+, Cu/HfO2/Si-n+, Au/Ti/TiO2/SiOx/Si-n+ structures. The conventional time series techniques were applied 25 to model the VSET and VRESET of these devices; to do so, autocorrelation functions and partial autocorrelation functions were obtained for all the types of RRAMs. Autoregressive models were obtained for Cu/HfO2/Si-n+ devices and the autocorrelation function between cycles was high, showing an important inertia between resistive switching cycles. The better diffusion of Cu ions in the dielectric is behind this behaviour. For Ni/HfO2/Si-n+ devices more complex models are needed and autocorrelation data show less inertia between resistive switching cycles. Finally, for Au/Ti/TiO2/SiOx/Si-n+ devices a significant correlation can be observed for more distant cycles in the reset voltage description. A physical explanation has been developed in connection with the correlation results found. The models obtained can be used to forecast the values of set and reset voltages in a resistive switching series if previous values are known. The information obtained in this context can be employed in modelling and in the characterization of RRAM variability. 32 [28] M. B. González, J. M. Rafí, O. Beldarrain, M. Zabala, and F. Campabadal, “Analysis of the switching variability in Ni/HfO2-based RRAM devices,” IEEE Trans. Device Mater. Reliab., vol. 14, no. 2, pp. 769–771, 2014. [29] N. Xiao, M.A. Villena, B. Yuan, S. Chen, B. Wang, M. Elias, Y. Shi, F. Hui, X. Jing, A. Sheuermann, K. Tang, P.C. McIntyre, M. Lanza, "Resistive Random Access Memory Cells with a Bilayer TiO2/SiOXInsulating Stack for Simultaneous Filamentary and Distributed Resistive Switching", Advanced Functional Materials, 27, p. 1700384, 2017. [30] P. J. Brockwell, R. A. Davis, “Introduction to Time Series and Forecasting”, Second Edition, Springer, 2002. [31] S.Bisgaard, M.Kulahci, “Time series analysis and forecasting by example”, Wiley, 2011. [32] G.U. Yule, “On a method of investigating periodicities in disturbed series, withreference to Wolfer’s Sunspot Numbers”, Philosophical Transactions of the Royal Society of London, Series A, 226, pp. 267–298, 1927. [33] D. C. Montgomery, C. L. Jennings, M. Kulahci, “Introduction to Time Series Analysis and Forecasting”, Wiley Series in Probability and Statistics, Wiley, 2015. [34] M.A. Villena, J.B. Roldán, F. Jiménez-Molinos, J. Suñé, S. Long, E. Miranda, M. Liu, "A comprehensive analysis on progressive reset transitions in RRAMs", Journal of Physics D: applied physics, 7, 205102, 2014. [35] A. Coghlan, “A little book of R for time series”, Release 0.2. In: https://media.readthedocs.org/pdf/a-little-book-of-r-for-time-series/latest/a-little-bookof-r-for-time-series.pdf, 2018. 33 [36] Yang, Y. et al. Electrochemical dynamics of nanoscale metallic inclusions in dielectrics. Nat. Commun. 5, 4232, 2014. [37] Yuanyuan Shi, Xianhu Liang, Bin Yuan, Victoria Chen, Haitong Li, Fei Hui, Zhouchangwan Yu, Fang Yuan, Eric Pop, H.-S. Philip Wong, Mario Lanza ,”Electronic synapses made of layered two-dimensional materials“, Nature Electronics 1, 458–465, 2018. [38] U. Celano, L.Goux, R. Degraeve, A.Fantini, O. Richard, H. Bender, M.Jurczak, W.Vandervorst, “Three-Dimensional Observation of the Conductive Filament in Nanoscaled Resistive Memory Devices”, Nanoletters, 14, pp. 2401-2406, 2014. [39] S. Kim, C. Du, P. Sheridan, W. Ma, S. Choi, W. D. Lu, “Experimental Demonstration of a Second-Order Memristor and Its Ability to Biorealistically Implement Synaptic Plasticity”, NanoLetters, 15, 2203, 2015. [40] https://nano.stanford.edu/stanford-rram-model [41] X. Guan, S. Yu, and H-S.P Wong, “A SPICE Compact Model of Metal Oxide Resistive Switching Memory with Variations”, IEEE Elec. Dev. Lett., vol. 33, pp. 1405-1407, 2012.