Holistic Variability Analysis in Resistive Switching Memories Using a Two-Dimensional Variability Coefficient

We present a new methodology to quantify the variability of resistive switching memories. Instead of statistically analyzing few data points extracted from current versus voltage (I–V) plots, such as switching voltages or state resistances, we take into account the whole I–V curve measured in each RS cycle. This means going from a one-dimensional data set to a two-dimensional data set, in which every point of each I–V curve measured is included in the variability calculation. We introduce a new coefficient (named two-dimensional variability coefficient, 2DVC) that reveals additional variability information to which traditional one-dimensional analytical methods (such as the coefficient of variation) are blind. This novel approach provides a holistic variability metric for a better understanding of the functioning of resistive switching memories.


Supplementary Note 1: Functional coefficient of variation calculation
We describe here the methodology for this new functional variability analysis in a formal manner to allow the reader to understand the steps followed in the calculations shown in the main manuscript.
The property of continuity in quadratic mean guarantees the continuity of the covariance function [Todorovic1992], which is crucial in many of the functional techniques, e.g., to obtain the spectral decomposition of the covariance operator required in Functional Principal Component Analysis (FPCA). In the context of resistive random access memories, FPCA based on Karhunen-Loève expansion was employed to describe the stochastic evolution of the main features related to reset curves in [AguileraMorillo2019, Aguilera2021, Ruiz-Castro2021].
The two main issues found in practice for the proposed variability analysis are the following: 1. The curves are not defined in the same domain because each curve has a different set/reset point or a different return point in the ramped voltage signal employed. This aspect is very important in the context of functional data analysis; therefore, "normalization" of the curve domain is needed. Since we are dealing with I-V curves measured under RVS, we have to use a normalized voltage.
2. The curves are observed in a discrete approach (non continuous), taking into consideration a finite grid of voltage values { 1 , 2 , … , } with being the voltage point where the i th curve is cut. Notice that the voltage values are ranked in increasing order. For the I-V curve domain normalization process (on the abscissa axis), the registration in the interval [0,1] is achieved by just dividing * = / with = 1, 2, … , (see, e.g., [Aguilera2021]). Hereinafter, the I-V curves defined in the interval [0,1] (see Figure S1) will be denoted as { ( * ): = 1, … , ; * ∈ [0,1]}. See that is the number of curves measured in the resistive switching series. For the Pt/TiO2/Ti devices equals 3900 and for the Au/Ti/TiO2/Au devices equals 749.
Once the curves are normalized, the second step of the analysis is performed. To do so, it is usual to consider a finite number of basis functions { 1 ( * ), … , ( * )} to obtain the real functional form of the sample I-V curves we are employing for the variability study [Ramsay2005, Acal2022]. In this manner, the normalized I-V curves can be expressed as follows, where are the basis coefficients (random variables with finite variance). Generalized cross validation can be used to determine the dimension of the basis of functions [Craven1978], whereas the choice of basis type depends on the I-V sample curves features. The most useful basis systems are Fourier functions for periodic data, B-Spline basis for non-periodic and smooth paths and wavelets basis for curves with a strong local behaviour. Given that B-Splines bases are considered in the current manuscript, a brief summary is given next: Let 0 < ⋯ < be a partition of knots (points) in the interval [0,1] in the abscissa of the I-V curve. More knots can be added to the partition as − < ⋯ < −2 < −1 < 0 < ⋯ < < +1 < +2 < ⋯ < + . Then, a B-spline basis of degree can be iteratively defined by A deep review about these and other bases, as well as their handling in the statistical software R can be checked in [Ramsay2005, Ramsay2009]. In particular, a cubic B-splines basis ( = 3) is a suitable option for the set/reset curves, given that these curves do not usually exhibit too much noise. Then, the B-spline basis coefficients for each curve are estimated by the least squares method as, The experimental I-V curves that we have employed in this study, for two different technologies, were reproduced by a cubic B-spline basis with 24 functions for each curve after normalization ( Figure S2 shows that this is an optimum choice for a correct accuracy). Figure S2. Normalized I-V curve and the fitting obtained with a cubic B-spline basis. As the number of functions in the basis increases, the accuracy of the fit improves significantly.
We performed the rebuilding process for all the curves in the RS series obtained in the laboratory. The results are shown in Figure S3 for the set and reset curves after the normalization process. Several classic tools used in exploratory analysis (the classical unidimensional statistical analysis commonly employed in conventional studies) are generalized to the functional framework (more details about descriptive statistics of functional data can be found in [Shang2015]). We have summarized some of them below. See that we change the domain from a one-dimensional dataset (e.g. reset voltages) to a two-dimensional dataset (e.g. the I-V curves for the reset processes).
The correlation plots corresponding to the r function described in Equation 7 are presented in Figure S4. They show in general low autocorrelation among the values of the curves at different normalized voltages, although the obtained data are nonuniform, see higher autocorrelation values at normalized voltages close to 1 for the set processes. Figure S4. a) Contour plots corresponding to the correlation surfaces for the set processes in Pt/TiO2/Ti devices, b) Contour plots corresponding to the correlation surfaces for the reset processes Pt/TiO2/Ti devices. c) Contour plots corresponding to the correlation surfaces for the set processes in the Au/Ti/TiO2/Au devices, d) Contour plots corresponding to the correlation surfaces for the reset processes in Au/Ti/TiO2/Au devices.
The functional coefficient of variation (FCV) can be considered to assess the variability of a set of curves, which is the case here. In this work, the classical coefficient of variation (CV = / ) is employed for the one-dimensional dataset and the FCV for the two-dimensional dataset. We are going to analyse here variability accounting for the whole set of I-V curves. This is a much accurate methodology than the one-dimensional statistical study of the set/reset voltages or any other parameter because of we are taking the complete curves (I-V curves) into account.
proposed a multivariate FCV. This new measure provides a single value representing the total variability in the whole domain (a single-point functional coefficient of variation that we name two-dimensional variability coefficient (2DVC) for convenience). By considering the basis expansion of the I-V curves defined in Equation (1), they adapt the multivariate coefficient of variation of Albert-Zhang type (see [Albert2010] for more details) as follows, where = ∫ ( * ) ( * ) 1 0 is the matrix (order × ) of inner products between basis functions, and ̅ are the covariance matrix (order × ) and the mean vector (order × 1) of basis coefficients, respectively. To determine the 2DVC value, the parameters and ̅ can be estimated by means of classical and robust estimators [Albert2010, Davies1987, Rousseeuw1985].
The main advantage of robust estimators in comparison with the classic ones is that the formers have a better performance when there are many outlying observations in the sample. The presence of outliers might have an important influence in the modeling of functional data. Multiple approaches have been developed with the goal of identifying functional outliers [Hyndman2007, Febrero2008]. Finally, the R-package rainbow contains functions for visualizing functional data and identifying outlier curves from different perspectives [Hyndman2010].

Calculation of the functional coefficient of variation for the bivariate case
In this section we consider the I-V curves measured for the forward ramp and for the reverse ramp. In this respect, we would have a function with two different current values for the same voltage. This cannot be considered in the simple functional data analysis. Therefore, we adapt the theoretical background with the bivariate case (the manner to consider the two I-V curve sections).
The previous methodology is only useful for the situation in which one functional variable is considered (the voltage (corresponding to the X axis) sweeps between 0V to Vset for the set process and 0V to Vramp for reset). However, as described above, the experimental I-V curves have different sections corresponding to the forward and reverse voltage ramps, both for the set and reset processes. Hereinafter, it is supposed that the forward voltage ramp denotes a functional variable and the reverse voltage ramp another variable. This situation is known as a functional bivariate case.
Let 1 ( * ), … , ( * ) be a functional sample from a bivariate random process denoted as ( * ) = ( 1 ( * ), 2 ( * )) , * ∈ [0,1]. The * ( * ) variable would correspond to the device current at voltage * . Note that ( * ) = ( 1 ( * ), 2 ( * )). We assume that the process ( * ) belongs to the Hilbert space 2 2 [0,1] with the usual characteristics defined above. Taking the basis expansion into account, each component of ( * ) can be expressed as follows, ∈ ℕ may be different for each functional variable. In matrix notation, Equation 10 adopts the following expression,  can be extended to the bivariate case (i.e., accounting for the I-V curve sections corresponding for the forward and reverse voltage ramps at once) as follows (see [Krzysko2019]), Making use of the previous developments, we could calculate the 2DVCf and 2DVCr coefficients described in the manuscript by means of Equation 9, and 2DVCt coefficient by means of Equation 13.

Supplementary Note 2: Compact modelling, a new variability implementation
We have used the Stanford Model (SM) [Guan2012, Chen2015, Jiang2016] to describe variability, accounting for the results of the methodology developed in the supplementary note 1. This model assumes filamentary conduction, and it calculates the evolution of a conductive nanofilament (CNF) in the dielectric; in particular, the CNF gap, g, i.e. the distance between the CNF tip and the electrode. In this manner, the device resistance evolution is calculated (the implementation is given in Verilog-A [Stanford model]) and the device operation can be described both in the quasi-static and transient operation regimes [Guan2012, Chen2015a, Jiang2016]. In the Chua's modeling formalism, the state variable is the CNF gap [Chua1971]. The abrupt shape of the Au/Ti/TiO2/Au device I-V curves suggests filamentary conduction; other TiO2-based devices have also been found to show filamentary conduction which is mainly due to the formation and rupture of nanometric paths (conductive nanofilaments) with a high concentration of oxygen vacancies that show ohmic-like charge conduction features [Lee2017, Carta2016].
where T stands for the CNF temperature [Guan2012]. The CNF gap variation (the model state variable) is calculated as follows, where ( ) accounts for random variations, ( ) a Gaussian distributed noise with zeromean and root mean square of unity [Jiang2016] generated at each time step. The dynamics of the CNF gap is given in Equation 16.
in which tox stands for the dielectric thickness, Eg is the activation energy for vacancy generation in the set process, Em is the energy barrier for oxygen ion migration that controls the reset process, v0 is the velocity dependent on the attempt-to-escape frequency, a0 the hopping site distance and VRRAM the applied voltage, which drops mainly at the CNF gap, γ is the field local enhancement factor that accounts for the polarizability of the material [Jiang2016].
In Figure S5a we plot the experimental curves from Figure 2a and the mean curve as obtained from Equation 4. Notice that after calculating the mean curve we have de-normalized it to plot it along the original experimental curves. A correct mean curve cannot be obtained without the previous normalization process, since we would sum currents limited by the compliance current (after the set event) with current values corresponding to voltages below the set voltages. The SM has been tuned to fit the mean curves obtained in Figure 4b and 4e. The parameter set employed is given in  Table 1 are in line with ref. [Roldan2023]. It can be seen in Figure S5b that the variability is not well modeled since the variation of the curves is low and the set and reset values are not well reproduced. Taking into consideration the limitations of this model for our devices, we modified the SM variability implementation to obtain curves in the interval marked with the dashed lines shown in Figures 4b and 4e in the main manuscript.
To do so, in addition to account for Equations 14 and 15, we included variations of the I0 parameter by generating random numbers in the [0,3] interval (uniform distribution) to redefine I0 parameter ( 0 ′ ), as described in Equation 17, 0 ′ = 0 + ( 0 * ) Furthermore, we allowed variations above g max for the gap in the set process that lead to higher variability in the HRS. The results obtained with the new implementation are shown in Figure  S5c. As can be seen, the variability corresponds much more accurately to the one characterized by the new statistical technique developed here. Figure S5. a Experimental I-V curves for the Au/Ti/TiO2/Au devices and mean (red line) curves (see Figures 4b, 4e in the main manuscript) for the set and reset processes. b SM simulated curves along with the mean (red line) curves, including the standard variability implementation in the SM. c Experimental and SM simulated I-V curves for the set and reset processes including the new model implemented to allow the variability characterized by the statistical procedure developed in this manuscript.