Robust Simulation of a TaO Memristor Model

This work presents a continuous and differentiable approximation of a Tantalum oxide memristor model which is suited for robust numerical simulations in software. The original model was recently developed at Hewlett Packard labs on the basis of experiments carried out on a memristor manufactured in house. The Hewlett Packard model of the nano-scale device is accurate and may be taken as reference for a deep investigation of the capabilities of the memristor based on Tantalum oxide. However, the model contains discontinuous and piecewise differentiable functions respectively in state equation and Ohm’s based law. Numerical integration of the differential algebraic equation set may be significantly facilitated under substitution of these functions with appropriate continuous and differentiable approximations. A detailed investigation of classes of possible continuous and differentiable kernels for the approximation of the discontinuous and piecewise differentiable functions in the original model led to the choice of near optimal candidates. The resulting continuous and differentiable DAE set captures accurately the dynamics of the original model, delivers well-behaved numerical solutions in software, and may be integrated into a commercially-available circuit simulator.


Introduction
The memristor and its applications represent one of the most interesting fields of research.The circuit theoretic introduction of the memristor traces back to 1971, when L. Chua [1] presented a new fundamental two-terminal element where a nonlinear relation links the time integrals of voltage and current, namely the flux across and the charge through the device respectively.Thanks to the crucial influence of the history of dynamics of the memristor on its current behavior, this device exhibits memory capability.This is the reason why L. Chua coined the term memory resistor, abbreviated as memristor, to denote it.Before the pioneering work from Chua, circuit theory classified only three fundamental electrical bipoles, namely the resistor, the inductor, and the capacitor.As a result, the advent of the memristor [2] represented a sensational breakthrough in electrical engineering.Nonetheless, the topic of memristors [3] mostly attracted the attention of theoretically-inclined researchers for the following three decades, until, in 2008, a group of Hewlett Packard (HP) engineers, supervised by S. Williams [4], uncovered the emergence of fingerprints [5] of memristive behavior [6] in nature, while studying the switching property of a double-layer nano-scale film based on Titanium dioxide and sandwiched between two metallic contacts [7].Since then the numbers of academic papers concerning the memristor has been increasing exponentially.The practical applications of the memristor are various.The most interesting field for the industry concerns the design of non-volatile memories.Real memristors have nanoscale dimension, consume negligible power, and may switch at very high speed.Furthermore a two-dimensional memristor array1 may be laid directly on top of CMOS circuitry, supporting the continuing nano-electronics industry trend to reduce integrated circuit (IC) area, and leveraging upon the compatibility between the well-established CMOS technology and memristor manufacturing process.Memristors are much more than "mere" memory cells.They also exhibit computing capability as demonstrated in numerous studies [10].They may pave the way to the introduction of novel computer architectures where processing unit and memory are no longer physically separated, thus clearly standing out from the classical Von Neumann computing approach, and paving the way to the era of memory intensive computing [11].Moreover, the memristor displays a close resemblance to a neural synapse.In fact the history-dependent conductance of the memristor may mimic the plasticity of the synaptic weight, which is finely tuned upon changes in the ionic flow through the synapse.Memristors may support various synaptic rules at the origin of neural learning, including habituation, long-term potentiation (LTP) [12], Hebb's and anti-Hebb's rules [13], spike-timing-dependent-plasticity [14]- [15] (based on either pairs or triplets of spikes), and even spike-rate-dependent-plasticity.The availability of CMOS circuits for the emulation of neurons, and the possibility to stack over them multiple layers of memristor arrays, may lead to the development of a bio-inspired computing engine with massively-parallel signal processing power, capable to reproduce the complex dynamics of the neural networks of the human brain.Finally, there exist memristors capable to amplify infinitesimal fluctuations of energy, e.g. the Niobium dioxide-based nanostructure manufactured at HP labs [16].The local activity of these memristors and the complex dynamics which may emerge in circuits employing them (for example the oscillatory behavior of a simple locally active memristor circuit built at NaMLab is investigated in [17]) may lead to the development of new circuits [18] which could outperform state-of-the-art electronic systems or complement their functionalities.This wide plethora of interesting applications urgently requires the developments of accurate mathematical models in order to uncover the full potential of memristors in the electronics of the future.An accurate mathematical model for a TaO memristor nano-device, fabricated at HP labs, was recently introduced [19] (Section 3 revisits its formulation, after a preliminary review of the latest memristor classification, given in Section 2).In view of the excellent performance of the HP nano device as non-volatile memory cell [20] the study of the relative model is timely.In this paper investigations are focused on the derivation of a mathematical representation of the original model suited for robust numerical simulations in software.In particular, past investigations have revealed the emergence of convergence issues in numerical integration of memristor differential algebraic equations (DAE) containing discontinuous and/or piecewise differentiable functions [21].Therefore, the computation of a computer-based solution to the original TaO memristor model may be facilitated through the introduction of suitable continuous and differentiable approximations to the discontinuous and piecewise differentiable functions appearing in state equation and Ohm's based law respectively.A thorough investigation of a set of candidate kernels led to the selection of the most appropriate approximating functions (see Section 4 for details), and resulted into an accurate mathematical description suited for robust numerical integration in software, as demonstrated in Section 5.The model may be integrated in circuit simulators available on the market, e.g.LTspice [22].A circuit implementation of the continuous and differentiable DAE set, coded in LTspice version IV, and presented in Section 6 in the form of a netlist, may be of interest to the circuit designer eager to explore the design opportunities the HP TaO resistance switching memory [23] may open up in the field of nano electronics.

Brief Review on Memristor Classes
A n th -order extended memristor is defined as [24] where u ∈ R and y ∈ R are input and output respectively (if one is the voltage v m across the memristor, the other is the current i m through it), x ∈ R n denotes the state vector, and f(•, •) : R n × R → R n stands for the vector field governing the state evolution with time.Importantly, H(•, •) : R n × R → R is a scalar state-and input-dependent function satisfying the constraint H(x, 0) = 0 irrespective of x, and representing the memristance M(x, u) or memductance W (x, u) under current or voltage input respectively.The memristor model analysed in Section 3 falls into the class of voltage-controlled extended memristors.Other exemplars of the extended memristors, based upon a suitable combination between a static nonlinear two-port and a dynamic one-port were presented in [25].A subclass of the class of memristors described by equations ( 1)-( 2), known as generic memristors [24], is expressed by the following DAE set: Here, under current (voltage) input the memristance (memductance) is a function of the state only.An example of generic memristors are the Sodium and Potassium ion channel memristors [26].With reference to first-order systems, a very important subclass of generic memristors consists of ideal memristors, where the state and Ohm's based equations assume the following form: Here the state is either the charge q m = t −∞ i m (t )dt through or the flux ϕ m = t −∞ v m (t )dt across the memristor, therefore the state evolution function is simply expressed by current or voltage input respectively.An ideal memristor may alternatively be described through the nonlinear constitutive relationship between charge and flux, namely under current input, and under voltage input.From the class of generic memristors a further subclass of memristors, known as ideal generic memristors, may be identified.Without loss of generality, referral to first-order systems is made here.The DAE set of these circuit elements may be cast as The name of these circuit elements originates from their equivalence to ideal memristors under suitable conditions.In fact, referring to the current-controlled case (the voltage input scenario descends by duality), the equations of an ideal generic memristor (i.e. ( 9)-( 10) with H(x) = M(x)) may be reduced to the constitutive equation of an ideal memristor (i.e. ( 7)) provided there exists a one-to-one differentiable function x = x(q m ) such that It follows that an infinite number of ideal generic memristor siblings [24] descends from each ideal memristor.As a result, the class of ideal generic memristors includes all the ideal memristors as a subset.

Model
The mathematical model of the voltage-controlled TaO x -based memristor nano-device fabricated at HP Labs and reported in [19] is expressed by the following DAE set which falls into the class of equations ( 1)-( 2), where the system order is n = 1, the input and output are respectively given by u = v m and y = i m , while H(x, u) = W (x, v m ) is the state and input-dependent memductance function, expressed by The memristor state x represents the volume fraction of the conductive channel.As a result, it has a limited domain of existence, namely the closed set [0, 1].The lower and upper bound of this set respectively refer to the fully-insulating and fully-conductive state of the nanostructure.In ( 13) step(•) stands for the step function, defined as

2
, where sign(•) denotes the sign function.The values of the model parameters [19] are given in Tab. 1.We shall use these values in all the investigations described in this manuscript.
Tab. 1. Values of the parameters in the model equations of the HP TaO x -based memristor.

Model Approximation
This section presents the key result of this research contribution.Recent investigations have highlighted the negative impact the presence of discontinuous and/or piecewise differentiable functions may have on the convergence properties of a numerical solver to DAE sets pertaining to memristor models [21].As a result, the replacement of discontinuous and piecewise differentiable functions with appropriate continuous and differentiable approximations may facilitate the numerical integration of the models.With reference to the modelling equations of the HP TaO memristor, it is evident that two are the functions under our zooming lens, particularly the discontinuous step function in the state equation (13), and the piecewise differentiable modulus function in the memductance function (see (15)) appearing in the Ohm's based law (14).The modulus function, acting on the memristor voltage v m and denoted as |v m | in (15), may be replaced by differentiable approximations falling into the following class: where ρ is a parameter defining the concavity of the proposed kernels at v m = 0 V. Figure 1 shows illustrative plots of modulus function and its differentiable approximations for values of ρ set to 10 1 , 10 2 , and 10 3 .In order to pick up a suitable value for ρ we focus on the Ohm's based law (14) only.In fact we want to check what is the effect of this parameter in modelling a scenario where the input across the device is applied for a very short period of time and takes rather small values, so that any state change over time may be neglected.As a result, under these conditions the state keeps approximately equal to its initial value x(0) = x 0 , and the dynamic equation ( 13) is unnecessary to study the current response to the voltage excitation.
The value of ρ has a major influence only when the state is close to the lower bound, where the exponential term in (15) dominates, and for values of the memristor voltage in the neighborhood of 0 V, where the function g ρ (v m ) most deviates from the modulus function for any ρ (see Fig. 1).Each plot in Fig. 2 shows a comparison between the currentvoltage curves descending from the original Ohm's based law (14) (adopting the modulus function in the memductance expression) and the proposed differentiable variant (which is calculated with the approximating function in (16) in place for the modulus function).The first (latter) curves are drawn through solid (dashed) lines.The parameter of the differentiable kernel g ρ (•) is set to 10 1 in plot (a), 10 2 in plot (b), and 10 3 in plot (c).Four state values close to the lower bound are considered in each plot (see the caption in Fig. 2 for details).For better visualization purposes only a negative memristor voltage range is shown.The approximated numerical results best agree with the ones pertaining to the original model for ρ = 10 3 (see plot (c)).The characteristic of Fig. 3 is obtained using (14) with |v m | in plot (a) and with g ρ (v m ) with ρ = 10 3 in plot (b), and sweeping the memristor initial state x 0 in the closed set Ψ = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1}.The agreement between the two pictures is a sign that our choice is near optimal.Thus in the remaining part of the paper we keep this value for ρ.
Further, the step function with argument v m , denoted in (13) as step(v m ), may be approximated by continuous kernels falling within the following class: where k controls the slope of the proposed kernel around v m = 0 V. Fig. 4 shows exemplary graphs of f k (v m ) for values of k ∈ {30, 40, 50} together with the shape of the step function.
-0.The excitation drives the series connection between the memristor and a resistor of value R = 70.1 Ω.The initial memristor state x 0 is set to 0.065.
Fig. 6(a) shows the numerical simulation result obtained by replacing step(v m ) in state equation ( 13) with function f k (v m ), as given in (17), for k ∈ {30, 40, 50} (here g ρ (v m ) with ρ = 10 3 was used in place for |v m | in the memductance function (15) appearing in Ohm's based law (14)).The curve referring to k = 50 better matches the characteristic observed using the original model equations ( 13)-( 14) with the discontinuous and piecewise differentiable functions (see Fig. 6(b)).Thus the value for k in( 17) is set to 50 in the study to follow.

Approximation Accuracy
Let us confirm the appropriateness of the approximations chosen for the discontinuous and piecewise differentiable functions.Sweeping the period T of the above mentioned asymmetric triangular input waveform in {1, 10 −2 , 10 −4 , 10 −6 , 10 −8 } s, and keeping all the other parameters of the voltage source v as well as x 0 as in the simulation of Fig. 6, the memristor steady state current-voltage curves observed in the circuit of Fig. 5(a) are reported in Fig. 7. Plots (a) and (b), respectively refer to the original model equations ( 13)-( 14), and to their approximations where f k (v m ) with k = 50 and g ρ (v m ) with ρ = 10 3 respectively replace step(v m ) and |v m |.Let us consider yet another simulation scenario.Here the same circuit of Fig. 5(a) with R = 70.1 Ω is taken into exam, but the voltage waveform exciting the system is a symmetric square wave of amplitude 0.55 V (see Fig. 8(a) for the case where the input period is equal to 1 s).The initial condition is kept equal to x 0 = 0.065.Distinct values are given to the period T of this input signal in order to investigate the frequency dependence of the memristor response.In plots (a) and (b) of Fig. 9 original model and its continuous and differentiable approximation are respectively adopted to reproduce the reduction in the lobe area of the memristor steady state pinched hysteresis current-voltage loop observed under increasing frequency for the square wave excitation under exam.The numerical results shown in plots (a) and (b) in each of Figs. 7 and 9 agree with each other, providing a further proof for the accuracy of the proposed continuous and differentiable approximation to the original discontinuous and piecewise differentiable model equations.Another excitation signal of interest is the sine wave.Driving the circuit of Fig. 5(a) with a signal of this kind, with amplitude equal to 0.55 V (Fig. In the investigations on the nonvolatile memory capability of the nanostructure, pulse train-based excitations are typical.Before concluding this section, it is therefore opportune to evaluate the efficacy of the proposed continuous and differentiable approximation to the original HP DAE set in a scenario where voltage source v delivers a pulse train in the circuit of Fig. 5(a) with R = 70.1 Ω.Let the input voltage waveform, evolving with time as depicted in Fig. 11, consist of a succession of 4 pulses.The first two have positive polarity, while the amplitude of the second two is negative.In the pair of positive (negative)-polarity pulses, the first one increases (decreases) the state of the device (here the initial state x 0 is kept equal to 0.1), while the second one allows an indirect measurement of the memductance value through the reading of the small memristor current it gives rise to.The amplitude of the read pulses V 0r is set to 0.1 V, while that of the write pulses V 0w is swept from 0.8 V to 1.1 V with constant increment 0.1 V.The width of read and positivepolarity write pulses is set to 5 ns, while the negative-polarity write pulses are applied for a period of time equal to 20 ns.The rise and fall time is taken as 1 ns for all pulses.Fig. 12 reports the time waveform of the memristor state resulting from the pulse train excitation of Fig. 11.A strong dependence of the dynamic behavior of the state upon polarity, amplitude, and duration of a pulse may be observed.The current flowing through the nano device evolves in time as illustrated in Fig. 13.Plots (a) and (b) in both Figs. 12 and 13 refer to the numerical integration of original and approximated models respectively.Their accord further validates the accuracy of the proposed continuous and differentiable approximation to the discontinuos and piecewise differentiable DAE set expressed by ( 13)-( 14).

LTspice Implementation
The continuous and differentiable model is suited for integration into commercially-available software packages with integrated circuit design emphasis, e.g., LTspice IV [22].In this section we present a simple circuit implementation of the approximated TaO x memristor model.Within the sub-circuit, the Ohm's based law is implemented through the use of a voltage-dependent current source G res , which couples the two access terminals and defines the current flowing between them.The latter represents the memristor current i m given in (14).The two control voltages for this current source are the potential difference between the access terminals plus and minus of the sub-circuit, modeling the memristor voltage v m , and the voltage drop across a linear capacitor C int (i.e.voltage at node Y ), which models the memristor state x.The state equation ( 13) is modelled by the charging rate of this capacitor driven by yet another voltage-dependent current source G Y placed in parallel to it.This current source, controlled by the same voltage inputs as G res , defines the current flowing through the capacitor, which is specified by the memristor state evolution function, reported on the right hand side of (13).A large auxiliary resistance R aux is also inserted in parallel to the capacitor to prevent convergence issues in the circuit simulator.The initial memristor state x 0 is modelled by the initial condition on the capacitor voltage.The netlist of the circuit implementation of Fig. 14(b), coded according to the guidelines specified in the LTspice IV guide [27], is reported below.Here the continuous and differentiable approximations to the discontinuous and piecewise differentiable functions of the original model, proposed in Section 4, are defined through the use of suitable functions.The netlist is ready to use, may be easily adapted to similar circuit simulators, and may be of interest to the circuit designer eager to explore the opportunities the promising TaO memristor may open up in the field of nanoelectronics.

Conclusions
This paper proposes an accurate continuous and differentiable approximation of a TaO memristor model [19] suited for robust numerical simulations in software.Discontinuous and piecewise differentiable functions respectively appear in state equation and Ohm's based law of the original Hewlett Packard model of the nano-device.A detailed investigation of classes of possible continuous and differentiable kernels for an accurate approximation of the discontinuous and piecewise differentiable functions of the original model led to the selection of the most appropriate candidates, resulting in a differential algebraic equation set which facilitates numerical integration in software.The proposed continuous and differentiable mathematical description is ideally suited for integration into commercially-available circuit simulators.The code of a circuit implementation of the approximated model was also presented.It may be directly used in LTspice version IV to investigate the full potential of the Hewlett Packard nanostructure in electronics.This work complements other important contributions on techniques for robust numerical simulation of memristor models [28]- [29].

Fig. 3 .Fig. 4 .
Fig. 3. Memristor current-voltage curves under excitation scenarios unable to trigger any state change from the initial condition.Here the initial state is swept from 0 to 1 in uniform steps of 0.1.For better visualisation purposes, only a negative memristor voltage range is shown.Plots (a) and (b) are respectively obtained using the original memductance expression with the modulus function and its differentiable approximation where |v m | is replaced by g ρ (v m ), as given in equation (16), with ρ = 10 3 .The two plots agree for all the state values.

Fig. 5 .
Fig. 5. (a) Test circuit with series resistor.(b) Time waveform of the asymmetric triangular voltage waveform of period T = 1 s driving the circuit in plot (a).

Fig. 6 .
Fig. 6.Steady state pinched hysteresis loops in the memristor current-voltage plane under the asymmetric triangular excitation of the circuit of Fig. 5(a), as illustrated in Fig. 5(b).(a) Numerical result using continuous transition function f k (v m ) in place for the step function in (13) for values of k in {30, 40, 50}.The modulus function |v m | is approximated with function g ρ (v m ), as given in (16), with ρ = 10 3 .(b) Numerical result using step and modulus functions, and its agreement with the curve for k = 50 in plot (a).

Fig. 7 .
Fig. 7. Frequency dependence of the memristor steady state current-voltage curve under periodic excitation of the circuit of Fig. 5(a) where R = 70.1 Ω and v is an asymmetric triangular voltage waveform with maximum and minimum voltages respectively equal to 0.8 V and −1.2 V. Numerical results from the original model (a) and from its continuous and differentiable approximation (k = 50, ρ = 10 3 ).

Fig. 8 .
Fig. 8. Time evolution of symmetrical waveforms at the input to the circuit of Fig. 5(a).(a) Square wave.(b) Sine wave.

Fig. 9 .
Fig. 9. Frequency dependence of the memristor steady state current-voltage curve under square wave excitation of the circuit of Fig. 5(a) where R = 70.1 Ω and v is a symmetric square voltage waveform with amplitude 0.55 V. Numerical results from the original model (a) and from its continuous and differentiable approximation (k = 50, ρ = 10 3 ).

4 Fig. 10 .Fig. 11 .
Fig. 10.Memristor steady state pinched hysteresis loops in the current-voltage plane under sine wave excitation of the circuit of Fig. 5(a).The resistance R is set to 70.1 Ω, the initial memristor state x 0 to 0.1, while the amplitude of the input voltage v is chosen as 0.55 V.The different curves pertain to distinct input periods (see the legend).Plots (a) and (b) respectively refer to original model and to the proposed continuous and differentiable approximation.In the latter k = 50, and ρ = 10 3 .
Figures 14(a)-(b) respectively show memristor circuit symbol and electronic realization of the proposed continuous and differentiable DAE set.The memristor is implemented as a sub-circuit with two access terminals (plus and minus).

3 Fig. 12 .
Fig. 12. Memristor state response to the pulse excitation of Fig. 11.Plots (a) and (b) respectively refer to original HP TaO model[19] and to its approximation as proposed in this paper.

Fig. 13 .
Fig. 13.Current flowing through the memristor in response to the pulse train excitation of Fig. 11.The results obtained through numerical integration of original and approximated models are respectively shown in plots (a) and (b).