Takacs Model of Hysteresis in Mathematical Modeling of Memristors

In this paper, the mathematical modeling of memristor via Takacs model of hysteresis is presented along with a modification of this model tailored to describe the asymmetric hysteresis loop and first order reversal curves. In particular, it is shown that there is a class of differential equations of the Duhem model of hysteresis where every member of the class could play a role of the state equation of memristor. Within this class of Duhem differential equations, there are two distinct subclasses: one corresponding to the Takacs model and the other one corresponding to the state equations of the memristor model with the Biolek window function of various degrees p. These two subclasses have a non-empty intersection, which contains the state equation of the memristor model with the Biolek window function for p = 1. To demonstrate the proposed approach, three examples are presented.


Introduction
The hysteresis is a nonlinear phenomenon, which emerges in numerous fields, including physics, electronics, materials sciences, biology, mechanics, economy, etc. [1]. It was first observed in ferromagnetic materials, and later in smart materials such as piezoelectric materials, electroactive polymers, magnetostrictive materials, shape memory alloys, etc. [2].
Due to the long history of the subject and the omnipresence of hysteresis phenomena, there exist a considerable number of models and published books on the modeling of hysteresis [3]- [12]. The main advantage of mathematical modeling is in its ability to describe experimental data via analytical expressions, which in turn provides a simplified and more efficient analysis of systems exhibiting hysteretic behavior. There is no model of hysteresis loop (HL) capable to grasp all necessary features of the object to be modeled [12].
The appearance of pinched hysteresis loop (PHL) in the voltage-current plane is one of the fingerprints of memristive behavior [31], [32], [33]. A comprehensive overview of fingerprints and the history of the subject are presented in [34]. Odd-symmetric self-crossing PHL is recognized as the signature of ideal memristor [35]. For various degrees of degeneracies, the classification of PHLs on self-crossing and touching has been further refined in [36]. Dependence of the type of PHL on the frequency content of the state variable is analyzed in [37]. Under the assumption that constant charge can be delivered to an ideal memristor within the half-period, the influence of the frequency of sinusoidal excitation on the area of PHL is analyzed in [38]. Apart from memristor, memcapacitor and meminductor, the other nonlinear elements in Chua's table also exhibit PHLs in the appropriate planes, as shown in [39]. Non-memristive elements having PHLs in the voltage-current plane are considered in [40]. Physical interpretation of lobe area (LA) of PHL is discussed in [41]. The computation of LA via the time-domain integration in the voltage-current plane has been studied in [42], [43] and [44]. Also, the computation of LA from the memristance vs. state map of the ideal and ideal generic memristors is considered in [45].
A simple integrator-multiplier model of PHL has been introduced in [46]. Theory of Lissajous figures were applied to the creation of PHL and experimentally verified in [47]. A graphical modeling of PHL that exhibits fingerprints of a memristor is presented in [48].
Mathematical model for the major and minor PHLs of memdiode (diode with memory) is proposed in [49] and then modified in [50] by introducing the rate-dependent state variable. In both models, the branches of hysteron are logistic curves, expressed via exponential functions. The simple relation that exists between the exponential and hyperbolic tangent function implies that the branches of logistic hysteron can also be expressed via the tangent hyperbolic function.
Hyperbolic-type memristor with memductance expressed via the hyperbolic tangent function of state variable is proposed in [51]. Although not explicitly reported, HL appears in the memductance-flux plane of this type of memristor. (As usual, by flux/charge we meant the time integral of voltage/current.) Similarly, HL in the memductance-flux plane of TaO x memristor can be observed in simulations with the model described in [52].
In addition, HL has been reported for spintronic memristor in memristance-flux plane [53], for ferroelectric memristor (along with the so-called reversal curves) in memristance-voltage plane [54], and for quantum point contact memristor in resistance-voltage plane [55].
HLs are measured for spintronic synapses for ANNs in memristance-current plane [56], [57]. Both real and imaginary parts of the admittance of RRAM devices described in [58], [59] exhibit hysteretic behavior with the change of the driving voltage. The coexistence of hysteretic memristive and memcapacitive behavior is analyzed in [60], [61]. Gas discharge lamps are the subclass of memristors having non-crossing PHLs [31] and inverse hysteresis in the flux-charge plane [62].
In this paper, we relate the classical hysteresis theory (Takacs and Duhem models) to the mathematical modeling of memristor. According to our best knowledge, the first usage of the Takacs model of hysteresis in mathematical modeling of memristor is proposed here. We also propose the modified Takacs model, in order to be able to describe asymmetric HL and corresponding first order reversal curves. We also provide two examples of the usage of Takacs model in fitting experimentally obtained data of some of fabricated memristors. Additionally, we show that the differential equation corresponding to the Takacs model belongs to a class of differential equations in the Duhem model of hysteresis. In particular, we show that a class of differential equations in the Duhem model can be used as state equations of memristors. Moreover, we show that the state equation of the memristor model with Biolek window function [63] also belongs to that class of Duhem differential equations. Accordingly, hysteresis appears in the state-charge plane for the current-controlled memristor and it is rate-independent. In the special case when p = 1 in the Biolek window, we show that the steady state solution of the state equation can be expressed in terms of the Takacs model.
The rest of the paper is organized as follows: Takacs model of HL is briefly described in Sec. 2. In Sec. 3, classical Takacs model is extended to the case of asymmetric HL with the first order reversal curves. In Sec. 4, it is shown that a class of differential equations of Duhem model of hysteresis could be used as the state equation of memristor. It is also shown that the state equation of the memristor model with the Biolek window function, as well as the differential equation corresponding to the Takacs model, both belong to the considered class of differential equations of Duhem model. Three examples, including spintronic memristor [53], ferroelectric memristor [54] and memristor model with the Biolek window function [63] are presented in Sec. 5. Conclusions are provided in Sec. 6.

Takacs Model of Hysteresis
This section provides a brief description of the Takacs model (also called T(x) model) of HL [11], enabling us to model major and minor HLs, inverse hysteresis, first-order and higher-orders reversal curves, demagnetization spiral, etc.
Takacs model is based on the T(x) function, which is a linear combination of a hyperbolic tangent and a linear function: In order to describe branches of the symmetric HL, the hyperbolic part of T(x) is shifted in horizontal direction (to the right and left by a 0 ) and in vertical direction (up and down by b 1 ), see Fig. 1. The ascending branch f + for the increasing values of x is described by whereas the descending branch ffor the decreasing values of x is described by The branches have two common points at the tips of HL (see Fig. 1). Since hyperbolic tangent is an odd function, we can assume, without the loss of generality, that For periodic driving signal x = x(t), with zero mean and -X m  x  X m , the tips occur at x = X m . Symmetric HL is closed when f + (X m ) = f -(X m ). This relation can be solved for b 1 : In the context of magnetic materials, parameter A 0 in (2) and (3) is related to reversible magnetization, and can be used to "skew" the HL. In what follows, coefficient A 0 is set to zero. For A 0 = 0, the horizontal asymptotes of f + (x) and f -(x) are B 0 + b 1 and B 0 -b 1 , respectively. Notice that HL can be easily shifted up or down by adding the same constant to both f + and f -.
Setting A 0 to zero in (2)-(3) and solving the resultant equations for x provide the branches of the inverse hysteresis: where artanh() denotes the inverse hyperbolic tangent function. Examples of HLs and corresponding inverse HLs are presented in Fig. 2(a) and Fig. 2(b), respectively.
The largest HL that can be achieved in the system or material is by definition the major HL (see e.g. [11]). When the driving signal x(t) is interrupted and reversed, the direction of f(x) is also reversed, see Fig. 3. Corresponding return path is often called the first order reversal curve (FORC) (see e.g. [11]). Furthermore, when the driving signal is stopped and reversed at X r on the ascending branch, -X m < X r < X m , and  returned to -X m (negative saturation), the corresponding return path f r-(x) can be described by [11]   The ascending branch (2) and down-going return path (8) have two common points: one corresponding to x = X r and the other to x = -X m . Therefore, System (9) can be solved for B down and c down : where b 1 is given by (5) and Dually, when the point of reversal is on the descending branch, the up-going return path f r+ (x) can be described by The descending branch (3) and up-going return path (12) have two common points: one corresponding to x = X r and the other to x = X m . By duality, unknown parameters can be obtained from (10)-(11) by replacing B down with B up , c down with c up , X m with -X m and X r with -X r .

Asymmetric HL
To model asymmetric major HL and corresponding FORCs we provide here a modification of the Takacs model. This type of asymmetric HL and corresponding FORCs are reported in [54] for ferroelectric memristors and in [60] for multilayered metal-oxide structures.
Asymmetric major HL can be modeled via a modified version of the Takacs model by using two sets of parameters: {C 0+ , a 0+ , B 0+ , b 1 as } parameters for the ascending branch and {C 0-, a 0-, B 0-, b 1 as } parameters for the descending branch: For the periodic driving signal x = x(t), with zero dc component and -X m  x  X m , the tips of HL occur for x =  X m , implying that as as as System of equations (15) can be solved for B 0-and b 1 as as 0 When the driving signal is stopped and reversed at X r on the ascending branch of asymmetric HL, -X m < X r < X m , and returned to -X m , the corresponding return path can be described by Since the ascending branch (13) and return path (18) have two common points at x = -X m and x = X r , it follows that as as as as Substitution of (13) and (18) where b 1 as is given by (16),  and  by (17), and An example of asymmetric major HL along with FORCs is presented in Fig. 4.
Dually, when the point of reversal is on the descending branch of asymmetric HL, the up-going return path f r+ (x) can be described by   as as as Notice that the descending branch (14) and return path (22)

State Equations and Hysteresis
In this section, we show that a class of differential equations of the Duhem model of hysteresis could play a role of state equations in the mathematical models of memristors. Since the Duhem model describes rate independent hysteresis, it follows that memristors with state equations of Duhem type exhibit rate-independent hysteresis in the state-charge (state-flux) plane for current-controlled (voltage-controlled) memristor. We also show that the state equation for memristor model with the Biolek window function belongs to the considered class of differential equations of Duhem model. Additionally, we demonstrate that the differential equation that corresponds to the Takacs model also belongs to the same class.
Let us consider a class of differential equations of the Duhem model of hysteresis in x-U plane. The corresponding equation is of the form where x stands for a variable, x(0) = x 0 is the initial condition, and u is the first derivative of U with respect to time, i.e.
max/min denotes the maximum/minimum function, and g(x) and h(x) are continuous functions. We consider a class of Duhem equations (23) in which g and h are functions of the variable x only, and hence independent of U. Equations (23), (24) represent the Duhem model of hysteresis, which is rate-independent, e.g. [5].
Both formation and disruption of filamentary like conductive channels across the insulating film is described by the differential equation of type (23) in [49]. In the model described by Equation (12) of [49], U corresponds to voltage and u to the first derivative of voltage, respectively.
It is easy to show that (23) can be rewritten as where H() is the Heaviside step function: H(u) = 1 for u  0, and H(u) = 0 for u < 0. According to the Duhem model, hysteresis will appear in the x U  plane, where d .
Additionally, (25) can be rewritten as In the context of memristors, (27) is state equation in terms of x (state variable) and u (the first derivative of U ). In the same context, u is the driving voltage or driving current. According to the Duhem model, the solution of (27)-(28) exhibits the rate-independent hysteresis in the x  charge plane for current-controlled memristors or x  flux plane for voltage-controlled memristors.
Particularly, substitution of where k is constant and p positive integer, along with u i  into (27) In what follows we derive the differential equation that corresponds to the Takacs model of hysteresis. We show that this Takacs differential equation is of type (27) and therefore can also be used in the modeling of memristors. For a particular choice of parameters, the Takacs differential equation coincides with the state equation of the Biolek model for 1, p  as it is presented in Sec. 5.3.
Let us assume that the state variable with hysteresis, shifted in the horizontal direction by d 0 and in the vertical direction by X 0 , is described by Recall that it can be assumed, without loss of generality, that B 0 > 0. Since the co-domain of the tangent hyperbolic function of finite argument is an open interval ( 1, 1),  it follows from (32) that Solving (33) we obtain the following range for the state variable with hysteresis described by (32), providing that Taking into account that the first derivative of hyperbolic tangent function satisfies the following identity the first derivative of (32) with respect to the time can be expressed as This expression can be rewritten in the form (27)- (28) providing that Expression (32) describes both the accommodation, as the transient part of the process, and closed HL as the corresponding limit cycle. Thus, the steady state solution of (37), as a limit cycle, is independent of initial conditions [64].
For the prescribed B 0 , C 0 , b 1 and X 0 , the remaining parameters a 0 and d 0 of closed HL (32) can be determined from the coordinates of the tips (0, x min ) and (U max , x max ), where min(U) = 0 and max(U) = U max . Appendix contains the full account of the arguments for the derivation of the following relations where  = tanh(C 0 U max ). The substitution of (42) into (40), (41), followed by the substitution of resulting relations into (32), provides the closed-form expression for the steady state solution of differential equation corresponding to the Takacs model. This result is also confirmed by the numerical simulations. The solution of the differential equation (37) for the prescribed initial condition, which includes transients and steady state, is presented in Fig. 5(a), while the steady state solution only is presented in Fig. 5(b).

Examples of Modeling HL
In this section, three examples of modeling HL are presented. Section 5.1 is related to spintronic memristor [53], Section 5.2 to ferroelectric memristor [54], and Section 5.3 to the model of memristor with the Biolek window function [63].

Spintronic Memristor
HL for spintronic memristor in the memristance-flux plane is reported in [53]. The branches of HL are described by where R + is the branch corresponding to positive voltage branches of (43) can be rewritten as On the other hand, asymmetric HL (Sec. 3), shifted horizontally by d 0 and vertically by R 0 , can be described by Comparison of (45) and (46) provides the following parameter values: Thus, model (43) can be rewritten in terms of modified Takacs model proposed in Sec. 3. In the case of classical hysteresis, the maximum/minimum input corresponds to the maximum/minimum output. In this example, the maximum resistance corresponds to the minimum flux and vice versa. This can be easily handled with both Takacs model and modified Takacs model by using negative values for parameters C 0+ and C 0-. HL described by (46)-(47) is presented in Fig. 6.

Ferroelectric Memristor
Let us consider voltage-controlled ferroelectric memristor, with two orders of magnitude of OFF-to-ON resistance ratio and tunable intermediate states, as presented in [54]. The result of measurements in memristancevoltage plane is depicted in Fig. 7(a).
In this case study, we model the FORCs related to the multilevel resistance states displayed in Fig. 7(a). The FORCs begin on the ascending branch and end at the lower left corner, where memristance is equal to R ON . We estimate R ON  0.18 M. The minimum voltage V min  -6 V and the maximum voltage V max  4.5 V, which can be read on the horizontal axis in Fig. 7(a), imply that the descending branch of major HL is missing. Thus, we adopt V m = 6 V, and estimate that the set of reversal points on the ascending branch is V r  {4.5, 4.05, 3.7, 3.55, 3.35, 3.2} V. Likewise, we assume that R OFF  68 M.
Here, we use the model of asymmetric HL with corresponding FORCs as described in Sec. 3. Additionally, we shift the model vertically by R 0 . The ascending branch for increasing voltage can be written as   as and the descending branch as   as Additionally, FORCs for decreasing voltage have the form Parameter R 0 is set to (R OFF + R ON )/2, i.e. R 0 = 34.08 M, while B 0+ and B 0-are calculated as (see (16)- (17)) Parameters C 0+ , a 0+ , C 0-, and a 0-are used for fitting the major HL. Their values are: C 0+ = 1.45 V -1 , a 0+ = 4.25 V, C 0-= 3.3 V -1 , and a 0-= 3 V. Parameter b 1 as can be obtained by substituting X m = V m into (16). Furthermore, parameters as down

B
and as down c can be calculated from (20) using that X m = V m and X r = V r . The result of modeling is shown in Fig. 7(b).

Memristor Model with Biolek Window
As presented in Sec. 4, the solution of the state equation of the current-controlled memristor with the Biolek window function exhibits rate-independent hysteresis in the state-charge plane.
In this subsection, we provide parameters for which the differential equation of the Takacs model (37) is reduced to the state equation with the Biolek window for p = 1. Since the steady state solution of (37) is given by (32), it immediately follows that the same set of parameters provides the closed-form expression for the steady state solution of the state equation with the Biolek window for p = 1. Additionally, we demonstrate that, for the steady state solution of state equation with the Biolek window for p  1, the maximum value x max and the minimum value x min of the state variable satisfy the relation x max + x min = 1. This relation can be used to estimate transient time as well as a numerical error in the course of finding the steady state solution.
Since the memristance of the model with the Biolek window is represented as a linear function of state variable, off off on it is easy to validate that the existence of the hysteresis in the state-charge plane implies the existence of hysteresis in the memristance-charge plane. State variable vs. charge and memristance vs. charge for memristor model with the Biolek window for p  {1, 2, 10}, are presented in Fig. 8(a) and Fig. 8(b), respectively. Simulations were performed using the Runge-Kutta-Fehlberg 4(5) method.
Notice that it has been recently proved that the memristor model with the Biolek window is characterized by a single stable fixed point [65], [66], meaning that the model is globally asymptotically stable. Fig. 9 presents HLs in the memristance-charge plane and corresponding PHLs for three different driving oddsymmetric periodic current waveforms. For prescribed Q max = max(Q) and given frequency f, the corresponding amplitudes for sinusoidal, triangular and rectangular currents are equal to Q max f, Q max 4f, and Q max 2f, respectively. HLs in the memristance-charge plane in Fig. 9(a) coincide for all three waveforms, which is in accordance with the fact that the Duhem model describes the rate-independent hysteresis.
By comparison of (37) with (30), (31), we conclude that the substitution of B 0 = 1, C 0 = k, b 1 = -1/2, X 0 = 1/2 and u = i into (37) provides the state equation with the Biolek window for p = 1:  Thus, the Takacs model of hysteresis is closely related to the model of memristor with the Biolek window function for p = 1.
Furthermore, the substitution of B 0 = 1, C 0 = k, b 1 = -1/2, X 0 = 1/2, u = i and U = Q into (32) yields Since the Gauss Hypergeometric function 2 F 1 is expressed via an infinite series, its simulation is usually slow and it is recommendable to use an efficient algorithm for solving differential equations, instead. In such case, relation (55) can be used to verify whether the steady state solution is reached in numerical simulations and/or to estimate the numerical errors of the solutions.

Conclusion
According to our best knowledge, this paper provides the first usage of the Takacs model of hysteresis in the mathematical modeling of memristor. A modified version of the Takacs model tailored to describe asymmetric HL and corresponding FORCs, which are observed in some experiments, is also proposed.
In particular, it is presented that a class of differential equations of the Duhem model of hysteresis coincides with a class of state equations of memristors. In this context, hysteresis appears in the state-charge (state-flux) plane for current-controlled (voltage-controlled) memristors. It is also confirmed that the state equation of the memristor model with the Biolek window function, as well as the differential equation that corresponds to the Takacs model, both belong to the class of differential equations related to the Duhem model. Since the Duhem model describes rateindependent hysteresis, it follows that hysteresis in the appropriate planes of memristors is also rate-independent.
Notice that in the Takacs model, the sigmoidal function other than the tangent hyperbolic function is also possible.
We believe that the idea behind this paper can be extended to the modeling of memcapacitor, meminductor and other mem-devices that exhibit HLs in the appropriate planes.
In order to demonstrate that (39) is the only solution of (72) let us introduce an auxiliary function It is easy to observe that the right hand side of (74) is equal to zero at x = X 0 . Furthermore, the sign of the right hand side is the same as the sign of b 1 for x < X 0 and opposite for X 0 < x. Therefore, for b 1 > 0 (b 1 < 0) W T aux (x) is monotonically increasing (decreasing) for x < X 0 and monotonically decreasing (increasing) for X 0 < x. These facts further imply that W T aux (x) attains the same value for some mutually distinct x 1 and x 2 , providing that x 1 < X 0 and X 0 < x 2 . Furthermore, the identity W T aux (x) = W T aux (2X 0 -x) (see (73)) implies that x 2 = 2X 0 -x 1 . Consequently, (39) is the only solution of (72).