Stochastic circuit breaker network model for bipolar resistance switching memories

We present a stochastic model for resistance switching devices in which a square grid of resistor breakers plays the role of the insulator switching layer. The probability of breaker switching between two fixed resistance values, ROFF\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_\mathrm{OFF}$$\end{document} and RON\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_\mathrm{ON}$$\end{document}, is determined by the corresponding voltage drop and thermal Joule heating. The breaker switching produces the overall device resistance change. Salient features of all the switching operations of bipolar resistance switching memories (RRAMs) are reproduced by the model and compared to a prototypical HfO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {HfO}_2$$\end{document}-based RRAM device. In particular, the need of a forming process that leads a fresh highly insulating device to a low resistance state (LRS) is captured by the model. Moreover, the model is able to reproduce the RESET process, which partially restores the insulating state through a gradual resistance transition as a function of the applied voltage and the abrupt nature of the SET process that restores the LRS. Furthermore, the multilevel capacity of a typical RRAM device obtained by tuning RESET voltage and SET compliance current is reproduced. The manuscript analyses the peculiar ingredients of the model and their influence on the simulated current–voltage curves and, in addition, provides a detailed description of the mechanisms that connect the switching of the single breakers and that of the overall device.


Introduction
Resistance switching devices, able to change their resistance when a sufficiently high voltage is applied to their two terminals, have been the subject of a widespread interest in last decade mainly as a storage class memory elements with nonvolatile retention (named, indeed, RRAMs, standing for resistance switching random access memories) [1][2][3].Alternative applications of such devices have been pursued more recently thanks to the acknowledgement that RRAMs can be considered memristive devices [4].In particular, alternative computation schemes employing RRAMs have been recently proposed like the neuromorphic [5][6][7][8] and the nonvolatile stateful logic [9,10].In both applications, software simulations of large systems composed of a huge number of RRAMs are usually employed: compact models exist which are able to simulate the average behavior of resistance switching devices [11,12].On the other hand, RRAMs are known to show large operation variability with complex statistical behavior [13,14], which has been recently seen in a positive light, e.g., for the generation of random numbers [15]; for extending the number of accessible resistance states in electronic synapses [16]; and for stochastic computing [17].
As a matter of fact, the variability of RRAM switching is intrinsic to its operation principle: indeed, the most mature class of RRAM devices is the one that bases its operation on the inherent stochastic processes of formation and disruption of nanometric filamentary regions where metallic ions or defects, which locally reduce the stoichiometry of the insulating layer, activate the electric conduction from one electrode to the opposite [18][19][20][21][22][23].As the filament dimensions shrinks, according to both device and power scaling requirements [24,25], variability raises as a consequence.For this reason, applications basing their functionality on device variability are becoming more and more appealing for future high-performance and low-power computation schemes [14,26] and compact models capturing the operation variability of RRAM devices are gaining importance.
In this context, we propose a RRAM stochastic modeling in which the insulating material of a bipolar RRAM is assumed as constituted by a square network of bistable resistor breakers.The update rule of the resistance of the breakers is stochastic according to switching probabilities that depend on the temperature and voltage drop on each breaker.Such model is in agreement with the general interpretation of thermo-electrochemical processes at the base of bipolar resistance switching [24,27].
Models based on random circuit breaker (RCB) networks have been initially implemented with some success for describing the operation of unipolar RRAMs [28][29][30][31][32].Such models have been used to describe the filamentary nature of unipolar switching, and the stochastic aspect was only relegated to the initialization of the breaker network.Few attempts have been made for describing bipolar RRAMs with RCB networks: The insertion of an interface layer is described in Ref. [33], while a more complex modeling based on the solution of RCB network taking into account ionic drift and diffusion in a Monte Carlo scheme has been developed by Li et al [34], explaining transmission electron microscopy imaging of conductive filaments in HfO 2 -based RRAMs.In general, indeed, Monte Carlo simulations that take into account defect generation, diffusion and recombination possibly in combination with complex conduction mechanisms, like the multi-phonon trap-assisted tunneling, revealed very powerful in describing physical mechanisms of RRAM switching [18,19,35,36].
The novelty of the present work is to show a model based on the solution of RCB networks explaining bipolar resistance switching and which does not require neither the assumption of interface layers with different properties with respect to the bulk of the insulating layer, nor the treatment of complex conduction mechanisms and of ionic drift and diffusion.The simplicity of the simulation allows a low computational load for each device simulation, which is not possible for complex physical phenomenological models that have been developed for describing the switching mechanisms [24,37,38].
The proposed model is benchmarked to experimental results obtained on bipolar HfO 2 -based RRAM devices.In particular, the simulation of the following peculiar aspects of bipolar RRAM operation will be investigated: (i) the need of a preliminary forming step that brings the device from the initial resistance to a low resistance state (LRS); (ii) the possibility of partially resetting the device to a high resistance state (HRS) with a lower resistance value than the initial preforming one, thanks to the application of a voltage with opposite polarity with respect to forming; (iii) the occurence of a SET transition from HRS to LRS at lower absolute value of the voltage than that of forming but same polarity; (iv) the different switching dynamics of RESET and SET, featuring, respectively, gradual LRS to HRS decreases and almost vertical HRS to LRS transitions; (v) multilevel operation through tuning of RESET stop voltage and SET compliance current.
The manuscript is organized as follows: Sect. 2 describes the experimental details; Sect. 3 deals with the modeling formulation and gives an overall description of the modeling technique.Simulation results for the forming, RESET and SET operations are presented in Sect. 4 together with a qualitative comparison with experimental data.Sections 4.1-4.2are dedicated to the detailed explanation of the RESET and SET processes.SET and forming can be considered equivalent processes, and the detailed description of the SET process is sufficient to clarify the simulation results.Section 4.3 deals with some general aspects of the model that allow understanding the overall model operation and, finally, the conclusion follows.

Experimental details
Experimental data reported in the paper have been acquired on 50 nm Pt/5.5 nm HfO 2 /40 nm TiN devices and on 50 nm Pt/35 nm NiO/180 nm W devices with dimension of 40 × 40 µm 2 .Samples have been fabricated according to Refs.[39] and [40,41], respectively.Current-voltage characteristics have been measured applying a voltage to the device Pt top electrodes and keeping the bottom ones grounded through two source measuring units of a Keysight B1500 framework.Current ramps are used for forming operations; current compliances are imposed with the instrument in the SET operations.

Modeling formulation
A RRAM device is modeled as a square grid of m × n nodes connected by vertical and horizontal breaker resistors, as shown in Fig. 1a.Throughout the document, simulations are carried out with 20 × 30 matrices.The voltage is applied simultaneously to all nodes in the top line and all nodes in the bottom line are grounded (Fig. 1a) to emulate the equipotential top and bottom electrodes of the devices.The breakers can assume two resistance values, R OFF and R ON .The breaker switching is produced by temperature raise due to Joule heating and electric field; in particular, the switching probabilities of the ith breaker in an infinitesimal time interval dt from OFF to ON states and vice versa are defined as where f is the attempt frequency; E a is the activation energy; T i is the temperature of the ith breaker in the grid; and k B is the Boltzmann's constant.p i,OFF→ON and p i,ON→OFF differ in the sign of the energy lowering/raising terms of the ith breaker E i that takes into account the barrier reduction because of the local electric field and that is given by where α is the asymmetry factor which is normalized with respect to the number of vertical breakers in a column, through the factor m − 1; V i is the voltage drop on the ith breaker and q is the unit charge.The result is that, on average, a negative voltage applied to the device produces negative values of the terms V i and E i that favor the switching of the breakers from OFF to ON states.Correspondingly, a positive voltage applied to the device will promote the ON to OFF switching of the breakers.The quantity V i can be evaluated with only one calculation step as a solution of the Kirchoff's laws in the matricial form [42].
The temperature of the ith breaker, T i , is evaluated according to the Fourier's equation for heat exchange with a thermal bath [29]: where c and a are the heat capacitance and thermal conductivity of the breakers, respectively; T b is the thermal bath temperature; and R i and I i are the resistance of ith breaker and the current flowing through it.As discussed in Ref. [24], the transient time for heating can be estimated in the order of tens of ps and, therefore, the steady state solution of Eq. 3 provides a good description of the system.Hence, the temperature of the ith breaker is given by where the thermal conductivity a is considered proportional the breaker electric conductance through the parameter β (i.e., a = β/R i ) in agreement with the Wiedemann-Franz law.The quantity T i is governed by the power dissipated locally by each breaker and globally on the temperature of the bath, which is linked to the overall power dissipated by the RRAM device, P D = V D • I D , according to the following Eq.5: where T b0 is the room temperature and γ is the thermal resistance of the bath.In summary, the solution of the Kirchoff's laws allows the evaluation of the voltage drop and temperature of each breaker through Eqs.4-5 and finally of the switching probabilities of Eq. 1.Such calculations can be done in one step with no need of checking the selfconsistency of the solutions, which lowers the computational load and speeds the simulations up.
The response of the system to a pulse with voltage V D and time width tw is simulated according to the flowchart reported in Fig. 1b.At first, the resistor grid is initialized.For what concerns the forming simulation, all the breakers are initially in the OFF state except few of them which are randomly initialized in the ON state to emulate role of the intrinsic defectivity of insulating layers [21,23,43,44].The fraction of ON breakers with respect to the total amount is controlled by the parameter f ON .In the simulations reported in the present paper, the breakers are initialized to the ON state randomly within the matrix.Such simulated defect configuration can find an experimental counterpart in the benchmarked HfO 2 -based RRAMs in a situation in which localized defects are present in the oxide films, e.g., local non-stoichiometries in amorphous insulators [21,23] or presence of doping species [39,45].The initial configuration of the simulations of the present paper, thus, differs from a condition in which defects are aggregated, e.g., as at the grain boundaries of a polycrystalline film [22,23,46].The initial condition for RESET and SET simulations is loaded from the final output of previous simulations, e.g., initial condition for RESET process is loaded from a forming simulation, while the initial condition for a SET process is loaded from a RESET simulation.
After the initialization, the Kirchoff's laws are solved and the overall device resistance is evaluated.According to many experimental observations for forming and SET operations, the current must be limited in order to prevent the insulator breakdown [1][2][3]20].For this reason, for forming and SET, the overall current flowing through the device is compared to a compliance current I C , which determines the ending ('No') or the continuation ('Yes') of the simulation (Fig. 1b).For the RESET operation, the current compliance is not needed and not implemented.If the simulation continues, the switching probabilities for breakers in the OFF state to switch to the ON state ( p i,OFF→ON ) and vice versa ( p i,ON→OFF ) are evaluated for a determined simulation time step dt.The breakers to be updated, either from OFF to ON or from ON to OFF states, are drawn according to their switching probabilities.If the simulated time exceeds the time width of the voltage pulse, the simulation ends ('No'); otherwise, the grid is solved again ('Yes') and the simulation proceeds until the voltage pulse runs out or, in case of forming and SET, until the compliance current is reached (Fig. 1b).
As a remark, it is worth stressing the distinction between the resistance of the breakers, which are simulation parameters indicated by either R OFF or R ON , and the resistance of the RRAM device R D , which can switch between a LRS and HRS.The resistance values of LRS and HRS are the results of the simulation; they depend on the applied voltage, the time width of the voltage pulse and the current compliance, and because of the stochasticity of the model, they vary from a simulation to another one for the same settings.
Finally, it must be pointed out that the current-voltage characteristics of the simulated device are linear in both the HRS and LRS states, because R OFF and R ON are fixed and constant parameters.Of course, the implementation of a voltage dependent resistance of the breakers would burden the computational load of the model that would require a self-consistent solution of the Kirchoff's equations in the breaker matrix.A refinement that can improve the model accuracy without impacting the computational efficiency is the introduction of a field enhancement factor as discussed for example by P. Y. Chen and S. Yu [11].

Simulation results
The model described above is run with the parameters listed in Table 1, and the representative results are summarized in Fig. 2. Figure 2a reports the simulated forming, RESET and SET operations driven by voltage staircases (time width tw of each step is 1 ms).Forming operation (for negative voltage) is visible as a sudden vertical transition up to the reaching of the compliance current I C .The subsequent RESET gradually brings the device to the HRS.As usual for filamentary RRAMs, the forming operation is only partially reversible [1][2][3]; indeed, the resistance of the HRS is lower than the preforming resistance.The SET operation occurs as a steep current increase at a negative voltage that is lower than that of the forming transition.For comparison, experimental measurements performed on Pt/HfO 2 /TiN devices are reported in Fig. 2b with forming at high negative voltages, gradual RESET process for positive voltages and vertical SET transition for low negative voltages.For obtaining a more precise current control during forming, a Voltage (V) Stop Voltage (V) current-driven sweep is performed, which results in a forming transition that appears horizontal, i.e., at constant current, in Fig. 2b.On the contrary, the simulation is performed with a voltage controlled sweep, in which the forming transition is vertical, i.e., at constant voltage (Fig. 2a).It can be noticed that the simulations capture the salient features of the experimental switching operations.Both simulated and experimental results agree on the fact that RESET transition is gradual and the higher the RESET stop voltage, the higher the HRS resistance.On the other hand, increasing the RESET stop voltages reduces the rate of resistance increase and finally leads to device breakdown, which can be recognized as a sharp transition toward low resistances as shown in the insets of Fig. 2a, b for simulated and experimental results, respectively.Finally, for both simulated (Fig. 2a) and experimental results (Fig. 2b), the SET transitions occur in correspondence of the vertical current increase toward the set I C value for voltage values that are lower than those of forming but with the same polarity.
The simulated multilevel programming operation is reported in Fig. 2c in comparison with the experimental one in Fig. 2d.Starting from a HRS, the resistance is continuously lowered through a sequence of SET operations with increasing I C (colorscale lines for negative voltages in Fig. 2c, d for simulated and experimental results) up to reaching the LRS.The sequence of SET operations constitutes a partialization of the SET operation that leads directly from HRS to LRS indicated as the black envelop line in both Fig. 2c, d [47].The multilevel programming operation in the RESET process is obtained by tuning the stop value of the voltage sweep (colorscale lines for positive voltages in Fig. 2c, d), which corresponds to a parceling of the gradual RESET transition as indicated by the black envelop line that leads directly from the LRS to the HRS [47].Though the simulation parameters have not been calibrated to fit the experimental data reported in Fig. 2, it is interesting to notice the good agreement between the resistance obtained after a RESET operation, R H RS , at different stop voltages for simulated (line) and experimen- tal results (symbols) as shown in Fig. 2e and the agreement between the resistance obtained after SET operations, R L RS , performed with different compliances I C for simulated (line) and experimental results (symbols) as shown in Fig. 2f.
In the following subsections, the simulations of the RESET and SET processes are analyzed in details to evidence the way in which the probabilities of breaker switching are affected by the voltage drops on them and their temperatures.A following subsection is dedicated to the roles of the competing breaker switching from ON to OFF and from OFF to ON states and the term of heating of the thermal bath γ , which significantly affects the simulation results and hides an interesting physical meaning.

RESET process
In this section, the details of the RESET process are described.Figure 3a shows the forward sweep of the RESET operation.The initial condition (LRS) is obtained from the forming simulation of Fig. 2a.According to Fig. 3a, as the voltage is increased, the current raises up to point 'A', where the RESET process initiates and the current starts decreasing.Point 'B' corresponds to a situation in which the RESET process is almost completed (close to the HRS).The maps of the breaker resistances in point 'A', at the beginning of the RESET process, and 'B', close to its end, are reported in Fig. 3b and c, respectively.Since the lateral dimension is comparable with the vertical dimension of the cell size (30 against 20 breakers, see Table 1), all the device volume is affected by the forming process, i.e., black segments, corresponding to breakers in the ON state, are visible throughout the entire map (please refer to the colorbar on the right of Fig. 3c).Simulations performed on larger breaker matrix result in filamentary-shaped conductive paths extending from top to bottom interface as expected from RRAM devices.Of course, increasing the grid dimension weighs the computational load down and a trade-off between simulation time and desired accuracy must be looked for.
The following panels (d-i) of Fig. 3 report the analysis of the quantities that influence the switching probabilities of the breakers p i,OFF←→ON , i.e., the energy lowering/raising terms E i and the thermal energy terms k B T i .These quantities are evaluated for all the breakers in the grid and are analyzed in the following, taking into account their average values to catch the overall trend as a function of the applied voltage.Furthermore, the fraction of breakers in the grid for which such quantities assume values above some specific thresholds is also evaluated in order to get pieces of information about localized processes inside the grid.
The analysis of the energy lowering/raising terms E i that enter Eq. 1 is reported in Fig. 3d-f: Fig. 3d shows the fraction of breakers with E i < −0.05 eV (solid line) and with E i > 0.05 eV (dashed line). 1 It must be reminded that positive values of E i favor OFF to ON switching, while negative values favor ON to OFF switching according to Eq. 1.The fraction of breakers with E i > 0.05 eV (dashed line) is larger than those with E i < −0.05 eV (solid line), though they differ by only about a factor 2. This justifies the fact that most breakers switch from ON to OFF, thus producing an overall resistance increase.In the maps of Fig. 3e, f, the distribution of E i values in the grid is reported, according to the colorbar at the right side of Fig. 3f.Negative values (favouring OFF to ON switching and SET) are hardly visible and they are sparse, while positive E i values are localized in special region, where switchings from ON to OFF state produce an interruption of the conductive paths between top and bottom electrodes, as the RESET proceeds.For instance, a visual correspondence exists between yellow regions in Fig. 3e, f (high positive values of the energy lowering terms) and white regions in Fig. 3b, c (breakers switched to OFF state).
The second factor that affects the switching probabilities of Eq. 1 is given by the quantities k B T i , which are analyzed in Fig. 3g-i.The mean value of all k B T i (k B T i , dashed line) and the fraction of breakers with k B T i > 0.06eV (solid line) are reported in Fig. 3g, left and right axes, respectively.It is interesting to notice that the average temperature remains almost constant from point 'A' on, because the RESET process keeps on reducing the current flowing through the device, while the voltage increases.However, the heat dissipation becomes localized in the regions in which the current flow from top to bottom electrodes is interrupted.Indeed, the k B T i map in Fig. 3h in correspondence of point 'A' is rather uniform while the one in Fig. 3i in correspondence of point 'B' shows high-temperature regions where large positive voltages drop (compare with Fig. 3f) and breakers are in the OFF state 1 the fraction of breaker is evaluated with respect of the total amount of breakers in the grid, i.e., (m − 1) • n vertical and m • (n − 1) horizontal breakers.
(compare with Fig. 3c).This picture is coherent with the fact that, even though the average temperature remains constant, the fraction of breakers with a significant thermal increase (k B T i > 0.06 eV) still increases with increasing voltage (Fig. 3g, solid line, left axis).
The terms E i and k B T i influence the probabilities of breaker switching from OFF to ON and from ON to OFF, which are described in Fig. 3j-l and m-o, respectively.The average switching probabilities ( p i,OFF←→ON , dashed lines) and the fraction of breakers with p i,OFF←→ON > 10 −4 (solid lines) are shown in Fig. 3j and m for OFF to ON switching and vice versa, respectively.It can be noted that both average probabilities and fraction of breaker with p i,OFF←→ON > 10 −4 are slightly larger for ON to OFF switching (Fig. 3m) compared those related to OFF to ON switching (Fig. 3j), as expected for a RESET process.This fact is mostly determined by the lowering/raising terms E i .Both probabilities on average remain almost constant from point 'A' on.However, from the probability maps, it can be noticed that the switching probabilities are high in special regions at point 'A' (Fig. 3k, n), while in point 'B' the probabilities assume almost uniform and low values throughout the whole grid (Fig. 3l, o), which marginally affects the overall device resistance and determines a slowing down of the RESET process.
In summary, the RESET process occurs because the ON to OFF probability is larger on average moving from point 'A' to 'B.' Furthermore, the larger rate of resistance increase around point 'A' with respect to point 'B' is due to the fact that the switching probability in 'B' assumes uniform values in the grid while in 'A' concentrates in regions where the breaker switching significantly influences the overall device resistance.The raise of the switching probability is mostly due to the energy lowering/raising terms E i at the beginning of the RESET process (around point 'A') where the temperature is nearly uniform in the grid.

SET process
In this section, the SET process is described in a similar manner as for the RESET.The initial condition (HRS) is the one obtained from the RESET process of Fig. 3.In Fig. 4, all the quantities that govern the device switching from HRS to LRS are analyzed.
Figure 4a shows the current-voltage characteristic.It initially raises proportionally to the applied voltage up to point 'A' and afterward a rapid current increase is observed through point 'B' and until the compliance current is reached (LRS).The maps of the resistance values of the breakers in points 'A' and 'B' are reported in Fig. 4b, c.At point 'A', some elongated regions extending from top left to bottom right can be recognized, which limit the conduction of the current through the device and are the result of the previous RESET process.Despite the slight voltage difference of points 'A' and 'B', most of the white regions at point 'A' (OFF state, Fig. 4b) are filled with black segments, i.e., breakers in ON state, at point 'B' (Fig. 4c).
The analysis of the energy lowering/raising terms E i that enter Eq. 1 is reported in Fig. 4d-f: Fig. 4d shows the fraction of breakers with E i < −0.05 eV (solid line) and with E i > 0.05 eV (dashed line).The fraction of breakers with E i < −0.05 eV is larger than those with E i > −0.05 eV.This evidence justifies the fact that most breakers switch from OFF to ON, thus producing an overall resistance decrease.In the maps of Fig. 4e, f, the distribution of E i values in the grid is reported, according to the colorbar at the right side of Fig. 4f.In Fig. 4e, positive values (favouring ON to OFF switching and RESET) are hardly visible and they are sparse, while negative E i values are localized in the regions where accumulations of OFF breakers exist (please compare the blue regions in Fig. 4e with white regions in Fig. 4b).At point 'B' in Fig. 4f, instead, mixed negative and positive E i regions are visible that contribute to both OFF to ON and ON to OFF breaker switchings, still with a slight prevailing of negative E i values that allow the following step toward the device LRS in Fig. 4a after point 'B'.
The second factor that influences the switching probabilities of Eq. 1 is given by the quantities k B T i , which are analyzed in Fig. 4g-i.The mean value of all k B T i (k B T i dashed line) and the fraction of breakers with k B T i > 0.06 eV (solid line) are reported in Fig. 4g, left and right axes, respectively.Because of the low applied voltage, the k B T i values remain low until the SET transition is triggered.During the SET transition, the connection from top to bottom electrodes is instated and the Joule heating is almost uniform throughout the volume of the device as visible in the k B T i map of Fig. 4i.
The terms E i and k B T i influence the probability of breaker switching from OFF to ON and from ON to OFF, which are described in Fig. 4j-l and in Fig. 4m-o, respectively.The  4j and m for OFF to ON and for ON to OFF switching, respectively.It can be noted that, before the SET transition (point 'A'), both p i,OFF←→ON and fraction of breakers with p i,OFF←→ON > 10 −4 are slightly larger for OFF to ON switching (Fig. 4m compared to Fig. 4j), as expected for a SET process.This fact is mostly determined by the lowering terms E i .For voltages lower (as absolute values) than that of point 'A', the switching probabilities are very low as visible from the probability maps at point 'A' of Fig. 4k-n, but the switchings of few breakers produce an avalanche phenomenon that results in a large variation of the overall device resistance.Such simulated process is in line with the interpretation of the SET operation in filamentary RRAMs as triggered by a threshold switching event [24].At point 'B' (Fig. 4l), a significant raise of the switching probabilities from OFF to ON, p i,OFF→ON , with respect to point 'A' (Fig. 4k).On the contrary, there is no significant difference for the ON to OFF switching probabilities between point 'A' and 'B' (Fig. 4n, o, respectively), which remain very low in both cases.In summary, the SET process initiates at voltages producing low values of thermal heating and barrier lowering but as soon as the process is triggered both terms explode and cause a fast device switching from HRS to LRS.

Roles of competing breaker switching and thermal bath heating
In this section, two aspects of the model are analyzed: (i) the implementation of a competing breaker switching either from their OFF to ON state or from their ON to OFF state and (ii) the inclusion of a heating mechanism of the thermal bath.
In the previously described simulations, the competing switching of the breakers consists in the fact that for every operation both p i,OFF→ON and p i,ON→OFF are evaluated together with their relative update of breaker resistance, according to the flowchart of Fig. 1b.The overall process of resistance increase or decrease is governed by the voltage polarity through the different sign of the lowering/raising terms E i in the argument of the exponential of Eq. 1.Such model ingredient is fundamental for taking into account secondary effects like the possibility of device failure through breakdown during RESET operation when the current is not limited (e.g., see inset of Fig. 2a) [48], or partial resistance decrease (increase) during RESET (SET), which is one of the sources of device variability.On the other hand, such competition between tendencies to SET and RESET processes is rarely taken into account in models that aims at providing a description of the different phenomenology of the switching mechanisms for SET and RESET [24,37].Only few models explicitly take such competing switching tendencies into account, e.g., as Monte Carlo models that evaluate the probabilities of generation, diffusion and recombination of oxygen vacancies that affect the electron transport through the oxide of a RRAM device [18,35,36].
The mechanism of heating for the thermal bath is governed by the parameter γ (Eq.5), and it might sound a technical aspect, but it hides an important physical meaning that affects the RESET dynamics as it will become more evident in the following.
Figure 5 shows representative simulations of forming, RESET and SET processes with different combinations of competing breaker switching and thermal bath heating.Figure 5a reports simulations performed in the same condition as those used for Fig. 2a and with parameters listed in Table 1. Figure 5b reports the results of the simulations performed neglecting the heating of the thermal bath (γ = 0): the forming and RESET voltages shift to higher voltages (as absolute values) compared to simulation of Fig. 5a.Indeed, the RESET stop voltage has been extended to 3 V in order to obtain a resistance change, to compensate for the low temperature obtained with γ = 0.However, because of the high RESET stop voltage, the probabilities p i,OFF→ON become important and determine the decrease of the overall device resistance in the voltage range 2.5-3 V, as indicated by the arrows in Fig. 5b.Since, at the end of the process, the overall resistance change is low, a very low voltage with the opposite polarity is sufficient to trigger a SET transition.
Figure 5c reports representative simulations in which the competing breaker switching is not taken into account, i.e., the probabilities p i,ON→OFF and p i,OFF→ON are not evaluated in forming/SET processes and in RESET process, respectively.The parameter γ is fixed at the same value as in Table 1.By comparing Fig. 5a, c, no significant difference can be noticed.However, the inclusion of the competing breaker switching allows simulating secondary effects as device breakdown (as in the inset of Fig. 2a) or additional variability, e.g., small resistance increases and decreases during the gradual RESET process.
Finally, Fig. 5d reports representative simulations in which both the competing breaker switching and the thermal bath heating are neglected.The interesting aspect that differentiates the results of Fig. 5d from the other simulations in Fig. 5 is that the RESET occurs as a sharp transition from a low to a high resistance state in case γ = 0. Since competing switching from OFF to ON state is neglected, the overall device resistance can only increase in the RESET process, contrarily to the RESET of Fig. 5b, where it first increases and then increases as previously discussed.The sharp resistance increase during RESET operation is typical of unipolar devices and is considered to be due to a mechanism of self-accelerated filament dissolution [49].According to such a model, the conductive filament responsible for the low resistance state of a unipolar device is dissolved because of temperature raise: As soon as the filament reaches the critical temperature, it starts dissolving and reducing its size, which increases the current density and increases the Joule heating.In this way, the temperature in the filament increases in a positive feedback loop which promotes a fast filament dissolution [49] .This process is justified if the heat exchange between the filament and the remnant insulating film is negligible and/or the filament is highly conductive [50].Furthermore, in unipolar devices, SET and RESET mechanisms are distinct and do not usually compete for the similar voltages values, which is coherent with the fact that in simulations reported in Fig. 5d the competing switching is neglected.For comparison, a representative RESET operation of a NiO-based unipolar device is reported in the inset Fig. 5d, together with the relative SET operation at the same voltage polarity.
Abrupt RESET transitions are often reported also for conductive bridge RRAMs (CBRAMs) employing low thermal conductivity materials (like the chalcogenides) in which the low resistance state is due to the formation of a conductive filament constituted by dense accumulation of metallic ions [51].In agreement with such experimental observations, neglecting the heating of the thermal bath (γ = 0) is mathematically equivalent to considering only the local contribution to the temperature raise, given by the quantity in Eq. 4 specific for each breaker and neglecting every variation of the bath temperature.Furthermore, also for CBRAMs showing sharp RESET, the RESET voltages are very different from SET voltages as absolute values [1,51], which justifies the fact that competing breaker switching is not important.
In general, also for oxide RRAMs, abrupt RESET transition has been sometimes found and associated with peculiar filament configurations or compositions [50,52].With respect to such experimental evidences, Fig. 5 demonstrates that the implementation of a competing breaker switching from OFF to ON or vice versa and of the thermal bath heating through the parameter γ allow tuning the RESET dynamics from abrupt to gradual.As a matter of fact, such expedients were not taken into account in previous versions of RCB models, which were dedicated to the simulation of unipolar devices with sharp RESET transitions, in agreement with the above considerations [29].
As an assessment of the presented simulation results, a comparison of the obtained temperature values with literature data is reported here.The RESET process is characterized by a temperature of 630 K as a mean value over the breaker matrix which reaches a maximum at point 'A' of about 700 K for some specific breaker.Conversely, the SET process is characterized by correspondingly lower temperatures of 560 K on average which reaches a maximum at point 'B' of about 580 K.Such temperature values are in line with literature results which are highly variable and range for RESET from 500 to 800 K [53][54][55][56] and for SET from 310 to 800 K [53,[55][56][57] for finite element, as well as, for Monte Carlo simulations.

Conclusions
Finally, in the present paper, a stochastic random circuit breaker model is proposed for the description of the bipolar operation of RRAM devices.The model considers that the breakers, arranged in a squared grid, possess a probability of switching which depends on the temperature and the electric field of the breakers themselves, in agreement with standard general interpretation of bipolar resistance switching effects.Temperature and voltage drops are evaluated in one step so that the computational load of the model is low and allows the simulation of a large number of devices or large number of resistance switching cycles.
The simulations are compared to experimental results obtained on a prototypical HfO 2 -based device and reproduce their salient features: the need of a forming step and alternate voltage polarities to SET and RESET the devices; the occurrence of RESET as a gradual transition and SET as an abrupt one.The evolution of SET and RESET processes has been described in detail to evidence the roles of temperature and voltage drops in determining the probabilities of breaker switching.Furthermore, the aspects related to competing breaker switching from OFF to ON and from ON to OFF and to the evaluation of the temperature of the thermal bath are described as additional degrees of freedom of the proposed model.In general, the model can be adapted to various material systems and can evaluate the time-dependent dynamics of the switching instead of the voltage dependent description reported in the present paper.
Given its simplicity and computational lightness, the model may result useful for the simulation of complex nanosystem architectures, e.g., alternative computing systems, in which many RRAMs, or even memristive devices, have to be simulated together with their stochastic operations.

Fig. 1 a
Fig. 1 a Sketch of a resistor grid with m × n nodes: voltage V D is applied to the nodes in the top line and the nodes in the bottom line are grounded.The overall current through the device is I D ; b flowchart of the RCB simulation for a fixed voltage value applied for a time interval equal to tw

Fig. 3
Fig. 3 RESET process: a I-V curve where points 'A' and 'B' are identified; b, c maps of the breaker resistances R i in points 'A' and 'B', respectively; d fraction of breakers with E i < −0.05 eV (solid line) and E i > 0.05 eV (dashed line) as a function of voltage (right axis); e, f maps of E i at points 'A' and 'B', respectively; g average temperature of the breakers k B T i (dashed line, left axis) and fraction of breakers with k B T i > 0.05 eV (solid line, right axis); h, i maps of k B T i in points 'A' and 'B', respectively; j average probability of switching from OFF

Fig. 4
Fig. 4 SET process: a I-V curve where points 'A' and 'B' are identified; b, c maps of the breaker resistances R i in points 'A' and 'B', respectively; d fraction of breakers with E i < −0.05 eV (solid line) and E i > 0.05 eV (dashed line) as a function of voltage (right axis); e, f maps of E i at points 'A' and 'B', respectively; g average temperature of the breakers k B T i (dashed line, left axis) and fraction of breakers with k B T i > 0.06eV (solid line, right axis); h, i maps of k B T i in points 'A' and 'B', respectively; j average probability of switching from OFF

Fig.
Fig. Simulation of forming, RESET and SET operations in different conditions: a simulation performed according to the model described above; b simulation in which the heating of the thermal bath is neglected, γ = 0; c simulation in which competing switching is neglected, i.e., forming and SET neglect possible ON to OFF switching and RESET neglects OFF to ON switching; d simulation in which both competing switching and heating of the thermal bath are neglected.The inset shows representative experimental results obtained on a Pt/NiO/W unipolar device: the downward arrow indicates the abrupt RESET transition

Table 1
The table reports the