Mathematical Modeling of Memristors

The memristor has quite a reputation as a missing circuit element. It is a powerful candidate for next-generation applications after being first implemented in HP’s laboratories. At this point, mathematical models were needed for the analysis of the memristor, and a lot of studies were done on this subject. In this chapter, mathematical modeling and simulations of the memristor device have been emphasized. Firstly, linear drift and nonlinear drift models have been described on the basic HP model. The window functions used in the nonlinear drift model have been widely examined. Different from HP model, the Simmons tunnel barrier and the threshold adaptive memristor model (TEAM) have been also mentioned. As a result, the most widely used modeling techniques have been described in detail.


Introduction
In the circuit theory, it refers to the existence of three basic circuit elements that define connections between basic circuit parameters such as current (i), voltage (v), charge (q), and magnetic flux (φ). These are resistor, inductor, and capacitor. However, a circuit element that determines the relationship between the charge and the magnetic flux is not defined. The fourth fundamental circuit element representing this relation was firstly presented by Chua in mathematical terms in 1971 with the name of the memristor (memory + resistor) [1]. In 2008, a group of researcher from HP laboratories announced that they were physically producing memristor [2]. In Figure 1, the relationship between fundamental circuit elements and basic circuit parameters is shown.
The relationship between current, voltage, charge, and flux for the memristor is given by: where M(q) has the unity of resistance and W(φ) has the unity of conductance [1].
Memristor has properties that are different from other fundamental circuit elements and can only be seen in a memristor such as nonvolatile memory effect, passivity, and pinched hysteresis loop.
When Eqs. (1) and (2) are opened, Eqs. (5) and (6) can be written as follows: Eqs. (5) and (6) show that the memristance value is related to the history of the current passing through the memristor. That is, when the current passing through the memristor is cut off, it remains at the value of the memristance value. The memristance value starts to change from the last value when it provides the current again to memristor. In other words, the memristor has a nonvolatile memory effect. On the other hand, memristor is not an element that stores energy [3,4].
Memristor is similar to resistor with memory. It shows a nonlinear resistance characteristic that the charge parameter is state variable [5].
Another distinguishing feature of the memristor is that the I-V change shows the pinched hysteresis loop characteristic. A memristor fed by a bipolar periodic signal always exhibits a pinched hysteresis I-V characteristic that passes through the origin. As the frequency of the excitation signal increases, the hysteresis lobe area decreases monotonically. When the frequency tends to infinity, the pinched hysteresis loop shrinks toward a single-valued function [6,7].
The memristor, with being a passive circuit element, has a unique ability to remember the state of resistance that it possesses by maintaining the relationship between voltage and current time integrals. Due to these features, they are being nominated for many different applications such as resistive memories, soft computing, neurocomputing, etc.

HP memristor model
Memristor-based applications require a suitable model for analysis and simulation of the system. When looking at the literature, the HP memristor model where the memristor mechanism based on the drift of oxygen vacancies is widely used. The memristor model developed by HP Lab is composed of Pt/TiO 2 /Pt structure as shown in Figure 2. Here, the TiO 2 layer in which the one side doped with positive charged-rich oxygen vacancies (TiO 2Àx ) is placed between two platinum layers [2].
The doped part of the TiO 2 layer exhibits a low resistance behavior, while the undoped part exhibits a high resistance behavior. As a result of the appropriate excitation on this structure, Figure 2. Structure of memristor reported by HP Lab [2].
Mathematical Modeling of Memristors http://dx.doi.org/10.5772/intechopen.73921 the ionic drift between the doped part and the undoped part results in a dynamic change in the width of the doped region. That is, the width of the doped region is taken as a state variable. As the width of the doped region approaches zero (w!0), the memristor goes to a high resistance state (HRS), and as the width of the doped region approaches D (w!D), the memristor goes to a low resistance state (LRS) as shown in Figure 3 [3].
Since the memristor's dimensions are very small (in few nm), it causes a change in the doped region even with a small stimulation applied. Thus the resistance of the memristor varies between HRS and LRS [3].

Linear drift model
In the model known as the linear drift model, the relation between the current and the voltage of the memristor is defined by the following equation where R on and R off are the values of the resistance for w(t) = D and w(t) = 0, respectively [2,8,9].
From Eq. (7), the value of memristance can be expressed by As shown in Figure 3, the state of change of the memristor resistance is represented by x(t) value in Eq. (8). The speed of movement of the boundary between the doped layer and undoped layer is expressed as dx/dt with Eq. (10) [12,13]: where μ v is the average drift mobility of the charges. If Eq. (10) is taken integral for time, the following expression is derived: If Eq. (11) is put in its place in Eq. (9), the following memristance expression is achieved: This expression can be written as if the left side of the total expression is neglected because R on < < R off [2].
Through this model, the characteristics of the memristor can be observed by the simulations as shown in the ongoing part. The following values are used for the simulations performed in this section. μ v = 10 À14 m 2 s À1 V À1 , D = 10 nm, initial value of w is 3 nm, input signal V input = V o .sin (ωt) where V o = 1 V and f = 1 Hz (ω = 2πf), R on = 100 Ω, and R off = 160 kΩ. Figure 4 shows the change of the current and voltage of the memristor with time for the given parameter values. The pinched hysteresis loop in the I-V plane shown in Figure 5 is one of the fingerprint characteristics of the memristor. Figure 6 shows the relationship between state variable and memristance. This indicates that the memristance depends on the state variable x. This figure also shows that the state variable is limited between 0 and 1. In Figure 7, the change of memristance with applied voltage is seen.
As shown in Figure 8, as the frequency increases, the I-V pinched hysteresis loops become narrower. As the frequency increases toward infinity, the I-V characteristic seems to be a linear resistance characteristic. Figure 9 shows the variation of the I-V characteristic for different amplitude values of the excitation signal.

Nonlinear drift model and window functions
The linear drift model supposes that the state variable (x) of the memristor is proportional to the charge flowing through the memristor. This proportion is acceptable to the interface between the electrodes and the interface between the doped and undoped parts of the memristor. The position of the doped part changes with the applied input signal. Furthermore, the linear drift model assumes that the vacancies have the freedom to move along the all length of the memristor. These assumptions made in the HP model have been greatly simplified, neglecting some basic laws. The reported literature shows that the drift of vacancies is not linear in the region near the boundary interfaces. The reason is that even a small excitation signal can create a large electric field causing nonlinear drift of the vacancies near the boundary interfaces in the memristor. Another problematic situation related to the linear drift model is that the state variable (x) never reaches zero, indicating that oxygen vacancies are not present in the memristor. Similarly, the doped region cannot cover the entire length of the memristor, because there will be no undoped part and the memristor will not work in this way [3,14,15].  In order to provide nonlinearity for the boundary problems mentioned above, functions called window function are introduced. This function is implemented by rearranging the expression Eq. (10) as The function f(x) should have zero at the limits of the memristor (x = 0 and x = 1) and maximum value at the middle of the memristor (x = 0.5) [11]. An effective window function should satisfy the following conditions for modeling of nonlinearity [16]:  • The function should take account of the boundary situation at the top and bottom electrodes of the memristor.
• The function should provide nonlinear drift across the entire active area of the memristor.
• The function should ensure linkage between the linear and nonlinear drift models.
• The function should be scalable in the interval of f max (x) can be obtained such that 0 ≤ f max (x) ≤ 1.
• The function should include the control parameter to set the model.
Many different window functions are proposed as a result of the studies carried out in order to provide these criteria.

Joglekar's window function
Joglekar's window function can be given as where p is the control parameter which changes the flatness of the f(x) curve around its maximum value at x = 0.5 and is a positive integer [17].
In Figure 10, the change of window function proposed by Joglekar for different p values is shown. The characteristic of this function is similar to the rectangular window function by increasing p value, and the nonlinear drift effect is reduced. The disadvantage of Joglekar's window function is the cling situation of the state variable at the boundaries, and it is difficult to change the window function due to the zero value at both boundaries. That is, the nonlinear drift problem is solved, but the boundary lock is not taken into account. When memristor arrives in R on or R off terminal condition, this state will be maintained forever due to zero value taken from the window function [13,14].

Biolek's window function
Biolek has introduced a window function that provides a solution for model errors (the cling situation of the state variable at the boundaries) of Joglekar's window function. Biolek's window function is expressed as follows: where p is positive integer and i is the memristor current [18]. Figure 11 shows the variation of the window function proposed by Biolek for different p values. The proposed window function by Biolek depends not only on the state variable but also on the current flow through the memristor. Thus, the problem of boundary lock is resolved. However, this window function does not include the scalability factor, so the maximum value of the window function cannot be set to a lower or greater value [13].

Prodromakis' window function
Prodromakis' window function is  where p and j are a positive real number [16]. In Figure 12, the change of Prodromakis' window function for j = 1 and different p values is shown. Figure 13 shows the variation of the Prodromakis' window function for p = 10 and different j values.
Prodromakis proposes a solution for the scalability problem in the aforementioned by the presented window function. Prodromakis' window function provides a connection to the linear dopant drift model for sufficiently large values of p. However, the model built by Prodromakis still contains the problem of boundary lock [13].

Zha's window function
This function has been introduced by Zha as a new window model so that boundary lock, scalability, and nonlinear effects can be met at the same time. Zha's window function is expressed as follows: where stp(i) is given in Eq. (17) and p and j are positive real numbers [13]. Zha's window function for j = 1 and different p values in Figure 14 is shown. Figure 15 shows the Zha's window function for p = 10 and different j values.

Comparison of window functions
In this section, the simulation results have been given over the nonlinear drift model of the window functions given in the previous sections. Figure 16 shows the change of window functions according to state variable x for p = 10 and j = 1 values. In Figure 17, I-V characteristics have been plotted for the window functions in Figure 16. In Figure 18, the change of the memristance of the memristor with the applied voltage for the window functions of Figure 16 has been presented. For simulations using these window functions, μ v = 10 À14 m 2 s À1 V À1 , D = 10 nm, initial value of w is 3.145 nm, input signal V input = V o .sin(ωt) where V o = 1.2 V and f = 1 Hz (ω = 2πf), R on = 100 Ω, and R off = 160 kΩ.

Exponential model
Even in the models described so far, the nonlinearity of the large electric field in the memristor is still not taken into consideration. In [19], an exponential model that accounts the effect has  been presented. In this model, the relation between the current of the memristor and the voltage is defined as follows: where β, α, χ, and γ are experimental fitting parameters. How the state variable can influence the current is determined by the n parameter [12,19]. According to Eq. (20), when the model is ON state, asymmetrical switching behavior is shown (sinh part). When the OFF state, the exponential part of the Eq. (20) has the dominant part of the current, which is similar to an ideal PN junction [12,14].
In this model, the differential equation of state variable is written as where a and m are fitting parameters. f(x) can be any window function [14].

Simmons tunnel barrier model
The models described so far were based on the HP model, which consisted of two regions, each of which was modeled as resistance. But Pickett presented another physical model of the memristor as an alternative to the HP model, consisting of a resistor and an electron tunnel barrier in series [20]. Figure 19 shows the memristor structure of the Simmons tunnel barrier model where w is the tunneling barrier and R s is the channel resistance.
w is the state variable of the model and can be written as where f off , f on , a off , a on , i off , i on , and b are fitting parameters [20]. The f on value has an amplitude order greater than f off . It is also effective in changing w in both parameters. The i on and i off parameters effectively limit the current threshold. In this model, a window function is not required, because a off and a on values force upper and lower bounds of x, respectively. Although this model is the most accurate model for a memristor, it is a nongeneric model that is defined for a particular type of memristor, which has a nonobvious relationship between current and voltage [14].

ThrEshold Adaptive Memristor (TEAM) model
TEAM model is a memristor model with several assumptions for analysis simplification and computational efficiency. These assumptions are as follows: 1. There is no change in the status variable for values below a certain threshold value.

2.
A polynomial relationship is established between the current of memristor and the internal state drift derivative instead of the exponential dependence.
Taking these assumptions into account, the derivation of the state variable is written as where k off (k off > 0,) k on (k on < 0),α off , and α on are constants, i off and i on are current thresholds, and w is the effective electric tunnel width. Dependency on state variable w by f off (w) and f on (w) functions is provided. These functions can be thought of as window functions to limit the state variable between w on and w off . If we assume that the memristance changes linearly with w as in Eq. (7), the relationship between current and voltage can be written as v t ð Þ ¼ R on þ R off À R on w off À w on w À w on ð Þ If we assume that the memristance changes exponentially with w, the relationship between current and voltage can be written as v t ð Þ ¼ R on e λ w off Àwon wÀwon ð Þ where λ is fitting parameter [21]:

Conclusion
This chapter describes the mathematical modeling and simulation of the memristor device.
Firstly, brief information about the historical development of the memristor has been given. The emergence of the memristor idea and the formation of mathematical theory have been mentioned. Then, information about the realization of the memristor as a physical element and the HP memristor model has been given.
In the memristor applications, the memristor device must be mathematically modeled correctly for analysis and simulation studies. For this reason, mathematical modeling and modeling methods of the memristor have been emphasized.