A Consistent Model for Short-Term Instability and Long-Term Retention in Filamentary Oxide-Based Memristive Devices

: Major challenges concerning the reliability of resistive switching random access memories based on the valence change mechanism (VCM) are short-term instability and long-term retention failure of the programmed resistance state, particularly in the high resistive state. On the one hand, read noise limits the reliability of VCMs via comparatively small current jumps especially when looking at the statistics of millions of cells that are needed for industrial applications. Additionally, shaping algorithms aiming for an enlargement of the read window are observed to have no lasting e ﬀ ect. On the other hand, long-term retention failures limiting the lifetime of the programmed resistance states need to be overcome. The physical origin of these phenomena is still under debate and needs to be understood much better. In this work, we present a three-dimensional kinetic Monte Carlo simulation model where we implemented di ﬀ usion-limiting domains to the oxide layer of the VCM cell. We demonstrate that our model can explain both instability and retention failure consistently by the same physical processes. Further, we ﬁ nd that the random di ﬀ usion of oxygen vacancies plays an important role regarding the reliability of VCMs and can explain instability phenomena as the shaping failure as well as the long-term retention failure in our model. Additionally, the results of the simulations are compared with experimental data of read noise and retention investigations on ZrO 2 -based VCM devices.


■ INTRODUCTION
Due to the sharply rising demand on portable electronic devices such as smartphones and being concomitant with that, the need for research on highly scaled non-volatile memories strongly increases. 1−3 Due to their promising operation features like scalability, non-volatility, reliability, complementary metal−oxide−semiconductor (CMOS) compatibility, and short switching times, 4 resistive switching random access memories (ReRAMs) are considered to succeed flash technology in the future. 5 Additionally, ReRAMs are considered promising candidates for several neuromorphic and machine learning applications. 6,7 One extensively studied type of ReRAM is bipolar switching valence change-based memories (VCMs) that typically consist of a transition metal oxide layer like HfO 2 , ZrO 2 , or Ta 2 O 5 sandwiched between two metal electrodes. 8−10 By applying an external voltage, oxygen can be extracted from the initially insulating oxide layer at the chemically active electrode. The so-called oxygen vacancies (V O ·· ) left behind in the oxide build a conductive filament significantly changing the resistance of the device. Depending on the polarity of the externally applied voltage, this filament can be ruptured and rebuilt consecutively switching the VCM device between a high resistive state (HRS) and a low resistive state (LRS). 11 These resistance states are non-volatile, can be read out non-destructively, and are used to store the information as logical "0" and "1".
Regardless of the numerous great properties of ReRAMs mentioned above, there are issues concerning the stability of the programmed states warding them from commercial breakthrough so far. 12−14 One challenge is short-term instability, also called read noise, limiting the read window, which separates the HRS and LRS. Here, program-verify or shaping algorithms were proposed to overcome these limitations 15 but have been reported to not have a lasting effect. 13 Instead, the shaped resistance (or read current) distributions always revert back toward their intrinsic statistics. These intrinsic distributions are observed to be very stable on short time scales. 16 However, on long time scales (>1 month) as required for typical memory applications, significant resistance changes may occur that are referred to as retention failure. 17 Both failure mechanisms are so far not fully understood, and to our knowledge, a physical simulation model has not yet been demonstrated to consistently cover both. Puglisi et al. explained a random telegraph noise via charge trapping at slow defects but did not discuss retention effects. 18 On the other hand, Aldana et al. stated that unbalanced V O ·· generation and recombination are responsible for the LRS retention failure 19 but the read current instability was not discussed. In the work of Huang et al., 20 the state instability was connected to Brownian-like hopping in a kinetic Monte Carlo (KMC) model. Furthermore, they used a second, clearly different analytical simulation model to explain retention phenomena mainly by V O ·· diffusion and recombination.
In this work, we will introduce a physical explanation for the processes taking place in VCM devices that are responsible for the instability and retention failures mentioned before, whereas the switching kinetics have been demonstrated in a previous work. 21 Therefore, we present a KMC-based simulation model based on the work of Abbaspour et al. 21 that we extended by the implementation of domains limiting the diffusion of V O ·· in the oxide layer. This extension is motivated by the work of Schie et al. 22 Using molecular dynamics simulations, they found different diffusion regimes for V O ·· in HfO 2 , suggesting that V O ·· may diffuse easily over short distances within certain domains. Within sub-diffusive regions between these domains, the diffusion is expected to require a higher activation energy. 22 These different diffusion regimes are a possible explanation for the high discrepancy in the time of occurrence of the different failure mechanisms. With the implementation of these diffusion-limiting domains, we can accurately show and explain several phenomena concerning the abovementioned instability and retention of the HRS of bipolar VCM devices. Both failure mechanisms can be explained by the same physical origin based on the random configuration and movement of V O ·· . In addition, we present experimental results obtained for VCM ReRAM devices based on ZrO 2 , which is considered nearly identical to HfO 2 with respect to its physicochemical properties. 23 Using these exemplary devices, we demonstrate typical short-term instability and long-term retention failures and compare our simulation model to the experimental findings.

■ SIMULATION MODEL
The simulation model used in this work is based on the work of Abbaspour et al. 21 and represents a typical VCM ReRAM cell. An oxide layer, e.g., HfO 2 or ZrO 2 , sandwiched between two metal electrodes, is modeled by a 6 × 6 × 6 nm 3 lattice structure with a spacing of 0.5 nm between the individual lattice points. Oxygen vacancies (V O ·· ) that are of major importance for the current transport in VCM memory devices 8 can be placed at the lattice points of the oxide structure. There, they locally reduce the lattice cations and induce defect states that act as electron traps, which are responsible for electron transport through the initially non-conducting oxide. 24, 25 At a given V O ·· configuration in the oxide of the VCM device, an external voltage can be applied to the electronically active top electrode (Pt in our work), while the ohmic bottom electrode (Ti in our work) is ground. Solving the threedimensional Poisson equation leads to the potential landscape in the whole device. Here, the charge of the V O ·· is considered twice positive reduced by the electron occupation probability that is between zero and one. Additionally, the heat-flow equation is solved, leading to the local temperature distribution in the cell. Afterward, the current through the device is calculated using a trap-assisted tunneling (TAT) current solver since the HRS of VCM devices with comparably low V O ·· concentrations can be well described by the multiphonon tunneling of electrons via electron traps. 26 Here, two main processes are considered, which are the electron tunneling between the electrode and trap and the electron hopping between two traps. 27 The electron tunneling rate between the electrode and trap is calculated by where R 0 is the electrode-trap coupling factor. The tunneling probability T R is given by the Wentzel−Kramers−Brillouin (WKB) approximation with x 1 and x 2 being the electron positions in the traps and the electrode, the effective mass for tunneling in HfO 2 m* = 0.1 m 0 , the tunneling barrier Φ B , and the electron energy E e . The Fermi factor F gives the number of occupied states in the electrode in the case of electrode-to-trap tunneling and correspondingly the number of unoccupied states in the electrode in the case of trap-to-electrode tunneling and is calculated via with the Fermi−Dirac density function F(E), E t e and E t f denoting the energies of empty or filled traps, respectively, and ρ(E) being the density of states of electrons in the metal electrode. For the electron hopping between traps, the Abraham−Miller formula is consulted so that the hopping rate from one trap at position i to another trap at position j is calculated with the following equation where v 0 e denotes the electron vibration frequency, d ij is the distance between the two traps, a 0 is the attenuation length of the localized wave function of the defects, and V i and V j are the potentials at the trap sites.
The modeling of the system dynamics is based on a kinetic Monte Carlo (KMC) algorithm, where the probabilities of all possible processes in the device are calculated. Then, randomly but weighted by their probabilities, one of these processes is chosen and fulfilled. In our simulation model, three processes concerning V O ·· are possible. The first and central process in this work is the diffusion of V O ·· from one lattice site to a neighboring one. Due to redox reactions, oxygen exchange can take place at the ohmic (bottom) electrode-oxide interface, leading to the generation or recombination of a V O ·· at the interface. In contrast to other simulation models, 27−29 no generation of V O ·· in the bulk due to anti-Frenkel pair formation is possible as they are found to be not stable but recombine nearly immediately. 30, 31 The transition rates are calculated via with the characteristic vibration frequency v 0 , the energy barriers for the respective process E D/G/R , and the potential difference between the two states ΔΦ, where the local potential, the self-potential, and the image potential due to the Coulomb interaction are considered. After choosing one process, the time of the simulation is updated using with a random number r 1 between 0 and 1 and the sum of all transition rates R sum . All parameters used in our simulation model that are relevant for this work are given in Table 1.
Following the results of Schie et al. 22 who found different diffusion regimes via molecular dynamics simulations for V O ·· in HfO 2 , we introduced domains limiting the V O ·· diffusion in our model aiming for a simplified approach to model the different energy barriers for the V O ·· diffusion. This is consistent with their findings that a simple random distribution of energy barriers is not sufficient to describe the diffusion behavior. Therefore, the oxide layer of the VCM cell is divided into several cubic boxes, as exemplarily sketched in Figure 1. Inside these boxes, the V O ·· can move as described before, whereas for the V O ·· diffusion from box to box, the diffusion energy barrier E D is increased to 1.2 eV. For simplification, all boxes have the same diffusion barriers and the same size that can be varied for different simulations.
In general, all simulations in this work are based on a VCM cell programmed to the HRS. This initial state is modeled by a fixed number of V O ·· being placed randomly in a predefined filament region with a relatively large gap between the filament and active electrode. The filament boundaries and the number of V O ·· are the same for all cells of a simulation frame so that the initial variations from cell to cell only originate from the random placement of the V O ·· in the filament region. The initial filament always contains an integer number of cubic boxes in all three dimensions. Additionally, individual V O ·· can also be randomly placed in the layer of boxes above the filament with respect to the non-abrupt transition from the filament to the gap after the RESET process. It should be mentioned that our initial state represents the HRS at a sufficiently large time after the programming so that possible relaxation effects observed in a previous experimental work 16 can be neglected. A typical HRS distribution can be seen in Figure 1.

■ EXPERIMENTAL SECTION
The experimental results in this work were obtained for VCM cells with a lateral size of 7 μm × 7 μm and a stack of 30 nm Pt/5 nm ZrO 2 /20 nm Ta/30 nm Pt, arranged in a 32 × 1 cell line array structure. Onto the Pt bottom electrode, covering the whole substrate, 30 nm SiO 2 was deposited via electron beam evaporation. Using UV lithography, the line structure was transferred to the SiO 2 layer. Subsequently, 7 μm broad lines were etched down to the Pt layer. It may be noted that the SiO 2 layer does not participate in resistive switching but is used to separate the line arrays. After removal of the photoresist, 5 nm ZrO 2 was deposited via reactive RF sputtering (80 W) and another UV lithography step was applied to structure the top electrode. For the latter, 20 nm Ta was deposited via RF sputtering (12.5 W) and in situ covered by 30 nm Pt (80 W) to prevent oxidation. Further details about the device structure and test setup may be found in ref 32.
The experiments were conducted using a custom array tester by aixACCT Systems, which was connected to the array structure via a wedge card with 32 probes. The tester provides an arbitrary waveform generator, a virtual ground to record the resulting current, and a switchbox that maps the applied signal to the cells of the line array. Using the switchbox, different resistances can be put in series with the VCM cell to limit the switching current. During the initial forming operation, which is required to build the conducting filament, a series resistance of R S = 10 kΩ was used. For SET operations, R S was reduced to 1 kΩ and RESET operations were performed without series resistance. The Pt bottom electrode (active electrode) was connected to ground, and all voltages were applied to the Ta top electrode (ohmic electrode). However, throughout this paper, all voltages are given with reference to 0 V at the Ta top electrode.   To program the devices, read-verify algorithms were used. For the forming operation, a voltage ramp with a length of t RAMP = 20 ms toward a maximum voltage V max was applied. After each attempt, the resulting cell resistance was determined by a read pulse (50 μs at −0.2 V), and in the case of failed electroforming, the operation was repeated with increased V max . Here, V max ranged from −2.0 to −4.0 V. The SET operation was performed analogously with t ramp = 5 ms and V max ranging from −0.6 to −1.4 V. During RESET, the same algorithm was used but with a rectangular pulse (rise/fall time 5 μs) with a length of t pulse = 3 ms and V max ranging from 1.3 to 2.0 V. It may be noted that the typical voltages for a successful operation were −3.2 V for forming, −1 V for SET, and 1.8 V for RESET.

■ RESULTS AND DISCUSSION
In the first part of the results, we focus on instability effects in VCM cells programmed to the HRS, which is mainly observed at read operations at comparably low voltages and room temperature. In all simulations under these conditions, only the diffusion of V O ·· can be observed. Due to the small read voltage and the low temperature, no V O ·· generation or recombination at the ohmic electrode happens at the time scales of the experiment. For the same reason, the diffusion of the V O ·· only takes place inside the introduced cubic boxes and no box-tobox diffusion takes place.
Typically, VCMs programmed to the HRS show an intrinsic log-normal behavior of their current distribution. 12,13 As proposed in previous works, random configurations and the movement of V O ·· seem to be the origin of this statistical behavior. 16,33 Therefore, the model was extended, and the so far static investigations are enlarged to dynamic simulations, allowing us to observe the time evolution of single cells as well as the statistics of multiple cells.
In the first step, the statistical current distributions resulting from the new simulation model are compared with the experimental data. A box size of 1 × 1 × 1 nm 3 has been chosen here and in all following simulations. Fifty V O ·· are placed randomly inside a filament volume of 2 × 2 × 3 nm 3 as well as five additional V O ·· in the layer of boxes above. In all simulations, a read voltage of 0.35 V is applied. The corresponding current distribution shows the log-normal behavior that is typical for VCM cells in the HRS 34 and can be observed in the dark blue graph (Sim I t 1 ) of Figure 2. As in experiments, 16 this simulated distribution is in general very stable over time under the mentioned conditions of room temperature and low read voltages. This stability is exemplarily shown in Figure 2 by the nearly similar current distributions at t 1 = 0 s at the beginning and t 2 = 1 s at the end of a read operation. Therefore, it does not make a difference whether a continuous read voltage is applied or single-read pulses with a convenient waiting time in-between are used to determine the current. This stability could not be reached without the introduction of the diffusion-limiting boxes in our model ( Figure S1). Except for the general log-normal behavior, the simulated current distribution differs from the experimental data both in the mean current and especially in the width of the distribution. In real devices, as in our presented experimental data, naturally, there is a device-to-device variability due to the fabrication process and a cycle-to-cycle variability. 35 Adding a variability to the number of distributed V O ·· in the filament (42−52) and the filament size (2−3 × 2−3 × 3−4 nm 3 ) in the simulation model leads to a slightly shifted and much wider current distribution, which looks very similar to the experimental one and can also be found in Figure 2 (red, Sim II). This additional variation of parameters defining the initial state obviously better represents the experimental results. Nevertheless, this variation is left out in the following due to the focus of this work laying on the stochastic processes during read and baking. In addition, excluding these additional parameters leads to a huge reduction of computation time costs.
Although the experimental and simulated current distributions are stable over time, the current of a single cell changes over time. Following our previous work, 16 where we proposed random jumps of V O ·· to be the origin of the observed current jumps, we now investigate the dynamic behavior with the simulation model. In Figure 3a, the current evolution over time is presented for three exemplary experimental (gray) and simulated (colored) current traces each. In general, both, experimental and simulated traces, show discrete jumps of different heights as well as no great changes in the mean current. The traces mainly differ in the high-frequent noise with low amplitude that is observed in the experimental data but not in the simulations. This can be easily explained by the fact that there is no source of electronic noise implemented in our simulation model, which may be the origin of the highfrequent noise in the experimental data, often called the random telegraph noise (RTN) in the literature. 12,36 In an enlarged view of the simulated data, as exemplarily shown in Figure 3b for the red curve (S1), the highest jumps in the current can be identified. In most cases, these large current jumps are caused by vertical jumps of V O ·· at or close to the filament-gap interface. The next smaller current changes correspond to vertical V O ·· jumps farther from the interface, whereas horizontal V O ·· jumps usually only have a small effect on the read current. Thus, the results of the dynamic simulations with the new box model are in perfect agreement with our previous results, where we found V O ·· jumps from the filament to the gap and vice versa to be the most important to discuss read variability in the HRS of VCM ReRAM devices. 16 One common way to analyze and characterize noise is the power spectral density (PSD) plot. 37 In Figure 4, the PSD plots of the experimental and simulated current traces are shown. The magnitude of their PSD values is very similar, and both show a 1/f 2 dependency, which corresponds to a Lorentzian noise. This accordance is further evidence that our basic assumption of randomly moving V O ·· originating read noise phenomena in VCM ReRAMs is valid. As already mentioned before, electronic noise is not covered by our simulation model. In a previous work, it was found that the current jumps caused by electron trapping and detrapping are much smaller than the effect of V O ·· diffusion. 38 Therefore, the high-frequent jumps with very low amplitude, which can be seen in the lower right of the experimental data of Figure 4, are not covered by the simulations. Since we asserted before that the highest jumps in the current are the most problematic for the reliability of the devices, it is reasonable to focus only on the current jumps caused by V O ·· movement. A major problem of HRS instability in VCM ReRAMs is the so-called shaping failure. 15 With program-verify or shaping algorithms, where cells with high currents are reprogrammed or removed, an enlargement of the read window between the HRS und LRS shall be achieved. Unfortunately, these approaches are reported to have no lasting effect. 13 The reasons for this failure are so far not fully understood but shall also be explained by the random movement of V O ·· in our simulation model. Therefore, the time evolution of 100 cells in the HRS generated, as described for Sim I in Figure 2, is investigated. To enlarge the read window, all cells with a current value above 7.5 μA at the start time (0 ms) are removed. The resulting current distribution is shown in black in Figure 5. Subsequently, the distribution relaxes over time and the log-normal behavior emerges from the shaped distribution so that the shaping effect is gone after about 500 ms. Subsequently, the log-normal distribution is stable again. The relaxation time observed in our simulations furthermore fits nicely to experimental results of comparable cells in ref 16, where the recovery of the log-normal shape took about 750 ms.
To sum up our results concerning the investigations of the instability of VCM ReRAMs in the HRS at room temperature and under read conditions, the basic assumption of randomly moving V O ·· can explain different phenomena. It is responsible for the intrinsic log-normal current distribution, the read noise in individual cells, and the shaping failure. Here, the introduced diffusion-limiting boxes are of major importance. Without the boxes, the two scenarios that are presented in Figure S1 are possible depending on the diffusion energy barrier for V O ·· in the oxide layer. Assuming comparably high diffusion energy barriers, no V O ·· would move at room temperature without or under low applied voltages. Therefore, the read noise as well as the shaping failure could not be explained. Assuming lower diffusion energy barriers as in our simulations and as it is also done in other simulations 39 would not lead to stable distributions without the limiting boxes. V O ·· would randomly diffuse through the whole oxide layer, and no stable filament would be possible. Already after seconds, the current of all cells would decrease and no longer resemble the stability of real devices.
In the second part of the results, we want to focus on the retention failure typically occurring in VCM ReRAM devices at long time scales. With our newly introduced simulation model with diffusion-limiting boxes, we want to give a suggestion for   the processes leading to the change of the initially very stable current distributions at large time scales. Again, we propose the random movement of V O ·· to play a major role here. Originally, retention phenomena can be observed at very large time scales of several years, where the programmed resistance state of the cells changes without any external interference except the always present room temperature. Since it is impractical to perform experimental investigations with a duration of several years, "baking" experiments are used to simulate the thermally excited processes usually taking place over years in much shorter time scales of minutes up to a few hours. In the simulations, a comparable dodge has to be exerted by strongly increasing the temperature. Therefore, the long-time retention phenomena can be investigated with an admissible amount of computational resources. This strong increase in temperature leads to very short time scales on which the retention processes take place in the simulation but is necessary to depict as much V O ·· from box to box as possible.
To perform thermally accelerated retention experiments, the resistance or read current could be monitored by continuously reading at a low but constant voltage, while the sample was kept at elevated temperatures. This approach, however, has two disadvantages: On the one hand, it adds an additional voltage stress, which could affect the result. On the other hand, this limits the number of cells that can be tested within a reasonable measurement time and thus lowers the achievable statistics. As we demonstrated before, it is crucial to evaluate retention characteristics regarding the intrinsic statistics of the devices. 33 Thus, we prepared three samples as explained in the Experimental Section. On each sample, a minimum of 200 cells is electroformed, pre-cycled 20 times, and then programmed into the blue HRS distributions shown in Figure 6 (0 min). Afterward, the samples are stored in an oven at 150, 175, or 200°C. After the bake times stated in the legends of Figure 6, the samples are taken out of the oven and cooled down to room temperature and the state of all prepared cells is determined by a read pulse of 200 μs at 0.2 V.
The resulting read current distributions in Figure 6 are impacted by the bake in two ways: The distributions are (i) shifted toward the lower read current with increasing bake time, resulting in a decrease in the median current μ, and (ii) a significant broadening of the distributions is observed, i.e., increasing standard deviation σ. Both trends are extracted from the recorded data and displayed in Figure 7. Here, the impact of the thermal acceleration becomes clear as both effects are increased with the bake temperature. Meanwhile, a decrease in μ might be beneficial for the application, and as it drives the HRS state further away from the LRS, the latter reduces the read window due to the observed broadening, which we have also demonstrated on industrial VCM ReRAMs. 33 In  conclusion, the experimental retention results identify a broadening of the read current distribution as a major retention challenge. Thus, our model aims toward a deeper understanding of the underlying physical mechanisms, which should be compatible with the already discussed short-term instability.
In general, the simulations concerning the retention phenomena are very similar to the simulations concerning instability before. The size and the barriers of the diffusionlimiting boxes have not changed as well as the initial random placement of V O ·· modeling a VCM ReRAM in the HRS. As mentioned before, the temperature is highly increased to 1000 K to simulate the long-term retention behavior at more manageable time scales. Again, the V O ·· jump inside the diffusion-limiting boxes, as already observed in the instability simulations. Additionally, the strongly increased temperature allows V O ·· to overcome the barriers of the boxes so that also box-to-box jumps are possible. Due to the V O ·· concentration gradient created by the high number of V O ·· in the filament compared to the gap or surrounding region with no or only few V O ·· , the box-to-box jumps in total have a predominant direction that leads to structural changes of the V O ·· configuration in the oxide layer of the VCM cell. Generation and recombination processes of V O ·· at the metal/oxide interfaces, which would also be possible under these conditions, are neglected in our simulations concentrating only on the effect of V O ·· diffusion. We will demonstrate that this process alone is sufficient to model typical retention characteristics. In total, two main processes can be observed. On the one hand, V O ·· can jump toward the active electrode, thereby closing the current limiting gap between the filament and active electrode. This leads to higher currents flowing through the cell during readout. On the other hand, V O ·· can diffuse in the radial direction. Thus, the filament widens, and the read current becomes lower. Both processes take place randomly, have a different strong influence on the read current, and superimpose each other. Therefore, no clear differentiation between the processes can be done. In Figure 8 Again, the statistical behavior of multiple cells gives more information than looking at a single cell. Therefore, the baking simulations at 1000 K were performed for 50 cells each. The results of two exemplary baking simulations differing in their initially programmed HRS are presented in Figures 9 and 10. The cells of the baking simulation presented in Figure 9 have initially been programmed to an HRS with a comparatively wide filament (3 × 3 × 3 nm 3 ) containing 50 randomly placed V O ·· and three V O ·· in the box layer above. Here, the effects of the radial diffusion and the diffusion into the gap region are well  balanced. Thus, the mean current (Figure 9b) of the 50 cells is mainly constant except for the random statistical fluctuations due to the low number of simulated cells compared to the experimental data. On the contrary, the current distribution clearly widens due to both diffusion processes away from the initial filament. This widening can be further confirmed by looking at the increase in the standard deviation of the current distribution in Figure 9c. This increase in the standard deviation fits to the experimental results in Figure 7b as well as to the observations of other groups. 20,40 The cells of the baking simulation presented in Figure 10 have initially been programmed to a narrower filament (2 × 2 × 3 nm 3 ) containing 50 randomly placed V O ·· and five V O ·· in the box layer above. Here, the increase in the standard deviation of the current distribution can be observed as well as in Figure  10c. In contrast to the results of the simulations in Figure 9, the mean current is not constant, but clearly decreasing (Figure 10b), due to the dominant radial diffusion. The combination of both trends, the decreasing mean current and the increasing standard deviation, is in good agreement with the experimental results presented in this work.
In addition to the two presented examples, also, other boundary conditions are conceivable, where the diffusion of V O ·· into the gap is dominating, leading to an increase in the mean current of the devices during baking. These considerations fit to the fact that, in the literature, decreasing as well as increasing median currents have been observed as HRS retention failure for VCM ReRAM devices. 33,40,41 Thinking of entropy, the increase in the standard deviation of the current distributions in all simulations is also reasonable. Additionally, it should be mentioned that, in addition to the diffusion of V O ·· , the generation or recombination of V O ·· at the ohmic electrode could also play a role at large time scales or at the elevated temperatures in baking experiments. However, as shown before, the diffusion of V O ·· in our simulation model alone is sufficient to explain the general behavior (tilt and shift) of the current distribution during baking. The broadening of the read current distributions observed in retention experiments necessarily involves single cells with increasing read current. Existing models that neglect the diffusion of V O ·· could only explain a current increase by the stress-induced generation of oxygen vacancy/interstitial pairs in the gap region. 28,42 However, such pairs have been demonstrated to be unable to form without any externally applied field and even if would recombine immediately. 30,31 During retention, there is also no current flow that could stabilize the formation of anti-Frenkel pairs. 43 In contrast, we present a model that is based on the diffusion of V O ·· , which explains the experimentally observed retention effects as well as the short-term instability consistently.

■ CONCLUSIONS
We have presented a 3D KMC simulation model explaining several instability and retention failure mechanisms in filamentary VCM ReRAMs. Therefore, we introduced diffusion-limiting domains to the oxide layer of the device in this work. The main physical process originating from the presented instability and retention phenomena was found to be the random diffusion of V O ·· . The simulated results have been compared with experimental data of ZrO 2 -based devices.
The log-normal behavior of the HRS read current distribution can be explained by the randomness of the V O ·· configuration. Typical read noise behavior with current jumps at short time scales in single devices at the one hand and the overall stable current distribution of multiple devices on the other hand have been depicted by our simulation model. Furthermore, the shaping failure has been explained accurately in the same manner.