Stochastic Dynamics of Resistive Switching: Fluctuations Lead to Optimal Particle Number

Resistive switching is one of the foremost candidates for building novel types of non-volatile random access memories. Any practical implementation of such a memory cell calls for a strong miniaturization, at which point fluctuations start playing a role that cannot be neglected. A detailed understanding of switching mechanisms and reliability is essential. For this reason, we formulate a particle model based on the stochastic motion of oxygen vacancies. It allows us to investigate fluctuations in the resistance states of a switch with two active zones. The vacancies' dynamics is governed by a master equation. Upon the application of a voltage pulse, the vacancies travel collectively through the switch. By deriving a generalized Burgers equation we can interpret this collective motion as nonlinear traveling waves, and numerically verify this result. Further, we define binary logical states by means of the underlying vacancy distributions, and establish a framework of writing and reading such memory element with voltage pulses. Considerations about the dis- criminability of these operations under fluctuations together with the markedness of the resistive switching effect itself lead to the conclusion, that an intermediate vacancy number is optimal for performance.

Resistive switching is one of the foremost candidates for building novel types of non-volatile random access memories. Any practical implementation of such a memory cell calls for a strong miniaturization, at which point fluctuations start playing a role that cannot be neglected. A detailed understanding of switching mechanisms and reliability is essential. For this reason, we formulate a particle model based on the stochastic motion of oxygen vacancies. It allows us to investigate fluctuations in the resistance states of a switch with two active zones. The vacancies' dynamics is governed by a master equation. Upon the application of a voltage pulse, the vacancies travel collectively through the switch. By deriving a generalized Burgers equation we can interpret this collective motion as nonlinear traveling waves, and numerically verify this result. Further, we define binary logical states by means of the underlying vacancy distributions, and establish a framework of writing and reading such memory element with voltage pulses. Considerations about the discriminability of these operations under fluctuations together with the markedness of the resistive switching effect itself lead to the conclusion, that an intermediate vacancy number is optimal for performance.

I. INTRODUCTION
Resistive switching (RS) refers to a change in the resistance of a dielectric due to the action of an external electric field or an electric flux through the medium. Thereby, the resistance depends on the history of the field or flux passing through the system, hence it can be considered as a hysteretic effect. RS has been observed in a wide range of transition metal oxides (TMO's), such as manganites M nO(OH) [1,2], perovskites CaT iO 3 [3][4][5] and titanium dioxide T iO 2 [6][7][8][9]. Its basic layout is a two terminal device consisting of a TMO sandwiched between two electrodes, as sketched in Fig. 1. In general the strength of the effect, i.e. the ratio of the high and low resistance, increases with smaller system sizes [7]. different realizations are possible. For example, we have phase change memory, consisting of a chalcogenide glass that switches between an amorphous and a crystalline phase [13? ].
Recently, applications particularly in the semiconductor industry have taken up at a rapid pace. Among the most promising candidates is resistive random access memory (ReRam). Other interesting applications are the integration of logic in memory [10], enabling example concepts of neuromorphic computing [11,12]. ReRam is expected to provide highly scalable, fast, non-volatile and low cost memory [13][14][15]. A single such cell is toggled in between its high-and low-resistive state by application of an external voltage or current. Typical implementations for industrial use aim for high density and stack those elements into a 3d nanocrossbar, layered grids of wires with RS cells in between [15][16][17]. However, along with high integration and miniaturization, challenges of reliability due to sneak paths and fluctuations become ever more significant.
To address the problem of read failures and heat dissipation due to sneak paths, Linn et. al. [17] suggested a complementary resistive switch (CRS). Therein, two RS elements are combined anti-serially to one memory cell. Both its logical states have a high resistance, albeit with differing internal states of the constituent elements. Due to the high resistance, the current and associated energy dissipation through the memory cell and the occurence of sneak paths around it are drastically reduced. As such, the concept has been picked up by various works, see e.g. [6,[18][19][20][21][22]. Conceptually, this setup is similar to having a single element with two active switching zones, one at each electrode-TMO boundary, as it will be applied in this work.
Fluctuations can appear externally and internally. They have been studied for phenomenological memristor models, whose resistance is determined by a single internal scalar variable, which denotes the relative sizes of a high and low and high resistivity area. Internal fluctuations are incorporated by adding white noise to this variable, [23,24], with beneficial effects on the RS-effect, such such as increasing the contrast between the resistive states. External fluctuations were studied in the form of noisy impulses switching the states, and depending on the setup can either have a positive [25,26] or detrimental effect [24]. As yet, no study of fluctuations has been conducted for a particle based model.
Such an approach is essential to address many characteristics of resistive switching. From experimental observations it is known that the functionality of a RS device is determined by the electrode-TMO interface and the distribution of oxygen vacancies [27][28][29]. In this setting, a one dimensional lattice model in which a probability distribution of the vacancies evolves depending on the external voltage and the local resistance of each lattice site is proportional to its density of vacancies has been proposed in [30] for manganites, termed the voltage enhanced oxygen-vacancy migration model (VEOV-model) by the authors. The Schottky barriers are incorporated by enhancing the vacancies effect near the interfaces compared to the bulk, resulting effectively in a two active switching zones. Further investigation of this model showed shock-wave like behavior, [31], made plausible by the formulation of a generalized Burgers equation for the time evolution of the vacancy distribution, which was used to predict the commutation speed of an RS element.
of the vacancy distribution. Hence, our goal is to describe a complementary resistive switch based on this mesoscopic model for discrete particles. We elaborate on the requirements to implement such a device, namely as a large contrast between the resistive states, and the reliability against fluctuations. For very small system, the fluctuations will be driven by the inherent stochasticity in the motion of just a few oxygen vacancies. This lets us determine a lower limit for the possible level of miniaturization and we can predict an optimal system size for the resistive switch. We proceed as follows: In section II we will formulate the VEOV-model for discrete particles, i.e. the oxygen vacancies. It is governed by a multivariate master equation with nonlinear transition rates. This allows us to consider not only the system dynamics depending on the parameters of the external driving, but also to further examine the fluctuations, whose magnitude is determined by the number of vacancies in the system. Subsequently, we investigate the stochastic dynamics of the oxygen vacancies in section III. Hereby, we pay special attention to the hysteretic effects, which we will quantify by the area of the corresponding hysteresis loops, and on the minimum and maximum resistances the system visits within a period of a periodic driving (cf. section III A). Further, we introduce continuous space and derive a nonlinear continuity equation governing the evolution of the oxygen vacancy distribution (cf. section III B). This equation can also be considered a generalized Burgers' equation, and hence the dynamics of the oxygen vacancies are interpreted as nonlinear traveling waves. Also, we successfully numerically integrate said equation to compare it with the results gained in the discrete particle picture, and show how these wave processes affect the electric properties of the device.
In section IV, we define the logical states of the a resistive switch with two active zones in terms of the underlying particle distributions. This allows us to express the actions of the driving voltage in these terms and develop a framework of how to write and read information in such a system (cf. IV A). The resulting read operation is finally considered for its stability with regard to fluctuations (cf. section IV B). Together with the previous results, this leads to an estimate for the optimal particle number.

II. A STOCHASTIC VACANCY HOPPING MODEL FOR BIPOLAR RESISTIVE SWITCHING
We base our considerations on the phenomenological VEOV-model. It focuses on the oxygen vacancy defects together with the influence of the electrode-TMO interface on bipolar switching [30]. Justified by experimental measurements showing filamentary paths of high conductance [32], the memristor is modeled as a one dimensional lattice, i.e. a network of resistors in serial order.
In transition metal oxides the oxygen vacancies affect the resistance in a complex way. On the one hand, they provide dopants and thereby facilitate conductance [33]. On the other hand, they disrupt the TM-O-TM bonds. In order to minimize the energetic cost of the vacancies, their place is filled by free electrons from the conduction band, that are localized now in this place and have an energy level below the conduction band [34]. This resistance effect is more pronounced at the boundary of the TMO, because at a metal-insulator transition the TMO is depleted of freely moving charges, an effect that is also known as Schottky-barrier. Hence taking electrons of the conduction band leads to a higher impact on the local resistance here.
For our investigations, the disrupting influence of the vacancies is considered, as it occurs for example in manganites [1,2]. We assume that a specific lattice site i has N i vacancies, every lattice site can hold a total of N 0 vacancies (with, in general, N i N 0 ), and that the resistance of that each site is proportional to its density of vacancies, The depletion zones at the electrode interface are taken into account by further multiplication with a resistivity factor A i , where A i in the bulk is small against A i near the boundaries of the domain. As a result, the local resistance is given by from which we gain the total resistance of the lattice by summation over all sites, The analytical form of A i is given by its actual course can be seen in Fig. 2.b, where it is indicated by the dash dotted lines, for a discrete lattice of length L, with a length of the depletion zones x 0 , and with 100 lattice sites at the positions i(x) = floor(100 · x/L). Therein, the parameters used are k = 20, x 0 /L = 0.1 and A ↑ /A ↓ = 100. Let a configuration of lattice vacancies (N 1 , N 2 , · · · , N L ) be denoted by N , and by {N } the set of configurations that arises from N through a single particle hopping to a neighboring site. This results in the following master equation An external voltage at the electrodes will induce an electrical current I through the system. It can be used to determine the voltage drop over each lattice site, ∆V i . The oxygen vacancies are quasiparticles holding a charge of q dop , whose diffusive motion is biased by the local voltage drop ∆V i . Agitated, oxygen vacancies now move with a certain probability from one site to another, limited by the number of possible free spaces at their respective target sites. Also, the vacancy requires an activation energy E A to leave its place in the lattice, that needs to be weighted against its thermal energy k B θ, with the temperature θ and the Boltzmann constant k B . Together, we can write down the hopping rate The total rate W i→j of a particle hop from lattice site i to j is proportional to the number of vacancies at the starting place, W i→j = N i w i→j . This basic setup is visualized in Fig. (1). Without external voltage, the oxygen vacancies will be equidistributed. If the external voltage is changed according to some time protocol, directed motion of the vacancies appears. We will investigate the system for a sinusoidal voltage driving V (t) = V 0 sin(2πt/T ) with amplitude V 0 and period T , the latter setting the scale for all system times. The voltage drop over a lattice site is proportional to its resistance, hence Here, the electrical current I(t) is, just as the total resistance R, a nonlocal property that depends on the configuration of the entire lattice. Simulation results for this dynamics are obtained via the Gillespie algorithm [35][36][37][38], which provides a feasible performance since it scales with the transition rates, which themselves are exponentially dependent on the external voltage.
For convenience, we introduce reduced units in which all quantities will be given, Thereby R 0 ≡ 1.76Ω follows by aa resistivity weighted equidistribution of the vacancies, which is explained in Appendix B 1. Unless stated otherwise, we choose for the simulations parameters

III. GENERAL DYNAMICS
The time evolution of a system governed by eqs. (2-6) is depicted in Fig. 2. The electrical current I follows the voltage input, until at a certain threshold value above or below zero the resistance suddenly drops significantly and the current spikes. Hereby, the course of I and R is periodic within each half-cycle of the voltage driving, due to the symmetry of the device.
Looking at the time dependent vacancy occupation distribution Fig. 2.b for the positive halfperiod of the voltage driving, we see that this corresponds to vacancies moving from the left boundary region of A i to the bulk. The form of this motion reminds us to a dispersing wave, as also noted in [31], until it reaches the bulk. From there some vacancies finally reach the right boundary region. In the final panels the hysteresis loops in the I − V plane and a two-legged structure in the R − V plane are depicted, cf. [30]. It reveals that the wave reaching the bulk corresponds to a drop in the device resistance, and a corresponding spike in the current. As soon as the vacancies reach the right side, the resistance increases again.
As the transport from an accumulation on the sides to the bulk differs in its dynamics from the transport from the bulk to the side, we see the emergence of hysteresis loops in the I −V and R−V diagrams. Also we note that the voltage drop ∆V i depends on both the external voltage V and the N i -configuration of the entire lattice, and as a consequence so does w i→j . Hence the distribution of the vacancies plays the role of the inherent memory variable of this resistive switching device, much as the boundary positions w does in the HP-memristor [7]. While a transition of oxygen vacancies between lattice sites with different form factors A i exerts a global influence on all reaction rates, jumps between lattice sites with the same A i only affect local rates.
We further want to investigate the system dependence on the sinusoidal driving, particularly on its period T as well as its amplitude V 0 . To that end we study two different quantities. Firstly the maximum and minima of the resistance during a switching cycle, R min and R max , and their ratio R max /R min . Secondly, the area of the hysteresis loops in the RV-diagram, F hyst .
The results are shown in Fig. 3. Therein we see that for very small T the system has no time to react to the sudden changes of the voltage. The minimum and maximum resistances do not differ and the vacancy distribution approaches a state in which, besides fluctuation, every lattice site contributes equally to the resistance, n i ∝ 1/A i , R → 1. The corresponding hysteresis loops vanish. For longer periods, R max and R min start to spread, and the area of the hysteresis loops increases. This is mainly due to the increment in R max when the vacancies accumulate near the interface; however, when the wave reaches the bulk, R min reaches values considerably lower than in the equilibrium distribution also.
The picture for the amplitudes is somewhat different. Specifically, for very small V 0 , the system has a high resistance. The diffusive motion that is counteracting the external driving is dominant and pushes the vacancies into a uniform distribution. We have n i t = 10 −2 , R → n i t A i 17.5. From a certain threshold above V 0 75, the resistance suddenly drops and the switching dynamics sets in, with corresponding spread in R max and R min and hysteresis loops. Now, the external driving dictates the dynamics. With increasing amplitudes V 0 more vacancies can gather in the high resistance regions near the boundaries, which results in a larger maximal resistance R max . Hence there is a connection between amplitude and period so that the system shows meaningful switching behavior.
For the relevant time scale, we remark that the transition rates are proportional to the activation energy, w i→j ∝ exp(−E A ), cf. eq. (5). The larger E A is, the slower all reactions take place. The other element that sets the time scale of the dynamics, is the external voltage V (t). Since we assume a periodically driven system, the length of the period must increase exponentially to offset the delay due to an increased E A .

A. Effects of Fluctuations
The influence of stochasticity is investigated by changing the number of vacancies per lattice sites, while keeping the average vacancy density n i = N i /N 0 and the other parameters constant. This means, that for a reduced total number of vacancies i N i each individual one bears a larger fraction of the total resistance R and hence causes a larger change in the system properties if it hops in between lattice sites with different A i , while the electrical quantities do not differ in principal, since only the ratio N i /N 0 enters (cf. eqs. (2)(3)(4)(5)(6)). For large vacancy numbers on the other hand, each individual vacancy barely affects the system and its dynamics will approach the deterministic behavior described by the corresponding master's equation eq. (4). In particle numbers per lattice site N i , together with the maxima and minima of the resistance and their respective ratio. We see that for small N i , the spread in between R max and R min monotonically increases, which also implies an increase in the size of the corresponding hysteresis loops. Also, we see that the variation between individual cycles are larger for small average particle numbers, which can be seen in form of the larger statistical deviations in this case. These observations are confirmed in the particle R − V diagram, where we indeed see that both the size of the loops and their fluctuations are bigger for smaller N i , whereas for many particles fluctuations barely play a role and the results approach the values we gain by direct numerical integration of the corresponding master equation (4). Hence we conclude that the resistive switching effect becomes more pronounced with increasing fluctuations. In this sense, the results mirror investigations of the HP-memristor, where it was found that an additional noise increases the memristive response of the system [23].

B. Continuum Limit and Burgers-Like Equation for Wave Transit
In this part, we want to introduce a continuum limit of the master equation in a mean field approximation (MFA). Here, only the outline of the calculations is given. A detailed derivation is shown in the appendix A. For the time derivative of the mean occupation density at lattice site i, cf. eq. (1), we have This expression follows from the master equation (4) and is derived in the appendix. Now we decouple the two point and higher order correlations, resulting in Next, we introduce the lattice spacing between neighboring sites and define x := i . Formally, the continuum limit is taken by letting the lattice spacing become infinitesimal while the number of lattice sites N runs to infinity, in such a way that the product of both remains the constant lattice length lim →0,N →∞ N = L. The hopping rates that formerly depended upon the density of neighboring lattice sites now depend on those at an infinitesimal distance , which will obviously introduce a derivative. The whole approach is quite similar to the derivations of Burgers equations as found for molecular motors in the ASEP model [39]. From here, we will use the rescaled length x = x/L and drop the tilde. Let the averaged density profile be denoted by ρ(x), from the connection we obtain ρ(x, t) = n i (t). Recall that our previous choice of a lattice with 100 sites corresponds to = 0.01. The continuous voltage V (x) is introduced analogously to eq. (12). Entering this into the MFA hopping rate eq. (11) and developing the resulting expression to the second order yields the equation As the electrical current through the system is not position dependent, we express local voltage in the following way (cf. eq. (6)) thereby splitting it in a nonlocal and local component. Bearing in mind the continuity equation Here, the term in braces is the oxygen vacancy flux. Since ρ can be considered a probability density with the normalization 1 0 dxρ(x) = 1, this constitutes nonlinear continuity equation for the probability. By rescaling the time, the prefactor 2 exp(−E A ) can be eliminated. By further, expanding the hyperbolic function to the first order, we are left with This equation can be formally compared to the viscous Burgers equation. With the viscosity denoted by ν, it reads ∂ t u = ν∂ 2 x u − 1/2∂ x u 2 . So eq. (16) is a Burgers-like equation with a current and resistance profile dependent prefactor. Following the wave-like motion of the vacancies, this connection has also been noted in [31].
To validate the derived result, we aim at numerically integrating eq. (16) and compare the thereby obtained results with the original numerical results. The integration is done with Neumann boundary conditions, there is no density flux out of the system. The traveling wave solution is shown in Fig. 6. Similar to the results seen in Fig. 2.b, at the beginning of a period there is a surplus of vacancies stored at the left hand side of the RS-device, near x = 0, which increases significantly with the driving period T . Upon being agitated by a voltage, these vacancies set in motion in form of a wave with decreasing peak height, until it reaches the bulk. For T = 1 the distribution in the bulk is almost flat, the incoming wave forms a hump that moves further to the right with time. Meanwhile, before the second bulk interface around x = 0.7 there is initially a hump that flattens and finally disappears as the distribution shifts further to the right. For T = 10 on the other hand, the initial distribution in the bulk looks like an inclined plane. The incoming wave rushes through the bulk here, reversing the orientation of the inclined plane on its way through it.
As we can see, the traveling wave reproduces the previous results with good principal agreement, as shown exemplary in the R − V for two different driving periods in Fig. 7 superimposed on the corresponding diagrams gained by integration of the particle picture. For all driving periods T , we see some alteration of the two-legged structure.

IV. LOGICAL STATES
As we have seen in the previous sections, the resistance of the system depends on how many vacancies are stored within the bulk versus in one of the two boundary regions. Principally, this configuration is similar to two serial switches with counter oriented polarities, and one active switching zone each. To that end, we think of the lattice as cut in half and each site constituting one resistor that can be either in a high resistive or low resistive state, depending on the distribution of vacancies. This also determines the logical state of the device, as indicated in Fig. 8. A gathering of vacancies near one of the electrodes, i.e. in an area with A i 0, leads to a large resistance, whereas in the bulk vacancies only marginally contribute to the resistance, since here A i 0. Hence, a gathering of vacancies near the first electrode-TMO interface corresponds to a setting R

A. Reading and Writing Operations
Now that we have defined the logical states, we can express the actions of the driving voltage in these terms and develop a framework of how to write and read information in our system.
In Fig. 2 we have seen how the positive half-cycle of the external voltage driving acts upon an initial vacancy accumulation near the left interface. Namely, by pushing it to the right interface. During its course, the resistance drops when the vacancy wave hits the bulk, but returns to its initial value once the voltage pulse ends and enough vacancies have accumulated near the right interface. In terms of the logical states, this corresponds an initial 0-state, onto which the positive voltage pulse acts. Afterwards, the system is in a 1-state. Hence the process is the write 1 operation. It is again displayed in Fig. 9.a. In Fig. 9.b the antagonistic process is shown, where the negative voltage half-cycle acts on a 1-state, resulting in a 0-state. Mirror reversed, the vacancies move from the right interface to the left one. Hence the negative voltage pulse performs the opposing write 0 operation. During its course, the resistance has the same dynamics. We note that the initial states here and for the read operation are prepared by several switching cycles 0 → 1 → 0 → . . . .
Let us turn to the reading operation. The state of the system is determined by sending a reading pulse through it and monitoring the reaction. The reading pulse itself is is the positive voltage half-cycle, albeit with half the amplitude of the writing signal V read = V 0 /2.
If the reading pulse acts on an initial 0-state, the resistance drops to its minimum level, R min . For the oxygen vacancy distributions this we see that the initial vacancy accumulation near the left interface dissolves and moves to the bulk. But unlike for the write case, there is no subsequent buildup of vacancies near the right interface. Hence when the reading pulse has passed, the system remains in the low resistive state R (1) low ⊕ R (2) low . The process is shown in Fig. 9.c. On the other hand, if the system is prepared in the 1-state, the reading pulse leads to a slight increase in the resistance above the value R max that is obtained for the switching cycle 1 → 0 → 1. For the distribution, there are virtually no vacancies near the left interface that can be transported to the bulk, while the accumulation near the right electrode already exists and only slightly increases in size. Hence in this case we have only a small effect on the resistance; cf. Fig.  9.d.
This means that two logical states are distinguished by their differing reactions to the reading pulse. Only when it acts upon the 0-state there is a distinct drop in the resistance. Obviously, the reading pulse has changed the initial distribution, hence a reversely oriented reset pulse need to restore the original configuration in either case. . Discriminability and associated probability of a reading error as a function of the average particle number for various numbers of read-reset-cycles. The inset denotes the optimal reading and reset voltages, i.e., those for which the reading resistance drops to a minimum and the reset resistance is close to the value reached in the writing cycle, respectively. Other parameters as in Fig. 2.

B. Clarity and Stability of the Reading Operation
As noted in section III A, the resistance drop increases with lower particle numbers and so do the statistical fluctuations associated with it. To be able to correctly determine the state of the device, we must assure that the respective resistances after a read pulse and the corresponding reset pulse do not overlap. This property is demanded for a single read-reset cycle as well as for several.
In Fig. 10 various such cycles are shown for two different average particle densities N i . As we can see, for N i = 0.2, the possible outcomes for the resistances fall almost perfectly onto discrete levels, one of which is overlapping. Here, we cannot determine in which state the system actually is. For N i = 178 on the other hand, the ratio of R reset /R read is much smaller, but so are the fluctuations accompanying them. Hence both states can be distinguished easily.
In order to quantify the effect, we employ the measure of the discriminability. Originally defined as the difference between the signals with and without noise [40], we adapt it to our current situation by setting with the standard variations of the respective resistances denoted by σ. As such, it is closely related to the signal to noise ratio, which simply is its square: SN R = d 2 . We approximate that the results fall onto a normal distribution, with peaks at R reset and R read and the standard deviation σ. This is justified since we are interested only in an estimate of the error. In Fig. 11 this measure is depicted as a function of the average particle number N i . With increasing N i , the discriminability of the output improves, since fluctuations largely vanish. Moreover, it drops with the number of read-reset cycles, a spread which is more pronounced for large N i . In view of Fig. 10, this can be attributed to the transient of the average resistance values. The optimal V read/reset (shown in the inset of Fig. 11) was determined after a writing operation, 1 0 → read, and apparently differs slightly from the value that is obtained after a read/reset cycle, read reset → read, thereby increasing the interval of the possible outcomes. To interpret the discriminability measure, we check on how many values will be falsely attributed. The measurements of the read and reset resistances fall onto two normal distributions whose peaks are σ · d apart. All measured values x ≤ R read + σd/2 =: x mid are then assigned to the distribution of R read , all larger values to R reset . Hence the probability of an erroneous attribution is given by those reset values smaller than R read + σd/2 plus the reading values larger than it. Expressed as a formula where erf(x) denotes the error function. The corresponding results are plotted in Fig. 11. We conclude that for an average vacancy occupation per lattice site of N i ≥ 10 we are secure independently of the number of cycles. Hence, the unambiguity of the reading establishes a lower boundary for the optimal particle number, while the clarity of the resistive switching effect is emphasized with lower particle numbers. This favors a intermediate range of the optimal average particle number of about N optimal = i N i = 100 · N i = 1000 particles. We note that the number of vacancies scales with the size of the device, hence the risk of erroneous read-outs limits the possible degree of miniaturization. In a rough estimate, we find that this translates into a length of about 18nm, assuming a cubical element of lanthanum-strontium manganite La 1−x Sr x MnO 3−x/2 with a vacancy density of x = 0.3 [41].
In actual systems, there is intrinsic asymmetry between the upper and lower electrodes, which can be reflected by different heights of the left and right flank of the resistivity profile A i in our system. This will result in two different 'leg' lengths in the resistance over voltage diagrams, but does not change the stochastic fluctuations, i.e. the results for the discriminability apply to this scenario as well.

V. CONCLUSION
Resistive switching is a topical and high current interest field in condensed matter physics and a prototype of a nonlinear stochastic switching event. Its main application, the non volatile resistive random access memory, is a promising candidate for future memory technologies. The scaling to small elements however, leads to an increasing role of fluctuations. It is therefore required, to have an understanding of the switching mechanism and the reliability under the influence of these fluctuations. Here we develop a particle based mesoscopic model based on the distribution of oxygen vacancies. Fluctuations caused by the stochasticity of their motion both enhance the resistive contrast and reduce the reliability of the resistive switch. These counteracting tendencies enable us to predict a limit to miniaturization and an optimal system size.
We have introduced a setup to describe a complementary resistive memory switch based on a discrete particle hopping model. Hereby, the spatial distribution of oxygen vacancies plays the vital part in determining the state of the system, which switches between a high and low resistive state. The vacancies' dynamics is given by a master equation with transition rates depending on the vacancy distribution and an external voltage driving.
The application of voltage pulses led to an accumulation of vacancies near the electrode-TMO interface. An antagonistic pulse then sets these accumulations in motion and they collectively wander through the system, thereby affecting the resistance. By formulating the problem in a particle picture, we gained a tool to vary the accompanying fluctuations: less vacancies, more fluctuations. We looked at the characteristics that define an RS element, such as the spread between the high and low resistance, and discovered that they become more pronounced with increasing fluctuations.
The nature of this collective vacancy motion could be elucidated by deriving a nonlinear continuity equation in a MFA from the master equation in continuum space. It has the structure of a Burgers equation with an additional nonlocal factor, that also incorporates the influence of the driving. We succeeded at numerically integrating this equation to find good agreement with the results gained by particle simulations. Hence we interpreted the motion of the vacancies as a nonlinear traveling wave.
Further, we defined binary logical states in terms of the underlying particle distributions. By linking the actions of voltage pulses to switches between these states, we have established the reading and writing operations necessary to run a memory element. Interestingly, investigations in the stability and discriminability of these operations let us gain a lower limit for the possible number of vacancies in the system, a quantity that is connected with the possible level of miniaturization. Together with the finding of enhanced RS-characteristics for fewer vacancies, this results in an optimal performance for about 1000 oxygen vacancies, which corresponds to a device length of about 18nm.
Since only hopping between neighboring sites is allowed, the number of possible transistions is vastly reduced. By N = (· · · , N i , N i+1 , · · · ) we denote a configuration, that differs from N only at the i-th and i+1-th position, where it takes the values N i and N i+1 respectively. The remaining transitions of the second sum of eq. (A3) can now be written explicitly, yielding We seek to collect the terms in the brackets by performing a substitution N i ±1 = N i , such that the probability function depends on a configuration N that is not altered by subtraction or addition of a particle. This substitution can either affect N k or not. First, the terms of eq. (A4) without particle hop at the k-th position: Since the summation is done over all possible configurations in both sums, we can rename N to N and see that all summands vanish. In the remaining terms the hopping involve the particles at the k-th position, and hence substitutions affect the prefactor N k in eq. (A4), Relabeling the marked quantities N to unmarked ones shows that all expressions with prefactor N k cancel each other out. This leaves us only with The transition rates W can now be connected to the reaction rates eq. (5) by considering the proportionality with the number of particles N i at an individual lattice site, With the further introduction of the density of filled oxygen vacancies n i := N i /N 0 , eq. (A7) becomes In this part, we introduce the continuum limit of the master equation in a mean field approximation (MFA). The MFA now decouples the two point and higher order correlations, resulting in Next, we introduce the lattice spacing between neighboring sites and define x := i . The continuum limit is taken by letting the lattice spacing become infinitesimal while the number of lattice sites N runs to infinity, in such a way that the product of both remains the constant lattice length lim →0,N →∞ N = L. Let the averaged density profile be denoted by ρ(x, t), from the relation we obtain ρ(x = i , t) = n i (t). The hopping rates that formerly depended upon the density of neighboring lattice sites now depend on those at an infinitesimal distance , which will obviously introduce an derivation. The continuum counterpart of the voltage drop eq. (6) further introduces a sign depending on whether we consider forward ( = +| |) or backward ( = −| |) fluxes, for it is a directed quantity. Entering this into the MFA hopping rate eq. (A10) yields The density is expanded up to the second order in ρ according to ρ(x+σ , t) = ρ(x, t)+σ ∂ x ρ(x, t)+ 1 2 2 ∂ 2 x ρ(x, t)+O( 3 ) and the resulting terms are collected with regard to the sign in the exponential, With the definitions of the hyberbolic functions follows In the last recast, we have also disregarded the 2 order in the sinh term. As the electrical current through the system is not position dependent, we express local voltage in the following way (cf. eq. (6)) thereby splitting it in a nonlocal and local component. Bearing in mind the continuity equation ∂ τ ρ + divj = 0, we want to put the spatial derivation in from of the right hand side of eq. (A15).
Since the inner derivatives of the hyperbolic functions introduce a further order in , we disregard these terms and the hyperbolic functions commute with the derivation operator. We gain from which we can easily read the expression for the diffuse and directed oxygen vacancy flux, j = j dir + j diff . Finally, by further expanding the hyberbolic functions to the first order, we obtain the expression Here we specify our assumptions for the form of the resistance profile A(x). While Rozenberg et. al. [30] used a step profile, for simulation purposes it is feasible to smoothen out the regions between the steps. Also, from a physical point of view, a smoother conduit seems reasonable, since both regions are essentially of the same material. The analytical form was given by eq. (3). The higher k is, the more rapid the transition from high to low resistivity occurs. For all simulations, we choose the parameter set k = 20, L = 1 and x 0 /L = 0.1. Obviously, the limit k → ∞ gives a step profile for A(x). For the discretized simulation, the number of lattice sites is set to N LS = 100.
Given a device of length L that has a certain total resistance R, the resistance profile cannot be independent of the number of lattice sites we choose, but scales with it. If for some simulations for example due to numerical expenditure, we need to reduce the number of lattice sites, A i scales accordingly. In fact, we should think of the resistance profile that is introduced in eq. (2) as a resistance per length We then choose A ↓ and A ↑ according to the resistances we want to accomplish, and the number of lattice sites, whereas the length D plays no role. Concerning the dynamics, the relative differences in A i are more important than its absolute value, whose influence can just as well be scaled with other parameters, such as the strength of the driving signal, V 0 . Further calculations are done with A ↑ = 100Ω and A ↓ = 1Ω. The initial distribution of vacancies is chosen in such a way, that all the lattice site contribute equally to the resistance. This choice minimizes the necessary thermalization time of the system. As the individual sites' resistance is proportional to A i (cf. eq. (2)), the individual vacancies are placed at site i according to the probability In effect, all sites then contribute with for the choice of parameters of the simulations, cf. Fig. 2, this yields the normalization R 0 1.76Ω.

Numerical integration of the Burgers Equation
The generalized Burgers equation is a nonlinear flux transport equation, which can be cast into weak form on multiplication by a suitable test function, ψ and integrating by parts The neglected boundary terms ensure that the natural boundary conditions are enforced: no flux out of the domain. We discretise solve equation (B4) using piecewise quadratic finite elements for the unknown density ρ. The nonlocal term I(t) is handled by treating it as an independent additional variable with the defining equation .

(B5)
This equation is assembled in the same integration loop over the elements used to construct the discrete nonlinear residuals from equation (B4). The non-linear residuals are assembled and solved by the C++ library oomph-lib [42]. Time is advanced by an implicit second-order backward difference method (BDF2). The resulting matrix system is solved by a direct solver. Adpativity in space and time is required to achieve accurate results and it has been confirmed that the results do not change if the error tolerances are reduced.