We present a unique compact model for oxide memristors based upon the concentration of oxygen vacancies as state variables. In this model, the increase (decrease) in oxygen vacancy concentration is similar in effect to the reduction (expansion) of the tunnel gap used as a state variable in existing compact models, providing a mechanism for the electronic current to increase (decrease) based upon the polarity of the applied voltage. Rate equations defining the dynamics of state variables are obtained from simplifications of a recent paper in which electronic processes (i.e., electron capture/emission) were combined with atomic processes (i.e., Frenkel-pair generation/recombination, diffusion) stemming from the thermochemical model of dielectric breakdown. Central to the proposed model is the effect of the electron occupancy of oxygen vacancy traps on resistive switching dynamics. The electronic current is calculated considering Ohmic, band-to-band, and bound-to-band contributions. The model includes uniform self-heating with Joule heating and conductive loss terms. The model is calibrated using experimental current–voltage characteristics for $HfO2$ memristors with different electrode materials. Though a general model is presented, a delta-shaped density of states profile for oxygen vacancies is found capable of accurately representing experimental data while providing a minimal description of bound-to-band transitions. The model is implemented in Verilog-A and tested using read/write operations in a $4\xd74$ 1T1R nonvolatile memory array to evaluate its ability to perform circuit simulations of practical interest. A particular benefit is that the model does not make strong assumptions regarding filament geometry of which scant experimental-evidence exists to support.

## I. INTRODUCTION

Several physics-inspired compact models have been proposed for circuit simulation of oxide memristors^{1–8}—key elements for emerging applications in non-volatile memories, unconventional computing, and biologically inspired computing alike.^{9–13} These are largely based on assumptions regarding the geometry of a “conductive filament,” described either in terms of a one-dimensional tunnel gap^{14–18} or the size/shape of the filament approximated as cylindrical^{19} or rectangular.^{20,21} Electronic conduction processes—derived from changes in filament geometry—are usually modeled with drift/diffusion,^{22–26} hopping,^{27} trap-assisted and band-to-band tunneling,^{17} Poole–Frenkel emission, site-percolation,^{28–31} or interfacial redox reactions.^{19} Additionally, conduction models sometimes incorporate simple thermal models (e.g., uniform Joule heating) since local temperature is purported to have a significant impact on electrical characteristics, particularly retention^{32} and multibit operation.^{33}

The scope of existing tunnel-gap models can be understood simply by taking any or all of the nominal conduction processes present in metal–insulator–metal devices (e.g., direct and Fowler–Nordheim tunneling, thermionic emission, Poole–Frenkel emission, Ohmic conduction, ionic conduction, and space-charge limited conduction^{34}) and replacing the insulator thickness parameter with a dynamic variable (i.e., the tunnel gap) ranging from zero to a few nanometers.^{17} The dynamics of the tunnel gap are based on the diffusion kinetics of oxygen ions presumed to be rate-limiting—approximated in one dimension along the oxide thickness using the Mott–Gurney drift velocity.^{35} A positive voltage drives ions toward the top electrode, reducing the tunnel gap; a negative voltage drives ions toward the bottom electrode, increasing the tunnel gap. Since this implicitly assumes Frenkel pairs recombine (generate) in the presence (absence) of oxygen ions, the thermal barriers for recombination/generation must be smaller than that of oxygen ion diffusion wherever the model is applied so as not to become rate-limiting. The zero-field thermal barrier for Frenkel-pair generation can be 1.4 eV (at grain boundary sites^{36}), 2.2 eV (at gettering electrode surfaces^{37}), and several electronvolts in the bulk (3.6 ,^{36} 6.8 eV^{37}), all of which are much larger than that of oxygen ion diffusion ($\u2248$0.7 eV^{32}). Consequently, tunnel-gap models are adequate for describing reset or set operations,^{17,35,38} particularly when the top electrode can be regarded as an oxygen reservoir with permeable interface and the electric field is large enough (tunnel gap small enough) that the thermal barrier for Frenkel-pair generation is reduced and no longer rate-limiting. A more complex or entirely new framework would be required to incorporate a physical description of forming self-consistently within this model.^{39} Given that oxide memristors are flux driven,^{40} a description that takes into account the entire electrical history of the device (forming, reset, and set) is critical from a self-consistent physical modeling perspective.

It is important to note that all compact models are phenomenological, involving many simplifications to enable circuit simulation at various levels of complexity and computational efficiency. Ultimately, the accuracy and speed of the device compact model are selected based on the needs of the circuit designer rather than the physical accuracy alone. For example, simple piecewise linear models may be desired when the focus is on developing larger neuromorphic systems,^{41} whereas predictive physics-based models may be desired for testing smaller memory array architectures.^{42} Ideally, memristor compact models would stem from a common physical framework in which complexity can be selectively refined at the circuit level, similar to the different levels of the well-known SPICE models for conventional bulk MOSFETs (square law, bulk charge, etc.).

In this work, a compact model for oxide memristors is presented that is based on the concentration of oxygen vacancies and their electron occupancy as opposed to one-dimensional filament geometry. The shift from a tunnel gap to a concentration state variable may be more physically appropriate, given that the formation of a one-dimensional defect chain is statistically unlikely considering other competing factors, such as vacancy diffusion, thermophoresis,^{32} and the increased thermodynamic stability of defect clusters,^{43} expected to lead to additional conduction pathways. These, along with the presence of grain-boundaries and other imperfections, which may reduce the zero-field formation enthalpy of oxygen vacancies, suggest a variety of filament shapes with varying connectivity as seen in comprehensive modeling approaches^{39,44,45} and experimental observations.^{46}

In our model, the oxygen vacancy concentration is split into two state variables to allow the electron occupancy of the vacancies to vary between occupied and unoccupied depending on electron capture and emission processes. This is based on recent KMC modeling approaches^{45} and reflects *ab initio* calculations showing increased thermodynamic stability of negatively charged oxygen vacancies under conditions of electron-injection^{43} as well as experimental observations of a negative space charge associated with suspected filament regions.^{46} This has immediate implications on retention since a larger thermal emission barrier would result in a more stable low-resistance state and, therefore, longer retention time as we show by comparison to experimental data of memristors having different electrode materials. The compact model is implemented in Verilog-A, using simplifying assumptions made to our previous model.^{45} In particular, a key approximation is that the field-dependence of the microstate transitions is approximated using the average electric field as opposed to the local electric field. Since electron capture and emission processes are taken into account, the model utilizes additional parameters associated with the electrodes (e.g., work function, effective mass) and those of the oxygen vacancy defects (e.g., trap energy level/ionization energy, capture cross section, and thermal capture barrier) that are intimately connected to the resistive switching dynamics. The model also considers temperature effects and, most importantly, does not assume a particular filament shape, gap size between filament tip, or any other restrictive geometrical constraints that are not firmly supported by experimental data. Each of these features makes this compact model attractive due to its predictive modeling capability and pairing ability with ongoing experimental investigations to identify and observe filament geometries.^{32,46} Being generally based upon a concentration of defects as opposed to a particular filament geometry, the model is open to different interpretations of conduction (e.g., drift/diffusion, trap-assisted-tunneling, or percolation) as we show.

Last, the treatment of oxygen vacancies as an electron trap allows the material-specific density of states to be specified to provide increased accuracy/reduced efficiency (continuous energy dependence) or decreased accuracy/increased efficiency (discrete trap level) based on the particular needs of the circuit designer. The density of states is more intimately connected to materials chemistry than is a tunneling gap, allowing greater flexibility when comparing model predictions to experimental data (e.g., defect spectroscopy^{47}). For circuit-simulation, a delta density of states profile is least mathematically complex to implement, requiring no numerical integration at each time step. Thus, we extensively evaluate the compact model (implemented in Verilog-A) based on a delta-profile density of states. The model is validated using independent experimental datasets having different electrode materials and tested using a simple $4\xd74$ 1T1R array representative of a small nonvolatile memory array to evaluate the model’s capability to perform basic circuit simulations.

## II. COMPACT-MODEL FRAMEWORK

### A. General assumptions and notation

In Table I, we define the notation for symbols used throughout the model as well as nominal material parameters for $HfO2$, which we evaluate herein as the oxide material. In Fig. 1, the device structure and definitions of energy levels are depicted. Although the model is presented using parameters for $HfO2$, it is generally applicable to transition metal oxides of interest for memristors since it is based upon a general thermochemical model of dielectric breakdown,^{48,49} common to many insulating oxides and/or polar solids in general.

Parameter . | Description . | Value and unit . |
---|---|---|

q | Elementary charge | 1.602 176 × 10^{−19} C |

k | Boltzmann constant | 1.380 649 × 10^{−23} J K^{−1} |

$\u210f$ | Reduced Planck constant | 1.054 571 × 10^{−34} J s |

T_{0} | Ambient temperature | 300 K |

ɛ_{r} | Relative permittivity | 30 |

R_{0} | Attempt-to-escape frequency | 10 THz |

$Ea,gen0$ | Activation energy for Frenkel-pair generation | 7.3 eV |

$Ea,rec0$ | Activation energy for Frenkel-pair recombination | 0.15 eV |

E_{ion} | Ionization energy of oxygen vacancy defects | 2.0−3.0 eV |

σ_{0} | Capture cross section of oxygen vacancies | 1 × 10^{−14} cm^{2} |

Φ^{TE} | TE work function ($Ti$) | 4.33 eV |

Φ^{BE} | BE work function ($TiN$) | 4.50 eV |

$meTE$ | TE electron effective mass ($Ti$) | 3.2 |

$meBE$ | BE electron effective mass ($TiN$) | 2.0 |

$neTE$ | TE electron concentration ($Ti$) | 5.67 × 10^{22} cm^{−3} |

$neBE$ | BE electron concentration ($TiN$) | 2.88 × 10^{22} cm^{−3} |

$meox$ | Oxide electron effective mass ($HfO2$) | 0.1 |

$p\u2192$ | Oxide molecular dipole moment ($HfO2$) | 15 × 10^{−10} C m |

E_{g} | Oxide bandgap energy ($HfO2$) | 5.9 eV |

χ | Oxide electron affinity ($HfO2$) | 2 eV |

μ_{eff} | Effective band mobility ($HfO2$) | 1 cm^{2} V^{−1} s^{−1} |

c_{p,ox} | Oxide specific heat capacity ($HfO2$) | 120 J kg^{−1} K^{−1} |

C_{add.} | Addenda heat capacity (device) | 1.45 × 10^{−9} J K^{−1} |

ρ_{ox} | Oxide mass density ($HfO2$) | 9800 kg m^{−3} |

κ_{ox} | Oxide thermal conductivity ($HfO2$) | 1 W m^{−1} K^{−1} |

N_{sites} | Site density ($HfO2$) | 4.38 × 10^{19} cm^{−3} |

y_{t} | Trap center ($HfO2$) | 0.5t_{ox} |

Parameter . | Description . | Value and unit . |
---|---|---|

q | Elementary charge | 1.602 176 × 10^{−19} C |

k | Boltzmann constant | 1.380 649 × 10^{−23} J K^{−1} |

$\u210f$ | Reduced Planck constant | 1.054 571 × 10^{−34} J s |

T_{0} | Ambient temperature | 300 K |

ɛ_{r} | Relative permittivity | 30 |

R_{0} | Attempt-to-escape frequency | 10 THz |

$Ea,gen0$ | Activation energy for Frenkel-pair generation | 7.3 eV |

$Ea,rec0$ | Activation energy for Frenkel-pair recombination | 0.15 eV |

E_{ion} | Ionization energy of oxygen vacancy defects | 2.0−3.0 eV |

σ_{0} | Capture cross section of oxygen vacancies | 1 × 10^{−14} cm^{2} |

Φ^{TE} | TE work function ($Ti$) | 4.33 eV |

Φ^{BE} | BE work function ($TiN$) | 4.50 eV |

$meTE$ | TE electron effective mass ($Ti$) | 3.2 |

$meBE$ | BE electron effective mass ($TiN$) | 2.0 |

$neTE$ | TE electron concentration ($Ti$) | 5.67 × 10^{22} cm^{−3} |

$neBE$ | BE electron concentration ($TiN$) | 2.88 × 10^{22} cm^{−3} |

$meox$ | Oxide electron effective mass ($HfO2$) | 0.1 |

$p\u2192$ | Oxide molecular dipole moment ($HfO2$) | 15 × 10^{−10} C m |

E_{g} | Oxide bandgap energy ($HfO2$) | 5.9 eV |

χ | Oxide electron affinity ($HfO2$) | 2 eV |

μ_{eff} | Effective band mobility ($HfO2$) | 1 cm^{2} V^{−1} s^{−1} |

c_{p,ox} | Oxide specific heat capacity ($HfO2$) | 120 J kg^{−1} K^{−1} |

C_{add.} | Addenda heat capacity (device) | 1.45 × 10^{−9} J K^{−1} |

ρ_{ox} | Oxide mass density ($HfO2$) | 9800 kg m^{−3} |

κ_{ox} | Oxide thermal conductivity ($HfO2$) | 1 W m^{−1} K^{−1} |

N_{sites} | Site density ($HfO2$) | 4.38 × 10^{19} cm^{−3} |

y_{t} | Trap center ($HfO2$) | 0.5t_{ox} |

Below, we summarize the main assumptions used in deriving the compact model:

- $\u2219$
The oxide is single phase (i.e., no polymorphism) with material properties of the bulk (e.g., effective mass, permittivity, bandgap, electron affinity). In practice, this implies that the material properties are homogeneous and is not necessarily tied to the bulk crystalline or amorphous state.

- $\u2219$
The local electric field can be approximated with the spatially averaged field.

- $\u2219$
Oxygen vacancies behave as electron traps (positive when empty $VO2+$, negative when occupied $VO2\u2212$, accommodating 4 electrons

^{43}), defined in terms of a continuous electronic density of states to specify the ionization energies/defect levels participating in electron capture/emission. - $\u2219$
The interface tunneling transmission probabilities for the top and bottom electrodes are approximately equal.

- $\u2219$
Diffusion of ionic point defects is not considered.

- $\u2219$
Trap-to-trap coupling is not considered.

### B. State variable dynamics

The state of the oxide is represented using three species: unoccupied oxygen vacancies $NVO2+$, occupied oxygen vacancies $NVO2\u2212$, and empty states $N0$; each of these is taken to be volume-averaged quantities (i.e., concentrations). These are depicted in Fig. 1(a). At any point in time, the sum of these is constant and equal to the total density of accessible sites, which is nominally equal to the atomic density of $HfO2$ ($Nsites\u22482.77\xd71022cm\u22123$),

We note that, in general, $Nsites$ can be regarded as a parameter that establishes an upper bound for the state variables $NVO2+$, $NVO2\u2212$, and $N0$. When $Nsites$ is equal to the atomic density, the intrinsic state of the memristor is reachable in the absence of external current-compliance. As will be shown, reducing $Nsites$ below the atomic density has a similar effect as forcing a current-compliance. A value of 4.38 × 10^{19} cm^{−3} was used throughout this work based on agreement to experimental data.

From (1), since $Nsites$ is a constant, the net rate of change in each species must sum to zero,

These rates can be expressed in terms of the following physical transitions: Frenkel-pair generation $Rgen$ and recombination $Rrec$, electron capture $Rc$, and emission $Re$ in which capture/emission of electrons can take place between either the top electrode (TE) or the bottom electrode (BE).

The rate of change in electron occupancy of oxygen vacancy traps due to electron capture and emission is given by Eq. (3),

In terms of $ft$ and the density of states $g(Et)$, the occupied, unoccupied, and total vacancy concentration can be calculated using Fermi–Dirac statistics,

Trap-to-trap coupling is neglected, which would otherwise give rise to a separate Fermi level for trap states $Eft$ and corresponding Fermi function $ft$ for use in calculating trap concentrations. Instead, the number density of traps is used to compute the electron-occupancy probability,

Using these definitions, integration of Eq. (3) over the density of trap states $g(Et)$ provides the rate equation for the occupied vacancy concentration $NVO2\u2212$,

where $g^(Et)$ is the normalized density of states profile for the oxygen vacancy levels,

Equation (9) simply states that unoccupied vacancy states $NVO2+$ are filled by either electrode (TE/BE), and occupied states $NVO2\u2212$ are emptied to either electrode (TE/BE). The energy integral allows for the density of states to take on any energy dependence (e.g., delta, uniform, Gaussian).

To account for defect creation, which initially gives rise to the formation of Frenkel pairs consisting of oxygen ions and unoccupied oxygen vacancies (i.e., $O2\u2212$, $VO2+$), atomic processes are incorporated by assuming that the empty states in our model $N0$ are increased (decreased) due to recombination (generation) of Frenkel pairs,

Using Eq. (2), along with Eqs. (9) and (11), the remaining rate equation for the unoccupied vacancies is given by Eq. (12),

Together, Eqs. (9), (11), and (12) define a dynamic relationship between *atomic* and *electronic* processes, which we have evaluated numerically in a recent paper in two dimensions, taking into account diffusion, Frenkel-pair generation and recombination, and electron capture and emission.^{45} For simplicity, the compact model developed here does not consider diffusion; diffusion of charged point defects would produce a spatially and time-dependent electric field, which would complicate compact-model development.

### C. Transition rates

The rates $Rgen$ and $Rrec$ are atomic processes describing the creation and removal of Frenkel pairs of oxygen ions and vacancies (i.e., $O2\u2212$ and $VO2+$), whereas the rates $Rc$ and $Re$ are electronic processes describing the capture and emission of electrons by oxygen vacancies to either electrode.

The rates of Frenkel-pair generation/recombination can be described according to classical transition state theory, commonly applied in kinetic Monte Carlo simulations of oxide memristors,

The thermal barriers for these transitions are electric field dependent (the local field can promote charge separation or recombination), which we model as symmetric raising/lowering of the thermal barrier with field $E$,

These equations are implemented such that the maximum rate cannot exceed the exponential prefactor $R0$, occurring for high positive/negative voltages, otherwise leading to negative thermal barriers.

The capture rates $RcTE/BE$ can be described using a classical or quantum approach. The latter uses Fermi’s golden rule for the rate of interaction $Rinter.$ between the bound state (within the oxide) and the band state (within the electrode). For simplicity, the interaction potential can be confined to a box-shaped volume $Vt$ in which the interaction potential energy is set to the defect level referenced to the conduction band $Et=Ec\u2212Eion$. For 2D or 3D modeling, in which spatially dependent features are integral, a quantum approach may be desired for its improved accuracy. A complete derivation of both is provided in the supplementary material of our previous work^{45} and elsewhere^{50,51} as implemented in standard device simulators.^{52}

Equation (17) shows the resulting quantum model for an elastic transition,

An additional factor is included to account for the tunneling probability $Ptun.$, using the Wentzel–Kramers–Brillouin (WKB) approximation for a triangular barrier of uniform field.^{34} Here, $f(E)TE/BE$ is the Fermi–Dirac function, defined separately for each electrode; $yt$ is the position of the trap relative to the electrode from which an electron is captured; $meTE/BE$ is the electron effective mass of the electrode; $Vt$ is the trap volume; and $y0TE/BE$ is a WKB tunneling parameter.^{34}^{,}$Vt$ and $y0$ can be expressed using a box-shaped potential for the trap state and a triangular barrier (uniform field) approximation for the conduction band as^{50,51}

where $\Delta ETE/BE=\Phi TE/BE\u2212\chi $ is the conduction band offset at each electrode interface, $\Phi TE/BE$ is the electrode work function, and $\chi $ is the electron affinity of the oxide; $y0$ has two values, one for each electrode.

Equation (17) assumes that the interaction between bound states and band states is elastic and that the electrodes can be approximated using a bulk parabolic band dispersion.

Here, for simplicity, we implement a classical model for the rate of interaction in terms of a capture cross section, thermal velocity, carrier concentration, and thermal barrier. The classical model has the same phenomenological aspects of the quantum model, though it differs quantitatively regarding the dependence on the trap level and applied voltage.

Equation (20) shows the resulting classical model,

where $Rc,0TE/BE$ and $EtfTE/BE$ are defined as

Regarding the $\u2213$ in front of the field-dependent terms in Eq. (20), the $\u2212$ is used for the top electrode (TE) and the $+$ is used for the bottom electrode (BE). In either model, the emission rate to either electrode is calculated from the capture rates by the detailed balance criterion,

For the classical model, this leads to the following:

Regarding the $\xb1$ in front of the field-dependent terms in Eq. (25), the $+$ is used for the top electrode (TE) and the $\u2212$ is used for the bottom electrode (BE).

The voltage dependence of these expressions enters via symmetric raising/lowering of the quasi Fermi energies from the equilibrium Fermi level (taken as zero, $EfTE,0=EfBE,0=0$) by amount $0.5qV$ such that the total difference between the Fermi energies is $EfBE\u2212EfTE=qV$,

We note that, in general, transitions will be field-dependent; a field-dependence introduces a time-dependence for the rates since the local electric field strength is dependent upon the spatial arrangement of charged defects—which varies with time. For this reason, numerical modeling approaches require the solution to Poisson’s equation at each time step as has been done elsewhere.^{44,45,53,54} Toward compact-model development, as a simplification, everywhere, the local field $E$ is replaced with the average field $E\xaf=Vtox$. Consequently, the rates depend on the applied voltage but not on time.

### D. Electronic conduction

The current density is computed from the following components: band-to-band tunneling between electrodes $Jt$, Ohmic conduction within the bands $JOhmic$, and bound-to-band tunneling from trapped electrons (i.e., trap-assisted tunneling) $Jtat$. These are shown diagrammatically in Fig. 1(b),

The band-to-band tunneling component $Jt$ accounts for Fowler–Nordheim tunneling through the thin, insulating portion of the $HfO2$,^{34}

Alternatively, this can be refined with image-force correction using the general Simmons formula^{55} as in Ref. 42.

The Ohmic contribution is included based on charge neutrality considerations and considerations of information derived from the prior literature. In our previous work, unoccupied oxygen vacancies were modeled as donors, and occupied oxygen vacancies were modeled as acceptors.^{45} The charge state of the vacancy changed each capture/emission event since dynamic updating of the charge state of the system for each discrete transition is possible using a KMC approach. Throughout these processes, charge neutrality can be satisfied using the space charge of oxygen ions and an additional compensating positive charge to balance the additional negative charge created due to the capture of electrons injected from the electrodes.

To illustrate, we consider a process leading to electron capture by an oxygen vacancy ($\u2217$ denotes a neutral site). Initially, Frenkel-pair generation leads to neutral pairs of oxygen ions $O2\u2212$ and oxygen vacancies $VO2+$,

After electron capture, negative oxygen ions remain, positive oxygen vacancies $NVO2+$ are removed, and negative oxygen vacancies $VO2\u2212$ are created—producing a net negative space charge equal to $\u22124e$ per capture event (assuming four electrons are captured^{43}),

Thus, by assuming that holes equal in concentration to $4NVO2\u2212$ are always present, bulk charge neutrality can be maintained within the oxide,

This assumes that the oxygen ions $O2\u2212$ are present throughout all processes behaving as negative fixed charges. However, oxygen ions can be removed or introduced into the oxide via gettering processes at the electrode interfaces, establishing a spatial dependence, and our model does not incorporate the oxygen ion concentration as a state variable. Instead of introducing a new state variable, we assume that an electron concentration equal to $NVO2+$ and a hole concentration equal to $NVO2\u2212$ are always present within the bands to ensure charge neutrality. We note that, since the oxygen vacancy state is a deep level state, the binding energy is stronger and the electron/hole states are not expected to behave as weakly bound, fully delocalized states with their respective band mobility, especially at high defect concentrations. Instead, we assume a single effective mobility $\mu eff$ parameter, with a drift current proportional to the total vacancy concentration,

Since this is a bulk component, it tends to become more important as the oxide thickness increases (reduced tunneling probability), as one might expect for highly defective insulators (e.g., $SiNx$); for thin films ($<5nm$), the tunneling current dominates. Note, we model the Ohmic conduction as a simple phenomenological model using a constant effective mobility parameter. If desired, this can be further modified to be temperature and field-dependent using a Poole–Frenkel mobility model as in disordered media,^{34}

where $\mu 0$, $E0$, $\beta $, and $\gamma $ are fitting parameters.^{52} This is expected to be more significant for oxide memristors consisting of smaller bandgap oxides (e.g., $TiO2$) presenting a smaller thermal barrier for local bound-to-band emission. However, it should be noted that, since Poole–Frenkel conduction is an additional emission component (to the bands within the oxide), the rate equation framework for the capture/emission process would need to be self-consistently modified if included since all capture/emission processes effect the electron occupancy of the oxygen vacancies. For simplicity, we do not consider such processes.

The trap-assisted-tunneling component is computed considering the current produced within a differential volume element centered about a discrete trap located within the oxide at $(xt,yt,zt)$.^{56} In one dimension, we can write

where $R$ is the net transition rate per unit volume, defined in terms of the electronic transition rates and density of states for traps $g(Et)$, according to previous work,^{45}

Each electronic transition rate is position-dependent due to the tunneling factors. For a symmetric energy barrier at each electrode interface, the net transition rate will be maximum and sharply peaked for defects situated near the center of the oxide $yt\u22480.5tox$ (the centroid of a uniform defect profile), regardless of the defect distribution. If the integral in Eq. (35) can be approximated by its peak value, the spatial dependence can be eliminated. Replacing $yt$ with $0.5tox$ yields the following:

This approximation effectively samples a narrow spatial range of the defect concentration (near the center of the oxide), regardless of the specific spatial distribution, since these will contribute more to the total integral in Eq. (35).

### E. Self-heating

Self-heating effects are approximated by assuming uniform Joule heating where the oxide volume dissipates heat via conduction to the external device temperature at $T0$,

Starting with an initial temperature $T(0)=T0$, equal to the device temperature, the temperature at the next time step is given in terms of the instantaneous temperature $T(t)$, thickness $tox$, specific heat capacity $cp,ox$, mass density $\rho ox$, thermal conductivity $\kappa ox$ of the oxide, current density $J$, and ambient temperature $T0$,

Equation (39) is the numerical solution to the time-dependent heat equation in one dimension, considering Joule heating and conductive heat losses to the electrodes [Eq. (38)]. Factor $2$ is due to the assumption that heat dissipates from the center of the device over a distance equal to half the oxide thickness ($\Delta x=0.5tox$). The model is approximate, and, in practice, the quantity $2\kappa oxtox$ is the specific thermal conductance and is a fitting parameter. Use of self-heating creates internal temperature as an additional state variable for each memristor and, since temperature is not intrinsically bound, requires a window function to prevent temperatures from reaching absurd values. However, dynamic inclusion of self-heating may not always be necessary. For example, in oxide memristors having thin dielectrics, primary electronic conduction involves coupling traps within the dielectric to nearby electrodes via tunneling. Tunneling current, by contrast, produces a boundary condition for the current density, as opposed to a local current density within the oxide. Joule heating is, therefore, expected to occur near the electrodes—where tunneling electrons dissipate their excess kinetic energy into the contacts, which behave as thermal reservoirs. Thus, for simplicity, if the electrode can be treated as an *efficient* thermal reservoir, the dynamic effects of Joule heating might be approximated simply by observing how the model responds to static changes to the ambient temperature.

Close inspection of Eq. (39) shows that the “device” is taken as the oxide volume. Since the oxide volume constitutes a very small fraction of the thermal system (see the SEM image within Ref. 57), a more realistic inclusion of self-heating considers that the oxide presents a negligible thermal mass in comparison with the electrodes and surrounding device structure, which is much more efficient at conducting heat due to the lower thermal conductivity of the oxide. Without such a correction, inclusion of thermal models leads to absurdly high temperatures and/or nonsensical fitting parameters. For example, using the device geometry and thermal resistance ($Rth=\Delta x\kappa A$) reported in one comprehensive modeling approach of $HfO2$ memristors,^{39} the effective thermal conductivity of the $HfO2$ evaluates to 100 W m^{−1} K^{−1}, which is $\u2248100$ times larger than the measured value.^{58} Equation (39) is, therefore, modified to include a correction factor that accounts for the addenda heat capacity due to the bulkier elements in thermal contact with the oxide so that temperature effects can be modeled self-consistently using material parameters of the $HfO2$,

where $CF$ is a correction factor that takes into account the additional heat capacity of the system $Cadd$,

It is simple to show [by adding additional thermal masses to the system defined by Eq. (38)] that the correction factor is defined as the ratio of the addenda heat capacity to the heat capacity of the oxide $Cox=Atox\rho oxcp,ox$—the addenda behaving as a heat sink; in this expression, the electrode area is $A$ and the oxide volume is $Atox$,

Table I includes an estimate of the addenda heat capacity used in this work.

### F. Numerical implementation

Equations (9), (11), and (12) are solved iteratively by choosing an initial total vacancy concentration [Eq. (6)], and the system is updated at each time step. For convergence, the time step must be chosen to be smaller than the fastest transition rate. Since no rate can exceed $R0$, the maximum time step expected to converge is

This time step practically limits the total time duration of voltage segments. Since memristors can be programmed using nanosecond pulses, this is not particularly restrictive. We note that Eq. (1) allows one of the three state variables to be eliminated such that it is only necessary to track the occupied and unoccupied vacancy concentration with time (i.e., $NVO2\u2212$ and $NVO2\u2212$). The complexity of the compact model (and its computation speed) depends on the choice of density of states profile for the oxygen vacancy traps as will be shown.

The implementation of the model is as follows:

- $\u2219$The equilibrium Fermi level is set to zero ($EfTE,0=0$), and nominal band offsets and energy levels are defined with respect to the zero energy,(44)$Vfb=\Phi TE\u2212\Phi BE,$(45)$EcTE/BE,0=\u2212(3neTE/BE\pi 2\u210f322meTE/BE3/2)2/3,$(46)$Et=\Phi BE\u2212\chi \u2212Eion\u2212|Vfb|toxyt,$(47)$\Delta ETE/BE=\Phi TE/BE\u2212\chi .$
- $\u2219$A density of states profile is chosen to represent the oxygen vacancy distributions for all memristors being simulated. The energy levels of oxygen vacancies can be equivalently represented in terms of an ionization energy (relative to the conduction band) or the trap level (relative to the equilibrium Fermi level). The density of states profile is normalized such that the integral (over all energies) is equal to one,Therefore, the creation of additional defects scales all energy levels equally without modification to the shape of the density of states profile.(48)$\u222bg^(Et)dEt=1.$
- $\u2219$The initial concentration of oxygen vacancies is set. This value would represent the concentration of defects present in as-prepared oxide films. A nominal value used is 1 × 10
^{10}cm^{−3},(49)$NVO2+(0)=NVO2\u2212(0)=1\xd71010cm\u22123,$(50)$N0(0)=Nsites\u2212NVO2+(0)\u2212NVO2\u2212(0).$ - $\u2219$
The voltage vs time profile is defined. In this work, we consider positive and negative linear ramps applied in sequence to assess forming, reset, and set operations.

- $\u2219$The Fermi level, Fermi functions, and transition rates are defined as functions of the instantaneous voltage and the trap level. These functions are called at each time step as the voltage changes,(51)$EfTE=EfTE,0\u22120.5qV,$(52)$EfBE=EfBE,0+0.5qV,$(53)$EcTE=EcTE,0\u22120.5qV,$(54)$EcBE=EcBE,0+0.5qV,$(55)$fTE/BE=1exp\u2061(Et\u2212EfTE/BEkT)+1.$
- $\u2219$
The state variables are updated using a forward Euler method to solve Eqs. (9), (11), and (12).

- $\u2219$
The current components are computed at each time step from the state variables.

- $\u2219$
Calculated quantities are saved at specified intervals and written to a file.

- $\u2219$
A MATLAB script is used to implement device equations for testing/validation, and Verilog-A code is implemented for circuit simulation in HSPICE, Cadence, and other simulation tools.

## III. RESULTS AND DISCUSSION

### A. General model behavior

#### 1. Density of states profile

In this section, we show how the total current density and current components vary with choice of density of states (DOS) profile for the oxygen vacancy states. Since oxygen vacancy states are electronically coupled to the electrodes, energy levels are necessarily broadened into continuous profiles. The ionization energy of a vacancy is also expected to vary with structural or chemical in-homogeneity, also leading to continuous profiles. For circuit simulation, the choice of DOS profile is important since it has an effect on the computational speed; any continuous DOS profile requires numerical integration at each time step for the calculation of capture/emission rates. The simplest case—the delta function DOS profile—is of particular interest for circuit simulation since this model can be implemented without numerical integration at each time step as will be shown. These considerations present an important trade-off between model complexity and computational efficiency, which can be selected by the needs of the circuit designer.

Figure 2 shows the current density vs voltage for the forming operation. Three different DOS profiles are compared, defined in terms of the ionization energies relative to the conduction band minimum $Eion$: delta profile (situated at 2.9 eV), uniform profile (throughout bandgap), and Gaussian profile (mean = 2.9 eV, standard deviation = 0.33 eV). Figures 2(a), 2(c), and 2(e) correspond to an oxide thickness of 10 nm, whereas Figs. 2(b), 2(d), and 2(f) correspond to an oxide thickness of 5 nm.

For larger oxide thicknesses (i.e., 10 nm), the Ohmic conduction model dominates the total current and is independent of the choice of the DOS profile. This is expected since tunneling current (band to band and trap-assisted) decreases exponentially with oxide thickness and bulk contributions become more dominant. Ohmic conduction, by contrast, is a local bulk component dependent upon the total defect concentration as opposed to the defect energy levels and/or their positions relative to the electrodes.

As the oxide thickness reduces (i.e., 5 nm), the importance of tunneling components increases, particularly during the transition from the pre-forming high-resistance state to the post-forming low-resistance state. Furthermore, the voltage dependence of the trap-assisted tunneling current $Jtat$ becomes increasingly linear as the DOS profile is broadened from a delta profile to that of a uniform and Gaussian profile. This is to be expected since for small voltages, the current is proportional to the emission rate of a single trap state (if a delta profile is assumed), which is exponentially dependent on voltage. A continuous DOS allows for emission from multiple states at a given voltage, broadening the current–voltage dependence computed using Eq. (37). In our previous work, a tunneling current was only considered, which produced linear current–voltage characteristics due to the varied contribution of defects in terms of their position and/or energy level.^{45} We also note that the TAT current saturates for a single energy level [as seen in Figs. 2(a) and 2(b)] when the applied voltage is large enough that the thermal barrier to emission becomes zero (or negative). Beyond this voltage, the emission rate is equal to the pre-exponential factor,

Importantly, these results suggest that, provided that the oxide thickness is not too small where tunneling significant occurs, the Ohmic conduction dominates the total current independent from the DOS profile. This allows a delta DOS profile to be used, which is more computationally efficient, despite being less representative of the physical reality as previously discussed. Importantly, setting the DOS profile as the delta profile,

Defining a net capture $Rc=fBERcBE+fTERcTE$ and emission rate $Re=(1\u2212fBE)ReBE+(1\u2212fTE)ReTE$, the rate equations can be succinctly and intuitively expressed,

The expression for the empty states $N0$ remains unchanged [Eq. (11)]. Together, these represent a condensed, intuitive, and more computationally efficient model that can represent the entirety of memristor characteristics as we show in Fig. 3. Here, the model capability is compared with and without self-heating, showing forming, reset, and set operations for a $10nm$ oxide thickness. Self-heating has the effect of increasing the extent to which the device resets to a high-resistance state. The higher temperature during set shifts the set voltage toward lower values [Fig. 3(c)]. We have assumed that the entire forming, reset, and set operations occur in sequence without delay [Fig. 3(b)]. If there is time between programming cycles (above $\u2248$ 1 ms), the device will cool back to the ambient temperature $T0$ [Fig. 3(c)] during “idling,” and the cumulative-effects of self-heating would be diminished. The time constant for recovery to the ambient temperature can be calculated from setting the Joule-heating term in Eq. (38) to zero, obtaining the familiar Newton’s cooling law result,

Using the parameter values in Table I, $\tau \u224860\mu s$. Assuming that the temperature has recovered within a total time equal to $3\tau \u22480.18ms$, the results shown in Fig. 3(c) are sensible. The exact rate of recovery is expected to differ due to differences in memristor and electrode materials and device geometry.

#### 2. State variable dynamics

In this section, we show how the state variables change with time. The essential aspects of the memristor model will be shown using the forming operation, to avoid restating similar trends that occur from simply changing the polarity/magnitude of voltage. Importantly, during forming, the current exhibits a sharp increase from a low value to a high value at a certain voltage, and the current–voltage characteristics exhibit hysteresis and zero-crossing. In our model, an increase in current stems from an increase in oxygen vacancy concentration—through defect generation and electron capture. By contrast, a decrease in current stems from a decrease in vacancies—through recombination of unoccupied vacancies and electron emission from occupied vacancies.

Figure 4 shows the evolution of the state variables with time for DC and transient conditions for a forming operation from 0 to 5 V. As shown in Figs. 4(a) and 4(b), the theoretical DC solution does not show hysteresis; the absence of hysteresis is evidenced by the fact that the state variables and the current density are symmetric with respect to the applied voltage vs time waveform. This observation is common among physics-based memristor models, which has led to debates regarding whether these devices can, strictly speaking, be classified as memristors. However, as argued by Wang and Roychowdhury,^{59} since the time scales needed to reach the theoretical DC solution are astronomically large, such arguments are not restrictive in practice and may be only of academic interest.

Furthermore, the lack of DC hysteresis suggests that the low-resistance state of a memristor is metastable; the high-resistance state can be recovered by simply raising the device temperature. The DC current is negligible for small voltages because vacancies are unstable to recombination and/or electron emission in the infinite time (i.e., DC) or the infinite temperature limit. This has been shown experimentally by Kumar *et al.*^{32}—oxygen diffusion ultimately leading to recombination of Frenkel pairs at vacancy sites with full reset occuring over a time scale of minutes at temperatures of $\u2248$550 K. In their work, the thermal barrier to recovery of the pre-forming high-resistance state is 0.7 eV, which is comparable to the thermal barrier to oxygen ion diffusion.^{45} Assuming an attempt frequency of 10 THz, at a temperature of 550 K, the time scale for diffusion would be $\u2248$0.2 *μ*s, which is much faster than the time scale observed for reset. Our model also considers electron emission as a parallel process whereby a conductive filament becomes more susceptible to dissolution [see Eq. (9)] since it has been shown that electron-occupied vacancies are more stable to recombination.^{43} Using nominal parameters, our model calculates an emission rate of $\u2248$1 mHz to an assumed $TiN$ electrode from a trap with ionization energy of 3 eV evaluated at low positive voltages $\u2248$0.1 V. The time scale of this process is of the order of minutes and is more consistent with the experimental data in Ref. 32 than a recombination process acting alone. It would be interesting to repeat this experiment for different electrode work functions to see if the thermal barrier is work function dependent and, therefore, if an electron-emission process has a rate-limiting effect on the recovery process.

Transient simulations, by contrast, do show hysteresis. Ramped conditions favor electron capture as opposed to emission as a precursor for recombination, which is a comparatively much slower process for deep traps; electron capture forms a negatively charged vacancy, which is more stable to recombination.^{43} This is shown in Figs. 4(c) and 4(d) for a ramp rate of 10 V ms^{−1}. This ramp rate is much slower, by many orders of magnitude, than ramp rates used during programming cycles in practice (e.g., write operations), which is of the order of $>1Vns\u22121$. As shown in Fig. 4(f), the emission rate remains lower than the capture rate (by many orders of magnitude) at low voltages until $\u2248$1.3 V. At this voltage, the Fermi level of the top electrode is approximately aligned with the trap level at $Eion\u22483eV$, increasing the emission rate to become comparable to the capture rate. When the capture and emission rates are both high—for positive voltages—the rate of recombination is low, preventing reset; for negative voltages, the recombination rate is high, leading to reset. Eventually, the generation rate becomes substantial at large positive voltages, causing the vacancy concentration to increase. On the reverse sweep, the emission rate falls off, stabilizing a large concentration of occupied vacancies and a corresponding low-resistance state. Thus, for all practical purposes, the model can be regarded as capable of reproducing hysteresis and would accurately predict the vanishing of hysteresis at temperatures exceeding $\u2248$550 K and for pulses exceeding minutes in duration.

In both cases (DC and transient), the transition from a pre-forming high-resistance state to a post-forming low-resistance state is coincident with a reduction in the concentration of empty states $N0$ and a corresponding increase in vacancies $NVO2+$ and $NVO2\u2212$. Qualitatively, this is similar in effect to a reduction in the tunneling gap according to previous models.^{17} On the forward sweep, due to the large difference in the rate of electronic processes and atomic processes, the increase in $NVO2\u2212$ tracks the exponential increase in $NVO2+$ since electron capture is energetically favorable from the bottom electrode for positive applied voltage $V=VTE\u2212VBE>0$. On the reverse sweep, the concentration of occupied vacancies is stable due to a sizeable thermal barrier to emission; the barrier is stable since the trap level is lower than the Fermi energy of either a electrode (i.e., $Et\u2212Ef<0$). The new value of $NVO2+$ that is reached once the voltage at the top electrode returns to zero is based on a balance between the rate of recombination of Frenkel pairs $Rrec$ and the combined rate of electron capture $Rc$ according to Eq. (12). This can be seen from the increase in the concentration of empty states and a decrease in unoccupied vacancies toward the end of the pulse $0.9ms<t<1ms$.

We note that a particular strength of our approach is that, based on Eq. (1), the state variables are bound between the atomic density $Nsites$ and zero, provided the rates are finite. As pointed out by Wang and Roychowdhury,^{59} the use of a physical description based on a tunneling gap as a state variable requires the use of a window function to prevent the gap from diminishing below zero or extending beyond the oxide thickness. Using the oxygen vacancy concentrations as state variables eliminates the need for such a window function.

#### 3. Predictive capability of model

Figure 5 shows the predictive capability of the model for variable temperatures (a) and (b), initial vacancy concentration (c) and (d), and activation enthalpy (e) and (f). Pulse durations used were 10 , 1 , and 10 *μ*s for forming, reset, and set, respectively.

Temperature was varied from 300 to 1000 K. Temperature changes affect the onset forming and set voltage and the resistance level following reset. At higher temperatures, the voltage required for forming is reduced; this indicates an increase in the rate of forming kinetics as predicted by the thermochemical model of dielectric breakdown;^{48} this has also been confirmed experimentally by measuring the time required to form at various temperatures.^{60} Importantly, a higher temperature increases the rate of electron emission from occupied vacancies, creating more unoccupied vacancies, which are available for recombination with oxygen ions, therefore increasing the rate of reset.

The initial vacancy concentration was varied from 1×10^{10} to 1 × 10^{18} cm^{−3}. An increase in vacancy concentration has the effect of modifying the high-resistance state of the oxide film prior to forming only; as expected, a higher defect concentration leads to a higher current density and a lower resistance level; this parameter can be adjusted to reflect differences in preparation techniques used to deposit oxide thin films, which will change due to processing parameters.

Last, the activation enthalpy of Frenkel-pair generation is modified from 4 to 8.5 eV. As expected from Eq. (16), the voltage required for forming linearly increases with an increase in activation enthalpy. Again, this follows from the thermochemical model of dielectric breakdown;^{48,49} this parameter can be adjusted to match forming voltages extracted from experimental datasets. To match set voltages, a separate activation enthalpy is used (not shown).

In Sec. III B, we illustrate a choice of these parameters using an experimental dataset.

### B. Model validation and testing

#### 1. Comparison to experimental data

Having demonstrated the general behavior of the model, we now utilize this framework to develop a circuit-compatible compact model in Verilog-A based upon a delta DOS profile.

First, we compare our compact model (using a delta DOS profile) with two different $HfO2$-based memristors having different electrode materials in which the $HfO2$ layer was identically processed using atomic layer deposition at 350 °C.^{60} This dataset was selected based on the availability of processing and characterization data, which provides details regarding the structure of the $HfO2$ layer in addition to the electrode interfaces.^{60–62} The memristors have identical geometry (area = 1.25 × 10^{-13} m^{2}, oxide thickness = 10 nm) with an electrode thickness of 25 nm. The devices have symmetric electrode configurations $TiN$/$HfO2$/$TiN$ [shown in Fig. 6(a)] and $Pt$/$HfO2$/$Pt$ [shown in Fig. 6(b)]. Apart from an apparent negative differential resistance in the as-prepared (pre-forming) state present in both but particularly pronounced in the $Pt$/$HfO2$/$Pt$ device (likely due to initial charging of interface traps^{63,64}), the devices are well-behaved. We note that other devices prepared using the same process show similar forming voltages and reduced current magnitudes without NDR behavior^{65} and are, therefore, within the large statistical spread of these devices^{60} and qualitatively identical. Electrode parameters were used from the prior literature reports for $TiN$^{66–68} and $Pt$^{69–71} thin films.

The simulated current–voltage characteristics for these two $HfO2$ memristors are compared with experimental data reported in Ref. 60 in Figs. 6(a) and 6(b). As evidenced by the smaller current in the pre-forming state, the $Pt$ electrode device is fit using a smaller initial oxygen vacancy concentration. The authors report a higher oxygen vacancy concentration in $HfO2$ films fabricated on $TiN$ (and $Ti$) as opposed to films grown on $Pt$, attributed to the formation of a titanium suboxide ($TiOx$) interfacial layer. This is sensible from a thermodynamics perspective given the large formation enthalpy of $TiO2$ and the associated gettering behavior of $Ti$.^{72} Importantly, from a modeling perspective, this justifies the increase in the initial oxygen vacancy concentration in $HfO2$ films deposited on $TiN$.

The forming and set voltages are higher for the $Pt$ electrode device than for the $TiN$ electrode device. The authors attribute the reduction in forming voltage to a reduced formation enthalpy of Frenkel pairs within the vacancy-rich interfacial layer adjacent to the $TiN$ electrodes. The lower formation enthalpy may reflect the fact that nearest-neighbor pairs and clusters of oxygen vacancies are more stable than point defects;^{43} vacancy clusters ultimately produce what is regarded as a “conductive filament” and are present after initial forming. For the forming data shown in Fig. 6(a), the forming voltage is on the higher end of the statistical variation reported for $Ti$/$TiN$ based electrodes, which can be much lower by a couple of volts.^{61} Although this explanation is consistent with theoretical investigations in the formation of vacancy clusters, it should be noted that the formation of a positive space charge (e.g., unoccupied oxygen vacancies) near the anode/TE or a negative space charge (e.g., negative oxygen vacancies) near the cathode/BE would be expected to lead to a local field enhancement within the oxide—both reducing the forming voltage. Additionally, the stoichiometric changes (when corrected for improvements to polarizability due to crystallization^{73}) can possibly lead to changes in forming and set voltages according to the thermochemical model [Eq. (16)]

Thus, either a reduction in formation enthalpy or an increase in polarization (e.g., relative permittivity or molecular dipole moment) could be used to describe the differences in transition voltages in these datasets. Here, to describe the forming/set operations, the activation energy for Frenkel-pair generation was reduced. We note that the transition voltages are sweep rate dependent; we have used linear voltage ramps over the indicated voltage ranges with a total time interval per segment of 10 *μ*s for forming and set and 1 *μ*s for reset operations in both devices corresponding to ramp rates between 0.1 and 1.2 MV s^{−1}, consistent with the reports for this dataset.^{61}

The authors also report the presence of different crystalline phases between the two $HfO2$ films—films grown on $Ti$ or $TiN$ are predominantly orthorhombic, whereas films grown on $Pt$ are predominantly monoclinic. The monoclinic phase is the least polarizable among the polymorphs of $HfO2$ with a relative permittivity of 18 and 29 for the cubic phase.^{74} The dominant orthorhombic phase observed on $TiN$ is, therefore, consistent with a reduced forming/set voltage relative to monoclinic films grown on $Pt$ without the need for reducing the activation enthalpy. The oxygen vacancy formation energy in both phases is similar at the same Fermi level position;^{75} therefore, the ionization energy is kept the same when fitting the datasets. Only a slight increase in capture cross section (from 1 × 10^{−14} to 2 × 10^{−14} cm^{2}) was needed to fit the $Pt$ devices.

In sum, for both datasets, our compact model demonstrates a close match with the experiments in all operating conditions (forming, set, and reset) of a memristor. The many physical considerations emphasize the need for detailed characterization in order to leverage the model’s predictive capability and justify changes to particular parameters over others. Table II lists the parameter values adjusted to achieve good agreement between the modeled and measured data.

. | $NVO2+(0)$ (cm^{−3})
. | N_{sites} (cm^{−3})
. | E_{ion} (eV)
. | σ_{0} (cm^{2})
. | $Ea,genforming$ (eV) . | $Ea,genset$ (eV) . |
---|---|---|---|---|---|---|

$TiN$/$HfO2$/$TiN$ | 5 × 10^{11} | 4.38 × 10^{19} | 2.957 | 1 × 10^{−14} | 7.35 | 1.90 |

$Pt$/$HfO2$/$Pt$ | 1.61 × 10^{10} | 4.38 × 10^{19} | 2.957 | 1.7 × 10^{−14} | 7.91 | 4.56 |

. | $NVO2+(0)$ (cm^{−3})
. | N_{sites} (cm^{−3})
. | E_{ion} (eV)
. | σ_{0} (cm^{2})
. | $Ea,genforming$ (eV) . | $Ea,genset$ (eV) . |
---|---|---|---|---|---|---|

$TiN$/$HfO2$/$TiN$ | 5 × 10^{11} | 4.38 × 10^{19} | 2.957 | 1 × 10^{−14} | 7.35 | 1.90 |

$Pt$/$HfO2$/$Pt$ | 1.61 × 10^{10} | 4.38 × 10^{19} | 2.957 | 1.7 × 10^{−14} | 7.91 | 4.56 |

#### 2. Simulation of a 4 × 4 1T1R nonvolatile memory array

Finally, to demonstrate the usability of our model and prove that our model can be self-consistently coupled with other circuit elements, we simulate a $4\xd74$ 1-transistor 1-memristor (1T1R) array [shown in Fig. 7(a)]. In our simulation, we use the 65 nm DGXFET models available in the IBM CMOS 10LPe process for the access transistors and our compact model calibrated with the experiment^{60} for the memristor-based memory elements. In memristor-based memory systems, two resistance levels—HRS (high-resistance state) and LRS (low-resistance state)—are used to define the memory states (logic “0” and logic “1,” respectively). Here, we demonstrate the write and read operations in the (2, 2) cell of the array. To access the cell, we fist apply the suitable voltage to the corresponding WL $(WL2)$ to turn the access transistor ON [see top panels of Figs. 7(b)–7(e)]. Then, to write into (read from) the accessed cell, we apply the suitable write voltage (read voltage) to the BL $(BL2)$. The remaining WLs and BLs are kept at ground.

Figures 7(b) and 7(c) show that with the suitable biases at WLs and BLs, the resistance of the accessed cell switches from one state to another: HRS $\u2192$ LRS during write “1” and LRS $\u2192$ HRS during write “0” operations. During the read operation, for the same amount of voltage at the BL, we observe two levels of current based on the stored memory state in the cell–high current for logic “1” and low current for logic “0” states. A simple current sense amplifier^{76,77} with a suitable reference can be used for the sensing purpose thanks to the distinct current levels corresponding to two memory states.

## IV. CONCLUSIONS

In summary, a new compact model for oxide memristors was presented based on the use of oxygen vacancy concentration as a state variable. The theoretical model is based upon a recent paper, which combined rates of atomic processes (e.g., Frenkel-pair generation and recombination, diffusion of point defects) with those of electronic processes (electron capture and emission) in a kinetic Monte Carlo simulation approach. The compact model was validated using bulk parameters for $HfO2$, though it is of general utility to a wide range of single phase, insulating oxides or polymorphs, which can be described using “effective” bulk-averaged parameters. The dynamic evolution of the state variables—the concentrations of occupied and unoccupied oxygen vacancies—was shown to be capable, both qualitatively and quantitatively, of reproducing the essential switching characteristics of oxide memristors. In particular, the increase (decrease) in oxygen vacancy concentration is qualitatively similar in effect to the more familiar reduction (expansion) of the tunnel gap, which has been used in existing compact models. Key strengths of this approach include the following: the state variables are bound between zero and the atomic density of the oxide without the use of a window function—needed for preventing tunnel gaps from becoming negative or exceeding the oxide thickness; the model makes few assumptions regarding the “filament” geometry, being chiefly based upon the number density of oxygen vacancy defects as opposed to their geometry; the model provides an intuitive description of resistive switching that is consistent with retention improvements observed in lower work-function metal electrodes (e.g., $Ti$/$TiN$) compared to $Pt$ due to an increased emission barrier, which stabilizes the low-resistance state; the model can be easily refined via inclusion of image-force correction or Poole–Frenkel modeling of high-resistance states without loss of generality.

As a test of its practical utility, the compact model was implemented in Verilog-A, verified using well-defined experimental datasets with different electrode geometries, and tested in circuit simulation using a $4\xd74$ one-transistor, one-resistor (1T1R) memory cell, representative of a small, nonvolatile memory array. We anticipate that this model will serve as an alternate description of memristor switching in terms of the concentration of defects as opposed to the state of the filament. Such a description will be particularly useful as new and ongoing experimental data emerge to suggest new or confirm existing geometrical descriptions of conductive filament shapes.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Ethics Approval

Ethics approval is not required.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*2015 International Conference on Memristive Systems, MEMRISYS 2015*(IEEE, 2016).

*IEEE International Conference on Design and Test of Integrated Micro and Nano-Systems, DTS 2019*(IEEE, 2019).

*HPCA—16 2010 The Sixteenth International Symposium on High-Performance Computer Architecture*(IEEE, 2010).

*Proceedings of IEEE IEDM*(IEEE, 2008), pp. 1–4.

*Proceedings of IEEE International Memory Workshop*(IEEE, 2012), pp. 1–4.

_{2}-based resistive switching memories,” in

*Proceedings of IEEE International Memory Workshop*(IEEE, 2011).

*Proceedings of International Reliability Physics Symposium*(IEEE, 2010), pp. 792–798.

*Proceedings of International Reliability Physics Symposium IEEE*(IEEE, 2013), pp. 5A.6.1–5A.6.8.

*International Conference on Simulation of Semiconductor Processes and Devices (SISPAD)*(IEEE, 2014), pp. 41–44.

*IEEE International Symposium on Circuits and Systems (ISCAS)*(IEEE, 2017).

_{2}/Si-n

^{+}-based RRAMs

_{2}memristors: Efficient device-circuit Co-design for neuromorphic systems

_{2}-based resistive random access memory by in situ TEM studies

_{2}using random telegraph noise signals from scanning tunneling microscopy

*Proceedings of the 2018 12th Spanish Conference on Electron Devices (CDE)*, edited by J. Mateos and T. Gonzalez (IEEE, New York, 2018).

*2012 International Electron Devices Meeting*(IEEE, 2012), pp. 20.1.1–20.1.4.

*2011 International Electron Devices Meeting*(IEEE, 2011), pp. 28.7.1–28.7.4.

*2017 75th Annual Device Research Conference (DRC)*(IEEE, 2017), pp. 1–2.

*et al.*, “