Mathematical Modeling of Memristors

2.1. 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₂/Pt structure as shown in Figure 2. Here, the TiO₂ 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₂ 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, 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].

2.2. 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⁻¹⁴ m²s⁻¹ V⁻¹, 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.

2.3. 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.

2.3.1. Joglekar’s window function

Joglekar’s window function can be given as

fx=1−2x−12pE15

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].

2.3.2. 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:

fx=1−x−stp−i2pE16

stpi=1,i≥00,i<0E17

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].

2.3.3. Prodromakis’ window function

Prodromakis’ window function is

fx=j1−x−0.52+0.75pE18

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].

2.3.4. 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:

fx=j1−0.25x−stp−i2+0.75pE19

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.

2.3.5. 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⁻¹⁴ m²s⁻¹ V⁻¹, 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Ω.

2.4. 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].

2.5. 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].

2.6. 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

If we assume that the memristance changes exponentially with w, the relationship between current and voltage can be written as

where λ is fitting parameter [21]: