Nonlinear effects in memristors with mobile vacancies

Because local concentration of vacancies in any material is bounded, their motion must be accompanied by nonlinear effects. Here, we look for such effects in a simple model for electric field-driven vacancy motion in a memristor, solving the corresponding nonlinear Burgers’ equation with impermeable nonlinear boundary conditions exactly. We find non-monotonous relaxation of the resistance while switching between the stable (‘on’ and ‘off’) states, and qualitatively different dependencies of switching time (under applied current) and relaxation time (under no current) on the memristor length. Our solution can serve as a useful benchmark for simulations of more complex memristor models.

Memristors were proposed by Leon Chua [1] as another (originally missing) building block for electric circuits. Their main feature is hysteresis in the current-voltage characteristic or the possibility of having (and switching between) different resistive states. Nowadays, memristors became an important part of developing information storage techniques, such as ReRAM [2]; in-memory computing [3]; neuromorphic computations [4,5]; and other applications [6].
Among different types of memristors, the ones in which oxygen vacancy movement processes play the key role in formation of resistive states have attracted substantial attention [7][8][9][10]. Such states are formed due to much higher mobility of vacancies than that of the metal cations in many transition metal oxides [11]. A model for them can be formulated in terms of the position of the interface between the vacancy-rich and the vacancy-depleted regions [2] moved by the electric field. It can be further generalized by describing its input-output relationship using Bernoulli ordinary differential equation [12]. A more detailed description in terms of spatial distribution of the vacancy concentration was developed in many recent modeling works [13][14][15][16][17][18][19], where the underlying kinetic equations were solved numerically. Inherent non-linearity in these models can be expressed via the nonlinear diffusion equation, predicting a prominent nonlinear effect -formation of vacancy concentration shock waves [20]. The purpose of this paper is to look for other essentially nonlinear effects (absent in the linear approximation) in the nonlinear vacancy-diffusion induced switching behavior of memristors. Here we study a simpler (but exactly solvable) memristor model in strongly nonlinear regime, reducing to Burgers' equation for vacancy concentration with impermeable nonlinear boundary conditions, which is formulated next.
Consider a thin film of a material (with charged mobile vacancies) sandwiched between two metals. The co-ordinate x is counted in the direction perpendicular to the film, which has thickness d. Assume that memristive interfaces at x = 0 and x = d are impermeable for vacancies, so that their total number in the film is conserved. The state of such a memristor at a time t is described by the instantaneous local vacancy concentration C(x, t). Because the number of vacancies is also locally conserved, C(x, t) obeys the continuity equation where ∂ t is the time derivative, and ∇·J(t, x) = ∂ x J(t, x) in the considered one-dimensional context when J = {J, 0, 0}. Suppose that each of the vacancies lives in a periodic potential with the distance a between its minima, separated by energy barriers of the height U A . Interaction of the vacancy electric charge q with the local electric field E makes the potential skewed, setting a preferred direction for vacancy jumps. The probability to overcome the energy barrier and move forward → or backward ← into the neighbouring energy minimum [21] can be expressed as where ν is the attempt frequency, k B is the Boltzmann constant and T is the absolute temperature. From Ohm's law E = ρ 0 I, where the resistivity ρ 0 is assumed here to be a constant and I is the electric current density. The vacancies can only make a jump if 1) they are present at the original energy minimum and 2) there is free space for them (e.g. a movable oxygen atom in the case of oxygen vacancies) at the location of neighboring energy minimum. The corresponding joint probability is c(1−c), where the normalized mobile vacancy concentration c = (C − C min )/(C max − C min ) is defined assuming that C min ≤ C ≤ C max , so that 0 ≤ c ≤ 1. The value of C min is the concentration of immobile vacancies and C max is the maximum concentration of vacancies, determined by the chemical composition of the film's material.
Summarizing, we can express the vacancy drift current due to the electric current I as where D = a 2 ν/2 exp (−U A /(k B T )) is the diffusivity. Another contribution to the vacancy current is due to diffusion and can be described by the Fick's law J diff = −D∂ x c. Substituting the total current J = J drift + J diff into (1), renormalizing the coordinate ξ = x/d, time τ = tD/d 2 and the vacancy current j = Jd/D we arrive at the nonlinear Burgers' equation for the dimensionless vacancy concentration c(τ, ξ) [20]: where p = 2(d/a) sinh (aqρ 0 I/(k B T )) = const. It reduces to the canonical Burgers' equation for the function c = 1 − 2c, but we will solve it directly for c in the present form. It is interesting that the external force (electric current) enters the equation (4) as a coefficient before the nonlinear term, but not as a separate term in the right hand side.
To solve the equation (4) first introduce the antiderivative function noting that u(τ, 1) = r is the total number of vacancies per unit of the film area (filling ratio). Integrating (4) over ξ produces the equation for u: After the substitution u = (1/p) log h (known as the Hopf-Cole substitution [22,23] or the Molenbroek-Chaplygin hodograph method [24]) nonlinear terms in the equation for h(τ, ξ) are canceled and we recover the linear diffusion equation: Its solution can be represented as a sum of a particular solution (satisfying inhomogeneous boundary conditions) and the general solution for the homogeneous boundary conditions of the corresponding type. Because the number of vacancies in the film is conserved, u(τ, 0) = 0 and u(τ, 1) = r, the boundary conditions are of Dirichlet type: h(τ, 0) = 1 and h(τ, 1) = e pr . Guessing the particular solution P and finding the general solution for the problem with homogeneous boundary conditions by separation of variables we get h n e −τ n 2 π 2 sin nπξ p , where the semicolon in function arguments separates (sometimes omitted) constant parameters p and r. One can verify directly that for any set of the Fourier coefficients h n the corresponding c = ∂ ξ u satisfies the original Burgers' equation (4) and has exactly zero vacancy current j = p c (1 − c) − ∂ ξ c at the boundaries: j(τ, 0) = j(τ, 1) = 0 at all times. This exact solution depends on the total number of vacancies in the film r.
The coefficients h n and the value of r can be uniquely determined from the initial conditions. Given c(0, ξ) = c 0 (ξ) we can compute u(0, ξ) = u 0 (ξ) = ξ 0 c 0 (ξ) dξ and using orthonormality of sin nπξ obtain h n = 2 1 0 e pu0(ξ) − P (ξ; p, r) p e −pξ/2 sin nπξ dξ (7) with r = u 0 (1). It is worth noting that under such definition of Fourier coefficients h n they attain finite limit at p → 0 and consequently P → 1. In this limit the derivative c = ∂ ξ u of (6) becomes the cosine series solution of the Cauchy problem for the linear (with p = 0) diffusion equation (4) with homogeneous linear Neumann-type boundary conditions (no vacancy current j = j diff = 0 at either boundary), which relaxes into the uniform state lim τ →∞ c(τ, ξ; 0, r) = r.
where ı = √ −1 and Im denotes taking imaginary part of a complex number. Substituting (8) into (6), we obtain evolution of this initial bump plotted in Fig. 1. As one can see, the external current causes the initial bump to relax into a stable configuration at τ → ∞. The analytical expression for this configuration follows from (6): c st (ξ; p, r) = (1/p)∂ ξ log P , which depends on the current p and the total number of vacancies in the initial state r = u 0 (1). This stable profile forms as the result of the competition between the directed jumps due to the applied current, trying to push the vacancies against the boundary, and undirected jumps due to the diffusion, trying to even-out the vacancy distribution. Except for the value of r, all information about the initial conditions is lost in the final state.
For a given magnitude of |p| there are two distinct stable configurations for positive p > 0 and negative p < 0 current directions. We will call them "on" and "off" states respectively: c on = c st (ξ; |p|, r) and c off = c st (ξ; −|p|, r). These states are mirror-symmetric: c off (ξ; |p|, r) = c on (1 − ξ; |p|, r) and P (1 − ξ; p, r) = e pr P (ξ; −p, r). Switching between them is the main operating mode of the memristor.
To consider the switching it is necessary to express the c off state in the Fourier basis, corresponding to the positive direction of the current p > 0. Using (7) we get where g = (e pr − 1)/(e p − e pr ) > 0 and B(z, a, b) = z 0 z a−1 (1 − z) b−1 dz is the incomplete Euler's beta function. Similarly to (8), the second expression for S was obtained by extending the integrand into the complex plane and making the substitution ξ = −(2/p) log η, which turns it into a product of a rational function in η and a power η α that can be rewritten in terms of beta functions. Equations (9) and (6) make it possible to compute the evolution of the memristor state during the switching. But what about its resistance ?
Since it was assumed from the start that the local resistivity ρ 0 is constant (which is the main leading order contribution), the change of the resistance in the present model may only come from the interfaces (10) where κ i stands for the surface resistivity of the interface i and Dirak's delta functions δ(x) are assumed to be left handed, sitting just outside the range x ∈ [0, d].
There are several possible resistivity change mechanisms (disruption of the crystal structure by defects, valence change of the Mn ions, formation of Shottky barrier, etc). Without delving into details, let us assume a phenomenological series expansion where the constants κ i,0 and k i,1 are defined by the material of the contact i, the material of the main body of the memristor (including its immobile vacancies) and the contact type. Integrating (10) from x = 0 − 0 to x = d + 0 (covering the delta functions) we get for the total resistance where R 0 = (dρ 0 + κ 1,0 + κ 2,0 )/A, k 1 = k 1,1 /A, k 2 = k 2,1 /A and A is the area of the contact. The difference between the resistances of the stable states ∆R = R on − R off is then It is only nonzero if the contacts are different (k 1 = k 2 ) and is maximized at r = 1/2 when ∆R| r=1/2 = ∆R max = (k 2 − k 1 ) tanh(p/4). When the memristor is switched, its resistance changes by ∆R from R off to R on or back. Let us now consider kinetics of this process. Fig. 2 shows the time evolution of the normalized resistance during the switching. It is interesting that R(t) in the figure is strictly monotonous only for r = 1/2,  (14). The dashed lines correspond to r = 0.1, 0.3, 0.7, 0.9 for each value of p, they illustrate that the limit is independent on r.
for other values of r it drops below R off or overshoots the value of R on during the switching process. This happens because either (for r < 1/2) the vacancies depart from the ξ = 0 contact and move as a soliton throughout the film before they start accumulating at ξ = 1; or (for r > 1/2) the vacancies start accumulating at ξ = 1 before they had time to depart from ξ = 0. This is a strictly nonlinear effect and for small values of p 1 the time evolution of the resistance is monotonous for all r. For larger p the range of r values around r = 1/2, corresponding to the monotonous evolution, progressively shrinks. For applications a large value of ∆R is desirable, which, as follows from (13), implies both large value of p and the optimal filling r = 1/2. In this case checking the monotonicity of the resistance relaxation under the applied electric current can be a useful tool for optimizing the filling ratio. It can indicate, based on measurements of a single sample, whether its filling ratio is above or below the optimal value. Evolution of the resistance under the applied current is, basically, a relaxation process towards the equilibrium configuration c on (or c off for p < 0). It is typical that such processes approach equilibrium according to the exponential law ∝ e −τ /τR , where τ R is the relaxation time. This is also the case for the present model. Fig. 3 shows the logarithmic derivative of time evolution of the resistance. At large times these curves become horizontal, which means that relaxation is exponential, their limiting value at t → ∞ is equal to −1/τ R . It is also worth noting that memristors with optimal filling r = 1/2 are the first to reach exponential relaxation regime.
It is not difficult to compute the relaxation time analytically from (6), (12) by neglecting all but the first terms in the Fourier series (as they are exponentially small, compared to the first). This gives which is independent on r and the initial distribution of vacancies. This simple formula contains two key characteristics of the memristor: at p > 0 it gives an estimate of the duration of the current pulse, necessary to switch the memristor into "on" (or "off") state; at p = 0 it gives an estimate of the memristor state lifetime τ 0 = 1/π 2 under the influence of thermal fluctuations. Practical considerations may require to introduce factors before these times, e.g. extending the current pulse duration to be several times τ R to ensure that bit is completely written, or consider the retention time to be several times less than τ R to ensure that bit can still be read reliably. Nevertheless τ R is a convenient simple estimate for both these times.
In seconds the relaxation time (14) is equal to where the last approximate equality assumes a k B T /(qρ 0 I). It has an interesting feature, which is entirely due to the nonlinear term in (4): at zero driving current the relaxation time (bit lifetime in this case) scales quadratically with film thickness d, but at large driving current it (bit writing time in this case) saturates for large d. It means that increasing the thickness (memristor length) to πa/ sinh[(qaρ 0 I)/(k B T )] will increase the bit lifetime quadratically, but the writing time will not be significantly increased. This suggests that vacancy-based memristors in strongly non-linear regime (with large p) may have promising applications for long term information storage. Let us also remark on the case when local resistivity ρ 0 = ρ 00 (1 + γc) is weakly γ 1 dependent on the vacancy concentration c. In the first order of perturbation theory over γ this only leads to rescaling of the constant R 0 by 1 + γr in (12). There is no further impact on the evolution of the resistance because the total number of vacancies r = const in a contact with closed boundaries.
The main limitation of the present consideration is that it does not take into account Joule heating of the film due to the applied current and phase changes in the material. From (15) it follows that the temperature increase reduces the value of p and effectively counteracts the effect of writing current. Thus, in practical applications the maximum achievable value of p is limited. This limit can be controlled by selection of the memristor's material.
In conclusion, we have considered and solved exactly a simple analytical vacancy-migration model for memristors. Its kinetics is governed by nonlinear Burger's equation with conserved number of vacancies and no vacancy current at the boundary. There are two substantially nonlinear effects: non-monotonous relaxation of resistance under the applied current (which can be used for optimizing the initial number of vacancies in the memristor) and the saturation of the memristor switching time for increasing film thickness (which decouples the bit writing time from the bit lifetime). We hope that both these effects can be useful in development of memristor-based memory applications and that the analytical solution (6) can become a benchmark for more complex memristor simulations.