TCAD Modeling of Resistive-Switching of HfO2 Memristors: Efficient Device-Circuit Co-Design for Neuromorphic Systems

In neuromorphic computing, memristors (or “memory resistors”) have been primarily studied as key elements in artificial synapse implementations, where the memristor provides a variable weight with intrinsic long-term memory capabilities, based on its modifiable resistive-switching characteristics. Here, we demonstrate an efficient methodology for simulating resistive-switching of HfO2 memristors within Synopsys TCAD Sentaurus—a well established, versatile framework for electronic device simulation, visualization and modeling. Kinetic Monte Carlo is used to model the temporal dynamics of filament formation and rupture wherein additional band-to-trap electronic transitions are included to account for polaronic effects due to strong electron-lattice coupling in HfO2. The conductive filament is modeled as oxygen vacancies which behave as electron traps as opposed to ionized donors, consistent with recent experimental data showing p-type conductivity in HfO x films having high oxygen vacancy concentrations and ab-initio calculations showing the increased thermodynamic stability of neutral and charged oxygen vacancies under conditions of electron injection. Pulsed IV characteristics are obtained by inputting the dynamic state of the system—which consists of oxygen ions, unoccupied oxygen vacancies, and occupied oxygen vacancies at various positions—into Synopsis TCAD Sentaurus for quasi-static simulations. This allows direct visualization of filament electrostatics as well as the implementation of a nonlocal, trap-assisted-tunneling model to estimate current-voltage characteristics during switching. The model utilizes effective masses and work functions of the top and bottom electrodes as additional parameters influencing filament dynamics. Together, this approach can be used to provide valuable device- and circuit-level insight, such as forming voltage, resistance levels and success rates of programming operations, as we demonstrate.


INTRODUCTION
In recent years, memristor devices have shown great potential for neuromorphic computing due to their resistive-switching dynamics and electrical behavior resembling that of biological synapses (Chua, 1971;Xia and Yang, 2019;Strukov et al., 2008). Memristors are resistors with memory whose resistance level can be controlled either through an applied voltage (i.e., flux-linkage) or current (i.e., charge-fluence). Oxide memristors tend to be voltage-controlled, having a metal-oxide-metal device structure identical to a parallel-plate capacitor. Since the oxide thickness tends to be thin (∼ 2-5 nm) (Pi et al., 2019) and the switching speed can be very fast (<1 ns) (Choi et al., 2016), a small amount of energy is required for programming resistance states. With these unique features, in addition to non-volatility, they have shown the great promise for building energy and area efficient memristive crossbar arrays (1T1R arrays) to form neural networks for a wide range of applications including robotics, computer vision, and speech recognition Li et al., 2018b,a;Hu et al., 2016). A 1T1R crossbar array ( Figure 1A) offers added benefits due to the use of a transistor in each resistive RAM (RRAM) memory cell. The transistor plays a major role in mitigating the sneak current path and programming disturbance associated with resistive (i.e., 1R) crossbar arrays (Manem et al., 2012;Yao et al., 2015). Furthermore, the transistor's gate terminal in the 1T1R cell allows for better control over the current through the memristive device. It also provides more resilience to the switching voltage magnitude and attains better uniformity (Liu et al., 2014).
Electroforming, a one-time forming or initialization process, is often required in transition metal oxide (TMO) memristors (Strukov et al., 2008), which have been widely used for memristive crossbar arrays, including those used in neuromorphic systems. However, for forming, the voltage often needs to be higher than the nominal supply voltage of modern CMOS processes, which causes significant design and integration challenges (Amer et al., 2017b,a). Memristors with high forming voltages require dedicated circuitry capable of tolerating such high voltage levels for executing the in-field forming. Furthermore, the area constraints associated with the in-field forming circuitry undermines the density benefits of the crossbars. In addition, such transistors are generally large compared to the regular devices to accommodate these high forming voltages ( Figure 1B). For example, for a 65 nm CMOS process (Beckmann et al., 2016;Amer et al., 2017a) used to prepare 80 nm × 80 nm memristor areas, the minimum length of the transistor used in the 1T1R cell could be as much as 0.5 μm to endure forming voltages up to 3.3 V. In contrast, the minimum length of the regular transistor used for peripheral circuitry is 60 nm with the nominal voltage 1.2 V or 1.0 V, depending on the process. Thus, researchers have focused on lowering the forming voltages to a level of operation that allows for better exploitation of memristive crossbar density (Govoreanu et al., 2011;Koveshnikov et al., 2012;Huang et al., 2013a;Chen, 2013;Kim et al., 2016;Amer et al., 2017c).
The electroforming process, in addition to reset and set operations, can be simulated from a condensed set of rate equations that define all possible changes of state of the system using Kinetic Monte Carlo (KMC). Rate equations used to model filament dynamics are typically based on the following physical transitions: oxygen ion (i.e., O 2− ) and vacancy (i.e., V 2+ O ) diffusion, and the generation and recombination of Frenkel pairs {V 2+ O , O 2− }. These have been implemented, successfully, by several authors for both 2D and 3D filaments with fitting capability to experimental data (Sementa et al., 2017;Aldana et al., 2018Aldana et al., , 2020Loy et al., 2020). The large difference in diffusion rates of oxygen ions and oxygen vacancies tends to favor filament growth along preexisting vacancy sites or positions in which the local electric field is high (e.g., grain boundaries or point defects) according to the thermochemical model of dielectric breakdown (McPherson et al., 2003). As a consequence, resulting filaments obtained from previous KMC approaches consist of positively charged oxygen vacancies resulting from the repetition of: 1) breaking Hf-O bonds and 2) the formation of Frenkel pairs consisting of nearly stationary oxygen vacancies and relatively diffuse, oxygen ions at interstitial sites. In other words, forming/set operations are thus determined by the local electric field-producing dendritic filament growth-whereas reset is determined by the coincidence of oxygen ion diffusion and recombination.
Despite growing evidence to the contrary, few modelling approaches allow the charge state of the oxygen vacancy (i.e., +2) to change during forming. In effect, oxygen vacancies are modeled as fixed charges for the purpose of determining filament evolution, yet as electron traps for the purpose of calculating current, which is based on nonlocal multiplephonon trap-assisted tunneling in which trap occupancy dynamics is fundamentally important. This inconsistency greatly limits the utility of existing approaches to provide increased physical insight into HfO x switching behavior-beyond that which existing compact models already provide (Bianchi et al. (2020); Yu and Wong (2010); Huang et al. (2013b); Guan et al. (2012b);Jiang et al. (2014))-which can be extended to device design, circuits and systems-level refinements (e.g., reducing forming voltage). The assumption of a static positive charge contradicts experimental evidence showing p-type conductivity in highly defective HfO x films (Hildebrandt et al., 2011)-suggesting that oxygen vacancies are deep acceptor-like traps (>3 eV from the conduction band edge). Moreover, ab-inito calculations have shown the thermodynamic stability of neutral and negatively charged vacancy states increases in conditions of electron injection (i.e., current flow) due to electron capture (Bradley et al., 2015). This is consistent with experimental work using in-situ TEM electron holography and EELS in which oxygen-vacancy filaments were observed, spatially, as regions of negative spacecharge (Li et al., 2017). Unlike previous KMC approaches, together these observations are, in fact, self-consistent with the physical assumptions used to model current flow in HfO 2 based memristive devices, in which, conduction occurs through a nonlocal, trap-assisted tunneling process involving electron capture and emission-appropriate for insulators having point defects (e.g., TaO x , HfO x , ZrO x , NbO x ).
Frontiers in Nanotechnology | www.frontiersin.org October 2021 | Volume 3 | Article 734121 Here, using a simple 2D model, we show that, in addition to the conventional set of rate equations (i.e., Frenkel pair generation/recombination and diffusion) filament evolution in HfO x can be modeled self-consistently as the result of bandto-trap electron capture and emission processes between the electrodes and oxygen vacancies ( Figure 2). In this way, the conductive filament consists of occupied oxygen vacancy electron traps, which lower their energy upon electron capture due to strong coupling between ionized defects and the lattice in HfO x (Huang-Rhys factor, S 17) as depicted in Figure 2. The primary benefit in this approach is that additional parameters associated with the electrodes (e.g., work function, effective mass) and those of the oxygen vacancy states (e.g., trap energy level, capture cross-section, thermal barrier and binding energy) are intimately linked to resistive switching behavior, as we show. Not only do these additional parameters provide more depth in terms of physical insight and modeling capability, they are readily accessible experimentally or through ab-initio estimates. Using TCAD Sentaurus (Synopsys, 2019), we demonstrate that the common forming, reset and set characteristics can be successfully reproduced and visualized. In particular, we show that certain regions within the filament have a negative potential-stemming from a negative space charge due to electron capture. This is consistent with recent experimental work describing the filament as a negative potential synapse (Li et al., 2017). Next, we couple our device model with a phenomenological compact model to bring in physics-based insights to the circuit-level simulation of a memristor-based synapse topology. Finally, to underscore the unique strength of our model, we investigate a device-circuit co-design strategy powered by Monte-Carlo simulations with different levels of initial oxygen vacancy volume faction.

Material Specifications and Device Geometry
A complete list of parameters used to specify the HfO 2 layer are provided in Table 1. Of these, the parameters related to hafniumoxygen bond energy and polarization were obtained from the thermochemical model (McPherson et al., 2003), assuming a 100% monoclinic phase composition. This makes clear the assumptions regarding crystalline phase of the HfO 2 thin film and the activation enthalpy required for breaking the hafniumoxygen bond-which differs due to differences in polarization and the number of bonds required to be broken about the Hf central atom. It should be noted that, in practice, a mixture of monoclinic and tetragonal phases are present in varying ratios-the monoclinic phase has a nominal breakdown field of 6.7 MV cm −1 whereas the tetragonal phase has a breakdown field of 3.9 MV cm −1 (McPherson et al., 2003). Empirically, the breakdown field is found to vary between 3 and 5 MV cm −1 (Sire et al., 2007) for atomic layer deposition (ALD) grown HfO 2 films on TiN of comparable thickness used in this work.
As indicated in Figure 3A, the memristor is modeled as a twodimensional top-electrode (TE)/HfO x /bottom-electrode (BE)  structure on a square grid. The HfO x thickness is 5 nm and the device width is 40 nm. The Ti (TE) and TiN (BE) electrodes are modeled as ideal, Ohmic contacts with a 0 Ω series resistance. The work function of the Ti and TiN layers were set to 4.33 and 4.5 eV respectively. The electron effective masses of the Ti and TiN layers were set to 3.2 and 2.0 respectively according to literature (Lima et al., 2012). Following the lattice-gas model, a grid point represents the smallest physical unit considered by this simulation, capable of representing either an empty "site" (for diffusion or the formation of a Frenkel pair), a positively charged oxygen vacancy (V 2+ O ), a negatively charged oxygen ion interstitial (O 2− ) or a negatively charged oxygen vacancy (V 2− O ) which also represents the conductive filament. Thus, field-independent transitions (e.g. Frenkel pair recombination) occur over nearest-neighbor distances whereas field-dependent transitions (e.g. Frenkel pair generation, ion diffusion) interact over many grid points through the screened Coulomb potential.
The initial state of the system can be defined by randomizing the location of oxygen vacancies and oxygen ions (needed to ensure charge-neutrality)-to represent an amorphous film. Alternatively, since it is known that ALDdeposited HfO 2 thin films exhibit a columnar grain 2 | Summary of transitions rates modelled using Kinetic Monte Carlo procedure and their parameters. Lattice coordinates are listed as relative positions to a given lattice point at (i, j) following the convention of the lattice gas model (Jansen, 2012

Transition Reaction Parameters
Oxygen Vacancy Diffusion Frontiers in Nanotechnology | www.frontiersin.org October 2021 | Volume 3 | Article 734121 morphology (≈8 nm grain size (Ho et al., 2003)), oxygen vacancies can be placed along grain boundaries due to the reduced formation energy-to represent a polycrystalline film.
The initial concentration of vacancies is determined by a variable volume fraction parameter which we nominally set to 0.0002 (i.e., 0.02 at. %). Our focus here is to demonstrate the key differences and advantages of our physical model incorporating additional electronic transitions to the constitutive rate equations describing filament dynamics and its implementation in TCAD Sentaurus.

Filament Evolution Under Voltage Stress
Filament evolution during forming, set and reset operations are described using a simple set of rate equations corresponding to the following physical processes, as outlined in Table 2: • Electron capture or emission by oxygen vacancies • Oxygen vacancy and ion diffusion • Frenkel pair generation and recombination (isolated bulk and nearest neighbor pairs) • Oxygen gettering by Ti These processes are implemented via a classical KMC selection algorithm applied to a chosen initial state of the system (top electrode, bottom electrode, oxygen ion, oxygen vacancy and filament), and updated in time according to Poisson statistics and the time scale of each selected mechanism. The details of the dynamical monte carlo algorithm and the meaning of simulation time (Fichthorn and Weinberg, 1991), and its application to the formation of 2D/3D conductive filaments (metallic and oxygen vacancy) has been discussed elsewhere (Sementa et al., 2017;Aldana et al., 2018Aldana et al., , 2020Loy et al., 2020). Here, we provide a minimal outline of the essential aspects, assumptions and parameters of the rate equations we've implemented. In particular, we highlight the physical assumptions that establish consistency between phenomenological models of filament evolution and models of electric conduction-needed for efficient and accurate device-circuit co-design.

Filament Changes Mediated by Electron Capture/Emission
We model filament precipitation (dissolution) as the net result of Frenkel pair generation (recombination) and electron capture (emission) by oxygen vacancy electron traps. For simplicity, we couple oxygen vacancy traps to the conduction band in the bottom-electrode, which permits straightforward evaluation of rate equations within the electron trap picture. Conventionally, capture and emission rates are defined in terms of a thermally activated capture-cross section, a tunneling coefficient, and a field-dependent trap barrier that depends on the relative energy difference between the trap level and the Fermi level of the bottom electrode. In other words, the forming, set and reset operations are described as band-to-trap (or trap-to-band) electronic transitions within the Wentzel-Kramers-Brillouin (WKB) approximation Eqs 1-3.
These expressions are derived in the Supplementary Material. Here, y t represents the y-coordinate of the trap relative to the electrode, y 0 is a parameter related to the wavefunction overlap between the electronic state in the trap and in the electrode. Within the WKB approximation, after applying the triangular barrier approximation for the bands we have the following: The quantity Φ be − χ ox represents the conduction band offset between the bottom electrode and the oxide, in terms of the work function of the bottom electrode Φ be and the electron affinity of the oxide χ ox . In Eqs 1-3, it is assumed that the electric field has a symmetric influence on transition rates, that is, the barrier lowering in the forward direction is equal and opposite that of the reverse in order to maintain steady-state equilibrium. The case statements in Eqs 1, 2 exist since the trap may be higher (i.e., E tf > 0) or lower (E tf < 0) than the Fermi level in the bottom electrode. Parameters which depend on the bottom electrode are the electron concentration, n be and the effective-mass, m * n,be , which enters through the thermal velocity (Eq. 5): We note that the time scale of electron capture and emission depends on the product of the carrier concentration in the electrodes, the thermal velocity and the capture cross section of oxygen vacancies as shown in Eqs 1-3 through a common exponential prefactor. Using values listed in Table 1, the ratio of the exponential prefactors for electronic (σ 0 n B Ev th ) and atomic processes (f ph ) is evaluated to be 8.25, so electronic processes are expected to occur much faster than atomic ones. However, the rate of electronic transitions also depends on the local electric field and the position of the trap relative to the electrode and so the above estimate only reflects those traps that are close to the bottom electrode. Therefore, the relative rates of electronic and Frontiers in Nanotechnology | www.frontiersin.org October 2021 | Volume 3 | Article 734121 atomic processes are expected to differ (generally reducing) as one moves from the BE to the TE under forming/set and from the TE to the BE under reset due to the change in voltage polarity.
The initial system consists of the electrodes and an initial concentration of positively charged, unoccupied oxygen vacancies. It should be noted that, although we have assumed a +2 charge state, the +2 oxygen vacancy is unstable in the presence of interstitial oxygen and/or conditions of electron injection. This is supported by ab-initio calculations suggesting that Frenkel pairs stabilize through the formation of neutral and/ or negatively charged oxygen vacancies-facilitated by electron capture (Bradley et al., 2015). Experimentally, this is supported by the observation of p-type conductivity in highly defective HfO x films, suggesting that oxygen vacancies interact and stabilize as deep acceptors (Hildebrandt et al., 2011) as opposed to shallow donors. Thus, the formation of conductive oxygen vacancy filaments can be regarded as thermodynamically driven by the increase in binding energy due to electron capture of single vacancies and expected to subsequently stabilize due to defect aggregation (i.e. filament growth). We account for these effects by lowering the energy level of an oxygen vacancy, E t upon electron capture, by an amount equal to the increase in binding energy of neutral and negatively charged vacancies (≈0.2-0.9 eV) as predicted by ab-initio calculations (Bradley et al., 2015;Sementa et al., 2017). In general, this is a lattice relaxation process involving the emission/absorption of multiple phonons as described by several authors (Englman and Jortner, 1970;Henry and Lang, 1977;Nasyrov et al., 2004;Nasyrov and Gritsenko, 2011). Here, for simplicity, we assume the lattice relaxation coincides with electron capture (or emission). The assumed electron capture and emission processes are illustrated in Figure 2.
The significance of this model of filament growth is that it is more consistent with the nonlocal trap-assisted tunneling processes associated with electron conduction, which we later implement in TCAD Sentaurus to calculate the current-voltage characteristics as a more rigorous extension of these assumptions.

Oxygen Ion and Vacancy Diffusion
The rate of diffusion of oxygen ion and oxygen vacancy species is described as an Arrhenius Eq. 6. Here, the important parameters are the thermal barrier, ionic charge, and jump distance for each diffusing species.

Frenkel Pair Generation
Oxygen vacancy formation is achieved through the production of Frenkel pairs, requiring the breaking of metal-oxygen bonds. A thermochemical description of dielectric breakdown exists (McPherson et al., 2003), in which the activation energy for breakdown is lowered by the local electric field projection along a polarizable bond axis. This expression is common in describing dielectric breakdown in thin insulators and is common in oxidereliability studies Eq. 7.
Within this model, breakdown is expected to begin at an electric field that lowers the effective thermal barrier to zero.
The value of E a,g determines the minimum voltage needed for forming, and is therefore an important consideration for the design of memristor circuits, as previously discussed. Using typical values for HfO 2 (p 0 11 × 10 −10 Cm, E a,g 4.5 eV, ϵ r 21), the breakdown field |E bd | 5.3 MV cm −1 . This value corresponds to a nominal forming voltage roughly equal to half the HfO 2 thickness when measured in nanometers (i.e., V form ≈ 2.5 V for a 5 nm film). Empirically, E a,g is found to reduce with increased volume fraction of oxygen vacancies, as is commonly observed in highly defective HfO x films, providing an empirical means for reducing forming voltage through the controlled introduction of defects. For example, choice of precursor (Hazra et al., 2019) and reaction time (Hazra et al., 2020) for the atomic-layer deposition of HfO x films have a profound influence on forming voltage and pre-forming high-resistance levels. Furthermore, recent work have incorporated a model having a large (68%) reduction in the activation energy of forming oxygen vacancies at the Ti/HfO 2 interface as opposed to the bulk (Xu et al., 2020). This may be anticipated, since the net energy cost of breaking Hf-O bonds is lowered by the large driving force of oxidation (gettering) in the Ti (Stout and Gibbons, 1955). These factors imply that E a,g is a spatially-dependent parameter, essential to filament formation dynamics. According to ab-initio work, which also includes the effect of electron injection, we assume a 36% reduction in the activation barrier to Frenkel pair generation in the vicinity of existing nearest-neighbor Frenkel defect pairs (Bradley et al., 2015). Phenomenologically speaking, this accounts for the accelerating effect that point defects have on dielectric breakdown, as is well-known from oxide reliability studies. Additionally, this attempts to account for the formation of stable clusters of oxygen vacancies in regions where the binding energy is high and may be a potential source of retention failure in addition to existing theories based on oxygen diffusion (Raghavan et al., 2015;Kumar et al., 2017).

Frenkel Pair Recombination
As previously discussed, upon formation, charged Frenkel pairs in HfO 2 are unstable, requiring additional electrons from the conduction band to neutralize the oxygen vacancy and prevent rapid recombination. It is therefore expected that the thermal barrier to recombination is small relative to other processes, producing a rapid recombination rate, which we model using a simple field-independent Arrhenius Eq. 9. We use a value of 0.2 eV, according to previous work (Larcher et al., 2012), though this value becomes most relevant during reset operations, when the recombination rate of Frenkel pairs and electron emission (filament precipitation) becomes comparable for deep level vacancy states.

Simulation Domain and Defect Modeling
Device geometry is defined in Synopsys TCAD Sentaurus with mesh refined using a maximum element size of 1 Å. Concentration profiles for each species (i.e., oxygen ions, unoccupied/occupied vacancies) were defined as point defects having Gaussian shape with decay length of 3 Å, corresponding to the minimum ion jump distance and grid spacing in our KMC model. Positions for each species were obtained from the output of the KMC simulation at each voltage step. Oxygen ions are modeled as negative fixed charges, with concentration as a parameter chosen to compensate the charge density of oxygen vacancies. Unoccupied oxygen vacancies are modeled as donors located 3 eV below the conduction band edge. Occupied oxygen vacancies are modeled as acceptors located 0.2-0.9 eV below the donor level, depending on the binding energy parameter (ΔE t ). As mentioned previously, the energy level difference between unoccupied/occupied oxygen vacancies reflects the increase in binding energy upon electron capture due to the large lattice coupling of ionized vacancies in HfO x . To model effects due to disorder, the energy levels of oxygen vacancies were defined having Gaussian energy broadening (σ 0.33 eV) consistent with similar approaches (Jiménez-Molinos et al., 2002).

Electric Current
Electric current was calculated using an electron barriertunneling model that couples each trap to the conduction band of the top and bottom electrodes through nonlocal, multiphonon-assisted inelastic and elastic transitions. Steadystate conditions were assumed. The rate of inelastic electron capture is described in terms of a maximum transition rate multiplied by the WKB tunneling probability T i,j | Ψ(y i ) Ψ(yj) | 2 and the phonon transition probability M i,j for a transition between two states denoted i and j located at y i and y j .
Sentaurus uses the asymptotic (large order) approximation to the conventional expression for I l (z), the modified Bessel function of order l contained within M i,j , and is therefore appropriate when the number of phonons emitted during a transition is large (Schenk and Heiser, 1997). Under this approximation, the capture rate for an electron in the conduction band of an electrode at y 0 to a trap located at y t can be written as: The emission rate is then computed from the capture rate using the principle of detailed balance.
Parameters for the model were defined as follows. The Huang-Rhys factor, S, was set to 17, the phonon energy to 0.07 eV, and the electron effective mass of HfO x was set to the band mass 0.1 according to similar reports (Guan et al., 2012a). Nonlocal tunneling paths were considered, extending outwards from each electrode towards the opposite electrode along the oxide thickness. Electrodes were treated as Ohmic, with a variable work function with nominal values defined to mimic realistic device structures having a Titanium top-electrode (Ψ TE 4.33 eV, m * n,te 3.2) and titanium nitride bottom-electrode (Ψ BE 4.55 eV, m * n,be 2.0) (Lima et al., 2012). Trap volumes were estimated according to the effective mass of electrons in HfO x and trap level relative to the conduction band edge (Palma et al., 1997;Jiménez-Molinos et al., 2002).

HSPICE Transient Simulations
To better utilize the insights provided by the device model, we couple the KMC + TCAD framework with a phenomenological compact model (Verilog-A) (Amer et al., 2017c) for the memristor to facilitate circuit-level simulations (in HSPICE). This compact model assumes a piecewise linear currentvoltage (I-V) behavior in pre-forming and post-forming states of the memristor. The compact model considers the resistance behavior of the memristor (during the forming process) as follows: Here, V M and R M are the voltage and resistance of the memristor (respectively). Clearly, forming voltage (V forming ), pre-forming (R Pre−forming ) and post-forming (R Post−forming ) resistance levels are the three necessary parameters for the forming operation of the memristor. This piecewise linear compact model considers ohmic Frontiers in Nanotechnology | www.frontiersin.org October 2021 | Volume 3 | Article 734121 current-voltage relation before and after forming. We extract the parameters for this compact model using extensive analysis powered by the KMC + TCAD framework. We simulate the DC I-V characteristics for 25 device instances with an idealized ramp voltage across the memristor and formulate distributions of device resistance and forming voltage. We simulate the forming operation for a memristor-based synapse circuit (shown in Figure 4A) in HSPICE to obtain the time dynamics of the voltage across the memristor. This synapse circuit can control all the operations of a memristor such as forming, set, reset and read. In this work, we only utilize the forming portion of the circuit. Therefore, the connections corresponding to other three operations are greyed out. Here, M p1 transistor controls the forming of the memristor and M n1 transistor is used to set the compliance current limit during the forming process.
To simulate the forming process in HSPICE, we calibrate the compact model with the relevant parameters (forming voltage, preforming HRS, and post-forming LRS), extracted from the KMC + TCAD simulations. Note, the ideal approach to calibrate the compact model would be to use iterations between the circuit level simulations and KMC + TCAD simulations. The approach would include running the KMC + TCAD simulations with the results from the circuit level simulation and then again simulating the circuit with the KMC + TCAD results. But, this would require a compact model that could capture the non-linear behavior observed in the KMC + TCAD simulations. Here, to calibrate the model for circuit simulations, we choose the mean values of the pre-forming HRS (1.2 MΩ), and post-forming LRS (1.3 kΩ) to set up the compact model. As for the forming voltage, we choose the maximum value of the corresponding distribution (2.4 V), to capture the worst-case scenario. For the transistors, we use the DGXFET NMOS/PMOS models for the IBM 65 nm 10LPe process.
Figures 4B-D illustrate the simulated transient characteristics of the synapse circuit (only the forming process), bearing the signature of the circuit-level interactions of the memristor. We govern the forming process with appropriately designed control signals ( Forming and Set). We first turn on the M n1 transistor by applying an appropriate gate voltage (Set). Then, we turn on the M p1 transistor to control the forming of the memristor. The voltage across the memristor (V M ) gradually increases when the forming operation begins ( Figure 4C). The memristor takes the prime share of the supply voltage (V DD ) when the two series transistors (M p1 and M n1 ) turn ON. Subsequently, after successful forming, the resistance of the memristor drastically reduces and so does the voltage across it. This circuit topology serves as the baseline for the Monte-Carlo simulations (discussed later). Figure 5 shows the filament formation dynamics for a single device simulation at a constant voltage of 3 V applied to the top electrode and 0 V applied to the bottom electrode. For this simulation, the initial defect volume fraction was chosen to correspond with 1 oxygen vacancy within the simulation domain. New defects tend to form in the vicinity of preexisting defects due to: 1) a lower formation energy in the presence of pre-existing Frenkel pairs; and 2) a higher electric field in the vicinity of a charged filament. As time progresses, filament growth proceeds towards the top electrode, where the electric field is highest. Once the electrodes have bridged, additional filament growth occurs along the width of the top electrode. This occurs due to: 1) a higher lateral electric field, 2) the effect of oxygen gettering by the Ti top electrode, which readily removes oxygen ion interstitials, and 3) screening of the electric field seen by points closer to the bottom electrode by charged vacancies near the top electrode.

Filament Formation Dynamics
Since bond breaking is modeled as a statistical process, it is necessary to evaluate the behavior of more than one device. Thus, we further investigate forming dynamics by assessing the statistical distribution of forming times for repeated device simulations. Here, we focus on the forming time-the time required to form, which we estimate by halting the simulation once a predefined number of filament states are formed. We use a volume fraction of 0.1, which should correspond to a concentration of the order of ≈1 × 10 21 cm −3 . This is comparable to what has been observed experimentally in HfOx films having high oxygen vacancy concentrations (Hildebrandt et al., 2011), reported as high as 6 × 10 21 cm −3 .
We note that since the forming process is modeled using the thermochemical model of dielectric breakdown, the forming time is synonymous with the time-dependent-dielectric-breakdown (TDDB), and is therefore an experimental observable which is straightforward to measure and, most importantly, can be used to validate kinetic monte carlo simulation models. Previous authors have compared the effect of top and bottom electrodes on the  forming time for nearly stoichiometric 10 nm thick HfO 2 thin films deposited by atomic layer deposition (Lorenzi et al., 2013;Cagli et al., 2011). Figure 6 shows a Weibull distribution of forming times obtained from our simulation and those of Lorenzi et al., 2013 for a TiN/HfO 2 /Pt device at an electric field of 5 MV cm −1 and 6 MV cm −1 . Simulation results at both values of electric field show reasonable quantitative agreement to experiment. At the lower electric field of 5 MV cm −1 , there appears to be a deviation from a Weibull distribution, however the forming times are similar in magnitude. These results illustrate the similarity between forming time and TDDB, and provide a potential route towards empirical validation of simulation models-from which the activation energy of Frenkel pair generation, E a,g and the fieldacceleration factor, c can be derived (McPherson and Mogul, 1998).
This is especially important to establish agreement to experimental results, since HfOx films can exhibit mixed crystalline phases and variable oxygen content depending on deposition conditions-both of which are expected to modify E a,g and c. Figure 7 shows the complete current-voltage (IV) characteristics for a forming, reset and set programming cycle. In order to obtain IV characteristics, snapshots of the state of the system are taken at the end of each voltage step shown in Figure 3B. In this case, the time step is 1 µs and the voltage step is 0.1 V. In Figure 7A, intermediate filament states are shown at different stages of programming. It can be seen that, at forming, a large volume fraction of the device consists of oxygen vacancy filament with a structure that extends laterally near the top electrode. The spatial extent of the filament (i.e., volume fraction of vacancies) will ultimately be controlled by the flux linkage (forming voltage × time), which is set by the compliance current in practice.

Current-Voltage Characteristics
Several of the well-known aspects of oxide memristors are captured by this simulated result shown in Figure 7B: 1) the highresistance initial state of the as-prepared thin-film; 2) A low resistance state after forming; 3) A gradual reset behavior with programmable analog high-resistance levels; and 4) A lowresistance state following set of the order of kilo-Ohms. The ability to quantify and visualize the increase in current upon set (potentiation) and the decrease in current upon reset (depression) is a key benefit of Figure 7B incorporating TCAD Sentaurus for modeling synaptic behavior. Figure 8 shows key aspects of the filament electrostatics. In Figure 8A, the electrostatic potential and x and y components of the electric field throughout the device are shown. In regions where the filament bridges the top and bottom electrode, the voltage drop across this region is large enough such that the filament region has a net negative potential near the bottom electrode. This agrees well with experimental results which relied on in-situ TEM electron holography measurements (Li et al., 2017), in which they described the filament as a "negative potential synapse." Our results show that this stems directly from the negative space-charge associated with the filament, assumed to be due to electron capture (V 2+ O → V 2− O ), and supported by both experimental and theoretical insights. We show this explicitly, by comparing two different line plots-outside the filament and within the filament-in Figures 8B,C. A dashed line at a potential of zero is added as a visual aid, clearly indicating a negative potential within the filament. This is also reflected in the energy band diagram in Figure 8C, which shows a negative curvature as expected for a negative space charge.

Monte Carlo Analysis of Synapse Forming Circuit
Finally, we use the unique capability of the KMC + TCAD model to investigate a device-circuit co-design strategy. We test memristor characteristics for different levels of initial oxygen vacancy volume faction (x), a design variable that can be easily controlled during the fabrication process. Figure 9A shows the pre-forming HRS for different values of x for 25 devices (each) obtained from the KMC + TCAD simulations. Considering the circuit-level scenario of a synapse, the variations in the characteristics of multiple transistors need to be superposed with the inherent device-level variations of the memristor. To account for all these variations concurrently, we run 1000-point (3σ) Monte-Carlo simulations for the forming circuit using the data obtained from the KMC + TCAD framework. We utilize the dependence of the pre-forming HRS on the initial oxygen vacancy volume faction and the threshold voltage variation of the PMOS and NMOS transistors to set the input distributions for the Monte-Carlo simulations. Figure 9A shows the pre-forming HRS for different values of x for 25 devices (each) obtained from the KMC + TCAD simulations. Without any loss of generality, we run the Monte-Carlo simulation for two values of x (0.02 and 0.04%). To incorporate the threshold voltage variation of the transistors, we use a gaussian distribution with a mean value equal to the nominal threshold voltage (0.65 V for 65 nm DGXFET transistors) and standard deviation of 20 mV (shown in the table of Figure 9B). We also run the Monte-Carlo simulations with different levels of current compliance, controlled by applying appropriate gate voltage (Set) to M n1 . Figures 9C,D show the scatter plots for the average current through the memristor and average power of the forming circuit, respectively. Each of these metrics have been reported for different levels of compliance currents (different levels of Set). To ensure a fair comparison, we allow a constant time for the forming process for all cases. Naturally, we observe that for lower compliance limits, many instances of the Monte-Carlo simulations exhibit "unsuccessful" forming. Note, a compliance limit may lead to a different level of post-forming LRS and hence might be treated as successful, if the post-forming LRS is known during the design stage. However, if such changes in the post-forming LRS occurs dynamically and randomly, those will lead to read/sense failure. Therefore, we simplify our analysis by tagging such cases as 'Forming Failure'. Higher compliance limit allows most memristor instances to successfully form and hence leads to a larger average current level ( Figure 9C). If the memristor can successfully form, it goes to post-forming LRS (1.3kω in our simulation). Otherwise, it remains in the preforming HRS which is much larger compared to LRS. Therefore, the increase in the value of Set increases the average current of the memristor due to the increase in the FIGURE 8 | (A) Filament electrostatics after forming, including the electrostatic potential, and electric field components. (B) Comparison of electrostatic potential inside and outside of a filament region as indicated in the red lines in (A). The filament region shows a region of negative potential and can be likened to a "negative potential synapse." (C) Comparison of the energy band diagram inside and outside of the filament. Within the filament, a negative curvature indicates that the filament is a negative space-charge region.
Frontiers in Nanotechnology | www.frontiersin.org October 2021 | Volume 3 | Article 734121 number of formed memristors. Since, the average power of the forming circuit is very closely related to the memristor current, the average forming power shows the same trend like the average memristor current ( Figure 9D). The initial oxygen vacancy volume factions lead to similar results, with different levels of mean value and standard deviation for the memristor current and forming power ( Figures 9C,D). Figure 9A shows that the preforming HRS for x 0.02% has a larger standard deviation compared to that for x 0.04%. Therefore, the memristor current and forming power obtained from the Monte-Carlo simulation show larger standard deviation for x 0.02% compared to the case of x 0.04%. But only for Set 1.5 V, the smallest value of pre-forming HRS of the memristors that cannot form becomes comparable to the value of the post-forming LRS. Therefore, the effect of x on the standard deviation of the memristor current and forming power gets suppressed. Figures 9E,F show the histogram plot for Monte-Carlo results of memristor current and forming power shown in the scatter plots ( Figures 9C,D).
Based on these Monte-Carlo simulations, we correlate the compliance limits (controlled by the Set pulse) with the pre- forming HRS (R Pre−forming ) level of the memristor. Figure 9G illustrates the combinations of Set and R Pre−forming that lead to successful forming (and vice versa). Clearly, for a given compliance limit, the pre-forming HRS level of a memristor needs to be higher than a critical threshold (R form−TH ), illustrated (in Figure 9G) as a line separating the successful and unsuccessful forming cases. Figure 9H shows similar trends for a different initial oxygen vacancy volume faction. Our analysis shows a pathway to optimize the synapse circuit by correlating the material and circuit-level design knobs.

CONCLUSION
The ability to design and implement fast, scalable and robust neuromorphic systems relies heavily upon our fundamental understanding of memristor switching. Oxide memristors, envisioned for RRAM-based neuromorphic systems, exhibit changes in resistance state through multiple synergistic effects involving electronic and atomic degrees of freedom, often modelled as separate influences. One of the main purposes of this work was to establish a more direct connection between the two in order to: 1) provide a unified view of filament evolution and electronic conduction; 2) to implement this description within a state-of-the-art TCAD framework for modeling electric conduction; and 3) gain circuit-level insight.
Here, we have argued the use of a simple model of filament evolution that makes explicit use of Fermi-Dirac statistics, coupling the rate of defect generation and recombination to electronic transitions associated with conduction and lattice relaxation. By combining Synopsys TCAD Sentaurus with Kinetic Monte Carlo simulations of filament evolution, we have shown the ability to quantify both the common and subtle aspects of resistive-switching behavior of HfO x memristors. Quasi-static snapshots of the device state-consisting of positive/negative oxygen vacancies, and oxygen ions-were taken at various voltages to obtain IV characteristics under stepped voltage ramp conditions. Electric conduction in oxygen vacancy filaments is modeled as trap-to-band transitions between occupied and unoccupied electronic states assisted by multiphonon absorption and emission. According to Fermi-Dirac statistics, such processes are expected to occur within a band of energies in the vicinity of the Fermi level wherein both occupied and unoccupied states are probable. Thus, the occupancy of a trap and its relation to the Fermi level is fundamentally related to transition rates associated with electronic conduction. The use of TCAD Sentaurus provides a powerful framework for modeling these and other conduction processes as well as visualizing filament electrostatics, as we've shown. In particular, we have obtained results that are consistent with experimental observations of a negative space charge and potential associated with a vacancy-rich filament. Our approach will enable the more efficient evaluation of memristor device behavior and circuit performance, stemming from physics-based modeling, having a direct impact and benefit on the fields of neuromorphic computing, memory design and dynamical systems.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
AZ conceived of and designed the Kinetic Monte Carlo algorithm and TCAD Sentaurus modeling framework. SA performed the transient simulations, and the circuit-level Monte Carlo analysis. ZW assisted with device-level simulations and parameter influences. RW designed the synapse circuit. AA analyzed the circuit simulations and Monte Carlo data. AZ AA, and GR jointly supervised the work and analyzed the results. All authors contributed to writing the article and approving the submitted version.

FUNDING
This material is based in part on research sponsored by Air Force Research Laboratory under agreement number FA8750-19-1-0025. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of Air Force Research Laboratory or the U.S. Government.