# Nonlinear elasticity under moderate to strong compression

B.L.N.. Kennett
Research School of Earth Sciences  The Australian National University
Canberra  ACT 2601  Australia
<EMAIL_ADDRESS>
###### Abstract

The strain-energy formulation of nonlinear elasticity can be extended to the
case of significant compression by modulating suitable strain energy terms by
a function of relative volume. For isotropic materials this can be
accomplished by the product of representations of shear, in terms of the
invariants of the Seth-Hill family of strain measures, and a function of
volume. The incremental shear modulus under pressure is determined by this
function, but nonlinear effects are retained for large strains. Suitable
functional forms can be derived from existing equations of state for moderate
to strong compression. For anisotropic materials, a similar development can be
made directly with strain energy terms depending directly on the Seth-Hill
strain tensors. Shear aspects can be emphasised by exploiting the
equivoluminal components of the strain tensors. Such formulations may be
helpful for materials under the conditions prevailing in the Earth’s interior.

###### keywords:

Compression, Shear Modulus, Strain Energy, Equations of State

##  1 Introduction

Many applications of nonlinear elasticity are concerned with extensional
environments with emphasis on shear properties, a useful review is provided by
Mihai & Goriely (2017). Compressibility has commonly been neglected in the
nonlinear case, but has been recognised to be significant in studies of soft
tissues (e.g. Beex, 2019). In contrast, in the study of the properties of
materials at high pressures the emphasis has been on the development of
equations of state for the bulk modulus. Improved experimental and
computational procedures mean that incremental shear properties from a
compressed state have become accessible, and so a full constitutive equation
is needed (Kennett, 2017). Large shears tend to be suppressed as pressure
increases, but can be significant in the Earth’s lithosphere.

Many of the formulations of shear properties are based on the superposition of
functions of members of the Seth-Hill strain tensors (Seth, 1964; Hill 1968)
and their associated conjugate stresses, which provide extensions of Hooke’s
law. The members of this suite of strain measures are characterised by the
exponent on the principal stretches. We here show how standard nonlinear
strain energy formulations can be adapted to carry shear properties into the
compressional regime, with the aid of an auxiliary function of density
modulating a deviatoric term.

For Earth materials, a semi-empirical linear relationship between the
incremental shear modulus, the bulk modulus and pressure can be used to
specify suitable functional forms for the auxiliary function. By this means a
shear modulus distribution can be associated with existing equations of state
to provide a full constitutive equation.. This strain energy formulation is
simplest in the isotropic case, but can be adapted to the anisotropic case by
using the full strain tensors rather than their invariants.

##  2 Isotropic materials under pressure

We consider a deformation from a reference state (unstressed) described by
coordinates $\boldsymbol{\xi}$ to a current state described by coordinates
$\mathbf{x}$. The relation between the states is provided by the deformation
gradient tensor $\mathbf{F}=\partial\mathbf{x}/\partial\boldsymbol{\xi}$, and
$J=\det\mathbf{F}=V/V_{0}$ is then the ratio of a volume element in the
current state ($V$) to that in the reference state ($V_{0}$). We introduce a
strain energy $W(\mathbf{F})$ depending on deformation, which specifies the
constitutive equation for a material.

In terms of $\mathbf{F}$ and the Green strain
$\mathbf{E}=\frac{1}{2}(\mathbf{F}^{T}\mathbf{F}-\mathbf{I})=\frac{1}{2}(\mathbf{C}-\mathbf{I})$,
the components of the stress tensor $\boldsymbol{\sigma}$ are given by

$J\sigma_{ij}=F_{ik}\frac{\partial W}{\partial
F_{jk}}=F_{ik}F_{jl}\frac{\partial W}{\partial E_{kl}},$ (1)

where we use the Einstein summation convention of summation over repeated
suffices.

The deformation gradient $\mathbf{F}$ can be written in terms of a stretching
component and a rotation in two ways

$\mathbf{F}=\mathbf{R}\mathbf{U}=\mathbf{V}\mathbf{R}$ (2)

where $\mathbf{U}^{2}=\mathbf{F}^{T}\mathbf{F}=\mathbf{C}$ and
$\mathbf{V}^{2}=\mathbf{F}\mathbf{F}^{T}=\mathbf{B}$. The matrices
$\mathbf{U}$, $\mathbf{V}$ have the same eigenvalues, the principal stretches
$\lambda_{1},\lambda_{2},\lambda_{3}$, but the principal axes vary in
orientation by the rotation $\mathbf{R}$.

The Seth-Hill class of strain measures take the form:

$\displaystyle\mathbf{E}_{q}(\mathbf{U})=\begin{cases}\frac{1}{q}(\mathbf{U}^{q}-\mathbf{I})&\text{if}\
q\neq 0,\\\ \ln{\mathbf{U}}&\text{if}\ q=0,\end{cases}$ (3)

where $\mathbf{I}$ is the identity tensor. The Green strain is thus
$\mathbf{E}_{2}$. All the members of this class of strain measures take the
same form for infinitesimal deformation.

The separation between volumetric deformation and shear-type deformation,
which is equivoluminal, can be achieved by working with $J$ and the normalised
deformation gradient $\mathbf{F}^{*}=J^{-1/3}\mathbf{F}$, so that
$\det\mathbf{F}^{*}=1$.

For an isotropic medium, the strain energy $W$ can be represented as a
function of invariants of strain measures (e.g. Spencer, 1980). Useful
invariants of $\mathbf{U}$, $\mathbf{V}$ are

$J=\lambda_{1}\lambda_{2}\lambda_{3}=\det\mathbf{U},$ (4)

a purely hydrostatic term, representing changes in volume, and the set

$L_{q}=J^{-q/3}[\lambda_{1}^{q}+\lambda_{2}^{q}+\lambda_{3}^{q}]=J^{-q/3}\mathrm{tr}\\{\mathbf{U}^{q}\\},\quad
q\neq 0,$ (5)

which concentrate on the deviatoric aspects of deformation. Note that
$\frac{1}{q}\\{L_{q}-3\\}$ corresponds to the trace of the equivoluminal part
of the Seith-Hill tensors, evaluated in terms of
$\mathbf{U}^{*}=J^{-1/3}\mathbf{U}$.

For an isotropic medium the principal axes of the stress tensor
$\boldsymbol{\sigma}$ align with those of $\mathbf{V}$, $\mathbf{B}$ (the
Eulerian triad), whereas the principal axes of $\mathbf{U}$, $\mathbf{C}$ and
$\mathbf{E}$ are rotated by $\mathbf{R}$ (the Lagrangian triad). In terms of
the principal stretches we can recast (1) in the form of an expression for the
$r$th principal stress

$\sigma_{r}=\frac{1}{J}\lambda_{r}\frac{\partial{W}}{\partial\lambda_{r}},\qquad\textrm{no
sum on }r,$ (6)

whilst recognising the rotation between the principal directions of the
elements on the left- and right-hand sides of the equation (6).

Many of the formulations for nonlinear shear given by Mihai & Goreiley (2017)
can be expressed as a linear combinations of the $L_{q}$ invariants, with
constant coefficients. Under compression Kennett (2017) has shown that it is
possible to associate a shear component to existing equations of state,
linking pressure and volume, by introducing a deviatoric term modulated by a
function of volume into the strain energy. The specific form used in Kennett
(2017) was derived from that for a neo-Hookean solid in terms of $L_{2}$, but
can be generalised to allow for a more complex shear behaviour.

Consider a strain energy function ${W}$ as a function of stretch invariants
$J$, $\\{L_{q}\\}$ with two independent volume terms $\Phi(J)$ and $\Psi(J)$:

${W}=\Phi(J)+\Psi(J)\sum_{q}a_{q}\frac{1}{q}\\{L_{q}-3\\},\quad\textrm{with}\
\ \sum_{q}a_{q}=1,$ (7)

incorporating a direct volume dependence in $\Phi(J)$ and a deviatoric
component in the second term. As noted above this is equivalent to an
expansion in terms of equivoluminal Seth-Hill tensors.
For purely hydrostatic compression:
$\lambda_{1}=\lambda_{2}=\lambda_{3}=\hat{\lambda}$, $J=\hat{\lambda}^{3}$ and
$\sum_{q}a_{q}\frac{1}{q}\\{L_{q}-3\\}=\sum_{q}a_{q}\frac{1}{q}\\{\hat{\lambda}^{-q}3\hat{\lambda}^{q}-3\\}=0$,
so that the deviatoric term
$\sum_{q}a_{q}\frac{1}{q}\\{L_{q}-3\\}\,\Psi(J)=0$.

For the strain energy (7) with both compressional and deviatoric components,
the $r$th principal stress takes the form:

$\sigma_{r}=\frac{\partial\Phi}{\partial J}+\frac{\partial\Psi}{\partial
J}\sum_{q}a_{q}\frac{1}{q}\\{L_{q}-3\\}+\frac{1}{J}\Psi(J)\sum_{q}a_{q}J^{-q/3}\\{\lambda_{r}^{q}-\textstyle{\frac{1}{3}}\displaystyle[\lambda_{1}^{q}+\lambda_{2}^{q}+\lambda_{3}^{q}]\\}.$
(8)

The full stress tensor $\boldsymbol{\sigma}$ can therefore be written as

$\boldsymbol{\sigma}=\mathbf{R}\left\\{\left[\frac{\partial\Phi}{\partial
J}+\frac{\partial\Psi}{\partial
J}\sum_{q}a_{q}\frac{1}{q}\\{L_{q}-3\\}\right]\mathbf{I}+\frac{1}{J}\Psi(J)\sum_{q}a_{q}J^{-q/3}\left\\{\mathbf{U}^{q}-\textstyle{\frac{1}{3}}\displaystyle\mathrm{tr}\\{\mathbf{U}^{q}\\}\mathbf{I}\right\\}\right\\}\mathbf{R}^{T}.$
(9)

For the purely hydrostatic case, the deviatoric terms vanish and the stress
tensor reduces to

$-p\mathbf{I}=\frac{\partial\Phi}{\partial J}\mathbf{I}.$ (10)

The incremental elastic moduli about this hydrostatically compressed state can
be extracted from the stress tensor (9) by making a first order expansion with
$\lambda_{r}=\hat{\lambda}(1+e_{r})$, so that
$J=\hat{\lambda}^{3}[1+\mathrm{tr}\\{\mathbf{e}\\}]+O(e^{2})$. In this case
the $r$th principal stress takes the form

$\sigma_{r}=-p+J\frac{\partial^{2}\Phi}{\partial
J^{2}}\,\mathrm{tr}\\{\mathbf{e}\\}+\frac{1}{J}\Psi(J)\left(\sum_{q}qa_{q}\right)[e_{r}-\textstyle{\frac{1}{3}}\displaystyle\mathrm{tr}\\{\mathbf{e}\\}].$
(11)

The representation of the principal stress in terms of the bulk modulus $K$
and shear modulus $G$ is

$\sigma_{r}=-p+K\mathrm{tr}\\{\mathbf{e}\\}+2G\left(e_{r}-\textstyle{\frac{1}{3}}\displaystyle\mathrm{tr}\\{\mathbf{e}\\}\right),$
(12)

and thus we identify the incremental moduli as:

$K=J\frac{\partial^{2}\Phi(J)}{\partial J^{2}},\qquad
G=\frac{1}{2J}\Psi(J)\sum_{q}qa_{q}.$ (13)

The shear properties for incremental strain are thus determined by $\Phi(J)$,
but for finite strain will be modulated by the nature of the sum over the
stretch invariants. Thus allows a wide variety of behaviour to be captured.
The choice of the functions of volume $\Phi(J)$ and $\Psi(J)$ depends on the
desired properties under pressure.

A more general representation of incremental properties about an initial
stress state in terms of the stretches $\\{\lambda_{i}$} has been provided by
Destrade & Ogden (2011), This treatment allows the possibility of non-
hydrostatic scenario, but reduces to (12) for a state of pure compression.

In applications to Earth materials, a number of different formulations have
been developed for equations of state through the strain energy term $\Phi(J)$
(see, e.g., Kennett, 2017). Formulations for shear are much less common, and
the most common form employed is the Birch-Murnaghan development in terms of
powers of Eulerian strain (Stixrude & Lithgow-Bertelloni, 2005). The
limitations of this approach for high pressures have been well documented by
Stacey & Davis (2004), who advocate instead the Keane equation of state for
bulk modulus, but this does not have an associated shear modulus. Kennett
(2017) has shown how the semi-empirical relation

$G=aK-bp,$ (14)

can be used to produce an effective representation of shear-properties under
pressure. Note that with the formulation above (.10,2.13) this means that
$\Psi(J)$ is related to the derivatives of $\Phi(J)$, with
$\partial\Phi/\partial J$ from pressure $p$ and $\partial^{2}\Phi/\partial
J^{2}$ from $K$. In terms of the bulk and shear moduli at zero pressure
($K_{0}$, $G_{0}$) and their pressure derivatives ($K^{\prime}_{0}$,
$G^{\prime}_{0}$)

$a=\frac{G_{0}}{K_{0}},\quad
b=\Big{(}\frac{G_{0}}{K_{0}}\Big{)}K_{0}^{\prime}-G_{0}^{\prime}.$ (15)

Equations (14, 15) provide a good representation of experimental results for
minerals, as illustrated in Figure 1 for the isotropic properties of MgO.

Figure 1: Illustration of the linear dependence of $G/K$ on $p/K$ to 20 GPa
pressure, for the adiabatic bulk modulus $K$ of periclase (MgO) using data
from Jackson & Niesler (1982), Sinogeikin & Bass (2000), and Zha et al.
(2000).

Figure 2: Illustration of the linear dependence of $G/K$ on $p/K$ for a
pyrolite composition lower mantle mineral assemblage from Gréaux et al.
(2019). The open symbols indicate the zone where some residual garnet may be
present.

The linear relation also provides a good description of the properties of
mineral assemblages. We illustrate the results for the Earth’s lower mantle
using the model developed by Gréaux et al. (2019) in Figure 2. The dominant
minerals are bridgmanite and ferropericlase, and some residual majorite garnet
is present at the top of the lower mantle where there is a slight deviation
from the linear trend. Although the linear form works well over a large range
of pressures (up to 140 GPa for the lower mantle), Burakovsky et al. (2004)
suggest that (14) should be modified with a slowly-varying pressure dependence
for $b$ to allow a match to the expectation for infinite pressure.

The combination of the strain energy development (7) with the identification
of moduli (13) and the relation (14) provides a flexible way of extending
nonlinear elastic effects to a compressed state, whilst retaining complex
shear behaviour for finite strains.

##  3 Anisotropic materials under pressure

Many natural materials such as wood and tissues show distinct anisotropy in
their properties. Most minerals have significant anisotropy, and its only in
aggregate that much of the Earth appears to have nearly isotropic properties.
For the description of the behaviour of materials under strong compression it
is therefore desirable to be able to provide a full description of the
anisotropic behaviour.

In the general anisotropic situation the directions of the principal stretches
do not remain constant and so their variation has to be taken into account in
any formulation. Even so it is possible to express the strain energy $W$ in
terms of $J$ and the normalised deformation gradient $\mathbf{F}^{*}$ as
$W(\mathbf{F})=W^{*}(J,\mathbf{F}^{*})$. With this representation the Cauchy
stress tensor $\boldsymbol{\sigma}$ is given by

$\boldsymbol{\sigma}=\frac{1}{J}\mathbf{F}\frac{\partial
W}{\partial\mathbf{F}}=\frac{\partial W^{*}}{\partial
J}\mathbf{I}+\frac{1}{J}\left[\mathbf{F}^{*}\frac{\partial
W^{*}}{\partial\mathbf{F}^{*}}-\textstyle{\frac{1}{3}}\displaystyle\mathrm{tr}\left\\{\mathbf{F}^{*}\frac{\partial
W^{*}}{\partial\mathbf{F}^{*}}\right\\}\mathbf{I}\right].$ (16)

Hence in asimilar way to the isotropic case above we can make a separation
into pressure dependence and shear deformations.

For orthotropic materials, Latorre & Montáns (2017) have split the strain
energy into an isotropic and a specifically anisotropic part. Such an approach
may well be suitable for small compressions, but when we want to include
strong compression we should allow for volume dependence of all the
components. Building on the approach used for isotropy we look to combine a
volumetric component with equivoluminal term modulated by a function of
relative volume. We use the equivoluminal equivalents of the Seth-Hill
measures (e.g., Miehe & Lambrecht, 2001)

$\displaystyle{\mathbf{E}}^{*}_{q}(\mathbf{U}^{*})=\begin{cases}\frac{1}{q}(\mathbf{U}^{*q}-\mathbf{I})&\text{if}\
q\neq 0,\\\ \ln{\mathbf{U}^{*}}&\text{if}\ q=0.\end{cases}$ (17)

and construct a strain energy function

$W^{*}(J,\mathbf{F}^{*})=\Phi(J)+\sum_{q}S_{q}(J)W_{q}(\mathbf{E}^{*}_{q}).$
(18)

Then in terms of the equivoluminal strain
$\mathbf{E}^{*}=\mathbf{E}_{2}^{*}=\textstyle{\frac{1}{2}}\displaystyle(\mathbf{U}^{*2}-\mathbf{I})$,

$\mathbf{F}^{*}\frac{\partial
W^{*}}{\partial\mathbf{F}^{*}}=\mathbf{F}^{*}\frac{\partial}{\partial\mathbf{E}^{*}}\sum_{q}S_{q}(J)W_{q}(\mathbf{E}^{*}_{q})\mathbf{F}^{*T}=\mathbf{F}^{*}\sum_{q}\frac{\partial
S_{q}(J)}{\partial\mathbf{F}^{*}}W_{q}(\mathbf{E}^{*}_{q})+\mathbf{F}^{*}\sum_{q}S_{q}(J)\frac{\partial
W_{q}(\mathbf{E}^{*}_{q})}{\partial\mathbf{E}^{*}}.$ (19)

The derivative of the functions of relative volume

$\frac{\partial S_{q}(J)}{\partial\mathbf{F}^{*}}=\frac{\partial
S_{q}(J)}{\partial J}\frac{\partial J}{\partial\mathbf{F}^{*}}=J\frac{\partial
S_{q}(J)}{\partial J}\mathbf{F}^{*-T}.$ (20)

Thus, the shear component of the Cauchy stress expression

$\mathbf{F}^{*}\frac{\partial
W^{*}}{\partial\mathbf{F}^{*}}=J\mathbf{I}\sum_{q}\frac{\partial
S_{q}(J)}{\partial
J}W_{q}(\mathbf{E}^{*}_{q})+\mathbf{F}^{*}\sum_{q}S_{q}(J)\frac{\partial
W_{q}(\mathbf{E}^{*}_{q})}{\partial\mathbf{E}_{q}^{*}}\boldsymbol{\mathsf{P}}_{q}.$
(21)

The fourth order projection tensor
$\boldsymbol{\mathsf{P}}_{q}=\partial\mathbf{E}^{*}_{q}/\partial\mathbf{E}^{*}$
is detailed in the Appendix.

The contributions to the stress (16) from the functions $\\{S_{q}(J)\\}$ are
purely hydrostatic, and the shear dependence comes from the choices made for
$\\{W_{q}(\mathbf{E}^{*})\\}$. The evolution of the stress tensor under
pressure and the consequent elastic properties can be evaluated by a
perturbation treatment around a state of pure compression as in Section 2. As
before the bulk modulus $K$ is given by $J\partial^{2}\Phi/\partial J^{2}$.

A simple form for the individual strain energy terms is quadratic:

$W_{q}=\textstyle{\frac{1}{2}}\displaystyle\mathbf{E}_{q}:\boldsymbol{\mathsf{Z}}_{q}:\mathbf{E}_{q},$
(22)

where $\boldsymbol{\mathsf{Z}}_{q}$ is a fourth-order stiffness tensor with 21
independent components, and $:$ denotes the double inner product, so that
$\mathbf{A}:\mathbf{B}=A_{ij}B_{ji}$. With a sum of a number of Seth-Hill
contributions a variety of deformation styles can be produced (e.g., Beex,
2019). In this case for hydrostatic stress $\mathbf{E}^{*}_{q}$ vanishes, and
as in the treatment of isotropic elasticity in Section 2 a perturbation
treatment about a hydrostatic state simplifies significantly to leave a shear
contribution specified by the $S_{q}(J)$.

In this anisotropic development we have introduced separate functions of
relative volume for each order of the Seth-Hill strain measures, but
simplified forms may be preferable. If the anisotropic properties are
consistent with increasing pressure, a suitable strain energy formulation for
a material under moderate compression would be

$W=\Phi(J)+S_{A}(J)\left[\mathbf{E}^{*}_{-1}:\boldsymbol{\mathsf{A}}:\mathbf{E}^{*}_{-1}\right],$
(23)

in terms of the equivoluminal component of the Almansi strain
${\mathbf{E}}^{*}_{-1}=\mathbf{I}-J^{1/3}\mathbf{U}^{-1}$ (Seth-Hill element
of order -1) with a purely volumetric function $S_{A}(J)$. The fourth order
tensor $\boldsymbol{\mathsf{A}}$, specifies the shear properties at the rest
state $J=1$. The choice of $\Phi(J)$ can be taken from formulations for
equation of state, and then $S_{A}(J)$ can be associated in a similar way to
the treatment of the shear modulus above.

For materials such as MgO whose anisotropy varies strongly with increasing
pressure (Karki et al., 1997) we need to add an additional term to the strain
energy, e.g.,

$S_{B}(J)\left[\mathbf{E}^{*}_{-2}:\boldsymbol{\mathsf{B}}:\mathbf{E}^{*}_{-2}\right],$
(24)

in terms of the equivoluminal Eulerian strain
${\mathbf{E}}^{*}_{-2}=\mathbf{I}-J^{2/3}\mathbf{U}^{-2}$ (Seth-Hill element
of order -2) . The function $S_{B}(J)=0$ at $J=1$, and can be tuned to
represent the variations in anisotropy with pressure.

##  4 Conclusion

We have shown how it is possible to develop formulations of nonlinear
elasticity that can accommodate large shear and high compression, by the
introduction of a shear function as a function of volume modulating a
deviatoric term. For Earth materials, the functional dependence of the shear
properties can be guided by the semi-empirical linear relation between shear
modulus, bulk modulus and pressure.

For anisotropy a similar development can be made with functions of volume
combined with strain energies depending on the the equivoluminal components of
the Seth-Hill family of strain tensors. The flexibility of the development
provides a means of representing a wide range of isotropic and anisotropic
scenarios suitbale for conditions in the Earth’s interior.

The formulation developed in this work has been oriented toward situations
with moderate to strong compression, but could also be used for strong
expansion with a switch in the style of strain measures employed. For
compression, the deviatoric component is best represented using measures
depending on strain exponent $q<0$, but in tension $q>0$ is to be preferred
(Beex, 2019).

## Appendix

The normalised stretch tensor $\mathbf{U}^{*}$ can be written in terms of its
eigenvalues, the normalised stretches $\lambda_{i}^{*}=J^{-1/3}\lambda_{i}$,
and their associated orthogonal eigenvectors $\mathbf{n}_{i}$ as:

$\mathbf{U}^{*}=\sum_{i=1}^{3}\lambda_{i}^{*}\mathbf{n}_{i}\mathbf{n}_{i},$
(A.25)

in terms of the dyadic product of the eigenvectors. The projection tensors
$\boldsymbol{\mathsf{P}}_{q}$ introduced in (19) depend on the evolution of
strain (Miehe & Lambrecht, 2001; Beex, 2019) and can also be written in terms
of the eigen-quantities:

$\boldsymbol{\mathsf{P}}_{q}=2\frac{\partial\mathbf{E}^{*}_{q}}{\partial\mathbf{E}}=\sum_{i=1}^{3}d^{\\{q\\}}_{i}\mathbf{n}_{i}\mathbf{n}_{i}\mathbf{n}_{i}\mathbf{n}_{i}+\sum_{r=1}^{3}\sum_{j\neq
i}^{3}\vartheta^{\\{q\\}}_{ij}\big{(}\mathbf{n}_{i}\mathbf{n}_{j}\mathbf{n}_{i}\mathbf{n}_{j}+\mathbf{n}_{i}\mathbf{n}_{j}\mathbf{n}_{j}\mathbf{n}_{i}\big{)}.$
(A.26)

The coefficients $d_{i}$ depend on the stretches and the order of the strain
element $q$

$d^{\\{q\\}}_{i}=\lambda_{i}^{*(q-2)}.$ (A.27)

For three distinct stretches,

$\vartheta^{\\{q\\}}_{ij}=\frac{1}{q}\frac{\lambda_{i}^{*q}-\lambda_{j}^{*q}}{\lambda_{i}^{*2}-\lambda_{j}^{*2}}.$
(A.28)

When two stretches are equal
$\lambda^{*}_{a}=\lambda^{*}_{b}\neq\lambda^{*}_{c}$,

$\vartheta^{\\{q\\}}_{ab}=\textstyle{\frac{1}{2}}\displaystyle d_{a}.$ (A.29)

For the hydrostatic case, $\lambda^{*}_{1}=\lambda^{*}_{2}=\lambda^{*}_{3}=1$,
the coefficient $\vartheta^{\\{q\\}}_{ij}=\frac{1}{2}$.

Second derivative projection operators can be defined in a similar way in
terms of the stretches and their associated eigenvectors, but now involve
sixth-order tensors (Miehe & Lambrecht, 2001; Beex, 2019).

## References

  * [1]
  * [2] eex L.A.A. 2019. Fusing the Seth–Hill strain tensors to fit compressible elastic material responses in the nonlinear regime. Int. J. Mech. Sci. 163 105072
  * [3]
  * [4] urakovsky L., Preston D.L., Wang Y., 2004. Cold shear modulus and Grüneisen parameter at all densities, Solid State Commun. 132 151–156.
  * [5]
  * [6] estrade M., Odgen R.W., 2013. On stress-dependent elastic moduli and wave speeds, IMA J. Applied Mathematics, 78, 965–977.
  * [7]
  * [8] réaux S., Irifune T., Higo Y., Tange Y., Arimoto T., Liu Z., Yamada A., 2019. Sound velocity of CaSiO3 perovskite suggests the presence of basaltic crust in the Earth’s lower mantle. Nature 565 218-221.
  * [9]
  * [10] ill R., 1968. On constitutive inequalities for simple materials—I. J. Mech. Phys. Solids 16 229–242.
  * [11]
  * [12] ackson I., Niesler H., 1982. The elasticity of periclase to 3GPa and some geophysical implications. In: Akimoto, S., Manghnani, M.H. (Eds.), High-pressure Research in Geophysics, pp. 93–113, Centre of Academic Publications, Japan.
  * [13]
  * [14] arki B.B., Stixrude L., Clark S.J., Warren M.C., Ackland G.J., Crain J., 1997. Structure and elasticity of MgO at high pressure. American Mineralogist, 82, 51–-60.
  * [15]
  * [16] ennett B.L.N., 2017. Towards constitutive equations for the deep Earth. Phys. Earth Planet. Inter. 270, 40–45.
  * [17]
  * [18] atorre M., Montáns F.J., 2017. WYPIWYG hyperelasticity without inversion formula: application to passive ventricular myocardium. Comput. Struct. 185 47–-58.
  * [19]
  * [20] iehe C., Lambrecht M., 2001. Algorithms for computation of stresses and elasticity moduli in terms of Seth-Hill’s family of generalized strain tensors. Commun. Numer. Methods Eng. 17, 337–353.
  * [21]
  * [22] ihai L.A., Goriely A., 2017. How to characterize a nonlinear elastic material? A review on nonlinear constitutive parameters in isotropic finite elasticity. Proc. R. Soc. A 473 20170607\.
  * [23]
  * [24] eth B.R., 1964. Generalized strain measure with application to physical problems. In: Second-Order Effects in Elasticity, Plasticity and Fluid Dynamics, Reiner M, Abir D (eds). Pergamon Press: Oxford, 162–172.
  * [25]
  * [26] inogeikin S.V., Bass J.D., 2000. Single-crystal elasticity of pyrope and MgO to 20 GPa by Brillouin scattering in the diamond cell. Phys. Earth Planet. Int. 120 43–-62.
  * [27]
  * [28] pencer A.J.M., 1980. Continuum Mechanics, Longman.
  * [29]
  * [30] tacey F.D., Davis P.M., 2004. High pressure equations of state with applications to the lower mantle and core. Phys. Earth Planet. Inter. 142 137–184.
  * [31]
  * [32] tixrude L., Lithgow-Bertelloni C., 2005. Thermodynamics of mantle minerals - I. Physical Properties. Geophys. J. Int., 162, 610–632.
  * [33]
  * [34] ha C.-S., Mao H.-K., Hemley R.J., 2000. Elasticity of MgO and a primary pressure scale to 55 GPa. PNAS 97 13494–-13499.


Contributed equally to this work Contributed equally to this work

# Dynamic Behaviors and Training Effects in TiN/Ti/HfOx/TiN Nanolayered
Memristors with Controllable Quantized Conductance States: Implications for
Quantum and Neuromorphic Computing Devices

Min-Hsuan Peng  Department of Physics, National Taiwan Normal University,
Taipei 116, Taiwan  Ching-Yang Pan  Department of Physics, National Taiwan
Normal University, Taipei 116, Taiwan  Hao-Xuan Zheng  Department of Physics,
National Sun Yat-Sen University, Kaohsiung 804, Taiwan  Ting-Chang Chang
Department of Physics, National Sun Yat-Sen University, Kaohsiung 804, Taiwan
Pei-hsun Jiang  Department of Physics, National Taiwan Normal University,
Taipei 116, Taiwan<EMAIL_ADDRESS>
###### Abstract

Controllable quantized conductance states of TiN/Ti/HfOx/TiN memristors are
realized with great precision through a pulse-mode reset procedure, assisted
with analytical differentiation of the condition of the set procedure, which
involves critical monitoring of the measured bias voltage. An intriguing
training effect that leads to faster switching of the states is also observed
during the operation. Detailed analyses on the low- and high-resistance states
under different compliance currents reveal a complete picture of the
structural evolution and dynamic behaviors of the conductive filament in the
HfOx layer. This study provides a closer inspection on the quantum-level
manipulation of nanoscale atomic configurations in the memristors, which helps
to develop essential knowledge about the design and fabrication of the future
memristor-based quantum devices and neuromorphic computing devices.

Keywords: HfO2, filament, resistive random-access memory (RRAM), memristor,
oxygen vacancy, resistive switching, conductance quantization, training effect

##  1 Introduction

Memristors with high scalability, low power consumption, and multilevel
switching is one of the promising candidates for artificial synapses in
neuromorphic computing to replace the conventional von-Neumann architecture 1,
2, 3. Linear conductance of a memristor in a wide voltage range is pursued for
the purpose of implementing vector-matrix multiplication in conductance
programming for memristor arrays 4, 5, 6. Nonlinear memristor dynamics due to
intrinsic conduction mechanisms 7, 8, 9 are therefore one of the key
challenges to build memristor-based dot-product engines.

In the mean time, quantized conduction in memristor-based devices is being
introduced to this field for its great potential for high-density data storage
through multilevel switching, and for analog synaptic weight update in
effective training of the artificial neural networks 10, 11. While
implementations of quantum neuromorphic computing platforms with quantum
memristors are being proposed 12, 13, realization of these architectures
remains difficult because conductance quantization of the memristors suffers
significant instability in terms of endurance and tuning accuracy, which
includes large half-widths in the histogram of the quantized conductance 14,
15. The occurrence of conductance quantization seems unstable and random even
when an optimal condition is used in the measurement 16, 17, 18. The mechanism
that guarantees conductance quantization in a memristor remains a mystery. The
detailed atomic dynamics of the conductive filaments in the memristors is not
yet fully explored and therefore demands more research.

In this letter, we have made in-depth investigations on conductance
characteristics of bipolar TiN/Ti/HfOx/TiN valence-change memristors (VCMs),
aiming to look into the instability issue of the quantized conductance. The
dynamics of the set procedure is observed to be decisive for the quantization
performance in the reset procedure. Detailed analyses on the low-resistance
state (LRS) have also been conducted to explicitly explore the electrical
characteristics associated with the nanoscale atomic structure of the
conductive filament in the HfOx layer. With better understanding of the atomic
dynamic behaviors of the conductive filament, we are able to perform a precise
control of the quantized conductance states of the memristors.

##  2 Device Fabrication and Measurement Methods

Fig. 1 shows the device layout of the TiN/Ti/HfOx/TiN memristor. A 300-nm SiO2
layer was grown via wet oxidation on a lightly doped p-type Si(100) substrate.
The SiO2 layer serves as an insulating layer between the Si substrate and the
bottom electrodes, which were formed by depositing TiN (50 nm)/Ti (50 nm)
using radio-frequency (rf) sputtering. After that, a 10-nm switching layer of
HfO2 was formed with atomic layer deposition, followed by a deposition of TiN
(40 nm)/Ti (10 nm) layer as top electrodes. Lithography and inductively-
coupled-plasma (ICP) etching were then used to define the active cell area.
Then, a low-temperature oxide (LTO) SiO2 layer was deposited, followed by
lithography and ICP etching to create via-hole structures in LTO. Finally,
AlCu/TaN contacts for the electrodes were formed via lithography and rf
sputtering. The electrical characteristics are measured at room temperature
using Keithley 2400 and Agilent B1500A. The bias voltages are applied to the
top electrodes as the bottom electrodes are grounded during electrical
measurements.

Figure 1: (a) Schematic drawing (not to scale) and (b) the TEM image of the
cross-sectional view of the TiN/Ti/HfOx/TiN memristor.

##  3 Results and Discussion

###  3.1 Measurements with the DC Voltage Sweep Mode

Figure 2: Currents in a log scale as functions of voltage with the compliance
current for the set procedure $I_{\mathrm{c}}=135$ $\upmu$A (blue curve) and
$60$ $\upmu$A (red curve), respectively. A stepwise feature is observed during
the reset of the curve with $I_{\mathrm{c}}=60$ $\upmu$A, and its blown-up
view in a linear scale is shown in the inset. The stepwise feature is absent
from the curve with $I_{\mathrm{c}}=135$ $\upmu$A.

Electrical measurements are performed on several devices with the same
structure as described in Section 2, and the results are found to be similar
and reproducible. The data presented in this paper are measured from a device
with a cell area of 0.36 $\upmu$m2. The forming procedure before the
electrical measurements for each device is described in detail in Supporting
Information Section LABEL:sec:forming. Two representative current-vs.-voltage
($I$–$V$) curves in the dc voltage sweep mode of the memristor operation are
shown in Fig. 2, one with signatures of conductance quantization in the reset
procedure and one without. A compliance current ($I_{\mathrm{c}}$) is applied
during each set procedure to prevent permanent breakdown. The $I$–$V$ curve
with $I_{\mathrm{c}}=135$ $\upmu$A (blue curve) shows the standard electrical
characteristics with a steep drop in $I$ at $-0.85$ V in the reset procedure
after the current reaches a maximum, whereas the one with $I_{\mathrm{c}}=60$
$\upmu$A (red curve) exhibits a series of descending steps (boxed with dashed
lines) starting at a smaller bias voltage of $-0.38$ V. The steps correspond
to the quantized conductance of a conducting channel in the switching layer,
which reveals the nanoscale atomic-level reaction of a conductive filament
consisting of oxygen vacancies in the switching layer 14, 19, 20, 21, 22.
During the reset procedure as the current is slowly switched from LRS to the
high-resistance state (HRS), the quantum point contact of the filament that
touches the negatively-charged top Ti electrode thins further and further
because the oxygen vacancies of the filament are gradually removed under the
voltage stress through recombination of the oxygen vacancies of the filament
and oxygen ions from the Ti layer 21. After the last oxygen vacancy in contact
is removed, breaking the circuit established by the filament, a Schottky
barrier is created in the conduction 23.

Figure 3: Representative examples of conductance quantization during the reset
procedure with $I_{\mathrm{c}}=60$ $\upmu$A. Conductance plateaus occur at
half-integer multiples of $2e^{2}/h$.

The corresponding conductance of the quantum point contact in the reset
procedure of the $I_{\mathrm{c}}=60$ $\upmu$A curve in Fig. 2 is calculated
and expressed in terms of the conductance quantum $G_{0}=2e^{2}/h$ in Fig.
3(a), along with other examples of conductance quantization of the device
shown in Figs. 3(b)–3(d). In each set-and-reset cycle, the voltage is swept at
a rate of $\Delta V=\pm 10$ mV per 0.5 second for each data point, except for
the reset procedure from $-0.35$ V to $-1$ V , during which the sweep rate is
decreased to $\Delta V=-2$ mV per 0.5 second to gently process the switching
of the quantized conductance states. The resistance in series with the atomic
point contact must be taken into consideration to precisely extract the
quantized conductance of the point contact 19, 24. It can be seen that the
conductance plateaus occur at some of the half-integer multiples of $G_{0}$.
Fig. 4 summarizes the numbers of counts of respective values of the
conductance plateaus collected from 330 $I$–$V$ curves using various
$I_{\mathrm{c}}$ (listed in Table 1). The histogram clearly demonstrates the
tendency of the device to yield quantized conductance. The appearance of half-
integer multiples of $G_{0}$ instead of merely integer multiples is not
universal in atomic point contacts. It is suggested to be caused by the
chemical potential difference between the two carrier reservoirs across the
filament 25, rearrangement of the atomic contact configuration 26, or possible
weak magnetism from oxygen vacancies that may lift spin degeneracy 14.

Figure 4: Histogram of the values of conductance plateaus retraced from 330
reset operations using the dc voltage sweep mode. More information about these
330 operations is listed in Table 1.

The compliance current $I_{\mathrm{c}}$ for the set procedure plays an
important role in search of the quantized conductance states of a memristor.
The set procedures have been performed with different $I_{\mathrm{c}}$ from
165 to 40 $\upmu$A, with 30 set-and-reset cycles completed for each
$I_{\mathrm{c}}$. The statistics of the results from a memristor with a cell
area of 0.36 $\upmu$m2 are listed in Table 1. (Statistics from other devices
with different cell areas show the similar behaviors; see Supporting
Information Section LABEL:sec:cell. Temperature dependence of the electrical
characteristics is also studied, as shown in Supporting Information Section
LABEL:sec:temp.) All the curves that exhibit conductance quantization in the
reset procedure (see the row “w/ Quant.”) belong to the fair-set group
(definition in the next paragraph). The chance of observing conductance
quantization in the reset procedure stays zero for $I_{\mathrm{c}}=165$
$\upmu$A to 105 $\upmu$A, and then gradually increases to $50\%$ as
$I_{\mathrm{c}}$ is gradually decreased to 70 $\upmu$A, and then reaches the
maximum $67\%$ when $I_{\mathrm{c}}=60$ $\upmu$A. The percentage then
decreases to $33\%$ as $I_{\mathrm{c}}$ is decreased further to 40 $\upmu$A.
This is an $I_{\mathrm{c}}$ so small that a good set (definition in the next
paragraph) can barely be acquired. With the highest yield of conductance
quantization, $I_{\mathrm{c}}=60$ $\upmu$A is considered to be the optimal
condition for later operation of the memristor for controlling the quantized
conductance states.

$I_{\mathrm{c}}$ ($\upmu$A) | 165 | 150 | 135 | 120 | 105 | 90 | 80 | 70 | 60 | 50 | 40
---|---|---|---|---|---|---|---|---|---|---|---
a) Good set | 30 | 30 | 30 | 30 | 30 | 23 | 20 | 11 | 5 | 4 | 0
w/ SCLC | 19 | 18 | 20 | 17 | 20 | 11 | 7 | 4 | 3 | 0 | 0
b) Fair set | 0 | 0 | 0 | 0 | 0 | 7 | 10 | 18 | 21 | 15 | 17
w/ Quant. | 0 | 0 | 0 | 0 | 0 | 6 | 10 | 15 | 20 | 12 | 10
c) Poor set | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 4 | 8 | 10
d) Set failure | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3
Table 1: Numbers of counts of various set conditions using different
$I_{\mathrm{c}}$. (a) Good sets. Bottom row: good sets with SCLC in LRS. (b)
Fair sets. Bottom row: fair sets with quantized conductance in the reset
procedure. Notice that the fair-set category is the only category in which
current quantization in the reset procedure can be observed. (c) Poor sets.
(d) Set failures.

The electrical characteristics of the set-and-reset cycles can be generally
classified into five categories, as illustrated in Fig. 5 with representative
examples. (More examples are presented in Supporting Information Section
LABEL:sec:set.) The blue dashed curves are plotted against the programed bias
voltage provided by the voltage source, whereas the red solid curves are
plotted against the measured bias voltage ($V_{\mathrm{m}}$). The only
discrepancy between them lies in the set procedure when the resulting current
abruptly jumps high to hit $I_{\mathrm{c}}$. The five categories are as
follows:

  1. 1.

Good set with space-charge-limited current (SCLC) (Figs. 5(a) and 5(b)): The
conduction in LRS after the set procedure is ohmic with a resistance of
2.0–4.0 k$\Omega$ (mostly around $3$ k$\Omega$), followed by a significant
slope increase (boxed with dashed lines) at a negative voltage ($-0.48$ V in
the representative example), known as the SCLC feature. The current then
gradually increases to a maximum before it drops abruptly into HRS in a reset
procedure (at $-$$0.85$ V in Fig. 5(a) and $-$$0.61$ V in Fig. 5(b)). A good
set procedure is featured with an $I$–$V_{\mathrm{m}}$ curve hitting
$I_{\mathrm{c}}$ at only one point for a stay, indicating a very stable
$V_{\mathrm{m}}$.

  2. 2.

Good set without SCLC (Fig. 5(c)): This has similar features with the previous
category, except that it tends to undergo the reset process at smaller
voltages ($\sim$0.18 V smaller on average), and the SCLC signature is missing.
(At the spot in the dashed box, it seems almost entering the SCLC regime
especially when compared with Fig. 5(b), but the filament fails to hold for
it. More discussion about the SCLC is presented in the following contexts and
in Supporting Information Section LABEL:sec:goodset.) Tiny fluctuations are
observed at larger negative bias voltages in LRS, before the reset procedure
takes place (at $-$$0.59$ V in this example).

  3. 3.

Fair set (Fig. 5(d)): The conduction in LRS after the set is ohmic with a
resistance of 4.8–7.5 k$\Omega$ (mostly around $6$ k$\Omega$), which is about
twice larger than those in the good-set cases. Tiny fluctuations are usually
observed at larger bias voltages in the ohmic state. The reset procedure
starts at an even smaller bias voltage ($-$$0.38$ V in this example), but with
a progressive phase that is prolonged to a more negative bias voltage. A fair
set procedure is featured with an $I$–$V_{\mathrm{m}}$ trace frequently
wiggling left and right along the horizontal compliance line. This is the only
category in which current quantization in the reset procedure can be observed.
Notice that $I_{\mathrm{c}}=60$ $\upmu$A in both Figs. 5(c) and 5(d). In very
few cases where quantization is absent (not presented in Fig. 5), the reset
process exhibits chaotic noise-like fluctuations similar to those found in
Fig. 5(e).

  4. 4.

Poor set (Fig. 5(e)): The conduction in LRS is no longer ohmic, but follows
the Schottky-emission equation, with resistance considerably higher than that
in the fair-set cases. A poor-set procedure is also featured with an
$I$–$V_{\mathrm{m}}$ trace wiggling along the compliance ceiling, and even
occasionally falling off and rising back to the ceiling before entering LRS.

  5. 5.

Set failure (Fig. 5(f)): The $I$–$V$ curve cannot enter LRS after the set
procedure.

Figure 5: The electrical characteristics of the set-and-reset cycles can be
classified into five categories: (a)(b) good set with SCLC, (c) good set
without SCLC, (d) fair set, (e) poor set, and (f) set failure. The blue dashed
curves are plotted against the programed bias voltage provided by the voltage
source, whereas the red solid curves are plotted against the measured bias
voltage.

The Schottky-emission equation is shown as follows:

${I}=aA^{*}{T}^{2}\exp\left[\frac{-(\phi_{\mathrm{B}}-\sqrt{{q^{3}V}/{4\pi\epsilon
d_{\mathrm{s}}}})}{kT}\right],$

where $a$ is the effective cross section of the filament, $A^{*}$ is the
effective Richardson constant, $q$ is the carrier charge, $T$ is the
temperature, $\phi_{\mathrm{B}}$ is the energy barrier height,
$d_{\mathrm{s}}$ is the effective switching thickness, and $\epsilon$ is the
permittivity, which is $\sim$$25\epsilon_{0}$ for HfO2. A fit to the LRS curve
in Fig. 5(e) yields $\phi_{\mathrm{B}}=0.36$ eV and $d_{\mathrm{s}}=2.4$ nm,
which reveals the characteristics of the device structure with the filament
grown between the bottom TiN electrode and the positively-biased top Ti
electrode to a point where it is only a tiny gap ($\sim$$2.4$ nm) away from
completion of the connection. The value of $\phi_{\mathrm{B}}$ is about half
of that of our previous similar HfOx memristors 23, which may imply a larger
amount of impurities or defects in the switching layer.

Shown in the inset of Fig. 5(a) is a log-scale blown-up view of the spot where
the $I$–$V$ slope changes due to SCLC, which is known to be mostly detected in
memristors with $\phi_{\mathrm{B}}\lesssim 0.3$ eV 27. The $I$–$V$
characteristic is ohmic until $V$ is swept to a more negative value of $-0.48$
V, where the device enters the trap-filling regime 28, 29, 30 with the slope
in the log scale prominently increased to $>2$, until the $I$–$V$
characteristic changes again to follow the Mott–Gurney law of SCLC 31 starting
at $V=-0.52$ V:

$I=\frac{9}{8}a\epsilon\mu\frac{V^{2}}{d_{s}^{3}},$

where $\mu$ is the electron mobility, which is $\sim$200 cm2/Vs in HfO2 32,
33. The $V^{2}$ law of SCLC only holds for a limited range of bias voltage.
After $V$ is swept to $-0.64$ V, the characteristic becomes $I\propto V^{x}$
with $1<x<2$, until the reset procedure starts. This may be interpreted with a
negative field dependence of the mobility of the space charges due to
positional and energetic disorder in the material 34, 35. At lower electric
fields, the most energetically favorable paths for percolative hopping
transport will proceed via randomly oriented jumps. However, with increasing
electric field, charge carriers are forced to make less energetically
favorable jumps in the direction of the field, leading to a reduced mobility.

It is generally believed that a larger $I_{\mathrm{c}}$ leads to a more
compact structure of the oxygen vacancies 36 or a larger diameter 37 of the
conductive filament. Therefore, a large $I_{\mathrm{c}}$ in our experiment
such as 135 $\upmu$A results in an $I$–$V$ curve with standard characteristics
of a fairly strong filament, as shown in Fig. 5(a). For $I_{\mathrm{c}}=60$
$\upmu$A, however, manifolds of set conditions and electrical characteristics
of the cycles can be observed (see Table 1), despite it being the set
$I_{\mathrm{c}}$ with the highest yield of conductance quantization in the
reset procedure. Figs. 5(c) and 5(d) both have $I_{\mathrm{c}}=60$ $\upmu$A,
and the key difference between them is the stability of $V_{\mathrm{m}}$ in
the set procedure. One may forecast conductance quantization in the later
reset procedure only when a “wiggling”, unsteady $V_{\mathrm{m}}$ is detected
in the set procedure.

Some previous works have tried to determine the size of the conductive
filaments in Ti/HfO2/TiN memristors either through TEM material analyses 38 or
through theoretical simulations 21, 39. From the electrical characteristics
and the TEM images or simulation results provided in these works, a filament
in the 10-nm thick HfO2 layer of our device with a resistance of around 3
k$\Omega$ to 6 k$\Omega$ may be roughly estimated to be only $\lesssim 3$ nm
in diameter. This explains why it has been extremely difficult to observe a
filament in cross-sectional TEM images of our devices. Electrical stresses on
a narrow filament structure that lead to conductance evolution in the order of
$G_{0}$, on the other hand, have also been studied in several simulation works
14, 40, 20, 41. However, a precise prediction of the quantized conductance
value as a function of the atomic evolution of the filament structure is not
available yet, nor is there a concise conclusion on the numerical values of
the stress voltage to optimize the chance of observing conductance
quantization. More experimental and theoretical research are necessary to
unveil the detailed mechanism of the conductance quantization, and our work
presents a step forward toward understanding the filament evolution.

Although multiple growths of filaments or branches composed of oxygen
vacancies are possible in the device, the conduction is believed to be
contributed by a single dominant filament because there is only one LRS in
each $I$–$V$ cycle (i.e., multiple filaments would have been resulted from
multiple set procedures in sequence in an $I$–$V$ cycle, exhibiting state
switching between multiple LRSs). Once a filament is established, the current
flows mostly through the connected filament, and further filament growth will
be suppressed owing to reduced electric field 42. It has been found from the
TEM images of SiO2-based planar devices that, in the LRS, there exists only
one completed filament accompanied with a few incomplete ones 43.

A good set without SCLC (Fig. 5(c)) may be regarded as an intermediate state
between the state with SCLC (Figs. 5(a) and 5(b)) and the state with
conductance quantization (Fig. 5(d)) in the sense of the resultant filament
strength. For a filament robust enough to stand a higher negative voltage, the
device can enter the SCLC-dominating regime until an abrupt drop of the
current occurs upon the reset procedure when the filament is ruptured under a
much higher voltage. A filament that exhibits a progressive reset, on the
other hand, features a relatively unstable figure showing unsteady
$V_{\mathrm{m}}$ in the set procedure, possibly with the oxygen vacancies at
the tip of the filament moving around among multiple metastable states to
establish or dismiss an ohmic contact. This instability is revisited in the
reset procedure starting around $V=-0.4$ V in a more distinct fashion, i.e.,
the conductance quantization, which again reveals nanoscale atomic-level
movements. As $I_{\mathrm{c}}$ is lowered to 40 $\upmu$A, with examples shown
in Figs. 5(e) and 5(f), 43% of the $I$–$V$ curves exhibit Schottky emissions
or set failures. From the facts above, $I_{\mathrm{c}}=60$ $\upmu$A is the
critical compliance current in our experiment that is just large enough to
build an ohmic filament, and yet simultaneously small enough to allow atomic-
level behaviors to be unveiled in the memristor operation. Our experiment is
the first of its kind to classify the different signatures of
$I$–$V_{\mathrm{m}}$ (the measured voltage) for a better understanding of the
atomic-level dynamics in a memristor.

Electrochemical metallization memristors may behave differently from VCMs in
the reset procedure in accordance with the filament resistance. For example,
Celano et al. 44 have found from Cu/Al2O3/TiN memristors that Cu filaments
with larger diameters and lower LRS resistance tend to exhibit progressive
resets, whereas Cu filaments with smaller diameters and higher LRS resistance
are inclined to undergo abrupt ruptures, which is opposite to our findings on
the HfOx-based devices. The different behaviors can be interpreted with the
filament growth and dissolution dynamic scenario being governed by the ion or
vacancy mobility and diffusivity, and the redox rate 45, 46. The activation
energies for diffusion of oxygen vacancies in HfO2 ($\sim$0.7 eV for a mobile
charged (2+) oxygen vacancy 22 and $\sim$3 eV for a neutral one 47) are much
higher than those of Cu ions in amorphous Al2O3 ($\sim$0.3 eV for mobile ones
48 and $\sim$0.9 eV for less mobile ones 49). As a very narrow oxygen-
deficient filament undergoes a reset procedure with extremely limited current
and thus limited Joule power, the oxygen vacancies, which have significantly
low diffusivity, may leave the filament one by one slowly and discretely,
allowing us to observe the step-wise conductance. The tendency to exhibit
progressive resets at lower $I_{\mathrm{c}}$ and consequently with higher LRS
resistance is typical of memristors based on the VCM mechanism 50, 51.

It is not clear yet why set-and-reset cycles with the same $I_{\mathrm{c}}$
(60 $\upmu$A) can lead to different set conditions for the same memristor.
Possible overshooting of the current is considered at the compliance point as
the memristor is quickly switched from HRS to LRS. The $I$–$V$ curves of our
device exhibit a fairly linear relation in LRS, as shown in Figs. 5(c) and
5(d). This ohmic behavior indicates that possible parasitic capacitance, and
hence current overshooting, are quite limited in our device 52. However,
because the conductance characteristics of the memristor are markedly affected
by the dynamic behaviors in the nanoscale, possible small overshooting of the
current even to a minimal extent may matter. Since it is difficult to directly
detect the probable variations of this minimal overshooting, monitoring
$V_{\mathrm{m}}$ becomes the only practical and effective method to determine
the set condition right away. The cause of the different resultant set
conditions under the same $I_{\mathrm{c}}$ may also involve the instant
internal state of the memristor in the atomic scale being affected dynamically
by the second-order state variables (the temperature decay for example)
present in the structure 53. With these nanoscale uncertainties in the system,
it seems that the most easily observable signal that reveals the multiple
metastable states of a point contact in the filament, and thus the potential
to yield quantized conductance states in the reset procedure, is the wiggling
$V_{\mathrm{m}}$ at the set procedure. Monitoring $V_{\mathrm{m}}$ therefore
becomes the critical method to track the qualities of the device fabrication
and measurements.

###  3.2 Measurements with the Pulse-Mode Reset Procedure

Distinguishing fair sets from the others allows us to achieve a high success
rate of control of the conductance quantization of the memristor using a
pulse-mode reset procedure. A typical example is demonstrated in Fig. 6(a).
The reset process is preceded by a set procedure using $I_{\mathrm{c}}=60$
$\upmu$A that exhibits a fair-set condition (i.e., with wiggling
$V_{\mathrm{m}}$) as depicted in Fig. 5(d). Voltage pulses with fixed width of
0.1 second and fixed value of $-0.35$ V are used to control the atom-by-atom
evolution. The pulse width and amplitude are chosen from amongst multiple
tests to achieve the optimal result, that is, to stimulate switching to the
next conductance state with a minimal average number of pulses. The pulse
value $-0.35$ V, which is very close to the onset voltage of the reset process
observed in dc voltage sweeps with a fair set (Fig. 5(d)), is speculated to be
a favorable value for our device to activate recombination between oxygen ions
and vacancies through oxygen migration by providing a proper electric field
and a local temperature enhancement due to Joule heating 21, 54. After each
pulse, the current is read at $-0.01$ V for 5 seconds, from which the
conductance of the point contact is computed and then presented in units of
$G_{0}$. It can be seen that the conductance decreases stepwise from $9G_{0}$
to $0.5G_{0}$ in steps of $0.5G_{0}$ with great precision, with an average
standard deviation of only $\sim$$0.014G_{0}$ for the quantized plateaus. The
width of each conductance plateau falls within 20 seconds, which corresponds
to 1 to 4 voltage stimuli before stepping down to the next plateau.

Figure 6: (a)(b) Two representative examples of controllable quantized
conductance states at integer multiples of $0.5G_{0}$ using a pulse-mode reset
procedure after a fair set with $I_{\mathrm{c}}=60$ $\upmu$A. (a) is taken
right after changing $I_{\mathrm{c}}$ from 70 $\upmu$A to 60 $\upmu$A. (b) is
taken after three successive cycles with fair sets. (c) Average total number
of exceptional steps, which is the sum of spontaneous steps with $|\Delta
G|\leq G_{0}$ and stimulated steps with $|\Delta G|=G_{0}$, as a function of
the order of successive fair-set cycles. Also displayed are the fastest and
the slowest training sequences from the data.

Fig. 6(b) presents another example of the pulse-mode measurement of the same
memristor after a fair set using the same $I_{\mathrm{c}}$ ($60$ $\upmu$A).
This set of data is from a cycle taken after three successive cycles with fair
sets using $I_{\mathrm{c}}=60$ $\upmu$A. (The quantized conductance data of
all the successive cycles of this sequence are presented in Supporting
Information Section LABEL:sec:training.) In this 4th cycle, not every step in
the reset procedure is $0.5G_{0}$ in height and stimulated by a voltage pulse.
A few exceptions are found between $6G_{0}$ and $2G_{0}$, where some of the
conductance drops have magnitudes of $1G_{0}$, and some occur spontaneously
without a pulsed voltage stimulus. This leads to a faster switching from a
high-conductance state to the lowest one. The occurrence probability of these
exceptions roughly increases with the number of repetitions of the
measurements under the same $I_{\mathrm{c}}$ ($60$ $\upmu$A), implying a
possible training effect. This training effect can be removed by applying a
higher $I_{\mathrm{c}}$ ($70$ $\upmu$A for example) to the set procedure
before returning to the original $I_{\mathrm{c}}$ to regain well-controlled
state switchings in steps of $0.5G_{0}$ without exceptions. The training
effect may be seen from the statistics shown in Fig. 6(c), where the average
number of exceptions (i.e., the number of spontaneous steps with $|\Delta
G|=0.5G_{0}$ or $G_{0}$, plus the number of stimulated steps with $|\Delta
G|=G_{0}$) is plotted against the order of successive fair-set cycles,
collected from 10 sequences in which no steps larger than $G_{0}$ are
involved. (In other words, a very minor group of sequences interrupted by a
reset procedure with steps larger than $G_{0}$ are not included in Fig. 6(c)
for simplicity.) For example, the total number of this kind of exceptional
steps is 4 for Fig. 6(b). Also displayed in Fig. 6(c) are the fastest (dashed
line) and the slowest (solid line) training sequences from the data. It is
rare for the fair set to appear continuously for more than 5 cycles, except
for a “slow learner” like the solid trace in Fig. 6(c), which has 7 successive
fair sets. Some may hope for an ideal memristor with controllable quantized
conduction that always generates a conductance decrease of $0.5$$G_{0}$ in
response to each voltage stimulus, but no realization of such a memristor has
been reported up to date. It is possible that the conductance is governed not
only by the external stimuli but also by its instant internal state, known as
the property of a second-order memristor. In fact, memristors with this kind
of instability on intermediate conductance states are being proposed as
candidate neuromorphic computing devices that can naturally emulate the
temporal behaviors, including sequence learning, of biological synapses 55,
56. More research is necessary to accommodate or even take advantage of the
second-order behaviors of the memristors for constructing practical
neuromorphic computing architectures.

Set conditions other than fair sets generally do not yield controllable
quantized conductance states in the pulse-mode reset procedure. For example,
most of the time a reset process preceded by a good set requires a larger
stimulating voltage, but only to lead to an abrupt conductance drop that
brings the device directly to HRS. In very few cases, a reset process preceded
by a good set that enters LRS at a relatively small $V_{\mathrm{m}}$ (similar
to that in Fig. 5(c)) can exhibit a few quantized conductance plateaus, but is
far away from accessing a complete set of integer multiples of $0.5$$G_{0}$.
Therefore, to efficiently repeat the operation of the quantum-level
manipulation, a pulse-mode reset procedure is executed only when a fair set is
detected. The ability of the nanoscale atomic structure of a filament to
switch among multiple metastable states upon the set process, as implied by
the wiggling $V_{\mathrm{m}}$ during a fair set (Fig. 5(d)), may be a
necessary feature for a memristor to permit excellent realization and
modulability of the quantized conductance states. Future device fabrications
and characterizations are encouraged to incorporate $V_{\mathrm{m}}$
measurements to analyze the set condition. Memristors that guarantee fair sets
under certain $I_{\mathrm{c}}$’s should be favorable.

For comparison, Xue et al. 57 have found from a Pt/HfOx/ITO memristor that
switching of the quantized conductance states needs to be stimulated by
extremely long (20-second) pulses, and becomes even more insensitive to
voltage stimuli at lower conductance, which they attributed to the lower
current and thus a lower power available for modulating the filament. In
contrast to their findings, our devices are more efficient in that they are
sensitive to short voltage stimuli throughout the whole reset procedure. This
indicates the high modulability of a very narrow filament even with a very low
current, which points to the criticality of the current density and the local
temperature enhancement based on heat transfer around the constriction (i.e.,
the narrowest point) of the filament during the reset procedure 21, 54. On the
other hand, there are other previous studies on conductance quantization of
memristors that also demonstrate state switching upon short pulsed voltage
stimuli, but with much lower precision of the quantized conductance values 14.
As analytical differentiation of the $I$–$V_{\mathrm{m}}$ characteristics
(Fig. 5) is employed in our experiment during device selection and
measurements, the precision and signal-to-noise ratios of the quantized
conductance are significantly improved in our experiment compared to previous
reports, and thereby brings the study of memristors closer to practical
application in neuromorphic computing.

##  4 Conclusions

In summary, we report on controllable quantized conductance states of
TiN/Ti/HfOx/TiN memristors in a pulse-mode reset procedure with significantly
improved precision. The high controllability and precision are realized
through analytical diagnoses of the set conditions of the fabricated devices.
The $I$–$V_{\mathrm{m}}$ characteristics of the set procedure can be
classified into good, fair, and poor conditions, and only those with fair sets
(i.e., with a “wiggling”, unstable measured voltage $V_{\mathrm{m}}$ at the
compliance current) can permit quantized conductance states in the reset
procedure. Controlled conductance decrease from $9G_{0}$ to $0.5G_{0}$ in
steps of $0.5G_{0}$ is successfully observed in pulse-mode reset procedures
that are preceded by a fair set with an optimal compliance current ($60$
$\upmu$A). A training effect that leads to a faster state switching is found
in the operation, which is regarded as a candidate mechanism for temporal
sequence learning. Our experiment is the first of its kind to point out the
importance of monitoring the measured bias voltage to track the qualities of
the device fabrication and measurements for the research of conductance
quantization of memristors. Our study unveils a full spectrum of the dynamic
behaviors under different set conditions to provide an overview of the
mechanisms of the conductive filament, from a strong ohmic structure with
space-charge-limited current (SCLC), to that without SCLC, then to a
relatively unstable configuration that supports quantized conduction as well
as ohmic conduction, and then to a nano-gapped channel with Schottky emission.
This allows a better understanding of the dynamics of the nanoscale atomic-
level structures in the memristors, which should promote the progress of
future design and fabrication of the memristors for neuromorphic computing and
quantum information processing algorithm.

## Associated Content

Supporting Information available: more examples for each set category;
information pertaining to the forming procedure, SCLC statistics, and
operation using a pulse-mode set procedure; measurements of devices with
different cell areas; temperature-dependent measurements down to 4.2 K; data
of a complete training sequence.

## Acknowledgements

We acknowledge the groups of Prof. Shu-Fen Hu and Prof. Yann-Wen Lan for
assisting us with the measurements. We thank Dr. Chih-Yang Lin for helpful
discussion. This study is sponsored by the Ministry of Science and Technology
of Taiwan under Grant No. MOST 109-2112-M-003-009 and MOST 110-2112-M-003-019.

## References

  * Chua 1971 Chua, L. Memristor-The missing circuit element. _IEEE Trans. Circuit Theory_ 1971, _18_ , 507–519
  * Yang et al. 2013 Yang, J. J.; Strukov, D. B.; Stewart, D. R. Memristive devices for computing. _Nat. Nanotechnol._ 2013, _8_ , 13–24
  * Milo et al. 2020 Milo, V.; Malavena, G.; Monzio Compagnoni, C.; Ielmini, D. Memristive and CMOS Devices for Neuromorphic Computing. _Materials (Basel)._ 2020, _13_ , 166
  * Li et al. 2015 Li, B.; Gu, P.; Shan, Y.; Wang, Y.; Chen, Y.; Yang, H. RRAM-Based Analog Approximate Computing. _IEEE Trans. Comput. Des. Integr. Circuits Syst._ 2015, _34_ , 1905–1917
  * Merced-Grafals et al. 2016 Merced-Grafals, E. J.; Dávila, N.; Ge, N.; Williams, R. S.; Strachan, J. P. Repeatable, accurate, and high speed multi-level programming of memristor 1T1R arrays for power efficient analog computing applications. _Nanotechnology_ 2016, _27_ , 365202
  * Sun et al. 2018 Sun, Z.; Ambrosi, E.; Bricalli, A.; Ielmini, D. Logic Computing with Stateful Neural Networks of Resistive Switches. _Adv. Mater._ 2018, _30_ , 1802554
  * Pickett et al. 2009 Pickett, M. D.; Strukov, D. B.; Borghetti, J. L.; Yang, J. J.; Snider, G. S.; Stewart, D. R.; Williams, R. S. Switching dynamics in titanium dioxide memristive devices. _J. Appl. Phys._ 2009, _106_ , 074508
  * Menzel et al. 2011 Menzel, S.; Waters, M.; Marchewka, A.; Böttger, U.; Dittmann, R.; Waser, R. Origin of the Ultra-nonlinear Switching Kinetics in Oxide-Based Resistive Switches. _Adv. Funct. Mater._ 2011, _21_ , 4487–4492
  * Strachan et al. 2013 Strachan, J. P.; Torrezan, A. C.; Miao, F.; Pickett, M. D.; Yang, J. J.; Yi, W.; Medeiros-Ribeiro, G.; Williams, R. S. State Dynamics and Modeling of Tantalum Oxide Memristors. _IEEE Trans. Electron Devices_ 2013, _60_ , 2194–2202
  * Lim et al. 2018 Lim, S.; Sung, C.; Kim, H.; Kim, T.; Song, J.; Kim, J.-J.; Hwang, H. Improved Synapse Device With MLC and Conductance Linearity Using Quantized Conduction for Neuromorphic Systems. _IEEE Electron Device Lett._ 2018, _39_ , 312–315
  * Chen et al. 2020 Chen, Q.; Han, T.; Tang, M.; Zhang, Z.; Zheng, X.; Liu, G. Improving the Recognition Accuracy of Memristive Neural Networks via Homogenized Analog Type Conductance Quantization. _Micromachines_ 2020, _11_ , 427
  * Marković and Grollier 2020 Marković, D.; Grollier, J. Quantum neuromorphic computing. _Appl. Phys. Lett._ 2020, _117_ , 150501
  * Pfeiffer et al. 2016 Pfeiffer, P.; Egusquiza, I. L.; Di Ventra, M.; Sanz, M.; Solano, E. Quantum memristors. _Sci. Rep._ 2016, _6_ , 29507
  * Li et al. 2015 Li, Y.; Long, S.; Liu, Y.; Hu, C.; Teng, J.; Liu, Q.; Lv, H.; Suñé, J.; Liu, M. Conductance Quantization in Resistive Random Access Memory. _Nanoscale Res. Lett._ 2015, _10_ , 420
  * Xue et al. 2019 Xue, W.; Gao, S.; Shang, J.; Yi, X.; Liu, G.; Li, R. Recent Advances of Quantum Conductance in Memristors. _Adv. Electron. Mater._ 2019, _5_ , 1800854
  * Ohno et al. 2011 Ohno, T.; Hasegawa, T.; Tsuruoka, T.; Terabe, K.; Gimzewski, J. K.; Aono, M. Short-term plasticity and long-term potentiation mimicked in single inorganic synapses. _Nat. Mater._ 2011, _10_ , 591–595
  * Yi et al. 2016 Yi, W.; Savel’ev, S. E.; Medeiros-Ribeiro, G.; Miao, F.; Zhang, M.-X.; Yang, J. J.; Bratkovsky, A. M.; Williams, R. S. Quantized conductance coincides with state instability and excess noise in tantalum oxide memristors. _Nat. Commun._ 2016, _7_ , 11142
  * Xie et al. 2020 Xie, Z.; Gao, S.; Ye, X.; Yang, H.; Gong, G.; Lu, Y.; Ye, J.; Liu, G.; Li, R.-W. Magnetism modulation and conductance quantization in a gadolinium oxide memristor. _Phys. Chem. Chem. Phys._ 2020, _22_ , 26322–26329
  * Lv et al. 2015 Lv, H.; Xu, X.; Sun, P.; Liu, H.; Luo, Q.; Liu, Q.; Banerjee, W.; Sun, H.; Long, S.; Li, L.; Liu, M. Atomic View of Filament Growth in Electrochemical Memristive Elements. _Sci. Rep._ 2015, _5_ , 13311
  * Long et al. 2013 Long, S.; Lian, X.; Cagli, C.; Cartoixà, X.; Rurali, R.; Miranda, E.; Jiménez, D.; Perniola, L.; Liu, M.; Suñé, J. Quantum-size effects in hafnium-oxide resistive switching. _Appl. Phys. Lett._ 2013, _102_ , 183505
  * Dirkmann et al. 2018 Dirkmann, S.; Kaiser, J.; Wenger, C.; Mussenbrock, T. Filament Growth and Resistive Switching in Hafnium Oxide Memristive Devices. _ACS Appl. Mater. Interfaces_ 2018, _10_ , 14857–14868
  * Capron et al. 2007 Capron, N.; Broqvist, P.; Pasquarello, A. Migration of oxygen vacancy in HfO2 and across the HfO2/SiO2 interface: A first-principles investigation. _Appl. Phys. Lett._ 2007, _91_ , 192905
  * Syu et al. 2013 Syu, Y.-E.; Chang, T.-C.; Lou, J.-H.; Tsai, T.-M.; Chang, K.-C.; Tsai, M.-J.; Wang, Y.-L.; Liu, M.; Sze, S. M. Atomic-level quantized reaction of HfOx memristor. _Appl. Phys. Lett._ 2013, _102_ , 172903
  * Tappertzhofen et al. 2015 Tappertzhofen, S.; Linn, E.; Menzel, S.; Kenyon, A. J.; Waser, R.; Valov, I. Modeling of Quantized Conductance Effects in Electrochemical Metallization Cells. _IEEE Trans. Nanotechnol._ 2015, _14_ , 505–512
  * Mehonic et al. 2013 Mehonic, A.; Vrajitoarea, A.; Cueff, S.; Hudziak, S.; Howe, H.; Labbé, C.; Rizk, R.; Pepper, M.; Kenyon, A. J. Quantum Conductance in Silicon Oxide Resistive Memory Devices. _Sci. Rep._ 2013, _3_ , 2708
  * Krishnan et al. 2017 Krishnan, K.; Muruganathan, M.; Tsuruoka, T.; Mizuta, H.; Aono, M. Highly Reproducible and Regulated Conductance Quantization in a Polymer-Based Atomic Switch. _Adv. Funct. Mater._ 2017, _27_ , 1605104
  * Malliaras and Scott 1999 Malliaras, G. G.; Scott, J. C. Numerical simulations of the electrical characteristics and the efficiencies of single-layer organic light emitting diodes. _J. Appl. Phys._ 1999, _85_ , 7426–7432
  * Kao and Hwang 1981 Kao, K. C.; Hwang, W. _Electrical Transport in Solids_ ; Pergamon: New York, 1981
  * Lampert and Mark 1970 Lampert, M. A.; Mark, P. _Current Injection in Solids_ ; Academic: New York, 1970
  * Mark and Helfrich 1962 Mark, P.; Helfrich, W. Space‐Charge‐Limited Currents in Organic Crystals. _J. Appl. Phys._ 1962, _33_ , 205–215
  * Mott and Gurney 1940 Mott, N. F.; Gurney, R. W. _Electronic Processes in Ionic Crystals_ ; Oxford University Press, 1940
  * Negara et al. 2009 Negara, M. A.; Cherkaoui, K.; Hurley, P. K.; Young, C. D.; Majhi, P.; Tsai, W.; Bauza, D.; Ghibaudo, G. Analysis of electron mobility in HfO2/TiN gate metal-oxide-semiconductor field effect transistors: The influence of HfO2 thickness, temperature, and oxide charge. _J. Appl. Phys._ 2009, _105_ , 024510
  * Özben 2010 Özben, E. D. _Carrier mobility in advanced channel materials using alternative gate dielectrics_ ; Forschungszentrum Jülich, 2010
  * Bässler 1993 Bässler, H. Charge Transport in Disordered Organic Photoconductors a Monte Carlo Simulation Study. _Phys. status solidi_ 1993, _175_ , 15–56
  * Bhattarai et al. 2018 Bhattarai, G.; Caruso, A. N.; Paquette, M. M. Steady-state space-charge-limited current analysis of mobility with negative electric field dependence. _J. Appl. Phys._ 2018, _124_ , 045701
  * Zahoor et al. 2020 Zahoor, F.; Azni Zulkifli, T. Z.; Khanday, F. A. Resistive Random Access Memory (RRAM): an Overview of Materials, Switching Mechanism, Performance, Multilevel Cell (mlc) Storage, Modeling, and Applications. _Nanoscale Res. Lett._ 2020, _15_ , 90
  * Yu 2016 Yu, S. Resistive Random Access Memory (RRAM). _Synth. Lect. Emerg. Eng. Technol._ 2016, _2_ , 1–79
  * Privitera et al. 2013 Privitera, S.; Bersuker, G.; Butcher, B.; Kalantarian, A.; Lombardo, S.; Bongiorno, C.; Geer, R.; Gilmer, D.; Kirsch, P. Microscopy study of the conductive filament in HfO2 resistive switching memory devices. _Microelectron. Eng._ 2013, _109_ , 75–78
  * Nardi et al. 2012 Nardi, F.; Larentis, S.; Balatti, S.; Gilmer, D. C.; Ielmini, D. Resistive Switching by Voltage-Driven Ion Migration in Bipolar RRAM—Part I: Experimental Study. _IEEE Trans. Electron Devices_ 2012, _59_ , 2461–2467
  * Li et al. 2017 Li, Y.; Zhang, M.; Long, S.; Teng, J.; Liu, Q.; Lv, H.; Miranda, E.; Suñé, J.; Liu, M. Investigation on the Conductive Filament Growth Dynamics in Resistive Switching Memory via a Universal Monte Carlo Simulator. _Sci. Rep._ 2017, _7_ , 11204
  * Zhong et al. 2016 Zhong, X.; Rungger, I.; Zapol, P.; Heinonen, O. Oxygen-modulated quantum conductance for ultrathin HfO2-based memristive switching devices. _Phys. Rev. B_ 2016, _94_ , 165160
  * Kwon et al. 2010 Kwon, D.-H.; Kim, K. M.; Jang, J. H.; Jeon, J. M.; Lee, M. H.; Kim, G. H.; Li, X.-S.; Park, G.-S.; Lee, B.; Han, S.; Kim, M.; Hwang, C. S. Atomic structure of conducting nanofilaments in TiO2 resistive switching memory. _Nat. Nanotechnol._ 2010, _5_ , 148–153
  * Yang et al. 2012 Yang, Y.; Gao, P.; Gaba, S.; Chang, T.; Pan, X.; Lu, W. Observation of conducting filament growth in nanoscale resistive memories. _Nat. Commun._ 2012, _3_ , 732
  * Celano et al. 2015 Celano, U.; Goux, L.; Belmonte, A.; Opsomer, K.; Degraeve, R.; Detavernier, C.; Jurczak, M.; Vandervorst, W. Understanding the Dual Nature of the Filament Dissolution in Conductive Bridging Devices. _J. Phys. Chem. Lett._ 2015, _6_ , 1919–1924
  * Yang and Lu 2013 Yang, Y.; Lu, W. Nanoscale resistive switching devices: mechanisms and modeling. _Nanoscale_ 2013, _5_ , 10076
  * Ielmini 2016 Ielmini, D. Resistive switching memories based on metal oxides: mechanisms, reliability and scaling. _Semicond. Sci. Technol._ 2016, _31_ , 063002
  * Duncan et al. 2016 Duncan, D.; Magyari-Kope, B.; Nishi, Y. Filament-Induced Anisotropic Oxygen Vacancy Diffusion and Charge Trapping Effects in Hafnium Oxide RRAM. _IEEE Electron Device Lett._ 2016, _37_ , 400–403
  * Pandey 2018 Pandey, S. C. Atomistic mechanisms of ReRAM cell operation and reliability. _Mater. Res. Express_ 2018, _5_ , 014005
  * Sankaran et al. 2012 Sankaran, K.; Goux, L.; Clima, S.; Mees, M.; Kittl, J. A.; Jurczak, M.; Altimime, L.; Rignanese, G.-M.; Pourtois, G. Modeling of Copper Diffusion in Amorphous Aluminum Oxide in CBRAM Memory Stack. _ECS Trans._ 2012, _45_ , 317–330
  * Wouters et al. 2019 Wouters, D. J.; Menzel, S.; Rupp, J. A. J.; Hennen, T.; Waser, R. On the universality of the I – V switching characteristics in non-volatile and volatile resistive switching oxides. _Faraday Discuss._ 2019, _213_ , 183–196
  * 51 Volos, C.; Pham, V.-T. Mem-elements for Neuromorphic Circuits with Artificial Intelligence Applications; Academic Press, 2021; p. 392.
  * Ambrogio et al. 2016 Ambrogio, S.; Milo, V.; Wang, Z.; Balatti, S.; Ielmini, D. Analytical Modeling of Current Overshoot in Oxide-Based Resistive Switching Memory (RRAM). _IEEE Electron Device Lett._ 2016, _37_ , 1268–1271
  * Kim et al. 2015 Kim, S.; Du, C.; Sheridan, P.; Ma, W.; Choi, S.; Lu, W. D. Experimental Demonstration of a Second-Order Memristor and Its Ability to Biorealistically Implement Synaptic Plasticity. _Nano Lett._ 2015, _15_ , 2203–2211
  * Roldán et al. 2021 Roldán, J. B.; González-Cordero, G.; Picos, R.; Miranda, E.; Palumbo, F.; Jiménez-Molinos, F.; Moreno, E.; Maldonado, D.; Baldomá, S. B.; Moner Al Chawa, M.; de Benito, C.; Stavrinides, S. G.; Suñé, J.; Chua, L. O. On the Thermal Models for Resistive Random Access Memory Circuit Simulation. _Nanomaterials_ 2021, _11_ , 1261
  * Mikheev et al. 2019 Mikheev, V.; Chouprik, A.; Lebedinskii, Y.; Zarubin, S.; Matveyev, Y.; Kondratyuk, E.; Kozodaev, M. G.; Markeev, A. M.; Zenkevich, A.; Negrov, D. Ferroelectric Second-Order Memristor. _ACS Appl. Mater. Interfaces_ 2019, _11_ , 32108–32114
  * Marrone et al. 2019 Marrone, F.; Zoppo, G.; Corinto, F.; Gilli, M. Second Order Memristor Models for Neuromorphic Computing. 2019 IEEE 62nd Int. Midwest Symp. Circuits Syst. 2019; pp 315–318
  * Xue et al. 2020 Xue, W.; Li, Y.; Liu, G.; Wang, Z.; Xiao, W.; Jiang, K.; Zhong, Z.; Gao, S.; Ding, J.; Miao, X.; Xu, X.; Li, R. Controllable and Stable Quantized Conductance States in a Pt/HfOx/ITO Memristor. _Adv. Electron. Mater._ 2020, _6_ , 1901055