Limit behaviour of Eringen’s two-phase elastic beams

Preprint of the article published in
European Journal of Mechanics - A/Solids
89, August–October 2021, 104315

Marzia Sara Vaccaro,
Francesco Paolo Pinnola,
Francesco Marotti de Sciarra,
Raffaele Barretta

https://doi.org/10.1016/j.euromechsol.2021.104315

ⓒ 2021. This manuscript version is made available under the CC-BY-NC-ND 4.0
license http://creativecommons.org/licenses/by-nc-nd/4.0/

###### Abstract

In this paper, the bending behaviour of small-scale Bernoulli-Euler beams is
investigated by Eringen’s two-phase local/nonlocal theory of elasticity.
Bending moments are expressed in terms of elastic curvatures by a convex
combination of local and nonlocal contributions, that is a combination with
non-negative scalar coefficients summing to unity. The nonlocal contribution
is the convolution integral of the elastic curvature field with a suitable
averaging kernel characterized by a scale parameter. The relevant structural
problem, well-posed for non-vanishing local phases, is preliminarily
formulated and exact elastic solutions of some simple beam problems are
recalled. Limit behaviours of the obtained elastic solutions, analytically
evaluated, studied and diagrammed, do not fulfill equilibrium requirements and
kinematic boundary conditions. Accordingly, unlike alleged claims in
literature, such asymptotic fields cannot be assumed as solutions of the
purely nonlocal theory of beam elasticity. This conclusion agrees with the
known result which the elastic equilibrium problem of beams of engineering
interest formulated by Eringen’s purely nonlocal theory admits no solution.

###### keywords:

Purely nonlocal elasticity , Local/nonlocal mixture , Nanobeams , Well-
posedness , Constitutive boundary conditions

##  1 Introduction

Over the past decades, growing interest in nanotechnological applications has
required a deep knowledge of mechanical behavior of micro- and nano-structures
conceived as basic components of smaller and smaller electro-mechanical
devices. Among the wide variety of micro/nano-electro-mechanical systems there
are biosensors (Soukarié et al., 2020), actuators (Lu et al., 2019),
resonators (Chorsi and Chorsi, 2018), DNA switches (Chao et al., 2020), porous
structures (Jankowski et al., 2020; Malikan et al., 2020b), valves/pumps
(Banejad et al., 2020), advanced and composite structures (Sedighi et al.,
2017; Żur et. al, 2020), energy harvesters (Ghayesh and Farokhi, 2020).

Micro/nano-structures exhibit technically significant size effects which
therefore have to be taken into account to accurately model and design new-
generation small-scale devices (Ghayesh and Farajpour, 2019). To this end,
nonlocal continuum mechanics (Polizzotto, 2001, 2002; Bažant and Jirásek,
2002; Borino et al., 2003) has been conveniently exploited by the scientific
community to simulate complex scale phenomena in place of atomistic approaches
which are computationally expensive.

One of the first theories of nonlocal elasticity is the strain-driven model
introduced by Eringen (1983) in the wake of seminal contributions by Rogula
(1965); Kröner et al. (1967); Rogula (1982). According to the strain-driven
model, the nonlocal stress field is the convolution integral between elastic
strain and an averaging kernel depending on a characteristic length.
Originally exploited by Eringen to deal with screw dislocations and surface
waves, when applied to structural problems involving bounded domains,
Eringen’s strain-driven model leads to ill-posed structural problems due to
incompatibility between nonlocal and equilibrium requirements (Romano et al.,
2017a).

An effective strategy to get well-posed problems consists in formulating a new
nonlocal model in which the roles of stresses and strains are swapped (Romano
et al., 2017b). The stress-driven nonlocal approach has been successfully
applied to a wide range of problems of nanotechnological interest (see e.g.
Barretta et al., 2018b; Roghani, 2020; Zhang, 2020; Barretta et al., 2018c;
Pinnola, 2020; Oskouie et al., 2018; Barretta et al., 2019).

Another possible strategy to overcome ill-posedness of strain-driven model
consists in formulating a two-phase local/nonlocal model, first introduced in
(Eringen, 1972, 1987) and recently applied in (Pisano and Fuschi, 2003;
Khodabakhshi and Reddy, 2015; Wang et al., 2016; Eptaimeros et al., 2016;
Barretta et al., 2018a). In particular, the mixture strain-driven model is
based on a convex combination of local and nonlocal phases by means of a
mixture parameter, where the nonlocal phase is the Eringen’s integral
convolution. Based on two parameters, the two-phase theory is able to model a
lot of structural problems. However, ill-posedness is eliminated only for
positive mixture parameters (Romano et al., 2017b).

Indeed, for vanishing local phases, Eringen’s purely nonlocal strain-driven
law is recovered, and hence inconsistencies are expected to occur for
structural nonlocal solutions of the limiting elastic equilibrium problem.

In this regard, inaccurate conclusions on the matter are still present in
current literature (see e.g. Mikhasev and Nobili, 2020; Mikhasev, 2021) by
improperly defining limiting cases of free vibration problems as solutions of
elastodynamic problems based on purely nonlocal strain-driven continuum model
which admits no solution (Romano et al., 2017a; Fernández-Sáez and Zaera,
2017; Fathi and Ghassemi, 2017; Vila et al., 2017; Zhang, 2017; Zhu and Li,
2017; Karami et al., 2018; Barretta et al., 2018b; Pisano et al., 2021).

Inconsistencies of limiting solutions of structural problems based on
Eringen’s two-phase theory will be definitely provided in the present
research. The plan is the following. The two-phase integral model for slender
beams is recalled in Sec.2 and the equivalent differential problem is also
provided. The purely nonlocal strain-driven law is then derived and discussed
in Sec.3. Exact structural solutions based on two-phase theory are given in
Sec.4 for simple applicative cases. Analytical limiting solutions, presented
in Sec.5, are shown to be in contrast with equilibrium and kinematic boundary
conditions. Consequently, such fields cannot be assumed as solutions of
structural problems formulated according to Eringen’s purely nonlocal strain-
driven model which is thus confirmed as inapplicable to nanocontinua of
technical interest.

##  2 Two-phase local/nonlocal theory for elastic beams

Let us consider a slender straight beam under flexure of length $\,L\,$. Beam
and bending axes are denoted by $\,x\,$ and $\,y\,$ respectively. According to
Bernoulli-Euler theory, the total curvature field $\,\chi:[0,L]\mapsto\Re\,$
associated with the transverse displacement field $\,v:[0,L]\mapsto\Re\,$ is
expressed by111 The symbol $\,\partial_{x}^{n}\,$ denotes $\,n$-times
differentiation along the beam axis $\,x\,$.

$\chi=\partial_{x}^{2}v=\chi^{el}+\chi^{nel}$ (1)

with $\,\chi^{el}\,$ elastic curvature and $\,\chi^{nel}\,$ all other non-
elastic curvature fields. Stress fields are described by bending moments
$\,M:[0,L]\mapsto\Re\,$ which have to fulfill the differential equation of
equilibrium

$\partial_{x}^{2}M=q$ (2)

with $\,q:[0,L]\mapsto\Re\,$ transversely distributed loading.

According to the two-phase elasticity model (Eringen, 1972), the bending
moment $\,M\,$ is convex combination of the source field $\,s:=K\chi^{el}\,$
and of the convolution between the source $\,s\,$ and a suitable averaging
kernel $\,\phi_{c}\,$

$M(x)=\alpha\>s(x)+(1-\alpha)\>\int_{0}^{L}{\phi_{c}(x,\xi)\>s(\xi)\>d\xi}$
(3)

where $\,K\,$ is the local elastic bending stiffness, i.e. second moment of
the field of Euler-Young elastic moduli $\,E\,$ on beam cross section.

The relation in Eq.(3) is a Fredholm integral equation of the second kind
(Tricomi, 1957; Polyanin and Manzhirov, 2008) in the unknown source field
$\,s\,$, with $\,0\leq\alpha\leq 1\,$ mixture parameter and $c>0$ nonlocal
scale parameter. The purely nonlocal law
$\displaystyle\,M(x)=\int_{0}^{L}{\phi_{c}(x,\xi)\>K\chi^{el}(\xi)\>d\xi}\,$
is got by setting $\,\alpha=0\,$, while for $\,\alpha=1\,$ the purely local
relation $\,M=K\chi^{el}\,$ is recovered.

The averaging kernel is assumed to be the bi-exponential function

$\phi_{c}(x)=\frac{1}{2c}\exp\biggl{(}-\frac{|x|}{c}\biggr{)}\,$ (4)

fulfilling symmetry, positivity and limit impulsivity (Eringen, 1983).

An equivalent differential formulation (Romano et al., 2017a) of the two-phase
model Eq.(3) can be got by observing that the special kernel in Eq.(4) is the
Green’s function of the linear differential operator $\,\mathcal{L}_{x}\,$
defined by

$\mathcal{L}_{x}:=1-c^{2}\,\partial_{x}^{2}$ (5)

Now, rewriting Eq.(3) as follows

$\big{(}M-\alpha\>K\chi^{el}\big{)}(x)=(1-\alpha)\>\int_{0}^{L}{\phi_{c}(x,\xi)\>K\chi^{el}(\xi)\>d\xi}$
(6)

and applying the differential operator $\,\mathcal{L}_{x}\,$ to Eq.(6), we get

$\mathcal{L}_{x}\big{(}M-\alpha\>K\chi^{el}\big{)}(x)=(1-\alpha)\>K\chi^{el}(x)$
(7)

which is the differential equation equivalent to the two-phase model in
Eq.(3).

Indeed, denoting by $\,\boldsymbol{\delta}\,$ the Dirac unit impulse, since
$\,\mathcal{L}_{x}\phi_{c}(x,\xi)\,=\boldsymbol{\delta}(x,\xi)\,$, we may
write

$\mathcal{L}_{x}\big{(}M-\alpha\>K\chi^{el}\big{)}(x)=(1-\alpha)\int_{0}^{L}{\mathcal{L}_{x}\phi_{c}(x,\xi)\>K\chi^{el}(\xi)\>d\xi}=(1-\alpha)K\chi^{el}(x)$
(8)

As proven in (Romano et al., 2017a), the special kernel in Eq.(4) satisfies
the following homogeneous boundary conditions

$\left\\{\begin{aligned} &\mathcal{B}_{0}\phi_{c}\,|_{0}=0\\\
&\mathcal{B}_{L}\phi_{c}\,|_{L}=0\end{aligned}\right.$ (9)

where $\,\mathcal{B}_{0}:=1-c\,\partial_{x}\,$ and
$\,\mathcal{B}_{L}:=1+c\,\partial_{x}\,$ are differential operators defined at
the boundary. By applying $\,\mathcal{B}_{0}\,$ and $\,\mathcal{B}_{L}\,$ to
Eq.(6) we get

$\left\\{\begin{aligned}
&\mathcal{B}_{0}\big{(}M-\alpha\>K\chi^{el}\big{)}\,|_{0}=0\\\
&\mathcal{B}_{L}\big{(}M-\alpha\>K\chi^{el}\big{)}\,|_{L}=0\end{aligned}\right.$
(10)

which are the constitutive boundary conditions associated with Eq.(7).

Finally, the equivalent differential problem in Eqs.(7)-(10) can be explicitly
formulated as follows

$\frac{M(x)}{c^{2}}-\,\partial_{x}^{2}M(x)=\frac{(K\chi^{el})(x)}{c^{2}}-\alpha\>\partial_{x}^{2}(K\chi^{el})(x)$
(11) $\left\\{\begin{aligned}
&\partial_{x}M(0)-\frac{1}{c}\,M(0)=\alpha\>\bigg{(}\partial_{x}(K\chi^{el})(0)-\frac{(K\chi^{el})(0)}{c}\bigg{)}\,\\\
&\partial_{x}M(L)+\frac{1}{c}\,M(L)=\alpha\>\bigg{(}\partial_{x}(K\chi^{el})(L)+\frac{(K\chi^{el})(L)}{c}\bigg{)}\end{aligned}\right.$
(12)

##  3 Limiting case: purely nonlocal elasticity

For vanishing mixture parameter $\,\alpha=0\,$, the purely nonlocal elasticity
model is recovered, that is

$M(x)=\int_{0}^{L}{\phi_{c}(x,\xi)\>(K\chi^{el})(\xi)\>d\xi}$ (13)

which is the Eringen’s purely nonlocal strain-driven law.

Remark. For any value of $\,c>0\,$, the integral convolution Eq.(13) admits a
unique solution (or no solution) if and only if the following constitutive
boundary conditions are satisfied by the bending field (Romano et al., 2017a)

$\left\\{\begin{aligned} &\partial_{x}M(0)-\frac{1}{c}\,M(0)=0\\\
&\partial_{x}M(L)+\frac{1}{c}\,M(L)=0\end{aligned}\right.$ (14)

Fulfillment of Eq.(14) ensures that differential equation (Eq.(15)) provides
the unique solution of the nonlocal problem (Romano et al., 2017a)

$\frac{M(x)}{c^{2}}-\,\partial_{x}^{2}M(x)=\frac{(K\chi^{el})(x)}{c^{2}}$ (15)

As proven in (Romano et al., 2017a), constitutive boundary conditions Eqs.(14)
are incompatible with equilibrium requirements, since they relate bending and
shearing fields through the characteristic nonlocal length $\,c\,$. The two-
phase model eliminates ill-posedness for $\,\alpha>0\,$. Indeed, left-hand
sides of Eq.(12) are known expressions (generally polynomial), unless of
$\,n\,$ integration constants, with $\,n\,$ redundancy degree. Hence, the laws
in Eq.(12) are non-homogeneous boundary conditions imposed on the elastic
curvature and its first derivative and no incompatibility arises with
kinematic boundary constraints since they act on displacement fields and its
first derivatives. However, as $\,\alpha\to 0\,$ the two-phase strain-driven
law will tend to an ill-posed theory and then inconsistencies of limiting
solutions are expected to occur. Predicted singularities for $\,\alpha\to 0\,$
will be investigated in detail in the sequel.

##  4 Case-studies and exact structural solutions

This section provides exact closed-form solutions of the elastic equilibrium
problem of inflected beams using Eringen’s two-phase model Eqs.(11)-(12).

The solution procedure is summarized as follows.

  * 1.

Step 1. Solution of the differential equilibrium condition in Eq.(2) equipped
with standard natural boundary conditions to get the bending moment as
functions of $n$ integration constants, with $n$ standing for redundancy
degree. For statically determinate beams, the equilibrated bending field is
univocally determined by equilibrium requirements.

  * 2.

Step 2. Inserting the obtained bending field $\,M\,$ in the equivalent
differential problem governed by Eqs.(11)-(12) to get the nonlocal elastic
curvature field $\,\chi^{el}\,$.

  * 3.

Step 3. Detection of nonlocal displacement field $v$ by solving Eq.(1) with
prescription of $\,3+n\,$ standard essential boundary conditions.

The solution procedure above is applied to solve the nonlocal elastostatic
problems below. The non-elastic curvature $\,\chi^{nel}\,$ is assumed to
vanish, so that, by virtue Eq.(1), total and elastic curvatures are coincident
$\,\chi=\chi^{el}\,$.

Cantilever under concentrated couple at free end.

Let us consider a cantilever of length $\,L\,$ and uniform bending stiffness
$\,K\,$ under a concentrated couple $\,\mathcal{M}\,$ applied at free end.

The bending moment field $\,M\,$ is obtained by the differential equilibrium
equation $\,\partial_{x}^{2}M=0\,$ equipped with standard natural boundary
conditions $\,M(L)=\mathcal{M}\,$ and $\,\partial_{x}M(L)=0\,$, so that
Eq.(11) becomes

$\frac{\mathcal{M}}{c^{2}}=\frac{K\chi^{el}(x)}{c^{2}}-\alpha\>K\partial_{x}^{2}\chi^{el}(x)$
(16)

supplemented with the constitutive boundary conditions

$\left\\{\begin{aligned}
&-\mathcal{M}\,=\alpha\>K\bigg{(}c\,\partial_{x}\chi^{el}(0)-\chi^{el}(0)\bigg{)}\,\\\
&\mathcal{M}\,=\alpha\>K\bigg{(}c\,\partial_{x}\chi^{el}(L)+\chi^{el}(L)\bigg{)}\end{aligned}\right.$
(17)

Solving the differential problem in Eqs.(16)-(17) provides the elastic
curvature $\,\chi^{el}\,$. Also, Eq.(1) with prescription of standard
essential boundary conditions $\,v(0)=0\,$ and $\,\partial_{x}v(0)=0\,$, gives
the nonlocal displacement field

$\begin{split}v(x)=&\,\frac{\mathcal{M}x^{2}}{2K}-\frac{\mathcal{M}}{\psi_{1}K}\,c\,(\alpha-1)\,e^{-\frac{x}{c\sqrt{\alpha}}}\big{(}(x-c\sqrt{\alpha})\,e^{\frac{L+x}{c\sqrt{\alpha}}}+\,c\,\sqrt{\alpha}\,e^{\frac{L}{c\sqrt{\alpha}}}\\\
&-e^{\frac{x}{c\sqrt{\alpha}}}(c\,\sqrt{\alpha}+x)+c\,\sqrt{\alpha}e^{\frac{2x}{c\sqrt{\alpha}}}\big{)}\end{split}$
(18)

with

$\begin{split}\psi_{1}=(\sqrt{\alpha}+1)\,e^{\frac{L}{c\sqrt{\alpha}}}+\sqrt{\alpha}-1\end{split}$
(19)

Cantilever under concentrated force at free end.

Let us consider a cantilever of length $\,L\,$ and uniform bending stiffness
$\,K\,$ under a concentrated force $\,F\,$ applied at free end.

The equilibrated bending moment field $\,M\,$ is defined by
$\,\partial_{x}^{2}M=0\,$ with prescription of standard natural boundary
conditions $\,M(L)=0\,$ and $\,\partial_{x}M(L)=-F\,$, so that Eq.(11) becomes

$\frac{F(L-x)}{c^{2}}=\frac{K\chi^{el}(x)}{c^{2}}-\alpha\>K\partial_{x}^{2}\chi^{el}(x)$
(20)

equipped with the constitutive boundary conditions

$\left\\{\begin{aligned} &-Fc-
FL\,=\alpha\>K\bigg{(}c\,\partial_{x}\chi^{el}(0)-\chi^{el}(0)\bigg{)}\,\\\
&-Fc\,=\alpha\>K\bigg{(}c\,\partial_{x}\chi^{el}(L)+\chi^{el}(L)\bigg{)}\end{aligned}\right.$
(21)

The unknown elastic curvature $\,\chi^{el}\,$ is obtained by solving the
differential problem in Eqs.(20)-(21). Then, Eq.(1) with prescription of
standard essential boundary conditions provides the following nonlocal
displacement field

$\begin{split}v(x)=\,&\frac{Fx^{2}\,(3L-x)}{6K}+\frac{(\alpha-1)\,c\,F}{\psi
K}\,e^{-\frac{x}{\sqrt{\alpha}c}}\,(\sqrt{\alpha}\,\psi_{1}\,c^{2}\,(e^{\frac{x}{\sqrt{\alpha}c}}-1)\,(e^{\frac{L}{\sqrt{\alpha}c}}+e^{\frac{x}{\sqrt{\alpha}c}})\\\
&-\psi_{2}\,L\,x\,e^{\frac{x}{\sqrt{\alpha}c}}+\,\psi_{3}\,c)\end{split}$ (22)

with

$\begin{split}&\psi=(\sqrt{\alpha}+1)^{2}\,e^{\frac{2L}{\sqrt{\alpha}c}}-(\sqrt{\alpha}-1)^{2}\\\
&\psi_{2}=\sqrt{\alpha}+(\sqrt{\alpha}+1)\,e^{\frac{2L}{\sqrt{\alpha}c}}-1\\\
&\psi_{3}=-(\sqrt{\alpha}+1)\,\sqrt{\alpha}\,L\,e^{\frac{2L}{\sqrt{\alpha}c}}+(\alpha-\sqrt{\alpha})\,L\,e^{\frac{2x}{\sqrt{\alpha}c}}-2\sqrt{\alpha}\,x\,e^{\frac{L+x}{\sqrt{\alpha}c}}+\\\
&(\sqrt{\alpha}+1)(\sqrt{\alpha}L-x)\,e^{\frac{2L+x}{\sqrt{\alpha}c}}-(\sqrt{\alpha}-1)\,e^{\frac{x}{\sqrt{\alpha}c}}(\sqrt{\alpha}L+x)\end{split}$
(23)

Simply supported beam under uniformly distributed loading.

Let us consider a simply supported beam of length $\,L\,$ and uniform bending
stiffness $\,K\,$ under uniformly distributed loading $\,q\,$.

Bending moment field $\,M\,$ is determined by differential equilibrium
equation $\,\partial_{x}^{2}M=q\,$ equipped with natural boundary conditions
$\,M(0)=0\,$ and $\,M(L)=0\,$, so that Eq.(11) becomes

$\frac{qx(x-L)}{2c^{2}}-\,q=\frac{K\chi^{el}}{c^{2}}-\alpha\>K\partial_{x}^{2}\chi^{el}$
(24)

supplemented with constitutive boundary conditions

$\left\\{\begin{aligned}
&-c\,\frac{qL}{2}\,=\alpha\,K\>\bigg{(}c\,\partial_{x}\chi^{el}(0)-\chi^{el}(0)\bigg{)}\,\\\
&c\,\frac{qL}{2}\,=\alpha\,K\>\bigg{(}c\,\partial_{x}\chi^{el}(L)+\chi^{el}(L)\bigg{)}\end{aligned}\right.$
(25)

The elastic curvature $\,\chi^{el}\,$ is obtained from Eqs.(24)-(25); then,
Eq.(1) together with essential boundary conditions $\,v(0)=0\,$ and
$\,v(L)=0\,$, provides the nonlocal displacement field

$\begin{split}v(x)=\,&\frac{qx\,(L^{3}-2Lx^{2}+x^{3})}{24K}-\frac{(\alpha-1)\,c^{2}\,q}{2\psi_{1}K}\,e^{-\frac{x}{\sqrt{\alpha}\,c}}\,(-e^{\frac{x}{\sqrt{\alpha}\,c}}\,(2\,\alpha^{3/2}\,c^{2}+L\,(\sqrt{\alpha}c\\\
&-\sqrt{\alpha}\,x+x)+(\sqrt{\alpha}-1)\,x^{2})+\,e^{\frac{L+x}{\sqrt{\alpha}c}}\,(-2\,\alpha^{3/2}\,c^{2}+L\,(-\sqrt{\alpha}\,c+\sqrt{\alpha}\,x\\\
&+x)-(\sqrt{\alpha}+1)\,x^{2})+\sqrt{\alpha}\,c\,e^{\frac{L}{\sqrt{\alpha}c}}\,(2\,\alpha\,c+L)+\sqrt{\alpha}\,c\,(2\,\alpha\,c+L)\,e^{\frac{2x}{\sqrt{\alpha}\,c}})\end{split}$
(26)

Doubly clamped beam under uniformly distributed loading.

Let us consider a doubly clamped beam of length $\,L\,$ and uniform bending
stiffness $\,K\,$ under uniformly distributed loading $\,q\,$.

The bending moment field $\,M\,$ obtained by the differential equilibrium
equation $\,\partial_{x}^{2}M=q\,$ is a function of two integration constants,
$\,a_{1}\,$ and $\,a_{2}\,$. Hence, Eq.(11) in terms of displacement field
$\,v\,$ becomes

$\frac{qx^{2}}{2}+a_{1}x+a_{2}-c^{2}q=K\partial_{x}^{2}v-\alpha\>c^{2}\,K\partial_{x}^{4}v$
(27)

supplemented with constitutive and kinematic boundary conditions

$\left\\{\begin{aligned} &c\,a_{1}-a_{2}=\alpha
K\,(c\,\partial_{x}^{3}v(0)-\partial_{x}^{2}v(0))\\\
&c\,(qL+a_{1})+\frac{qL^{2}}{2}+a_{1}L+a_{2}=\alpha
K\,(c\,\partial_{x}^{3}v(L)+\partial_{x}^{2}v(L))\\\ &v(0)=0\\\
&\partial_{x}v(0)=0\\\ &v(L)=0\\\ &\partial_{x}v(L)=0\end{aligned}\right.$
(28)

By solving the differential problem in Eqs.(27)-(28) the nonlocal displacement
field is given by

$\begin{split}v(x)=\,&\frac{q\,x^{2}\,(L-x)^{2}}{24K}+\frac{c\,(\alpha-1)\,q}{12\,\psi_{4}\,K}(12\,c^{2}+6c\,L+L^{2})\,e^{-\frac{x}{\sqrt{\alpha}c}}(\sqrt{\alpha}\,c\,L\,e^{\frac{L}{\sqrt{\alpha}c}}\\\
&-e^{\frac{L+x}{\sqrt{\alpha}c}}(\sqrt{\alpha}\,cL-
Lx+x^{2})+e^{\frac{x}{\sqrt{\alpha}c}}(-\sqrt{\alpha}\,cL-
Lx+x^{2})+\sqrt{\alpha}\,c\,L\,e^{\frac{2x}{\sqrt{\alpha}c}})\end{split}$ (29)

with

$\psi_{4}=\,L\,(-\sqrt{\alpha}-(\sqrt{\alpha}+1)\,e^{\frac{L}{\sqrt{\alpha}c}}+1)+2\,(\alpha-1)\,c\,(e^{\frac{L}{\sqrt{\alpha}c}}-1)$
(30)

##  5 Inconsistency of limiting solutions

As enlightened by theoretical outcomes in Sect.3, no solution can be found to
Eq.(3) for a vanishing mixture parameter, i.e.: $\,\alpha=0\,$. Indeed, for
any value $\,c>0\,$, constitutive boundary conditions in Eq.(12) can be
satisfied by equilibrated bending moments only for a strictly positive mixture
parameter. Hence, inconsistencies of two-phase solutions are expected to occur
as $\,\alpha\to 0\,$.

In the following, analytical limiting responses of case-studies in Sec.4 are
provided. Non-dimensional variables are used in parametric plots by
introducing the non-dimensional abscissa $\,\xi=x/L\,$, the nonlocal parameter
$\,\lambda=c/L\,$ and the following non-dimensional limiting displacement
fields

$\bar{v_{l}}=v_{l}\frac{K}{FL^{3}}\quad
or\quad\bar{v_{l}}=v_{l}\frac{K}{\mathcal{M}L^{2}}\quad
or\quad\bar{v_{l}}=v_{l}\frac{K}{qL^{4}}$ (31)

Cantilever under concentrated couple at free end.

Computing the limit as $\,\alpha\to 0\,$ of Eq.(18), we get the following
limiting displacement field

$v_{l}(x)=\frac{\mathcal{M}}{K}\bigg{(}\frac{x^{2}}{2}+c\,x\bigg{)}$ (32)

The displacement field in Eq.(32) is clearly in contrast with essential
boundary conditions prescribed by constraints. Indeed, except for the local
case $\,\lambda=0^{+}\,$, kinematic inconsistency is apparent in parametric
plots of displacement (Fig.1) and rotation fields (Fig.2) obtained from
derivation of Eq.(32).

Figure 1: Cantilever under concentrated couple at free end: displacement
$\,\bar{v_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 2: Cantilever under
concentrated couple at free end: rotation $\,\bar{\varphi_{l}}\,$ versus
$\,\xi\,$ for $\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$.

Double derivation of Eq.(32) leads to the limiting elastic bending curvature
$\,\chi_{l}\,$ that inserted into the integral convolution in Eq.(3) (with
$\,\alpha=0\,$) should provide the equilibrated bending moment. Moreover,
since the redundancy degree is zero, equilibrated bending interaction is the
uniform field univocally determined by equilibrium conditions. As clearly
shown by parametric plots in Fig.3, the limiting elastic nonlocal curvature
$\,\chi_{l}\,$ provides a bending moment which is not equilibrated; indeed, it
is not uniform and not equal to the applied couple, except for the asymptotic
local bending moment which is equilibrated for $\,\xi\in\,]0,1[\,$. A further
derivation leads to shear force field which is not vanishing (see Fig.4),
except for $\,\lambda=0^{+}\,$ with $\,\xi\in\,]0,1[\,$.

Figure 3: Cantilever under concentrated couple at free end: bending moment
$\,\bar{M_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 4: Cantilever under
concentrated couple at free end: shear force $\,\bar{T_{l}}\,$ versus
$\,\xi\,$ for $\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$.

Cantilever under concentrated force at free end.

Computing the limit as $\,\alpha\to 0\,$ of Eq.(22), we get the following
limiting displacement field

$v_{l}(x)=\frac{F}{K}\bigg{(}L\,\frac{x^{2}}{2}-\frac{x^{3}}{6}+c^{2}\,x+L\,c\,x\bigg{)}$
(33)

Kinematic inconsistencies are clearly shown in parametric plots of
displacement and rotation fields as a function of $\,\lambda\,$, except for
the local case $\,\lambda=0^{+}\,$ (see Figs.5-6). Double derivation of
Eq.(33) leads to the limiting bending curvature that inserted into the
integral convolution in Eq.(3), for $\,\alpha=0\,$, provides a bending moment
field which is not compatible with differential and boundary equilibrium
requirements.

Figure 5: Cantilever under concentrated force at free end: displacement
$\,\bar{v_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 6: Cantilever under
concentrated force at free end: rotation $\,\bar{\varphi_{l}}\,$ versus
$\,\xi\,$ for $\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$.

Indeed, since the redundancy degree is zero, the equilibrated bending moment
is the linear field univocally determined by equilibrium conditions. Instead,
as clearly shown by parametric plots in Fig.7, the limiting elastic bending
curvature $\,\chi_{l}\,$ provides a bending moment which is not equilibrated
since it is not linear and does not vanish at free-end of the cantilever,
except for the asymptotic local bending moment which is equilibrated for
$\,\xi\in\,]0,1]\,$. A further derivation leads to the shear force field which
is not equal to the value of the applied force (see Fig.8) and hence, it is
not uniform (except for $\,\lambda=0^{+}\,$ with $\,\xi\in\,]0,1[\,$), i.e.:
the emerging distributed loading is not vanishing.

Figure 7: Cantilever under concentrated force at free end: bending moment
$\,\bar{M_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 8: Cantilever under
concentrated force at free end: shear force $\,\bar{T_{l}}\,$ versus $\,\xi\,$
for $\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$.

Simply supported beam under uniformly distributed loading.

Computing the limit as $\,\alpha\to 0\,$ of Eq.(26), we get the following
limiting displacement

$v_{l}(x)=\frac{q}{K}\bigg{(}\frac{x^{4}}{24}-\frac{Lx^{3}}{12}-\frac{c^{2}x^{2}}{2}+\frac{L^{3}x}{24}(12\lambda^{2}+1)\bigg{)}$
(34)

Parametric plots of displacement and rotation fields as function of
$\,\lambda\,$ are shown in Figs.9-10. Double derivation of Eq.(34) leads to
the limiting bending curvature that inserted into the integral convolution in
Eq.(3) (for $\,\alpha=0\,$) provides the static fields shown in Figs.11-12.
The bending moment field is clearly incompatible with natural boundary
conditions and the shear force field is not linear, i.e.: the emerging
distributed loading is not uniform (see Fig.13). Hence, the interaction fields
for $\,\lambda>0\,$ are not compatible with differential and boundary
equilibrium requirements.

Figure 9: Simply supported beam under uniformly distributed loading:
displacement $\,\bar{v_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 10: Simply supported
beam under uniformly distributed loading: rotation $\,\bar{\varphi_{l}}\,$
versus $\,\xi\,$ for $\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure
11: Simply supported beam under uniformly distributed loading: bending moment
$\,\bar{M_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 12: Simply supported
beam under uniformly distributed loading: shear force $\,\bar{T_{l}}\,$ versus
$\,\xi\,$ for $\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 13:
Simply supported beam under uniformly distributed loading: emerging loading
$\,\bar{q_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$.

Doubly clamped beam under uniformly distributed loading.

The limiting displacement field as $\,\alpha\to 0\,$ of Eq.(29) is given by

$v_{l}(x)=\frac{q}{K}\left(\frac{x^{4}}{24}-\frac{Lx^{3}}{12}-\frac{c^{2}x^{2}}{2}+\frac{c\left(12\lambda^{2}+6\lambda+1\right)L^{2}x}{12(2\lambda+1)}+\frac{L^{2}x^{2}}{24(2\lambda+1)}\right)$
(35)

The displacement field in Eq.(35) is in contrast with essential boundary
conditions prescribed by constraints. These kinematic inconsistencies are
clearly shown in parametric plots of displacement and rotation fields (see
Figs.14-15). Double derivation of Eq.(33) leads to the limiting nonlocal
elastic curvature that is used to compute the integral convolution in Eq.(3),
for $\,\alpha=0\,$, to get the bending moment field depicted in Fig.16. A
further derivation leads to shear force field which is not linear (see Fig.17)
and the emerging distributed loading is not uniform and not equal to the
applied one (see Fig.18). Hence, limiting solutions of the elastostatic
problem do not satisfy kinematic compatibility and equilibrium requirements,
except for the asymptotic local displacement and rotation fields
$\,\forall\xi\in\,[0,1]\,$ and for the asymptotic local bending and shearing
fields for $\,\xi\in\,]0,1[\,$.

Figure 14: Doubly clamped beam under uniformly distributed loading:
displacement $\,\bar{v_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 15: Doubly clamped
beam under uniformly distributed loading: rotation $\,\bar{\varphi_{l}}\,$
versus $\,\xi\,$ for $\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure
16: Doubly clamped beam under uniformly distributed loading: bending moment
$\,\bar{M_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 17: Doubly clamped
beam under uniformly distributed loading: shear force $\,\bar{T_{l}}\,$ versus
$\,\xi\,$ for $\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$. Figure 18:
Doubly clamped beam under uniformly distributed loading: emerging loading
$\,\bar{q_{l}}\,$ versus $\,\xi\,$ for
$\,\lambda\in\\{0^{+},0.1,0.2,0.3,0.4,0.5\\}\,$.

##  6 Closing remarks

The bending behaviour of elastic beams has been investigated by Eringen’s two-
phase integral theory. Exact closed-form solutions of simple structural
problems of applicative interest have been preliminarily provided in terms of
mixture and nonlocal parameters. Then, relevant asymptotic fields (as the
mixture parameter tends to zero) have been evaluated and studied. It has been
proven that such fields are in contrast with kinematic boundary conditions and
equilibrium requirements for any nonlocal parameter value. Accordingly,
limiting responses of Eringen’s two-phase formulation cannot be used as
solutions of the ill-posed Eringen purely nonlocal structural problem.
Therefore, inaccurate outcomes contributed in current literature (see e.g.
Mikhasev and Nobili, 2020; Mikhasev, 2021), regarding the well-posedness of
limiting elastodynamic problems based on Eringen’s two-phase continuum theory,
should be amended consequently.

Acknowledgment \- Financial support from the MIUR in the framework of the
Project PRIN 2017 - code 2017J4EAYB _Multiscale Innovative Materials and
Structures (MIMS)_ ; University of Naples Federico II Research Unit - is
gratefully acknowledged.

## References

  * Banejad et al. (2020) Banejad, A., Passandideh-Fard, M., Niknam, H., Mirshojaeian Hosseini, M.J., Mousavi Shaegh, S.A., 2020. Design, fabrication and experimental characterization of whole-thermoplastic microvalves and micropumps having micromilled liquid channels of rectangular and half-elliptical cross-sections. Sensors and Actuators A 301, 111713.
  * Barretta et al. (2018a) Barretta, R., Fabbrocino, F., Luciano, R., Marotti de Sciarra, F., 2018a. Closed-form solutions in stress-driven two-phase integral elasticity for bending of functionally graded nano-beams. Physica E: Low-dimensional Systems and Nanostructures 97, 13–30.
  * Barretta et al. (2018b) Barretta, R., Čanadija, M., Feo, L., Luciano, R., Marotti de Sciarra, F., Penna, R., 2018b. Exact solutions of inflected functionally graded nano-beams in integral elasticity. Composites Part B 142, 273-286.
  * Barretta et al. (2018c) Barretta, R., Diaco, M., Feo, L., Luciano, R., Marotti de Sciarra, F., Penna, R., 2018c. Stress-driven integral elastic theory for torsion of nano-beams. Mechanics Research Communications 87, 35-41.
  * Barretta et al. (2019) Barretta, R., Marotti de Sciarra, F., Vaccaro, M.S., 2019. On nonlocal mechanics of curved elastic beams. International Journal of Engineering Science 144, 103140.
  * Bažant and Jirásek (2002) Bažant, Z.P., Jirásek, M., 2002. Nonlocal integral formulation of plasticity and damage: survey of progress. Journal of Engineering Mechanics - ASCE, 128, 1119-1149.
  * Borino et al. (2003) Borino, G., Failla, B., Parrinello, F., 2003. A symmetric nonlocal damage theory. International Journal of Solids and Structures 40, 3621-3645.
  * Chao et al. (2020) Chao, M., Wang, Y., Ma, D., Wu, X., Zhang, W., Zhang, L., Wan, P., 2020. Wearable mxene nanocomposites-based strain sensor with tile-like stacked hierarchical microstructure for broad-range ultrasensitive sensing. Nano Energy 78, 105187.
  * Chorsi and Chorsi (2018) Chorsi, M.T. , Chorsi, H.T., 2018. Modeling and analysis of MEMS disk resonators. Microsystem Technologies 24 (6), 2517–2528 .
  * Eringen (1972) Eringen, A. C., 1972. Linear theory of nonlocal elasticity and dispersion of plane waves. International Journal of Engineering Science 10(5), 425-435.
  * Eringen (1983) Eringen, A. C., 1983. On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves. Journal of Applied Physics 54, 4703.
  * Eringen (1987) Eringen, A. C., 1987. Theory of nonlocal elasticity and some applications. Res Mechanica 21, 313-342.
  * Eptaimeros et al. (2016) Eptaimeros, K.G., Koutsoumaris, C.C., Tsamasphyros, G.J., 2016. Nonlocal integral approach to the dynamical response of nano-beams. Int J Mech Sci 115–116, 68–80 .
  * Fathi and Ghassemi (2017) Fathi, M., Ghassemi, A. (2017). The effects of surface stress and nonlocal small scale on the uniaxial and biaxial buckling of the rectangular piezoelectric nanoplate based on the two variable-refined plate theory. Journal of the Brazilian Society of Mechanical Sciences and Engineering 39 3203-3216.
  * Fernández-Sáez and Zaera (2017) Fernández-Sáez, J., Zaera, R., 2017. Vibrations of Bernoulli-Euler beams using the two-phase nonlocal elasticity theory. International Journal of Engineering Science 119, 232-248.
  * Ghayesh and Farajpour (2019) Ghayesh, M.H., Farajpour, A., 2019. A review on the mechanics of functionally graded nanoscale and microscale structures. International Journal of Engineering Science 137, 8-36.
  * Ghayesh and Farokhi (2020) Ghayesh, M.H., Farokhi, H., 2020. Nonlinear broadband performance of energy harvesters. Int J Eng Sci 147, 103202
  * Jankowski et al. (2020) Jankowski, P., Żur, K.K., Kim, J., Reddy, J.N., 2020. On the bifurcation buckling and vibration of porous nanobeams. Composite Structures 250, 112632.
  * Karami et al. (2018) Karami, B., Shahsavari, D., Janghorban, M., Li, L., 2018. Wave dispersion of mounted graphene with initial stress. Thin-Walled Structures 122, 102-111.
  * Kiani (2021) Kiani, K., Żur, K.K., 2021. Vibrations of double-nanorod-systems with defects using nonlocal-integralsurface energy-based formulations. Composite Structures 256, 113028.
  * Khodabakhshi and Reddy (2015) Khodabakhshi, P., Reddy, J.N., 2015. A unified integro-differential nonlocal model. International Journal of Engineering Science 95, 60-75.
  * Kröner et al. (1967) Kröner, E., 1967. Elasticity theory of materials with long range cohesive forces. International Journal of Solids and Structures 3(5), 731-742.
  * Lu et al. (2019) Lu, P., Veletić, M., Laasmaa, M., Vendelin, M., Louch, W.E., Steinar Halvorsen P., Bergsland, J., Balasingham, I., 2019. Multi-nodal nano-actuator pacemaker for energy-efficient stimulation of cardiomyocytes. Nano Communication Networks 22, 100270.
  * Malikan et al. (2020b) Malikan, M., Eremeyev, V.A., Żur, K.K., 2020. Effect of axial porosities on flexomagnetic response of in-plane compressed piezomagnetic nanobeams. Symmetry 12(12), 1935.
  * Mikhasev and Nobili (2020) Mikhasev, G., Nobili, A., 2020. On the solution of the purely nonlocal theory of beam elasticity as a limiting case of the two-phase theory. International Journal of Solids and Structures 190, 47–57.
  * Mikhasev (2021) Mikhasev, G., 2021. Free high-frequency vibrations of nonlocally elastic beam with varying cross-section area. Continuum Mechanics and Thermodynamics. https://doi.org/10.1007/s00161-021-00977-6
  * Oskouie et al. (2018) Oskouie, M. F., Ansari, R., Rouhi, H., 2018. Bending of Euler-Bernoulli nanobeams based on the strain-driven and stress-driven nonlocal integral models: a numerical approach. Acta Mechanica Sinica. Doi: 10.1007/s10409-018-0757-0
  * Pinnola (2020) Pinnola, F.P., Vaccaro, M.S., Barretta, R., Marotti de Sciarra, F., 2020. Random vibrations of stress-driven nonlocal beams with external damping. Meccanica URL https: //doi.org/10.1007/s11012-020-01181-7.
  * Pisano and Fuschi (2003) Pisano, A.A., Fuschi, P., 2003. Closed form solution for a nonlocal elastic bar in tension. International Journal of Solids and Structures 40, 13-23.
  * Pisano et al. (2021) Pisano, A.A., Fuschi, P., Polizzotto, C., 2021. Integral and differential approaches to Eringen’s nonlocal elasticity models accounting for boundary effects with applications to beams in bending. ZAMM - Journal of Applied Mathematics and Mechanics. https://doi.org/10.1002/zamm.202000152
  * Polizzotto (2001) Polizzotto, C., 2001. Nonlocal elasticity and related variational principles. International Journal of Solids and Structures 38, 7359-7380.
  * Polizzotto (2002) Polizzotto, C., 2002. Thermodynamics and continuum fracture mechanics for nonlocal-elastic plastic materials. European Journal of Mechanics A/Solids 21, 85-103.
  * Polyanin and Manzhirov (2008) Polyanin, A.D., Manzhirov A.V., 2008. Handbook of integral equations. 2nd ed. Boca Raton, FL: Chapman & Hall/CRC.
  * Roghani (2020) Roghani, M., Rouhi, H., 2020. Nonlinear stress-driven nonlocal formulation of timoshenko beams made of fgms. Continuum Mechanics and Thermodynamics. https://doi.org/10. 1007/s00161-020-00906-z
  * Rogula (1965) Rogula, D., 1965. Influence of spatial acoustic dispersion on dynamical properties of dislocations. Bulletin de l’Académie Polonaise des Sciences, Séries des Sciences Techniques 13, 337-343 .
  * Rogula (1982) Rogula, D., 1982. Introduction to nonlocal theory of material media. In D. Rogula (Ed.), Nonlocal theory of material media, CISM courses and lectures (vol. 268, pp. 125-222). Wien: Springer.
  * Romano et al. (2017a) Romano, G., Barretta, R., Diaco, M., Marotti de Sciarra, F., 2017a. Constitutive boundary conditions and paradoxes in nonlocal elastic nano-beams. International Journal of Mechanical Sciences 121, 151-156.
  * Romano et al. (2017b) Romano, G., Barretta, R., Diaco, M., 2017b. On nonlocal integral models for elastic nano-beams. International Journal of Mechanical Sciences 131-132, 490-499.
  * Sedighi et al. (2017) Sedighi, H.M., Malikan, M., Valipour, A., Żur, K.K., 2020. Nonlocal vibration of carbon/boron-nitride nano-hetero-structure in thermal and magnetic fields by means of nonlinear finite element method. Journal of Computational Design and Engineering 7(5), 591–602.
  * Soukarié et al. (2020) Soukarié, D., Ecochard, V., Salomé, L., 2020. Dna-based nanobiosensors for monitoring of water quality. Int. J. Hyg. Environ. Health 226, 113485.
  * Tricomi (1957) Tricomi, F.G., 1957. Integral Equations. Interscience, New-York, USA. Reprinted by Dover Books on Mathematics, 1985.
  * Vila et al. (2017) Vila, J., Fernández-Sáez, J., Zaera, R., 2017. Nonlinear continuum models for the dynamic behavior of 1D microstructured solids. International Journal of Solids and Structures 117, 111-122.
  * Wang et al. (2016) Wang, Y., Zhu, X., Dai, H., 2016. Exact solutions for the static bending of Euler-Bernoulli beams using Eringen two-phase local/nonlocal model. AIP Advances 6(8), 085114. Doi:10.1063/1.4961695
  * Zhang (2017) Zhang, Y., 2017. Frequency spectra of nonlocal Timoshenko beams and an effective method of determining nonlocal effect. International Journal of Mechanical Sciences 128-129, 572-582.
  * Zhang (2020) Zhang, J.Q., Qing, H., Gao, C.F., 2020. Exact and asymptotic bending analysis of microbeams under different boundary conditions using stress-derived nonlocal integral model. Z. Angew. Math. Mech. 100(1), e201900148.
  * Zhu and Li (2017) Zhu, X., Li, L., 2017. Closed form solution for a nonlocal strain gradient rod in tension. International Journal of Engineering Science 119, 16-28.
  * Żur et. al (2020) Żur, K.K., Arefi, M., Kimc, J., Reddy, J.N., 2020. Free vibration and buckling analyses of magneto-electro-elastic FGM nanoplates based on nonlocal modified higher-order sinusoidal shear deformation theory. Composites Part B 182, 107601.


# AD-DMKDE: Anomaly Detection through Density Matrices and Fourier Features
††thanks: Citation: Joseph et al., AD-DMKDE: Anomaly Detection through Density
Matrices and Fourier Features.

Oscar A. Bustos-Brinez , Joseph A. Gallego-Mejia , Fabio A. González
MindLab
Universidad Nacional de Colombia
Bogotá, Colombia
<EMAIL_ADDRESS>0000-0003-0704-9117
0000-0001-8971-4998 0000-0001-9009-7288

###### Abstract

This paper presents a novel density estimation method for anomaly detection
using density matrices (a powerful mathematical formalism from quantum
mechanics) and Fourier features. The method can be seen as an efficient
approximation of Kernel Density Estimation (KDE). A systematic comparison of
the proposed method with eleven state-of-the-art anomaly detection methods on
various data sets is presented, showing competitive performance on different
benchmark data sets. The method is trained efficiently and it uses
optimization to find the parameters of data embedding. The prediction phase
complexity of the proposed algorithm is constant relative to the training data
size, and it performs well in data sets with different anomaly rates. Its
architecture allows vectorization and can be implemented on GPU/TPU hardware.

_K_ eywords density matrix $\cdot$ random features $\cdot$ anomaly detection
$\cdot$ quantum machine learning.

##  1 Introduction

An anomaly can be broadly defined as an observation or datum that deviates
significantly from the patterns of the data set from which it originates, in
one or more features. In most cases, data are generated by complex processes
and allow for different types of measurements, so an anomaly may contain
valuable information about anomalous behaviors of the generative processes, or
elements that are impacting the generation or measurement stages [1]. Then,
recognizing this type of data (which can be referred to as unusual, atypical,
unexpected, malicious or rare, depending on the scenario), discerning real,
meaningful anomalies from normal noisy samples (known as ‘outliers’) and
identifying the unusual processes that originate them are the main objectives
of Anomaly Detection (AD) [2]. Methods and algorithms that perform AD are key
in various applications such as bank fraud detection, identification of abrupt
changes in sensors, medical diagnostics, natural sciences, network security
intrusion detection, among many others [21].

However, classic AD methods face significant challenges that limit their
performance and range of application [14]. In most cases, anomalies are much
more sparse and scarce than normal data, and are not identified a priori as
anomalies, which makes the use of supervised classifiers difficult.
Furthermore, the boundaries between “normal" and “anomaly" regions are very
subjective and depend on the particularities of each case. Data may also be
contaminated by noise, so that both types of data may be intermingled,
blurring the very idea of separation. Many of the AD algorithms address some
of these difficulties, but at the cost of being vulnerable to others; however,
the combination of classical methods with a good mathematical foundation and
deep learning algorithms is one of the most promising paths [11], although it
also presents shortcomings on the detection of complex types of anomalies and
fewer possibilities of explanation.

The core idea of the proposed method is to use random Fourier features to
approximate a Gaussian kernel centered in each training sample and then use a
density matrix to summarize these kernels in an efficient and compact way. The
density matrix is then used to estimate the density of new samples, and those
whose density lies below a certain threshold are classified as anomalies. The
method uses optimization to obtain the best random Fourier feature embedding
(a process called “Adaptive Fourier Features"), and is able to calculate the
best threshold using percentiles from a validation data set.

The outline of the paper is as follows: in Section 2, we describe anomaly
detection and density estimation in more depth, and present the benchmark
methods to which we will compare our algorithm. In Section 3, we present the
novel method, explaining the stages of the algorithm and how it uses Fourier
features and density matrices. In Section 4, we systematically compare the
proposed algorithm with state-of-the-art anomaly detection algorithms. In
Section 5, we state the conclusions of this work and sketch future research
directions.

##  2 Background and Related Work

###  2.1 Anomaly Detection

The main mechanism of anomaly detection algorithms is the construction of a
model that determines a degree of “normality" for the data points, and then
detects anomalies as points that deviate from this model. The importance of
anomaly detection lies in the fact that recognizing and understanding
anomalous behavior allows one to make decisions and predictions. For example,
anomalous traffic patterns on a network could mean that sensitive data is
being transmitted across it, or that there is an intruder attempting to attack
a system [3]. Anomalous user behavior in payment transactions could signify
credit card fraud and give some clues as to how to avoid it in the future
[23]. Anomaly detection should be distinguished from outlier detection, the
purpose of which is data cleaning: to locate outliers (mainly normal data
affected by noise) and subsequently remove them from the data set.

Let $r$ be the ratio of anomaly samples with respect to all samples. When $r$
is high, the most common approach to anomaly detection is to use a supervised
classification method. However, when $r$ is low (which is usually the case),
the best approach is to use semi-supervised or unsupervised AD algorithms.
According to [16], classical AD algorithms can be classified into six main
types: classification-based approaches (such as the well-known one-class SVM),
probabilistic models (such as histograms and density estimation methods like
KDE [8]), clustering models, information-based methods, spectral analysis
(where dimensionality reduction methods such as PCA are found), and distance-
based methods (including isolation forest and nearest neighbors).

In recent years, these classical models have been progressively replaced by
deep learning-based algorithms, showing high performance in anomaly detection
tasks [25]. These algorithms can efficiently model the inherent features of
complex data through their distributed and hierarchical architectures, i.e.,
they can perform implicit feature engineering, especially when a large amount
of data is available [20]. For this reason, they are commonly used for anomaly
detection in scenarios involving very large data sets with high-dimensional
instances, such as speech recognition, image/video analysis, or natural
language processing [14].

###  2.2 Anomaly Detection Baseline Methods

In this paper, we select 11 state-of-the-art anomaly detection methods that
represent a complete sample of the most common types of anomaly detection
methods. All algorithms consider the proportion of anomalies in the data as a
necessary parameter for finding threshold values; other parameters for
specific algorithms are mentioned in later sections.

We selected five well-known methods that are based on classic mathematical
approaches to anomaly detection. These methods include One Class SVM [24], a
kernel-based algorithm that encloses normal data in a boundary that leaves
anomalies outside of it; Covariance Estimator, that finds the smallest
ellipsoid that wraps normal data; Isolation Forest [12], that tries to
separate points using decision trees, and those who are easier to isolate are
the outliers; Local Outlier Factor (LOF) [4], based on a distance measure from
each point to a set of its neighbors; and K-nearest neighbor (KNN) [18], that
makes an scoring based on the distance from the $k$-th nearest neighbor.

We also decided to include methods developed in the last decade that do not
use neural networks into their architectures, but instead take other
approaches. These type of methods include SOS [7], that builds an affinity
matrix between all points that acts as a similarity measure; COPOD [10], that
builds univariant probability functions fro each dimension and then joins them
in a unified multivariant function that models the true distribution of the
data; and LODA [15], that combines classifiers in low dimensions and uses
histograms to detect anomalies.

Finally, we complete the overview of baseline methods with some of the
proposals that are based on the use of neural networks as their central
elements. These three models include VAE-Bayes [9], built around a variational
autoencoder that builds a latent space from data, where the probability
distribution of the data can be retrieved by using Bayesian assumptions; DSVDD
[22], in where a neural network is used to transform data into a latent space
where normal data can be encompasssed into a hypersphere whose features are
neurally optimized; and LAKE [13], that includes a variational autoencoder to
reduce dimensionality in a way that preserves data distribution, and then
performs KDE in this space and separates anomalies using a threshold.

##  3 Anomaly Detection through Density Matrices and Fourier Features (AD-
DMKDE)

González et al. [5] proposed a new algorithm based on the joining of density
matrices and random Fourier features for neural density estimation, called
“Density matrix for Kernel Density Estimation" (DMKDE), that has its core in
kernel density estimation. The DMKDE algorithm consists of a random Fourier
feature mapping, a density matrix calculation, and a quantum density
estimation phases. This method, whose original purpose was only density
estimation, is used as the foundation for the novel method presented here.

AD-DMKDE enhances the random Fourier feature mapping from DMKDE through the
use of optimization, that allows to find better values for the parameters of
the mapping function. This process, called Adaptive Fourier features, can be
seen in more detail in subsection 3.1. After that, the density matrix is built
from the mappings of the training samples (see subsection 3.2), and the
density of new samples is calculated through a process known as Quantum
Density Estimation, inspired in quantum measurement theory from quantum
mechanics (see subsection 3.3). Then, as final step, a threshold is obtained
through a cross-validation process over a validation data set. This threshold
is calculated using the rate of the anomalies in the validation set and is
explained in the Subsection 3.4. The test phase of the algorithm comprises
applying the already-built Fourier feature mapping over the test set, Quantum
Density Estimation and threshold comparison to make the decision whether each
test point is anomalous or not. Figure 1 shows a summary of the method.

Figure 1: AD-DMKDE algorithm. Each section corresponds to a different stage of
the method. The parameters of the mapping/embedding and the density matrix do
not change after they are calculated.

###  3.1 Random Fourier features and Adaptive Fourier features

In the training phase, each data point in the training data set is implicitly
transformed into a higher dimension. The Random Fourier feature mapping, as
first stated in [17], applies the following function
$\text{rff}(\bf{x}_{\textit{i}})=\sqrt{2}\cos(\bf{w}^{\textit{T}}\bf{x}_{\textit{i}}+\bm{b})$
where $\bf{w}$ is a vector of $iid$ from a $\mathcal{N}(\bm{0},\bm{I}_{D})$
and $\bf{b}$ is a $iid$ sample from a $\text{Uniform}(0,2\pi)$. According to
Bochner’s theorem, the kernel $k$ can be approximated as an expectation,
$k(\bm{x},\bm{y})\simeq\mathbb{E}_{\bm{w}}[\langle\phi_{\text{rff}}(\bm{y}),\phi_{\text{rff}}(\bm{x})\rangle]$.
So, since we expect to approximate a Gaussian kernel to perform KDE, the
parameters of the kernel, in this case the $\sigma$ (standard deviation)
parameter, is one of the major influences over the selection of $\bf{w}$ and
$\bf{b}$.

This mapping is improved through AFF (Adaptive Fourier features), a mechanism
that, before of the RFF mapping, performs an optimization of the values of
$\bf{w}$ and $\bf{b}$ by building a neural network that adjust their values
through a learning process. After the application of the RFF formula with the
optimized values of $\bf{w}$ and $\bf{b}$, the obtained results are normalized
to obtain the final AFF embedding $\phi_{\text{aff}}(\bf{x}_{\textit{i}})$
that will be the representation of each $\bf{x}_{\textit{i}}$.

###  3.2 Density matrix calculation

Once the AFF embedding has been applied to all training samples, the
transformed training data correspond to the set
$\\{\phi_{\text{aff}}(\bf{x}_{\textit{1}}),\phi_{\text{aff}}(\bf{x}_{\textit{2}}),...,\phi_{\text{aff}}(\bf{x}_{\textit{n}})\\}$,
being $n$ the number of train samples. For each sample, the algorithm computes
its corresponding pure-state density matrix $\textsl{R}_{i}$, given by the
outer product of each $\phi_{\text{aff}}(\bf{x}_{\textit{i}})$:

$\textsl{R}_{i}=\left[\phi_{\text{aff}}(\bf{x}_{\textit{i}})\right]^{T}\phi_{\text{aff}}(\bf{x}_{\textit{i}})$

The mixed-state density matrix that comprehends all the training samples is
calculated as $\textsl{R}=1/n\sum_{i=1}^{n}\textsl{R}_{i}$. Although R
contains all the information of the embedded train samples, its size does not
depend on $n$; instead, it only depends on the size of the embedding. This is
a major improvement over other methods, making it possible to compactly
represent (and store) very large data sets using only a single matrix whose
size is much smaller in general than the training set. Other advantage of the
use of the density matrix representation is the fact that the calculation of R
only has to occur once in the entire training phase.

###  3.3 Quantum Density Estimation

With the calculation of R, we can estimate the density for new samples, in
particular those of the validation and test sets. Given a validation sample
$\bf{x_{\textit{j}}}$, we need to transform it by using the exact same mapping
and normalization we previously defined in order to obtain
$\phi_{\text{aff}}(\bf{x_{\textit{j}}})$, and then calculate

$\hat{f}(\bf{x}_{\textit{j}})=[\phi_{\text{aff}}(\bf{x}_{\textit{j}})]^{\textit{T}}\,\textsl{R}\>[\phi_{\text{aff}}(\bf{x}_{\textit{j}})]$

as an estimate of the density of $\bf{x_{j}}$. This process is repeated for
every validation sample, in order to build the validation density set
$\hat{F}_{val}=\left\\{\hat{f}(\bf{x}_{\textit{j}})\>|\>\bf{x}_{\textit{j}}\in\mathit{X_{val}}\right\\}$.

###  3.4 Threshold calculation and prediction phase

To determine which of the samples are anomalies or not, we need to find a
threshold $\theta$ that can be used as a sort of boundary between these two
types of data points. Therefore, we use the anomaly rate of the data set
(either an a priori known value or the proportion of anomalies we expect to
find). We take the estimates from the validation set and find the percentile
that corresponds to the anomaly rate we defined above to calculate the
threshold $\theta:=q_{(\text{anomaly rate \%})}(\hat{F}_{\mathit{val}})$.

As a final step, we estimate the density of the test samples (by mapping them
and applying density estimation with R), and then use the threshold $\theta$
to discriminate them. So, the classification of a test sample $\bf{x_{k}}$ is
given by:

$\hat{y}(\bf{x}_{\textit{k}})=\left\\{\begin{matrix}\text{
`normal'}&\textit{if }\hat{f}(\bf{x}_{\textit{k}})\geq\theta\\\ \text{
`anomaly'}&\textit{otherwise}\end{matrix}\right.$

##  4 Experimental Evaluation

###  4.1 Experimental Setup

For our experiments, we compared our proposed method with all the baseline
algorithms listed above (see section 2.2). To run the first four methods, we
used the Python implementation provided by the scikit-learn library. KNN, the
shallow methods, VAE and DVSDD were executed through the implementation
provided by the PyOD Python library [26]. For the LAKE algorithm [13], the
implementation we used comes from the Github repository of its
authors111https://github.com/1246170471/LAKE, although we had to correct some
issues in the original code, particularly the way to split the test data set
to include both normal and anomalous samples.

To handle the inherent randomness that AD-DMKDE can present in some stages
(particularly in AFF neural network training), we decided to select a unique,
invariant value for every random seed that affect the behavior of the method.
All the experiments were carried out on a machine with a 2.1GHz Intel Xeon
64-Core processor with 128GB RAM and two NVIDIA RTX A5000 graphic processing
units, that run Ubuntu 20.04.2 operating system.

####  4.1.1 Data sets.

We used twenty public data sets to evaluate the performance of our proposed
algorithm in performing anomaly detection. We chose data sets with a wide
variety of characteristics, such as their size or anomaly rate, to test the
robustness of our algorithm in multiple scenarios. The main characteristics of
all the data sets can be seen in Table 1. The source for all the data sets was
the ODDS Library of Stony Brook University [19].

Table 1: Main features of the data sets. Data Set | Instances | Dimensions | Outlier Rate
---|---|---|---
Arrhythmia | 452 | 274 | 0,146
Cardio | 2060 | 22 | 0,2
SpamBase | 3485 | 58 | 0,2
Thyroid | 3772 | 36 | 0,0247
KDDCUP | 10000 | 118 | 0,1934
Glass | 214 | 9 | 0,042
Lympho | 148 | 18 | 0,04
Ionosphere | 351 | 33 | 0,359
Letter | 1600 | 32 | 0,0625
MNIST | 7603 | 100 | 0,092
Musk | 3062 | 166 | 0,0317
OptDigits | 5216 | 64 | 0,0288
PenDigits | 6870 | 16 | 0,0227
Pima | 768 | 8 | 0,349
Satellite | 6435 | 36 | 0,3164
SatImage | 5803 | 36 | 0,0122
Shuttle | 10000 | 9 | 0,0715
Vertebral | 240 | 6 | 0,125
Vowels | 1456 | 12 | 0,03434
WBC | 378 | 30 | 0,0556

####  4.1.2 Experiment Configuration.

For each data set, the following configuration was chosen: the data set was
split in a stratified way (keeping the same proportion of outliers in each
subset) by randomly taking 30% of the samples as the test set, and from the
remaining samples, again randomly taking 30% for validation and the remaining
70% as the training set. This splitting was performed only once per data set,
so that all algorithms worked with exactly the same data partitions.

When implementing the Scikit-Learn and PyOD algorithms, the default
configurations defined in these libraries were used, modifying only the
parameters referred to the data set features mentioned in the previous
section. The only exception to this rule was the VAE algorithm, whose
architecture, based on autoencoders, needed to be slightly modified when using
low-dimensional data sets, to ensure that the encoded representation was
always of lower dimensionality than the original sample. Also, as mentioned
above, LAKE algorithm required some corrections to the original code published
by its authors; however, the basic internal structure, the internal logic of
the algorithm and the functions on which it is based remained the same as in
its original configuration.

Our AD-DMKDE algorithm was implemented from the original code of DMKDE,
presented in [6], with some modifications. Since the main goal of DMKDE is
density estimation, the separation of data as normal or outlier was performed
by adding a function that calculates the classification threshold using the
rate of anomalies in each data set, a known value in all cases (see table 1).

####  4.1.3 Evaluation Metrics.

The main metric we chose to determine the performance of the algorithms was
the F1 score (with weighted average), a well-known and very common metric used
when testing classification algorithms. However, the calculation of other
common metrics, such as accuracy and the exclusive F1-Score to the anomalous
class, was also implemented. For each algorithm, a set of parameters of
interest was selected in order to perform a comprehensive search for the
combinations of values that gave the best result for each data set. This
search was performed by training only with the training set and reporting the
results of the validation set. The selection of the best combination of
parameters was performed by choosing the model with the best average F1 score.
Once the best parameters were found, these same values were used to run each
algorithm once again, this time combining the training and validation sets for
the training phase and reporting the results on the test set. These final
results were the values we used to compare the overall performance of all
algorithms.

###  4.2 Results and discussion

Table 2: F1 Score for all classifiers over all data sets. The first and second best values are marked in bold and underlined, respectively. Data Set | OCSVM | iForest | Cov. | LOF | KNN | SOS | COPOD | LODA | VAE-B | DSVDD | LAKE | AD-DMKDE
---|---|---|---|---|---|---|---|---|---|---|---|---
Arrhythmia | 0.813 | 0.821 | 0.818 | 0.804 | 0.861 | 0.773 | 0.844 | 0.798 | 0.856 | 0.864 | 0.909 | 0.911
Cardio | 0.804 | 0.752 | 0.756 | 0.702 | 0.753 | 0.739 | 0.750 | 0.717 | 0.783 | 0.735 | 0.772 | 0.831
Glass | 0.916 | 0.931 | 0.931 | 0.925 | 0.900 | 0.907 | 0.916 | 0.848 | 0.900 | 0.907 | 1.000 | 0.974
Ionosphere | 0.765 | 0.710 | 0.876 | 0.830 | 0.817 | 0.784 | 0.736 | 0.510 | 0.714 | 0.674 | 0.993 | 0.959
KDDCUP | 0.744 | 0.838 | 0.975 | 0.789 | 0.932 | 0.716 | 0.785 | 0.770 | 0.803 | 0.755 | 1.000 | 0.984
Letter | 0.893 | 0.897 | 0.909 | 0.930 | 0.910 | 0.911 | 0.895 | 0.899 | 0.897 | 0.897 | 0.838 | 0.927
Lympho | 0.934 | 1.000 | 0.934 | 1.000 | 1.000 | 0.934 | 0.962 | 0.923 | 1.000 | 1.000 | 1.000 | 1.000
MNIST | 0.881 | 0.881 | 0.841 | 0.886 | 0.909 | 0.864 | 0.868 | 0.866 | 0.895 | 0.880 | 0.959 | 0.911
Musk | 0.958 | 0.931 | 0.997 | 0.958 | 0.991 | 0.951 | 0.964 | 0.954 | 0.984 | 0.992 | 0.639 | 1.000
OptDigits | 0.952 | 0.952 | 0.952 | 0.949 | 0.952 | 0.953 | 0.949 | 0.955 | 0.951 | 0.952 | 0.819 | 0.981
PenDigits | 0.963 | 0.967 | 0.966 | 0.961 | 0.965 | 0.963 | 0.966 | 0.966 | 0.966 | 0.963 | 0.994 | 0.994
Pima | 0.592 | 0.624 | 0.559 | 0.615 | 0.644 | 0.636 | 0.615 | 0.597 | 0.632 | 0.679 | 0.740 | 0.758
Satellite | 0.681 | 0.757 | 0.813 | 0.634 | 0.716 | 0.597 | 0.732 | 0.709 | 0.761 | 0.761 | 0.841 | 0.845
SatImage | 0.984 | 0.998 | 0.991 | 0.979 | 0.998 | 0.980 | 0.994 | 0.997 | 0.996 | 0.996 | 0.946 | 1.000
Shuttle | 0.933 | 0.992 | 0.968 | 0.900 | 0.973 | 0.891 | 0.990 | 0.973 | 0.981 | 0.978 | 0.990 | 0.998
SpamBase | 0.702 | 0.794 | 0.714 | 0.702 | 0.719 | 0.719 | 0.799 | 0.699 | 0.741 | 0.738 | 0.885 | 0.816
Thyroid | 0.953 | 0.958 | 0.986 | 0.949 | 0.953 | 0.949 | 0.953 | 0.958 | 0.958 | 0.961 | 0.803 | 0.967
Vertebral | 0.750 | 0.778 | 0.817 | 0.796 | 0.810 | 0.817 | 0.817 | 0.817 | 0.817 | 0.819 | 0.807 | 0.904
Vowels | 0.950 | 0.942 | 0.941 | 0.951 | 0.969 | 0.954 | 0.943 | 0.929 | 0.952 | 0.943 | 1.000 | 0.979
WBC | 0.942 | 0.941 | 0.949 | 0.957 | 0.942 | 0.913 | 0.970 | 0.947 | 0.957 | 0.957 | 0.897 | 0.961

The results for all data sets and all algorithms (baseline methods and our
proposed method) are shown in Table 2. The best value for each data set is
highlighted in bold, and the second best value is underlined. In the F1 Score
table, a noticeable difference appears between most of the baseline methods
and our proposed method, with LAKE being a notable exception; however, AD-
DMKDE shows the best performance in the majority of the data sets, followed by
LAKE and LOF methods.

When we consider the influence of the characteristics of the data set on the
performance of our method, we see an interesting pattern. The outlier rate of
the data sets does not seem to have any influence over the performance of AD-
DMKDE, since the method performs well when there are a low number of anomalies
(like in Musk or Shuttle), and also when there are many of them (like in
Cardio or Pima). The number of dimensions have a weak influence, since the
method has the best performance more frequently in data sets with less than
twenty dimensions (such as in Vertebral, Pima or Pendigits), but also in half
of the data sets with more than 100 dimensions (Arrhythmia and Musk). Finally,
when we look at the number of samples, we see that AD-DMKDE stands out more
for data sets with few instances (less than 1000) , and also for many
instances (more than 5000), but its performance is not as good in the
intermediate data sets (the ones that have between 1000 and 5000 instances).
In summary, AD-DMKDE is a method that can perform over the average state-of-
the-art anomaly detection methods, standing out in many different scenarios
and being capable of affording data sets with both low and high values in its
number of dimensions, number of samples and outlier rates.

##  5 Conclusion

In this paper, we presented a novel approach for anomaly detection using
density matrices from quantum mechanics and random Fourier features. The new
method AD-DMKDE was systematically compared against eleven different anomaly
detection algorithms using F1 Score as main metric. AD-DMKDE show better than
state-of-the-art performance over twenty anomaly detection data sets, being
notably superior than classic algorithms and comparable to deep learning-based
methods. The performance of the method does not seem to be affected by the
anomaly rate of the data sets or the size of the data sets, but it seems to
perform better for data sets with low dimensionality of the samples. AD-DMKDE
do not have huge memory requirements due to that the method constructs a
single density matrix in all training phase whose size is defined by the
embedding. So it is possible to resume very large data sets in relatively
small matrices, which makes AD-DMKDE a scalable solution in those cases. Also,
the method allows easy interpretability of the results, because each data
point labeled as “anomaly" can be understood as a sample lying in low density
regions with respect to normal data. As future work, we will continue to
further develop the main concepts of AD-DMKDE, building new algorithms based
on the combination of Fourier features and density matrices with deeper neural
networks, such as autoencoders and variational autoencoders, whose reduction
power can be coupled to the new method in order to better handle data sets
with higher dimensionalities.

## References

  * [1] Aggarwal, C.C.: Outlier analysis second edition (2016)
  * [2] Blázquez-García, A., Conde, A., Mori, U., Lozano, J.A.: A review on outlier/anomaly detection in time series data. ACM Computing Surveys (CSUR) 54(3), 1–33 (2021)
  * [3] Bouyeddou, B., Harrou, F., Kadri, B., Sun, Y.: Detecting network cyber-attacks using an integrated statistical approach. Cluster Computing 24, 1435–1453 (2021)
  * [4] Breunig, M.M., Kriegel, H.P., Ng, R.T., Sander, J.: Lof: identifying density-based local outliers. In: Proceedings of the 2000 ACM SIGMOD international conference on Management of data. pp. 93–104 (2000)
  * [5] Gallego, J.A., González, F.A.: Quantum adaptive fourier features for neural density estimation (2022). https://doi.org/10.48550/ARXIV.2208.00564
  * [6] González, F.A., Gallego, A., Toledo-Cortés, S., Vargas-Calderón, V.: Learning with density matrices and random features. Quantum Machine Intelligence 4(2) (2022)
  * [7] Janssens, J., Huszar, F., Postma, E., van den Herik, H.: Stochastic outlier selection (2012)
  * [8] Kalair, K., Connaughton, C.: Anomaly detection and classification in traffic flow data from fluctuations in the flow–density relationship. Transportation Research Part C: Emerging Technologies 127 (2021), cited By 0
  * [9] Kingma, D.P., Welling, M.: Auto-encoding variational bayes (2014)
  * [10] Li, Z., Zhao, Y., Botta, N., Ionescu, C., Hu, X.: Copod: Copula-based outlier detection. pp. 1118–1123 (11 2020)
  * [11] Liu, F., Huang, X., Chen, Y., Suykens, J.A.K.: Random Features for Kernel Approximation: A Survey on Algorithms, Theory, and Beyond (apr 2020)
  * [12] Liu, F.T., Ting, K.M., Zhou, Z.H.: Isolation forest. In: 2008 eighth ieee international conference on data mining. pp. 413–422. IEEE (2008)
  * [13] Lv, P., Yu, Y., Fan, Y., Tang, X., Tong, X.: Layer-constrained variational autoencoding kernel density estimation model for anomaly detection. Knowledge-Based Systems 196 (2020)
  * [14] Pang, G., Shen, C., Cao, L., Hengel, A.V.D.: Deep learning for anomaly detection: A review. ACM Computing Surveys (CSUR) 54(2), 1–38 (2021)
  * [15] Pevný, T.: Loda: Lightweight on-line detector of anomalies. Machine Learning 102, 275–304 (2016)
  * [16] Prasad, N.R., Almanza-Garcia, S., Lu, T.T.: Anomaly detection. Computers, Materials and Continua 14, 1–22 (2009)
  * [17] Rahimi, A., Recht, B.: Random features for large-scale kernel machines. In: Proceedings of the 20th International Conference on Neural Information Processing Systems. p. 1177–1184. NIPS’07, Curran Associates Inc. (2007)
  * [18] Ramaswamy, S., Rastogi, R., Shim, K.: Efficient algorithms for mining outliers from large data sets. p. 427438. Association for Computing Machinery (2000)
  * [19] Rayana, S.: ODDS library (2016), http://odds.cs.stonybrook.edu
  * [20] Rippel, O., Mertens, P., Konig, E., Merhof, D.: Gaussian anomaly detection by modeling the distribution of normal data in pretrained deep features. IEEE Transactions on Instrumentation and Measurement (2021)
  * [21] Ruff, L., Kauffmann, J.R., Vandermeulen, R.A., Montavon, G., Samek, W., Kloft, M., Dietterich, T.G., Muller, K.R.: A unifying review of deep and shallow anomaly detection. Proceedings of the IEEE 109, 756–795 (2021)
  * [22] Ruff, L., Vandermeulen, R., Goernitz, N., Deecke, L., Siddiqui, S.A., Binder, A., Müller, E., Kloft, M.: Deep one-class classification. vol. 80. PMLR (2018)
  * [23] Santosh, T., Ramesh, D.: Machine learning approach on apache spark for credit card fraud detection. Ingenierie des Systemes d’Information 25, 101–106 (2020)
  * [24] Schölkopf, B., Platt, J.C., Shawe-Taylor, J., Smola, A.J., Williamson, R.C.: Estimating the support of a high-dimensional distribution. Neural computation (2001)
  * [25] Zhang, C., Liu, J., Chen, W., Shi, J., Yao, M., Yan, X., Xu, N., Chen, D.: Unsupervised anomaly detection based on deep autoencoding and clustering. Security and Communication Networks 2021 (2021)
  * [26] Zhao, Y., Nasrullah, Z., Li, Z.: Pyod: A python toolbox for scalable outlier detection. Journal of Machine Learning Research 20 (2019)


# Diophantine equation for the Rice–Mele model:
Topological aspect of filling numbers and associated spatial pump

Koichi Asaga  Takahiro Fukui  Department of Physics, Ibaraki University, Mito
310-8512, Japan

###### Abstract

We introduce a long-period generic spatial modulation into a typical model of
the Thouless pump, namely, the Rice–Mele (RM) model, to examine the lattice
analog of the fermion charge in quantum field theory. We derive a Diophantine
equation relating the fermion charge and the pumped charge, which leads to the
one-dimensional (1D) analog of the Streda formula in the quantum Hall effect
(QHE). This formula implies that an adiabatic change of the periodicity of the
spatial modulation yields a spatial charge pump such that the rightmost charge
is pumped to the right by the Chern number compared with the leftmost charge.
This causes a change in the length of the fermion chain by an integer, thus
providing the opportunity for direct measurement of the Streda formula in 1D
systems.

††preprint: preprint

The discovery of solitons in polyacetylene by Su, Schrieffer, and Heeger had a
great impact on various fields in physics Su et al. (1979). In particular,
fractional fermion charges on solitons attracted considerable attention not
only in condensed-matter physics Takayama et al. (1980); Thouless (1983) but
also in quantum field theory Goldstone and Wilczek (1981); Niemi and Semenoff
(1986). Let us consider $(1+1)$-dimensional fermions coupled with two scaler
fields Goldstone and Wilczek (1981)

$\displaystyle\bar{\psi}(\varphi_{1}+\varphi_{2}i\gamma_{5})\psi\propto\bar{\psi}e^{i\theta\gamma_{5}}\psi,$
(1)

where $\varphi_{j}$ and, hence, the parameter $\theta$ are assumed to be
slowly varying in space and time. Then, the expectation value of the current
$j^{\mu}\equiv\bar{\psi}\gamma^{\mu}\psi$ is given by Goldstone and Wilczek
(1981)

$\displaystyle\langle
j^{\mu}(x,t)\rangle=\frac{1}{2\pi}\epsilon^{\mu\nu}\partial_{\nu}\theta.$ (2)

When $\theta$ is a function of $x$ with the boundary condition
$\theta(+\infty)=\theta(-\infty)+\pi$, the fermion charge $\nu$ becomes
$\nu=-1/2$. This was discussed in terms of the fractional fermion charge in
background soliton fields Goldstone and Wilczek (1981); Niemi and Semenoff
(1986). Generically, Eq. (2) implies that when $\theta$ is a periodic function
of $x$ with $\theta(x+L/2)=\theta(x-L/2)+2\pi p$, the fermion charge in one
period becomes

$\displaystyle\nu\equiv\int_{-L/2}^{L/2}\langle j^{0}(t,x)\rangle dx=-p.$ (3)

The aforementioned fractional charge is regarded as the special case
$L=\infty$ and $p=1/2$. On the other hand, when $\theta$ is a function of $t$
with period $T$, $\theta(t+T)=\theta(t)\pm 2\pi$, the pumped charge $c$
defined by integrating $\langle j^{1}\rangle$ with respect to $t$ over one
period becomes $c=\pm 1$. This is currently known as a Thouless pump Thouless
(1983). Lattice Wilson–Dirac fermions indeed exhibit nontrivial Thouless pumps
Fukui and Fujiwara (2017). Thus, Eq. (2) implies that the fermion charge and
the pumped charge have the same origin.

Generically, a Thouless pump is characterized by the Chern number of a gapped
ground state, which is the change of the 1D polarization as a function of $t$.
Bulk polarization is defined as the Berry phase Vanderbilt and King-Smith
(1993); King-Smith and Vanderbilt (1993); Marzari et al. (2012) at a fixed
$t$. For a system with boundaries, it corresponds to the center of mass (CM)
Wang et al. (2013), and its quantized shift was observed experimentally
Nakajima et al. (2016); Lohse et al. (2016). Recently, quasi-periodic
potentials have been introduced to study disorder effects in the Thouless
pumping Marra and Nitta (2020); Nakajima et al. (2021).

In this study, we generalize the RM model Rice and Mele (1982); Xiao et al.
(2010); Wang et al. (2013) to include long-period spatial modulations and
discuss the fermion charge with topological origin. We show that the fermion
charge and the pumped charge obey a Diophantine equation similar to that for
the Harper equation. This leads to a one-dimensional analog of the Streda
formula in the QHE Streda (1982). Namely, a change of the fermion charge is
due to a change in the periodicity of the spatial modulation. For systems with
charge conservation, a change of the fermion charge induces a change in the
length of the fermion chain. We refer to these phenomena as spatial charge
pumping, which provides a novel opportunity for experimental observation of
topological pumping. It also allows a direct measurement of the Streda formula
in the QHE.

The RM model consists of alternating hoppings as well as staggered potentials.
Let $H(t)=\sum_{j}c_{j}^{\dagger}\hat{\cal H}_{j}(t)c_{j}$ be the Hamiltonian.
Then, the Hamiltonian operator $\hat{\cal H}_{j}$ is written as

$\displaystyle\hat{\cal
H}_{j}(t)=\frac{t_{0}+{\color[rgb]{0,0,0}\sigma^{j}}\delta_{j}(t)}{2}\hat{\delta}+\hat{\delta}^{*}\frac{t_{0}+{\color[rgb]{0,0,0}\sigma^{j}}\delta_{j}(t)}{2}+{\color[rgb]{0,0,0}\sigma^{j}}\Delta_{j}(t),$
(4)

where ${\color[rgb]{0,0,0}\sigma}=-1$, and $\hat{\delta}$ and
$\hat{\delta}^{*}$ are forward and backward shift operators, respectively,
acting on the right, that is, $\hat{\delta}f_{j}=f_{j+1}$ and
$\hat{\delta}^{*}f_{j}=f_{j-1}$. The alternating hopping and staggered
potential are defined by

$\displaystyle\delta_{j}(t)=\delta_{0}\left[\cos 2\pi(t/T+pj/q)+r_{h}\right],$
$\displaystyle\Delta_{j}(t)=\Delta_{0}\left[\sin 2\pi(t/T+pj/q)+r_{m}\right],$
(5)

where $q$ is assumed to be an even integer. These are periodic functions of
$t$ and $j$, as $\delta_{j}(t+T)=\delta_{j}(t)$ and
$\delta_{j+q}(t)=\delta_{j}(t)$, and likewise for $\Delta_{j}(t)$. This model
will be referred to as the generalized RM (gRM) model. When $p=0$, it reduces
to the conventional RM model.

First, we show the relationship between the fermion charge $\nu$ and the
pumped charge $c$. We temporarily assume that $q$ is sufficiently large
compared with $p$. Then, the Hamiltonian depends on $j$ smoothly through
$(p/q)j$ but rapidly through $(-)^{j}$ . To describe the latter, we introduce
a 2-unit cell including two neighboring sites $(2j-1,2j)$ and the
corresponding local charge operator $\rho_{j}\equiv
c_{2j-1}^{\dagger}c_{2j-1}+c_{2j}^{\dagger}c_{2j}$. From the Heisenberg
equation $i\hbar\partial_{t}\rho_{j}=[H(t),\rho_{j}]$, the continuity equation
reads

$\displaystyle\partial_{t}\rho_{j}+({\color[rgb]{0,0,0}I_{j}-I_{j-1}})=0,$ (6)

where the local current operator is defined by
${\color[rgb]{0,0,0}I_{j}}=(i/\hbar)t_{2j}(c_{2j+1}^{\dagger}c_{2j}-\mbox{H.c})$
with $t_{j}=\left(t_{0}+\sigma^{j}\delta_{j}\right)/2$. We note that $t_{2j}$
depends on $\xi\equiv t/T+2pj/q$ smoothly, so that the current operator is
denoted as ${\color[rgb]{0,0,0}I_{j}}(\xi)$. For sufficiently large $q$, we
have $t_{2j-2}\sim t_{2j}-(2p/q)\partial_{\xi}t_{2j}$. Thus, Eq. (6) becomes

$\displaystyle\partial_{t}\rho_{j}+{\color[rgb]{0,0,0}I_{j}(\xi)-I_{j-1}(\xi)}+(2p/q)\partial_{\xi}{\color[rgb]{0,0,0}I_{j-1}}(\xi)=0,$
(7)

where $\xi$ is regarded as a fixed parameter independent of $j$. We define
averaged expectation values $\sum_{j}\langle\rho_{j}\rangle/(N_{\rm
s}/2)\equiv\rho(\xi)$ and,
$\sum_{j}\langle{\color[rgb]{0,0,0}I_{j}}\rangle/(N_{\rm
s}/2)\equiv{\color[rgb]{0,0,0}I(\xi)}$, where $j$ runs over all $N_{\rm s}/2$
unit cells along the chain. Then, we have

$\displaystyle\partial_{t}\rho(\xi)+(2p/q)\partial_{\xi}{\color[rgb]{0,0,0}I(\xi)}=0.$
(8)

It follows from $\partial_{t}=\partial_{\xi}/T$ that

$\displaystyle\rho(\xi)/T+(2p/q){\color[rgb]{0,0,0}I(\xi)}=\mbox{const.}$ (9)

For a large $q$, let us define the fermion charge in the $q$-unit cell, $\nu$,
which is referred to as a filling number:

$\displaystyle\nu$
$\displaystyle=\sum_{j=1}^{q/2}\rho(\xi{\color[rgb]{0,0,0}=t/T+2pj/q})=\frac{q}{2p}\int_{0}^{p}d\xi\rho(\xi)$
$\displaystyle=-T\int_{0}^{p}d\xi{\color[rgb]{0,0,0}I(\xi)}=-\int_{0}^{pT}{\color[rgb]{0,0,0}I(\xi)}dt=-pc,$
(10)

where we have neglected the constant term in Eq. (9), and $c$ denotes the
pumped charge, which is in fact the Chern number Thouless (1983).

|  |
---
Figure 1:  (a) Spectrum at $t=0$ as a function of $\phi=p/q$. The dashed lines
denote the Fermi energies along which the Chern numbers in (b) and (c) are
computed. (b) and (c) show the Chern numbers as a function of $\phi$ at fixed
Fermi energies $0$ and $-0.7$, respectively. The parameters used are
$t_{0}=1$, $\delta_{0}=\Delta_{0}=0.5$, $r_{h}=0$, and $r_{m}=0.2$.

Equation (10) provides a method of computing the Chern number from the filling
number. To this end, it is convenient to define the fermion charge density
$n\equiv\nu/q$ at a given time, for instance, $t=0$. Then, Eq. (10) reads

$\displaystyle\frac{\partial n}{\partial\phi}=-c,\quad\phi=p/q.$ (11)

This serves as the Streda formula in the QHE Streda (1982). It should be noted
that $\nu$ and, hence, $n$ are given at a fixed time. Therefore, one can
compute the Chern numbers without using any time dependence. In Fig. 1 (a),
the spectrum at $t=0$, and in (b) and (c), the Chern numbers computed by Eq.
(11) at fixed Fermi energies, are shown as functions of $\phi$. It can be seen
that when the Fermi energy lies within a gap, the left-hand side of Eq. (11)
is stably constant, from which the Chern number of the gap can be obtained.
The Streda formula (11) is in sharp contrast to the conventional method of
computing Chern numbers using Bloch wave functions in the $k$-$t$ space
Thouless (1983).

Interestingly, Eq. (11) provides an alternative opportunity to observe
topological pumps by controlling $\phi$ in Eq. (5) with $t$ fixed. We assume
that there is a chain with $N$ fermions confined in a harmonic trap. We also
assume that the harmonic trap is so weak that all fermions occupy the bands up
to the $\nu$th, and there is a large gap above the $\nu$th band. Then, the
length of the chain is $L/a=N/n$, where $a$ is the lattice constant, and
$n=\nu/q$ is the charge density. Now, we consider an adiabatic insertion of
$\Delta\phi$ into this system through Eq. (5). This induces a charge-density
change $n\rightarrow n^{\prime}=n+\Delta n$ by Eq. (5), where $\Delta
n=-c\Delta\phi$. Owing to the charge conservation of the present system, the
charge-density change implies a length change $L\rightarrow
L^{\prime}=L+\Delta L$, and hence, we have
$nL/a{\color[rgb]{0,0,0}=n^{\prime}L^{\prime}/a}=(n+\Delta n)(L+\Delta L)/a$,
which yields

$\displaystyle\Delta L/L=(c/n)\Delta\phi.$ (12)

This equation implies that the length of the fermion chain becomes longer by
the amount $\Delta L$ because of the adiabatic insertion of $\Delta\phi$. This
is one of our main results. Let us apply this to the conventional RM model.
When $\phi=0$, the RM model yields $c=1$ at $n=1/2$. As above, we assume a
chain of $N$ fermions of length $L/a=2N$ in a harmonic trap. The gap is also
assumed to be so large that all fermions occupy the lower band. This was
indeed realized in the experiment Nakajima et al. (2016). We note that in the
case of $\phi=0$, we can choose any $q$; thus, let us set $q=L/a$. Then, an
adiabatic $\Delta\phi=p/q$ insertion yields $\Delta L/a=2p$, implying that the
length of the chain becomes longer by $2p$. For a more generic choice of $q$,
see below Eq. (14). We can interpret this as a spatial charge pump: The charge
at the right end is pumped $2p$ to the right compared with the charge at the
left end. This causes the change in the length of the system. On the other
hand, in the Thouless pump, each charge is pumped in the same direction by the
same amount under time evolution, so that the length of the system does not
change, but the CM is shifted instead. In this sense, we refer to the Thouless
pump as a temporal pump to distinguish it from the aforementioned spatial pump
owing to the adiabatic change of $\phi$.

|
---|---
|
|
Figure 2:  (a) and (b) show energy spectra as functions of $p$ with $q=100$
fixed. The length $L/a$ of the system is chosen as $L/a=q$, and open boundary
conditions are imposed on both ends. Dashed lines denote the same Fermi
energies as in Fig. 1, which are located in gaps with $c=1$ and $2$; (c) and
(d) show the number of states below the Fermi energies specified in (a) and
(b), respectively; (e) and (f) show the CM of the occupied states in (a) and
(b), respectively. The black dots indicate the values of the jumps at the
discontinuity of the CM.

We can also consider the spatial pump from the perspective of the bulk-edge
correspondence. Specifically, we consider the same system mentioned above,
that is, the system of density $n$ and size $L$ that includes $N=nL/a$
fermions, on which the open boundary conditions are imposed. Now, we assume
that the system is in contact with a particle reservoir. It follows from Eq.
(11) that the adiabatic $\Delta\phi=\Delta p/q$ ($q=L/a$) insertion causes a
change of the total fermion charge,

$\displaystyle\Delta N=-c\Delta p.$ (13)

Such a change is due to the spectral flow of the edge states across the Fermi
energy when $p$ is varied. In Figs. 2 (a) and (b), the spectra of the systems
of size $L/a=q=100$ as functions of $p$ are shown. It is clearly seen that the
number of states below the Fermi energy only decreases by $c$ for every
$\Delta p=1$ because of the edge states across the Fermi energy. Furthermore,
these edge states are localized at the right end, as the CM exhibits
discontinuous changes by $-1/2$ Hatsugai and Fukui (2016), as shown in Figs. 2
(e) and (f). Here, the CM is defined to be $-1/2$ at $j=1$, and $+1/2$ at
$j=L/a$.

So far, we have discussed the spatial pump based on the filling number (10)
and the resultant Streda formula (11). To apply Eq. (10) to the lattice model
more precisely, some corrections and definitions are required. First, when
$p=0$, $\nu$ is the number of occupied bands for the RM model; completely
unoccupied (0), half-filled (1), and fully occupied (2). In the gRM model in
Eq. (5), these filling numbers are translated into $0$, $q/2$, and $q$. This
should be included in Eq. (10), probably as the constant term in Eq. (9)
neglected in Eq. (10). Second, the gRM model has $q$ bands generically, so
that Eq. (10) is expected to hold at each gap separately. Let us focus on the
$i$th gap. The filling number $\nu_{i}$ is the number of bands below the $i$th
gap, and the pumped charge $c_{i}$ is the Chern number of the gap, that is,
the sum of the Chern numbers of all bands below the $i$th gap. With these
corrections and definitions, we propose the following relationship between
$\nu_{i}$ and $c_{i}$:

$\displaystyle\nu_{i}=s_{i}q/2-pc_{i},$ (14)

where $s_{i}$ denotes a certain integer. This equation is a Diophantine
equation similar to that Thouless et al. (1982) for the Harper equation Harper
(1955). Using this, Eq. (12) provides a more precise length of the system as
follows: Let us begin with a model including $\phi_{0}$, and control the
fermion number $N$ such that all fermions occupy the bands below the $i$th
gap, implying the charge density $n_{i}=s_{i}/2-\phi_{0}c_{i}$. This may be
possible when the $i$th gap is sufficiently large. Then, the initial length of
the chain is $L_{0}/a=N/n_{i}$. We assume that we change
$\phi_{0}\rightarrow\phi$ and obtain $n_{i^{\prime}}=s_{i^{\prime}}/2-\phi
c_{i^{\prime}}$. The length becomes $L/a=N/n_{i^{\prime}}$. It then follows
that

$\displaystyle L/L_{0}=(s_{i}-2c_{i}\phi_{0})/(s_{i}-2c_{i}\phi),$ (15)

provided that $c_{i^{\prime}}=c_{i}$ and $s_{i^{\prime}}=s_{i}$, implying that
all fermions continue to occupy the bands up to the same topological gap
characterized by $c_{i}$ and $s_{i}$. For the conventional half-filled RM
model discussed below Eq. (12), we have $L/L_{0}=1/(1-2\phi)$, where $L_{0}$
is the length of the system when $\phi_{0}=0$. For a small $\phi$, this
relation is approximated by $L=2\phi L_{0}$, and if we choose $q$ in
$\phi=p/q$ as $q=N/2(=L_{0}/a)$, the previous result can be reproduced.

|
---|---
Figure 3:  Chern numbers of gaps as functions of the Fermi energy $\epsilon$
for the model with $q=22$ and (a) $p=1$, (b) $p=-1$. Horizontal colored thick
lines indicate gaps, and the number above the $i$th gap is $\nu_{i}$. The
green, red, and cyan lines correspond to the gaps satisfying Eq. (14) with
$s_{i}=0,1$, and $2$, respectively. The parameters of the model are same as
those in Fig. 1.

To verify the validity of the formula (14), the filling number $\nu_{i}$ as
well as the Chern numbers $c_{i}$ Thouless (1983) are shown in Fig. 3 in the
case of $p/q=\pm 1/22$. Here, the Chern numbers are those defined by the Berry
curvature of the Bloch states in the $k$-$t$ space Thouless et al. (1982);
Thouless (1983), calculated using the method in Fukui et al. (2005). It can be
seen that the system has $q$ separated bands: Each band has Chern number $\pm
1$, and as a result, the Chern number of the gap as a function of energy
exhibits a step-like behavior similar to that of the Landau levels in the QHE.
Such a step-like structure is divided into three regions by the van Hove
singularities, which carry large Chern numbers Hatsugai et al. (2006). Each
region is characterized by $s_{i}$: $s_{i}=0$ $(2)$ series is associated with
particle (hole) states in the fully-unoccupied (fully-occupied) vacuum,
whereas $s_{j}=1$ series is associated with the half-filled ground state of
the RM model.

|
---|---
Figure 4:  Same as in Fig. 3 but for $p=3$. In (a), the orange lines and the
black line correspond to $s_{i}=3$ and $s_{j}=-1$, respectively. Other colors
have the same interpretation as in Fig. 3.  It can be seen that the Chern
number of each band, calculated by the difference of the gap Chern numbers
$c_{i}-c_{i-1}$, is $-4,7,-4,7,-4,-4,\cdots$, starting from the 1st band; (b)
is the same system as (a) but uses a larger threshold gap: Here, the threshold
gap implies that any two bands are regarded as degenerate if the minimum
direct gaps between them are smaller than the threshold gap.

Other examples are shown in Fig. 4. When $p=3$, each band has large Chern
numbers, namely, 7 or $-4$, as seen in Fig. 4 (a), and gaps with various
$s_{i}$ coexist in the same energy region. Nevertheless, if tiny gaps are
ignored and regarded as degenerate bands, as shown in Fig. 4 (b), three series
of gaps separated by van Hove singularities can be seen, as in the case of
$p=1$ in Fig. 3. These examples as well as various other numerical studies
involving small $q$ and/or large $p$ demonstrate that the Diophantine equation
(14) is valid generically for the gRM model.

The effect of the spatial pump can be understood as a change of the filling
number $\nu_{i}$ for nonzero $\phi$: In Figs. 3 and 4,  it can be seen that
the the gaps at zero energy have filling numbers $\nu=10$, $12$, and $8$ for
$p=1$, $-1$, and $3$, respectively, as well as the same Chern number $c=1$.
This gap corresponds to the half-filled gap with $\nu=11$, when $\phi=0$
($p=0$ and $q=22$). Therefore, in all these cases, we have the filling-number
shift $\Delta\nu=\pm 1,-3=-c\Delta p$.  This is the effect of the spatial pump
from the bulk perspective.

As is the case with the RM model Wang et al. (2013), the present system would
be realized in cold atoms using the optical superlattice Nakajima et al.
(2016); Lohse et al. (2016)

$\displaystyle V(x,\phi,t)=$ $\displaystyle-v_{1}\cos^{2}\frac{\pi
x}{a}-v_{2}\cos^{2}\left[\left(\frac{1}{2}+\phi\right)\frac{\pi
x}{a}+\frac{\pi t}{T}\right]$
$\displaystyle{\color[rgb]{0,0,0}-v_{3}\cos^{2}\left(\frac{\pi
x}{2a}+\theta\right),}$ (16)

where $\theta$ controls the parameters $r_{h}$ and $r_{m}$ in Eq. (5). An
adiabatic change of $t$ with $\phi=0$ fixed yields the temporal Thouless pump
already observed Nakajima et al. (2016); Lohse et al. (2016), whereas an
adiabatic change of $\phi$ with $t$ fixed would yield the spatial pump
proposed in this paper. Here, $\phi$ could be controlled by the change in the
tilted angle between the optical superlattices Marra and Nitta (2020);
Nakajima et al. (2021). In particular, the latter could provide a direct
measurement of the Streda formula in 1D systems.

In conclusion, we discussed the topological aspect of the filling number or
fermion charge of the gRM model. We derived a Diophantine equation relating
the filling number and the pumped charge, and the resultant Streda formula.
These predict a spatial pump, which could be observed as a change in the
length of the fermion chain. A simpler model can be considered if
${\color[rgb]{0,0,0}\sigma=1}$ is chosen in Eq. (4). This model is quite
analogous to the conventional QHE, as Eq. (14) is modified into
$\nu_{i}=s_{i}q-pc_{i}$, which is the same as the Diophantine equation for the
Harper equation. Such systems with ${\color[rgb]{0,0,0}\sigma=\pm 1}$ would
offer a promising platform for studying more intimate relationships between
topological pumps and the QHE.

TF would like to thank Y. Takahashi for fruitful discussions. This work was
supported in part by Grants-in-Aid for Scientific Research Numbers 17H06138
from Japan Society for the Promotion of Science.

## References

  * Su et al. (1979) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Physical Review Letters 42, 1698 (1979).
  * Takayama et al. (1980) H. Takayama, Y. R. Lin-Liu, and K. Maki, Physical Review B 21, 2388 (1980).
  * Thouless (1983) D. J. Thouless, Physical Review B 27, 6083 (1983).
  * Goldstone and Wilczek (1981) J. Goldstone and F. Wilczek, Physical Review Letters 47, 986 (1981).
  * Niemi and Semenoff (1986) A. J. Niemi and G. W. Semenoff, Phys. Rep. 135, 99 (1986).
  * Fukui and Fujiwara (2017) T. Fukui and T. Fujiwara, Physical Review B 96, 205404 (2017).
  * Vanderbilt and King-Smith (1993) D. Vanderbilt and R. D. King-Smith, Physical Review B 48, 4442 (1993).
  * King-Smith and Vanderbilt (1993) R. D. King-Smith and D. Vanderbilt, Physical Review B 47, 1651 (1993).
  * Marzari et al. (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Reviews of Modern Physics 84, 1419 (2012).
  * Wang et al. (2013) L. Wang, M. Troyer, and X. Dai, Physical Review Letters 111, 026802 (2013).
  * Nakajima et al. (2016) S. Nakajima, T. Tomita, S. Taie, T. Ichinose, H. Ozawa, L. Wang, M. Troyer, and Y. Takahashi, Nat Phys 12, 296 (2016).
  * Lohse et al. (2016) M. Lohse, C. Schweizer, O. Zilberberg, M. Aidelsburger, and I. Bloch, Nat Phys 12, 350 (2016).
  * Marra and Nitta (2020) P. Marra and M. Nitta, Physical Review Research 2, 042035 (2020).
  * Nakajima et al. (2021) S. Nakajima, N. Takei, K. Sakuma, Y. Kuno, P. Marra, and Y. Takahashi, Nature Physics (2021).
  * Rice and Mele (1982) M. J. Rice and E. J. Mele, Phys. Rev. Lett. 49, 1455 (1982).
  * Xiao et al. (2010) D. Xiao, M.-C. Chang, and Q. Niu, Reviews of Modern Physics 82, 1959 (2010).
  * Streda (1982) P. Streda, Journal of Physics C: Solid State Physics 15, L717 (1982).
  * Hatsugai and Fukui (2016) Y. Hatsugai and T. Fukui, Physical Review B 94, 041102 (2016).
  * Thouless et al. (1982) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Physical Review Letters 49, 405 (1982).
  * Harper (1955) P. G. Harper, Proceedings of the Physical Society. Section A 68, 874 (1955).
  * Fukui et al. (2005) T. Fukui, Y. Hatsugai, and H. Suzuki, Journal of the Physical Society of Japan 74, 1674 (2005).
  * Hatsugai et al. (2006) Y. Hatsugai, T. Fukui, and H. Aoki, Physical Review B 74, 205414 (2006).