A Physical Description of the Variability in Single‐ReRAM Devices and Hardware‐Based Neuronal Networks

Resistive switching devices are candidates to be used as weight for the future hardware realization of neuronal networks. Herein, these resistive switches base on oxygen vacancy migration and can be toggled between different resistance states. Internal processes lead to a variability over time in these states. To understand this, a dynamical model is developed, which links the internal physical mechanism to the current response. Thereby, the electronic current is calculated by tunneling processes over oxygen vacancies. The model allows to understand the physical transport properties and the current–voltage dependence. Integrating ionic hopping processes clarifies the variability in the electronic conduction. The derived model does not require empirical input parameters and only depends on the properties of the materials. This allows to use literature values and independently compare them to electrical measurements. The model predicts the experiments correctly. Thereupon, the physical properties of the conduction are investigated with respect to the transport over discrete energy levels, which change under the hopping of oxygen vacancies. Afterwards, the current level and variability are simulated for increasing oxygen vacancy concentrations to investigate analog switching, which is required for neuronal network layers. Finally, the variability is tested with multiple interacting cells in a neuronal layer.

certain findings.This generalization requires a precise and dynamical transport description, which is still missing.For this step, the current must be described quantitatively correct.The missing description hinders the design of modern VCM-based neuronal computing networks.For the exploration of the complexity of many interacting cells, a model with high predictive power is required.This model has to cover the current understanding of VCM devices and only rest on material properties to be universally applicable.
Another important issue for the application is the intrinsic variability in transition metal oxide-based ReRAM cells, which leads to current variation over time.First experiments by Schnieders et al show that the variability of the resistive states depends on the type of conduction. [35]For a deeper analysis on the variability, a model is required, which keeps track of the position of oxygen vacancies and connects this to current modification.Ideally, this model is free from empirical parameters and should only depend on the material properties of the oxide and the electrode metal.Under these circumstances, the model enables material engineering to optimize the device performance.In particular, the model needs to predict the variability correctly to be of use in the simulation of neuronal networks, which are physically realized by means of ReRAM weights.The variability of these weights influences the robustness of the network. [36]his study presents such a model, which is only based on the material properties.A dynamical variability allows the oxygen vacancies to jump to different positions and to directly map this local configuration change to the current transport.After presenting the model, the predictive power of the model is quantitatively tested using the experimental data of an HfO 2 /Pt reference system.By these means, the multilevel switching is investigated under the perspective of variability.The here-developed model is then applied to a neuronal network, by means of the underlying matrix multiplication, to test the influence of the variability on this highly interesting application.The influence on the performance of the network is measured.The simulation results provide a physical emulation for the considered network and it allows to explore certain physical limits.

Physical Principles
The development and analysis of a dynamical and scalable model requires a short review of the electronic conduction in VCM devices based on the physical insights of the atomistic ab initio simulations in ref. [37].The conduction is based on tunneling through the Schottky barrier.The tunneling length is determined by the depletion zone of the Schottky barrier. [37]Changes in the vacancy concentration vary this depletion zone and, therefore, the tunneling distance.Thus, the switching is invoked by changing the length of the depletion zone, due to a modification in the oxygen vacancy concentration.Two different tunneling processes, called type 1 and type 2, can be identified based on the energetic gap between vacancy state and CB minimum (CBM).These two types correlate with the two species of oxygen vacancy states found by Linderälv et al. [38] In their work, s-type or a p-type oxygen vacancy defect states were identified.The s-type vacancy category is identified by a deep defect state with respect to the CBM.The p-type defect states lead to a shallow defect state.Note, that both are donator-like defect states.[41][42][43] At high defect concentration, the defect state orbitals start to overlap.As a result, the sharp defect states broaden over an energy range to form a sub-band.The schematic density of states (DOS) is sketched in Figure 1c.In total, the number of states in the sub-band is fixed to the number of defect states.This is small compared to the number of states within the CB.
Type 1 conduction is found in oxides with a shallow defect level.Here, the electrons cross the Schottky depletion zone by tunneling from or into the CBM, as sketched in Figure 1a. [44,45]he shallow defect states allow an easy thermal excitation into the CB.The CB contains a high number of states at each energy level.This allows the electrons to tunnel from continuous energetic positions.This is visualized by the spectral current, which describes the current fraction at each energy level (current energy distribution), as shown in the right part of Figure 1a.Previous simulations show that for the type 1 conduction the spectral current forms a peak, such that the most of the current arise at a tiny energy range around the peak maximum. [44]This is shown as a green line in Figure 1a.Two mechanisms are responsible for the peak formation.At energies lower than the peak energy, the tunneling probability is limited by the long tunnel distance.At energies higher than the peak energy, the current is limited by the reduced thermal excitation of electrons into higher energetic levels.
The type 2 conduction is seen for oxides with a large gap between defect state and CBM.Due to the lower probability of thermal excitation into the CB, the main conduction mechanism is the tunneling directly between vacancy defect states and metal electrode.Thus, the Schottky barrier is overcome by a trap-tointerface limited transport.The position of these states is found at the Fermi edge.The type 2 conduction is sketched in Figure 1b.The type 2 conduction is limited by interface-to-trap transport, as suggested in other studies.[37,46].The number of conducting states in type 2 devices is limited to the number of oxygen vacancies, see DOS in Figure 1c.These are significantly less compared to the tunneling into the CB in case of the type 1 conduction mechanism.Figure 1b shows the tunneling represented by the single dashed orange lines, which connect the defect states to the electrode.The spectral current is given by single points at the energetic positions of the defect states.The amount of current over each spot is obtained by the combination of thermal excitation and tunneling probability, identical to type 1 conduction.Thus, it is expected, that the single spots in the spectral current are distributed on a comparable peak as for the type 1 conduction.This is indicated by the dashed green line.If the position of an oxygen vacancy changes, the tunneling length is also modified and the position of the corresponding spot in the spectral current is altered.In the sketch in Figure 1b, this scenario is plotted by black arrows.Thus, a relocation of an oxygen vacancy induces an abrupt current jump on the macroscopic scale as only a small number of vacancies take the major part in the conduction.
It is expected for type 1 conduction that the same small fluctuations occur due to single jumps of oxygen vacancies.The current response to this fluctuation is low, because of the many available tunneling states in the CB, which average out this fluctuation.Furthermore, the impact on the screening length, which determines the tunneling distance, is only small.As a result, the current strongly depends on the position of these vacancies (see black arrows in Figure 1b).Thus, a model has to be developed, to prove the concepts suggested in Figure 1b.Up to now, no comprehensive model has been developed that allows an analysis of the spectral current transport for type 2 conduction under dynamically applied voltage.Currently, the theoretical description is limited to the static ab initio simulation level at zero bias. [37]his problem is solved using the upcoming model, which leads to a nonzero bias description of the current and a scalability to multiple cells.This also allows to integrate the model in typically used circuit design simulators.Furthermore, it connects the application to the physics.
Based on this short review, the application of this model must be understood in the sense of VCM-type resistive switching devices.As a conclusion, the application within neuronal networks of the presented model must be understood in the scope of VCM devices.The variability of other types of electronic devices, like in studies [47,48], cannot been proofed within this study.Thus, the generalization is limited to VCM devices.

Schottky Contact Modeling
The defect chemistry of oxygen vacancies suggests that oxygen vacancies can be treated as twofold donors, which is true for both types of defect states. [31,32]Consequently, the donor concentration N D in the model is given by whereby [V O ] is the oxygen vacancy concentration.As introduced, these donors screen the Schottky contact at the active electrode.
During the switching process, the concentration of oxygen vacancies in the corresponding depletion zone is changed.The modeling of the Schottky depletion L dp zone is text book knowledge and can be described by with the build-in potential ψ Bi ðVÞ, the dielectric constant ε r , the vacuum permittivity ε 0 , the elementary charge e, and the voltage V (cf.ref. [49]).The build-in potential is given by the difference in the work-function between metal and the oxide plus the applied bias Here, ϕ M is the work function of the metal and ϕ Ox of the oxide.The voltage V is applied to the metal contact.The work function in the oxide is a result of the electron affinity χ Ox plus the difference between the CB and the Fermi level E t,depth = E c À E f .In the depletion zone, the CB has a parabolic band bending.Outside the depletion zone the bands are considered to be flat.This is sketched in Figure 2 and can be written as equation Here, the Schottky barrier is placed at position z = 0 m.Also, the defect states follow the band bending of the CB.Here, it is assumed that the vacancy states are distributed over a sub-band with an energy width of E SB;width .Note that it is named sub-band here, but it could be any kind of distribution of single vacancies due to their interaction.Even for a fixed defect state without subband, the model presented in this section can be applied in the limit of E SB;width !0 eV.However, due to the high concentration observed in VCM cells, a sub-band or energetic distribution is realistic and a valid assumption.The lower edge of the of the sub-band follows the parabolic dependence of the CB according to The position of the traps relevant for the transport is near the lower edge of the sub-band according to Equation (5).In Figure 2 the energetic landscape is sketched consisting of the CB and the lower sub-band edge.For the scope of this model the energy reference level is set to the oxide Fermi level of 0 eV.Obviously, the concentration is a macroscopic value for the portion of particles per volume.As we are considering small volumes, the number of vacancies can be small to approach realistic vacancy concentrations.As only integer numbers are physically reasonable, a discretization is required.For this reason, one needs to recalculate the number of oxygen vacancies N V in the investigated system via with the lattice constant a Ox and an integer number k.A switch is the switching area.For filamentary switching cells, this is the lateral extension of the filament at the active electrode.The sum of the depletion zone and k ⋅ a Ox defines the simulation region.Thereby the minimum simulation range is given by k = 0.The meaning of k becomes more relevant for the hopping of vacancies later on.
Here, it just defines the simulation range.Note that the ceil operation in Equation ( 6) rounds the number of oxygen vacancies to the next higher integer.The influence of rounding to the next integer is small but required to avoid the (unrealistic) case of zero vacancies in the depletion zone.
The oxygen vacancies are at equidistance distributed over the simulation range.Therefore, the position of the i-th oxygen vacancy is given by whereby r 0,i is between 0 to 1 and calculated by i/N V .From the position z t,i the energetic position E t,i of the i-th vacancy is calculated by E t,i = E SB,LE (z t,i ) according to Equation ( 5).In Figure 2 these traps are marked by green circles.Using the described positions and the energetic landscape presented in this section, the transport over each vacancy is calculated.

Current Calculation
[52][53][54] The principle of these simulations is that the oxygen vacancy-induced traps are filled with a certain probability.If a trap is filled it can release an electron into the electrode with a certain tunneling probability.If a trap is empty, it can capture an electron, which tunnels from the electrode into the trap.Vice versa, a filled state at the electrode interface can release an electron and an unoccupied state can capture an electron.Another way to describe the tunneling processes is by the coupling/attempt frequency R 0 of the electrons in the trap for the direct tunneling into the electrode.Those are often described as the coupling factor.The number of states in the metal which couple to the trap states is approximately given by the atom states at the interface and is called N(E).
The described process can be modeled by the tunneling rate from trap i to the electrode that is given by Free states at the metal electrode (8)   and the tunneling rate from the electrode into the trap i given by Here, f(E,E f,m ) is the Fermi distribution, E f , m the Fermi level in the metal, and Tr the tunneling probability.The number of free/ occupied states in the metal is given by the factors of occupation probability and N(E) as marked in Equation ( 8) and ( 9).The current for the two tunneling directions is then given by and Here, p t,i is the occupation probability of the trap i.The net current I t,i of the trap i is given by the difference between in and out flow The total current I of the system is obtained by the summation over all trap currents

The Occupation Probability
When a voltage is applied to a heterojunction, it is typically considered that the full voltage drop occurs at the interface.This leads to a split of the quasi-Fermi level at the oxide/semiconductor to metal interface.As a result, a constant quasi-Fermi level in the oxide E f,Ox and a constant quasi-Fermi level in the metal E f , m can be defined.Under a bias voltage V both Fermi levels are separated by eV.The oxide Fermi level E f,Ox is chosen as the reference energy at 0 eV; for E f,m this yields As a result, the oxide is in a quasiequilibrium and the occupation probability p t,i of the ith trap is calculated by

The Attempt Frequency
The current calculation requires an estimation of the frequency R 0 .The idea is to describe the electronic wave in the trap as a single-particle problem.The trap is considered as an isolated electrostatic potential, like in a quantum well.For this approximation, the trap potential must be sufficiently local.The trap potential is formed by the surrounding unit cell and the crystal.][57][58][59] This leads to the conclusion that the trap potential is sufficiently local to apply the approximation of a single-particle problem.Of course, the potential leading to the localized eigenstate has a complex unknown and materialdependent form.In this trap potential it is assumed that the electron oscillates in the ground state.For this approximation the electron must be localized for a sufficiently long time, which is higher than the attempt frequency.The trap potential is approximated as a box potential with a width of the oxide lattice constant a Ox .The ground state is given by with the electron mass m e and the Planck constant ħ.In this system the electron k-vector (momentum) is given by Then the velocity is calculated by v e = ħk/m e .The ground state can be understood as a standing wave which is generated by a left and right propagating wave.It gets reflected at the box potential sites.The frequency of reflection of the outgoing wave is interpreted as attempt frequency, which leads to the expression

Tunneling Probability
To calculate the tunneling, the popular Wentzel-Kramers-Brillouin (WKB) approximation is applied.From this approximation, the transmission coefficient can be calculated by defining the tunnel barrier for a particle at a specific energy.The transmission depends exponentially on the integration over the tunneling barrier.The tunneling conditions are independent for each vacancy position.Using the energy landscape derived in the Equation ( 1)-( 7), the tunneling position is defined by z t,i and the tunneling energy E t,i for the i-th vacancy is calculated.The second tunneling edge is found at the interface at z = 0 for each vacancy.
In a first-order approximation the barrier is treated as trapezoidal, which is sketched as red line for one selected oxygen vacancy in Figure 2. In equational form the trapezoid is described by The transmission according to the WKB approximation is then given by with the electron tunneling mass m t .Inserting equation (19) into equation (20) and executing the integration over the barrier yield Please note that according to Equation ( 19) the parameters a and b in Equation ( 21) are functions of V, z t , E t (meaning a(V, z t , E t ) and b(V, z t )).For reasons of readability, the dependence is not specifically mentioned.

Hopping Transport
The oxygen vacancy concentration is the state variable of the resistive switching cell.According to Equation ( 6), N V is conserved within one resistive state.However, there are manifold options to place the N V vacancies into the simulation volume of A switch ⋅ ðL DP ðVÞ þ k ⋅ a ox Þ, which varies as a function of voltage and time t.Due to thermal fluctuations, oxygen vacancies can hop from their initial positions to other lattice sites, see Equation (7).The hopping probability follows an Arrhenius-type law with the migration barrier E B,Vo for oxygen vacancies.The hopping rate of the vacancies is calculated by multiplication with the attempt frequency, which is given by the phonon frequency υ Ph Here, k B is the Boltzmann constant and T the temperature.The average number of jumps P Av in one simulation time step t step is given by P Av is a real number, but obviously the number of jumps must be an integer.Therefore, P Av should be interpreted as jump probability and transformed to a physically meaningful integer.For this goal P Av is split into P total is the total number of full jumps and is calculated by a floor operation The fractional part P frac is the rest probability that another vacancy jump occurs on top to the P total jumps.The jumps of the i-th oxygen vacancy P V,i are rounded to the next higher integer with a probability of P frac and reduced to the next lower integer (including zero) with probability of 1 À P frac .This leads to the formulation The P V,i jumps could represent jumps back-and forward with equal probability and partly amortize the total displacement.Thereby r i 's are random numbers in the range 0 to 1.The net displacement DP V,i of the i-th vacancy is given by the random walk DP V;i ¼ with the random number q k,i ∈ fÀ1, 1g.The jump distance of the vacancies is given by the lattice constant, such that the updated position of the oxygen vacancy i after the time step t step is According to Equation ( 28), oxygen vacancies could hop out of the simulation range.This could theoretically reduce the vacancy concentration, but this is physically not realistic, because also vacancies from outside of the simulation range may hop in.As a result, vacancies, which leave the system according to Equation ( 28), are replaced.Therefore, the new vacancies are placed on a random position in the simulation range.

Model Results and Quantification
The model described by Equation ( 1)-( 28) connects the physical insights from the ab initio simulations with the compact model level.Therefore, the model is independent from a specific device and does not require empirical input data.Consequently, the model could be fed with literature data for the material properties of the simulated device.Afterward, the predictive power of the model is tested by comparing it to independent experimental data for the high resistive state (HRS).This prevents that other conduction mechanisms compete with the resistance of the cell.For this proof HfO 2 is taken as relevant material with type 2 conduction.As a typical high-work-function electrode material Pt is used.The literature data are listed in Table 1 and the tunneling mass is approximated as the free electron mass.First a triangular current-voltage (I-V ) sweep is modeled and compared to HfO 2 /Pt-based reference system.The sweep data are taken from the published work of Hardtdegen et al. in ref. [60].The measurements were done with a sweep rate of 0.01 V/15 ms and measured in time steps t step of 15 ms.The measured resistive switching cycles are plotted in blue color in Figure 3.The device area is (100 nm) 2 .Thus, a filament radius in the range of 3-20 nm is reasonable.To model the I-V sweep an oxygen vacancy content of 3% and a filament radius of r fil = 10 nm are assumed.The percentage of oxygen vacancies in this manuscript refers to the maximum of possible oxygen which could be removed from HfO 2 .This is taken from the concentration of oxygen in the monoclinic cell of HfO 2 , here approximated by the eight oxygen atoms in the volume 1/a 3 HfO 2 .The filament radius and the vacancy concentration cannot be taken from literature.But these two parameters can be estimated to some Extracted from the calculated DOS for reduced HfO 2 and metallic Pt.
extent from the device size.Furthermore, transmission electron microscopy can be applied to visualize the single-device filaments, which are consistent to our simulation parameters. [61]he calculated I-V sweep using the same sweep rate and time steps as in the experiment is plotted in orange on top of the experiment in Figure 3.The calculated I-V curve is in good quantitative and qualitative agreement with the I-V characteristic of the HRS.Also, the variability and current jumps observed at this time scale are accurately predicted.As a result, the literature values as input parameters lead to the experimentally obtained current values over the full voltage range.The LRS of the cells could be determined by other internal mechanisms like the contact resistance.As those mechanisms are unclear the model might not be fully applicable in the LRS.
To trace the variability of a resistive state, read noise measurements are often applied, where one resistance state is written into the cell.Afterward, this resistance state is read out over time.This allows to visualize the current fluctuation over time.For an HfO 2 /Pt system, such a measurement is taken from the work of Schnieders et al. in ref. [35].The data has been measured with a sampling of t step = 0.13 ms for 2 s.To model this, the HRS vacancy concentration is kept at 2.5% as in the sweep, but the radius has been extended to 20 nm.The experimental noise current trace (NCT) measurement under a voltage signal of À0.25 V applied to the Pt electrode is redrawn as the black line in Figure 4. Using the same voltage and sampling time as in the experiment, three different NCTs are calculated.In Figure 4, the calculated NCTs are plotted in different colors in addition to the experiment.The analysis of the plot reveals a similar behavior as in the experiment for the current change over time.This analysis is further supported quantitatively by comparing the average and the standard deviation.The average current in the experiment is 7.36Â10 À6 A and the standard deviation 5.12Â10 À7 A. The simulations lead to average currents of: 6.61Â10 À6 , 6.74Â10 À6 , and 6.63Â10 À6 A. The simulated standard deviations σ are 4.55Â10 À7 , 6.31Â10 À7 , and 4.75Â10 À7 A. Hence, the model reproduces the experimental current fluctuation over time, such as the statistical parameters.
In conclusion, the model is able to reproduce the experimental results based on the literature values for the material properties.Therefore, the model can be used to investigate the internal physical processes responsible for the variability.

Model Analysis and Discussion
The predictive power of the model allows to explore the raised questions and assumptions of the introduction.To this end, the parameters of the result section are used (Table 1), with a filament radius of 10 nm and an applied voltage of À0.25 V.The resistive state is approximated with 2.5% occupation of oxygen vacancies.The introduced spectral current as energetically resolved current is an important measure.The oxygen vacancies in our model yield the dominating pathway for electron transport, which means that the current distribution is discrete at the oxygen vacancy energy levels.This can be visualized by plotting the current I t,i of the i-th vacancy over its energetic position E t,i .The E t,i À I t,i spectral current plot is on the right side of the Figure 5a-c at three different times.At the same time, on the left side, the vacancy positions z t,i and energy E t,i are plotted.The color code is utilized in the figure to distinguish and identify the vacancies over time, due to the hopping-induced position change.These three figures are snapshots for the NCT in Figure 5.The color code represents the time.In addition to the three snapshots, all states are visualized in a video format in the supplementary information.
To analyze the current distribution, a continuous dashed line is plotted, which connects the single-oxygen-vacancy currents on the right side of Figure 5a-c.This line can be understood as hypothetical continuous spectral current (current energy distribution).The form of the spectral current is a peak.The reason is the same as for type 1 conduction.In contrast to type 1 in the here-investigated type 2 conduction the current cannot flow at each point of the continuous hypothetical spectral current.The peak form rather describes the hypothetical current for the vacancies at a specific energy position (E t,i z t,i ).Thus, when an vacancy hops to an updated (E t,i z t,i) pair, the updated current over this vacancy is determined by the hypothetical spectral current.The hypothetical continuous distribution is a helpful tool to analyze the physical processes determining the current for each vacancy.The hypothetical current distribution forms a peak at a specific energy.For energies higher and lower than the peak the current is limited.This has two reasons as discussed in the introduction: For higher energies a thermal excitation of the electrons is required as the tunneling energies lay higher than the Fermi level.Consequently, the upper peak side stems from a thermal excitation limit.For energies below the Fermi level the tunneling distance and barrier increase, which leads to a decreased tunneling probability.Consequently, the lower peak is tunneling limited.These two limitations shape the peak-like distribution at the metal Fermi level.The broadening of these peaks gives the small energy range, where the most of the current is transported.Thus, most of the current is transported over the few vacancies, which are energetically located near the maximum of the peak.Therefore, the number of oxygen vacancies relevant for the current is very limited.
Even under a conserved concentration the configuration of oxygen vacancies is variable and they hop to different positions.The snapshots at different time stamps in Figure 5 are transferred to a video, which visualizes the change in the vacancy position (see motion of colored spots in supplementary information).The current transport is sensible to the few vacancies located in the energy range of the spectral current peak.The hopping and configuration of these few vacancies have a strong impact on the current, which leads to the fluctuations in the NCTs.
The NCT analysis is transferred to the previous voltage sweep simulation shown in Figure 3.The sweep in the applied voltage moves the position of the metal quasi-Fermi level.For further analysis, three reference voltages are chosen at À1 V, %0 V (À0.005V), 1 V. Snapshots at these voltages are shown in Figure 6.Also, for this sweep a video is given in the supplementary material.In the sweep, the current of the defect states forms a peak distribution for the same reasons as in the NCT analysis.The width of this peak is the transport-relevant energy range.The peak forms for negative voltages at the metal Fermi level and for positive voltages at the quasi-Fermi-level of the oxide.Thus, the peak forms at the energetically higher quasi-Fermi level.For the voltage near to zero the peak height reduces.The variability in the sweep simulated in Figure 6 has the same reason as for the NCT simulation.The permanent reconfiguration of oxygen vacancies due to the natural hopping leads to different currents.

Concentration-Controlled Resistance State
The switching in these cells is invoked by a voltage-stimulated change in the oxygen vacancy concentration.Thereby, different oxygen vacancy concentrations can be programmed.Thus, multilevel states in this manuscript are represented by a change in the concentration, while keeping the other parameters constant.Therefore, the current traces for 2000 concentrations between 2.5% oxygen vacancies to 8.75% are shown in Figure 7a.
The mean current I mean for various concentrations is plotted in Figure 7b.The results reveal a nonlinear increase of the current with increasing concentration.The conclusion is that when an application requires a linear spacing of resistance states the concentration must be changed in a nonlinear form.The plot might give the false impression that more conductive states show a higher variability compared to the high resistive states.However, rather the opposite is true.Plotting the relative error σ/I mean as inset in Figure 7b reveals that more conductive states show less relative variability.
In the measurements, the R HRS /R LRS ratio is often in the range 10-20 for HfO 2 .In contrast, the resistance change modeled here is higher than a factor of 70, which indicates that the considered concentration range considered is too broad.For the fit to the experimental NCT measurements, a concentration of 2.5% of oxygen vacancies was suitable (see Section 3).In the simulation shown in Figure 7, the concentration must exceed values of 4.43%, 5.05%, 5.55% for an R HRS /R LRS ratio of 10, 15, 20, respectively.

Statistical Influence on Neuronal Network Application
Based on the single-device behavior, an application-near mechanism of multiple interacting systems is investigated.The previous section has been focused on the single device, but does not scale to the variability of multiple interacting cells.The qualitative analysis reveals that the models correctly predict important measures such as variability.The implementation of the device physics is in accordance to the experiment and the ab initio results. [37,44]This conclusion allows us to use the model for the study of an ensemble of cells at the same time.
To emulate a real computing operation, the matrix-vector multiplication is tested.This is a fundamental operation used for the implementation of neuronal networks.To realize a matrix multiplication a dense layer is trained using FLUX in Julia. [62,63]The dense layer has no bias and a linear activation function θðxÞ ¼ x.Thus, the free adjustable parameters are the layer weights.To scale down the calculation resources, the MNIST data set is reduced by cutting the left/right/lower/upper three pixels and skipping every second pixel.This leads to 12 Â 12 pixel image, which is digitalized by a gray scale threshold of 0.1.The pixels are used as input vector to the dense layer with a size of 144 values.The index of the output vector with a size of 10 represents the handwritten numbers of the MNIST dataset.This leads to 1440 weights.An example is plotted in Figure 8a.The layer is trained in software with 64 Bit weights and plotted in a color code in Figure 8b (the color map legend for all figures is given in Figure 8i).The quality of the layer is measured by the false  prediction rate R err of 17.3% or a right prediction rate (1 À R err ) of 82.7%.The index of the maximum of the output vector is taken as predicted value.To investigate the influence of a lower number of bit states, the 64 bits are mapped to 1-8 Bit by rounding up to the next higher level for positive values and to the next lower level for negative values.The 1 bit representation and 4 bit representation are plotted as example in Figure 8c.The discretization of the weights to a lower number of bits turns out to reduce the performance.For the 1 bit representation the false prediction rate R err increases to %35%.As shown in Figure 8h, the false The cell is plotted for 1 bit and 4 bit.h) Rate of the false prediction for the software dense layer (blue dots), error rate for the mean conductance (green stars), and the overtime variation visualized as box plots.Each box plot is the statistical reference to the false prediction over time for the number of bits.Two explicit error rates over time are plotted for the 1 and 4 bit as insets.i) Color map for all dense layer plots in Figure 8.
prediction R err converges to the 64 bit limit at around 4-5 bit representations of the weight.[66][67] The matrix learnt above is a pure software solution.In the following this matrix will be rebuilt by the VCM cells.This rebuild will be compared to the previous computer solution, which is therefore referenced as software realization.
A possible hardware realization based on VCM cells of such a dense layer is sketched in Figure 8d, see refs.[68-71].For our simulation approach, we use the conductance (=inverse cell resistance) for the representation of the weights of the dense layer.In this setup the input vector contains voltages and the output vector currents.To realize negative values, each weight is represented by two VCM cells, whereby one cell is positively biased and the other one negatively.Thus, the current could be negative or positive depending on the conductance balance of the two cells, which represent one weight.The pixels are transferred to a voltage.Therefore, black pixels are identified as digital "0" with a voltage of V D = AE0 V and the white pixels as digital "1" with a read voltage of V D = AE0.25 V.The predicted hand written number is taken by the index of the maximum of the current output vector.For this proof of concept, an R HRS /R LRS ratio of the simulated cells is given by a factor of 10.
The VCM realization causes a variability in theoretically equal weights due to the cell-to-cell variability.In Figure 8g the mean values for one cell are shown as example.The mean values of this example are used to realize one weight, but the mean values vary from cell to cell.Another issue for the variability is the modification over time.The current trace over time in Figure 8g shows that the example weight starts to vary.This time variation is found for each weight represented by two VCM cells.As a result, the predictive power in terms of R err becomes a function over time R err (t).
Investigating the sample cell in Figure 8g reveals that there is no overlap over time for the 1 bit resolution.Thus, the two possible weights are distinguishable at all times.This changes when 4 Bits are stored in the cell.Here, the states start to overlap.Thus, the information of one weight at one time is not unique.If a fixed information should be stored in the bit states of the cell this raises a problem, because the information is not clearly defined.However, in the application of a neuronal network, it might be sufficient that the states are only statistically distinguishable, as the correctly predicted output of a neuronal network is a statistical measure.Thus, the variation of the individual network weights might be of secondary importance.To analyze this principle, the software reference network will be emulated by VCM cells using the presented model.
The network presented here has 144 weights per line.This requires 288 cells due to the positive and negative values.These 288 cells are simulated using different sets of random numbers and randomly permutated for each of the ten lines.The simulation parameters for the cells are taken according to the analysis in Figure 5, which is based on HfO 2 material parameters.The conductance reference level is given by the mean of the negative HRS conductance, which is modulated by an oxygen vacancy concentration of 10%.In Figure 8g, this level is plotted by the lower mean level.The R HRS /R LRS ratio is adjusted to 10.Thus, the LRS reference level is defined by ten times the HRS mean value (current level I HRS,mean ).The device is set into the LRS by increasing the oxygen vacancy concentration stepwise with a concentration step size of 0.094% of vacancies.When the concentration is high enough that the mean current exceeds ten times I HRS,mean the LRS is reached.In the same way the intermediate levels are reached to realize the multilevel states.Depending on the number of bits the R HRS /R LRS range is split into the corresponding number of states, which are the reference levels for the bit states.Afterward the bit states are adjusted in the same way.Starting from the HRS reference level the oxygen vacancy concentration is increased until it exceeds the reference level.The mean currents for the 1 bit and 4 bit states are adjusted and plotted in Figure 8g.The same procedure is done for the negatively biased cell, whereby the reference level remains the positive mean of the HRS.Therefore, also the HRS must be adjusted by increasing the oxygen vacancy concentration.Finally, the positively and negatively biased cells must be balanced to realize the weights of the dense layer.To represent negative weights the positively biased cell is set to the HRS, and the negative biased cell is put to the intermediate state or LRS.The intermediate state or the LRS represents the value of the dense layer.For positive values this procedure is applied vice versa.For example, in a 1 bit representation the software realized reference matrix in Figure 8c must be represented by the cells.For all negative weights the negative biased cell is set to the LRS and the positive cell to the HRS, which lead to a net negative current.In contrast for all positive weights the negatively biased cell is set to the HRS and the positive cell into the LRS.For 4 bit the cells are not set to the LRS; instead, they are balanced to a corresponding intermediate resistance state.The same is valid for the other bit representations of the VCM constructed matrix.
In summary, the described procedure maps the weights of the learned dense layer to the conductance of the resistive switching cells by increasing the number of oxygen vacancies.
Finally, the layer is represented by VCM cells, and the weights are connected to the natural variability of these cells.In Figure 8e the mean current of each weight is plotted for the 1, 2, and 4 bit case.The cell-to-cell variability leads to varying weights, in contrast to the fixed software reference.The variation over time is plotted in Figure 8f, where a snapshot of the weights is plotted at three times for the 1 bit case.In the supplement material a video of this is given including all time steps.This simulated hardware realization of the dense layer allows to calculate the error rate R err (t) as time-dependent quantity.In advance the mean values of the current states are used to calculate a time-independent R err .The error rate from the mean values is plotted as green stars in Figure 8h together with the software implementation at 64 bits resolution (blue dots), to visualize the difference.The time dependences for the 1 and 4 bit error rates R err (t) are plotted as insets.This analysis reveals that R err (t) scatters over time but shows an autocorrelation behavior.This autocorrelation becomes clear by comparing the error rates spaced by one time step, R err (t) and R err (t þ t step ).These two neighboring points scatter in a smaller range compared to the total fluctuation over time.This is caused due to the minor vacancy reconfiguration in one time step.Therefore, for future analysis this model allows to determine the critical time when the current response is not correlated to the starting point anymore.In this case the simulation can be accelerated, as the single steps of vacancies do not require a tracking of the ionic hops.The individual hops can be replaced by a suitable distribution function.To compare the scattering over time to the software and medial reference, the time scattering is plotted in a box plot format.The variation over time is given for all bits and plotted together with the mean (green stars) and software error rate (blue dots) in Figure 8h.
The median of the box plot is given by the orange line.Interestingly it fits closely to the mean value of the error rate of the simulated hardware realization.Consequently, the mean values could be applied in larger networks.Increasing the number of bits improves the predictive power of the network.Thereby, the hardware realization exactly follows the improvement found for the software reference.The error rate starts at high values, which reduces significantly until it converges to the performance limits.The point of convergence follows the software analysis.Thus, the limit is approximately reached using 4-5 bits.This limit is identical for software and hardware emulation.The performance limit of the hardware emulation, however, is notably higher than for the software reference.In conclusion, for the here tested classification problem, a number of 4-5 bits is sufficient.Therefore, it is sufficient to operate the VCM cells in 4-5 bit mode.The simulation reveals that different states variates over time and therefore, could not be clearly distinguished at each time.As a consequence of these variations, the performance limit of the hardware emulation is higher than for the software reference, leading to a performance gap.This performance gap is marked in Figure 8h.Clearly, this is a drawback for the hardware realization.Whether this drawback will be significant depends on the individual application of the individual network.
From the upper analysis the advantages and disadvantages of the VCM-based matrix multiplication are contrasted.First, the broadening of the boxplot analysis emphasizes that the outcome is time dependent and nondeterministic.Thus, only certain statistical values of the network are reproducible for a certain time span.Second, the software solution reaches a better performance limit than the VCM implementation.The reason for this performance gap is the intrinsic variability of type 2 ReRAM device.This performance gap, however, seems small compared to the advantages of the hardware realization presented.The total number of floating-point operations required in the software solution is massive, as the underlying matrix multiplication is based on a product of many floating point numbers and, on top, the sum for each line must be calculated.Additional features like a more complex activation function would further increase the number of floating-point operations.Those operations can slow down the calculation due to sequential execution.Furthermore, the calculation has a high-energy requirement.The hardware realization allows a massive parallelization as the matrix multiplication and sum over the lines is executed in one parallel step.This parallelization also significantly reduces the power consumption of the network as the massive number of floating point operations are surplus.As most neuronal networks are not perfect and have a statistical error quote, the additional error quote from the ReRAMs might be negligible.Thus, it is likely that most performance drawbacks will be of minor importance for the future hardware realizations of neuronal networks.In the future the developed model needs to be applied to more advanced neuronal networks, with higher number of layers.

Summary and Conclusion
The presented study gives a short physical review of the two dominating conduction mechanisms in VCM-based ReRAM devices.It has been explored that one of the conduction mechanisms reveals a higher variability, due to small fluctuations.To test this hypothesis, the study presents a dynamical model describing this conduction type 2. The model is based on a trap to electrode tunneling over individual oxygen vacancies.The total number of oxygen vacancies is calculated from the concentration.These vacancies are randomly placed in the simulation range.The simulation range is calculated from the Schottky depletion zone, which determines the tunneling barrier for electrons.The physical description of this process presented here enables to calculate the current over each individual vacancy.This current is highly dependent on the (energetic-)position of the vacancy.Experimentally, observed current fluctuations can be explained by relocations of oxygen vacancies over time, as the vacancy undergoes thermal hopping processes over the migration barrier.
The developed model only relies on the material properties.Using HfO 2 /Pt as a reference material system with corresponding literature values for the materials, the calculated current traces and their variability can be compared to experimental literature data.The model predicts the current evolution over time and the current response to a quasi-static voltage sweep qualitatively correct.This predictive power allows to explain the experimental current by the transport over distinguishable vacancies.The highest portion of current flows over a couple of vacancies in a specific energy range.The fraction and position of oxygen vacancies located in the energy range are most important for conduction.Thus, the current is extremely sensitive to relocations of this small fraction of vacancies.Furthermore, the mean current and its standard deviation over time for several oxygen vacancy concentrations are analyzed to investigate multilevel behavior.For the application-near analysis a dense layer of a neuronal network is simulated.The dense layer is trained in software and used as reference.Afterward, a hardware realization of this layer is emulated by simulation of multiple VCM-resistive switching cells.Those VCM cells act as weights of the dense layer.This leads to the variation of the weights over time and from cell to cell.The performance is then compared to the reference software dense layer.Two main conclusions can be drawn.First, the false prediction rate varies over time and is not deterministic.However, the variation is significantly smaller than the error itself for all tested bit states.Second, the variability of the cells leads to a performance gap, which describes the offset to the theoretical optimum of the software solution.This performance gap cannot be mitigated by increasing the resolution of the devices.This accuracy offset is overcompensated by the benefit in throughput, as the software realization performs the computing operations sequentially, while the VCM-based dense layer enables fully parallelized operations.
The model and conclusions derived within this manuscript pave the way for further investigations.Thus, the influence of more advanced networks must be investigated.The presented model is written in a script language and is, therefore, easily transferable and also applicable to simulators, which allow scripting language inputs.

Figure 1 .
Figure 1.Sketches of the different types of conduction mechanisms that involve oxygen vacancy defect states.a) Type 1 conduction is based on tunneling into the CB.The sketch shows the CB of the oxide as a black line and the shallow defect states as a gray line.The green box represents the active metal electrode.Upon applied voltage V, the Fermi level, represented as dashed line, splits.The tunneling path is represented by the orange dashed arrow.The continuous green line describes the energetically resolved current.b) Type 2 conduction is characterized by trap states, which are deep in the bandgap.The single red spots describe individual oxygen vacancies.Instead of a broad tunneling arrow, several single-orange tunneling arrows are plotted.The energetic distribution is plotted in a continuous manner, as a dashed green line.The tunneling can only occur at the single spots of oxygen vacancies on this dashed line.The hopping of a vacancy is plotted by the black arrow.This changes the (energetic) position and therefore, the current transport at the position of the dashed line.c) Sketch of the DOS for oxides with significant oxygen vacancy concentration.The red range describes the defect states, which are smaller and less broad compared to the blue sketched CB states.The separation of the red and blue zone describes the defect states.d) Representative sketch of the I-V sweep measurements of the two types of conduction sketched in (a,b), taken from ref. [37] (adapted under the terms of a Creative Commons Attribution CC-BY 4.0 License.Copyright 2022, The Authors, published by ACS).

Figure 2 .
Figure 2. Picture of the energetic landscape used in the simulation model.The green circles represent the oxygen vacancies.The red line describes the approximated tunneling barrier of one exemplary oxygen vacancy.The dashed orange line is the tunneling path for electrons from and into the exemplary vacancy.The defect states for the tunneling are built from the surrounding metal states, which are near to the vacancy position.

Figure 3 .
Figure 3. Simulation of an I-V sweep of an HfO 2 /Pt ReRAM device in HRS and comparison to switching data of Hardtdegen et al. (ref.[60]).The simulation does not include a change in the vacancy concentration.Therefore, the orange line does not follow the SET process.

Figure 4 .
Figure 4. NCT simulation in comparison to the experimental data.The experimental data are taken from ref. [35].

Figure 5 .
Figure 5. Analysis of an NCT at different times.The color of the current spots in subfigure (d) describes the time evaluation.The subfigures (a-c) are snapshots of the oxygen configuration at the three different times steps 0, 0.5, and 1 s.On the left side the position and energy of each oxygen vacancy are plotted.The colors of the dots represent one specific vacancy.At each time the colored spots moved to different positions.The right side of the plot describes the current over the individual oxygen vacancies.The dashed line is a hypothetical energetically resolved current distribution, which forms a peak.The colored spots describe the current over the specific vacancies.

Figure 6 .
Figure 6.Analysis of the sweep simulated in Figure 5 at three different voltages.

Figure 7 .
Figure 7. Resistance and variability for different oxygen vacancy concentrations.

Figure 8 .
Figure 8. a) Digitalization of and transformation to voltages of one example MNIST-figure (left site).b) A software trained dense layer consisting of 144 pixels as input and ten digits as output (right site).The weights are stored in 64 bit and color coded according to (i).c) Software mapping of the dense layer to 1 bit and 4 bit.d) Sketch of the hardware realization of a dense layer of ReRAM devices using positive and negative values.The weights are represented by the conductance G nm of the individual ReRAM cells.e) Simulation of the hardware realization sketched in (d) whereby the mean values of the current are plotted.This way the 1, 2, and 4 bit state representation is plotted.f ) 1 bit time evolution of the hardware realization.g) Example for the current simulation over time and the mean values.The cell is plotted for 1 bit and 4 bit.h) Rate of the false prediction for the software dense layer (blue dots), error rate for the mean conductance (green stars), and the overtime variation visualized as box plots.Each box plot is the statistical reference to the false prediction over time for the number of bits.Two explicit error rates over time are plotted for the 1 and 4 bit as insets.i) Color map for all dense layer plots in Figure 8.

Table 1 .
Material parameters for the simulation.