# Conversion between electron spin and microscopic atomic rotation

Masato Hamada  Department of Physics, Tokyo Institute of Technology, 2-12-1
Ookayama, Meguro-ku, Tokyo 152-8551, Japan  Shuichi Murakami  Department of
Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo
152-8551, Japan  TIES, Tokyo Institute of Technology, Ookayama, Meguro-ku,
Tokyo 152-8551, Japan

###### Abstract

We theoretically investigate the microscopic mechanism of conversion between
the electron spin and the microscopic local rotation of atoms in crystals. In
phonon modes with angular momenta, the atoms microscopically rotate around
their equilibrium positions in crystals. In a simple toy model with phonons,
we calculate the spin expectation value by using the adiabatic series
expansion. We show that the time-averaged spin magnetization is generated by
the microscopic local rotation of atoms via the spin-orbit interaction. On the
other hand, in the system with a simple vibration of atoms, time-averaged spin
magnetization becomes zero due to the time-reversal symmetry. Moreover, the
magnitude of the time-averaged spin magnetization depends on the inverse of
the difference of instantaneous eigenenergy, and we show that it becomes
smaller in band insulators with a larger gap.

##  I Introduction

The conversion between the magnetization and the mechanical rotation, such as
Einstein-de Haas effect Einstein and de Haas (1915) and Barnett effect Barnett
(1915), has been known since more than a hundred years ago, and this has been
focused on in spintronics recently. These effects have been understood by
means of the spin-rotation coupling Matsuo  _et al._ (2013, 2011a, 2011b),
which means that the mechanical rotation acts on electron spins as an
effective magnetic field. In spintronics, mechanical generation of spin
current has been theoretically proposed or experimentally reported by several
means, such as the surface acoustic wave Matsuo  _et al._ (2013); Kobayashi
_et al._ (2017), the twisting vibrational mode in carbon nanotubes Hamada  _et
al._ (2015), and the flow of the liquid metals Takahashi  _et al._ (2016), via
the spin-rotation coupling. However, microscopic mechanism of conversion
between electron spins and the rotational motion of atoms in solids is not
well known.

Recently, an angular momentum of phonons, which is a microscopic local
rotation of atoms in crystals, has been formulated Zhang and Niu (2014), and
various phenomena related with it have been energetically investigated, such
as the phonon Hall effect Strohm  _et al._ (2005); Inyushkin and Taldenkov
(2007); Sugii  _et al._ (2017); Kagan and Maksimov (2008); Zhang  _et al._
(2010), spin relaxation Garanin and Chudnovsky (2015); Nakane and Kohno
(2018); Streib  _et al._ (2018), the orbital magnetic moment of phonons
Juraschek and Spaldin (2019), and conversion between phonons and magnons
Holanda  _et al._ (2018); Guerreiro and Rezende (2015). In a system without
inversion symmetry, the phonons have chirality Zhang and Niu (2015). In
particular, at the valleys in momentum space, the angular momentum of phonons
is quantized, and such phonons are called chiral phonons. Chiral phonons can
be excited by a circularly polarized light via intervalley scattering of
electrons, and this has been experimentally observed in monolayer tungsten
diselenide Zhu  _et al._ (2018). Methods of generation of the phonon angular
momentum have been proposed, such as a magnetic field Zhang and Niu (2014);
Juraschek  _et al._ (2017), a Coriolis force in rotating frame Wang  _et al._
(2015); Liu  _et al._ (2017), a circularly polarized light Zhang and Niu
(2015); Zhu  _et al._ (2018), an infrared excitation Juraschek and Spaldin
(2019); Nova  _et al._ (2017), a temperature gradient in the system without
inversion symmetry Hamada  _et al._ (2018), and an electric field in the
magnetic insulator Mas . Here, one question arises whether electron spins
couple with a microscopic local rotation of atoms in crystals. Because the
microscopic local rotation is quite different from the macroscopic rotation
studied in the Einstein-de Haas effect and the Barnett effect, it is a crucial
issue whether one can apply the spin-rotation coupling to a coupling between
the electron spins and a microscopic local rotations.

In this paper, we theoretically propose microscopic mechanism of conversion
between an electron spin and a microscopic local rotation of atoms in
crystals. We first introduce a simple toy model, which is a two-dimensional
system forming the honeycomb lattice with a microscopic local rotation of
atoms. We set the rotation of atoms to be periodic in time, and we calculate
the spin expectation value. Here we assume that the rotational motion is much
slower than the characteristic energy scale of electrons, and we use the
adiabatic series expansion in two previous theories; Berry phase in adiabatic
process Trifunovic  _et al._ (2019) and the iteration of the unitary
transformation Berry (1987); Rigolin  _et al._ (2008) in order to calculate a
time-averaged spin expectation value. By using the adiabatic series expansion,
we show that the microscopic local rotation of atoms leads to a nonzero spin
magnetization even after averaging over one period in our toy model. Meanwhile
a simple vibration of atoms induces spin magnetization as a function of time,
but after time averaging it becomes zero. Moreover, we discuss dependence on
the parameters for time-averaged spin magnetization.

##  II Toy model with phonon angular momentum

We introduce a two-dimensional tight-binding model on the honeycomb lattice
with a microscopic local rotation of atoms as shown in Fig. 1(a). In Fig.
1(b), we show the first Brillouin zone. The two-dimensional honeycomb lattice
without inversion symmetry has chiral phonons at the valleys $K$ and
$K^{\prime}$ in the momentum space Zhang and Niu (2015). The $K$ and
$K^{\prime}$ phonons have the chirality $\pm 1$, and in these modes, the atoms
circularly rotate. On the other hand, at the zone center $\Gamma$, the
longitudinal and the transverse optical phonon modes are degenerate and these
modes do not have phonon angular momentum. However, we can get the circular
polarized phonon modes by superposition of these degenerate modes at the
$\Gamma$ point. Here, we focus on a microscopic local rotation of atoms
corresponding to the $\Gamma$-phonon by superposition of two optical
degenerate modes. In Fig. 1(a), the red and blue balls represent the A and B
atoms at the two sublattices respectively, and the atoms A and B circularly
rotate in the $xy$-plane, with their phases of rotations different by $\pi$.
The primitive vectors are written as $\bm{a}_{1}=a(1,0)$,
$\bm{a}_{2}=a(1/2,\sqrt{3}/2)$ with a lattice constant $a$, and the vectors
between the nearest neighbor atoms are written as
$\bm{d}_{1}=a_{0}(\sqrt{3}/2,1/2)$, $\bm{d}_{2}=a_{0}(-\sqrt{3}/2,1/2)$,
$\bm{d}_{3}=a_{0}(0,-1)$ with the nearest neighbor bond length
$a_{0}=a/\sqrt{3}$. We set the angular velocity of the atoms A and B to be
$\Omega$, and the displacement vectors of atoms A and B at time $t$ are
represented as

$\displaystyle\bm{u}_{\rm A}=u_{0}^{\rm A}(\cos\Omega t,\sin\Omega
t),~{}\bm{u}_{\rm B}=-u_{0}^{\rm B}(\cos\Omega t,\sin\Omega t),$ (1)

and the displacement vector from the atom A to the atom B is

$\bm{u}=\bm{u}_{\rm B}-\bm{u}_{\rm A}=-u_{+}(\cos\Omega t,\sin\Omega t)$ (2)

with $u_{+}=u_{0}^{\rm A}+u_{0}^{\rm B}$.

Next, we consider the following nearest-neighbor tight-binding Hamiltonian for
electrons, which is similar to the Kane-Mele model , representing a two-
dimensional topological insulator Kane and Mele (2005),

$\hat{H}_{0}=t_{0}\sum_{\braket{ij}}\hat{c}_{i}^{\dagger}\hat{c}_{j}+\lambda_{\rm
v}\sum_{i}\xi_{i}\hat{c}_{i}^{\dagger}\hat{c}_{i}+i\frac{\lambda_{\rm
R}}{a_{0}}\sum_{\braket{ij}}\hat{c}_{i}^{\dagger}(\bm{s}\times\bm{d}_{ij})_{z}\hat{c}_{j}.$
(3)

The first term is a nearest-neighbor hopping term and $t_{0}$ is a hopping
parameter. The second term is a staggered potential, which we include to break
inversion symmetry. $\xi_{i}$ is $\xi_{\rm A(B)}=+1(-1)$ for the A and B atoms
and $\lambda_{\rm v}$ is an on-site potential. The third term is a nearest-
neighbor Rashba term, which leads to spin-split band structure, and
$\lambda_{\rm R}$ is Rashba parameter, $\bm{d}_{ij}$ is the nearest neighbor
bond vector, and $s_{\alpha}$ is the Pauli matrices of electron spin. The
operators
$\hat{c}_{i}^{\dagger}=(c_{i\uparrow}^{\dagger},\hat{c}_{i\downarrow}^{\dagger})$,
$\hat{c}_{i}=(\hat{c}_{i\uparrow},\hat{c}_{j\downarrow})^{T}$ represent the
creation and annihilation operators for an electron at the $i$ site,
respectively. By using the Bloch wave function $\ket{\epsilon_{n}(\bm{k})}$,
this Hamiltonian $\hat{H}_{0}$ can be diagonalized, and the Bloch Hamiltonian
$\hat{H}_{0}(\bm{k})$ in the $\bm{k}$-space is represented as

$\displaystyle\hat{H}_{0}(\bm{k})$ $\displaystyle=t_{0}(F_{\rm
R}(\bm{k})\Gamma^{1}-F_{\rm Im}(\bm{k})\Gamma^{12})+\lambda_{\rm v}\Gamma^{2}$
$\displaystyle-\lambda_{\rm
R}(F_{1}(\bm{k})\Gamma^{4}+F_{2}(\bm{k})\Gamma^{3}+F_{3}(\bm{k})\Gamma^{24}-F_{4}(\bm{k})\Gamma^{23}),$
(4)

and

$\begin{split}F_{\rm R}(\bm{k})&=2\cos x\cos y+\cos 2y,\\\ F_{\rm
Im}(\bm{k})&=2\cos x\sin y-\sin 2y,\\\ F_{1}(\bm{k})&=\sqrt{3}\sin x\sin y,\\\
F_{2}(\bm{k})&=\cos x\cos y-\cos 2y,\\\ F_{3}(\bm{k})&=\sqrt{3}\sin x\cos
y,\\\ F_{4}(\bm{k})&=\cos x\sin y+\sin 2y,\end{split}$ (5)

with $x=\sqrt{3}k_{x}a_{0}/2$, $y=k_{y}a_{0}/2$. To express $H_{0}(\bm{k})$ in
a compact form, we introduced the five Dirac matrices:
$\Gamma^{1,2,3,4,5}=(\sigma_{x}\otimes s_{0},\sigma_{z}\otimes
s_{0},\sigma_{y}\otimes s_{x},\sigma_{y}\otimes s_{y},\sigma_{y}\otimes
s_{z})$, where the Pauli matrices $\sigma_{\alpha}$ and $s_{\alpha}$ represent
the sublattice index and the electron spin respectively, and
$\sigma_{0},s_{0}$ are unit matrices. We also introduced their ten commutators
$\Gamma^{ab}=[\Gamma^{a},\Gamma^{b}]/(2i)$.

Figure 1:  Schematic figures of (a) the honeycomb lattice with a microscopic
local rotation and (b) the first Brillouin zone.

Next, we consider a perturbation term due to the microscopic local rotation of
atoms. Since the nearest-neighbor bond lengths change with time $t$ by the
microscopic local rotations, we take these rotations as modulation of the
nearest-neighbor hopping and the Rashba parameters. We assume that the lattice
deformation $\bm{u}$ modifies the nearest-neighbor hopping parameter as
$t_{0}\to t_{0}+\delta t_{a}~{}(a=1,2,3)$ Sasaki and Saito (2008). We define
the Hamiltonian representing the modulation of the nearest-neighbor hopping
term as

$\hat{H}_{\rm h}=\sum_{\braket{ij}}\delta
t_{ij}(t)\hat{c}_{i}^{\dagger}\hat{c}_{j}.$ (6)

We assume that the modulation of the nearest-neighbor hopping parameter
$\delta t_{a}~{}(a=1,2,3)$ is proportional to the change of the nearest-
neighbor bond length, which is given by $\bm{u}(t)\cdot(\bm{d}_{a}/a_{0})$.
Thus, the modulation of the nearest-neighbor hopping parameter $\delta
t_{a}(t)$ is given by

$\delta t_{a}(t)=-\frac{\delta t_{0}}{a_{0}^{2}}\bm{u}(t)\cdot\bm{d}_{a}.$ (7)

Therefore, the Hamiltonian for the modulation of the hopping term for the wave
vector $\bm{k}$ is represented as

$\displaystyle\hat{H}_{\rm h}(\bm{k},t)$ $\displaystyle=-\delta
t(F_{1}(\bm{k})\Gamma^{1}+F_{3}(\bm{k})\Gamma^{12})\cos\Omega t$
$\displaystyle-\delta
t(F_{2}(\bm{k})\Gamma^{1}-F_{4}(\bm{k})\Gamma^{12})\sin\Omega t,$ (8)

with $\delta t=\delta t_{0}u_{+}/a_{0}$. Similarly, we consider the modulation
of the nearest-neighbor Rashba term. We assume that the Rashba term is
modulated as

$\hat{H}_{\rm R}(t)=i\lambda_{\rm
R}\sum_{\braket{ij}}\hat{c}_{i}^{\dagger}\left(\bm{s}\times\frac{\bm{d}_{ij}^{\prime}}{|d_{ij}^{\prime}|}\right)_{z}\hat{c}_{j},$
(9)

where $\bm{d}_{ij}^{\prime}$ is the nearest neighbor bond vector with the
modulation $\bm{u}(t)$. This can be expanded in Taylor series in powers of the
displacement vector $\bm{u}$

$\displaystyle\hat{H}_{\rm R}(t)$ $\displaystyle=\frac{i\lambda_{\rm
R}}{a_{0}}\sum_{\braket{ij}}\hat{c}_{i}^{\dagger}(\bm{s}\times\bm{d}_{ij})_{z}\hat{c}_{j}$
$\displaystyle+\frac{-i\lambda_{\rm
R}}{a_{0}^{3}}\sum_{\braket{ij}}(\bm{d}_{ij}\cdot\bm{u}(t))\hat{c}_{i}^{\dagger}(\bm{s}\times\bm{d}_{ij})_{z}\hat{c}_{j}$
$\displaystyle+\frac{i\lambda_{\rm
R}}{a_{0}}\sum_{\braket{ij}}\hat{c}_{i}^{\dagger}(\bm{s}\times\bm{u}(t))_{z}\hat{c}_{j}+O(\bm{u}^{2}).$
(10)

The first term is the Rashba term without the microscopic local rotation and
the second and the third terms are the modulated Rashba terms by the
displacement vector $\bm{u}$. The modulated Rashba term $H_{\rm R}(t)$ at the
wave vector $\bm{k}$ is represented as

$\displaystyle\hat{H}_{\rm R}(\bm{k},t)$ $\displaystyle=-\delta\lambda_{\rm
R}\Big{[}F_{3}(\bm{k})\Gamma^{23}+F_{1}(\bm{k})\Gamma^{3}$
$\displaystyle\hskip 28.45274pt-(2F_{\rm Im}(\bm{k})+3\cos x\sin
y)\Gamma^{24}$ $\displaystyle\hskip 28.45274pt+(2F_{\rm R}(\bm{k})+3\cos x\cos
y)\Gamma^{4}\Big{]}\cos\Omega t$ $\displaystyle-\delta\lambda_{\rm
R}\Big{[}(2F_{\rm Im}(\bm{k})+F_{5}(\bm{k}))\Gamma^{23}$ $\displaystyle\hskip
28.45274pt-(2F_{\rm R}(\bm{k})+F_{6}(\bm{k}))\Gamma^{3}$ $\displaystyle\hskip
28.45274pt-F_{3}(\bm{k})\Gamma^{24}-F_{1}(\bm{k})\Gamma^{4}\Big{]}\sin\Omega
t,$ (11)

and

$\begin{split}F_{5}(\bm{k})&=\cos x\sin y-2\sin 2y,\\\ F_{6}(\bm{k})&=\cos
x\cos y+2\cos 2y,\end{split}$ (12)

with $\delta\lambda_{\rm R}=\frac{\lambda_{\rm R}u_{+}}{2a_{0}}$. Therefore,
the Hamiltonian with periodic microscopic local rotation of atoms for the wave
vector $\bm{k}$ is represented as

$\hat{H}(\bm{k},t)=\hat{H}_{0}(\bm{k})+\hat{H}_{\rm h}(\bm{k},t)+\hat{H}_{\rm
R}(\bm{k},t).$ (13)

##  III Adiabatic series expansion

Now we calculate spin expectation values for the time-periodic Hamiltonian. In
general, the motions of atoms due to phonons are much slower than those of
electrons because the masses of atoms are much larger than those of electrons.
Therefore, in this section, in order to calculate the spin expectation values
we employ the adiabatic series expansion in two ways: Berry-phase method in
adiabatic process and iteration of unitary transformation. The use of the
adiabatic expansion is justified when the band gap is much larger than
$\hbar\Omega$

###  III.1 Berry-phase method in adiabatic process

In this section, we explain the Berry-phase method to calculate dynamical
response in adiabatic process proposed in Ref. Trifunovic  _et al._ (2019). We
will calculate an expectation values of an operator $\hat{X}(t)$. To this end,
we introduce the Hamiltonian $\hat{H}(t,\lambda)$ such that

$\displaystyle\hat{H}(t,\lambda)|_{\lambda=0}$ $\displaystyle=\hat{H}(t)$ (14)
$\displaystyle\partial_{\lambda}\hat{H}(t,\lambda)|_{\lambda=0}$
$\displaystyle=\hat{X}(t).$ (15)

For example, in static systems, when the parameter $\lambda$ is an external
magnetic field $B_{z}$, the quantity $\hat{X}$ is magnetization $M_{z}$. In
particular, when the magnetic field enters the Hamiltonian as a vector
potential, the quantity $\hat{X}$ is the orbital magnetization. On the other
hand, when the magnetic field enters the Hamiltonian as a Zeeman field,
$\hat{X}$ is the spin magnetization.

Next, we calculate the time-averaged expectation values of the quantity
$\hat{X}(t)$, under the time-periodic Hamiltonian $\hat{H}(t)$ with time
period $T$. Let $\ket{\epsilon_{n}(t)}$ be an instantaneous eigenstate of the
band $n$ at time $t$ with instantaneous eigenvalue $\epsilon_{n}(t)$, where
$n$ is the band index. We assume that time-dependence of the Hamiltonian is
slow enough so that the state hardly transits to other states, which means
that the gap $\Delta_{nm}$ between bands $n$ and $m$ satisfies
$\Delta_{nm}T\gg\hbar$ with time period $T$. We introduce the Neumann equation
$i\hbar\partial_{t}\hat{\rho}(t)=[\hat{H}(t),\hat{\rho}(t)]$ with the density
matrix $\hat{\rho}(t)$. By using the density matrix $\hat{\rho}(t)$, we define
the time-average of the expectation value $\hat{X}(t)$ as

$X\equiv\frac{1}{T}\int_{0}^{T}dt~{}{\rm tr}[\hat{\rho}(t)\hat{X}(t)].$ (16)

We suppose that in the absence of the time evolution the density matrix is
identical to $\ket{\epsilon_{n}(t)}\bra{\epsilon_{n}(t)}$ corresponding to the
band $n$. The density matrix $\hat{\rho}(t)$ can be expanded to the first
order with respect to $(\Delta_{nm}T)^{-1}$ as

$\hat{\rho}(t)=\ket{\epsilon_{n}(t)}\bra{\epsilon_{n}(t)}+\delta\hat{\rho}(t)+O((\Delta_{nm}T)^{-2}).$
(17)

The first term is the projector onto the instantaneous eigenstate in the
absence of the time evolution. The second term describes transitions to the
excited states through the time evolution. The matrix element between bands
$n$ and $m$ is given by

$\displaystyle\bra{\epsilon_{m}(t)}\hat{\rho}(t)\ket{\epsilon_{n}(t)}=i\hbar\frac{\braket{\epsilon_{m}(t)}{\partial_{t}\epsilon_{n}(t)}}{\epsilon_{m}(t)-\epsilon_{n}(t)},$
(18)

from the Neumann equation and Eq. (17). From Eq. (15), the matrix elements of
$\hat{X}(t)$ between bands $n$ and $m$ is written as

$\displaystyle\bra{\epsilon_{m}(t)}\hat{X}(t)\ket{\epsilon_{n}(t)}$
$\displaystyle=\bra{\epsilon_{m}(t)}\partial_{\lambda}\hat{H}(t,\lambda)|_{\lambda=0}\ket{\epsilon_{n}(t)}$
$\displaystyle=(\epsilon_{n}(t)-\epsilon_{m}(t))\braket{\epsilon_{m}(t)}{\partial_{\lambda}\epsilon_{n}(t,\lambda)}|_{\lambda=0}.$
(19)

Here we used the Sternheimer identity
$(\partial_{\lambda}\hat{H}(t,\lambda))\ket{\epsilon_{n}(t,\lambda)}=(\epsilon_{n}(t,\lambda)-\hat{H}(t,\lambda))\ket{\partial_{\lambda}\epsilon_{n}(t,\lambda)}+\partial_{\lambda}\epsilon_{n}(t,\lambda)\ket{\epsilon_{n}(t,\lambda)}$.
By combining these equations, Eq. (16) is rewritten as

$\displaystyle
X_{n}=\frac{1}{T}\int_{0}^{T}dt(\partial_{\lambda}\epsilon_{n}(t,\lambda)+B_{n}(t,\lambda))|_{\lambda=0},$
(20)

where
$B_{n}(t,\lambda)=i\hbar\partial_{t}\braket{\epsilon_{n}(t,\lambda)}{\partial_{\lambda}\epsilon_{n}(t,\lambda)}-i\hbar\partial_{\lambda}\braket{\epsilon_{n}(t,\lambda)}{\partial_{t}\epsilon_{n}(t,\lambda)}$
is the Berry curvature of a band $n$ in the $(t,\lambda)$ space. Moreover, in
the time-periodic systems, the instantaneous eigenstates can be chosen to
satisfy $\ket{\epsilon_{n}(t,\lambda)}=\ket{\epsilon_{n}(t+T,\lambda)}$.
Therefore, the time-average of the expectation value is rewritten as

$\displaystyle X_{n}$ $\displaystyle=X_{n}^{\rm inst}+X_{n}^{\rm geom},$ (21)
$\displaystyle X_{n}^{\rm inst}$
$\displaystyle\equiv\frac{1}{T}\int_{0}^{T}dt\bra{\epsilon_{n}(t)}\hat{X}(t)\ket{\epsilon_{n}(t)}$
$\displaystyle=\frac{1}{T}\int_{0}^{T}dt\partial_{\lambda}\epsilon_{n}(t,\lambda)|_{\lambda=0},$
(22) $\displaystyle X_{n}^{\rm geom}$
$\displaystyle\equiv-\frac{1}{T}\partial_{\lambda}\varphi_{n}(\lambda)|_{\lambda=0},$
(23) $\displaystyle\varphi_{n}(\lambda)$
$\displaystyle\equiv\int_{0}^{T}dt\braket{\epsilon_{n}(t,\lambda)}{i\hbar\partial_{t}\epsilon_{n}(t,\lambda)},$
(24)

where $X_{n}^{\rm inst}$ is the time-average of the expectation value of the
band $n$ for the instantaneous eigenvalues, $X_{n}^{\rm geom}$ is the
geometric contribution of the band $n$ due to the adiabatic time evolution,
and $\varphi_{n}(\lambda)$ is the Berry phase associated with the adiabatic
time evolution.

When the parameter $\lambda$ is the external magnetic field $B_{z}$, the time-
averaged expectation values $X$ corresponds to the orbital or spin
magnetization. $X_{n}^{\rm geom}$ in Eq. (23) is the geometrical
magnetization, and $X_{n}^{\rm inst}$ in Eq. (22) is the magnetization from
the snapshot of the eigenstates averaged over one period.

###  III.2 Iteration of unitary transformation

In this section, we explain the method of the adiabatic series expansion
proposed by Berry Berry (1987); Rigolin  _et al._ (2008), and then we show the
expectation value of each state in the second-order adiabatic approximation.

First, we consider the time-periodic Hamiltonian, $\hat{H}(t)$, and its time-
dependent Schrödinger equation, $i\hbar\frac{d}{dt}\psi(t)=\hat{H}(t)\psi(t)$,
where $\psi(t)$ is a wave function. This Hamiltonian is periodic in time with
a period $T=2\pi/\Omega$, $\hat{H}(t+T)=\hat{H}(t)$. Here, we introduce the
rescaled time $\tau=\Omega t$, and rewrite the functions of time as
$\psi(t)=\psi(\tau),\hat{H}(t)=\hat{H}(\tau),\hat{H}(\tau+2\pi)=\hat{H}(\tau)$.
Hence, the time-dependent Schrödinger equation is rewritten as

$i\hbar\Omega\frac{d}{d\tau}\psi(\tau)=\hat{H}(\tau)\psi(\tau).$ (25)

Next, we introduce the time-evolution operator $\hat{U}$, which is unitary. In
terms of the time-evolution operator $\hat{U}(\tau)$, the wave function is
written as $\psi(\tau)=\hat{U}(\tau)\psi(0)$. By substituting it into the
Schrödinger equation, we get

$[\hat{U}^{\dagger}(\tau)\hat{H}(\tau)\hat{U}(\tau)-i\hbar\Omega\hat{U}^{\dagger}(\tau)\partial_{\tau}\hat{U}(\tau)]\psi(0)=0$
(26)

Here, we introduce the instantaneous eigenstate $\ket{\varepsilon_{n}(\tau)}$
with the instantaneous eigenvalue $\varepsilon_{n}(\tau)$,
$\hat{H}(\tau)\ket{\varepsilon_{n}(\tau)}=\varepsilon_{n}(\tau)\ket{\varepsilon_{n}(\tau)}$,
where $n$ is the band index. In the following, we will calculate
$\hat{U}(\tau)$ in the form
$\hat{U}(\tau)=\hat{R}_{0}(\tau)\hat{R}_{1}(\tau)\hat{R}_{2}(\tau)\cdots$,
where $\hat{R}_{i}(\tau)~{}(i=1,2,3,\dots)$ is a unitary operator. As a
zeroth-order approximation of $\hat{U}(\tau)$, we define the unitary operator
$\hat{R}_{0}(\tau)$ as

$\displaystyle\hat{R}_{0}(\tau)=\sum_{n}\ket{\varepsilon_{n}(\tau)}\bra{\varepsilon_{n}(0)}.$
(27)

Then we rewrite Eq. (26),

$\displaystyle\hat{U}^{\dagger}[\hat{H}(\tau)-i\hbar\Omega\partial_{\tau}]\hat{U}$
$\displaystyle=\cdots\hat{R}_{2}^{\dagger}\hat{R}_{1}^{\dagger}[\hat{R}_{0}^{\dagger}\hat{H}(\tau)\hat{R}_{0}-i\hbar\Omega\hat{R}_{0}^{\dagger}\partial_{\tau}\hat{R}_{0}-i\hbar\Omega\partial_{\tau}]\hat{R}_{1}\hat{R}_{2}\cdots$
$\displaystyle=\cdots\hat{R}_{2}^{\dagger}\hat{R}_{1}^{\dagger}[\tilde{H}^{(0)}(\tau)+\Omega\tilde{V}^{(0)}(\tau)-i\hbar\Omega\partial_{\tau}]\hat{R}_{1}\hat{R}_{2}\cdots,$
(28)

and

$\displaystyle\tilde{H}^{(0)}(\tau)$
$\displaystyle=\hat{R}_{0}^{\dagger}\hat{H}(\tau)\hat{R}_{0}$
$\displaystyle=\sum_{n}\varepsilon_{n}(\tau)\ket{\varepsilon_{n}(0)}\bra{\varepsilon_{n}(0)},$
(29) $\displaystyle\tilde{V}^{(0)}(\tau)$
$\displaystyle=\hat{R}_{0}^{\dagger}(-i\hbar\partial_{\tau})\hat{R}_{0}.$ (30)

Note that time-dependence of $\tilde{H}^{(0)}(\tau)$ is contained only in the
eigenvalues $\varepsilon_{n}(\tau)$.

In the first-order approximation of the adiabatic series expansion,
$\Omega\tilde{V}^{(0)}$ is much smaller than $\tilde{H}^{(0)}(\tau)$, and we
drop the $\Omega\tilde{V}^{(0)}$ term. Here, to the first order, we put the
unitary operator $\hat{R}_{2}=\hat{R}_{3}=\dots=\hat{I}$ and calculate
$\hat{R}_{1}$. From Eq. (28), the solution of
$\tilde{H}^{(0)}(\tau)\hat{R}_{1}(\tau)-i\hbar\Omega\partial_{\tau}\hat{R}_{1}(\tau)=0$
is given by

$\hat{R}_{1}(\tau)\equiv\bar{R}_{1}(\tau)=\exp\left(\frac{1}{i\hbar\Omega}\int_{0}^{\tau}ds\tilde{H}^{(0)}(s)\right).$
(31)

In the first-order approximation, the state at $\tau$ is represented as

$\psi(\tau)=\hat{U}(\tau)\psi(0)=\hat{R}_{0}(\tau)\bar{R}_{1}(\tau)\psi(0)$
(32)

Next, we consider the second-order approximation of the adiabatic series
expansion. Here, we put
$\tilde{H}_{1}(\tau)=\tilde{H}^{(0)}(\tau)+\Omega\tilde{V}^{(0)}$, and we set
the instantaneous eigenvector $\ket{\tilde{\varepsilon}_{1,n}(\tau)}$ with the
instantaneous eigenvalues $\tilde{\varepsilon}_{1,n}(\tau)$,
$\tilde{H}_{1}(\tau)\ket{\tilde{\varepsilon}_{1,n}(\tau)}=\tilde{\varepsilon}_{1,n}(\tau)\ket{\tilde{\varepsilon}_{1,n}(\tau)}$.
We perform the perturbation expansion for $\tilde{H}_{1}(\tau),$ and
$\tilde{\varepsilon}_{1,n}(\tau)$ and $\ket{\tilde{\varepsilon}_{1,n}(\tau)}$
are represented to the first order of $\Omega$ as

$\displaystyle\tilde{\varepsilon}_{1,n}(\tau)$
$\displaystyle=\varepsilon_{n}(\tau)+\Omega\bra{\varepsilon_{n}(0)}\tilde{V}^{(0)}(\tau)\ket{\varepsilon_{n}(0)},$
(33) $\displaystyle\ket{\tilde{\varepsilon}_{1,n}(\tau)}$
$\displaystyle=\ket{\varepsilon_{n}(0)}+\Omega\sum_{m(\neq
n)}\ket{\varepsilon_{m}(0)}\frac{\bra{\varepsilon_{m}(0)}\tilde{V}^{(0)}(\tau)\ket{\varepsilon_{n}(0)}}{\varepsilon_{n}(\tau)-\varepsilon_{m}(\tau)}.$
(34)

Note that we assume that the instantaneous eigenvalues $\varepsilon_{n}(\tau)$
for the instantaneous eigenstates $\ket{\varepsilon_{n}(\tau)}$ have no
degeneracy. We define the unitary operator $\hat{R}_{1}(\tau)$ as

$\displaystyle\hat{R}_{1}(\tau)=\sum_{n}\ket{\tilde{\varepsilon}_{1,n}(\tau)}\bra{\tilde{\varepsilon}_{1,n}(0)}.$
(35)

From Eq. (28), one has

$\displaystyle\hat{U}^{\dagger}[H(\tau)-i\hbar\Omega\partial_{\tau}]\hat{U}$
$\displaystyle=\cdots\hat{R}_{2}^{\dagger}[\hat{R}_{1}^{\dagger}\tilde{H}_{1}(\tau)\hat{R}_{1}-i\hbar\Omega\hat{R}_{1}^{\dagger}\partial_{\tau}\hat{R}_{1}-i\hbar\partial_{\tau}]\hat{R}_{2}\cdots$
$\displaystyle=\cdots\hat{R}_{2}^{\dagger}[\tilde{H}^{(1)}(\tau)+\Omega\tilde{V}^{(1)}(\tau)-i\hbar\Omega\partial_{\tau}]\hat{R}_{2}\cdots,$
(36)

and

$\displaystyle\tilde{H}^{(1)}(\tau)$
$\displaystyle=\hat{R}_{1}^{\dagger}\tilde{H}_{1}(\tau)\hat{R}_{1}$
$\displaystyle=\sum_{n}\tilde{\varepsilon}_{1,n}(\tau)\ket{\tilde{\varepsilon}_{1,n}(0)}\bra{\tilde{\varepsilon}_{1,n}(0)},$
(37) $\displaystyle\tilde{V}^{(1)}(\tau)$
$\displaystyle=\hat{R}_{1}^{\dagger}(-i\hbar\partial_{\tau})\hat{R}_{1}.$ (38)

Note that time-dependence of $\tilde{H}^{(1)}(\tau)$ is contained only in the
eigenvalues $\tilde{\varepsilon}_{1,n}(\tau)$.

In the second-order approximation of the adiabatic series expansion, we assume
that $\Omega\tilde{V}^{(1)}(\tau)$ is of the order $O(\Omega^{2})$ and this
term is sufficiently smaller than $\tilde{H}^{(1)}(\tau)$. Hence, we drop the
$\Omega\tilde{V}^{(1)}(\tau)$ term. Therefore, to the second order, we put the
unitary operators $\hat{R}_{3}=\hat{R}_{4}=\dots=\hat{I}$ and calculate
$\hat{R}_{2}$. From Eq. (36), the solution of
$\tilde{H}^{(1)}(\tau)\hat{R}_{2}(\tau)-i\hbar\Omega\partial_{\tau}\hat{R}_{2}(\tau)=0$
is given by

$\hat{R}_{2}(\tau)\equiv\bar{R}_{2}(\tau)=\exp\left(\frac{1}{i\hbar\Omega}\int_{0}^{\tau}ds\tilde{H}^{(1)}(s)\right).$
(39)

In the second-order approximation, the state at time $\tau$ is represented as

$\psi(\tau)=\hat{U}(\tau)\psi(0)=\hat{R}_{0}(\tau)\hat{R}_{1}(\tau)\bar{R}_{2}(\tau)\psi(0).$
(40)

Similarly, we can iterate this expansion $n$ times, and then the state at
$\tau$ is given as

$\psi(\tau)=\hat{R}_{0}(\tau)\hat{R}_{1}(\tau)\cdots\bar{R}_{n}(\tau)\psi(0).$
(41)

Then, the expectation value of an operator $\hat{X}$ in the second order is
given by

$\displaystyle\braket{X(\tau)}_{n}$
$\displaystyle=\bra{\psi_{n}(\tau)}\hat{X}\ket{\psi_{n}(\tau)}$
$\displaystyle=\bra{\psi_{n}(0)}\bar{R}_{2}^{\dagger}\hat{R}_{1}^{\dagger}\hat{R}_{0}^{\dagger}\hat{X}\hat{R}_{0}\hat{R}_{1}\bar{R}_{2}\ket{\psi_{n}(0)},$
(42)

where $\ket{\psi_{n}(0)}$ is the initial state. First, the unitary operator
$\bar{R}_{2}(\tau)$ is rewritten as

$\displaystyle\bar{R}_{2}(\tau)=$
$\displaystyle\exp\left(\frac{-i}{\hbar\Omega}\int_{0}^{\tau}ds\sum_{n}\tilde{\varepsilon}_{1,n}(s)\ket{\tilde{\varepsilon}_{1,n}(0)}\bra{\tilde{\varepsilon}_{1,n}(0)}\right)$
$\displaystyle=\sum_{n}e^{\frac{-i}{\hbar\Omega}\int_{0}^{\tau}ds\tilde{\varepsilon}_{1,n}(s)}\ket{\tilde{\varepsilon}_{1,n}(0)}\bra{\tilde{\varepsilon}_{1,n}(0)},$
(43)

from the orthogonality relation
$\braket{\tilde{\varepsilon}_{1,n}(\tau)}{\tilde{\varepsilon}_{1,m}(\tau)}=\delta_{n,m}$.
By substituting Eqs. (27), (35), (43) into Eq. (42), the expectation value is
rewritten as

$\displaystyle\braket{X(\tau)}_{n}$
$\displaystyle=X_{nn}(\tau)+\Omega\sum_{m(\neq
n)}\left[\left(1-e^{\frac{i}{\hbar\Omega}\int_{0}^{\tau}ds[\tilde{\varepsilon}_{1,n}(s)-\tilde{\varepsilon}_{1,m}(s)]}\right)X_{nm}(\tau)\frac{\bra{\varepsilon_{m}(0)}\tilde{V}^{(0)}(\tau)\ket{\varepsilon_{n}(0)}}{\varepsilon_{n}(\tau)-\varepsilon_{m}(\tau)}+{\rm
c.c.}\right]$ $\displaystyle=X_{nn}(\tau)+\hbar\Omega\sum_{m(\neq
n)}\left[\left(1-e^{\frac{i}{\hbar\Omega}\int_{0}^{\tau}ds[\varepsilon_{n}(s)-\varepsilon_{m}(s)-\hbar\Omega(\gamma_{n}(s)-\gamma_{m}(s))]}\right)X_{nm}(\tau)\frac{\braket{\varepsilon_{m}(\tau)}{(-i\partial_{\tau})\varepsilon_{n}(\tau)}}{\varepsilon_{n}(\tau)-\varepsilon_{m}(\tau)}+{\rm
c.c.}\right],$ (44)

where
$X_{nm}(\tau)=\bra{\varepsilon_{n}(\tau)}\hat{X}\ket{\varepsilon_{m}(\tau)}$,
and
$\gamma_{n}(\tau)=\braket{\varepsilon_{n}(\tau)}{(-i\partial_{\tau})\varepsilon_{n}(\tau)}$
is Berry connection. The first term shows the instantaneous expectation value
at time $\tau$, and the second term is proportional to the inverse time period
$1/T$. The time-averaged expectation value $\bar{X}_{n}$ is given as

$\displaystyle\bar{X}_{n}$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}d\tau\braket{X(\tau)}_{n}$
$\displaystyle\simeq\frac{1}{2\pi}\int_{0}^{2\pi}d\tau\Bigg{[}X_{nn}(\tau)+\hbar\Omega\sum_{m(\neq
n)}X_{nm}(\tau)\frac{\bra{\varepsilon_{m}(\tau)}(-i\partial_{\tau})\ket{\varepsilon_{n}(\tau)}}{\varepsilon_{n}(\tau)-\varepsilon_{m}(\tau)}+{\rm
c.c.}\Bigg{]}.$ (45)

Since we have assumed $\tilde{\varepsilon}_{1,n}/\hbar\Omega\gg 1$, the phase
term in Eq. (44) rapidly oscillates, and this phase term has been neglected
after time averaging. This time-averaged expectation value is equal to Eqs.
(22) and (23) in the Berry-phase method in adiabatic process.

##  IV time-averaged spin magnetization

In this section, we calculate the time-averaged spin magnetization for our toy
model. The spin operator is expressed as

$\displaystyle\hat{S}_{\alpha}=\frac{\hbar}{2}\sigma_{0}\otimes
s_{\alpha}\equiv\frac{\hbar}{2}\hat{s}_{\alpha}~{}~{}(\alpha=x,y,z).$ (46)

From Eq. (15), the spin magnetization is represented as

$\displaystyle\hat{H}_{\rm B}=-\hat{\bm{m}}_{\rm
spin}\cdot\bm{B}=\frac{g\mu_{B}}{\hbar}\hat{\bm{S}}\cdot\bm{B}=\mu_{B}\bm{\hat{s}}\cdot\bm{B},$
(47)

with g-factor $g=2$ and the Bohr magneton $\mu_{\rm B}$. First, we calculate
the spin expectation values without phonons. We show the electronic energy
band structure for the Hamiltonian without phonons Eq. (4) in Fig. 2. In Fig
2, we assume that the Fermi energy $\varepsilon_{F}$ is zero, and the
temperature is zero, and then this toy model is a band insulator.

Figure 2:  Electronic energy band structure along the high-symmetry line in
the system without phonons. We set the parameter values to be $\lambda_{\rm
v}=0.2t_{0}$, and $\lambda_{\rm R}=0.4t_{0}$.

The staggered potential term breaks inversion symmetry, and the electron spin
is split by the Rashba spin-orbit interaction term. In the following we
consider the spin expectation values below the Fermi energy
$\varepsilon_{F}(=0)$.

We introduce the $\Theta$ operator $\Theta=i(\sigma_{0}\otimes s_{y})K$, where
$K$ is complex conjugate operator, and this corresponds to time-reversal
operator in static systems without phonons. The relations between the $\Theta$
operator and $\Gamma$ matrices are given by
$\Theta\Gamma^{a}\Theta^{-1}=\Gamma^{a}$,
$\Theta\Gamma^{ab}\Theta^{-1}=-\Gamma^{ab}$ with $a,b=1,2,3,4,5$. From Eq.
(4), the Hamiltonian $H_{0}(\bm{k})$ and its eigenstate
$\ket{\varepsilon_{n}(\bm{k})}$ with the eigenenergy $\varepsilon_{n}(\bm{k})$
is represented as $\Theta H_{0}(\bm{k})\Theta^{-1}=H_{0}(-\bm{k})$ and
$\Theta\ket{\varepsilon_{n}(\bm{k})}=\ket{\varepsilon_{n}(-\bm{k})}$. Since
the spin operator $\hat{S}_{\alpha}$ satisfies
$\Theta\hat{S}_{\alpha}\Theta^{-1}=-\hat{S}_{\alpha}$, the spin expectation
values $\braket{\hat{\bm{S}}(\bm{k})}_{n}$ is an odd function of the
wavevector $\bm{k}$ in our toy model without phonons, and the total spin
expectation values is zero due to cancellation between $\bm{k}$ and $-\bm{k}$.

Next, we show the spin expectation values with the microscopic local rotation
of atoms. By applying Eqs. (21), (22), and (23) to our toy model, the time-
averaged spin magnetization are given by

$\displaystyle\bm{M}_{\rm spin}$
$\displaystyle=\frac{1}{VT}\int_{0}^{T}dt\sum_{\bm{k}}\sum_{n}^{\rm
occ}\Bigg{[}\bra{\varepsilon_{n}(\bm{k},t)}\partial_{\bm{B}}\hat{H}_{\rm
B}|_{\bm{B}=0}\ket{\varepsilon_{n}(\bm{k},t)}$ $\displaystyle\hskip
56.9055pt+\sum_{m(\neq
n)}\frac{\bra{\varepsilon_{n}(\bm{k},t)}\partial_{\bm{B}}\hat{H}_{\rm
B}|_{\bm{B}=0}\ket{\varepsilon_{m}(\bm{k},t)}\braket{\varepsilon_{m}(\bm{k},t)}{(-i\hbar\partial_{t})\varepsilon_{n}(\bm{k},t)}}{\varepsilon_{n}(\bm{k},t)-\varepsilon_{m}(\bm{k},t)}+{\rm
c.c.}\Bigg{]}$ $\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}d\tau\sum_{n}^{\rm
occ}[\bm{m}_{n}^{\rm inst}(\tau)+\bm{m}_{n}^{\rm geom}(\tau)],$ (48)

and

$\displaystyle\bm{m}_{n}^{\rm inst}(\tau)$ $\displaystyle=\frac{\mu_{\rm
B}}{V}\sum_{\bm{k}}\hat{\bm{s}}_{nn}(\bm{k},\tau),$ (49)
$\displaystyle\bm{m}_{n}^{\rm geom}(\tau)$ $\displaystyle=\frac{\mu_{\rm
B}}{V}\sum_{\bm{k}}\sum_{m(\neq
n)}\Bigg{[}\frac{\hbar\Omega\hat{\bm{s}}_{nm}(\bm{k},\tau)A_{mn}(\bm{k},\tau)}{\varepsilon_{n}(\bm{k},\tau)-\varepsilon_{m}(\bm{k},\tau)}+{\rm
c.c.}\Bigg{]}$ (50)

where $\ket{\varepsilon_{n}(\bm{k},t)}$ is the instantaneous eigenstate with
instantaneous eigenenergy $\varepsilon_{n}(\bm{k},t)$ for the snapshot
Hamiltonian in momentum space $\hat{H}(\bm{k},t)$,
$\hat{\bm{s}}_{nm}(\bm{k},\tau)=\bra{\varepsilon_{n}(\bm{k},\tau)}\hat{\bm{s}}\ket{\varepsilon_{m}(\bm{k},\tau)}$
is the instantaneous matrix element for spin between bands $n$ and $m$ at time
$\tau=\Omega t$, and
$A_{mn}(\bm{k},\tau)=-i\braket{\varepsilon_{m}(\bm{k},\tau)}{\partial_{\tau}\varepsilon_{n}(\bm{k},\tau)}$
is the Berry connection defined by the instantaneous eigenstates between bands
$n$ and $m$. Here, the summation over $n$ represents that over the occupied
bands. $\bm{m}_{n}^{\rm inst}(\tau)$ is the instantaneous spin magnetization
for the energy band $n$ at time $\tau$, and $\bm{m}_{n}^{\rm geom}(\tau)$ is
the geometrical spin magnetization in the adiabatic process at time $\tau$ and
this term is proportional to the phonon frequency $\Omega$. We note that while
$\hat{s}_{nm}$ and $\hat{A}_{nm}~{}(m\neq n)$ depend on a gauge choice,
$m_{n}^{\rm geom}(\tau)$ is gauge invariant.

Here, we apply the $\Theta$ operator to our toy model with phonons. From Eq.
(13), the time-dependent Hamiltonian is rewritten as

$\displaystyle
H(\bm{k},\tau)=H_{0}(\bm{k})+H_{c}(\bm{k})\cos\tau+H_{s}(\bm{k})\sin\tau.$
(51)

Therefore, the time-dependent Hamiltonian satisfies

$\displaystyle\Theta H(\bm{k},\tau)\Theta^{-1}=H(-\bm{k},\tau),$ (52)

because $\Theta H_{c(s)}(\bm{k})\Theta^{-1}=H_{c(s)}(-\bm{k})$. Thus, within
our dynamical model, the operator $\Theta$ does not flip the sign of $\tau$,
and therefore, to avoid confusion we call it a $\Theta$ operation instead of
time-reversal operation. Hence the instantaneous eigenstates
$\ket{\varepsilon_{n}(\bm{k},\tau)}$ and the instantaneous eigenenergy
$\varepsilon_{n}(\bm{k},\tau)$ satisfy

$\displaystyle\Theta\ket{\varepsilon_{n}(\bm{k},\tau)}=\ket{\varepsilon_{n}(-\bm{k},\tau)},~{}~{}\varepsilon_{n}(\bm{k},\tau)=\varepsilon_{n}(-\bm{k},\tau).$
(53)

From these relations, the instantaneous matrix elements for spin
$\hat{\bm{s}}_{nm}(\bm{k},\tau)$ and the Berry connection
$A_{nm}(\bm{k},\tau)$ are given by

$\displaystyle\hat{\bm{s}}_{nm}(\bm{k},\tau)$
$\displaystyle=\bra{\varepsilon_{n}(\bm{k},\tau)}\hat{\bm{s}}\ket{\varepsilon_{m}(\bm{k},\tau)}$
$\displaystyle=\bra{\varepsilon_{n}(-\bm{k},\tau)}(-\hat{\bm{s}})\ket{\varepsilon_{m}(-\bm{k},\tau)}$
$\displaystyle=-\hat{\bm{s}}_{nm}(-\bm{k},\tau),$ (54)

and

$\displaystyle A_{nm}(\bm{k},\tau)$
$\displaystyle=\braket{\varepsilon_{n}(\bm{k},\tau)}{(-i\partial_{\tau})\varepsilon_{m}(\bm{k},\tau)}$
$\displaystyle=\braket{\varepsilon_{n}(-\bm{k},\tau)}{(i\partial_{\tau})\varepsilon_{m}(-\bm{k},\tau)}$
$\displaystyle=-A_{nm}(-\bm{k},\tau).$ (55)

By substituting these relations into Eq. (48), the instantaneous spin
magnetization $\bm{m}_{n}^{\rm inst}(\tau)$ becomes zero since its sum over
the first Brillouin zone cancels between $\bm{k}$ and $-\bm{k}$ at time
$\tau$. On the other hand, the geometrical spin magnetization $\bm{m}_{n}^{\rm
geom}(\tau)$ can be nonzero since it is an even function of the wave vector
$\bm{k}$ at time $\tau$. Therefore, only the geometrical spin magnetization
contributes to the time-averaged spin magnetization.

We calculate the geometrical spin magnetization $\bm{m}_{\alpha,n}^{\rm geom}$
for each band below the Fermi energy $\varepsilon_{F}(=0)$, and it is nonzero
as shown in Figs. 3(a) - (c).

Figure 3:  Geometrical spin magnetization for one cycle in the system with
microscopic rotation versus the rescaled time $\tau$ for (a) $x$, (b) $y$, and
(c) $z$ components. (d) Evolution of spin magnetization for one cycle in the
$xy$-plane. The small circles show the values at $\tau=0$ and the arrows show
the time-evolution direction. We set the parameters as $\lambda_{\rm
v}=0.2t_{0},\lambda_{\rm R}=0.4t_{0},\delta t=0.1t_{0}$, and
$\delta\lambda_{\rm R}=0.1\lambda_{\rm R}$. In (a) - (d), the blue line shows
the lowest band ($n=1$), and the red line shows the second lowest band
($n=2$), and the black line shows the sum of the bands below the Fermi energy.

In this case, the generators of the space-time symmetries are
$\\{C_{3z}|T/3\\}$ and $\\{m_{\tau}m_{x}\Theta|T/2\\}$, where
$\\{C_{3z}|T/3\\}$ is the threefold rotation around the $z$ axis, followed by
$T/3$ time translation, and $\\{m_{\tau}m_{x}\Theta|T/2\\}$ is a product of
the $\Theta$ operation, the mirror reflection with respect to the $yz$-plane,
and $m_{\tau}$ operation, which represents inversion of the sign of time
$\tau$, followed by $T/2$ time translation. We note that the dynamical
symmetry is defined using the discrete time-translation invariance Floquet
systems Morimoto  _et al._ (2017); Xu and Wu (2018). These symmetries
$\\{C_{3z}|T/3\\}$ and $\\{m_{\tau}m_{x}\Theta|T/2\\}$ are represented by
unitary operators $U_{1}=\exp(-i\frac{\pi}{3}(\sigma_{0}\otimes s_{z}))$ and
$U_{2}=(\sigma_{0}\otimes s_{x})\Theta$, respectively. Here, the time-
dependent Hamiltonian satisfies

$\displaystyle
U_{1}\hat{H}(\bm{k},\tau)U_{1}^{-1}=\hat{H}\left(R\bm{k},\tau+\frac{2\pi}{3}\right),$
(56) $\displaystyle
U_{2}\hat{H}(k_{x},k_{y},\tau)U_{2}^{-1}=\hat{H}(k_{x},-k_{y},-\tau+\pi),$
(57)

where $R$ is the three-fold rotational matrix around the $z$-axis,
$R\bm{k}=(-k_{x}/2-\sqrt{3}k_{y}/2,\sqrt{3}k_{y}/2-k_{x}/2)$. From these
relations, the instantaneous matrix elements for spin
$\hat{\bm{s}}_{nm}(\bm{k},\tau)$ and the Berry connection
$A_{nm}(\bm{k},\tau)$ are given by

$\displaystyle\hat{\bm{s}}_{nm}(\bm{k},\tau)$
$\displaystyle=R^{-1}\hat{\bm{s}}_{nm}\left(R\bm{k},\tau+\frac{2\pi}{3}\right),$
(58) $\displaystyle A_{nm}(\bm{k},\tau)$
$\displaystyle=A_{nm}\left(R\bm{k},\tau+\frac{2\pi}{3}\right),$ (59)

and

$\displaystyle\begin{pmatrix}\hat{s}_{x,nm}(k_{x},k_{y},\tau)\\\
\hat{s}_{y,nm}(k_{x},k_{y},\tau)\\\
\hat{s}_{z,nm}(k_{x},k_{y},\tau)\end{pmatrix}$
$\displaystyle=\begin{pmatrix}-\hat{s}_{x,nm}(k_{x},-k_{y},-\tau+\pi)\\\
\hat{s}_{y,nm}(k_{x},-k_{y},-\tau+\pi)\\\
\hat{s}_{z,nm}(k_{x},-k_{y},-\tau+\pi)\end{pmatrix},$ (60) $\displaystyle
A_{nm}(k_{x},k_{y},\tau)$ $\displaystyle=A_{nm}(k_{x},-k_{y},-\tau+\pi).$ (61)

Hence, the geometrical spin magnetization satisfies

$\displaystyle\bm{m}_{n}^{\rm geom}(\tau)$
$\displaystyle=R^{-1}\bm{m}_{n}^{\rm geom}\left(\tau+\frac{2\pi}{3}\right),$
(62) $\displaystyle m_{x,n}^{\rm geom}(\tau)$ $\displaystyle=-m_{x,n}^{\rm
geom}(-\tau+\pi),$ (63) $\displaystyle m_{y(z),n}^{\rm geom}(\tau)$
$\displaystyle=m_{y(z),n}^{\rm geom}(-\tau+\pi).$ (64)

The geometrical spin magnetization of the $z$ components oscillates with a
$T/3$ period. On the other hand, the geometrical spin magnetization in the
$xy$-plane rotates, with its trajectory forming a triangle-like form for one
cycle as shown in Fig 3(d). From Eqs. (63), (64), at $\tau=\pi/2,3\pi/2$, the
geometrical spin magnetization for each band lie in the $yz$-plane. Moreover,
from the $\\{C_{3z}|T/3\\}$ symmetry, the geometrical spin magnetization in
the $xy$-plane rotates once within the period $T$, and its time average over
one period vanishes. In addition, when the direction of the microscopic local
rotation is reversed, the geometrical spin magnetizations are reversed.
Therefore, the microscopic local rotation works like an effective magnetic
field, and the spin magnetization along the rotational axis is proportional to
the phonon frequency.

##  V Discussion

In this section, we discuss how the microscopic local rotation couples with
electron spins. First, when the Rashba parameter $\lambda_{\rm R}$ is zero in
our toy model, the time-averaged spin magnetization is zero since the energy
bands are spin degenerate and the instantaneous matrix element for spin
$\hat{\bm{s}}_{nm}(\bm{k},\tau)$ becomes zero for arbitrary wave vector, time,
and bands. Therefore, within our model, the spin-orbit coupling is required
for spin polarization induced by microscopic atomic rotations.

Next we compare two cases: a simple vibration of atoms and a microscopic local
rotation of atoms. For this comparison, we consider a case where the atoms A
and B mutually vibrate with the displacement vector $\bm{u}$ given by

$\displaystyle\bm{u}=\bm{u}_{\rm B}-\bm{u}_{\rm A}=-u_{+}(\cos\tau,0).$ (65)

Then the time-dependent Hamiltonian Eq. (13) is rewritten as
$H(\bm{k},\tau)=H_{0}(\bm{k})+H_{c}(\bm{k})\cos\tau$. We calculate the
geometric spin magnetization in this system with the simple vibration of
atoms, and it is shown in Fig. 4.

Figure 4:  Geometrical spin magnetization for one cycle in the system with the
simple vibration versus rescaled time $\tau$ for (a) $x$, (b) $y$, and (c) $z$
components. (d) Evolution of spin magnetization for one cycle in the
$xy$-plane. The small circles show the values at $\tau=0$ and the arrows show
the time-evolution direction. We set the parameters as $\lambda_{\rm
v}=0.2t,\lambda_{\rm R}=0.4t_{0},\delta t=0.1t_{0},$ and $\delta\lambda_{\rm
R}=0.1\lambda_{\rm R}$. In (a) - (d), the blue lines show the lowest band
$(n=1)$, and the red lines show the second lowest band $(n=2)$, and the black
lines show the sum of the below the Fermi energy.

In Figs. 4(a), (b), and (c), we show the $x,y$ and $z$ componets of
geometrical spin magnetization versus time $\tau$, respectively. In Fig. 4(d),
we show the evolution of the geometrical spin magnetization for one period in
the $xy$-plane. As is similar to the results with the microscopic local
rotation of atoms, the geometrical spin magnetization can be finite at time
$\tau$. Meanwhile the time-averaged spin magnetization vanishes.

We discuss the symmetry reason why the time-averaged spin magnetization
vanishes. In our toy model with a simple vibration of atoms, the Hamiltonian
and instantaneous eigenstates satisfy not only Eqs. (52) and (53) but also the
following relations,

$\displaystyle\Theta H(\bm{k},\tau)\Theta^{-1}$
$\displaystyle=H(-\bm{k},-\tau),$ (66)
$\displaystyle\Theta\ket{\varepsilon_{n}(\bm{k},\tau)}$
$\displaystyle=\ket{\varepsilon_{n}(-\bm{k},-\tau)},$ (67)

because $H(\bm{k},\tau)=H(\bm{k},-\tau)$. The instantaneous matrix element for
spin $\hat{\bm{s}}_{nm}(\bm{k},\tau)$ and the Berry connection
$A_{nm}(\bm{k},\tau)$ satisfies

$\displaystyle\hat{\bm{s}}_{nm}(\bm{k},\tau)$
$\displaystyle=-\hat{\bm{s}}_{nm}(-\bm{k},-\tau),$ (68) $\displaystyle
A_{nm}(\bm{k},\tau)$ $\displaystyle=A_{nm}(-\bm{k},-\tau).$ (69)

From Eqs. (54), (55), (68), and (69), the geometrical spin magnetization is an
even function of the wave vector $\bm{k}$ and at the same time an odd function
of time $\tau$, and therefore the geometrical spin magnetization has opposite
signs between $\tau$ and $-\tau$. In this case, the generators of the space-
time symmetry of the system are $\\{\Theta\\}$ and $\\{m_{x}|T/2\\}$, which is
the mirror reflection with respect to the $yz$-plane $m_{x}$ and the time
translation by half the period $T/2$. This symmetry $\\{m_{x}|T/2\\}$ is
represented by the unitary operator $U_{3}=\sigma_{0}\otimes s_{x}$. From
$U_{3}\hat{H}(k_{x},k_{y},\tau)U_{3}^{-1}=\hat{H}(-k_{x},k_{y},\tau+\pi)$, the
instantaneous matrix element for spin and the Berry connection satisfy the
following relations:

$\displaystyle\hat{s}_{x,nm}(k_{x},k_{y},\tau)$
$\displaystyle=\hat{s}_{x,nm}(-k_{x},k_{y},\tau+\pi),$ (70)
$\displaystyle\hat{s}_{y(z),nm}(k_{x},k_{y},\tau)$
$\displaystyle=-\hat{s}_{y(z),nm}(-k_{x},k_{y},\tau+\pi),$ (71) $\displaystyle
A_{nm}(k_{x},k_{y},\tau)$ $\displaystyle=A_{nm}(-k_{x},k_{y},\tau+\pi).$ (72)

By combining these relations for the $\Theta$ operator (Eqs. (68) and (69))
and symmetry operator $\\{m_{x}|T/2\\}$ (Eqs. (70) - (72)), the geometrical
spin magnetization $m_{\alpha,n}^{\rm geom}(\tau)$ satisfies

$\displaystyle m_{x,n}^{\rm geom}(\tau)$ $\displaystyle=-m_{x,n}^{\rm
geom}(-\tau)=-m_{x,n}^{\rm geom}(-\tau+\pi),$ (73) $\displaystyle
m_{y(z),n}^{\rm geom}(\tau)$ $\displaystyle=-m_{y(z),n}^{\rm
geom}(-\tau)=m_{y(z),n}^{\rm geom}(-\tau+\pi).$ (74)

At $\tau=0,\pi$, the geometrical spin magnetization becomes zero, and at
$\tau=\pi/2,3\pi/2$, the $x$ component of the geometrical spin magnetization
becomes zero, while its the $y$ and $z$ components are symmetric with respect
to $\tau=\pi/2,3\pi/2$. This symmetry analysis coincides the numerical results
shown in Fig. 4. Therefore, the time-averaged spin magnetization vanishes for
a simple vibration of atoms by time-reversal symmetry.

Next, we discuss an effect of the modulation of the Rashba term
$\delta\lambda_{\rm R}$.

Figure 5:  Geometrical spin magnetization for one cycle in the system with
microscopic rotation as $\delta\lambda_{\rm R}=0$ versus the rescaled time
$\tau$ for (a)$x$, (b)$y$, and (c)$z$ components. (d) Evolution of spin
magnetization for one cycle in the $xy$-plane. The small circles show the
$\tau=0$ and the arrows show the time-evolution direction. We set the
parameters as $\lambda_{\rm v}=0.2t_{0},\lambda_{\rm R}=0.4t_{0},\delta
t=0.1t_{0}$, and $\delta\lambda_{\rm R}=0$. In (a) - (d), the blue lines show
the lowest band ($n=1$), and the red lines show the second lowest band
($n=2$), and the black lines show the sum of the bands below the Fermi energy.

In our model, even without the modulation of the Rashba term
$\delta\lambda_{\rm R}=0$, the time-averaged spin magnetization along the
microscopic rotational axis becomes nonzero as shown in Fig. 5(c). Figures
5(a), (b), and (d) show that the geometrical spin magnetization in the
$xy$-plane isotropically rotates for one cycle. When the modulation of the
Rashba term is absent, the electron spin couples with the microscopic local
rotation indirectly via the Rashba spin-orbit term, and the magnitude of the
geometrical spin magnetization becomes smaller than that with the modulation
of the Rashba term (Fig. 3). Therefore, the microscopic local rotation of
atoms couples with the electron spins via the Rashba spin-orbit interaction,
and then the time-averaged spin magnetization along the rotational axis is
generated.

Next, we discuss how the spin magnetization depends on the gap. The
geometrical spin magnetization can be rewritten as

$\displaystyle\sum_{n}\bm{m}_{n}^{\rm geom}(\tau)$
$\displaystyle=\hbar\Omega\sum_{\bm{k}}\sum_{n}^{\rm occ}\sum_{m(\neq
n)}\left[\frac{\hat{\bm{s}}_{nm}(\bm{k},\tau)A_{mn}(\bm{k},\tau)}{\varepsilon_{n}(\bm{k},\tau)-\varepsilon_{m}(\bm{k},\tau)}+{\rm
c.c.}\right]$ $\displaystyle=\hbar\Omega\sum_{\bm{k}}\sum_{n}^{\rm
occ}\sum_{m}^{\rm
unocc}\left[\frac{\hat{\bm{s}}_{nm}(\bm{k},\tau)A_{mn}(\bm{k},\tau)}{\varepsilon_{n}(\bm{k},\tau)-\varepsilon_{m}(\bm{k},\tau)}+{\rm
c.c.}\right].$ (75)

In the second line, we divide the sum over all the bands $m$ to that over the
occupied bands and that over the unoccupied bands, and the former sum
vanishes. From Eq. (75), since the geometrical spin magnetization is
proportional to the inverse of the difference of energy
$[\epsilon_{n}(\bm{k},\tau)-\epsilon_{m}(\bm{k},\tau)]^{-1}$, the magnitude of
the geometrical spin magnetization becomes larger when the gap is small,
namely when the staggered potential $\lambda_{\rm v}$ is small, as shown in
Fig. 6(a). In Figs 6 (b), (c), (d) and (e), we show the the time-averaged
energy expectation value $\bar{E}$ and the time-averaged geometrical spin
magnetization $m_{n,z}^{\rm geom}$ on the high-symmetry line, where
$\bar{E}_{n}(\bm{k})=\frac{1}{2\pi}\int_{0}^{2\pi}d\tau\bra{\varepsilon_{n}(\bm{k},\tau)}\hat{H}(\bm{k},\tau)\ket{\varepsilon_{n}(\bm{k},\tau)}$
is the time-averaged energy expectation value and the time-averaged
geometrical spin magnetization $m_{z,n}^{\rm geom}(\bm{k})$ is given from Eq.
(48) as

$\displaystyle\bm{m}_{n}^{\rm geom}(\bm{k})$ $\displaystyle=\frac{\mu_{\rm
B}\hbar\Omega}{2\pi}\int_{0}^{2\pi}d\tau\sum_{m(\neq
n)}\left[\frac{\bm{\hat{s}}_{nm}(\bm{k},\tau)A_{mn}(\bm{k},\tau)}{\varepsilon_{n}(\bm{k},\tau)-\varepsilon_{m}(\bm{k},\tau)}+{\rm
c.c.}\right].$ (76)

Moreover, we show the dependence on the modulated hopping and the modulated
Rashba parameters of the geometrical spin magnetization in Fig 7. From Figs.
7(a) and (b), the geometrical spin magnetization is proportional to the square
of the modulated parameters $\delta t$ and $\delta\lambda_{\rm R}$ in the
lowest order. Thus the physical picture of the present effect as follows. Our
model is insulating in the static case, and the phonons dynamically excite
electrons across the gap. The excitation depends on spin via the Rashba term,
and as a result, the system acquires a spin magnetization. Figures 6 (b), (c),
(d) and (e) show the distribution of the contribution to the spin
magnetization in Eq. (76) as a color map. When the staggered potential
$\lambda_{\rm v}$ is small (Figs. 6(b) and (c)), the gap is small, and the
time-averaged geometrical spin magnetization around $K$ and $K^{\prime}$
points for $n=2$ mainly contributes to the spin magnetization. On the other
hand, when the staggered potential $\lambda_{\rm v}$ becomes larger (Figs.
6(d) and (e)), the time-averaged geometrical spin magnetization around $K$ and
$K^{\prime}$ points cancels between $n=1$ and $n=2$ bands, and then the spin
magnetization becomes small. Moreover, since the geometrical spin
magnetization is of the second order of the displacement by phonons, the spin
polarization is generated by interactions between the dynamically excited
states for electrons. This clearly shows that the dynamically excited
electrons across the gaps around $K$ and $K^{\prime}$ points generate spin
polarization. Therefore, as the band gap becomes larger, the geometrical spin
magnetization becomes smaller.

Figure 6:  Dependence on the staggered potential $\lambda_{\rm v}$ of
geometrical spin magnetization. (a) Geometrical spin magnetization with or
without modulated Rashba parameter $\delta\lambda_{\rm R}$ versus staggered
potential $\lambda_{\rm v}$. The blue lines show the geometrical spin
magnetization without modulated Rashba parameter $\delta\lambda_{\rm R}=0$,
and the red line show it with $\delta\lambda_{\rm R}=0.04t$. (b) - (e): Time-
averaged energy expectation values and the geometrical spin magnetization. For
(b), (c), (d) and (e), the colormap shows the time-averaged geometrical spin
magnetization on the high-symmetry line. We set the parameters as
$\lambda_{\rm R}=0.4t_{0},\delta t=0.1t_{0}$ and (b) $\lambda_{\rm
v}=0.2t_{0},\delta\lambda_{\rm R}=0$, (c) $\lambda_{\rm
v}=0.2t_{0},\delta\lambda_{\rm R}=0.04t_{0}$, (d) $\lambda_{\rm
v}=t_{0},\delta\lambda_{\rm R}=0$, and (e) $\lambda_{\rm
v}=t_{0},\delta\lambda_{\rm R}=0.04t_{0}$.  Figure 7:  Dependence on the
modulated hopping parameter $\delta t$ and the modulated Rashba parameter
$\delta\lambda_{\rm R}$ of the geometrical spin magnetization. (a) Geometrical
spin magnetization without the modulated Rashba parameter $\delta\lambda_{\rm
R}=0$ versus the modulated hopping parameter $\delta t^{2}$. (b) Geometrical
spin magnetization without the modulated hopping parameter $\delta t=0$ versus
the modulated Rashba parameter $\delta\lambda_{\rm R}^{2}$. We set the
parameters as $\lambda_{\rm v}=0.2t_{0},\lambda_{\rm R}=0.4t_{0}$. The blue
lines show the lowest band $(n=1)$, and the red lines show the second lowest
band $(n=2)$, and the black lines show the sum of the bands below the Fermi
energy.

##  VI Conclusion

In summary, we have theoretically show that a microscopic local rotation of
atoms induces electron spins in the system with Rashba spin-orbit interaction.
Using our toy model which is a two-dimensional tight-binding model on the
honeycomb lattice with Rashba spin-orbit interaction and phonons, we calculate
the time-averaged spin magnetization using the adiabatic approximation, which
is valid because the phonon frequency is much smaller than the frequency scale
from the electron bandwidth. The time-averaged spin magnetization have two
terms; the instantaneous spin magnetization and the geometrical spin
magnetization. The instantaneous spin magnetization is given by the
instantaneous eigenstates for the snapshot Hamiltonian, and this term vanishes
in the snapshot system with time-reversal symmetry. The geometrical spin
magnetization is represented in terms of the product of the Berry connection
and the instantaneous matrix element for spin, and it is proportional to the
phonon frequency $\Omega$. When the system has both the Rashba spin-orbit
interaction and the microscopic local rotation of atoms, the geometrical spin
magnetization becomes nonzero. The geometrical spin magnetization for each
band perpendicular to the rotational axis have a finite value as a function of
time, and its time average over one period is zero. Meanwhile the time-
averaged spin magnetization along the rotational axis has a finite value and
it is proportional to the phonon frequency $\Omega$. On the other hand, when
the motions of atoms are vibrations, the geometrical spin magnetization for
each band has a finite value as a function of time, but the time-averaged spin
magnetization vanishes due to the time-reversal symmetry. Therefore, we have
shown that the electron spins couple to the microscopic local rotation of
atoms via the spin-orbit interaction, and the time-averaged spin magnetization
along the rotational axis is proportional to the phonon frequency.

In the present paper, we studied phonons at the $\Gamma$ point. We expect that
phonons at general $\bm{k}$ points leads to a qualitatively similar result.
This is because the effect which we discuss is of the second order in the
phonon displacement, and so even for phonons with wavevector $\bm{k}$ can give
a uniform magnetization in this second order due to $\bm{k}+(-\bm{k})=0$.
Detailed behaviors of the spin magnetization sensitively depend on systems and
phonons considered, and are left as futrue problems.

The generation of spin polarization due to rotational motion looks similar to
the effect due to spin-rotation coupling discussed in Refs. Matsuo  _et al._
(2013, 2011a, 2011b). Nonetheless there is one important difference between
our work and the effects due to spin-rotation coupling. The effect in the
present requires a spin-orbit coupling, and the resulting spin polarization
depends on the spin-orbit coupling, whereas the spin-rotation coupling does
not depend much on spin-orbit coupling. Thus the result is the present papaer
cannot be attributed to the spin-rotation coupling in Refs. Matsuo  _et al._
(2013, 2011a, 2011b). We also note that a relateded phenomenon is studied for
an orbital magnetization in Ref. Trifunovic  _et al._ (2019). It is shown that
a geometric orbital magnetization appears in a periodic adiabatic process. It
does not require spin-orbit coupling. Meanwhile, our result shows that
generation of spin magnetization via rotations of phonons require spin-orbit
coupling.

###### Acknowledgements.

We thank Keiji Saito for useful discussions. This work was partly supported by
a Grant-in-Aid for Scientific Research on Innovative Area, “Nano Spin
Conversion Science” (Grant No. 26100006), by MEXT Elements Strategy Initiative
to Form Core Research Center (TIES), and also by JSPS KAKENHI Grant Numbers
JP17J10342 and 18H03678.

## References

  * Einstein and de Haas (1915) A. Einstein and W. J. de Haas, Phys. Ges. 17, 152 (1915).
  * Barnett (1915) S. J. Barnett, Phys. Rev. 6, 239 (1915).
  * Matsuo  _et al._ (2013) M. Matsuo, J. Ieda, K. Harii, E. Saitoh, and S. Maekawa, Phys. Rev. B 87, 180402(R) (2013).
  * Matsuo  _et al._ (2011a) M. Matsuo, J. Ieda, E. Saitoh, and S. Maekawa, Phys. Rev. Lett. 106, 076601 (2011a).
  * Matsuo  _et al._ (2011b) M. Matsuo, J. Ieda, E. Saitoh, and S. Maekawa, Phys. Rev. B 84, 104410 (2011b).
  * Kobayashi  _et al._ (2017) D. Kobayashi, T. Yoshikawa, M. Matsuo, R. Iguchi, S. Maekawa, E. Saitoh, and Y. Nozaki, Phys. Rev. Lett. 119, 077202 (2017).
  * Hamada  _et al._ (2015) M. Hamada, T. Yokoyama, and S. Murakami, Phys. Rev. B 92, 060409(R) (2015).
  * Takahashi  _et al._ (2016) R. Takahashi, M. Matsuo, M. Ono, K. Harii, H. Chudo, S. Okayasu, J. Ieda, S. Takahashi, S. Maekawa, and E. Saitoh, Nature Physics 12, 52 (2016).
  * Zhang and Niu (2014) L. Zhang and Q. Niu, Phys. Rev. Lett. 112, 085503 (2014).
  * Strohm  _et al._ (2005) C. Strohm, G. L. J. A. Rikken, and P. Wyder, Phys. Rev. Lett. 95, 155901 (2005).
  * Inyushkin and Taldenkov (2007) A. V. Inyushkin and A. N. Taldenkov, JETP Letters 86, 379 (2007).
  * Sugii  _et al._ (2017) K. Sugii, M. Shimozawa, D. Watanabe, Y. Suzuki, M. Halim, M. Kimata, Y. Matsumoto, S. Nakatsuji, and M. Yamashita, Phys. Rev. Lett. 118, 145902 (2017).
  * Kagan and Maksimov (2008) Y. Kagan and L. A. Maksimov, Phys. Rev. Lett. 100, 145902 (2008).
  * Zhang  _et al._ (2010) L. Zhang, J. Ren, J.-S. Wang, and B. Li, Phys. Rev. Lett. 105, 225901 (2010).
  * Garanin and Chudnovsky (2015) D. A. Garanin and E. M. Chudnovsky, Phys. Rev. B 92, 024421 (2015).
  * Nakane and Kohno (2018) J. J. Nakane and H. Kohno, Phys. Rev. B 97, 174403 (2018).
  * Streib  _et al._ (2018) S. Streib, H. Keshtgar, and G. E. W. Bauer, Phys. Rev. Lett. 121, 027202 (2018).
  * Juraschek and Spaldin (2019) D. M. Juraschek and N. A. Spaldin, Phys. Rev. Materials 3, 064405 (2019).
  * Holanda  _et al._ (2018) J. Holanda, D. S. Maior, A. Azevedo, and S. M. Rezende, Nature Physics 14, 500 (2018).
  * Guerreiro and Rezende (2015) S. C. Guerreiro and S. M. Rezende, Phys. Rev. B 92, 214437 (2015).
  * Zhang and Niu (2015) L. Zhang and Q. Niu, Phys. Rev. Lett. 115, 115502 (2015).
  * Zhu  _et al._ (2018) H. Zhu, J. Yi, M.-Y. Li, J. Xiao, L. Zhang, C.-W. Yang, R. A. Kaindl, L.-J. Li, Y. Wang, and X. Zhang, Science 359, 579 (2018).
  * Juraschek  _et al._ (2017) D. M. Juraschek, M. Fechner, A. V. Balatsky, and N. A. Spaldin, Phys. Rev. Materials 1, 014401 (2017).
  * Wang  _et al._ (2015) Y.-T. Wang, P.-G. Luan, and S. Zhang, New Journal of Physics 17, 073031 (2015).
  * Liu  _et al._ (2017) Y. Liu, Y. Xu, S.-C. Zhang, and W. Duan, Phys. Rev. B 96, 064106 (2017).
  * Nova  _et al._ (2017) T. F. Nova, A. Cartella, A. Cantaluppi, M. Först, D. Bossini, R. V. Mikhaylovskiy, A. V. Kimel, R. Merlin, and A. Cavalleri, Nature Physics 13, 132 (2017).
  * Hamada  _et al._ (2018) M. Hamada, E. Minamitani, M. Hirayama, and S. Murakami, Phys. Rev. Lett. 121, 175301 (2018).
  * (28) M. Hamada and S. Murakami (unpublished).
  * Trifunovic  _et al._ (2019) L. Trifunovic, S. Ono, and H. Watanabe, Phys. Rev. B 100, 054408 (2019).
  * Berry (1987) M. V. Berry, Proc. R. Soc. London. A 414, 31 (1987).
  * Rigolin  _et al._ (2008) G. Rigolin, G. Ortiz, and V. H. Ponce, Phys. Rev. A 78, 052508 (2008).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005).
  * Sasaki and Saito (2008) K.-i. Sasaki and R. Saito, Prog. Theor. Phys. Suppl. 176, 253 (2008).
  * Morimoto  _et al._ (2017) T. Morimoto, H. C. Po, and A. Vishwanath, Phys. Rev. B 95, 195155 (2017).
  * Xu and Wu (2018) S. Xu and C. Wu, Phys. Rev. Lett. 120, 096401 (2018).


#  NeuroBlend: Towards Low-Power yet Accurate Neural Network-Based Inference
Engine Blending Binary and Fixed-Point Convolutions

Arash Fayyazi  University of Southern CaliforniaLos AngelesCaliforniaUSA
<EMAIL_ADDRESS>,  Mahdi Nazemi  University of Southern CaliforniaLos
AngelesCaliforniaUSA<EMAIL_ADDRESS>,  Arya Fayyazi  University of Southern
CaliforniaLos AngelesCaliforniaUSA<EMAIL_ADDRESS>and  Massoud Pedram
University of Southern CaliforniaLos AngelesCaliforniaUSA<EMAIL_ADDRESS>
(2024)

###### Abstract.

This paper introduces NeuroBlend, a novel neural network architecture
featuring a unique building block known as the Blend module. This module
incorporates binary and fixed-point convolutions in its main and skip paths,
respectively. There is a judicious deployment of batch normalizations on both
main and skip paths inside the Blend module and in between consecutive Blend
modules. Additionally, we present a compiler and hardware architecture
designed to map NeuroBlend models onto FPGA devices, aiming to minimize
inference latency while maintaining high accuracy. Our NeuroBlend-20
(NeuroBlend-18) model, derived from ResNet-20 (ResNet-18) trained on CIFAR-10
(CIFAR-100), achieves 88.0% (73.73%) classification accuracy, outperforming
state-of-the-art binary neural networks by 0.8% (1.33%), with an inference
time of 0.38ms per image, 1.4x faster than previous FPGA implementation for
BNNs. Similarly, our BlendMixer model for CIFAR-10 attains 90.6%
accuracy(1.59% less than full precision MLPMixer), with a 3.5x reduction in
model size compared to full precision MLPMixer. Furthermore, leveraging DSP
blocks for 48-bit bitwise logic operations enables low-power FPGA
implementation, yielding a 2.5x reduction in power consumption.

Binary neural networks, FPGA, Low power, MLPMixer

††journalyear: 2024††copyright: rightsretained††conference: Great Lakes
Symposium on VLSI 2024; June 12–14, 2024; Clearwater, FL, USA††booktitle:
Great Lakes Symposium on VLSI 2024 (GLSVLSI ’24), June 12–14, 2024,
Clearwater, FL, USA††doi: 10.1145/3649476.3658736††isbn:
979-8-4007-0605-9/24/06††ccs: Computing methodologies Machine learning††ccs:
Hardware Hardware accelerators††ccs: Hardware Reconfigurable logic and FPGAs

##  1\. Introduction

Deep neural networks (DNNs) have surpassed the accuracy of conventional
machine learning models in many challenging domains, including computer vision
(Huang et al., 2017) and natural language processing (NLP) (Devlin et al.,
2019). Recently, inspired by the successes in NLP, transformers (Vaswani et
al., 2017) are adopted by the computer vision community. Built with self-
attention layers, multi-layer perceptrons (MLPs), and skip connections,
transformers make numerous breakthroughs on visual tasks (Liu et al., 2021).
To reduce the transformer model complexity, MLPMixers (Tolstikhin et al.,
2021), which replace the multi-head self-attention module in transformers with
a two-layer spatial MLP, are introduced.

Unfortunately, many DNN-based inference engines have a high latency cost and
use enormous hardware resources, which, in turn, prevent their deployment in
latency-critical applications, especially on resource-constrained platforms.
The high latency and large hardware cost are due to the fact that practical,
high-quality deep learning models entail billions of arithmetic operations and
millions of parameters, which exert considerable pressure on both the
processing and memory subsystems.

Quantization has emerged as a promising model compression method where
parameters and/or activations are replaced with low-precision, quantized,
fixed-point values. Despite such a transformation, quantized models can match
the accuracy of full-precision models utilizing 32-bit floating-point (FP-32)
while requiring fewer data transfers and storage space. Early works on binary
neural networks (BNNs) (Hubara et al., 2016) and XNOR-net (Rastegari et al.,
2016) demonstrated the potential advantages of extreme quantization, i.e.,
binarization. BNNs are 1-bit quantized models where all weights and
activations are represented by two values, -1/0 and +1, significantly
decreasing the memory footprint. Additionally, to speed up the inference, the
multiplication/addition operations are switched out for less complex
operations like the XNOR logical operation and bit count (Rastegari et al.,
2016). However, this superior performance is achieved at the cost of a
significant accuracy drop in deep neural networks.

To address the issue of significant accuracy loss, some prior work (Liu et
al., 2018; Martínez et al., 2020; Nazemi et al., 2021; Liu et al., 2020; Zhang
et al., 2021) propose to modify the well-known architectures (e.g., ResNet) as
they show the network architecture can affect BNN performance. Although some
of these works can achieve improved accuracy, their proposed models cannot be
efficiently deployed on hardware platforms. For instance, ReActNet (Liu et
al., 2020) significantly improves the accuracy by activation shifting and
reshaping the MobileNet V1 (Howard et al., 2017) architecture at the cost of
increasing the parameters such that the total number of parameters is about 30
million more than the number of parameters in MobileNet V2 (Sandler et al.,
2018).

This paper presents a novel building block called Blend module, which utilizes
binary convolution on its main path and fixed-point convolution on its skip
path. The key contributions of this work can be summarized as follows.

  * •

We present NeuroBlend, a hardware-friendly neural network architecture with
binary activations, where all convolutional layers are computed using binary
multiply-accumulate (BMAC) operations (except on skip paths that utilize
fixed-point convolutions). On a ResNet-20-like architecture designed for the
CIFAR-10 dataset, NeuroBlend outperforms the state-of-the-art binary-based
implementation (Zhang et al., 2021) by 0.8% in top-1 accuracy.

  * •

We introduce a powerful and flexible compiler for mapping a given NeuroBlend
inference engine running on any dataset onto our optimized accelerator design
by converting the network model to a computational graph, scheduling the
graph’s execution, and optimizing its nodes by leveraging intrinsic fusions of
the convolution and batch normalization layers.

  * •

We present a flexible FPGA-based design that enables the DSP-based realization
of BMAC operations. The reconfigurability of DSPs for performing bitwise logic
operations is utilized to achieve a low-power implementation.

  * •

We apply our transformations, integrate our blocks to MLPMixers, and improve
the naively binarized MLPMixer models by 6%. To the best of our knowledge,
this is the first paper that presents a binary model of MLPMixers with
negligible accuracy drop.

##  2\. Preliminaries

(a)

(b)

Figure 1. Overview of MLPMixer model. (a) Overall MLPMixer and (b) mixing
block architectures.

In this section, we briefly introduce the MLPMixer, the conventional BNN
model, and discuss some recent BNN architectures that have significantly
improved accuracy.

###  2.1. MLPMixer

To further reduce the inductive biases introduced by CNNs, MLP-Mixer has
recently proposed a more straightforward solution that is fully based on
multi-layer perceptrons (MLPs) (Tolstikhin et al., 2021) (see fig. 1). The
basic layer in MLP-Mixer consists of two components: the channel-mixing block
and the token-mixing block. Each of these mixing blocks has the same units and
is shown in Fig. 1(b). In the channel-mixing block, the feature map is
projected along the channel dimension for the communications among various
channels, while the feature map is projected along the spatial dimension, and
communications among spatial locations are accomplished concurrently by the
token-mixing block.

###  2.2. Conventional BNN

BNNs have several properties that enable more efficient mapping to FPGAs
without affecting network accuracy. For implementing conventional BNN models
that are binarized (both weights and activations are quantized to 1-bit
values), the product of a binary weight and activation can be replaced with a
binary XNOR operation.

Furthermore, by assuming that an unset bit represents -1 and a set bit
represents +1, there are only two possible values of +1 and -1 for the result
of the XNOR operation and, thus, synapse input. Therefore, the summation of a
binary dot product can be implemented by a popcount operation that counts the
number of set bits instead of accumulation with signed arithmetic.
Furthermore, all BNN layers use batch normalization on convolutional or fully
connected layer outputs and then apply the sign function to determine the
output activation. For hardware implementation, the same output activation can
be computed via thresholding (Umuroglu et al., 2017). Lastly, Max-pooling on
binary activations can be implemented in hardware using the Boolean OR
function.

###  2.3. Prior Work on State-of-the-art BNN Architectures

The state-of-the-art BNN models achieved a high accuracy using blocks similar
to ResNet models (He et al., 2016). Martinez et al. (Martínez et al., 2020)
presented a strong baseline for BNN models which is based on the modified
ResNet block suitable for 1-bit CNNs. They used double skip connections and
PReLU (PReLU was first introduced in (Bulat et al., 2019)) as the activation
function. They follow the idea of using real-valued downsample layers proposed
in (Liu et al., 2018) that improves the accuracy significantly. They also
presented a method wherein the architectural gap between real and binary
networks is, step by step, bridged via a sequence of teacher-student pairs.

More recently, a more accurate BNN model is proposed called ReActNet (Liu et
al., 2020) to mitigate the precision gap between the binarized model and its
counterpart of real-valued. ReActNet took one step further and is based on
MobileNetV1 (Howard et al., 2017) architecture. It reaches a top-1 accuracy of
69.4% in the IMAGENET data set using 4.82G BMACS with a 4.6 MB model size. The
key block in ReActNet is a biased PReLU (BPReLU) activation function that
shifts and reshapes the feature maps between two convolutional layers. This
substantially improves the accuracy of the model.

To improve ReActNet (Liu et al., 2020), The authors in (Zhang et al., 2021)
have proposed FracBNN which employs a dual-precision activation scheme to
compute features with up to two bits, using an additional sparse binary
convolution layer. They have achieved MobileNetV2-level accuracy with
competitive model size.

##  3\. The Proposed Building Blocks

The proposed BNN model comprises two types of building blocks, as shown in
Fig. 2: one with no operations on the skip path and another with average
pooling, convolution, and batch normalization on its skip path. There are a
few differences between the presented building blocks compared to prior work.

First, both types of building blocks include a batch normalization layer as
their final output layer. This batch normalization layer, which does not
include any data-driven, trainable channel rescaling/shifting parameters,
ensures the output of each block is in a range centered around zero, and is,
therefore, amenable to fixed-point or binary quantization. Adding channel
rescaling/shifting to this batch normalization layer tends to reduce the
classification accuracy by a few percentage points.

Second, a PReLU activation is placed in the main path and before the addition
operation, which yields improvements in classification accuracy compared to
other activation functions such as ReLU (Martínez et al., 2020).

Last but not least, the building blocks are designed with the hardware
implementation cost in mind. For example, as shown in Fig. 4, the final batch
normalization layer of a building block can be fused into the next building
block, resulting in a thresholding operation in the main path and modified
convolutional weights in the skip path (see details in Sec. 4.1). A similar
fusing can be performed for batch normalization in the skip path, all leading
to reduced end-to-end inference latency. We extend our building blocks to the
MLPMixer model and propose a new Mixing block, as shown in Fig. 3.

(a) Normal block

(b) Downsample block

Figure 2. The proposed building blocks. The differences with respect to real-
to-binary blocks are highlighted in red.

##  4\. Proposed Accelerator Design

In this section, we describe the proposed hardware accelerator and the
associated compiler to perform optimization tailored to the employed
accelerator design.

###  4.1. Compiler Optimization

Our compiler performs optimizations tailored to the employed accelerator
design and fuses operations to reduce the hardware cost and generate an
intermediate graph. Finally, our compiler compiles the intermediate graph,
extracts the required parameters for the accelerator design, and generates a
static schedule.

All BNN layers proposed in this paper (cf. Fig 2) begin with a sign function
(SG) and end with a batch normalization (BN) operation. First, we move the
last BN block of a layer to the next layer and pass it through both feed-
forward and skip paths of the next layer (cf. Fig. 4(b). Moreover, the BN
block that comes before the SG block can be replaced by a thresholding (TH)
function in order to reduce the hardware cost (c.f. Fig. 4(c)). Using such a
technique, we can process the input feature map using an unsigned comparison
and avoid expensive operations such as multiplication that are required in the
BN block. Reference (Umuroglu et al., 2017) explains how the hardware cost of
a regular BN-SG block is reduced from 2 DSPs, 55 FFs, and 40 LUTs for separate
BN and SG computations to merely 6 LUTs for the TH block computations using
such a technique.

In addition, the BN block that is passed onto the skip path is also fused with
the CONV layer in the skip path if one exists. In summary, we will replace a
BN-CONV-BN sequence of layers on the skip path with a CONV layer (c.f. Fig.
4(c)). More details of such fusion are provided below. Note that the BN block
that is passed onto the skip path is passed through the average pooling block
in the inference (c.f. Fig. 4(b)). This simplification also reduces the
chances of encountering overflow/underflow even if we assign fewer bits than
required bits for summation because the weights/biases of the fused BN-CONV-BN
layer will help normalize this layer’s outputs. In other words, the output of
BMAC can be considered 16 bits, and no quantization is required.

Algorithm 1 shows the well-known batch normalization algorithm where the two
parameters $\gamma$ and $\beta$ are learned during the training process. Note
that $\epsilon$ is a small constant value used to ensure that division-by-zero
error is not encountered.

Algorithm 1 Batch normalization for activation $y$ in a mini-batch,
$y^{\prime}$ is the normalized result

Input $\mathcal{B}$ = {$y_{0},y_{1},...,y_{m}$}; $\gamma$; $\beta$
Output $y^{\prime}_{i}$ = $BN_{\gamma,\beta}(y_{i})$

1: $\mu_{\mathcal{B}}\leftarrow\frac{1}{m}\sum_{i=1}^{m}y_{i}$ // mini-batch
mean

2:
$\sigma_{\mathcal{B}}^{2}\leftarrow\frac{1}{m}\sum_{i=1}^{m}(y_{i}-\mu_{\mathcal{B}})^{2}$
// mini-batch variance

3:
$\hat{y_{i}}\leftarrow\frac{y_{i}-\mu_{\mathcal{B}}}{\sqrt{\sigma_{\mathcal{B}}^{2}+\epsilon}}$
// normalize

4: $y^{\prime}_{i}\leftarrow\gamma\hat{y_{i}}+\beta\equiv
BN_{\gamma,\beta}(y_{i})$ // scale and shift

Figure 3. The proposed Mixing block in this paper.

When the BN layer is fused into a subsequent convolutional layer, the fused
layer’s computations may be written as,

(1) $\small
y^{\prime}=\mathbold{w}(\gamma\frac{\mathbold{x}-\mathbold{\mu^{\prime}}_{\mathcal{B}}}{\sqrt{\mathbold{\sigma^{\prime}}_{\mathcal{B}}^{2}+\epsilon}}+\mathbold{\beta^{\prime}})+\mathbold{b}$

Hence, the new fused parameters can be calculated as:

(2)
$\small\begin{gathered}\mathbold{w^{\prime}}=\frac{\mathbold{\gamma^{\prime}}\mathbold{w}}{\sqrt{\mathbold{\sigma^{\prime}}_{\mathcal{B}}^{2}+\epsilon}},~{}~{}~{}\mathbold{b^{\prime}}=-\frac{\mathbold{\gamma^{\prime}}\mathbold{w}\mathbold{\mu^{\prime}}_{\mathcal{B}}}{\sqrt{\mathbold{\sigma^{\prime}}_{\mathcal{B}}^{2}+\epsilon}}+\mathbold{w}\mathbold{\beta^{\prime}}+\mathbold{b}\end{gathered}$

Note that $\mathbold{x}$ in Eqn. 1 is the output of previous layer, so the 4-D
$\mathbold{w}$ must be squeezed into a 2-D vector of size ($c_{\mathrm{in}}$,
$c_{\mathrm{out}}$ $\times$ $w_{\mathrm{k}}$ $\times$ $h_{\mathrm{k}}$).
Indeed, if we use zero padding in convolutions, we will have 0 values entering
the convolution after the BN. When we fuse the parameters, this must still be
the case. So, to correctly fold $\mathbold{w}$ into the 2-D vector, we have to
replace the default 0 values for padding by
$\mathbold{\mu^{\prime}}_{\mathcal{B}}-\frac{\mathbold{\beta^{\prime}}\sqrt{\mathbold{\sigma^{\prime}}_{\mathcal{B}}^{2}+\epsilon}}{\mathbold{\gamma^{\prime}}}$,
that is, we must apply the BN transformation to the padding. Hence, by
multiplying the fused weights $\mathbold{w^{\prime}}$ to the padding values,
we apply the inverse of the BN transformation.

Using algorithm 1 and proceeding with fused parameters as in the previous
case, the weights and bias of the resulting fused BN-CONV-BN block may be
expressed as:

(3)
$\small\begin{gathered}\mathbold{w^{\prime\prime}}=\frac{\mathbold{\gamma^{\prime\prime}}\mathbold{w^{\prime}}}{\sqrt{\mathbold{\sigma^{\prime\prime}}_{\mathcal{B}}^{2}+\epsilon}},~{}~{}~{}\mathbold{b^{\prime\prime}}=\mathbold{\beta^{\prime\prime}}+\mathbold{\gamma^{\prime\prime}}\frac{\mathbold{b^{\prime}}-\mathbold{\mu^{\prime\prime}}_{\mathcal{B}}}{\sqrt{\mathbold{\sigma^{\prime\prime}}_{\mathcal{B}}^{2}+\epsilon}}\end{gathered}$

(a)

(b)

(c)

Figure 4. (a) The proposed building block. (b) Distribute the BN from the
precedent NeuroBlend block to both main and skip paths. (c) Merge components
that result in the optimized NeuroBlend block. Note that the red and yellow
dash boxes show the precedent and current building blocks. The green dash
boxes show the potential components for merging.

###  4.2. Accelerator Architecture

We adopted a heterogeneous streaming architecture where each layer uses its
own hardware resources for this work.

Our Accelerator consists of the following types of operations: 1) $3\times 3$
BNN convolution, 2) $1\times 1$ Fixed-Point convolution, 3) Average pooling,
4) Max Pooling, 5) Linear layer, 6) Thresholding, 7) BN-PReLU/BN-PReLU-BN, and
8) Residual connection and summation.

We have designed and implemented an efficient accelerator that supports these
operations on FPGA devices. We separate our design into three domains: binary,
fixed-point, and joint domains. The BNN (FPNN) convolutions are placed in
binary (fixed-point) domains, whereas all other operations are in the joint
domain. Our hardware is designed by using Xilinx’s High-Level Synthesis (HLS)
tools. In the following, we will describe the hardware engines and the HLS
design techniques in detail.

####  4.2.1. BNN and FPNN Block

The BMAC/FPMAC of the proposed accelerator are mapped to DSPs. Each DSP block
(i.e., DSP48E2 in Xilinx FPGA) is capable of performing a 18 by 27-bit or a
48-bit bitwise logic operation including AND, OR, NOT, NAND, NOR, XOR, and
XNOR. We used this feature to perform 48 XNOR operations simultaneously as one
BMAC operation. Note that similar 48 XNOR operations on LUTs result in the
usage of 48 LUTs and 49 FFs. So, the input feature map and weights must be
packed into groups of size 48. For this purpose, We pack bits along the
channel dimension into 48-bit unsigned integers for concurrent access. We will
highlight the advantage of mapping logic operations to DSPs in Section 5.3.

Our design results in a balanced usage of hardware resources since LUTs are
mostly used for other operations while DSPs are used to do BMAC/FPMAC
operations. Note that the reconfigurability of DSPs in each clock cycle can be
used to switch between BMAC and FPMAC operations. This is very helpful in the
case of a homogeneous single execution engine design where we only use a fixed
set of resources for all layers, although this is not in the scope of the
present paper. Note that option can be used in single execution engine
architecture. In the streaming architecture, each unit has its own resource. A
2D array of DSPs with the size of 32 * 32 is designed to perform FPMAC
operations.

####  4.2.2. FP/B Joint Blocks

In addition to the accelerators for the convolution operations, which account
for the majority of the computations in a vision neural network, we also
implement hardware accelerators for other operations that must be performed in
the joint domain, including the max (average) pooling and thresholding blocks.
For instance, the TH unit compares each output activation from the previous
layer with a programmable threshold and outputs +1 (0) if the output
activation is greater (smaller) than the corresponding threshold. Since these
TH blocks only contain channel-wise parameters, their impact on the total
number of model parameters is negligible. Although we can achieve the maximum
concurrency by processing all output activations simultaneously, this would
require a lot of resources, including both memories (we must increase the
number of access lines by partitioning the data to several memory blocks) and
TH block. The performance gain may not be worth the cost of additional
resources needed. We use the greatest common divisor (GCD) of the height of
the systolic array in the FPNN block (e.g., 32) and the width of the array of
BMACs in the BNN block (e.g., 48) as the parallelism factor for all blocks in
the joint domain.

##  5\. Experimental Results

For evaluation purposes, we targeted a high-end Virtex® UltraScale+ FPGA
(Xilinx VU9P FPGA, which is available in the cloud as the AWS EC2 F1
instance). This FPGA contains approximately 2.5 million logic elements and
approximately 6,800 DSP units.We use Vitis 2020.2 for mapping to the FPGA and
set the target frequency to 340 Mhz. We also use the Vivado power report
provided by Xilinx to assess the power consumption of each design. Finally, we
evaluate our proposed method on a well-known CNN, i.e., ResNet-20 (He et al.,
2016) and MLPMixers (Tolstikhin et al., 2021) and a commonly used computer-
vision dataset for object recognition, i.e., the CIFAR-10 (Krizhevsky et al.,
2009) dataset.

###  5.1. Experimental setup

In the case of MLPMixers, the resolution of the input image is 32*32, and the
patch size that the experiments are based on is 4*4. So, we have S non-
overlapping image patches that are mapped to a hidden dimension $C$. $D_{S}$
and $D_{C}$ are tunable hidden widths in the token-mixing and channel-mixing
MLPs, respectively. The summary of design specifications is shown in Table 1.

Table 1. Summary of design specifications for MLPMixers used in this paper. The ”S” and ”B” (small and base) models scales down follow Tolstikhin et al. (Tolstikhin et al., 2021). The notation ”B/4” means the model of base scale with patches of resolution 4*4. Specification | S/4 | B/4 | 2S/4
---|---|---|---
Sequence length S | 64 | 64 | 64
Hidden size C | 128 | 192 | 256
Patch resolution | 4*4 | 4*4 | 4*4
MLP dimension $D_{C}$ | 512 | 768 | 1024
MLP dimension $D_{S}$ | 64 | 96 | 128
Number of layers | 8 | 12 | 8

###  5.2. ResNet-18 and ResNet-20

Table 2 demonstrates the superiority of NeuroBlend on ResNet-20 and ResNet-18
compared to other binary-based approaches on CIFAR-10 and CIFAR-100 datasets.
We can observe that NeuroBlend achieves a higher top-1 accuracy by 0.8% for
the CIFAR-10 dataset and 1.33% for the CIFAR-100 dataset compared to the
state-of-the-art BNNs and improves the accuracy by 1.5% to about 9% compared
to the other approaches. NeuroBlend-20 achieves the same model size as FracBNN
even while keeping the first and last layer full precision (i.e., 16-bit
fixed-point).

Table 2. Classification accuracy of ResNet-20-like and Resnet-18 models on
CIFAR-10 and CIFAR-100, respectively.

CIFAR-10 | CIFAR-100
---|---
Approach | Accuracy (%) | Approach | Accuracy (%)
DoReFa-Net (Zhou et al., 2016) | 79.3 | Bi-RealNet-18 (Liu et al., 2018) | 63.51
DSQ (Gong et al., 2019) | 84.1 | ReActNet-18 (Liu et al., 2020) | 70.76
ReActNet (Liu et al., 2020) | 85.8 | AresB-Net-18(Kim, 2021) | 71.98
IR-Net (Qin et al., 2020) | 86.5 | PresB-Net-18(Shin and Kim, 2022) | 72.30
FracBNN (Zhang et al., 2021) | 87.2 | Hyper-BinaryNet(Wang et al., 2021) | 72.40
NeuroBlend-20 | 88.0 | NeuroBlend-18 | 73.73

Table 3. Comparison between the hardware metrics of different types of
NeuroBlend-20 implementation on CIFAR-10.

Approach | DSP48E2 | LUT | FF | BRAM_18K |  $f_{target}$ | Power
---|---|---|---|---|---|---
DSP-based NeuroBlend-20 |  2240 (32%) |  10K (1%) |  30K (1%) |  912 (42%) |  342 MHz |  6.1 w
LUT-based NeuroBlend-20 |  1132 (17%) |  500K (50%) |  900K (34%) |  912 (42%) |  342 MHz |  15.3 w

###  5.3. Hardware cost and performance of NeuroBlend

In this section, we evaluate our NeuroBlend-20 model on the FPGA platform.
Compared to the BNN accelerator in FracBNN (Zhang et al., 2021), our design
achieves a higher frame rate (3846 FPS vs. 2807 FPS reported in FracBNN (Zhang
et al., 2021)) and higher working frequency (342 MHz vs. 250 MHz) while
yielding higher accuracy. Note that the FPGA used in this paper is a server-
class FPGA while authors in (Zhang et al., 2021) deployed their design on an
embedding FPGA with fewer resources. However, the frame rate comparison is
fair since they unroll the entire model, which is similar to what we achieved.
The unrolling is feasible because ResNet-20-based models on CIFAR-10 are very
compact and can be fitted into an FPGA to eliminate unnecessary transactions
between the logic blocks and the DDR memory.

Table 3 compares the hardware cost of two approaches for implementing the
NeuroBlend-20 on the CIFAR-10 dataset where the DSP-based NeuroBlend-20 is the
approach presented in Section 4.2.1. The LUT-based is a naive implementation
where logic operations are performed using LUTs. Note that DSP usage in LUT-
based implementation is for the first and last layers. Decreasing the LUT
usage results in power saving. Our measurements show that the proposed
implementation yields 2.5x lower power consumption.

###  5.4. MLPMixers

As table 4 shows, our model (i.e., BlendMixer) outperforms the naively
binarized version of MLPMixer (i.e., BinaryMixer) yet achieves a smaller
memory footprint (i.e., model size). Increasing the model size (e.g., model
B/4) reduces the accuracy drop due to binarization. Note that FPMAC operations
in both BinaryMixer and BlendMixer are due to the patch embedding block and
the last FC layer. When it comes to analyzing the performance, BlendMixer can
achieve 0.02ms (0.06ms) latency for a small (base) model, while MLPMixer can
process an image in 0.24ms (1.14ms).

###  5.5. Ablation study with normalization unit in mixing blocks

To further evaluate the influence of the normalization unit on the proposed
Mixing blocks, we performed an ablation study with different normalization
units, layer normalization as suggested in the paper, batch normalization over
channels, and batch normalization over patches. As shown in Table 5, the
Mixing block with batch normalization over channels outperforms other models.
Using batch normalization also reduces the memory footprint, requiring storing
fewer parameters. Note that the MLPMixer consists of the proposed mixing
blocks called BlendMixer.

###  5.6. Computation/memory cost vs accuracy in mixing blocks

In this section, we assess the trade-off between computation/memory complexity
and result accuracy. As demonstrated in Table 6, we can manipulate the type of
computation in FC layers of mixing blocks to improve the accuracy at the cost
of computation and memory complexity. The BlendMixer(BB/BB-2S/4) is simply the
widened version of BlendMixer S/4, and BlendMixer(BB/FPB-2S/4) is the same as
the previous architecture with the exception that the FC layer of the second
mixing block being calculated in FP format. The other models are named
accordingly. The BlendMixer(BB/FPB-2S/4) achieves a comparable accuracy (less
than 1% accuracy drop compared to MLPMixer S/4).

Table 4.  Comparison between the mapping metrics of BlendMixer with those of
MLPMixer (Tolstikhin et al., 2021) and BinaryMixer on CIFAR-10 dataset.

Model Spec | Model |  Precision (W/A) |  Model Size (MB) |  FPMAC ($\times 10^{8}$) |  BMAC ($\times 10^{8}$) |  Top-1 (%)
---|---|---|---|---|---|---
S/4 | MLPMixer | 16/16 | 2.5 | 0.76 | 0 | 92.25
BinaryMixer | 1/1 | 0.41 | 0.005 | 0.75 | 74.38
BlendMixer | 1/1 | 0.15 | 0.005 | 0.75 | 80.43
B/4 | MLPMixer | 16/16 | 8.41 | 3.83 | 0 | 92.93
BinaryMixer | 1/1 | 1.07 | 0.015 | 3.82 | 82.5
BlendMixer | 1/1 | 0.49 | 0.015 | 3.82 | 87.34

Table 5. Accuracy comparison between BlendMixer S/4 with different
normalization unit.

Normalization unit | LayerNormalization |  BatchNormalization (over channels) |  BatchNormalization (over patches)
---|---|---|---
Accuracy | 90.60 | 91.96 | 90.57

Table 6.  Comparison between different design choices for trading-off
computation cost for accuracy. Each layer in BlendMixer has two mixing blocks,
each of which contains two FCs. Symbols B and FP show the implementation of
FCs where B (FP) is binary (fixed-point) implementation.

Model |  Model Size (MB) |  FPMAC ($\times 10^{8}$) |  BMAC ($\times 10^{8}$) |  Top-1 (%)
---|---|---|---|---
BlendMixer(BB/BB-2S/4) | 0.58 | 0.021 | 6.03 | 88.39
BlendMixer(BB/FPB-2S/4) | 2.42 | 0.063 | 3.36 | 91.35
BlendMixer(BB/BFP-2S/4) | 2.42 | 0.063 | 3.36 | 91.28
BlendMixer(FPB/BB-2S/4) | 0.70 | 0.023 | 5.7 | 90.66
BlendMixer(BFP/BB-2S/4) | 0.70 | 0.023 | 5.7 | 90.53

##  6\. Conclusion

This paper introduces NeuroBlend, an innovative neural network design that
makes use of the Blend module, a new building block that performs binary and
fixed-point convolutions in the main and skip paths, respectively. Batch
normalizations are deployed intelligently on both the main and skip paths
within the Blend module as well as between adjacent Blend modules.
NeuroBlend-20, a descendant of ResNet-20 trained on the CIFAR-10 dataset,
achieves 88.0% classification accuracy (0.8 percent better than the state-of-
the-art binary neural network), yet 1.4x higher throughput. Finally,
BlendMixer, our new block inspired by MLPMixer blocks, outperforms the naively
binarized version of MLPMixer yet achieves a smaller memory footprint.

## References

  * (1)
  * Bulat et al. (2019) Adrian Bulat et al. 2019\.  Improved training of binary networks for human pose estimation and image recognition.  _CoRR_ abs/1904.05868 (2019).  arXiv:1904.05868
  * Devlin et al. (2019) Jacob Devlin et al. 2019\.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding. In _Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_.
  * Gong et al. (2019) Ruihao Gong et al. 2019\.  Differentiable Soft Quantization: Bridging Full-Precision and Low-Bit Neural Networks. In _International Conference on Computer Vision_.
  * He et al. (2016) Kaiming He et al. 2016\.  Deep Residual Learning for Image Recognition. In _Conference on Computer Vision and Pattern Recognition_. IEEE.
  * Howard et al. (2017) Andrew G. Howard et al. 2017\.  MobileNets: Efficient Convolutional Neural Networks for Mobile Vision Applications.  (2017).  arXiv:1704.04861
  * Huang et al. (2017) Gao Huang et al. 2017\.  Densely Connected Convolutional Networks. In _Conference on Computer Vision and Pattern Recognition_. IEEE Computer Society.
  * Hubara et al. (2016) Itay Hubara et al. 2016\.  Binarized Neural Networks. In _Advances in Neural Information Processing Systems_.
  * Kim (2021) HyunJin Kim. 2021.  AresB-Net: accurate residual binarized neural networks using shortcut concatenation and shuffled grouped convolution.  _PeerJ Computer Science_ 7 (2021), e454.
  * Krizhevsky et al. (2009) Alex Krizhevsky, Geoffrey Hinton, et al. 2009\.  Learning multiple layers of features from tiny images.  (2009).
  * Liu et al. (2018) Zechun Liu et al. 2018\.  Bi-Real Net: Enhancing the Performance of 1-Bit CNNs with Improved Representational Capability and Advanced Training Algorithm. In _European Conference on Computer Vision_.
  * Liu et al. (2020) Zechun Liu et al. 2020\.  ReActNet: Towards Precise Binary Neural Network with Generalized Activation Functions. In _European Conference on Computer Vision_.
  * Liu et al. (2021) Ze Liu et al. 2021\.  Swin Transformer: Hierarchical Vision Transformer using Shifted Windows. In _International Conference on Computer Vision_.
  * Martínez et al. (2020) Brais Martínez et al. 2020\.  Training binary neural networks with real-to-binary convolutions. In _8th International Conference on Learning Representations, ICLR,_.
  * Nazemi et al. (2021) Mahdi Nazemi et al. 2021\.  NullaNet Tiny: Ultra-low-latency DNN Inference Through Fixed-function Combinational Logic. In _29th IEEE FCCM_.
  * Qin et al. (2020) Haotong Qin et al. 2020\.  Forward and Backward Information Retention for Accurate Binary Neural Networks. In _Conference on Computer Vision and Pattern Recognition, CVPR_. Computer Vision Foundation / IEEE.
  * Rastegari et al. (2016) Mohammad Rastegari et al. 2016\.  XNOR-Net: ImageNet Classification Using Binary Convolutional Neural Networks. In _ECCV_. Springer.
  * Sandler et al. (2018) Mark Sandler et al. 2018\.  MobileNetV2: Inverted Residuals and Linear Bottlenecks. In _IEEE Conference on Computer Vision and Pattern Recognition, CVPR_.
  * Shin and Kim (2022) Jungwoo Shin and HyunJin Kim. 2022.  PresB-Net: parametric binarized neural network with learnable activations and shuffled grouped convolution.  _PeerJ Computer Science_ 8 (2022), e842.
  * Tolstikhin et al. (2021) Ilya O. Tolstikhin et al. 2021\.  MLP-Mixer: An all-MLP Architecture for Vision. In _Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems, NeurIPS_.
  * Umuroglu et al. (2017) Yaman Umuroglu et al. 2017\.  FINN: A Framework for Fast, Scalable Binarized Neural Network Inference. In _Proceedings of the ACM/SIGDA International Symposium on Field-Programmable Gate Arrays, FPGA_. ACM.
  * Vaswani et al. (2017) Ashish Vaswani et al. 2017\.  Attention Is All You Need.  _CoRR_ abs/1706.03762 (2017).  arXiv:1706.03762
  * Wang et al. (2021) Qi Wang, Nianhui Guo, Zhitong Xiong, Zeping Yin, and Xuelong Li. 2021.  Gradient matters: Designing binarized neural networks via enhanced information-flow.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 44, 11 (2021), 7551–7562.
  * Zhang et al. (2021) Yichi Zhang et al. 2021\.  FracBNN: Accurate and FPGA-Efficient Binary Neural Networks with Fractional Activations. In _FPGA_. ACM.
  * Zhou et al. (2016) Shuchang Zhou et al. 2016\.  DoReFa-Net: Training Low Bitwidth Convolutional Neural Networks with Low Bitwidth Gradients.  (2016).  arXiv:1606.06160