Non-Uniform Spline Quasi-Interpolation to Extract the Series Resistance in Resistive Switching Memristors for Compact Modeling Purposes

An advanced new methodology is presented to improve parameter extraction in resistive memories. The series resistance and some other parameters in resistive memories are obtained, making use of a two-stage algorithm, where the second one is based on quasi-interpolation on non-uniform partitions. The use of this latter advanced mathematical technique provides a numerically robust procedure, and in this manuscript, we focus on it. The series resistance, an essential parameter to characterize the circuit operation of resistive memories, is extracted from experimental curves measured in devices based on hafnium oxide as their dielectric layer. The experimental curves are highly non-linear, due to the underlying physics controlling the device operation, so that a stable numerical procedure is needed. The results also allow promising expectations in the massive extraction of new parameters that can help in the characterization of the electrical device behavior.


Introduction
Advanced mathematics can be applied to solve many engineering problems. In fact, centuries ago, many of the mathematical developments were born to solve specific problems that came up as the different engineering disciplines unfolded. However, many techniques developed currently in the context of research groups of mathematicians do not make their way to the fields of applied physics and of other technical subjects. In fact, a quick overview of the mathematical foundations taught in engineering colleges reveals content created mainly in the 19th century. In this manuscript, we present a multidisciplinary application where a clear and well-defined problem, contextualized in the field of electron devices (in particular, related to Resistive Random Access Memories, RRAMs), can be tackled with techniques based on state-of-the-art approximation theory tools, such as spline-quasi-interpolation and other procedures. These advanced tools come to the rescue when we are faced with complex numerical techniques. One of them is connected to accurate derivative calculation. If, for example, a procedure depends on an accurate derivative determination, a simple numerical calculation might not be enough, mostly if the measured device current curves show strong discontinuities, as is the case here with experimental current-voltage curves in RRAMs whose operation is based upon non-linear physical processes. In order to address the intricacies of the numerical procedure presented in this manuscript, we have proposed, in addition to the approximation theory approach, an algorithm to detect particular curve shapes.
There is currently a great interest in research and development activities devoted to RRAMs. These devices are compatible with conventional technologies in the electronics industry and present an exciting set of features for applications such as the possibility of being used as entropy sources for random number generation circuits, their suitability for non-volatile memories and their potential for neuromorphic computing circuits, such as hardware neural networks [1][2][3][4][5]. For the fabrication of neuromorphic circuits [1,3,4,6], these devices mimic biological synapses and allow for many compact circuits with lower power consumption than in previous technological approaches.
The experimental characterization of electron devices is essential in electronics. A great number of measurements are needed to extract parameters that describe a certain technology [7][8][9][10][11]. Later on, these parameters can be included in analytical expressions to calculate the device current, capacitances and other representative physical magnitudes that are required for their use in circuit design. The corresponding circuit simulation infrastructure is built upon compact models (a set of analytical expressions and technology-dependent parameters) that describe the devices. The parameter extraction routines and the analytical models can be complex; therefore, as stated above, advanced mathematics have to be employed to deal with modeling problems with guaranties.
In relation to the latter issue, it is important to highlight that several works have been presented on the subject. For instance, in [12], a numerical method based on smoothing splines was proposed to extract the threshold voltage in Metal Oxide Semiconductor Field Effect Transistors (MOSFETs). This method leads to the solution of a system of linear equations whose order increases with the number of intervals of the spline space considered. For this kind of device, a strategy relying on the use of a weighted essentially non-oscillatory method was described in [13]. MOSFET devices operate differently from RRAMs; however, some of the numerical issues that come up in the parameter extraction context are similar. A modeling problem for MOSFETs with a square shape was successfully tackled with advanced partial differential equation methodologies in [14]. For RRAMs, we can highlight the approach followed in [15], where state-of-the-art techniques for derivative calculation were employed.
We present here a mathematical technique to improve the extraction of the series resistance in RRAMs [16][17][18]. The series resistance is an essential parameter that needs to be obtained to represent the device correctly at the circuit level. A key issue in the corresponding extraction technique lies in the curve slope determination, and this can be performed by means of a derivative calculation. For this purpose, we study current curves (representing the set process of the device, where the device resistance drastically drops off due to the formation of an internal conductive filament) and go on with an in-depth modeling process based on non-uniform spline quasi-interpolation. We have employed related methodologies [19] based on the charge-flux domain; nevertheless, we have to go a step forward in these methodologies to address this new problem due to the strong non-linearity of our devices.
In Section 2, we shortly introduce the engineering problem and some details of the devices employed and measured. Section 3 is devoted to describing the mathematical foundations to deal with the numerical procedures needed. The final results are explained in Section 4, along with the main conclusions.

The Problem
The devices analyzed here were fabricated at the Institute of Microelectronics of Barcelona (CNM-CSIC) [18,20,21] and measured in our laboratory at the University of Granada. A scheme of the structures fabricated and measured is shown in Figure 1a. The devices were square-shaped with an area of 15 × 15 µm 2 . The bottom electrode (Wolfram) was grounded with a titanium layer below. A voltage ramp of (0.08 V/s) was applied to the top electrode (titanium with a TiN layer above) with a voltage step of 0.01 V. Hundreds of curves (resistive switching cycles) were measured in a successive manner. That is why our new technique has to be numerically stable and efficient, in order to be applied to many different current curves. For the measurements employed here, we employed 200 points per curve. This is reasonable since this means the use of voltage steps in the order of 0.01, where the main features of the devices under study are captured. We choose the curve with the maximum slope (in fact, it is a vertical region, shown in symbols) as the one that allows us to extract the correct series resistances, as depicted in [18]. Figure 1c shows an experimentally measured curve and those resulting when the variable corresponding to the X-axis is changed (from V applied to V modified ). The original applied voltage is corrected by the voltage in the series resistance; i.e., we reduce V applied by a factor calculated as current × R series (see Figure 1b). When we do so, we obtain a new voltage, V modified , the one we employ to replot the current curves for different series resistances (see Figure 1c). Among the set of curves obtained for different series resistances, we select the series resistance corresponding to the curve with the vertical section (plotted with symbols in Figure 1c) [18]. Figure 2 shows several original (black) and modified (red) curves, once the series resistance has been extracted. Observe that when the effects of the series resistance have been eliminated, the current versus modified voltage can show a section of negative slope (just after the V TS 2 points) [16][17][18]. We define this new parameter as the turning point of the curve for the region of the negative slope, when the curves turn left in the plot. The threshold set voltage (V TS ) is defined as the modified voltage corresponding to the vertical region in the current curve, see Figure 2. Once the vertical slope is detected and the corresponding series resistance is obtained, the algorithm to extract V TS and V TS 2 needs a smooth approximant that accurately represents the current-voltage curve.

Algorithms Description
Let V i and I i be the voltage and current applied, at time t i , respectively, i = 0, . . . , n. We assume here a non-uniform approach because, as a general rule, the measurement instruments do not always provide a uniform set of data. Therefore, for the sake of generality, a non-uniform quasi-interpolation scheme is followed. The proposed procedure for extracting the parameters of interest associated with the curve V applied − I is structured in two parts. In the first one, a procedure is established that provides the series resistance R series such that the modified data {(V i − R series I i , I i ), i = 0, . . . , n} shows a nearly vertical segment. It also provides the set of points that determine that segment. In the second part, the quasi-interpolation spline on non-uniform partitions is used to deal with the modified curve V modified − I data associated with the series resistence R series in order to extract other parameters of interest.

Estimating the Series Resistance
In this section, we will propose the algorithm to estimate R series and the threshold set voltage V TS . In the V modified − I curve, V TS corresponds to a vertical set of data for the value R series [18].
In order to understand the behavior of the V modified − I curve, firstly, we represent the measured values of a specific cycle (we consider a cycle as a set process followed by a reset process; therefore, the device resistance is changed from high to low and the other way around in each cycle that is repeated many times, 1000 in our case, in what is called a resistive switching series in the RRAM characterization process). A series resistance R = 0 Ω is applied to obtain the black plot in Figure 3, the original measurements. Then, the initial values are modified with R = 28.25 Ω to obtain the red plot with symbols. Two other curves are represented, for which R = 13 Ω (green plot) and R = 36 Ω (pink plot).
As shown in Figure 3a, the cycle starts and ends at (0, 0). We observe that the vertical segment is approximately in the center of the cycle, so that, for our measurements, the beginning and the end of the cycle are neglected [18], around 25% of the data in each part. The admissible part of the cycle will begin at index i 1 , the index of the first point in the 25% limit to end at i n , and the index of the last point in the 25% limit. The behavior of the curves modified by using a resistance R is quite similar, so this procedure is reasonable, see more details in [18]. To detect the vertical set of data, we work with I − V mod curves, where V mod := V applied − R × Current, so we will detect a horizontal set from which R series will result (see Figure 3b). We will consider constant segments formed by, at least, m points, with m ≥ 10.
To find if a given I − V mod cycle (with a fixed R value, starting from R = 0) has a horizontal subset of data, we compute the least square constant for data starting at I i 1 , V modi 1 and ending at I i m−1 , V modi m−1 , and the least square constant, V ls m+1 (i 1 ), for data starting at I i 1 , V modi 1 and ending at (I i m , V modi m ). Define also m max = m. Given a tolerance , andV ls (i 1 ) being the mean value of the computed V ls (i 1 ) starting If we have not found a constant segment in the previous process starting with i 1 , we increase this index by 1 and check again for a constant segment. The process ends when we find a constant segment or when the starting point i i is outside the valid data range, i.e., in the 25% ending data limit.
We stop increasing R when we find negative values of V mod . Now we have an interval [R min , R max ] of values of R that provide a constant segment of data in the I − V mod circle (as shown in Figure 3b, for the red curve) within the given tolerance . We have to select the R value that gives the optimal segment, the "more" constant one.
In order to estimate the optimal value R series of R, we first select the maximal interval (i st , i end ) among all the intervals found for the different values of R we have checked.
For this maximal interval, we will perform a bisection search of the value R series in [R min , R max ] applied to the mean square error to compute the value of R that provides the smallest E(R). In order to achieve that, we compute the values E(R min ), E(R max ), and we take a new interval [R min , R max ] as the one with one of the endpoints at the midpoint of the interval [R min , R max ] and the other endpoint as the value of R min or R max with smaller E(R). We end the bisection process when the length of the interval (or the value E(R)) is less than a certain tolerance. The application of the explained process to the cycle that appears in Figure 3a (with R = 0, black curve) provides the resistance R = 28.25 (red curve).
The algorithm is summarized, step by step, in the following lines: 1. Take R = 0.

2.
If any value for V mod (R) is negative, go to point 12.

3.
Take a point i st = i 1 starting 25% from the beginning of the data.

4.
Take an ending point i end = i m−1 , with fixed m ≥ 10. If i end is outside the 25% limit, go to stage 11. 5.
Compute the mean-square constant V ls i end that aproximates (I k (R), V modk (R)) with k = i st , . . . , i end . 6. Compute whereV ls is the mean value of the computed V ls starting at index i 1 .

7.
If i end > i m−1 and E ∞ < , the data with indexes between i st and i end form a constant segment. 8.
If i end = i m−1 or the segment (i st , i end ) is constant, increase i end = i end + 1 (ending when the value is outside the 25% limit) and go to point 5. 9.
If the segment (i st , i end−1 ) is constant, but E ∞ > , for R, the segment (i 1 , i end−1 ) is a potential candidate for the final horizontal segment. Save the values associated with R. 10. If the segment (i st , i end ) is not constant, increase i st = i st + 1 (ending when the value is outside the 25% limit) and go to stage 4. 11. Increase R = R + 1 and go to stage 2. 12. For all the values of R with constant segments, take the reference interval as the greatest one among the saved segment candidates. 13. Perform a bisection search between R min and R max applied to the mean square error given by Equation (1) to compute the value of R that provides the smallest E(R).

Non-Uniform Quasi-Interpolating Splines for Parameter Extraction
Once the resistance R that provides the appropriate modified voltages has been determined, the data {(V i − R I i , I i ), i = 0, . . . , n} should be approximated in order to extract other parameters of interest. The variable change introduced demonstrates the need of the non-uniform approach followed here evident. Therefore, now we propose a general method to construct high-precision spline approximants to a function f defined on an interval I := [a, b], from its values f (t i ) at non-equally spaced points t i such that a := t 0 < t 1 < t 2 < · · · < t n−1 < t n := b. (2) It will be applied to the voltage and current functions to define approximating splines from their values at times t i . Now, we introduce the space where the approximating splines will be defined and recall some results [22,23]. Definition 1. Let n, d ∈ N, and let ∆ be a knot sequence satisfying (2). The spline space S d (∆) of order d + 1 is the space of all C d−1 (I)−piecewise polynomial functions whose restrictions to each sub-interval [t i , t i+1 ], i = 0, . . . , n − 1 is a polynomial of degree at most d.
It is well-known that S d (∆) is a linear space of dimension n + d. To get a good basis of this space, an extended partition is needed. Among all options, we choose the one with multiple end knots (see Figure 4): (3) where [z 0 , . . . , z k ]g and (·) r + stand for the divided difference of g at z 0 , . . . , z k and the truncated power of degree r, respectively. Then, N i,d 1≤i≤n+d forms a basis for S d (∆) with Since monomials up to the degree d belong to the space S d (∆), they are expressed in terms of the B-splines. These representations are essential to define the approximating splines we want to use to approximate a given function. They are obtained as a consequence of Marsden's identity [23] (Th. 4.21).

Proposition 2.
For any non-negative integer r ≤ d, the monomial m r (t) := t r can be written as a linear combination of the basis of B-splines as follows: where θ (r) and the classical symmetric functions are In particular, from (4) and (5), we obtain θ As a consequence, the operator S : C(I) → S d (∆) defined as is exact on P 1 , the space of polynomials of degree less than or equal to 1. That Schoenberg operator provides an approximating spline for f , but only an approximation order equal to O h 2 is obtained, where h := max 0≤i≤n−1 h i and h i := t i+1 − t i . Note that, when the t i 's are equally spaced values obtained by dividing the interval I into n equal parts, then t i = a + i h with h = (b − a)/n and Therefore, for an even d, the Greville's abscissa θ (1) i is the mid-point of an interval.
However, when d is odd, θ (1) i is one of the knots. Since the function to be approximated is only known at the knots, in the general case, we will consider an odd degree d = 2k + 1, k ≥ 1 and construct operators Q d : C(I) → S d (∆) having the form and yielding an approximation order O h d+1 . This is achieved by imposing the exactness of Q d on the space P d of all polynomials of degree at most d.

Lemma 1. Q d is exact on P d if and only if
Proof. Q d is exact on P d if and only if Q d m r = m r , 0 ≤ r ≤ d. On the one hand, by (6) it holds that On the other hand, by (4) we have From now on we will assume that d is odd. Figure 5) so that the linear functionals µ i,d will be determined by the values of f at d + 1 knots lying in the support, i.e., Then, d − i + 1 additional knots are needed to define Analogously, for n − d ≤ i ≤ n − 1 we define Theorem 1. There exist unique values α i,j such that the quasi-interpolation Q d given by (6) and the linear functionals defined by (8)-(10) are exact on P d .
Proof. Let us suppose d + 1 ≤ i ≤ n. By (7), Q r is exact on P d if and only if µ i,d (m r ) = θ (r) The Vandermonde matrix of this system of linear equations is non-singular as the involved knots are pairwise distinct. Therefore, the claim follows in this case.
The proof in the case 1 ≤ i ≤ d (resp. n − d ≤ i ≤ n − 1) is similar, but now the linear functional given by (9) (resp. (10)) is involved, so that the matrix of coefficients is of Vandermonde type based on the knots t 0 , . . . , t d (resp. t n−d , . . . , t n ).
In the cubic case, for s + 1 ≤ i ≤ n, it holds and a straightforward calculation leads to the following result on the unique quasi-interpolation operator Q 3 defined from linear functionals based on point evaluations at the knots and yields the optimal approximation order. Proposition 3. The coefficients of the unique quasi-interpolation operator Q 3 exact on P 3 and given by (6), defined by the linear functionals in (8)- (10), are given by the following expressions: The resulting operator will be used to approximate the voltage-current curves. Note that after solving the system ensuring the exactness of the quasi-interpolation operator, an approximating spline is immediately available for the function being considered, and no specific system needs to be solved to determine it. Only the values of the function at the knots determine the coefficients of the linear combination of B-splines.

Results and Discussion
To show the effectiveness of the proposed procedure to treat experimental data corresponding to the considered non-linear devices, we will make use of two cycles. For the first one, we will apply the algorithm of series resistance extraction that produces a modified set that gives rise to a discrete curve with an almost vertical section. It has provided a resistance of R = 28.97 ohms. Both the original and the modified curve are shown in Figure 6a. For the second cycle (see Figure 6b), the algorithm gives a series resistance equal to R = 21.25 ohms.
Therefore, to estimate the points V TS 2 shown in Figure 6a, we have taken into account that the modified curves show different regions. At low voltages, the current is low because the conductive filament is not formed; therefore, the curve slope is also relatively low. At the onset of the low resistance state, the conductive filament is formed and there is a sudden current rise, as expected. This constitutes a second region where negative slopes could be seen. The transition of these operation regions is marked by V TS 2 . In some cases, the voltage in the series resistance increases considerably, and the modified voltage is reduced, even if the applied voltage goes on increasing.
For each cycle, the C 2 cubic quasi-interpolants QV mod and QI to voltage and current, respectively, are computed. The derivative of I as a function of V is estimated from the derivatives of QI and QV mod as functions of t to yield dI dV QI (t) QV mod (t) , and the latter is computed to determine the first point at which the slope is negative. It is an estimate of V TS 2 . For the cycle considered in Figure 6a, it has been obtained that the slope is negative at t = 7.69 and V TS 2 (0.43131V, 0.00166A). For the one in Figure 6b, t = 8.17 and V TS 2 (0.43158V, 0.00457A). These results are in good agreement with the results shown by the only information available, namely voltages and currents measured at non-uniformly spaced times.
See in Figure 6 that the quasi-interpolant is built correctly both for the original and modified curves. In addition, the vertical section of the curve is correctly detected to extract the corresponding series resistance. For the cycle 2 (resp. 106), the point V TS is attained at t = 11.18 (resp. t = 9.48) and V TS (0.27014V, 0.01770A) (resp. V TS (0.43319V, 0.0089A)).

Conclusions
A new methodology based on an approximation theory approach is presented to deal with the extraction of the series resistance and set transition voltage in resistive memories. A robust procedure is presented in order to develop an algorithm that can be programmed to analyze hundreds of consecutive current-voltage curves in an automatic way.