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×4 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.

Several physics-inspired compact models have been proposed for circuit simulation of oxide memristors1–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 gap14–18 or the size/shape of the filament approximated as cylindrical19 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 retention32 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 conduction34) 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 sites36), 2.2 eV (at gettering electrode surfaces37), and several electronvolts in the bulk (3.6 ,36 6.8 eV37), all of which are much larger than that of oxygen ion diffusion (0.7 eV32). 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 approaches39,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 approaches45 and reflects ab initio calculations showing increased thermodynamic stability of negatively charged oxygen vacancies under conditions of electron-injection43 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 spectroscopy47). 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×4 1T1R array representative of a small nonvolatile memory array to evaluate the model’s capability to perform basic circuit simulations.

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.

FIG. 1.

(a) Device structure, showing the physical processes and memristor components modeled herein. (b) Band diagram showing definitions of energy parameters (not to scale).

FIG. 1.

(a) Device structure, showing the physical processes and memristor components modeled herein. (b) Band diagram showing definitions of energy parameters (not to scale).

Close modal
TABLE I.

Summary of nominal parameters used in this work unless otherwise stated.

ParameterDescriptionValue and unit
q Elementary charge 1.602 176 × 10−19 C 
k Boltzmann constant 1.380 649 × 10−23 J K−1 
 Reduced Planck constant 1.054 571 × 10−34 J s 
T0 Ambient temperature 300 K 
ɛr Relative permittivity 30 
R0 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 
Eion Ionization energy of oxygen vacancy defects 2.0−3.0 eV 
σ0 Capture cross section of oxygen vacancies 1 × 10−14 cm2 
ΦTE TE work function (Ti4.33 eV 
ΦBE BE work function (TiN4.50 eV 
meTE TE electron effective mass (Ti3.2 
meBE BE electron effective mass (TiN2.0 
neTE TE electron concentration (Ti5.67 × 1022 cm−3 
neBE BE electron concentration (TiN2.88 × 1022 cm−3 
meox Oxide electron effective mass (HfO20.1 
p Oxide molecular dipole moment (HfO215 × 10−10 C m 
Eg Oxide bandgap energy (HfO25.9 eV 
χ Oxide electron affinity (HfO22 eV 
μeff Effective band mobility (HfO21 cm2 V−1 s−1 
cp,ox Oxide specific heat capacity (HfO2120 J kg−1 K−1 
Cadd. Addenda heat capacity (device) 1.45 × 10−9 J K−1 
ρox Oxide mass density (HfO29800 kg m−3 
κox Oxide thermal conductivity (HfO21 W m−1 K−1 
Nsites Site density (HfO24.38 × 1019 cm−3 
yt Trap center (HfO20.5tox 
ParameterDescriptionValue and unit
q Elementary charge 1.602 176 × 10−19 C 
k Boltzmann constant 1.380 649 × 10−23 J K−1 
 Reduced Planck constant 1.054 571 × 10−34 J s 
T0 Ambient temperature 300 K 
ɛr Relative permittivity 30 
R0 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 
Eion Ionization energy of oxygen vacancy defects 2.0−3.0 eV 
σ0 Capture cross section of oxygen vacancies 1 × 10−14 cm2 
ΦTE TE work function (Ti4.33 eV 
ΦBE BE work function (TiN4.50 eV 
meTE TE electron effective mass (Ti3.2 
meBE BE electron effective mass (TiN2.0 
neTE TE electron concentration (Ti5.67 × 1022 cm−3 
neBE BE electron concentration (TiN2.88 × 1022 cm−3 
meox Oxide electron effective mass (HfO20.1 
p Oxide molecular dipole moment (HfO215 × 10−10 C m 
Eg Oxide bandgap energy (HfO25.9 eV 
χ Oxide electron affinity (HfO22 eV 
μeff Effective band mobility (HfO21 cm2 V−1 s−1 
cp,ox Oxide specific heat capacity (HfO2120 J kg−1 K−1 
Cadd. Addenda heat capacity (device) 1.45 × 10−9 J K−1 
ρox Oxide mass density (HfO29800 kg m−3 
κox Oxide thermal conductivity (HfO21 W m−1 K−1 
Nsites Site density (HfO24.38 × 1019 cm−3 
yt Trap center (HfO20.5tox 

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

  • 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.

  • The local electric field can be approximated with the spatially averaged field.

  • Oxygen vacancies behave as electron traps (positive when empty VO2+, negative when occupied VO2, accommodating 4 electrons43), defined in terms of a continuous electronic density of states to specify the ionization energies/defect levels participating in electron capture/emission.

  • The interface tunneling transmission probabilities for the top and bottom electrodes are approximately equal.

  • Diffusion of ionic point defects is not considered.

  • Trap-to-trap coupling is not considered.

The state of the oxide is represented using three species: unoccupied oxygen vacancies NVO2+, occupied oxygen vacancies NVO2, 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 (Nsites2.77×1022cm3),

N0+NVO2++NVO2=Nsites.
(1)

We note that, in general, Nsites can be regarded as a parameter that establishes an upper bound for the state variables NVO2+, NVO2, 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 × 1019 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,

dN0dt+dNVO2+dt+dNVO2dt=0.
(2)

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),

dftdt=(1ft)(fBERcBE+fTERcTE)Rcft((1fBE)ReBE+(1fTE)ReTE)Re=(1ft)RcftRe.
(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,

NVO2=ftg(Et)dEt,
(4)
NVO2+=(1ft)g(Et)dEt,
(5)
NVO=NVO2++NVO2=g(Et)dEt.
(6)

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,

ft=NVO2NVO,
(7)
1ft=NVO2+NVO.
(8)

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,

dNVO2dt=1NVOg(Et)(NVO2+RcNVO2Re)dEt=g^(Et)(NVO2+RcNVO2Re)dEt,
(9)

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

g^(Et)=g(Et)g(Et)dEt=g(Et)NVO.
(10)

Equation (9) simply states that unoccupied vacancy states NVO2+ are filled by either electrode (TE/BE), and occupied states NVO2 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, 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,

dN0dt=RrecNVO2+RgenN0.
(11)

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

dNVO2+dt=RgenN0RrecNVO2+dNVO2dt=RgenN0RrecNVO2+g^(Et)(NVO2+RcNVO2Re)dEt.
(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.

The rates Rgen and Rrec are atomic processes describing the creation and removal of Frenkel pairs of oxygen ions and vacancies (i.e., O2 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,

Rgen=R0exp(Ea,genkT),
(13)
Rrec=R0exp(Ea,reckT).
(14)

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,

Ea,gen=Ea,gen0E|p|(εr+23),
(15)
Ea,rec=Ea,rec0+E|p|(εr+23).
(16)

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=EcEion. 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 work45 and elsewhere50,51 as implemented in standard device simulators.52 

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

RcTE/BE=Ptun.Rinter.=exp(yty0TE/BE)84π(meTE/BE)3/2×VtEion2EtEcTE/BEf(Et)TE/BE.
(17)

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 as50,51

Vt=(2meoxEion)3,
(18)
y0=342meTE/BEΔETE/BE,
(19)

where ΔETE/BE=ΦTE/BEχ is the conduction band offset at each electrode interface, ΦTE/BE is the electrode work function, and χ 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,

RcTE/BE=Rc,0TE/BEexp(yty0TE/BE)exp(Ea,tkBT){exp(qEytEtfkBT),EtfTE/BE>0,exp(qEytkBT),EtfTE/BE<0,
(20)

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

Rc,0TE/BE=σ0vthTE/BEneTE/BE,
(21)
EtfTE/BE=EtEfTE/BE,
(22)
vthTE/BE=3kTmeTE/BE.
(23)

Regarding the in front of the field-dependent terms in Eq. (20), the 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,

ReTE/BE=exp(EfTE/BEEtkT)RcTE/BE.
(24)

For the classical model, this leads to the following:

ReTE/BE=Rc,0TE/BEexp(yty0TE/BE)exp(Ea,tkBT){exp(±qEytkBT),Etf>0,exp(±qEyt+EtfkBT),EtfTE/BE<0.
(25)

Regarding the ± in front of the field-dependent terms in Eq. (25), the + is used for the top electrode (TE) and the 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 EfBEEfTE=qV,

EfTE=EfTE,00.5qV,
(26)
EfBE=EfBE,0+0.5qV.
(27)

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¯=Vtox. Consequently, the rates depend on the applied voltage but not on time.

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),

J=Jt+JOhmic+Jtat.
(28)

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

Jt=q38πhΔETE/BEE¯2exp(42meox(ΔETE/BE)33qE¯).
(29)

Alternatively, this can be refined with image-force correction using the general Simmons formula55 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 ( denotes a neutral site). Initially, Frenkel-pair generation leads to neutral pairs of oxygen ions O2 and oxygen vacancies VO2+,

VO2++O2.
(30)

After electron capture, negative oxygen ions remain, positive oxygen vacancies NVO2+ are removed, and negative oxygen vacancies VO2 are created—producing a net negative space charge equal to 4e per capture event (assuming four electrons are captured43),

VO2++O2VO2+O2.
(31)

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

VO2++O2VO2+O2+4h+.
(32)

This assumes that the oxygen ions O2 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 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 μeff parameter, with a drift current proportional to the total vacancy concentration,

JOhmic=qμeff(NVO2++NVO2)E¯.
(33)

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 

μeff=μ0exp(E0kT)exp(E(βTγ)),
(34)

where μ0, E0, β, and γ 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

dJtat=qR(y)dyJtat=q0toxR(y)dy,
(35)

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 

R(y)=(RcBEReTERcTEReBERcBE+RcTE+ReBE+ReTE)g(Et)dEt.
(36)

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 yt0.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:

Jtatqtoxg(Et)(RcBEReTERcTEReBERcBE+RcTE+ReBE+ReTE)dEt.
(37)

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).

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,

CdevicedTdtI(t)V(t)κoxAΔx(T(t)T0).
(38)

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 ρox, thermal conductivity κox of the oxide, current density J, and ambient temperature T0,

T(t+Δt)=T(t)+Δttoxcp,oxρox(J(t)V(t)2κoxtoxT(t)T0).
(39)

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 (Δx=0.5tox). The model is approximate, and, in practice, the quantity 2κ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=Δxκ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 100 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,

T(t+Δt)=T(t)+Δttoxcp,oxρoxCF(J(t)V(t)2κoxtoxT(t)T0),
(40)

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

Cdevice=Cox+Cadd.
(41)

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ρoxcp,ox—the addenda behaving as a heat sink; in this expression, the electrode area is A and the oxide volume is Atox,

CF=1+CaddCoxCaddAtoxρoxcp,ox.
(42)

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

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

Δt1R0=1×1013s.
(43)

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 and NVO2). 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:

  • 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,
    Vfb=ΦTEΦBE,
    (44)
    EcTE/BE,0=(3neTE/BEπ2322meTE/BE3/2)2/3,
    (45)
    Et=ΦBEχEion|Vfb|toxyt,
    (46)
    ΔETE/BE=ΦTE/BEχ.
    (47)
  • 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,
    g^(Et)dEt=1.
    (48)
    Therefore, the creation of additional defects scales all energy levels equally without modification to the shape of the density of states profile.
  • 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 × 1010 cm−3,
    NVO2+(0)=NVO2(0)=1×1010cm3,
    (49)
    N0(0)=NsitesNVO2+(0)NVO2(0).
    (50)
  • 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.

  • 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,
    EfTE=EfTE,00.5qV,
    (51)
    EfBE=EfBE,0+0.5qV,
    (52)
    EcTE=EcTE,00.5qV,
    (53)
    EcBE=EcBE,0+0.5qV,
    (54)
    fTE/BE=1exp(EtEfTE/BEkT)+1.
    (55)
  • The state variables are updated using a forward Euler method to solve Eqs. (9), (11), and (12).

  • The current components are computed at each time step from the state variables.

  • Calculated quantities are saved at specified intervals and written to a file.

  • 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.

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.

FIG. 2.

Current density vs voltage comparison for different density of states profiles and oxide thickness. Delta DOS profile (situated at 2.9 eV) for an oxide thickness of (a) 10 and (b) 5 nm; uniform DOS profile (throughout bandgap) for an oxide thickness of (c) 10 and (d) 5 nm; and Gaussian DOS profile (mean = 2.9 eV, standard deviation = 0.33 eV) for an oxide thickness of (e) 10 and (f) 5 nm. Nominal parameters for this simulation are as follows: T0=T(t)=300K, tox=10nm, NVO2+(0)=NVO2(0)=5×1011cm3, Nsites=4.38×1019cm3, Ea,gen=7.35eV (forming), and σ0=1×1014cm2.

FIG. 2.

Current density vs voltage comparison for different density of states profiles and oxide thickness. Delta DOS profile (situated at 2.9 eV) for an oxide thickness of (a) 10 and (b) 5 nm; uniform DOS profile (throughout bandgap) for an oxide thickness of (c) 10 and (d) 5 nm; and Gaussian DOS profile (mean = 2.9 eV, standard deviation = 0.33 eV) for an oxide thickness of (e) 10 and (f) 5 nm. Nominal parameters for this simulation are as follows: T0=T(t)=300K, tox=10nm, NVO2+(0)=NVO2(0)=5×1011cm3, Nsites=4.38×1019cm3, Ea,gen=7.35eV (forming), and σ0=1×1014cm2.

Close modal

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,

Rmax=neTE/BEvthTE/BEσ0exp(yty0TE/BE).
(56)

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,

g(Et)=NVOδ(EtEt0),
(57)

removes the integrals from Eqs. (9) and (37),

dNVO2dt=NVO2+(fBERcBE+fTERcTE)NVO2((1fBE)ReBE+(1fTE)ReTE),
(58)
Jtatqtox(NVO2+NVO2+)(RcBEReTERcTEReBERcBE+RcTE+ReBE+ReTE).
(59)

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

dNVO2dt=NVO2+RcNVO2Re,dNVO2+dt=RgenN0RrecNVO2+NVO2+Rc+NVO2Re.

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 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,

T(t)=T0+T(t0)exp(tt0τ);τ=toxCdevice2κoxA.
(60)

Using the parameter values in Table I, τ60μs. Assuming that the temperature has recovered within a total time equal to 3τ0.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.

FIG. 3.

(a) Nominal current–voltage characteristics computed by delta-shaped density of states profile, with (solid lines) and without self-heating (dashed lines) for an assumed TiN/HfO2/TiN device structure. (b) Voltage vs time waveform. (c) Calculated temperature vs time profile due to self-heating. Parameters for this simulation are as follows: T0=T(0)=300K, tox=10nm, NVO2+(0)=NVO2(0)=5×1011cm3, Nsites=4.38×1019cm3, Ea,gen=7.35eV (forming), Ea,gen=1.9eV (set), σ0=1×1014cm2.

FIG. 3.

(a) Nominal current–voltage characteristics computed by delta-shaped density of states profile, with (solid lines) and without self-heating (dashed lines) for an assumed TiN/HfO2/TiN device structure. (b) Voltage vs time waveform. (c) Calculated temperature vs time profile due to self-heating. Parameters for this simulation are as follows: T0=T(0)=300K, tox=10nm, NVO2+(0)=NVO2(0)=5×1011cm3, Nsites=4.38×1019cm3, Ea,gen=7.35eV (forming), Ea,gen=1.9eV (set), σ0=1×1014cm2.

Close modal

In Secs. III A 2 and III A 3, we apply this delta-profile DOS model to gain an intuitive understanding of state variable dynamics, parameter dependencies, and predictability, as well as making comparisons to experimental data.

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.

FIG. 4.

State variables and current density vs time for DC (a) and (b) and transient (c) and (d) simulations. Transition rates vs time (e) and voltage (f). Note: Rc and Re are the total capture and emission rates for both electrodes. Nominal parameters for this simulation are as follows: T0=T(t)=300K, tox=10nm, NVO2+(0)=NVO2(0)=5×1011cm3, Nsites=4.38×1019cm3, Ea,gen=7.35eV (forming), and σ0=1×1014cm2.

FIG. 4.

State variables and current density vs time for DC (a) and (b) and transient (c) and (d) simulations. Transition rates vs time (e) and voltage (f). Note: Rc and Re are the total capture and emission rates for both electrodes. Nominal parameters for this simulation are as follows: T0=T(t)=300K, tox=10nm, NVO2+(0)=NVO2(0)=5×1011cm3, Nsites=4.38×1019cm3, Ea,gen=7.35eV (forming), and σ0=1×1014cm2.

Close modal

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 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 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 1 mHz to an assumed TiN electrode from a trap with ionization energy of 3 eV evaluated at low positive voltages 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 >1Vns1. As shown in Fig. 4(f), the emission rate remains lower than the capture rate (by many orders of magnitude) at low voltages until 1.3 V. At this voltage, the Fermi level of the top electrode is approximately aligned with the trap level at Eion3eV, 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 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. 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 tracks the exponential increase in NVO2+ since electron capture is energetically favorable from the bottom electrode for positive applied voltage V=VTEVBE>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., EtEf<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.

FIG. 5.

(a) and (b) Effect of temperature: temperature was varied from 300 to 1000 K in 100 K increments. (c) and (d) Effect of initial vacancy concentration: initial vacancy concentration was varied logarithmically from 1 × 1010 to 1 × 1018 cm−3 in one decade increments. (e) and (f) Effect of activation enthalpy of Frenkel-pair generation: enthalpy was varied from 4 to 8.5 in 0.5 eV increments. The current density is shown in both logarithmic (top) and linear (bottom) scales.

FIG. 5.

(a) and (b) Effect of temperature: temperature was varied from 300 to 1000 K in 100 K increments. (c) and (d) Effect of initial vacancy concentration: initial vacancy concentration was varied logarithmically from 1 × 1010 to 1 × 1018 cm−3 in one decade increments. (e) and (f) Effect of activation enthalpy of Frenkel-pair generation: enthalpy was varied from 4 to 8.5 in 0.5 eV increments. The current density is shown in both logarithmic (top) and linear (bottom) scales.

Close modal

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×1010 to 1 × 1018 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.

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 m2, 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 traps63,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 behavior65 and are, therefore, within the large statistical spread of these devices60 and qualitatively identical. Electrode parameters were used from the prior literature reports for TiN66–68 and Pt69–71 thin films.

FIG. 6.

Comparison of current–voltage characteristics from a compact model and experimental data for two memristor devices with identical geometry (area = 1.25 × 10-13 m2, oxide thickness = 10 nm) having different electrode structures, with an identically processed HfO2 layer:60 (a) TiN/HfO2/TiN and (b) Pt/HfO2/Pt. The compact model (solid lines) shows a reasonable match with the experimentally observed characteristics (open circles).

FIG. 6.

Comparison of current–voltage characteristics from a compact model and experimental data for two memristor devices with identical geometry (area = 1.25 × 10-13 m2, oxide thickness = 10 nm) having different electrode structures, with an identically processed HfO2 layer:60 (a) TiN/HfO2/TiN and (b) Pt/HfO2/Pt. The compact model (solid lines) shows a reasonable match with the experimentally observed characteristics (open circles).

Close modal

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 crystallization73) can possibly lead to changes in forming and set voltages according to the thermochemical model [Eq. (16)]

Ebd=Ea,gen|p|(2+εr3).
(61)

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 cm2) 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.

TABLE II.

Summary of parameters used in fitting the data in Fig. 6.

NVO2+(0) (cm−3)Nsites (cm−3)Eion (eV)σ0 (cm2)Ea,genforming (eV)Ea,genset (eV)
TiN/HfO2/TiN 5 × 1011 4.38 × 1019 2.957 1 × 10−14 7.35 1.90 
Pt/HfO2/Pt 1.61 × 1010 4.38 × 1019 2.957 1.7 × 10−14 7.91 4.56 
NVO2+(0) (cm−3)Nsites (cm−3)Eion (eV)σ0 (cm2)Ea,genforming (eV)Ea,genset (eV)
TiN/HfO2/TiN 5 × 1011 4.38 × 1019 2.957 1 × 10−14 7.35 1.90 
Pt/HfO2/Pt 1.61 × 1010 4.38 × 1019 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×4 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 experiment60 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.

FIG. 7.

(a) Schematic of a 1T1R array where we perform the write and read operations in the (2, 2) cell. Top panels show the time dynamics of WL and BL biases for (a) write “1,” (b) write “0,” (c) read “1,” and (d) read “0” operations in the accessed cell. Bottom panels show the corresponding effects of the WL and BL biases on the accessed cell during (a) write “1,” (b) write “0,” (c) read “1,” and (d) read “0” operations.

FIG. 7.

(a) Schematic of a 1T1R array where we perform the write and read operations in the (2, 2) cell. Top panels show the time dynamics of WL and BL biases for (a) write “1,” (b) write “0,” (c) read “1,” and (d) read “0” operations in the accessed cell. Bottom panels show the corresponding effects of the WL and BL biases on the accessed cell during (a) write “1,” (b) write “0,” (c) read “1,” and (d) read “0” operations.

Close modal

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 LRS during write “1” and LRS 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 amplifier76,77 with a suitable reference can be used for the sensing purpose thanks to the distinct current levels corresponding to two memory states.

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×4 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.

The authors have no conflicts to disclose.

Ethics approval is not required.

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

1.
S.
Ambrogio
,
S.
Balatti
,
D. C.
Gilmer
, and
D.
Ielmini
, “
Analytical modeling of oxide-based bipolar resistive memories and complementary resistive switches
,”
IEEE Trans. Electron Devices
61
,
2378
2386
(
2014
).
2.
S.
Brivio
,
E.
Covi
,
A.
Serb
,
T.
Prodromakis
,
M.
Fanciulli
, and
S.
Spiga
, “Gradual set dynamics in HfO2-based memristor driven by sub-threshold voltage pulses,” in 2015 International Conference on Memristive Systems, MEMRISYS 2015 (IEEE, 2016).
3.
A.
Hardtdegen
,
C.
La Torre
,
F.
Cuppers
,
S.
Menzel
,
R.
Waser
, and
S.
Hoffmann-Eifert
, “
Improved switching stability and the effect of an internal series resistor in HfO2/TiOx bilayer ReRAM cells
,”
IEEE Trans. Electron Devices
65
,
3229
3236
(
2018
).
4.
K.
Mbarek
,
F.
Rziga
,
S.
Ghedira
, and
K.
Besbes
, “A practical HfO2-based OxRAM memristor model suitable for circuit design and simulation,” in IEEE International Conference on Design and Test of Integrated Micro and Nano-Systems, DTS 2019 (IEEE, 2019).
5.
F.
Rziga
,
K.
Mbarek
,
S.
Ghedira
, and
K.
Besbes
, “
An efficient Verilog—A memristor model implementation: Simulation and application
,”
J. Comput. Electron.
18
,
1055
1064
(
2019
).
6.
F.
Alonso
,
D.
Maldonado
,
A.
Aguilera
, and
J.
Roldán
, “
Memristor variability and stochastic physical properties modeling from a multivariate time series approach
,”
Chaos, Solitons Fractals
143
,
110461
(
2021
).
7.
C.
Bengel
,
F.
Cüppers
,
M.
Payvand
,
R.
Dittmann
,
R.
Waser
,
S.
Hoffmann-Eifert
, and
S.
Menzel
, “
Utilizing the switching stochasticity of HfO2/TiOx-based ReRAM devices and the concept of multiple device synapses for the classification of overlapping and noisy patterns
,”
Front. Neurosci.
15
,
661856
(
2021
).
8.
D.
Zhevnenko
,
F.
Meshchaninov
,
V.
Kozhevnikov
,
E.
Shamin
,
O.
Telminov
, and
E.
Gornev
, “
Research and development of parameter extraction approaches for memristor models
,”
Micromachines
12
,
1220
(
2021
).
9.
M. P.
Raj
and
G.
Kavithaa
, “
Memristor based high speed and low power consumption memory design using deep search method
,”
J. Ambient Intell. Hum. Comput.
12
,
4223
4235
(
2021
).
10.
N.
Du
,
H.
Schmidt
, and
I.
Polian
, “
Low-power emerging memristive designs towards secure hardware systems for applications in internet of things
,”
Nano Mater. Sci.
3
(
2
),
186
204
(
2021
).
11.
G.
Sun
,
Y.
Joo
,
Y.
Chen
,
D.
Niu
,
Y.
Xie
,
Y.
Chen
, and
H.
Li
, “A hybrid solid-state storage architecture for the performance, energy consumption, and lifetime improvement,” in HPCA—16 2010 The Sixteenth International Symposium on High-Performance Computer Architecture (IEEE, 2010).
12.
A.
Thomas
, “
Memristor-based neural networks
,”
J. Phys. D: Appl. Phys.
46-9
,
093001
(
2013
).
13.
Q.
Lu
,
F.
Sun
,
L.
Liu
,
L.
Li
,
Y.
Wang
,
M.
Hao
,
Z.
Wang
,
S.
Wang
, and
T.
Zhang
, “
Biological receptor-inspired flexible artificial synapse based on ionic dynamics
,”
Microsyst. Nanoeng.
6
,
84
(
2020
).
14.
L.
Chen
,
C.
Li
,
T.
Huang
,
Y.
Chen
,
S.
Wen
, and
J.
Qi
, “
A synapse memristor model with forgetting effect
,”
Phys. Lett. A
377
,
3260
3265
(
2013
).
15.
Y.
Zhang
,
X.
Wang
,
Y.
Li
, and
E. G.
Friedman
, “
Memristive model for synaptic circuits
,”
IEEE Trans. Circuits Syst. II
64
,
767
771
(
2017
).
16.
I.
Vourkas
,
A.
Batsos
, and
G. C.
Sirakoulis
, “
SPICE modeling of nonlinear memristive behavior
,”
Int. J. Circuit Theory Appl.
43
,
553
565
(
2015
).
17.
X.
Guan
,
S.
Yu
, and
H.-S. P.
Wong
, “
A SPICE compact model of metal oxide resistive switching memory with variations
,”
IEEE Electron Device Lett.
33
,
1405
1407
(
2012
).
18.
P.-Y.
Chen
and
S.
Yu
, “
Compact modeling of RRAM devices and its applications in 1T1R and 1S1R array design
,”
IEEE Trans. Electron Devices
62
(
12
),
4022
4028
(
2015
).
19.
M.
Bocquet
,
D.
Deleruyelle
,
H.
Aziza
,
C.
Muller
,
J.-M.
Portal
,
T.
Cabout
, and
E.
Jalaguier
, “
Robust compact model for bipolar oxide-based resistive switching memories
,”
IEEE Trans. Electron Devices
61
,
674
681
(
2014
).
20.
D.
Biolek
,
V.
Biolkova
, and
Z.
Biolek
, “
SPICE model of memristor with nonlinear dopant drift
,”
Radioengineering
18
,
210
214
(
2009
).
21.
S.
Kvatinsky
,
E. G.
Friedman
,
A.
Kolodny
, and
U. C.
Weiser
, “
TEAM: Threshold adaptive memristor model
,”
IEEE Trans. Circuits Syst. I
60
,
211
221
(
2013
).
22.
B.
Gao
,
S.
Yu
,
N.
Xu
,
L. F.
Liu
,
B.
Sun
,
X. Y.
Liu
,
R. Q.
Han
,
J. F.
Kang
,
B.
Yu
, and
Y.
Wang
, “Oxide-based RRAM switching mechanism: A new ion-transport recombination model,” in Proceedings of IEEE IEDM (IEEE, 2008), pp. 1–4.
23.
S.
Larentis
,
F.
Nardi
,
S.
Balatti
,
D. C.
Gilmer
, and
D.
Ielmini
, “
Resistive switching by voltage-driven ion migration in bipolar RRAM—Part II: Modeling
,”
IEEE Trans. Electron Devices
59
,
2468
2475
(
2012
).
24.
D. C.
Gilmer
,
G.
Bersuker
,
S.
Koveshnikov
,
M.
Jo
,
A.
Kalantarian
,
B.
Butcher
,
R.
Geer
,
Y.
Nishi
,
P. D.
Kirsch
, and
R.
Jammy
, “Asymmetry, vacancy engineering and mechanism for bipolar RRAM,” in Proceedings of IEEE International Memory Workshop (IEEE, 2012), pp. 1–4.
25.
L.
Vandelli
,
A.
Padovani
,
L.
Larcher
,
G.
Bersuker
,
D.
Gilmer
, and
P.
Pavan
, “Modeling of forming operation in HfO2-based resistive switching memories,” in Proceedings of IEEE International Memory Workshop (IEEE, 2011).
26.
F.
Nardi
,
S.
Larentis
,
S.
Balatti
,
D. C.
Gilmer
, and
D.
Ielmini
, “
Resistive switching by voltage-driven ion migration in bipolar RRAM—Part I: Experimental study
,”
IEEE Trans. Electron Devices
59
,
2461
2467
(
2012
).
27.
P.
Huang
,
X. Y.
Liu
,
B.
Chen
,
H. T.
Li
,
Y. J.
Wang
,
Y. X.
Deng
,
K. L.
Wei
,
L.
Zeng
,
B.
Gao
,
G.
Du
,
X.
Zhang
, and
J. F.
Kang
, “
A physics-based compact model of metal-oxide-based RRAM DC and AC operations
,”
IEEE Trans. Electron Devices
60
,
4090
4097
(
2013
).
28.
J.
Suñé
, “
New physics-based analytic approach to the thin-oxide breakdown statistics
,”
IEEE Electron Device Lett.
22
,
294
296
(
2001
).
29.
J.
Suñé
,
S.
Tous
, and
E. Y.
Wu
, “
Analytical cell-based model for the breakdown statistics of multilayer insulator stacks
,”
IEEE Electron Device Lett.
30
,
41359
1361
(
2009
).
30.
S.
Tous
,
E. Y.
Wu
, and
J.
Suñé
, “A compact analytic model for the breakdown distribution of gate stack dielectrics,” in Proceedings of International Reliability Physics Symposium (IEEE, 2010), pp. 792–798.
31.
S.
Long
,
X.
Lian
,
C.
Cagli
,
L.
Perniola
,
E.
Miranda
,
D.
Jiménez
,
H.
Lv
,
Q.
Liu
,
L.
Li
,
Z.
Huo
,
M.
Liu
, and
J.
Suñé
, “Compact analytical models for the SET and RESET switching statistics of RRAM inspired in the cell-based percolation model of gate dielectric breakdown,” in Proceedings of International Reliability Physics Symposium IEEE (IEEE, 2013), pp. 5A.6.1–5A.6.8.
32.
S.
Kumar
,
Z.
Wang
,
X.
Huang
,
N.
Kumari
,
N.
Davila
,
J.
Strachan
,
D.
Vine
,
A.
Kilcoyne
,
Y.
Nishi
, and
R.
Williams
, “
Conduction channel formation and dissolution due to oxygen thermophoresis/diffusion in hafnium oxide memristors
,”
ACS Nano
10
,
11205
11210
(
2016
).
33.
A.
Alexandrov
,
A.
Bratkovsky
,
B.
Bridle
,
S.
Savel’ev
,
D.
Strukov
, and
R. S.
Williams
, “
Current-controlled negative differential resistance due to Joule heating in TiO2
,”
Appl. Phys. Lett.
99
,
202104
(
2011
).
34.
S. M.
Sze
and
K. K.
Ng
,
Physics of Semiconductor Devices
(
John Wiley & Sons
,
2006
).
35.
S.
Yu
and
H.-S. P.
Wong
, “
A phenomenological model for the reset mechanism of metal oxide RRAM
,”
IEEE Electron Device Lett.
31
,
1455
1457
(
2010
).
36.
S.
Aldana
,
P.
Garcia-Fernandez
,
R.
Romero-Zaliz
,
M. B.
Gonzalez
,
F.
Jimenez-Molinos
,
F.
Gomez-Campos
,
F.
Campabadal
, and
J. B.
Roldan
, “
Resistive switching in HfO2 based valence change memories, a comprehensive 3D kinetic Monte Carlo approach
,”
J. Phys. D: Appl. Phys.
53
,
225106
(
2020
).
37.
X.
Xu
,
B.
Rajendran
, and
M. P.
Anantram
, “
Kinetic Monte Carlo simulation of interface-controlled hafnia-based resistive memory
,”
IEEE Trans. Electron Devices
67
,
118
124
(
2020
).
38.
Z.
Jiang
,
S.
Yu
,
Y.
Wu
,
J. H.
Engel
,
X.
Guan
, and
H.-S. P.
Wong
, “Verilog—A compact model for oxide-based resistive random access memory (RRAM),” in International Conference on Simulation of Semiconductor Processes and Devices (SISPAD) (IEEE, 2014), pp. 41–44.
39.
X.
Guan
,
S.
Yu
, and
H.-S. P.
Wong
, “
On the switching parameter variation of metal-oxide RRAM—Part I: Physical modeling and simulation methodology
,”
IEEE Trans. Electron Devices
59
,
1172
1182
(
2012
).
40.
L.
Chua
, “
Memristor—The missing circuit element
,”
IEEE Trans. Circuit Theory
18
,
507
519
(
1971
).
41.
S.
Amer
,
S.
Sayyaparaju
,
G. S.
Rose
,
K.
Beckmann
, and
N. C.
Cady
, “A practical hafnium-oxide memristor model suitable for circuit design and simulation,” in IEEE International Symposium on Circuits and Systems (ISCAS) (IEEE, 2017).
42.
M. D.
Pickett
,
D. B.
Strukov
,
J. L.
Borghetti
,
J. J.
Yang
,
G. S.
Snider
,
D. R.
Stewart
, and
R. S.
Williams
, “
Switching dynamics in titanium dioxide memristive devices
,”
J. Appl. Phys.
106
,
074508
(
2009
).
43.
S. R.
Bradley
,
A. L.
Shluger
, and
G.
Bersuker
, “
Electron-injection-assisted generation of oxygen vacancies in monoclinic HfO2
,”
Phys. Rev. Appl.
4
,
064008
(
2015
).
44.
S.
Aldana
,
P.
García-Fernández
,
A.
Rodríguez-Fernández
,
R.
Romero-Zaliz
,
M. B.
González
,
F.
Jiménez-Molinos
,
F.
Campabadal
,
F.
Gómez-Campos
, and
J. B.
Roldán
, “
A 3D kinetic Monte Carlo simulation study of resistive switching processes in Ni/HfO2/Si-n+-based RRAMs
,”
J. Phys. D: Appl. Phys.
50
,
335103
(
2017
).
45.
A.
Zeumault
,
S.
Alam
,
Z.
Wood
,
R. J.
Weiss
,
A.
Aziz
, and
G. S.
Rose
, “
TCAD modeling of resistive-switching of HfO2 memristors: Efficient device-circuit Co-design for neuromorphic systems
,”
Front. Nanotechnol.
3
,
734121
(
2021
).
46.
C.
Li
,
B.
Gao
,
Y.
Yao
,
X.
Guan
,
X.
Shen
,
Y.
Wang
,
P.
Huang
,
L.
Liu
,
X.
Liu
,
J.
Li
,
C.
Gu
,
J.
Kang
, and
R.
Yu
, “
Direct observations of nanofilament evolution in switching processes in HfO2-based resistive random access memory by in situ TEM studies
,”
Adv. Mater.
29
,
1602976
(
2017
).
47.
R.
Thamankar
,
N.
Raghavan
,
J.
Molina
,
F. M.
Puglisi
,
S. J.
O’Shea
,
K.
Shubhakar
,
L.
Larcher
,
P.
Pavan
,
A.
Padovani
, and
K. L.
Pey
, “
Single vacancy defect spectroscopy on HfO2 using random telegraph noise signals from scanning tunneling microscopy
,”
J. Appl. Phys.
119
,
084304
(
2016
).
48.
J.
McPherson
,
J.-Y.
Kim
,
A.
Shanware
, and
H.
Mogul
, “
Thermochemical description of dielectric breakdown in high dielectric constant materials
,”
Appl. Phys. Lett.
82
,
2121
2123
(
2003
).
49.
J. W.
McPherson
and
H. C.
Mogul
, “
Underlying physics of the thermochemical E model in describing low-field time-dependent dielectric breakdown in SiO2 thin films
,”
J. Appl. Phys.
84
,
1513
1523
(
1998
).
50.
F.
Jiménez-Molinos
,
F.
Gámiz
,
A.
Palma
,
P.
Cartujo
, and
J. A.
López-Villanueva
, “
Direct and trap-assisted elastic tunneling through ultrathin gate oxides
,”
J. Appl. Phys.
91
,
5116
5124
(
2002
).
51.
A.
Palma
,
A.
Godoy
,
J. A.
Jiménez-Tejada
,
J. E.
Carceller
, and
J. A.
López-Villanueva
, “
Quantum two-dimensional calculation of time constants of random telegraph signals in metal-oxide–semiconductor structures
,”
Phys. Rev. B
56
,
9565
9574
(
1997
).
52.
Synopsys, “TCAD Sentaurus,” (2019); see https://www.synopsys.com.
53.
S.
Aldana
,
P.
Garcia-Fernandez
,
R.
Romero-Zaliz
,
M. B.
Gonzalez
,
F.
Jimenez-Molinos
,
F.
Campabadal
,
F.
Gomez-Campos
, and
J. B.
Roldan
, “A kinetic Monte Carlo simulator to characterize resistive switching and charge conduction in Ni/HfO2/Si RRAMs,” in Proceedings of the 2018 12th Spanish Conference on Electron Devices (CDE), edited by J. Mateos and T. Gonzalez (IEEE, New York, 2018).
54.
L.
Larcher
,
A.
Padovani
,
O.
Pirrotta
,
L.
Vandelli
, and
G.
Bersuker
, “Microscopic understanding and modeling of HfO2 RRAM device physics,” in 2012 International Electron Devices Meeting (IEEE, 2012), pp. 20.1.1–20.1.4.
55.
J. G.
Simmons
, “
Generalized formula for the electric tunnel effect between similar electrodes separated by a thin insulating film
,”
J. Appl. Phys.
34
,
1793
1803
(
1963
).
56.
A.
Schenk
,
Advanced Physical Models for Silicon Device Simulation
(
Springer Science & Business Media
,
1998
).
57.
L.
Zhang
,
Y.-Y.
Hsu
,
F. T.
Chen
,
H.-Y.
Lee
,
Y.-S.
Chen
,
W.-S.
Chen
,
P.-Y.
Gu
,
W.-H.
Liu
,
S.-M.
Wang
,
C.-H.
Tsai
,
R.
Huang
, and
M.-J.
Tsai
, “
Experimental investigation of the reliability issue of RRAM based on high resistance state conduction
,”
Nanotechnology
22
,
254016
(
2011
).
58.
M. A.
Panzer
,
M.
Shandalov
,
J. A.
Rowlette
,
Y.
Oshima
,
Y. W.
Chen
,
P. C.
McIntyre
, and
K. E.
Goodson
, “
Thermal properties of ultrathin hafnium oxide gate dielectric films
,”
IEEE Electron Device Lett.
30
,
1269
1271
(
2009
).
59.
T.
Wang
and
J.
Roychowdhury
, “Well-posed models of memristive devices,” arXiv:1605.04897 (2016).
60.
P.
Lorenzi
,
R.
Rao
, and
F.
Irrera
, “
Forming kinetics in HfO2-based RRAM cells
,”
IEEE Trans. Electron Devices
60
,
438
443
(
2013
).
61.
T.
Cabout
,
J.
Buckley
,
C.
Cagli
,
V.
Jousseaume
,
J. F.
Nodin
,
B.
de Salvo
,
M.
Bocquet
, and
C.
Muller
, “
Role of Ti and Pt electrodes on resistance switching variability of HfO2-based resistive random access memory
,”
Thin Solid Films
533
,
19
23
(
2013
).
62.
C.
Jorel
,
C.
Vallée
,
E.
Gourvest
,
B.
Pelissier
,
M.
Kahn
,
M.
Bonvalot
, and
P.
Gonon
, “
Physicochemical and electrical characterizations of atomic layer deposition grown HfO2 on TiN and Pt for metal-insulator-metal application
,”
J. Vac. Sci. Technol. B
27
,
378
(
2009
).
63.
R.
Rofan
and
C.
Hu
, “
Stress-induced oxide leakage
,”
IEEE Electron Device Lett.
12
,
632
634
(
1991
).
64.
Y.
Chen
, “Reliability characterizations of ultra-thin gate oxides of MOSFETs,” Ph.D. thesis (University of Maryland, College Park, MD, 1998).
65.
C.
Cagli
,
J.
Buckley
,
V.
Jousseaume
,
T.
Cabout
,
A.
Salaun
,
H.
Grampeix
,
J. F.
Nodin
,
H.
Feldis
,
A.
Persico
,
J.
Cluzel
,
P.
Lorenzi
,
L.
Massari
,
R.
Rao
,
F.
Irrera
,
F.
Aussenac
,
C.
Carabasse
,
M.
Coue
,
P.
Calka
,
E.
Martinez
,
L.
Perniola
,
P.
Blaise
,
Z.
Fang
,
Y. H.
Yu
,
G.
Ghibaudo
,
D.
Deleruyelle
,
M.
Bocquet
,
C.
Müller
,
A.
Padovani
,
O.
Pirrotta
,
L.
Vandelli
,
L.
Larcher
,
G.
Reimbold
, and
B.
de Salvo
, “Experimental and theoretical study of electrode effects in HfO2 based RRAM,” in 2011 International Electron Devices Meeting (IEEE, 2011), pp. 28.7.1–28.7.4.
66.
L.
Lima
,
J.
Diniz
,
I.
Doi
, and
J.
Godoy Fo
, “
Titanium nitride as electrode for MOS technology and Schottky diode: Alternative extraction method of titanium nitride work function
,”
Microelectron. Eng.
92
,
86
90
(
2012
).
67.
M. N.
Solovan
,
V. V.
Brus
,
E. V.
Maistruk
, and
P. D.
Maryanchuk
, “
Electrical and optical properties of TiN thin films
,”
Inorg. Mater.
50
,
40
45
(
2014
).
68.
C. G. H.
Walker
,
J. A. D.
Matthew
,
C. A.
Anderson
, and
N. M. D.
Brown
, “
An estimate of the electron effective mass in titanium nitride using UPS and EELS
,”
Surf. Sci.
412–413
,
405
414
(
1998
).
69.
J. K.
Schaeffer
,
L. R. C.
Fonseca
,
S. B.
Samavedam
,
Y.
Liang
,
P. J.
Tobin
, and
B. E.
White
, “
Contributions to the effective work function of platinum on hafnium dioxide
,”
Appl. Phys. Lett.
85
,
1826
1828
(
2004
).
70.
G.
Fischer
,
H.
Hoffmann
, and
J.
Vancea
, “
Mean free path and density of conductance electrons in platinum determined by the size effect in extremely thin films
,”
Phys. Rev. B
22
,
6065
6073
(
1980
).
71.
H.
Hoffmann
and
G.
Fischer
, “
Electrical conductivity in thin and very thin platinum films
,”
Thin Solid Films
36
,
25
28
(
1976
).
72.
V. L.
Stout
and
M. D.
Gibbons
, “
Gettering of gas by titanium
,”
J. Appl. Phys.
26
,
1488
1492
(
1955
).
73.
S.
Venkataiah
,
S. J.
Chandra
,
U.
Chalapathi
,
C.
Ramana
, and
S.
Uthanna
, “
Oxygen partial pressure influenced stoichiometry, structural, electrical, and optical properties of DC reactive sputtered hafnium oxide films
,”
Surf. Interface Anal.
53
,
206
214
(
2021
).
74.
X.
Zhao
and
D.
Vanderbilt
, “
First-principles study of structural, vibrational, and lattice dielectric properties of hafnium oxide
,”
Phys. Rev. B
65
,
233106
(
2002
).
75.
J.
Wei
,
L.
Jiang
,
M.
Huang
,
Y.
Wu
, and
S.
Chen
, “
Intrinsic defect limit to the growth of orthorhombic HfO2 and (Hf, Zr)O2 with strong ferroelectricity: First-principles insights
,”
Adv. Funct. Mater.
31
,
2104913
(
2021
).
76.
A.
Aziz
,
X.
Li
,
N.
Shukla
,
S.
Datta
,
M.-F.
Chang
,
V.
Narayanan
, and
S. K.
Gupta
, “Low power current sense amplifier based on phase transition material,” in 2017 75th Annual Device Research Conference (DRC) (IEEE, 2017), pp. 1–2.
77.
M.-F.
Chang
,
S.-J.
Shen
,
C.-C.
Liu Chang
,
S.-J.
Shen
,
C.-C.
Liu
et al., “
An offset-tolerant fast-random-read current-sampling-based sense amplifier for small-cell-current nonvolatile memory
,”
IEEE J. Solid-State Circuits
48
,
864–877
(
2013
).