Fractional Marcus–Hush–Chidsey–Yakopcic current–voltage model for redox-based resistive memory devices

Abstract

We propose a circuit-level model combining the Marcus–Hush–Chidsey electron current equation and the Yakopcic equation for the state variable for describing resistive switching memory devices of the structure metal–ionic conductor–metal. We extend the dynamics of the state variable originally described by a first-order time derivative by introducing a fractional derivative with an arbitrary order between zero and one. We show that the extended model fits with great fidelity the current–voltage characteristic data obtained on a Si electrochemical metallization memory device with Ag–Cu alloy.

---

1 Introduction

Substantial research efforts have been dedicated to the development of electrically-controlled resistive switching in metal–insulator–metal (MIM) devices or memristors, going from new materials discovery^1–7 to modelling and simulation,^8–10 and design and applications.^3,11–13 With both memory and logic capabilities combined at the hardware level, in addition to long retention times¹¹ and high switching rates¹⁴ at relatively low energy consumption,^1,15 these devices are favorably seen as the next-generation building blocks for nonvolatile memories and neuromorphic computing applications.^11,12 In a typical memristor, the resistive switching is based on the electrically-stimulated change of cell resistance usually driven by internal ion redistribution, which actually depends not only on the applied excitation but also on the past history of the excitation.⁶ Physical mechanisms associated with these reversible transitions have been attributed to different effects including valence-change,¹⁶ electrochemical metallization,¹⁷ and phase change effects.¹⁸ They can be either abrupt (binary) or gradual (analogue), and evolve at different timescales, leading to rich and complex device behaviors in this seemingly simple device structure of just three layers.¹⁹ Furthermore, with the wide range of diversity in memristors materials and their morphologies, operating mechanisms, and manufacturing technologies there is an urgent need for the development of a general model capable of capturing accurately and effectively their complex nonlinear dynamics. This is crucial not only for the characterization and comparison between different memristor devices, but also for the investigation of larger scale memristor-based circuits and hybrid hardware architectures, and also to explore similar behaviors observed for instance in biological synapse systems.²⁰

While models at different size scales and thus with different degrees of physical details and computational complexity have been developed for memristors, including but not limited to ab initio,²¹ kinetic Monte Carlo, and finite element method models,²² in this work we focus on the circuit-level (compact) current–voltage behavior of memristors. From this point of view, memristors are generally described by the system of coupled equations:²³

i = G(x,v)·v, (1)

ẋ = F(x,v), (2)

where i = i(t) is the current through the device, v = v(t) is the applied voltage, and x = x(t) corresponds to a state variable or a group of state variables that quantify the internal dynamics of the device. These are, for example, width of doping region, concentration of vacancies in the gap region, and tunneling barrier width.⁸ State variables can not be observed from external electrical behavior.²⁴Eqn (1) follows the i–v curve of the resistive device in question with G(x,v) being the generalized conductance, whereas eqn (2) describes the dynamics of its internal state x based on its prehistory.²⁵ The actual state of a memristor can only be determined by solving eqn (1) and (2) self-consistently. Memristive systems as featured in terms of eqn (1) and (2) are known to possess a pinched hysteresis loop at the origin in the i–v plane in the response to any periodic voltage source.²⁶

Being versatile and modular enough it is the Yakopcic model^27–29 which is most often used to simulate the nonlinear i–v characteristic of wide range of memristors in response to sinusoidal and repetitive sweeping inputs. The model takes into account electron transmission effects, voltage threshold for state variable motion, and nonlinear velocity function for oxygen vacancies or dopant drift, considered to be the most relevant internal state information.²⁹ It follows on the steps of Strukov et al. work,³⁰ and describes the memristor as two resistors in series, one is undoped with high resistance and the other is doped with low resistance, characterized by electron transmission equations so that:²⁹

i(t) = h₁(v)x + h₂(v)(1 − x). (3)

Here, h₁ is used to model the behavior in the low-resistance state of the device, and h₂ captures its behavior in the high-resistance state. The two electron transmission equations are weighted and mixed by the state variable x which is set to take values between zero and one.²⁵ In memristive devices, it is the rate of change of the state variable x that is explicitly determined (2), and is given in the Yakopcic memristor model by the product of the two composite functions g(v) and f(x) such that:²⁹

ẋ(t) = g(v)f(x). (4)

An exponential dependency of the state change to the positive and negative regions of the input voltage v is modelled in terms ofincluding programming voltage thresholds u_p and u_n. The magnitude of state change for a voltage potential is defined with a_p and a_n. The second function f(x) is determined byfor v > 0, while for v < 0, it is defined asEffectively, this function introduces the nonlinear ion motion, as it becomes harder to change the state of the devices when the state variable approaches the boundaries. The degree of this nonlinearity is adjusted by x_p and x_n since the electrode metal used on either side of the dielectric film can react to the dopants differently. In eqn (6), w_p(x,x_p) is a windowing function that ensures f(x) equals zero when x(t) = 1, and in (7), w_n(x,x_n) keeps x(t) from becoming less than 0 when the current flow is reversed. These two functions can explicitly be written as w_p(x,x_p) = 1 + (x_p − x)/(1 − x_p) and w_n(x,x_n) = x/(1 − x_n).

Clearly, in (3), the functions h₁ and h₂ depend on the structure and type of the memristor under study. Several types of resistive switching memory devices can be classified as nanoionic-based electrochemical systems, wherein an ion conductor in the form of electron insulator layer is placed between two electrodes.^31,32 For the case of cation-migration-based electrochemical metallization memory cells, Ag or Cu are typically used as active electrodes, Pt or W as counter electrodes, and a variety of oxides or chalcogenides thin films as solid electrolytes. When a positive voltage is applied, the active electrode material is oxidized at the electrode–electrolyte interface leading to the release of metallic ions in the adjacent electrolyte, followed by drift and diffusion of these ions across the electrolyte, and then their deposition in filamentary-like metal structures at the counter electrode surface. Short-circuit occurs when the filament has grown sufficiently far to make an electronic contact with the opposite electrode, which defines the low-resistance state of the cell. When a negative voltage is applied, the cell returns back, in principle reversibly, to the high-resistance state.³¹ Anion-migration-based valence change cells, on the other hand, are formed by placing a metal oxide between for example Pt or TiN electrodes and another oxygen-affine, lower work function electrode. The low-resistance and high-resistance states are defined based on the electrochemical formation of oxygen-deficient, mixed ionic–electronic conducting filaments, and the nanoionic modification of the potential barrier between the tip of the filament and the electrode it faces.³¹ For these types of redox-based resistive memory cells, it is more appropriate to consider electron transfer theory associated with the kinetics of redox reactions to better describe their i–v characteristics. Furthermore, because the formation and rupture of the metallic filaments follow random paths, the possibility of charge trapping from one operation sequence to another, charge leakage, the dynamics of an internal state variable associated with these cells cannot be defined solely based on its immediate past, in other words via integer-order derivative as in (2). Taking into account the integral past is believed to be more representative for a proper mathematical description of the complexity and dissipative nature of these cells.

Motivated by these observations, we herein propose a circuit-level model for redox-based resistive memory devices, where the current eqn (1) is taken from the Marcus–Hush–Chidsey (MHC) theory^33–35 of heterogeneous electron transfer, while the state variable eqn (2) is taken from the Yakopcic generalized memristive model.²⁷ We consider the dynamics of the state variable with respect to time to be of fractional, non-integer, order. Mathematically, this adds an extra degree of freedom to the model that can be generically correlated to the non-perfect reversibility of the device when looking at it from one cycle to another. We fit the extended model to the experimental data obtained on a Si memristor with Ag–Cu alloy, shown in Fig. 1, as reported in ref. 36, and make a direct comparison with the superstatistics approach developed therein. A close inspection of numerical results unambiguously reveals that switching to the fractional derivative allows one to significantly improve the agreement between the theory and experimental data.

Fig. 1 A schematic showing the layer-by-layer structure of the memristor under test.

2 Memristor model

The generalized i–v relationship for the proposed memristor model is specified by eqn (1) with

h_j(v) = γ_j·h(δ_j·v), j = 1, 2, (8)

where δ₁, δ₂, γ₁, γ₂ > 0 are fitting parameters. As a rule, these parameters are material-specific and temperature-dependent, so that δ₁ and δ₂ can be viewed as magnitudes of the current conductivities, while γ₁ and γ₂ control the curvatures in the i–v curve relative to the applied voltage v. The function

h(v) = h₊(v) − h₋(v), (9)

is based on the MHC model for electron transfer described by the Gauss-Fermi integral,³⁵Here, the ± signs refer to the oxidative and reductive transition rate functions, λ is the dimensionless reorganization energy scaled to k_BT, while the integral over the dimensionless variable z accounts for the Fermi statistics of electron energies, distributed around the electrode potential. The prefactor β specifies the electronic coupling strength and the electronic density of states of the electrode. In eqn (10), λ and β are assumed to be fitting parameters, knowing that β is usually expressed as an exponential term itself that depends on the distance between the donor and acceptor of electrons. This, however, does not affect the generality of the proposed model. Finally, v in eqn (1)–(10) is actually the electrochemical overpotential defined as the difference between the equilibrium Nernst-potential of the metal and the actual electrode potential defined by the external power supply. We will consider the equilibrium potential to be negligible, so that the electrochemical potential is equal to the applied voltage on the device.

For the dynamics of the state variable x(t), we introduce a fractional time derivative in (4) as follows,

D^α_tx(t) = g(v)f(x), x(0) = x₀, (11)

where D^α_t is the fractional derivative operator of order α > 0 in the sense of Caputo,where , while in the denominator stands for the Gamma function, and x^(m)(τ) is the m-th order derivative. In our present study, the order of the derivative is assumed to be 0 < α < 1.

Non-integer α < 1 in (12) for the state variable eqn (11) indicates that its dynamics does not evolve without prior knowledge of all past information of its state. Or in other words, the right-hand side of (12) contains information about all the previous states of the physical system, representing the so-called memory trace.³⁷ The memory trace term increases with α decreasing away from one, and when we are close to the actual instant t at which the variable is evaluated. At the limiting case of α = 1, corresponding to first-order integer derivative, the memory trace part vanishes, and does not have any effect on the dynamics of the state variable.

Thus, the parameters of the proposed memristor model can be summarized into a single vector

p = (α, x_p, x_n, a_p, a_n, u_p, u_n, β, λ, γ₁, γ₂, δ₁, δ₂). (13)

3 Methods

3.1 Evaluation of the MHC integral

The Gauss–Fermi integrals like (10) can be evaluated numerically using the Gauss-Hermite quadrature (see, e.g., ref. 38). In practice,where n corresponds to the amount of sample points, q(z) is an arbitrary function, while z_k are the roots of the Chebyshev–Hermite polynomial H_n(z_k) = 0 with k = 1,…, n. For a given n ≥ 2 the Chebyshev–Hermite polynomial H_n(z) can be identified from recurrence relations

H_n+1(z) = 2zH_n(z) − 2nH_n−1(z), (15)

provided H₀(z) = 1 and H₁(z) = 2z. In (14), the coefficients c_k are given byThus, one can rewrite the MHC integrals (10) in the formIn our numerical simulations, the order of quadrature n = 25, which is deemed more than sufficient for our purpose.

3.2 Solution to the fractional differential equation

The nonlinear fractional differential equation (11) is solved numerically using the Adams-type predictor–corrector method.³⁹ For nonlinear fractional differential equations of the form

D^α_tx(t) = F(t,x), x^(k)(0) = x^(k)₀, (18)

where k = 0, 1,…, m − 1 and , this method can be described as follows. The approach is based on the fact that the initial value problem is equivalent to the Volterra integral equationWe assume that the choice of F(t,x) guarantees the existence of a unique solution in a certain interval 0 ≤ t ≤ T. We divide this interval into N equal pieces as specified by a uniform grid at the points t_n = hn, n = 0, 1,…, N and h = T/N. The basic idea is that using pre-calculated approximations x_h(t_j) ≈ x(t_j), j = 0, 1,…, n, we get the next time step approximation x_h(t_n+1) by means of eqn (19).

Replacing the integral on the right-hand side of eqn (19) by the product rectangle rule, we obtain

where F_j = F(t_j,x_h(t_j)) and b_k = (k + 1)^α − k^α, provided that 0 ≤ k ≤ n. The predicted value x^P(t_n+1) is determined by the fractional Adams–Bashforth method,To obtain a formula for the corrector, one uses the product trapezoidal quadrature formula to replace the integral in eqn (19), where nodes t_j are taken with respect to the weight function (t_n+1 − τ)^α−1. Using standard techniques from quadrature theory, we can write the integral on the right-hand side of eqn (19) aswhere a₋₁ = 1 and a_n = n^α+1 − (n − α)(n + 1)^α, while a_k = (k + 2)^α+1 − 2(k + 1)^α+1 + k^α+1 for k = 1,…, n − 1. We thus come to the corrector approximation, which can be thought of as the fractional variant of the one-step Adams–Moulton method,The numerical error of this method is shown to behave aswith p = min{2, 1 + α}. In practice, we first calculate and store the coefficients given by {b_k} and {a_k} of eqn (20) and (22) as arrays. After that, on each time step we calculate the predictor (21) and then use it to calculate the corrector (23). To speed up the calculations, we apply the fast Fourier transform algorithm to compute the convolutions on the right-hand sides of eqn (21) and (23).

3.3 Fitting method

Suppose the i–v curve is yielded by N measurements {(t_k,v_k,i_k)}^N_k=1, where v_k = v(t_k) and i_k = i(t_k). To fit the model as specified by eqn (1) and (11) to this data, we search for the set of fitting parameters (13) using the least squares method. This is done by applying the trust region reflective algorithm⁴⁰ to minimize the cost function,where i^mod(v,t,x,p) is the model current specified by the right-hand side of (1). The parameters are non-negative, and the fractional derivative order is bounded, 0 < α ≤ 1. Since the current (1) depends on the state variable x(t), for each p we self-consistently solve either the ordinary (4) or fractional (11) differential equation with respect to the state variable in 0 ≤ t ≤ T. As long as the evolution of x(t) is described in terms of the ordinary differential equation, we keep the parameter α = 1 excluded from the fitting parameters vector (13). Eqn (4) is solved numerically using the Runge–Kutta–Fehlberg method, while eqn (11) is addressed by a means of the Adams-type predictor–corrector method on condition that x(0) = 0. Once the state variable x(t) is calculated, we interpolate it at time steps t_k and evaluate the current i^mod(v,t,x,p) for specific points {(v_k,t_k,x(t_k))}^N_k=1. The fitted model is then evaluated and compared to the experimental data using the Normalized Root-Mean-Square Error (NRMSE),where i^mod_k = i^mod(v_k,t_k,x(t_k),p) is the evaluated model current.

4 Results and discussion

We fitted the memristor specified by eqn (1) and (11) combining the MHC-based state-controlled current–voltage relationship and the fractional Yakopcic state variable model to the i–v characteristic data of the electrochemical metallization memory device taken from ref. 36. The device is a Si memristor with Ag–Cu alloying conducting channels that was fabricated following the method of Yeon et al.⁴¹ A schematic of the fabricated memristor can be seen in Fig. 1. For the i–v measurements that were carried out on a BioLogic VSP-300 workstation, six successive sinusoidal voltage waveforms were applied across the two terminals of the device such that

v(t) = u₀ sin(2πft), (27)

with u₀ = 6 V and f = 1 Hz in the time course 0 ≤ t ≤ 6 s. The fitting is then performed on this time interval to the whole six cycles of switching.⁴² Note that the fitting method is implemented on C++ using the FFTW library⁴³ for Fast Fourier Transform and least-squares routine from the GSL.⁴⁴ The implemented framework allows a model with a single set of parameters to be fitted to multiple datasets at once, so that the main procedure returns a complete array of all fitting parameters and costs for each performed run. Noteworthy, our findings suggest the obtained value of α is robust to the variance of the input data.

The numerical simulations of the models with the ordinary and fractional derivatives fitted to the experimental data are presented in Fig. 2. Here, we show the numerical solution of the canonical memristor system (1) and (2) for the MHC–Yakopcic model in terms of i(t) and x(t), where the calculated current is compared to the measured one. As mentioned above, we considered the equilibrium Nernst-potential of the electrode to be zero, so that the electrochemical potential in eqn (1) and (2) is equal to the actual applied voltage on the device. The corresponding fitting parameters are provided in Table 1. As one can see, the MHC–Yakopcic model fits very well to the experimental data with NRMSE = 0.776 for the fractional order α = 0.677. Remarkably enough, NRMSE = 0.781 when the state variable evolves according to the ordinary differential equation, i.e., with α = 1. It is worth mentioning that averaging over cycles is typically done in literature to show a characteristic curve of the device under study. Meanwhile, the device history, or memory, effect is blurred when we try to fit averaged data. Instead, in our analysis, we fit the proposed model to the whole six cycles of switching as shown in Fig. 2. As the result, our model with α < 1 shows improvement in terms of the NRMSE against the same model with fixed α = 1. This improvement demonstrates that the fractional dynamics could be valuable enough for describing memristive switching, provided ageing and effects of degradation.

Fig. 2 Current i(t) and state variable x(t) as calculated by the MHC–Yakopcic model fitted to the data of the memristor under study, provided the dynamics of the state variable is described in terms of (a) the ordinary differential eqn (4) and (b) the fractional differential eqn (11). The experimental data for the current is depicted by the red points.

Table 1 Fitting results for the MHC–Yakopcic model with ordinary and fractional differential equation for the state variable

Parameters	Integer order	Fractional order
α	1.000	0.677
x p	0.587	0.577
x n	0.0	0.989
a p	0.068	0.064
a n	0.093	0.006
u p	2.373	4.883
u n	0.000	0.000
β	33.37	1.377
λ	28.27	17.40
γ 1	30.17	1.743
γ 2	3.663	2.567
δ 1	1.072	4.509
δ 2	1.597	2.315
NRMSE	0.781	0.776

---

In comparison, we evaluated the q-deformed memristor model recently reported in ref. 36 using the provided parameters and got NRMSE = 0.827. Note that this value is almost twice larger than the reported one in ref. 36 since the fit was performed to the data averaged over six cycles of switching. The q-deformed model was derived by taking into account gamma-distributed local spatial inhomogeneities in the device structure. This provided a noticeable improvement in the fitting of the i–v response of the same device under study here when compared to the currently used existing model (i.e., the Yakopcic model with MIM and Schottky electron transmission equations).³⁶ In practice, it was shown that by introducing a single additional parameter that is related to the q-deformed exponential, one can achieve better agreement with the experimental data for the device under study. This q-parameter being of pure mathematical nature can be however associated with the fractional order α.⁴⁵ In the meantime, as mentioned above, the kinetics of electrode reactions in redox-based electrochemical metallization memory cells should be rather described by electron transfer theory. That is, the MHC + Yakopcic model provides about 8% improvement in comparison to the q-deformed model (Fig. 3).

Fig. 3 Current–voltage characteristic calculated by the fitted MHC–Yakopcic model with fractional order dynamics in comparison with the experimental data. The data is averaged over six cycles of switching for both experimental and model i–v characteristic. The NRMSE, shown in inset, reaches its minimum at α = 0.677 as marked by the red dot. Remarkably, the NRMSE exceeds 0.827 obtained as a result of fitting the q-deformed model to the same data in the range of α < 0.64 and α > 0.71.

Here, because α ≠ 1 we may speak of an intrinsic memory embedded in our redox-based resistive memory device. Fractional dynamics is in fact very often observed in electrochemical devices and complex systems.^46–52 Interestingly enough, the saturation time of a typical memristor that is needed to bring it from the low resistance state to the high resistance state under the applied voltage is sensitive to the fractional order α.⁵³ Thus, the latter, in principle, can be identified from knowing the saturation time.

5 Conclusions

In this work we proposed a compact and accurate model for describing the electrical behavior of redox-based resistive memory devices in which (i) the state-controlled current–voltage equation is based on the MHC theory for electron transfer, and (ii) the dynamics of the state variable is assumed to follow fractional time derivatives of order α (0 < α < 1), with the latter adding a non-Markovian or memory trace term to the modeled dynamics. For the numerical solution to the MHC integral we used the Gauss–Hermite quadrature method and for the fractional differential equation of the state variable we used an Adams-type predictor–corrector technique. Goodness of fit to the experimental data is evaluated in terms of NRMSE, and indicates advanced capabilities of the proposed model when compared to recently reported ones. The developed numerical routine allows one to uniquely determine the value of α. It should be stressed that the NRMSE is used herein for evaluation to produce a normalized, dataset-independent quality metric for the tested model. This allows for a more standardized way of comparing models between this work and other studies. The obtained results, in connection to the electrochemical nature of the device under test, point out to necessity to take into consideration fractional dynamics, that could be of importance provided ageing and general degradation of the device, when describing the i–v characteristics of redox-based resistive memory cells. We showed that the proposed model with fractional order α < 1 dynamics provides 1% improvement in terms of the NRMSE in comparison to the model with fixed α = 1. It is worth mentioning that the q-number as introduced previously in the spirit of superstatistics approach³⁶ is can be linked to the fractional order α.⁴⁵

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

A. A. P. and D. Y. acknowledge the support from the Russian Ministry of Science and Higher Education Project no. 075-15-2021-607.

References

1. B. J. Choi, A. C. Torrezan, J. P. Strachan, P. G. Kotula, A. J. Lohn, M. J. Marinella, Z. Li, R. S. Williams and J. J. Yang, Adv. Funct. Mater., 2016, 26, 5290–5296.
2. C.-Y. Wang, C. Wang, F. Meng, P. Wang, S. Wang, S.-J. Liang and F. Miao, Adv. Electron. Mater., 2020, 6, 1901107.
3. Y. van De Burgt, A. Melianas, S. T. Keene, G. Malliaras and A. Salleo, Nat. Electron., 2018, 1, 386–397.
4. E. C. Ahn, H.-S. P. Wong and E. Pop, Nat. Rev. Mater., 2018, 3, 1–15.
5. V. K. Sangwan and M. C. Hersam, Nat. Nanotechnol., 2020, 15, 517–528.
6. S. Satapathi, K. Raj and M. A. Afroz, Phys. Rev. Appl., 2022, 18, 017001.
7. S.-G. Xu, P. Zhang and X. Zhang, Phys. Rev. Mater., 2021, 5, 024603.
8. N. V. Agudov, A. V. Safonov, A. V. Krichigin, A. A. Kharcheva, A. A. Dubkov, D. Valenti, D. V. Guseinov, A. I. Belov, A. N. Mikhaylov and A. Carollo, J. Stat. Mech.: Theory Exp., 2020, 2020, 024003.
9. W. Wang, M. Laudato, E. Ambrosi, A. Bricalli, E. Covi, Y.-H. Lin and D. Ielmini, IEEE Trans. Electron Devices, 2019, 66, 3802–3808.
10. K. Zhang, J. Wang, Y. Huang, L.-Q. Chen, P. Ganesh and Y. Cao, npj Comput. Mater., 2020, 6, 1–10.
11. J. J. Yang, D. B. Strukov and D. R. Stewart, Nat. Nanotechnol., 2013, 8, 13–24.
12. S. Kumar, X. Wang, J. P. Strachan, Y. Yang and W. D. Lu, Nat. Rev. Mater., 2022, 1–17.
13. H. Bao, Z. Hua, H. Li, M. Chen and B. Bao, IEEE Trans. Circuits Syst. I: Regul. Pap, 2021, 68, 4534–4544.
14. A. C. Torrezan, J. P. Strachan, G. Medeiros-Ribeiro and R. S. Williams, Nanotechnology, 2011, 22, 485203.
15. J. Zhou, F. Cai, Q. Wang, B. Chen, S. Gaba and W. D. Lu, IEEE Electron Device Lett., 2016, 37, 404–407.
16. H.-S. P. Wong, H.-Y. Lee, S. Yu, Y.-S. Chen, Y. Wu, P.-S. Chen, B. Lee, F. T. Chen and M.-J. Tsai, Proc. IEEE, 2012, 100, 1951–1970.
17. R. Waser, R. Dittmann, G. Staikov and K. Szot, Adv. Mater., 2009, 21, 2632–2663.
18. H.-S. P. Wong, S. Raoux, S. B. Kim, J. Liang, J. P. Reifenberg, B. Rajendran, M. Asheghi and K. E. Goodson, Proc. IEEE, 2010, 98, 2201–2227.
19. M. A. Zidan, J. P. Strachan and W. D. Lu, Nat. Electron., 2018, 1, 22–29.
20. A. Ascoli, F. Corinto, V. Senger and R. Tetzlaff, IEEE Circuits Syst. Mag., 2013, 13, 89–105.
21. B. Traoré, P. Blaise, E. Vianello, L. Perniola, B. De Salvo and Y. Nishi, IEEE Trans. Electron Devices, 2015, 63, 360–368.
22. S. Larentis, F. Nardi, S. Balatti, D. C. Gilmer and D. Ielmini, IEEE Trans. Electron Devices, 2012, 59, 2468–2475.
23. L. Chua, IEEE Trans. Circuit Theory, 1971, 18, 507–519.
24. Y. Shang, W. Fei and H. Yu, IEEE Trans. Circuits Syst. I: Regul. Pap., 2012, 59, 1906–1918.
25. T. Chang, S.-H. Jo, K.-H. Kim, P. Sheridan, S. Gaba and W. Lu, Appl. Phys. A: Mater. Sci. Process., 2011, 102, 857–863.
26. L. Chua, in Memristors and Memristive Systems, ed. R. Tetzlaff, Springer, New York, 2014, ch. 2, pp. 17–90.
27. C. Yakopcic, T. M. Taha, G. Subramanyam, R. E. Pino and S. Rogers, IEEE Electron Device Lett., 2011, 32, 1436–1438.
28. C. Yakopcic, T. M. Taha, G. Subramanyam and R. E. Pino, IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst., 2013, 32, 1201–1214.
29. C. Yakopcic, T. M. Taha, D. J. Mountain, T. Salter, M. J. Marinella and M. McLean, IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst., 2019, 39, 1084–1095.
30. D. B. Strukov, G. S. Snider, D. R. Stewart and R. S. Williams, Nature, 2008, 453, 80–83.
31. I. Valov, E. Linn, S. Tappertzhofen, S. Schmelzer, J. van den Hurk, F. Lentz and R. Waser, Nat. Commun., 2013, 4, 1771.
32. I. Valov, R. Waser, J. R. Jameson and M. N. Kozicki, Nanotechnology, 2011, 22, 254003.
33. C. E. D. Chidsey, Science, 1991, 251, 919–922.
34. R. A. Marcus, J. Chem. Soc., Faraday Trans., 1996, 92, 3905–3908.
35. Y. Zeng, R. B. Smith, P. Bai and M. Z. Bazant, J. Electroanal. Chem., 2014, 735, 77–83.
36. R. Konlechner, A. Allagui, V. N. Antonov and D. Yudin, Phys. A, 2023, 614, 128555.
37. W. Teka, T. M. Marinov and F. Santamaria, PLoS Comput. Biol., 2014, 10, e1003526.
38. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, New York, 2017.
39. K. Diethelm, N. J. Ford and A. D. Freed, Nonlinear Dyn., 2002, 29, 3–22.
40. M. A. Branch, T. F. Coleman and Y. Li, SIAM J. Sci. Comput., 1999, 21, 1–23.
41. H. Yeon, P. Lin, C. Choi, S. H. Tan, Y. Park, D. Lee, J. Lee, F. Xu, B. Gao, H. Wu, H. Qian, Y. Nie, S. Kim and J. Kim, Nat. Nanotechnol., 2020, 15, 574–579.
42. D. V. Prodan, Fitting framework for memristor electrical models, 2023, https://github.com/Dmitrii2209-skoltech/iv_fiting_w.
43. M. Frigo and S. Johnson, Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP '98 (Cat. no. 98CH36181), 1998, vol. 3, pp. 1381–1384.
44. GSL – GNU Scientific Library, https://www.gnu.org/software/gsl/.
45. R. Herrmann, Phys. A, 2010, 389, 4613–4622.
46. V. E. Tarasov, Fractional dynamics: applications of fractional calculus to dynamics of particles, fields and media, Springer Science & Business Media, 2011.
47. A. Allagui, H. Benaoum, A. S. Elwakil and M. Alshabi, IEEE Trans. Electron Devices, 2022, 69, 5792–5799.
48. A. Allagui and H. Benaoum, J. Electrochem. Soc., 2022, 169, 040509.
49. E. Hernández-Balaguera, B. Arredondo, G. del Pozo and B. Romero, Commun. Nonlinear Sci. Numer. Simul., 2020, 90, 105371.
50. D. Zhang, A. Allagui, A. S. Elwakil, A. M. Nassef, H. Rezk, J. Cheng and W. C. H. Choy, Org. Electron., 2019, 70, 42–47.
51. Y. F. Luchko, M. Rivero, J. J. Trujillo and M. P. Velasco, Comput. Math. Appl., 2010, 59, 1048–1056.
52. R. Metzler and J. Klafter, Phys. Rep., 2000, 339, 1–77.
53. S. F. Wang and A. Ye, Symmetry, 2020, 12, 437.

---