# Quantum theory of the magnetochiral anisotropy coefficient in ZrTe5

Yi-Xiang Wang<EMAIL_ADDRESS>School of Science, Jiangnan
University, Wuxi 214122, China  School of Physics and Electronics, Hunan
University, Changsha 410082, China  Fuxiang Li<EMAIL_ADDRESS>School of
Physics and Electronics, Hunan University, Changsha 410082, China

###### Abstract

Recent experiments performed the nonreciprocal magneotransport in ZrTe5 and
obtained a giant magnetochiral anisotropy (MCA) coefficient $\gamma^{\prime}$.
The existing theoretical analysis was based on the semiclassical Boltzmann
equation. In this paper, we develop a full quantum theory to calculate
$\gamma^{\prime}$. We reveal that the $xz$-mirror symmetry breaking term also
breaks the parity symmetry of the system, and leads to the mixed selection
rules and a nonvanishing second-order conductivity $\sigma_{xxx}$. The
calculations show that $\gamma^{\prime}$ decreases with the magnetic field,
survives only to weak impurity scatterings, and grows almost linearly with the
strength of the $xz$-mirror symmetry breaking. Our paper can provide a deeper
insight into the intrinsic nonreciprocal magnetotransport phenomena in the
topological semimetal material.

##  I Introduction

Symmetries in the microscopic dynamics are vital in determining the
macroscopic physical reponses to the external excitations. One example is that
the quantum anomalous Hall effect can emerge in a system with broken time-
reversal symmetry (TRS) Haldane ; C.Z.Chang2013 ; C.Z.Chang2023 , while the
quantum spin Hall effect is protected by the TRS M.Z.Hasan ; X.L.Qi . Another
example is that magnetochiral anisotropy (MCA) will occur in a system with
broken inversion symmetry Y.Tokura . The MCA means that under a magnetic field
$B$, the resistance of a three-dimensional (3D) crystal sample is
nonreciprocal when an electric current flows in the opposite directions. In
this case, the nonreciprocal resistance $R$ includes both the linear and
quadratic terms, and is written as

$\displaystyle R=\frac{V_{x}}{I_{x}}=R_{0}(1+\gamma BI_{x}),$ (1)

where $I_{x}$ is the electric current that we assume to flow along the $x$
direction, $V_{x}$ is the voltage, $R_{0}$ is the reciprocal resistance, and
the coefficient $\gamma$ characterizes the magnitude of MCA. Eq. (1) is valid
for both the inner-product-type nonreciprocal behavior Rikken2001  and the
vector-product-type nonreciprocal behavior Rikken2005 , for which the electric
field and magnetic field are parallel $\bm{B}\parallel\bm{I}$ and
perpendicular $\bm{B}\perp\bm{I}$, respectively.

On the other hand, the current density $j_{x}$ can be written as a function of
the electric field $E_{x}$,

$\displaystyle
j_{x}=j_{x}^{(1)}+j_{x}^{(2)}=\sigma_{xx}E_{x}+\sigma_{xxx}E_{x}^{2},$ (2)

in which $\sigma_{xx}$ and $\sigma_{xxx}$ are the first-order and second-order
longitudinal conductivities, respectively. If the current is expressed as
$I_{x}=L_{y}L_{z}j_{x}$ and the voltage as $V_{x}=L_{x}E_{x}$, with
$L_{x,y,z}$ denoting the size of the sample, we obtain
$\gamma=\frac{\sigma_{xxx}}{L_{y}L_{z}B\sigma_{xx}^{2}}$. Since $\gamma$
depends on the size and shape of the 3D system, we instead use the cross
section-independent coefficient

$\displaystyle\gamma^{\prime}=L_{y}L_{z}\gamma=\frac{\sigma_{xxx}}{B\sigma_{xx}^{2}},$
(3)

to characterize MCA T.Ideue .

The common bulk nonreciprocal response in a 3D crystal originates from the
relativistic effect and are quite small Rikken2005 . The experiments in
trigonal tellurium showed that MCA had an inner-product type with the
coefficient as large as $\gamma^{\prime}\sim 10^{-8}$ m2T-1A-1 Rikken2019 ,
where the inversion symmetry was broken by its helical crystal structure. The
recent magnetotransport experiments in 3D ZrTe5 at low temperatures reported
that MCA had a vector-product type, in which $\gamma^{\prime}$ could reach the
order of $10^{-7}$ m2T-1A-1 and was gigantic Y.Wang ; N.Wang . The microscopic
mechanism was identified as the Fermi surface deformation due to the magnetic
field-induced Zeeman effect T.Ideue ; Y.Wang ; P.He . The supporting
theoretical analysis adopted the semiclassical Boltzmann equation (SBE) and
the conductivities were calculated by the Drude formula H.Watanabe ; C.Wang ,
from which $\gamma^{\prime}\sim 10^{-11}$ m2T-1A-1 and was four orders of
magnitude smaller than the experimental value Y.Wang . Although the SBE based
on the wave packet dynamics was widely used in the nonlinear conductivity
studies I.Sodemann ; T.Morimoto ; Z.Z.Du ; M.X.Deng ; Y.Gao ; D.Sinha , a full
quantum theory should be developed and applied in order to understand the
underlying physics in the nonreciprocal magnetotransport.

In this paper, with the help of nonequilibrium density matrix and
thermodynamical averaging over current operators, we derive the conductivity
formula of $\sigma_{xx}$ and $\sigma_{xxx}$. Such a quantized theory has the
advantages that it can well deal with the band characters of the system, and
some other factors, such as the impurity scatterings. Within the topological
semimetal model describing 3D ZrTe5 H.Weng ; R.Y.Chen ; Z.G.Chen ; F.Tang ;
S.Galeski2021 , where the inversion symmetry is broken intrinsically, we make
a systematic calculation of $\gamma^{\prime}$ and analyze the effects of the
magnetic field, the impurity scatterings as well as the $xz$-mirror symmetry
breaking term.

Our main results are given as follows: (i) the intrinsic $xz$-mirror symmetry
breaking term also breaks the parity symmetry of the system, which results in
the mixed selection rules, $n\rightarrow n,n\pm 1,n\pm 2,\cdots$, with $n$ the
LL index, and plays a decisive role in determining $\sigma_{xxx}$; (ii)
$\gamma^{\prime}$ decreases with the magnetic field $B$, and is nonvanishing
at a weak $B$, which are consistent with the experimental observations Y.Wang
; N.Wang ; (iii) $\gamma^{\prime}$ is susceptible to the impurity scatterings,
and survives only to the weak impurity scatterings; and (iv) more importantly,
$\gamma^{\prime}$ grows linearly with the strength of the $xz$-mirror symmetry
breaking, which is understood by its role played in $\sigma_{xxx}$. It is
worth noting that the last two points are distinct from the SBE results Y.Wang
and the implications will be discussed. Our paper can help better understand
the MCA phenomenon in ZrTe5 experiments.

##  II The first- and second-order conductivities

When an electric field ${\bm{E}}=E{\bm{e}}_{x}$ acts on a system, the
electrostatic force will drive the the electrons deviate from its equilibrium
state. If we consider the electric field as a perturbation to the system, the
nonequilibrium density matrix up to the second order is written as L.Smrcka ;
G.D.Mahan

$\displaystyle\rho=\rho_{0}+\rho_{1}e^{st}+\rho_{2}e^{2st}+\cdots,$ (4)

where $\rho_{0}=1/[e^{\beta(\hat{H}_{0}-\mu)}+1]$ is the equilibrium Fermi-
Dirac distribution function with $\mu$ denoting the chemical potential, and
$\beta=\frac{1}{k_{B}T}$ with $k_{B}$ being the Boltzmann constant and $T$ the
temperature. $\rho_{1}$ and $\rho_{2}$ are the first-order and second-order
density matrix, respectively. The factors $e^{st}$ and $e^{2st}$ are added to
$\rho_{1}$ and $\rho_{2}$ with $s>0$ to ensure that they vanish at
$t\rightarrow-\infty$. After a direct calculation, $\rho_{1}$ and $\rho_{2}$
are obtained as ($\hbar=1$) (see Supplemental Material (SM) SM )

$\displaystyle\rho_{1}=\frac{1}{i}\int_{0}^{\infty}dte^{-st}[\hat{H}^{\prime}(-t),\rho_{0}],$
(5)
$\displaystyle\rho_{2}=\frac{1}{i^{2}}\int_{0}^{\infty}dte^{-2st}\int_{0}^{\infty}dt^{\prime}e^{-st^{\prime}}\big{[}\hat{H}^{\prime}(-t),[\hat{H}^{\prime}(-t^{\prime}),\rho_{0}]\big{]}.$
(6)

where
$\hat{H}^{\prime}(t)=e^{i\hat{H}_{0}t}\hat{H}^{\prime}e^{-i\hat{H}_{0}t}$ and
$\hat{H}^{\prime}=-e{\bm{E}}\cdot{\hat{\bm{r}}}$ describes the interaction
between the system and the electric field, with $\hat{\bm{r}}$ being the
position operator.

The measurable first-order and second-order current densities are calculated
through the thermodynamical averaging over the current density operator
$\hat{\bm{j}}=e\hat{\bm{v}}$, i.e.,
${\bm{j}}^{(1)}=\text{Tr}[\rho_{1}\hat{\bm{j}}]$ and
${\bm{j}}^{(2)}=\text{Tr}[\rho_{2}\hat{\bm{j}}]$. With the help of the Green’s
function, the first-order and second-order conductivities, $\sigma_{xx}$ and
$\sigma_{xxx}$, can be derived (see SM SM ). At zero temperature, their
expressions read

$\displaystyle\sigma_{xx}=\frac{e^{2}\eta^{2}}{\pi
V}\sum_{\bm{k}}\sum_{n}\sum_{n^{\prime}}\frac{|\langle\psi_{n}|\hat{v}_{x}|\psi_{n^{\prime}}\rangle|^{2}}{[(\mu-\varepsilon_{n})^{2}+\eta^{2}][(\mu-\varepsilon_{n^{\prime}})^{2}+\eta^{2}]},$
(7)

and

$\displaystyle\sigma_{xxx}=$
$\displaystyle\frac{2e^{3}}{V}\sum_{\bm{k}}\sum_{\varepsilon_{n}<\mu}\sum_{\varepsilon_{n^{\prime}}>\mu}\sum_{n^{\prime\prime}}$
$\displaystyle\text{Re}\big{[}\frac{\langle\psi_{n}|\hat{v}_{x}|\psi_{n^{\prime}}\rangle\langle\psi_{n^{\prime}}|\hat{v}_{x}|\psi_{n^{\prime\prime}}\rangle\langle\psi_{n^{\prime\prime}}|\hat{v}_{x}|\psi_{n}\rangle}{(\varepsilon_{n}-\varepsilon_{n^{\prime}}-i\eta)^{2}(\varepsilon_{n}-\varepsilon_{n^{\prime\prime}}+i\eta)^{2}}$
$\displaystyle-\frac{\langle\psi_{n}|\hat{v}_{x}|\psi_{n^{\prime\prime}}\rangle\langle\psi_{n^{\prime\prime}}|\hat{v}_{x}|\psi_{n^{\prime}}\rangle\langle\psi_{n^{\prime}}|\hat{v}_{x}|\psi_{n}\rangle}{(\varepsilon_{n}-\varepsilon_{n^{\prime}}+i\eta)^{2}(\varepsilon_{n^{\prime\prime}}-\varepsilon_{n^{\prime}}-i\eta)^{2}}\big{]},$
(8)

respectively. Here $V$ is the volume of the system, $\varepsilon_{n}$ is the
energy and $|\psi_{n}\rangle$ the corresponding eigenstate, $\hat{v}_{x}$ is
the velocity operator, and $\eta$ denotes the linewidth broadening that is
introduced to represent the impurity scatterings phenomenologically. For
simplicity, we take $\eta$ as a constant for all eigenstates in the system.

Note that the expression of $\sigma_{xx}$ is the same as those obtained from
the Kubo-Streda formula in the previous studies H.W.Wang ; Y.X.Wang2020 . In
$\sigma_{xxx}$, $|\psi_{n}\rangle$ and $|\psi_{n^{\prime}}\rangle$ represent
the initial and final state, respectively, while
$|\psi_{n^{\prime\prime}}\rangle$ denotes the middle state. Intuitively, since
$\sigma_{xxx}$ does not include the $\hat{v}_{y}$ term, it has no relation to
the Berry curvature and thus is distinct from the nonlinear Hall conductivity
C.Wang . Eqs. (7) and (8) relate the conductivities to the band structures of
the system directly, which are transparent to interpretation.

The presence of a magnetic field will drive the formation of the Landau levels
(LLs). In Eqs. (7) and (8), the nonvanishing velocity operator matrix element
determines the selection rules that will play the decisive roles in the
transport properties. In $\sigma_{xx}$, there are no restrictions on the
selection rules: they may take the conventional ones $n\rightarrow n,n\pm 1$
Y.Jiang ; Y.X.Wang2023 ; Z.Cai , with $n$ the LL index, or the unconventional
ones $n\rightarrow n\pm 2$ Y.X.Wang2022 . But in $\sigma_{xxx}$, there exist
multiplications of three matrix elements that form a closed loop, requiring
that the selection rules would not simply take $n\rightarrow n\pm 1,n\pm 2$;
they should take $n\rightarrow n$ or other types.

##  III Topological Semimetal Model and LLs

Figure 1: (Color online) The energy bands versus $k_{x}$ for $k_{y}=k_{z}=0$,
with the band index $(s\lambda)$ being labeled. In (a) when $m=0$, the system
lies in the nodal line semimetal phase, while in (b) when $m\neq 0$, the
system is a gapped insulator. We set the parameter $\xi=0$.

We use the topological semimetal model to describe the low-energy excitations
in 3D ZrTe5. In the four-component basis
$\big{(}|+\uparrow\rangle,|-\uparrow\rangle,|+\downarrow\rangle,|-\downarrow\rangle\big{)}^{T}$,
the Hamiltonian $\hat{H}_{0}(\bm{k})$ is written as H.Weng ; R.Y.Chen ;
Z.G.Chen ; F.Tang ; S.Galeski2021

$\displaystyle\hat{H}_{0}(\bm{k})=$ $\displaystyle
v(k_{x}\sigma_{z}\otimes\tau_{x}+k_{y}I\otimes\tau_{y})+v_{z}k_{z}\sigma_{x}\otimes\tau_{x}$
$\displaystyle+mI\otimes\tau_{z},$ (9)

where $\sigma$ and $\tau$ are the Pauli matrices acting on the spin and orbit
degrees of freedom, respectively. $v$ is the Fermi velocity in the $xy$ plane,
$v_{z}$ is the Fermi velocity in the $z$ direction, and $m$ denotes the Dirac
mass. In the 3D system, there exist the $xz$-mirror symmetry, $\hat{\cal
M}_{xz}^{-1}{\hat{H}}_{0}(k_{x},k_{y},k_{z})\hat{\cal
M}_{xz}={\hat{H}}_{0}(k_{x},-k_{y},k_{z})$, with the operator $\hat{\cal
M}_{xz}=i\sigma_{y}\otimes\tau_{z}$, the $xy$-mirror symmetry, $\hat{\cal
M}_{xy}^{-1}{\hat{H}}_{0}(k_{x},k_{y},k_{z})\hat{\cal
M}_{xy}={\hat{H}}_{0}(k_{x},k_{y},-k_{z})$, with $\hat{\cal
M}_{xy}=\sigma_{z}\otimes I$, and the inversion symmetry $\hat{\cal
I}^{-1}{\hat{H}}_{0}(\bm{k})\hat{\cal I}={\hat{H}}_{0}(-\bm{k})$, with
$\hat{\cal I}=I\otimes\tau_{z}$. In the following, we take the model
parameters $(v,v_{z})=(6,0.5)\times 10^{5}$ m/s, and $m=5$ meV.

In a recent magnetotransport experiment Y.Wang , the inversion symmetry of 3D
ZrTe5 was demonstrated to be broken intrinsically, which was attributed to the
staggered displacement of Te atom along the $y$ direction. Such inversion
symmetry breaking terms are represented by

$\displaystyle\hat{H}_{IB}=\Delta
I\otimes\tau_{x}+\xi\sigma_{x}\otimes\tau_{y}.$ (10)

in which the parameters $\Delta$ and $\xi$ denote the strengths of the $xz$\-
and $xy$-mirror symmetry breaking, respectivley. Clearly, the introduction of
$\hat{H}_{IB}$ breaks the inversion symmetry as $\hat{\cal
I}^{-1}{\hat{H}}_{IB}\hat{\cal I}=-{\hat{H}}_{IB}$.

Figure 2: (Color online) The LL spectra in (a)$-$(c), with the index
$ns\lambda$ being labeled, and the velocity operator matrix element
$|\langle\psi_{ns\lambda}|\hat{v}_{x}|\psi_{n^{\prime}s^{\prime}\lambda^{\prime}}\rangle|$
between the initial state and the final state for different $\Delta$ in (d),
with the initial state index chosen as $(3--)$. The parameters are set as
$(\Delta,\xi)=(0,0)$ in (a), $(\Delta,\xi)=(0,6)$ meV in (b),
$(\Delta,\xi)=(10,6)$ meV in (c), the magnetic field $B=1$ T and $\xi=6$ meV
in (d). In (a) and (b), each LL owns a definite even or odd parity that are
represented by the red solid or blue dashed lines, respectively. In (d), the
neighboring bars are shifted vertically for clarity.

The energies of $\hat{H}_{0}+\hat{H}_{IB}$ are obtained as

$\displaystyle\varepsilon_{s\lambda}(\bm{k})=$ $\displaystyle
s\Big{[}\big{(}\sqrt{v^{2}k_{x}^{2}+v_{z}^{2}k_{z}^{\prime
2}}+\lambda\sqrt{\Delta^{2}+\xi^{2}}\big{)}^{2}+v^{2}k_{y}^{\prime 2}$
$\displaystyle+m^{2}\Big{]}^{\frac{1}{2}},$ (11)

where the index $s=\pm 1$ denotes the conduction/valence band, and
$\lambda=\pm 1$ the two branches. The redefined wave vectors are
$vk_{y}^{\prime}=vk_{y}\text{cos}\phi-v_{z}k_{z}\text{sin}\phi$,
$v_{z}k_{z}^{\prime}=vk_{y}\text{sin}\phi+v_{z}k_{z}\text{cos}\phi$, and
tan$\phi=\frac{\xi}{\Delta}$, which is equivalent to rotate the $yz$ plane
around the $x$ axis by angle $-\phi$. In Fig. 1, the energy bands are plotted
for $k_{y}=k_{z}=0$. We see that the two $\lambda=1$ branches have a gap of
$2\delta_{1}=2\sqrt{m^{2}+\Delta^{2}+\xi^{2}}$; while another two $\lambda=-1$
branches cross at zero energy for $m=0$ [Fig. 1(a)] and are gapped for a
nonvanishing $m$ [Fig. 1(b)]. The former occurs when
$k_{x}=\pm(\Delta^{2}+\xi^{2}-v_{z}^{2}k_{z}^{\prime 2})^{\frac{1}{2}}/v$ and
$k_{y}=v_{z}k_{z}\text{tan}\phi/v$, indicating that the system lies in the
nodal line semimetal phase and the Fermi surface has a torus shape Y.Wang ;
D.Sinha ; while the latter means that the system is a gapped insulator.

Figure 3: (Color online) (a)$-$(b) The conductivities $\sigma_{xx}$ and
$\sigma_{xxx}$ versus the magnetic field $B$ for different linewidths $\eta$.
In (b), the left inset shows that when weak $\eta=0.1$ meV, $\sigma_{xxx}$
exhibits strong oscillations, and the right inset shows the magnetic field
$B_{p}$ and the corresponding conductivity $\sigma_{xxx}(B_{p})$ versus
$\eta$. (c) The MCA coefficient $\gamma^{\prime}$ at $\eta=1$ meV and $2$ meV.
We take the parameters $(\Delta,\xi)=(10,6)$ meV, and the chemical potential
$\mu=30$ meV.

The Landau bands develop when a magnetic field is applied along the $z$
direction, $\bm{B}=B\bm{e}_{z}$. To solve the LLs, we choose the Landau gauge
$\bm{A}=-By\bm{e}_{x}$, which is minimally coupled to the crystal momentum
through the Periels substitution $\bm{p}\rightarrow\bm{p}+e\bm{A}$. Then we
replace the momentum operators in $\hat{H}_{0}$ by the standard ladder
operators, $k_{+}\rightarrow\frac{\sqrt{2}}{l_{B}}\hat{a}^{\dagger}$ and
$k_{-}\rightarrow\frac{\sqrt{2}}{l_{B}}\hat{a}$, with the magnetic length
$l_{B}=\frac{1}{\sqrt{eB}}$. The total Hamiltonian $\hat{H}_{0}+\hat{H}_{IB}$
under a magnetic field is written as

$\displaystyle\hat{H}_{B}$
$\displaystyle=\begin{pmatrix}m&P\hat{a}+\Delta&0&v_{z}k_{z}-i\xi\\\
P\hat{a}^{\dagger}+\Delta&-m&v_{z}k_{z}+i\xi&0\\\
0&v_{z}k_{z}-i\xi&m&-P\hat{a}^{\dagger}+\Delta\\\
v_{z}k_{z}+i\xi&0&-P\hat{a}+\Delta&-m\end{pmatrix},$ (12)

in which the parameter $P=\frac{\sqrt{2}}{l_{B}}v$.

When $\Delta=0$, the LLs can be solved analytically. For the $n\geq 1$ LLs, we
use the trial wave function
$\psi_{n}=\big{(}c_{n1}|n-1\rangle,c_{n2}|n\rangle,c_{n3}|n\rangle,c_{n4}|n-1\rangle\big{)}^{T}$,
with $|n\rangle$ being the harmonic oscillator state defined as
$\hat{a}^{\dagger}\hat{a}|n\rangle=n|n\rangle$. The LL energies are obtained
as

$\displaystyle\varepsilon_{ns\lambda}(k_{z})=s\Big{[}v_{z}^{2}k_{z}^{2}+m^{2}+\big{(}P\sqrt{n}+\lambda\xi\big{)}^{2}\Big{]}^{\frac{1}{2}}.$
(13)

For the zeroth LLs, we use the trial wave function
$\psi_{0}=\big{(}0,c_{02}|0\rangle,c_{03}|0\rangle,0\big{)}^{T}$. The energies
are given as

$\displaystyle\varepsilon_{0s}(k_{z})=s\big{(}v_{z}^{2}k_{z}^{2}+m^{2}+\xi^{2}\big{)}^{\frac{1}{2}},$
(14)

which remain unchanged to the magnetic field.

The LL spectra are displayed in Figs. 2(a)$-$(c). When $\xi=0$, the
$\lambda=\pm$ branches are degenerate [Fig. 2(a)]. The nonvanishing $\xi$ will
break the twofold degeneracy of the $n\geq 1$ LLs [Fig. 2(b)]. At weak
magnetic fields, when $\varepsilon_{n^{\prime}++}>\varepsilon_{n+-}$ with
$n^{\prime}<n$, the LL crossings at higher energies are clearly seen, leading
to the separation of the $ns\pm$ LLs. In the quantum Hall measurements, the LL
crossings were also caused by the broken inversion symmetry and can induce
complicated filling factors with the varying carrier density and magnetic
field A.C.Lygo ; D.A.Kealhofer . At strong magnetic fields with
$P\sqrt{n}\gg\xi$, the $ns\pm$ LLs will rearrange and pair with each other.

##  IV Selection Rules

More importantly, when $\Delta=0$, the Hamiltonian $H_{B}$ owns the parity
symmetry with the operator $\hat{\cal
P}=(-1)^{\hat{a}^{\dagger}\hat{a}}\sigma_{z}\otimes\tau_{z}$ T.Devakul ;
Y.X.Wang2022 . This means that each LL carries a definite even or odd parity,
as represented by the red solid or blue dashed lines in Figs. 2(a) and 2(b).
For the velocity operator $\hat{v}_{x}=v\sigma_{z}\otimes\tau_{x}$, because
$\hat{\cal P}^{-1}\hat{v}_{x}\hat{\cal P}=-\hat{v}_{x}$, it has an odd parity.
Thus the nonvanishing velocity operator matrix element
$\langle\psi_{n}|\hat{v}_{x}|\psi_{n^{\prime}}\rangle$ requires that the
initial and final states own opposite parities, leading to the selection rule
$n\rightarrow n\pm 1$.

When $\Delta$ is finite, $\hat{H}_{B}$ needs to be solved with numerics: in
the infinite Hilbert space spanned by the harmonic oscillator state
$|n\rangle$, the LLs are obtained by diagonalizing $H_{B}$ after truncating
the Hilbert space at the cutoff $N_{c}$. In the calculations, we set
$N_{c}=200$ to obtain the well-convergent results. Now since $\hat{\cal
P}^{-1}(\Delta\tau_{x})\hat{\cal P}=-\Delta\tau_{x}$, the parity symmetry of
the system is broken and the LLs will not carry a definite parity. With a
finite $\Delta=10$ meV, the LL crossings remain at higher energies and their
movements are clearly seen [Fig. 2(c)]. Note that as $\hat{\cal
P}^{-1}(\xi\sigma_{x}\otimes\tau_{y})\hat{\cal
P}=\xi\sigma_{x}\otimes\tau_{y}$, a finite $\xi$ does not break the parity
symmetry of the system.

To find out the selection rules, in Fig. 1(d), we display the velocity
operator matrix element
$|\langle\psi_{ns\lambda}|\hat{v}_{x}|\psi_{n^{\prime}s^{\prime}\lambda^{\prime}}\rangle|$
between the initial state $|\psi_{ns\lambda}\rangle$ and the final state
$|\psi_{n^{\prime}s^{\prime}\lambda^{\prime}}\rangle$ for different $\Delta$.
The magnetic field is set as $B=1$ T and the initial state is chosen as
$(3--)$. We see that when $\Delta=0$, the matrix elements are nonvanishing
only for the final states $(2+\pm)$ and $(4+\pm)$, from which the selection
rules are inferred as the conventional ones $n\rightarrow n\pm 1$ Y.Jiang ;
Y.X.Wang2023 ; Z.Cai  and agree with the above parity symmetry analysis. When
$\Delta=2$ meV and the parity of each LL is broken, we find that besides
$n\rightarrow n\pm 1$, the selection rules can take $n\rightarrow n$ as well
as $n\rightarrow n\pm 2$. Further increasing $\Delta$, the LL eigenstate
$|\psi_{n}\rangle$ may involve the state $|n^{\prime}\rangle$ that is far away
from $|n\rangle$, leading to more selection rules. When $\Delta=20$ meV, the
selection rules can even take $n\rightarrow n+4$. Therefore, we suggest that
the broken parity symmetry by the $\Delta$ term can lead to the mixed
selection rules, $n\rightarrow n,n\pm 1,n\pm 2,\cdots$, which makes
$\sigma_{xxx}$ nonvanishing.

##  V MCA coefficient $\gamma^{\prime}$

Figure 4: (Color online) (a)$-$(b) The conductivities $\sigma_{xx}$ and
$\sigma_{xxx}$ versus $\Delta$ for different linewidths $\eta$. The inset in
(b) shows that when weak $\eta=0.1$ meV, $\sigma_{xxx}$ exhibits strong
oscillations. (c) The MCA coefficient $\gamma^{\prime}$ versus $\Delta$ at
$\eta=1$ meV and $2$ meV. We take the parameter $\xi=6$ meV, the magnetic
field $B=0.2$ T, and the chemical potential $\mu=30$ meV.

To obtain the MCA coefficient $\gamma^{\prime}$, we need to calculate the
conductivities $\sigma_{xx}$ and $\sigma_{xxx}$ by using Eqs. (7) and (8). The
results are plotted in Figs. 3(a) and 3(b) as functions of the magnetic field
$B$ for different linewidths $\eta$, in which the chemical potential $\mu$ is
fixed at $\mu=30$ meV to ensure that the 3D system owns a low carrier density
F.Tang ; E.Martino .

Firstly, we study $\sigma_{xx}$. In the limiting clean case of $\eta=0.1$ meV,
$\sigma_{xx}$ exhibits a series of oscillation peaks that sit on a descending
background. This is because under a weak magnetic field, the LLs are densely
distributed. With increasing $B$, the LLs will cross the chemical potential
$\mu$ one by one, resulting in the oscillation peaks of $\sigma_{xx}$
D.Shoenberg ; Y.X.Wang2023 . Meanwhile, the decreasing LL number around $\mu$
leads to the descending background. When $B>2.44$ T, the system enters the
extreme quantum limit, with all electrons confined to the zeroth LL Z.Cai ;
S.Galeski2022 ; W.Wu . As $\eta$ grows, in the quantum limit, $\sigma_{xx}$
increases; while in the quantum oscillation regime, $\sigma_{xx}$ exhibits a
nonmonotonous variation, which are consistent with the previous studies
Y.X.Wang2020 ; Y.X.Wang2023 .

Next we study $\sigma_{xxx}$. When $\eta=0.1$ meV, $\sigma_{xxx}$ exhibits
strong oscillations [the left inset of Fig. 3(b)], but will quickly drop to
zero at a relatively large magnetic field. The former behavior can also be
explained by the densely distributed LLs, while the latter is attributed to
the enlarged LL spacings. When $\eta$ grows, the neighboring oscillations will
merge together. For weak impurity scatterings, e.g., $\eta=1$ meV and $2$ meV,
we see that $\sigma_{xxx}$ first increases with $B$; after reaching its peak
value at the magnetic field $B_{p}$, $\sigma_{xxx}$ then decreases to zero. To
see this behavior more deeply, in the right inset of Fig. 3(b), we plot
$B_{p}$ and the corresponding conductivity $\sigma_{xxx}(B_{p})$ as functions
of $\eta$. We observe that $B_{p}$ increases with $\eta$ and is saturated when
$\eta\geq 4$ meV; while $\sigma_{xxx}$ decreases steadily to zero. Such a
dependence of $\sigma_{xxx}$ on the magnetic field is consistent with the
second-order longitudinal resistance observations in the experiment N.Wang ,
where temperature plays a role similar to the impurity scatterings. Further
increasing $\eta$ to $\eta\geq 5$ meV, $\sigma_{xxx}$ is completely
suppressed. These results tell us that $\sigma_{xxx}$ is evident only at a
weak magnetic field, and is susceptible to the impurity scatterings.

After obtaining $\sigma_{xx}$ and $\sigma_{xxx}$, we calculate the MCA
coefficient $\gamma^{\prime}$ by using Eq. (3). From the above analysis, we
see that $\gamma^{\prime}$ is dependent on the impurity scatterings, which is
different from the conclusion of SBE. In the framework of the SBE Y.Wang , it
was suggested that $\sigma_{xx}\sim\tau$ and $\sigma_{xxx}\sim\tau^{2}$, with
the relaxation time equal to the inverse linewidth, $\tau=\eta^{-1}$, leading
to the independence of $\gamma^{\prime}$ on $\eta$. Fig. 3(c) plots
$\gamma^{\prime}$ at $\eta=1$ meV and $2$ meV, for which $\sigma_{xxx}$ takes
a relatively large value. We observe that (i) $\gamma^{\prime}$ exhibits a
decreasing trend with the magnetic field $B$, which captures the qualitative
trend of the experimental results Y.Wang ; N.Wang ; and (ii) $\gamma^{\prime}$
reaches the order of $10^{-10}$ T-1A-1m2. Compared with the experimental value
$\gamma^{\prime}\sim 10^{-7}$ T-1A-1m2 Y.Wang ; N.Wang , although there is
still a large discrepancy, it has improved one order of magnitude over the SBE
result $\gamma^{\prime}\sim 10^{-11}$ T-1A-1m2 Y.Wang .

Since the $\Delta$ term plays a decisive role in determining the selection
rules, we study its influence on MCA. With the magnetic field set as $B=0.2$
T, in Figs. 4(a)$-$4(c), the conductivities, $\sigma_{xx}$ and $\sigma_{xxx}$,
and the MCA coefficient $\gamma^{\prime}$ are displayed as functions of
$\Delta$, respectively. In Fig. 4(a), when $\eta=0.1$ meV, $\sigma_{xx}$
exhibits certain oscillations with $\Delta$, which can be attributed to the
$\Delta$ term-driven LL crossing over the chemical potential. When $\eta\geq
1$ meV, the oscillation are smeared by the impurity scatterings. We also see
that $\sigma_{xx}$ increases slowly with $\Delta$, due to the almost unchanged
LL number around the chemical potential.

In Fig. 4(b) when $\eta=1$ meV, $\sigma_{xxx}$ shows oscillations that lie on
an increasing background. Compared to $\sigma_{xx}$, more oscillations are
visible in $\sigma_{xxx}$. When $\Delta=0$, the parity symmetry of the system
is preserved and the LLs carry definite odd or even parities Y.X.Wang2022 ;
T.Devakul . Correspondingly, the selection rules take the conventional ones
$n\rightarrow n\pm 1$ and thus $\sigma_{xxx}$ vanishes. At a finite $\Delta$,
the parity symmetry is broken and the resulting selection rules are mixed,
which may include $n\rightarrow n,n\pm 1,n\pm 2,\cdots$, and lead to a finite
$\sigma_{xxx}$. With increasing $\Delta$, on one hand, besides the initial and
final states, the middle state $|\psi_{n^{\prime\prime}}\rangle$ can also
cross the chemical potential [see Eq. (8)], resulting in more oscillations in
$\sigma_{xxx}$; on the other hand, since more selection rules are allowed,
$\sigma_{xxx}$ shows a pronounced enhancement with $\Delta$. When $\eta=2$
meV, $\sigma_{xxx}$ increases linearly with $\Delta$, and the oscillations
almost disappear.

In Fig. 4(c), we observe that when $\eta=1$ meV, $\gamma^{\prime}$ also
oscillates with $\Delta$ on a linear increasing background. When $\eta=2$ meV,
$\gamma^{\prime}$ increases linearly with $\Delta$,

$\displaystyle\gamma^{\prime}\sim\Delta.$ (15)

Such a linear dependence is distinct from the SBE results Y.Wang , where
$\gamma^{\prime}$ was found to be inversely proportional to $\Delta$,
$\gamma^{\prime}\sim\Delta^{-1}$. Here as the increasing $\Delta$ can induce
more selection rules, more electronic states that are far away from the Fermi
surface can paticipate in the magnetotransport and correspondingly,
$\gamma^{\prime}$ will increase. Thus we suggest that the linear relationship
is more physically reasonable.

##  VI Discussions and Conclusions

To summarize, in this paper, by developing a full quantum theory, we calculate
the MCA coefficient $\gamma^{\prime}$ in 3D ZrTe5 and explore the effects of
the magnetic field, the impurity scatterings and the mirror symmetry breaking.
We reveal the significance of the $xz$-mirror symmetry breaking in forming the
mixed selection rules, which in turn determines $\sigma_{xxx}$ and
$\gamma^{\prime}$. The results show that $\gamma^{\prime}$ can capture the
qualitative features of the experiments Y.Wang ; N.Wang , and exhibits a
linear dependence on the strength of the $xz$-mirror symmetry breaking.

Since the full quantum theory is based on the LL energies and wave functions,
it goes beyond the previous SBE framework and the Drude conductivity formula.
We expect that the derived $\sigma_{xxx}$ can be used to analyze the recently
reported MCA in the Weyl semimetal WTe2 T.Yokouchi , and the formula of
$\sigma_{\alpha\beta\gamma}$ can be further exteneded to study the nonlinear
Hall effect Q.Ma ; K.Kang .

##  VII Acknowledgments

This work was supported by the Natural Science Foundation of China (Grant No.
11804122, No. 11905054, No. 12275075), the National Key Research and
Development Program of Ministry of Science and Technology (Grant No.
2021YFA1200700), the Fundamental Research Funds for the Central Universities
from China, and the China Postdoctoral Science Foundation (Grant No.
2021M690970).

## References

  * (1) F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * (2) C. Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L. L. Wang, Z. Q. Ji, Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S. C. Zhang, K. He, Y. Wang, L. Lu, X. C. Ma, and Q. K. Xue, Science 340, 167 (2013).
  * (3) C. Z. Chang, C. X. Liu, and A. H. MacDonald, Rev. Mod. Phys. 95, 011002 (2023).
  * (4) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (5) X. L. Qi and S. C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (6) Y. Tokura and N. Nagaosa, Nat. Commun. 9, 3740 (2018).
  * (7) G. L. J. A. Rikken, J. Fölling, and P. Wyder, Phys. Rev. Lett. 87, 236602 (2001).
  * (8) G. L. J. A. Rikken and P. Wyder, Phys. Rev. Lett. 94, 016601 (2005).
  * (9) T. Ideue, K. Hamamoto, S. Koshikawa, M. Ezawa, S. Shimizu, Y. Kaneko, Y. Tokura, N. Nagaosa and Y. Iwasa, Nat. Phys. 13, 578 (2017).
  * (10) G. L. J. A. Rikken and N. Avarvari, Phys. Rev. B 99, 245153 (2019).
  * (11) Y. Wang, H. F. Legg, T. Bömerich, J. Park, S. Biesenkamp, A. A. Taskin, M. Braden, A. Rosch, and Y. Ando, Phys. Rev. Lett. 128, 176602 (2022).
  * (12) N. Wang, J. Y. You, A. Wang, X. Zhou, Z. Zhang, S. Lai, H. J. Tien, T. R. Chang, Y. P. Feng, H. Lin, G. Chang and W. B Gao, arXiv: 2202.07143 (2022).
  * (13) P. He, C. H. Hsu, S. Shi, K. Cai, J. Wang, Q. Wang, G. Eda, H. Lin, V. M. Pereira, and H. Yang, Nat. Commun. 10, 1290 (2019).
  * (14) H. Watanabe and Y. Yanase, Phys. Rev. X 11, 011001 (2021).
  * (15) C. Wang, Y. Gao, and D. Xiao, Phys. Rev. Lett. 127, 277201 (2021).
  * (16) I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * (17) T. Morimoto, S. Zhong, J. Orenstein, and J. E. Moore, Phys. Rev. B 94, 245121 (2016).
  * (18) Z. Z. Du, C. M. Wang, H. Z. Lu, and X. C. Xie, Phys. Rev. Lett. 121, 266601 (2018).
  * (19) M. X. Deng, G. Y. Qi, R. Ma, R. Shen, R. Q. Wang, L. Sheng, and D. Y. Xing, Phys. Rev. Lett. 122, 036601 (2019).
  * (20) Y. Gao and F. Zhang, Phys. Rev. B 103, L041301 (2021).
  * (21) D. Sinha and A. Taraphder, Phys. Rev. B 104, 245141 (2021).
  * (22) H. Weng, X. Dai, and Z. Fang, Phys. Rev. X 4, 011002 (2014).
  * (23) R. Y. Chen, Z. G. Chen, X. Y. Song, J. A. Schneeloch, G. D. Gu, F. Wang, and N. L. Wang, Phys. Rev. Lett. 115, 176404 (2015).
  * (24) Z. G. Chen, R. Y. Chen, R. D. Zhong, J. Schneeloch, C. Zhang, Y. Huang, F. Qu, R. Yu, Q. Li, G. D. Gu, and N. L. Wang, Proc. Natl. Acad. Sci. USA 114, 816 (2017).
  * (25) F. Tang, Y. Ren, P. Wang, R. Zhong, J. Schneeloch, S. A. Yang, K. Yang, P. A. Lee, G. Gu, Z. Qiao, and L. Zhang, Nature 569, 537 (2019).
  * (26) S. Galeski, T. Ehmcke, R. Wawrzyńczak, P. M. Lozano, K. Cho, A. Sharma, S. Das, F. Küster, P. Sessi, M. Brando, R. Küchler, A. Markou, M. König, P. Swekis, C. Felser, Y. Sassa, Q. Li, G. Gu, M. V. Zimmermann, O. Ivashko, D. I. Gorbunov, S. Zherlitsyn, T. Förster, S. S. P. Parkin, J. Wosnitza, T. Meng, and J. Gooth, Nat. Commun. 12, 3197 (2021).
  * (27) L. Smrcka and P. Streda, J. Phys. C: Solid State Phys. 10, 2153 (1977).
  * (28) G. D. Mahan, Many-Particle Physics, 3rd ed. (Plenum, New York, 2000).
  * (29) See Supplemental Material for detailed calculations of the density matrix and conductivities.
  * (30) H. W. Wang, B. Fu, and S. Q. Shen, Phys. Rev. B 98, 081202(R) (2018).
  * (31) Y. X. Wang and F. Li, Phys. Rev. B 101, 085201 (2020).
  * (32) Y. Jiang, J. Wang, T. Zhao, Z. L. Dun, Q. Huang, X. S. Wu, M. Mourigal, H. D. Zhou, W. Pan, M. Ozerov, D. Smirnov, and Z. Jiang, Phys. Rev. Lett. 125, 046403 (2020).
  * (33) Y. X. Wang and Z. Cai, Phys. Rev. B 107, 125203 (2023).
  * (34) Z. Cai and Y. X. Wang, Phys. Rev. B 108, 155202 (2023).
  * (35) Y. X. Wang and F. Li, Phys. Rev. B 106, 205102 (2022).
  * (36) A. C. Lygo, B. Guo, A. Rashidi, V. Huang, P. Cuadros-Romero, and S. Stemmer, Phys. Rev. Lett. 130, 046201 (2023).
  * (37) D. A. Kealhofer, L. Galletti, T. Schumann, A. Suslov, and S. Stemmer, Phys. Rev. X 10, 011050 (2020).
  * (38) T. Devakul, Y. H. Kwan, S. L. Sondhi, and S. A. Parameswaran, Phys. Rev. Lett. 127, 116602 (2021).
  * (39) E. Martino, I. Crassee, G. Eguchi, D. Santos-Cottin, R. D. Zhong, G. D. Gu, H. Berger, Z. Rukelj, M. Orlita, C. C. Homes, and A. Akrap, Phys. Rev. Lett. 122, 217402 (2019).
  * (40) D. Shoenberg, Magnetic Oscillations in Metals (Cambridge University Press, Cambridge, 1984).
  * (41) S. Galeski, H. F. Legg, R. Wawrzyńczak, T. Förster, S. Zherlitsyn, D. Gorbunov, M. Uhlarz, P. M. Lozano, Q. Li, G. D. Gu, C. Felser, J. Wosnitza, T. Meng, and J. Gooth, Nat. Commun. 13, 7418 (2022).
  * (42) W. Wu, Z. Shi, Y. Du, Y. Wang, F. Qin, X. Meng, B. Liu, Y. Ma, Z. Yan, M. Ozerov, C. Zhang, H. Z. Lu, J. Chu, and X. Yuan, Nat. Mat. 22, 84 (2023).
  * (43) T. Yokouchi, Y. Ikeda, T. Morimoto, and Y. Shiomi, Phys. Rev. Lett. 130, 136301 (2023).
  * (44) Q. Ma, S. Y. Xu, H. Shen, D. MacNeill, V. Fatemi, T. R. Chang, Andrés M. M. Valdivia, S. Wu, Z. Du, C. H. Hsu, S. Fang, Q. D. Gibson, K. Watanabe, T. Taniguchi, R. J. Cava, E. Kaxiras, H. Z. Lu, H. Lin, L. Fu, N. Gedik, and P. Jarillo-Herrero, Nature 565, 337 (2019).
  * (45) K. Kang, T. Li, E. Sohn, J. Shan and K. F. Mak, Nat. Mat. 18, 324 (2019).


# A Fully Memristive Spiking Neural Network with Unsupervised Learning

Peng Zhou, Dong-Uk Choi, Jason K. Eshraghian,
and Sung-Mo Kang P. Zhou, D. Choi, and S. M. Kang are with the Department of
Electrical and Computer Engineering, UC, Santa Cruz, CA, USA.J. K. Eshraghian
is with the Dept. of Electrical Engineering and Computer Science, University
of Michigan, MI, USA and the Dept. of Computer Science and Software
Engineering, University of Western Australia, Australia.

###### Abstract

We present a fully memristive spiking neural network (MSNN) consisting of
physically-realizable memristive neurons and memristive synapses to implement
an unsupervised Spike Timing Dependent Plasticity (STDP) learning rule. The
system is fully memristive in that both neuronal and synaptic dynamics can be
realized by using memristors. The neuron is implemented using the SPICE-level
memristive integrate-and-fire (MIF) model, which consists of a minimal number
of circuit elements necessary to achieve distinct depolarization,
hyperpolarization, and repolarization voltage waveforms. The proposed MSNN
uniquely implements STDP learning by using cumulative weight changes in
memristive synapses from the voltage waveform changes across the synapses,
which arise from the presynaptic and postsynaptic spiking voltage signals
during the training process. Two types of MSNN architectures are investigated:
1) a biologically plausible memory retrieval system, and 2) a multi-class
classification system. Our circuit simulation results verify the MSNN’s
unsupervised learning efficacy by replicating biological memory retrieval
mechanisms, and achieving 97.5$\%$ accuracy in a 4-pattern recognition problem
in a large scale discriminative MSNN.

###### Index Terms:

Neuromorphic computing, memristor, spiking neural network, unsupervised
learning, STDP.

##  I Introduction

Neuromorphic computing is guided by the rich neural dynamics present in the
brain, the highly parallelized nature of neural computation, and the sparse
encoding of data as spikes, in pursuit of optimizing memory and computation
for energy efficiency. The pervasive von Neumann architecture disaggregates
memory and computation which leaves much room for improvement for threads with
a deterministic set of instructions. This deficiency has spurred the
development of a variety of neuromorphic computing systems [1, 2, 3, 4, 5, 6,
7, 8, 9], which integrate spiking neurons and simplified synaptic models onto
a silicon substrate. In almost all instances, synaptic weights are stored in
random access memory (RAM), thus moving memory closer to processor. But memory
access and computation remain as two separate steps, which does not address
the cost of data communication and memory access, which impose the most
overhead in both neuromorphic and general purpose computing systems.

The memristor has been presented as an option to merging the computation and
memory substrates [10, 11, 12]. What was once a theoretical postulate proposed
by Chua in 1971 [10] has become a commercially available technology that can
be integrated in the back-end-of-the-line (BEOL) of modern CMOS processes [13,
14]. Their non-volatile retention capacity is often likened to synapses [15,
16, 17, 18, 19], and their threshold-switching characteristics occasionally
exploited as a spiking neuron model [20, 21, 22, 23]. For example, the
modulation of device resistance has been correlated to synaptic plasticity,
where achieving short-term plasticity (STP) and long-term plasticity (LTP)
using memristors is as simple as applying programming pulse trains [24]. The
benefits of memristive neurons and synapses arise from CMOS-compatibility,
high-density, nanoscale vertical integration, and their ability to directly
implement biological features, as opposed to requiring several arithmetic
steps as with transistor-only circuits [25]. The use of spiking neural
networks (SNNs) as opposed to modern deep learning paradigms has shown
significant energy and latency benefits as a result of activation sparsity,
spike-based representations of data, and event-based data processing [26, 27,
28, 29, 30], and can be harnessed using the growing infrastructure to support
SNN simulations [31, 32, 33, 34, 35].

Much of the prior work on memristive spiking neural networks (MSNNs) is
constrained to using either memristive synapses or memristive neurons. The
work in [36] demonstrated a fully integrated memristive system that emulated
both synapses and spiking neuron models to achieve simple pattern recognition
tasks. The result is demonstrated using an in-house fabricated diffusive
memristor, which is not readily accessible to the broader research community.
In this work, we introduce a fully memristive approach to designing a system
that shows the capacity to learn pattern recognition tasks using memristive
neurons and synapses. We achieve this using models of commercially available,
low-cost memristors [37]. We use a memristive neuron circuit which consists of
the minimum number of circuit elements necessary to achieve distinct
depolarization, hyperpolarization, and repolarization voltage phases. This
neuron is integrated in a multi-layer MSNN with memristive synaptic
interconnections, modelled using an equivalent SPICE model of a memristive
array to train the MSNN with the STDP learning rule [38]. We successfully
demonstrate the potential of MSNN to perform associative memory retrieval, as
well as pattern recognition achieving 97.5% accuracy on a multi-class
classification task. Our results validate homogeneous memristive systems that
are capable of neuromorphic learning without global error signals.

##  II Methods

###  II-A MIF model

In this paper, we apply the MIF neuron model [39]. The MIF circuit is shown in
Fig. 1. The proposed MIF circuit is characterized by the differential
equations below:

Figure 1: Pre-synaptic and post-synaptic memristive neurons with a trainable
memristive synapse interposed between the two. The neuron model is a
memristive integrate-and-fire (MIF) neuron circuit, consisting of two
memristors $M_{1}$ and $M_{2}$, connected to DC voltage sources $E_{\rm rest}$
and $E_{\rm reset}$, in parallel with a capacitor $C$. Voltage spikes
generated by the MIF neuron propagates through the synapse, and trigger an
input current to the post-synaptic MIF neuron, which in turn will generate
spikes.

$\frac{{\mathrm{d}}v}{\mathrm{d}t}=\frac{I-G_{1}(v-E_{rest})-G_{2}(v-E_{reset})}{C}$
(1a)
$\frac{{\mathrm{d}x_{1}}}{\mathrm{d}t}=\frac{1}{\tau_{1}}(\frac{1-x_{1}}{1+e^{\frac{v_{on1}-(v-E_{rest})}{k_{th}}}}-\frac{x_{1}}{1+e^{\frac{(v-E_{rest})-v_{off1}}{k_{th}}}})$
(1b)
$\frac{{\mathrm{d}x_{2}}}{\mathrm{d}t}=\frac{1}{\tau_{2}}(\frac{1-x_{2}}{1+e^{\frac{v_{on2}-(v-E_{rest})}{k_{th}}}}-\frac{x_{2}}{1+e^{\frac{(v-E_{rest})-v_{off2}}{k_{th}}}})$
(1c) $G_{1}=\frac{x_{1}}{R_{on1}}+\frac{1-x_{1}}{R_{off1}}$ (1d)
$G_{2}=\frac{x_{1}}{R_{on2}}+\frac{1-x_{2}}{R_{off2}}$ (1e)

where $G_{1}$ and $G_{2}$ are the memductances, $x_{1}$ and $x_{2}$ are a pair
of internal states, $\tau_{1}$ and $\tau_{2}$ are time constants governing the
rate of change in internal states, $k_{th}$ is the effective thermal voltage.
They are the characteristic variables of $M_{1}$ and $M_{2}$, respectively.
This system of equations mirrors several prominent SPICE memristor models, and
has been used to emulate the commercially available Knowm memristor [37].

To show the dynamics of the above system, we apply an alpha-shaped input
current, with the resultant memristor internal states and voltage waveforms
shown below in Fig. 2:

Figure 2: Simulation result of MIF receiving an alpha shape input current.
(a) Internal states $x_{1}$, $x_{2}$. (b) Voltage response $v$.

The alpha current is modeled by Eq. (2):

$\displaystyle\tau_{syn}\frac{{\mathrm{d}I}}{\mathrm{d}t}$ $\displaystyle=a-I$
(2a) $\displaystyle\tau_{syn}\frac{{\mathrm{d}a}}{\mathrm{d}t}$
$\displaystyle=-a+W_{j}\cdot\sum_{f}{\delta(t-t_{j}^{f})}$ (2b)

where $W_{j}$ is a weighted synaptic current generated between presynaptic
neuron $j$ and associated postsynaptic neuron. $\sum_{f}{\delta(t-t_{j}^{f})}$
indicates the total number of spikes emitted by presynaptic neuron $j$, and
$I$ is the input current.

The parameters used in this network are listed in Table I:

TABLE I: MIF circuit parameters Parameter | Value | Parameter | Value
---|---|---|---
$E_{rest}$ | 0 mV | $E_{reset}$ | 50 mV
$C$ | 100 pF | $k_{th}$ | 15 mV
$v_{off_{1}}$, $v_{off_{2}}$ | 0 mV |  $v_{on_{1}}$, $v_{on_{2}}$ | 100 mV
$R_{off_{1}}$, $R_{off_{2}}$ | 0.1 M$\Omega$ |  $R_{on_{1}}$, $R_{on_{2}}$ | 1 k$\Omega$
$\tau_{1}$, $\tau_{2}$ | 1 $\mu$s |  |

###  II-B Memristive Synapse with STDP

The Spike-timing-dependent plasticity (STDP) learning rule modulates synaptic
weights based on the time difference of pre- and postsynaptic spike arrivals
[38]. This type of learning rule can be be achieved by using memristors [40,
18], and can be verified using SPICE-level models [41]. As Fig. 1 shows, a
synaptic memristor is interposed between two neurons, where the pre- and
postsynaptic spikes will generate a voltage across the memristor that causes
the memristance to be updated. Moreover, the time difference between pre- and
postsynaptic spikes will modulate the changes in the memristance of the
synapse.

In the memristive synapse with STDP learning mechanism, the weight change will
increase rapidly and then decrease slowly when the time difference of a pre-
and post- synaptic spike increases from zero. The weight change of a
memristive synapse is determined by the voltage across the memristor. The
largest voltage will result in the largest weight change, and occurs when one
of the pre- or post- synaptic neuron reaches the threshold level, while the
other is at the reset voltage level. Therefore, the learning window of
memristive STDP will rise with a fast time constant, followed by a slow decay.
Thus, it is reasonable to model the memristive synapse behavior with an alpha
shape learning window as:

$\displaystyle W_{pre}(x)$
$\displaystyle=U_{pre}\cdot\frac{\mathinner{\\!\left\lvert\Delta
t\right\rvert}}{\tau_{pre}}\cdot e^{-\frac{\mathinner{\\!\left\lvert\Delta
t\right\rvert}}{\tau_{pre}}}\quad$ (3a) at $\displaystyle\quad
t_{post}\quad\text{for}\quad t_{pre}<t_{post}$ $\displaystyle W_{post}(x)$
$\displaystyle=U_{post}\cdot\frac{\mathinner{\\!\left\lvert\Delta
t\right\rvert}}{\tau_{post}}\cdot e^{-\frac{\mathinner{\\!\left\lvert\Delta
t\right\rvert}}{\tau_{post}}}\quad$ (3b) at $\displaystyle\quad
t_{pre}\quad\text{for}\quad t_{post}<t_{pre}$

where Upre, Upost, $\tau_{\rm pre}$, $\tau_{\rm post}$ are fitted constants,
and $\Delta$t=tpost-tpre. Generally, the curve in the learning window will
take an alpha shape in both the positive and negative planes of the $x$-axis.

This type of STDP can be modeled with a system of differential equations
defined in Eq. (4):

$\tau_{pre}\frac{{\mathrm{d}}A_{pre}}{\mathrm{d}t}=T_{pre}-A_{pre}$ (4a)
$\tau_{pre}\frac{{\mathrm{d}}T_{pre}}{\mathrm{d}t}=T_{pre}+U_{pre}\cdot\sum_{f}{\delta(t-t_{j}^{f})}$
(4b) $\tau_{post}\frac{{\mathrm{d}}A_{post}}{\mathrm{d}t}=T_{post}-A_{post}$
(4c)
$\tau_{post}\frac{{\mathrm{d}}T_{post}}{\mathrm{d}t}=-T_{post}+U_{post}\cdot\sum_{n}{\delta(t-t_{i}^{n})}$
(4d) $W_{j}\leftarrow W_{j}+A_{post}\quad\text{upon presynaptic spike}$ (4e)
$W_{j}\leftarrow W_{j}+A_{pre}\quad\text{upon postsynaptic spike}$ (4f)

where $f$ indicates the total numbers of spikes emitted by the presynaptic
neuron $j$, and $n$ defines the total number of spikes emitted by the
postsynaptic neuron $i$. Upre is typically positive, while Upost is usually
negative. $W_{j}$ is the weighted synaptic current. Upon arrival of a
presynaptic spike, $W_{j}$ immediately increases by the amount of Apost. Tpre
induces a rise by the amount of Upre, which then evolves according to Eq.
(4b). This in turn modulates Apre (4a). Upon arrival of a postsynaptic spike,
$W_{j}$ immediately increases by the amount of Apre. Tpost increases by the
amount of Upost, and then effects Apost in a similar manner to the impact from
a presynaptic spike.

An example of the temporal evolution of the learning rule is shown in Fig. 3,
where an update occurs if presynaptic spike arrival occurs before the
postsynaptic spike. Fig. 3(a) shows Tpre and Tpost instants, at which the
curve takes an immediate increase followed by an exponential decay. Fig. 3(b)
shows Apre and Apost as functions of time, which have alpha-shaped curves, and
Fig. 3(c) shows the fast weight increase when the presynaptic spike occurs
before the postsynaptic spike. When the presynaptic spike arrives at $t=5ms$,
$W_{j}$ increases by the amount of Apost at this time, which is zero. At
$t=10ms$, the postsynaptic spike arrives, $W_{j}$ increases by the amount of
Apre at the this time, which is around 0.3 $\mu A$.

Figure 3: An example of the proposed memristive STDP with $\tau_{\rm pre}$ =
$\tau_{\rm post}$=3 $\mu s$, Upre=1 $\mu A$, Upost=-1 $\mu A$. (a) Tpre and
Tpost in Eq. (4) are determined by a pre- and a post- synaptic spike,
respectively. (b) Apre and Apost in Eq. (4) are determined by a pre- and a
post- synaptic spike, respectively. (c) The weight is updated according to Eq.
(4).

###  II-C Type-1 MSNN for Memory Retrieval

The neural network architecture we implement with the fully memristive neuron
and synapse models consists of three layers [42]. The input layer applies a
Poisson spike train which encodes the input patterns via rate encoding. The
second layer is a MIF excitatory layer. The third layer is an inhibitory layer
and includes feedback connections to the excitatory layer. The architecture is
shown in Fig. 4. The input to the excitatory layer consists of excitatory
synapses with fixed weights between the input layer and the MIF excitatory
layer, with one-to-one connections. The output of the excitatory layer
includes trainable (via STDP) excitatory recurrent synaptic connections, in
addition to fixed-weight excitatory synapses between the MIF excitatory layer
and the inhibitory layer, also with one-to-one connections. The inhibitory
layer includes inhibitory feedback synapses with fixed weights between the
inhibitory layer and the MIF excitatory layer with one-to-(all-1) connections
(not including the single MIF neuron connected to the inhibitory neuron). In
this neural network, the number of neurons in all three layers are identical,
and STDP is enabled within the MIF excitatory layer. The input patterns
consist of 32$\times$32 pixels, such that the total number of neurons in each
layer is 1024. Each pattern generates a Poisson spike train of 35 $\mu s$
duration, followed by a 15 $\mu s$ refractory period without any input such
that each input pattern does not affect the network dynamics upon arrival of
the next input pattern. Each input spike will generate an alpha current where
$\tau_{\rm syn}=10ns$. The threshold of each MIF neuron is set to 25 $mV$. As
in the MIF neuron, there is no extra control circuit to force reset, and we
set a refractory of 3 $\mu s$ to prevent duplicating spikes during the
simulation. We find the weight increase has a more significant effect when
compared with the weight decrease, and it also enables the network to learn
faster. Therefore, we set U${}_{\rm pre}=10\mu A$ and U${}_{\rm post}=-0.1\mu
A$. Note that Upre,post refers to an internal state variable that is modulated
by the current-based neuron model, and is in dimensions of amperes. The
recurrent connections are initialized with randomly weighted synaptic currents
between 0 to 0.2 $\mu A$. The exact weight values between excitatory layer and
inhibition layer do not have a large effect, which enables the inhibitory
neurons to trigger a spike following an output from excitatory neurons. On the
other hand, the weights between input-to-excitatory layers and the inhibitory
feedback mechanism both need to be carefully tuned to be neither too weak, nor
too strong, to prevent complete suppression of downstream spikes, or excessive
reinforcement of firing. We chose the fixed input to the excitatory weight to
be 50 $\mu A$, and the fixed inhibitory layer feedback weight to be 20 $\mu
A$.

Figure 4: The memristive spiking neural network architecture, which consists
of the input layer, MIF excitatory layer, and the inhibitory layer. The blue
arrows denote the excitatory synaptic forward connections, the orange arrow
indicates excitatory synaptic recurrent connections, and the green arrow shows
the inhibitory synaptic forward connection.

###  II-D Type-2 MSNN for Pattern Recognition

The sequence of the MSNN for pattern recognition is identical to the memory
retrieval network, in that we use a Poisson input layer, an excitatory layer,
and an inhibitory layer. To improve unsupervised learning we modify the nature
of the synaptic connections. STDP-based learning is enabled between the input
layer and the excitatory MIF layer with all-to-all connections. As before,
each excitatory neuron connects to one inhibitory neuron and the inhibitory
layer uses inhibitory synapses with fixed weights connected to the excitatory
layer. In this network, no recurrent connections are present due to the lack
of temporal dependencies in the pattern recognition task. The number of
neurons between the excitatory and inhibitory layers must be identical, though
they no longer need to match the input layer. We use 1024 neurons in the input
layer, and 320 neurons each in the excitatory and inhibitory layers, where
STDP takes place between the input and the excitatory layers. We set U${}_{\rm
pre}=10nA$, U${}_{\rm post}=-0.1nA$, and the fixed inhibitory layer feedback
weight to 200 $\mu A$. The weighted synaptic current between input and
excitatory layers are randomly initialization between 0 and 12 $\mu A$ based
on empirical evaluation. Additionally, a random 0-2 $\mu s$ delay between
input layer and excitatory layer is assigned to avoid excessive simultaneous
spiking, thus mitigating the exploding weights problem. All other neuron and
input parameters are identical to the memory retrieval tasks.

##  III Results

###  III-A Type-1 MSNN for Memory Retrieval Result

We provided four different patterns to the MSNN successively for each
iteration, and simulated across multiple epochs. This is to coarsely emulate
how biological systems experience real-world data in a batch size of ‘1’ in
online learning systems. The Poisson input layer converts each pixel intensity
of the input pattern into a spiking probability of each neuron. We used four
different patterns, as shown in the top row of Fig. 5. The weight matrix is
updated with the STDP rule. The training phase consists of 20 epochs for all
patterns, corresponding to 80 iterations.

Figure 5: Top row: the four input spiking patterns (1:Square, 2: Cross, 3:
Diamond, 4: Triangle) applied to the input of the MSNN. Bottom row: the number
of spikes in the excitatory layer during the test stage, upon receiving these
four input patterns.

After training, the same four patterns are successively applied to the MIF.
The MSNN is expected to ‘recall’ the patterns as the connections have ideally
made associations in the form of synaptic memory. When one of the patterns is
applied to the MSNN again, other patterns are also recalled with less
intensity, which is related to ‘memory retrieval’, illustrated in the bottom
row of Fig. 5. In order to evaluate memory retrieval performance, we calculate
the percentage of overlapping spike count with incorrect patterns, with the
results shown in Fig. 6. We find that when an input pattern is provided to the
network, its resemblance to the spiking behavior with the target behavior is
at a maximum. This is depicted in Fig. 6, where the values along the diagonal
are the largest. As shown in the bottom row of Fig. 5, for each input pattern,
the network continues to recall features of the other three patterns, which is
a result of low-resistance pathways that cause overlap with other patterns,
but with less resemblance to the target pattern.

Figure 6: The resultant heatmap shows the memory retrieval resemblance for
four patterns, which are 1:Square, 2: Cross, 3: Diamond, 4: Triangle.

###  III-B Type-2 MSNN for Pattern Recognition Result

During training, we applied the four different patterns (Fig. 5) to the MSNN
successively, then accumulated and recorded the number of spikes in each
output MIF neuron for each pattern. During the training process, each MIF
neuron is randomly assigned a pattern and updated after every 5 epochs. This
allows each MIF neuron to be associated with a pattern such that it may be
used to encode one of four patterns during the test stage. After assignments
are updated, we use the current weights to test the accuracy of the network.
During the test stage, the four patterns are passed to the network for 40
epochs and the total number of spikes are counted for each pattern for all
output excitatory neurons. The neuron with the highest spike count is deemed
as the winner neuron for this input pattern. The accuracy evolution during
training is shown in Fig. 7, where a total accuracy of $97.5\%$ is attained at
the 70th iteration, shown by the dashed orange line. Early stopping is applied
here, as further training causes excessive reinforcement of high-conductance
pathways.

Figure 7: Accuracy across multiple iterations.

##  IV Discussion and Conclusion

In this paper, we proposed two different types of fully memristive neural
networks for unsupervised learning, based on a memristive integrate-and-fire
neuron model together with memristive synapses. The synapses and neuron models
are designed using SPICE-level memristor models to relate circuit-level
plausibility with the biological plausibility of spiking neurons. We
demonstrated memory retrieval and pattern recognition across the four input
patterns, and validated the potential of our fully-memristive approach across
both tasks. Future work will address more challenging tasks, such as testing
memory retrieval for incomplete patterns, and increasing dataset complexity
for the multi-pattern classification problems. Explorations of new learning
paradigms that rely on error-propagation may need to be integrated together
with the STDP update rule to enable the success of more complex tasks. The
proposed framework will render circuit-level implementations using a broad
class of memristors to perform neuromorphic computing.

## References

  * [1] M. Davies, N. Srinivasa, T.-H. Lin, G. Chinya, Y. Cao, S. H. Choday, G. Dimou, P. Joshi, N. Imam, S. Jain _et al._ , “Loihi: A neuromorphic manycore processor with on-chip learning,” _Ieee Micro_ , vol. 38, no. 1, pp. 82–99, 2018.
  * [2] S. B. Furber, F. Galluppi, S. Temple, and L. A. Plana, “The spinnaker project,” _Proceedings of the IEEE_ , vol. 102, no. 5, pp. 652–665, 2014\.
  * [3] B. V. Benjamin, P. Gao, E. McQuinn, S. Choudhary, A. R. Chandrasekaran, J.-M. Bussat, R. Alvarez-Icaza, J. V. Arthur, P. A. Merolla, and K. Boahen, “Neurogrid: A mixed-analog-digital multichip system for large-scale neural simulations,” _Proceedings of the IEEE_ , vol. 102, no. 5, pp. 699–716, 2014\.
  * [4] P. A. Merolla, J. V. Arthur, R. Alvarez-Icaza, A. S. Cassidy, J. Sawada, F. Akopyan, B. L. Jackson, N. Imam, C. Guo, Y. Nakamura _et al._ , “A million spiking-neuron integrated circuit with a scalable communication network and interface,” _Science_ , vol. 345, no. 6197, pp. 668–673, 2014\.
  * [5] J. Schemmel, D. Brüderle, A. Grübl, M. Hock, K. Meier, and S. Millner, “A wafer-scale neuromorphic hardware system for large-scale neural modeling,” in _2010 IEEE International Symposium on Circuits and Systems (ISCAS)_. IEEE, 2010, pp. 1947–1950.
  * [6] S. Moradi, N. Qiao, F. Stefanini, and G. Indiveri, “A scalable multicore architecture with heterogeneous memory structures for dynamic neuromorphic asynchronous processors (dynaps),” _IEEE transactions on biomedical circuits and systems_ , vol. 12, no. 1, pp. 106–122, 2017.
  * [7] L. A. Camuñas-Mesa, B. Linares-Barranco, and T. Serrano-Gotarredona, “Neuromorphic spiking neural networks and their memristor-CMOS hardware implementations,” _Materials_ , vol. 12, no. 17, p. 2745, 2019.
  * [8] S. Furber, “Large-scale neuromorphic computing systems,” _Journal of neural engineering_ , vol. 13, no. 5, p. 051001, 2016.
  * [9] G. Orchard, E. P. Frady, D. B. D. Rubin, S. Sanborn, S. B. Shrestha, F. T. Sommer, and M. Davies, “Efficient neuromorphic signal processing with loihi 2,” in _2021 IEEE Workshop on Signal Processing Systems (SiPS)_. IEEE, 2021, pp. 254–259.
  * [10] L. O. Chua, “Memristor-the missing circuit element,” _IEEE Transactions on Circuit Theory_ , vol. 18, no. 5, pp. 507–519, 1971.
  * [11] L. O. Chua and S. M. Kang, “Memristive devices and systems,” _Proceedings of the IEEE_ , vol. 64, no. 2, pp. 209–223, 1976.
  * [12] D. B. Strukov, G. S. Snider, D. R. Stewart, and R. S. Williams, “The missing memristor found,” _Nature_ , vol. 453, no. 7191, pp. 80–83, 2008.
  * [13] “eflash,” https://www.tsmc.com/english/dedicatedFoundry/technology/
specialty/eflash.

  * [14] “Mad200,” https://mycmp.fr/nvm-mad200/.
  * [15] S. H. Jo, T. Chang, I. Ebong, B. B. Bhadviya, P. Mazumder, and W. Lu, “Nanoscale memristor device as synapse in neuromorphic systems,” _Nano letters_ , vol. 10, no. 4, pp. 1297–1301, 2010.
  * [16] M. Rahimi Azghadi, Y.-C. Chen, J. K. Eshraghian, J. Chen, C.-Y. Lin, A. Amirsoleimani, A. Mehonic, A. J. Kenyon, B. Fowler, J. C. Lee _et al._ , “Complementary metal-oxide semiconductor and memristive hardware for neuromorphic computing,” _Advanced Intelligent Systems_ , vol. 2, no. 5, p. 1900189, 2020.
  * [17] F. Cai, J. M. Correll, S. H. Lee, Y. Lim, V. Bothra, Z. Zhang, M. P. Flynn, and W. D. Lu, “A fully integrated reprogrammable memristor–CMOS system for efficient multiply–accumulate operations,” _Nature Electronics_ , vol. 2, no. 7, pp. 290–299, 2019.
  * [18] T. Serrano-Gotarredona, T. Masquelier, T. Prodromakis, G. Indiveri, and B. Linares-Barranco, “STDP and STDP variations with memristors for spiking neuromorphic learning systems,” _Frontiers in neuroscience_ , vol. 7, p. 2, 2013.
  * [19] J. K. Eshraghian, X. Wang, and W. D. Lu, “Memristor-based binarized spiking neural networks: Challenges and applications.” _IEEE Nanotechnology Magazine_ , 2022.
  * [20] M. D. Pickett, G. Medeiros-Ribeiro, and R. S. Williams, “A scalable neuristor built with mott memristors,” _Nature materials_ , vol. 12, no. 2, pp. 114–117, 2013.
  * [21] J. K. Eshraghian, K. Cho, C. Zheng, M. Nam, H. H.-C. Iu, W. Lei, and K. Eshraghian, “Neuromorphic vision hybrid RRAM-CMOS architecture,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 26, no. 12, pp. 2816–2829, 2018.
  * [22] X. Zhang, W. Wang, Q. Liu, X. Zhao, J. Wei, R. Cao, Z. Yao, X. Zhu, F. Zhang, H. Lv _et al._ , “An artificial neuron based on a threshold switching memristor,” _IEEE Electron Device Letters_ , vol. 39, no. 2, pp. 308–311, 2017.
  * [23] Y. Zhang, Z. Wang, J. Zhu, Y. Yang, M. Rao, W. Song, Y. Zhuo, X. Zhang, M. Cui, L. Shen _et al._ , “Brain-inspired computing with memristors: Challenges in devices, circuits, and systems,” _Applied Physics Reviews_ , vol. 7, no. 1, p. 011308, 2020.
  * [24] W. Wang, W. Song, P. Yao, Y. Li, J. Van Nostrand, Q. Qiu, D. Ielmini, and J. J. Yang, “Integration and co-design of memristive devices and algorithms for artificial intelligence,” _Iscience_ , p. 101809, 2020.
  * [25] E. Covi, S. Brivio, A. Serb, T. Prodromakis, M. Fanciulli, and S. Spiga, “Analog memristive synapse in spiking networks implementing unsupervised learning,” _Frontiers in neuroscience_ , vol. 10, p. 482, 2016.
  * [26] J. K. Eshraghian, M. Ward, E. Neftci, X. Wang, G. Lenz, G. Dwivedi, M. Bennamoun, D. S. Jeong, and W. D. Lu, “Training spiking neural networks using lessons from deep learning,” _arXiv preprint arXiv:2109.12894_ , 2021\.
  * [27] M. R. Azghadi, C. Lammie, J. K. Eshraghian, M. Payvand, E. Donati, B. Linares-Barranco, and G. Indiveri, “Hardware implementation of deep network accelerators towards healthcare and biomedical applications,” _IEEE Transactions on Biomedical Circuits and Systems_ , vol. 14, no. 6, pp. 1138–1159, 2020.
  * [28] J. K. Eshraghian and W. D. Lu, “The fine line between dead neurons and sparsity in binarized spiking neural networks,” _arXiv preprint arXiv:2201.11915_ , 2022.
  * [29] P. Zhou, J. K. Eshraghian, D.-U. Choi, and S.-M. Kang, “Spiceprop: Backpropagating errors through memristive spiking neural networks,” _arXiv preprint arXiv:2203.01426_ , 2022.
  * [30] P. Zhou and S. Abbaszadeh, “Towards real-time machine learning for anomaly detection,” in _2020 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC)_. IEEE, pp. 1–3.
  * [31] M. Stimberg, R. Brette, and D. F. Goodman, “Brian 2, an intuitive and efficient neural simulator,” _Elife_ , vol. 8, p. e47314, 2019.
  * [32] M.-O. Gewaltig and M. Diesmann, “Nest (neural simulation tool),” _Scholarpedia_ , vol. 2, no. 4, p. 1430, 2007.
  * [33] T. Bekolay, J. Bergstra, E. Hunsberger, T. DeWolf, T. C. Stewart, D. Rasmussen, X. Choo, A. Voelker, and C. Eliasmith, “Nengo: a python tool for building large-scale functional brain models,” _Frontiers in neuroinformatics_ , vol. 7, p. 48, 2014.
  * [34] A. P. Davison, D. Brüderle, J. M. Eppler, J. Kremkow, E. Muller, D. Pecevski, L. Perrinet, and P. Yger, “Pynn: a common interface for neuronal network simulators,” _Frontiers in neuroinformatics_ , vol. 2, p. 11, 2009.
  * [35] M. L. Hines and N. T. Carnevale, “Neuron: a tool for neuroscientists,” _The neuroscientist_ , vol. 7, no. 2, pp. 123–135, 2001.
  * [36] Z. Wang, S. Joshi, S. Savel’ev, W. Song, R. Midya, Y. Li, M. Rao, P. Yan, S. Asapu, Y. Zhuo _et al._ , “Fully memristive neural networks for pattern classification with unsupervised learning,” _Nature Electronics_ , vol. 1, no. 2, pp. 137–145, 2018.
  * [37] T. W. Molter and M. A. Nugent, “The generalized metastable switch memristor model,” in _CNNA 2016; 15th International workshop on cellular nanoscale networks and their applications_. VDE, 2016, pp. 1–2.
  * [38] G.-q. Bi and M.-m. Poo, “Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type,” _Journal of neuroscience_ , vol. 18, no. 24, pp. 10 464–10 472, 1998.
  * [39] S. M. Kang, D. Choi, J. K. Eshraghian, P. Zhou, J. Kim, B.-S. Kong, X. Zhu, A. S. Demirkol, A. Ascoli, R. Tetzlaff, and L. O. Chua, “How to build a memristive integrate-and-fire model for spiking neuronal signal generation,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , 2021.
  * [40] G. Maranhão and J. G. Guimarães, “Low-power hybrid memristor-CMOS spiking neuromorphic STDP learning system,” _IET Circuits, Devices & Systems_, vol. 15, no. 3, pp. 237–250, 2021.
  * [41] C. Yakopcic, T. M. Taha, G. Subramanyam, and R. E. Pino, “Generalized memristive device spice model and its application in circuit design,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 32, no. 8, pp. 1201–1214, 2013.
  * [42] N. Diederich, T. Bartsch, H. Kohlstedt, and M. Ziegler, “A memristive plasticity model of voltage-based STDP suitable for recurrent bidirectional neural networks in the hippocampus,” _Scientific reports_ , vol. 8, no. 1, pp. 1–12, 2018.