# Electrostatic interactions in twisted bilayer graphene

Tommaso Cea  Imdea Nanoscience, Faraday 9, 28015 Madrid, Spain  Pierre A.
Pantaleón  Imdea Nanoscience, Faraday 9, 28015 Madrid, Spain  Niels R. Walet
Department of Physics and Astronomy, University of Manchester, Manchester, M13
9PY, UK  Francisco Guinea  Imdea Nanoscience, Faraday 9, 28015 Madrid, Spain
Donostia International Physics Center, Paseo Manuel de Lardizábal 4, 20018 San
Sebastián, Spain

(August 29, 2024)

###### Abstract

The effects of the long range electrostatic interaction in twisted bilayer
graphene are described using the Hartree-Fock approximation. The results show
a significant dependence of the band widths and shapes on electron filling,
and the existence of broken symmetry phases at many densities, either
valley/spin polarized, with broken sublattice symmetry, or both.

##  I Introduction

The discovery of a complex phase diagram in twisted bilayer graphene [1, 2],
which includes superconducting and insulating phases, has led to an extensive
study of the role of interactions in the electronic structure of this system
(see also[3, 4, 5]). The existence of narrow bands can be obtained from a
continuum model for a twisted bilayer graphene bilayer [6] (see also [7, 8]),
as reported in [9]. Similar results have been obtained using lattice based
tight binding models [10, 11]. The two techniques tend to give consistent
results. Twisted bilayer graphene forms a moiré pattern with the lattice
period $L=a/[2\sin(\theta/2)]$, where $a\approx 2.46\text{\AA}$ is the lattice
constant of graphene, and $\theta$ is the twist angle expressed in radians.
The reported narrow bands occur for twist angles $\theta\approx 1^{\circ}$,
and thus for a moiré lattice vector of $L\approx a/\theta\approx
15\,\text{nm}$. Typical widths of the bands nearest to the Fermi energy range
from $W\sim 1-5$ meV. The central bands accommodate an electron density of
order $-4\times L^{-2}\lesssim 4\times\rho\lesssim\sim L^{-2}$, where
$L^{-2}\sim 10^{12}\mathrm{cm}^{-2}$. The factor $4$ arises from the spin and
valley degeneracy of each band. The fact that $L\gg a$ allows us to classify
the strength of interactions in terms of their dependence on $a/L\sim
10^{-2}$, or, alternatively, in terms of the number of carbon atoms in the
unit cell, $N\sim(L/a)^{2}\sim 10^{4}$. n inhomogeneous charge distribution on
the scale of $L$ leads to electrostatic energies of the order $e^{2}/(\epsilon
L)\sim 10-50$ meV, where $\epsilon\sim 4-10$ is the dielectric constant of the
environment [12]. This energy scales as $L^{-1}\sim 1/\sqrt{N}$. The existence
of atomic orbitals modifies the electron-electron repulsion at length scales
comparable to the lattice constant, $a$. This effect is best described by
including an intra-atomic Hubbard term, $U_{H}\lesssim e^{2}/a\sim 2-6$ eV,
and other inter-atomic interactions. Two electrons in the central bands have a
probability $\sim 1/N^{2}$ of being on the same atom. Hence, the contribution
of a Hubbard term to the energy of the central bands per moiré unit cell is
$\sim U_{H}\times N\times 1/N^{2}\sim U_{H}/N$.

A similar analysis can be carried out for the electron-electron interactions
mediated by phonons. Longitudinal acoustic phonons are coupled to electrons by
the deformation potential, $D\approx 20-30$ eV. In Fourier space, this
coupling translates into an effective electron-electron interaction
$V_{e-e}^{ph_{l}}\sim D^{2}/(\lambda+2\mu)$, where $\lambda$ and $\mu$ are the
Lamé coefficients which describe the elastic properties of a graphene layer.
By averaging this coupling over the moiré scale, we obtain a coupling
$V_{e-e}^{ph_{l}}\sim D^{2}/(\lambda+2\mu)\times L^{-2}\sim 1/N$. Transverse
and optical phonons do not induce local changes in the chemical potential, and
couple to electrons only through a modification of the hopping parameters,
which can be described as an effective gauge field. The dimensionless
parameter which describes this coupling is $\beta=(a/t)\partial t/\partial
a\sim 3$, where $t$ is the nearest neighbor hopping. Transverse acoustic
phonons lead to an effective interaction in Fourier space
$\sim(t^{2}\beta^{2})/\mu$, which translates into an effective interaction
averaged over the moiré unit cell into
$V_{e-e}^{ph_{tr}}\approx(t^{2}\beta^{2})/(a^{2}\mu)\times L^{-2}\sim 1/N$.
For optical phonons, the effective electron-electron coupling defined in
Fourier space is $(t/a)^{2}\beta^{2}(\rho a^{2}\omega_{ph}^{2})^{-1}$, where
$\rho$ is the mass density, and $\omega_{ph}$ is the frequency of the relevant
phonon. In real space, averaging over the moiré unit cell, the previous
quantity defines an energy scale $V_{e-e}^{ph_{opt}}\sim(t^{2}\beta^{2})/(\rho
a^{2}\omega_{ph}^{2})\times L^{-2}\sim 1/N$. A summary of the different
interactions is shown in Table 1.

Interaction |  Energy scale | Dependence on the number of atoms per unit cell
---|---|---
Long range electrostatic interaction | $\displaystyle\frac{e^{2}}{\epsilon L}$ | $\displaystyle\frac{1}{\sqrt{N}}$
Intra-atomic interaction | $U_{H}\times\frac{a^{2}}{L^{2}}$ | $\displaystyle\frac{1}{N}$
Longitudinal acoustic phonons | $\displaystyle\frac{D^{2}}{(\lambda+2\mu)L^{2}}$ | $\displaystyle\frac{1}{N}$
Transverse acoustic phonons | $\displaystyle\frac{t^{2}\beta^{2}}{(a^{2}\mu)}\times\frac{a^{2}}{L^{2}}$ | $\displaystyle\frac{1}{N}$
Optical phonons | $\displaystyle\frac{t^{2}\beta^{2}}{\rho a^{2}\omega_{ph}^{2}}\times\frac{1}{L^{2}}$ | $\displaystyle\frac{1}{N}$
Table 1: Interactions (left), typical energies averaged over the moiré unit
cell (center), and dependence on the number of atoms in the unit cell (right).

The qualitative analysis presented here shows that the leading electron-
electron interaction in twisted bilayer graphene at small angles, $\theta\sim
1^{\circ}$, where the number of atoms in the unit cell is large, $N\sim
10^{4}\gg 1$, is the long range electrostatic interaction. In most crystalline
solids, the leading effects of this interaction can be captured by the
Hartree-Fock approximation. Moreover, twisted bilayer graphene at sufficiently
high temperatures shows a metallic behavior, with features consistent with the
existence of a Fermi surface. This metallic phase can be expected to be
amenable to a mean field, Hartree-Fock, description.

In the following, we describe the main features obtained so far in Hartree-
Fock analyses of twisted bilayer graphene. This approximation is standard in
quantum mechanics, and it gives useful insights in theoretical condensed
matter physics. We will focus mostly on studies carried out using the
continuum model for the electronic structure of twisted bilayer graphene [9],
although lattice-based tight-binding models have also been considered [10,
11]. Using the continuum approximation, the problem is defined in terms of
three quantities with dimensions of energy, which describe the non interacting
band, $v_{F}/L,g_{1}$ and $g_{2}$. The parameter $v_{F}$ is the Fermi
velocity, and $g_{1}$ and $g_{2}$ are averages of the interlayer hopping, see
below. Finally the interactions are described by the parameter
$e^{2}/\epsilon$, where $e$ is the electric charge, and $\epsilon$ is the
dielectric constant of the environment.

This brief review is focused on recent results on the electronic structure of
twisted bilayer graphene, obtained within the Hartree-Fock theory. It is worth
mentioning that the Hartree-Fock approximation has been extensively used for
the analysis of interactions in two dimensional materials [13], including the
energy splittings of flat Landau levels[14]. This paper does not attempt to
cover in full the already extensive literature devoted to this topic. We leave
out analyses beyond mean field, such as RPA calculations of the low energy
excitations [15, 16, 17, 18, 19, 20, 21, 22, 23]. The main conclusions, which
description a number of novel features almost unique to twisted bilayer
graphene, are consistent with most studies carried out to date.

The next section discusses the main features of the model. Then, we describe
the implementation of the Hartree-Fock approximation to the study of long
range interactions in the model. The following section presents the dependence
of the electronic bands on band filling and twist angle. We discuss next
possible broken symmetry phases. Finally, we highlight a number of novel
features induced by long range electrostatic interactions, treated within the
Hartree-Fock approximation,

##  II The Hartree-Fock approximation within the continuum model of TBG

Rotating two layers of graphene by a relative small angle, $\theta$, gives
rise to a moiré pattern. The period of the moiré,
$L=\frac{a}{2\sin(\theta/2)}$, dramatically increases by reducing $\theta$. We
describe the TBG within the low energy continuum model considered in Refs. [6,
9, 24, 25], which is meaningful for sufficiently small angles, so that an
approximately commensurate structure can be defined for any twist. The moiré
mini-BZ, resulting from the folding of the two BZs of each monolayer (see
Fig.1(a)), is generated by the two reciprocal lattice vectors:

$\boldsymbol{G}_{1}=\frac{2\pi}{L}\left(\frac{1}{\sqrt{3}},1\right)\text{ and
}\boldsymbol{G}_{2}=\frac{4\pi}{L}\left(-\frac{1}{\sqrt{3}},0\right),$ (1)

shown in green in Fig. 1(b).

Figure 1:  (a) Folding of the BZs of each of the monolayers in bilayer
graphene. The BZ of layer 1 (indicated by the red hexagon) is rotated by
$-\theta/2$, while that of the layer 2 (the blue hexagon) is rotated by
$\theta/2$. The small black hexagons represent the unit cells forming the
reciprocal moiré lattice. In the inset we show $K_{1,2}$ are the Dirac points
of the twisted monolayers, which identify the corners of the mini-BZ [25]. (b)
mini-BZ. $\boldsymbol{G}_{1,2}$ are the two reciprocal lattice vectors. The
blue line shows the high symmetry circuit in the mini-BZ used to compute the
bands shown in the following.

Let $K_{\xi}=\xi 4\pi(1,0)/3a$ be the two Dirac points of the unrotated
monolayer graphene, with $\xi=\pm 1$. For small twists, the coupling between
the $K_{+}$ and $K_{-}$ valleys of the two monolayers can be safely neglected,
as the interlayer hopping has a long wavelength modulation. The fermionic
field operators of the TBG are 4-component Nambu spinors:

$\displaystyle\Psi_{\xi\sigma}=\left(\psi_{\xi\sigma}^{A_{1}},\psi_{\xi\sigma}^{B_{1}},\psi_{\xi\sigma}^{A_{2}},\psi_{\xi\sigma}^{B_{2}}\right)^{T},$
(2)

where $A,B$ and the subscripts $1,2$ denote the sub-lattice and layer indices,
respectively, and $\sigma$ is the spin index. We introduce a relative twist
$\theta$ between the two monolayers by rotating the layer $1$ by $-\theta/2$
and the layer $2$ by $\theta/2$. Without loss of generality, we assume that in
the aligned configuration, at $\theta=0$, the two layers are $AA$-stacked. In
the continuum limit, the effective Hamiltonian of the TBG in a volume $\Omega$
can be generally written as [6, 9, 24, 25]:

$\displaystyle\hat{H}_{TBG}=\sum_{\xi\sigma}\int_{\Omega}\,d^{2}\boldsymbol{r}\Psi_{\xi\sigma}^{\dagger}(\boldsymbol{r})\begin{pmatrix}H_{\xi
1}&U_{\xi}(\boldsymbol{r})\\\ U_{\xi}^{\dagger}(\boldsymbol{r})&H_{\xi
2}\end{pmatrix}\Psi_{\xi\sigma}(\boldsymbol{r}),$ (3)

where

$\displaystyle H_{\xi l}=\xi\hbar v_{F}\left(-i\boldsymbol{\nabla}-\xi
K_{l}\right)\cdot\boldsymbol{\tau}^{\xi}_{\theta_{l}}$ (4)

is the Dirac Hamiltonian for the $\xi$ valley of layer $l$,
$v_{F}=\sqrt{3}ta/(2\hbar)$ is the Fermi velocity, $t$ is the hopping
amplitude between localized $p_{z}$ orbitals at nearest neighbors carbon
atoms, $\theta_{1,2}=\mp\theta/2$,
$\boldsymbol{\tau}^{\xi}_{\theta_{l}}=e^{i\xi\tau_{z}\theta_{l}/2}\left(\tau_{x},\xi\tau_{y}\right)e^{-i\xi\tau_{z}\theta_{l}/2}$,
$\tau_{i}$ are the Pauli matrices, and
$K_{l}=4\pi\left(\cos\theta_{l},\sin\theta_{l}\right)/(3a)$ are the Dirac
points of the two twisted monolayers corresponding to the $\xi=+$ valley,
which identify the corners of the mini-BZ shown in Fig. 1(a).
$U_{\xi}(\boldsymbol{r})$ is the inter layer potential, which is a periodic
function in the moiré unit cell. In the limit of small angles, its leading
harmonic expansion is determined by only three reciprocal lattice vectors [6]:
$U_{\xi}(\boldsymbol{r})=U_{\xi}(0)+U_{\xi}\left(-\boldsymbol{G}_{1}\right)e^{-i\xi\boldsymbol{G}_{1}\cdot\boldsymbol{r}}+U_{\xi}\left(-\boldsymbol{G}_{1}-\boldsymbol{G}_{2}\right)e^{-i\xi\left(\boldsymbol{G}_{1}+\boldsymbol{G}_{2}\right)\cdot\boldsymbol{r}}$,
where the amplitudes $U_{\xi}\left(\boldsymbol{G}\right)$ are given by:

$\displaystyle U_{\xi}(0)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}g_{1}&g_{2}\\\ g_{2}&g_{1}\end{pmatrix},$
$\displaystyle U_{\xi}\left(-\boldsymbol{G}_{1}\right)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}g_{1}&g_{2}e^{-2i\xi\pi/3}\\\
g_{2}e^{2i\xi\pi/3}&g_{1}\end{pmatrix},$ (5) $\displaystyle
U_{\xi}\left(-\boldsymbol{G}_{1}-\boldsymbol{G}_{2}\right)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}g_{1}&g_{2}e^{2i\xi\pi/3}\\\
g_{2}e^{-2i\xi\pi/3}&g_{1}\end{pmatrix}.$

In the following we adopt the parametrization of the TBG given in the Ref.
[25]: $\hbar v_{F}/a=2.1354$eV, $g_{1}=0.0797$eV and $g_{2}=0.0975$eV. The
difference between $g_{1}$ and $g_{2}$, as described in [25], accounts for the
inhomogeneous interlayer distance, which is minimum in the $AB/BA$ regions and
maximum in the $AA$ ones, or it can be seen as a model of a more complete
treatment of lattice relaxation [26]. In that latter reference it is also
explained that one may want to use the full dispersion of the graphene bands,
and it is shown that in general the Fermi surface shifts as a function of
twist angle. If we focus, e.g., on the $\xi=+$ valley, then the Hamiltonian of
Eq. (3) hybridizes states of layer $1$ with momentum $\boldsymbol{k}$ close to
the Dirac point with the states of layer $2$ with momenta
$\boldsymbol{k},\boldsymbol{k}+\boldsymbol{G}_{1},\boldsymbol{k}+\boldsymbol{G}_{1}+\boldsymbol{G}_{2}$.

In the absence of interactions, the band structure and the DOS per moiré unit
cell of the mini-bands at charge neutrality (CN) are shown in Fig. 2, for
$\theta=1.05^{\circ}$. The green dashed line identifies the Fermi energy. The
bands are computed along the high symmetry circuit of the BZ denoted by the
blue arrows in Fig. 1(b). Here we are showing only one valley, the opposite
one being related by the time-reversal symmetry, upon inverting
$\boldsymbol{k}$ to $-\boldsymbol{k}$. $A_{C}=L^{2}\sqrt{3}/2$ is the area of
the moiré unit cell and the DOS is normalized to 8, accounting for two bands
and four spin/valley flavors. As deeply studied in the past literature [6, 9],
these bands are characterized by a very narrow bandwidth, $\lesssim 5$ meV,
and by an almost vanishing Fermi velocity as compared to that of monolayer
graphene, thus pinning the DOS at the two van Hove singularities in Fig. 2(b).

Figure 2:  (a) mini-bands of the non-interacting TBG at CN, obtained for the
twist angle $\theta=1.05^{\circ}$ and computed along the high symmetry circuit
of the BZ denoted by the blue arrows in Fig. 1(b). The green dashed line
identifies the Fermi energy. (b): DOS per moiré unit cell, normalized to 8.

Next we introduce the Coulomb interaction, as described by the Hamiltonian

$\displaystyle\hat{H}_{C}=\frac{1}{2}\int_{\Omega}\,d^{2}\boldsymbol{r}d^{2}\boldsymbol{r}^{\prime}\delta\hat{\rho}(\boldsymbol{r})v_{C}(\boldsymbol{r}-\boldsymbol{r}^{\prime})\delta\hat{\rho}(\boldsymbol{r}^{\prime}),$
(6)

where
$\delta\hat{\rho}(\boldsymbol{r})\equiv\hat{\rho}(\boldsymbol{r})-\rho_{CN}(\boldsymbol{r})$
is the quantum operator associated to the density fluctuations,
$\hat{\rho}(\boldsymbol{r})=\sum_{\mu}\Psi^{\dagger}_{\mu}(\boldsymbol{r})\Psi_{\mu}(\boldsymbol{r})$
is the density operator, $\mu=(\xi,\sigma)$ are the generalized valley and
spin index, $\rho_{CN}(\boldsymbol{r})$ is the average density corresponding
to the non-interacting TBG at CN, and $v_{C}(\boldsymbol{r})$ is the Coulomb
potential. In the following we assume that the Coulomb interaction is actually
screened by a double metallic gate, and thus described by the Fourier
transform

$v_{C}(\boldsymbol{q})\equiv\int\,d^{2}\boldsymbol{r}v_{C}(\boldsymbol{r})e^{-i\boldsymbol{q}\cdot\boldsymbol{r}}=\frac{2\pi
e^{2}}{\epsilon|\boldsymbol{q}|}\tanh\left(|\boldsymbol{q}|d\right),$ (7)

where $e$ is the electron charge, $\epsilon$ the dielectric constant of the
environment and $d$ the distance of the sample from the gate.

The calculations to be discussed have been carried out in momentum space,
using the continuum model. The Fourier transform of the density operator is

$\displaystyle\hat{\rho}_{\boldsymbol{q}}$
$\displaystyle=e^{i\boldsymbol{q}\boldsymbol{r}}{\cal I}_{s}{\cal
I}_{\tau}{\cal I}_{\sigma}{\cal I}_{l}$ (8)

where $\boldsymbol{r}$ is the position operator in the continuum model, and
${\cal I}_{i}$ is a unit $2\times 2$ matrix in the spin, valley, sublattice,
and layer degrees of freedom. The mean field analysis described below can be
expressed in terms of expectation values of the operator in Eq. (8). Of
particular relevance are the expectation values (or form factors) of the type
$\langle\boldsymbol{k},m|\hat{\rho}_{\boldsymbol{q}}|\boldsymbol{k}+\boldsymbol{q},n\rangle$,
where the indices $m$ and $n$ stand for the central bands. Results in the
following sections have been calculated using the parameters $\epsilon=10$ and
$d=40\leavevmode\nobreak\ \text{nm}$, which are realistic values for most
experiments. At mean-field level, the Hamiltonian $\hat{H}_{C}$ is replaced by

$\displaystyle\hat{H}_{C}\rightarrow\hat{H}^{MF}_{C}=\hat{H}_{H}+\hat{H}_{F}+E_{0},$
(9)

where

$\displaystyle\hat{H}_{H}=\sum_{i\mu}\int_{\Omega}\,d^{2}\boldsymbol{r}\psi^{i,\dagger}_{\mu}(\boldsymbol{r})\psi^{i}_{\mu}(\boldsymbol{r})V_{H}(\boldsymbol{r})$
(10a) is the Hartree Hamiltonian and
$V_{H}(\boldsymbol{r})=\int_{\Omega}\,d^{2}\boldsymbol{r}^{\prime}v_{C}(\boldsymbol{r}-\boldsymbol{r}^{\prime})\left\langle\delta\hat{r}(\boldsymbol{r}^{\prime})\right\rangle$
is the local Hartree potential;
$\displaystyle\hat{H}_{F}=\sum_{ij\mu}\int_{\Omega}\,d^{2}\boldsymbol{r}\,d^{2}\boldsymbol{r}^{\prime}\psi^{i,\dagger}_{\mu}(\boldsymbol{r})V^{ij}_{F,\mu}(\boldsymbol{r},\boldsymbol{r}^{\prime})\psi^{j}_{\mu}(\boldsymbol{r}^{\prime})$
(10b) is the Fock Hamiltonian and
$V^{ij}_{F,\mu}(\boldsymbol{r},\boldsymbol{r}^{\prime})=-\left\langle\psi^{j,\dagger}_{\mu}(\boldsymbol{r}^{\prime})\psi^{i}_{\mu}(\boldsymbol{r})\right\rangle
v_{C}(\boldsymbol{r}-\boldsymbol{r}^{\prime})$ is the non-local Fock
potential.

The Hartree potential is absent in single layer graphene, and in aligned
bilayers, because the charge distribution is homogeneous, and compensated by
the external gate. On the other hand, the exchange (Fock) potential is present
in monolayer and bilayer graphene. Its effect is by no means negligible, as it
causes the renormalization of the Fermi velocity [27]. In very clean graphene
samples, and at charge neutrality, the Fermi velocity diverges logarithmically
as function of a cut-off defined by the ratio between the bandwidth and a low
energy energy cutoff set by the Fermi energy, or by the fluctuations in the
local charge due to disorder. Part of the effects of this exchange potential
is already included in the effective parameters which describe the bands of
graphene. When calculating the exchange potential in twisted bilayer graphene
this renormalization needs to be taken into account. Different procedures have
been proposed [28, 29], such as the subtraction of the exchange potential in a
monolayer, or the subtraction of the exchange potential due to the non
interacting central bands in twisted bilayer graphene. The exchange potential
per graphene unit cell in one monolayer scales as $\Sigma_{x}\propto
e^{2}\Lambda$, where $\Lambda\sim 1/a$ is a momentum cut-off comparable to the
lattice unit cell. Hence, the effect of the central bands in twisted bilayer
graphene can be expected to be of order $\sim e^{2}/L$ per moiré unit cell of
side $L$. In the following, we will consider that the contribution of the
central bands to the exchange energy in monolayer graphene can be neglected,
and that no subtraction is required. We have,

$\displaystyle E_{0}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left[\left\langle\hat{H}_{H}+\hat{H}_{F}\right\rangle\
+\int_{\Omega}\,d^{2}\boldsymbol{r}\rho_{CN}(\boldsymbol{r})V_{H}(\boldsymbol{r})\right]$
(10c)

as the zero point energy, which is required to avoid double counting of the
total energy at mean-field level. The mean-field Hamiltonian for the
interacting TBG is then:

$\displaystyle\hat{H}^{MF}=\hat{H}_{TBG}+\hat{H}_{C}^{MF},$ (11)

which we diagonalize self-consistently, by computing the quantum averages of
the Eqs. (10) by means of $\hat{H}^{MF}$ and iterating until convergence. It
is worth noting that this procedure is equivalent to minimize the GS energy of
$\hat{H}^{MF}$. In order to diagonalize $\hat{H}^{MF}$, we exploit the Bloch’s
theorem, by expressing the eigenfunctions in the basis of Bloch’s plane waves
defined on the moiré:

$\displaystyle\ket{\boldsymbol{k},\alpha,\mu}=\sum_{\boldsymbol{G}i}\phi_{\boldsymbol{k}+\boldsymbol{G},\alpha,\mu,i}\ket{\boldsymbol{k}+\boldsymbol{G},\mu,i},$
(12)

where $\boldsymbol{k}\in$mBZ, the $\boldsymbol{G}$’s are reciprocal lattice
vectors, $\alpha$ is the band index and
$\phi_{\boldsymbol{k}+\boldsymbol{G},\alpha,\mu,i}$ are numerical eigenvectors
normalized according to:
$\sum_{i\boldsymbol{G}}\phi^{*}_{\boldsymbol{k}+\boldsymbol{G},\alpha,\mu,i}\phi_{\boldsymbol{k}+\boldsymbol{G},\alpha^{\prime},\mu,i}=\delta_{\alpha\alpha^{\prime}}$.
Upon using the Eq. (12) to evaluate the Hartree and Fock potentials, the
matrix elements of the Eqs. (10) can be written in the Bloch’s basis as:

$\displaystyle\bra{\boldsymbol{k}+\boldsymbol{G},\mu,i}\hat{H}_{H}\ket{\boldsymbol{k}+\boldsymbol{G}^{\prime},\mu^{\prime},i^{\prime}}=\delta_{ii^{\prime}}\delta_{\mu\mu^{\prime}}\frac{v_{C}(\boldsymbol{G}-\boldsymbol{G}^{\prime})}{\Omega}\times$
$\displaystyle\times\sum_{\boldsymbol{k}^{\prime}\boldsymbol{G}^{\prime\prime}}\sum_{\alpha\mu^{\prime\prime}i^{\prime\prime}}\phi_{\boldsymbol{k}^{\prime}+\boldsymbol{G}^{\prime\prime}+\boldsymbol{G},\alpha,\mu^{\prime\prime},i^{\prime\prime}}\phi^{*}_{\boldsymbol{k}^{\prime}+\boldsymbol{G}^{\prime\prime}+\boldsymbol{G}^{\prime},\alpha,\mu^{\prime\prime},i^{\prime\prime}}\equiv\delta_{ii^{\prime}}\delta_{\mu\mu^{\prime}}V_{H}\left(\boldsymbol{G}-\boldsymbol{G}^{\prime}\right),$
(13a) where the sum over the band index, $\alpha$, runs over all the occupied
states counted from CN, and $V_{H}\left(\boldsymbol{G}\right)$ is noting but
the Fourier transform of the Hartree potential,
$V_{H}\left(\boldsymbol{r}\right)$, evaluated in $\boldsymbol{G}$.
$\displaystyle\bra{\boldsymbol{k}+\boldsymbol{G},\mu,i}\hat{H}_{F}\ket{\boldsymbol{k}+\boldsymbol{G}^{\prime},\mu^{\prime},i^{\prime}}$
$\displaystyle=$
$\displaystyle-\delta_{\mu\mu^{\prime}}\sum_{\boldsymbol{k}^{\prime}\boldsymbol{G}^{\prime\prime}\alpha}\frac{v_{C}(\boldsymbol{k}-\boldsymbol{k}^{\prime}-\boldsymbol{G}^{\prime\prime})}{\Omega}\phi_{\boldsymbol{k}^{\prime}+\boldsymbol{G}^{\prime\prime}+\boldsymbol{G},\alpha,\mu,i}\phi^{*}_{\boldsymbol{k}^{\prime}+\boldsymbol{G}^{\prime\prime}+\boldsymbol{G}^{\prime},\alpha,\mu^{\prime},i^{\prime}}$
$\displaystyle\equiv$
$\displaystyle\delta_{\mu\mu^{\prime}}\Sigma_{F}^{\mu}(\boldsymbol{k}+\boldsymbol{G},i;\boldsymbol{k}+\boldsymbol{G}^{\prime},i^{\prime}),$
(13b)

where $\alpha$ runs over all the occupied states above a given threshold. In
the present context, we set this threshold to the lowest energy of the mini-
bands in the middle of the spectrum, meaning that we are neglecting the
contribution of the bulk bands. However, including other bands might affect
quantitatively the results. Because the Eqs. (13) express the matrix elements
in terms of the energy levels and of the corresponding eigenfunctions, $\phi$,
they completely define the self-consistent problem.

The main contributions of the long-range interaction are expected to come from
small momenta. Therefore, we only consider the matrix elements of the Hartree
potential, Eq. (13), with $\boldsymbol{G}-\boldsymbol{G}^{\prime}$ belonging
to the first star of reciprocal lattice vectors:
$\pm\boldsymbol{G}_{1},\pm\boldsymbol{G}_{2},\pm\left(\boldsymbol{G}_{1}+\boldsymbol{G}_{2}\right)$.
Concerning the matrix elements of the Fock potential, Eq. (13), for each
external momentum $\boldsymbol{k}$ we truncate the sum over
$\boldsymbol{k}^{\prime}$ and $\boldsymbol{G}^{\prime\prime}$ so that:
$\boldsymbol{k}-\boldsymbol{k}^{\prime}-\boldsymbol{G}^{\prime\prime}$ belongs
to the BZ. We checked that including larger momenta affects negligibly the
results.

Finally, the energy of the GS, as following from the Eqs. (9)-(11), is given
by:

$\displaystyle
E_{GS}=\sum_{\boldsymbol{k}\alpha\mu}\varepsilon_{\boldsymbol{k}\alpha\mu}+E_{0},$
(14)

where $\varepsilon_{\boldsymbol{k}\alpha\mu}$ are the single-particle energies
and the sum over $\alpha$ runs over all the occupied states.

Figure 3: Charge density at the high symmetry points within the mBZ for the
lower band in Fig. 2. White hexagon is the real space unit cell.

##  III Evolution of the band structure as function of filling

Figure 4 shows the band structure and DOS of the mini-bands of the interacting
TBG at $\theta=1.05^{\circ}$, for the non-polarized GS at positive filling of
the conduction bands: $\nu=0,1,2$, corresponding to CN, one and two electrons
per moiré unit cell, respectively. The green dashed line identifies the Fermi
energy. Each band is four-fold degenerate, so that the spin/valley flavors are
equally occupied. At $\nu=0$, only the lower band is filled, the
$\mathcal{C}_{2}$ symmetry is broken and the Fermi surface (FS) is fully
gapped. For this choice of parameters, the width of the gap is $\sim 5$meV,
comparable to the overall bandwidth. At $\nu=1$, one quarter of the upper band
is filled, the $\mathcal{C}_{2}$ symmetry is still broken, but the FS exists.
At $\nu=2$, the $\mathcal{C}_{2}$ symmetry is completely restored.

Figure 4: Band structure and DOS corresponding to the non-polarized GS,
obtained within the Hartree-Fock approximation at the twist angle
$\theta=1.05^{\circ}$ and filling: $\nu=0(a),1(b),2(c)$. The green dashed line
identifies the Fermi energy.

The evolution of the band structure corresponding to the non-polarized
solution is shown in Fig. 5, for $-2.5\leq\nu\leq 2.5$. As it can be seen,
solutions breaking the $\mathcal{C}_{2}$ symmetry are not stable for $\nu>1.5$
and $\nu\leq-1.5$. Two features of the band structure are worth to be further
noticed: i) the lack of particle-hole symmetry, so that the bands are not
symmetric upon inverting $\nu$ to $-\nu$; ii) the bands are rigid at the
$\Gamma$ point of the BZ. This is a consequence of the fact that the charge
density of the TBG evaluated at the $\Gamma$ point of the Brillouin Zone is
almost homogeneous as compared to the other high symmetry points of the BZ
[30], as shown in Fig. 3. The resulting dependence of the band shape on
filling factor has been extensively discussed [12, 31, 32, 33, 34, 35, 28, 36,
37]. This significant dependence of the band shape on electron filling implies
the existence of a special effective electron-electron interaction, electron
assisted hopping [38, 39, 40, 41, 42].

The dependence of the Hartree-Fock potential on electron occupancy arises from
the variations in charge distribution in different regions of the Brillouin
Zone [30]. These variations are determined by form factors of the type:

$\displaystyle{\cal
M}_{\boldsymbol{G}}(\boldsymbol{k},\boldsymbol{k}+\boldsymbol{q})$
$\displaystyle=\left\langle\boldsymbol{k},m\left|\hat{\rho}_{\boldsymbol{q}+\boldsymbol{G}}\right|\boldsymbol{k}+\boldsymbol{q},n\right\rangle=\int_{\Omega}d^{2}\boldsymbol{r}u_{\boldsymbol{k},m}^{*}(\boldsymbol{r})e^{i\boldsymbol{G}\boldsymbol{r}}u_{\boldsymbol{k}+\boldsymbol{q},n}(\boldsymbol{r})$
(15)

where $\boldsymbol{G}$ is a reciprocal lattice vector, $\boldsymbol{k}$ and
$\boldsymbol{k}+\boldsymbol{q}$ are in the moiré Brillouin Zone, $m$ and $n$
are band indices, $u_{\boldsymbol{k},m}(\boldsymbol{r})$ is the periodic part
of the wavefunction associated to $|\boldsymbol{k},n\rangle$, and $\Omega$ is
the moiré unit cell. These from factors encode the complexity of the band
wavefunctions of twisted bilayer graphene. There are some simplifying
approximations made in the literature. These are very convenient because they
allow for simple analytical calculations, like the ”flat metric condition”
used in [43, 44], however, the neglect the dependence of the form factors in
Eq. (15) on the position $\boldsymbol{k}$ in the Brillouin Zone. Since the
variations of the form factors over the BZ are substantial, it is not
immediately clear how realistic this approximation is.

Figure 5: Evolution of the band structure in the Hartree-Fock approximation
upon varying the filling, for $-2.5\leq\nu\leq 2.5$. The green dashed line
identifies the Fermi energy.

Figure 6:  Band structure obtained at integer filling, $-3\leq\nu\leq 3$, for
three different polarized solutions: $2+2$ (a), $1+3$ (b) and $3+1$ (c). The
blue and magenta lines show the bands of higher and lower occupancy,
respectively. The green dashed line in each panel identifies the Fermi energy.
Figure 7:  Fermi energy as a function of the filling for the $1+3$ solution.
The continuous black line shows the value of $E_{F}$ as obtained from the
highest energy occupied states, while the dashed cyan line represents an
alternative definition of $E_{F}$, the derivative of the GS energy with
respect to $\nu$.

##  IV Broken symmetry solutions.

###  IV.1 Polarized solutions.

We consider first valley and/or spin polarized phases, motivated by the
”cascade” of phases found experimentally [45]. We define these phases using
two integers, $n_{+}$ and $n_{-}$, so that there are $n_{+}$ bands with high
occupancy, and $n_{-}$ bands with low occupancy. For a total occupancy $\nu$,
such that $-4\leq\nu\leq 4$, the occupancies of the two types of bands are
listed in table [2].

| $2+2$ | $3+1$ | $1+3$
---|---|---|---
$\displaystyle n_{-}$ | $\displaystyle 4+\frac{\nu-|\nu|}{2}$ | $\displaystyle 3+\frac{\nu-|\nu+2|}{2}$ | $\displaystyle 5+\frac{\nu-|\nu-2|}{2}$
$\displaystyle n_{-}$ | $\displaystyle\frac{\nu+|\nu|}{2}$ | $\displaystyle 1+\frac{\nu+|\nu+2|}{2}$ | $\displaystyle-1+\frac{\nu+|\nu-2|}{2}$
Table 2: Filling of the low occupancy bands, $n_{-}$, and the high occupancy
bands, $n_{+}$, as function of the total filling, $-4\leq n_{-}\leq 4$, see
Fig. [6].

Figure 6 shows the evolution of the band structure of the polarized solutions:
$2+2$ (a), $1+3$ (b) and $3+1$ (c), for integer fillings and $-3\leq\nu\leq
3$. The high and low occupancy sets of bands are represented by the blue and
magenta lines, respectively. Solutions breaking the $\mathcal{C}_{2}$ symmetry
occur at $\nu=\pm 2$ for the $2+2$ solution, $\nu=1$ for the $1+3$ solution
and $\nu=-1,3$ for the $3+1$ solution. Interestingly, the solution $1+3$ does
not break $\mathcal{C}_{2}$ at $\nu=-3$, in contrast to what is expected. We
argue that here the $\mathcal{C}_{2}$ symmetry breaking is prevented by the
small value of the interaction. Furthermore, it is worth noting that, in the
$1+3$ solution at $\nu=-2$, the empty low occupancy bands mostly stay below
the Fermi level. At this filling the Fermi energy, which we define as the
highest energy of the occupied states, decreases upon increasing the filling,
implying a negative compressibility. This behavior can be better seen in the
Fig. 7, showing $E_{F}$ as a function of $\nu$ and displaying a jump-like
discontinuity at $\nu=-2$. The cyan dashed line in the Fig. 7 represents the
value of $E_{F}$ as obtained from the derivative of the GS energy with respect
to $\nu$, which indeed matches quite well the curve of $E_{F}$ computed as
described above. However, this anomalous behavior is actually not very
meaningful in the present context, as the solution $1+3$ turns out to not be
stable close to $\nu=-2$, as emphasized by Figs. 5 and Fig. 6.

###  IV.2 Solutions at half filling.

The competition between phases at half filling has extensively been discussed
[28, 46, 29, 37]. Hartree-Fock calculations at the bands near half filling
have already been presented in Fig. [4] and Fig. [5]. At half filling it is
reasonable to expect that the Hartree potential vanishes (note that the
continuum model allows us only to compute differences in the charge
distribution at different fillings). This assumption is consistent with
lattice based tight binding calculations which include the entire $p_{z}$ band
[35]. It is worth noting that modifications of the moiré structure, such as
the one induced by heterostrains, can alter significantly the effect of
interactions [47].

We assume first that the leading mean field decoupling of the interaction is
of the type

$\displaystyle{\cal
H}_{MF}^{int}=\sum_{\boldsymbol{q}}V_{\boldsymbol{q}}\langle\hat{\rho}_{\boldsymbol{q}}\rangle\hat{\rho}_{\boldsymbol{q}}+h.c.$
(16)

where $\hat{\rho}_{\boldsymbol{q}}$ is the density operator, eq.(8). Within
this approximation, the only coupling between electrons of different flavors,
valley and spin, is through the Hartree potential, as the density operator,
Eq. (8) is diagonal in valley space.

The lack of a Hartree potential implies that the resulting Fock hamiltonian,
within the decoupling scheme in Eq.(16) is split into four independent
hamiltonians, one for each spin and valley. Each of these hamiltonians admits
solutions with the full symmetries of the non interacting Hamiltonian, and, at
least one solution where the equivalence between the two sublattices,
described by a ${\cal C}_{2}{\cal T}$ symmetry is broken. The latter solution
shows a gap between the valence and the conduction bands, and it has a lower
energy when the valence band is fully occupied and the conduction band is
empty. As in this phase the charge has a higher weight in one sublattice, we
can describe the corresponding Slater determinant in terms of the spin,
valley, and majority sublattice, $|Sp,V,S\rangle$, where
$Sp=\uparrow,\downarrow;V=K,K^{\prime};S=A,B$. In each of the possible
$2^{3}=8$ wavefunctions the occupied valence band is associated to a valley,
so that it carries a Chern number ${\cal C}=\pm 1$, depending on the
combination of valley and sublattice, $\tau_{z}\sigma_{z}$.

A possible ground state at half filling is:

$\displaystyle|GS\rangle_{\sigma}$
$\displaystyle\equiv|\uparrow,K,S_{\uparrow,K}\rangle\otimes|\downarrow,K,S_{\downarrow,K}\rangle\otimes|\uparrow,K^{\prime},S_{\uparrow,K^{\prime}}\rangle\otimes|\downarrow,K^{\prime},S_{\downarrow,K^{\prime}}\rangle$
(17)

where the sublattice labels $S_{Sp,V}$ can take the values $A,B$. As this
sublattice index can take two values in each component, there are 16 possible
wavefunctions of the type shown in Eq. (17). There are two wavefunctions with
Chern number ${\cal C}=\pm 4$, four wavefunctions with ${\cal C}=\pm 2$, and
six wavefunctions with ${\cal C}=0$.

Other Slater determinants which maximize the overlap between one electron
wavefunctions are superpositions of states with the same valley and spin
indices. These determinants are globally spin and valley polarized.

$\displaystyle|GS\rangle_{P}$
$\displaystyle\equiv|Sp_{1},V_{1},v\rangle\otimes|Sp_{1},V_{1},c\rangle\otimes|Sp_{2},V_{2},v\rangle\otimes|Sp_{2},V_{2},c\rangle$
(18)

where $Sp_{1},Sp_{2}$ are spin indices, $V_{1},V_{2}$ are valley indices, and
$v,c$ stand for valence and conduction bands. In the limit when the bands are
infinitely flat [48, 49] this valley/spin polarized wavefunction is
degenerate, within the Hartree-Fock approximation, with the sublattice
polarized wavefunction in Eq. (17).

The wavefunction in Eq. (17) shows four occupied valence bands, separated by a
gap from four unoccupied conduction bands, (the problem reduces to four copies
of the bands shown in Fig.[8]). The valley and spin polarized wavefunction in
Eq. (18) shows two occupied valence bands and two occupied conduction bands.
In this case, there is no energy gain in opening a gap between the valence and
conduction bands.We can expect, in both cases, that the overlap between one
electron wavefunctions is higher between states in the same band (valence and
conduction) than between states in different bands. If that is the case, the
wavefunction in Eq. (17) will have a lower energy, as no conduction band is
occupied.

At this stage, the Hartree-Fock approximation leaves undetermined the relative
phase between wavefunctions in different valleys and different sublattices.
Solutions with well defined phases between Slater determinants in different
valleys arise if the mean field hamiltonian allows for inter-valley matrix
elements, not included in the decoupling scheme described in Eq. (16).

Such terms can favor a inter-valley coherent (IVC) combination of the type
[29]:

$\displaystyle|GS\rangle_{IVC}$
$\displaystyle\equiv\left(|\uparrow,K,A\rangle+ie^{i\phi}|\uparrow,K^{\prime},B\rangle\right)\otimes\left(|\uparrow,K,B\rangle-
ie^{-i\phi}|\uparrow,K^{\prime},A\rangle\right)\otimes$
$\displaystyle\otimes\left(|\downarrow,K,A\rangle+ie^{i\phi}|\downarrow,K^{\prime},B\rangle\right)\otimes\left(|\downarrow,K,B\rangle-
ie^{-i\phi}|\downarrow,K^{\prime},A\rangle\right)$ (19)

The existence of a gap in the one particle spectrum of the independent valley
solutions shown in Fig.[8] implies that inter-valley terms will only appear
above a threshold determined by the ratio between the inter-valley
interactions and the gap. In order to check for the emergence of inter-valley
terms in the Hartree-Fock hamiltonian, we have carried out calculations
allowing for off-diagonal valley terms in the Hartree-Fock hamiltonian,
projecting the hamiltonian onto the two central bands. The effect of
intervalley biases in the initial wavefunctions tends to vanish in the
iteraction towards self consistency, and the bands with and without initial
intervalley bias agree within numerical accuracy. These results suggest that
the valleys remain decoupled.

Interactions at the atomic scale, will favor a phase in which the two
sublattices have equal occupation, over the phase where all the charge resides
in the same sublattice [50, 51, 52, 53]. Finally, the zeroth point energy of
excitations above the Hartree-Fock solutions discussed here can also favor
specific combinations of wavefunctions of the type shown in Eq. (17) with
different Chern numbers [22], see Eq. (19).

It is worth noting that, even in the absence of the Hartree potential, the
exchange term induces significant deformations in the band dispersion. The
electronic bands associated to the state in Eq. (17) for initially dispersive
and for the non dispersive bands resulting from the chiral model [48, 49] are
shown in Fig. 8.

Figure 8:  Hartree-Fock bands at a magic angle and at half filling, for a
state of the form Eq. (17). We show the self-consistent solutions for
initially a) dispersive bands, as in Fig. 2, and in b) for initially non-
dispersive bands (chiral model at a magic angle) [48, 49].

##  V Conclusions.

We have discussed the electronic properties of interacting electrons in
twisted bilayer graphene obtained using the mean field Hartree-Fock
approximation. This is usually the initial method to study interactions in
theoretical condensed matter physics, and it is also a reasonable starting
point to analyze effects beyond mean field theory. The Hartree-Fock method has
been applied to the study of the long range electrostatic interactions. Simple
scaling arguments suggest that this is the leading interaction between
electrons in twisted bilayer graphene.

Despite the simplicity of the method, the Hartree-Fock analysis of
interactions in twisted bilayer graphene shows a number of features not
commonly found in most materials studied in condensed matter physics. Among
them, we can emphasize:

  * •

A significant distortion of the band shapes as function of the electron
filling. This effect is mostly driven by the Hartree term, although it is also
significant at half filling, solely due to the exchange potential.

This filling dependence makes twisted bilayer graphene a system where a good
case can be made for the existence of electron assisted hopping.

  * •

The existence of a gapped phase at half filling, driven by the exchange term.

This phase is sublattice polarized in each flavor (spin and valley) sector.
The Hartree-Fock allows for many combinations of these separate sectors.

  * •

Away from half filling, the Hartree-Fock approximation leads to spin and/or
valley polarized phases. As at half filling, the Hartree-Fock approximation
respects the spin and valley degeneracy, so that the details of the
polarization cannot be resolved.

  * •

The energy differences between the polarized phases away from half filling are
small, of the order of a few meV per unit cell. At these scales, other
interactions, like the intra-atomic Hubbard repulsion, are likely to play a
role, and select the most stable ground state.

## Acknowledgements

This work was supported by funding from the European Commission, under the
Graphene Flagship, Core 3, grant number 881603, and by the grants NMAT2D
(Comunidad de Madrid, Spain), SprQuMat and SEV-2016-0686, (Ministerio de
Ciencia e Innovación, Spain). NRW is supported by STFC grant ST/P004423/1. We
are thankful to H. Ochoa ant to M. A. Cazalilla for helpful discussions, and
to M. I. Katsnelson for pointing out to the references [38, 39, 40].

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018a), arXiv:1802.00553 .
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Nature Materials 19, 1265 (2020).
  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Nature Physics 16, 725 (2020).
  * Wang  _et al._ [2020] J. Wang, X. Mu, L. Wang, and M. Sun, Materials Today Phys. 9, 100099 (2020).
  * Lopes dos Santos  _et al._ [2007] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Mele [2010] E. J. Mele, Phys. Rev. B 81, 161405 (2010).
  * Mele [2011] E. J. Mele, Phys. Rev. B 84, 235439 (2011).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Suárez Morell  _et al._ [2010] E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82, 121407 (2010).
  * Trambly de Laissardière  _et al._ [2010] G. Trambly de Laissardière, D. Mayou, and L. Magaud, Nano Lett. 10, 804 (2010).
  * Guinea and Walet [2018] F. Guinea and N. R. Walet, Proceedings of the National Academy of Sciences 115, 13174 (2018).
  * Kotov  _et al._ [2012] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro Neto, Rev. Mod. Phys. 84, 1067 (2012).
  * Goerbig [2011] M. O. Goerbig, Rev. Mod. Phys. 83, 1193 (2011).
  * Lewandowski and Levitov [2019] C. Lewandowski and L. Levitov, Proceedings of the National Academy of Sciences 116, 20869 (2019), https://www.pnas.org/content/116/42/20869.full.pdf .
  * Kozii  _et al._ [2019] V. Kozii, H. Isobe, J. W. F. Venderbos, and L. Fu, Phys. Rev. B 99, 144507 (2019).
  * Wu and Das Sarma [2020] F. Wu and S. Das Sarma, Phys. Rev. Lett. 124, 046403 (2020).
  * Sharma  _et al._ [2020] G. Sharma, M. Trushin, O. P. Sushkov, G. Vignale, and S. Adam, Phys. Rev. Research 2, 022040 (2020).
  * Khalaf  _et al._ [2020] E. Khalaf, N. Bultinck, A. Vishwanath, and M. P. Zaletel, arXiv e-prints , arXiv:2009.14827 (2020), arXiv:2009.14827 [cond-mat.str-el] .
  * Kumar  _et al._ [2020] A. Kumar, M. Xie, and A. H. MacDonald, arXiv e-prints , arXiv:2010.05946 (2020), arXiv:2010.05946 [cond-mat.str-el] .
  * Kozii  _et al._ [2020] V. Kozii, M. P. Zaletel, and N. Bultinck, arXiv e-prints , arXiv:2005.12961 (2020), arXiv:2005.12961 [cond-mat.str-el] .
  * Khalaf  _et al._ [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Sci. Adv. 7, eabf5299 (2021).
  * Cea and Guinea [2021] T. Cea and F. Guinea, Proc. Nat. Acad. Sci. (USA) (2021).
  * Lopes dos Santos  _et al._ [2012] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8, 031087 (2018).
  * Guinea and Walet [2019] F. Guinea and N. R. Walet, Phys. Rev. B 99, 205134 (2019).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Phys. Rev. Lett. 124, 097601 (2020).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Phys. Rev. X 10, 031034 (2020).
  * Rademaker and Mellado [2018] L. Rademaker and P. Mellado, Phys. Rev. B 98, 235158 (2018).
  * Cea  _et al._ [2019] T. Cea, N. R. Walet, and F. Guinea, Phys. Rev. B 100, 205113 (2019).
  * Rademaker  _et al._ [2019] L. Rademaker, D. A. Abanin, and P. Mellado, Phys. Rev. B 100, 205114 (2019).
  * Calderón and Bascones [2020] M. J. Calderón and E. Bascones, Phys. Rev. B 102, 155149 (2020).
  * Zhang  _et al._ [2020] Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Phys. Rev. B 102, 035136 (2020).
  * Goodwin  _et al._ [2020] Z. A. H. Goodwin, V. Vitale, X. Liang, A. A. Mostofi, and J. Lischner, Electronic Structure 2, 034001 (2020).
  * Liu and Dai [2021] J. Liu and X. Dai, Phys. Rev. B 103, 035427 (2021).
  * Liu  _et al._ [2021] S. Liu, E. Khalaf, J. Y. Lee, and A. Vishwanath, Phys. Rev. Research 3, 013033 (2021).
  * Schubin and Vonchowsky [1934] S. Schubin and S. Vonchowsky, Proc. R. Soc. Lond. A 145, 159 (1934).
  * Vonsovsky and Katsnelson [1979a] S. V. Vonsovsky and M. I. Katsnelson, Journ. Phys. C: Solid St. Phys. , 2043 (1979a).
  * Vonsovsky and Katsnelson [1979b] S. V. Vonsovsky and M. I. Katsnelson, Journ. Phys. C: Solid St. Phys. 12, 2055 (1979b).
  * Hirsch and Marsiglio [1990] J. E. Hirsch and F. Marsiglio, Phys. Rev. B 41, 2049 (1990).
  * Marsiglio and Hirsch [1990] F. Marsiglio and J. E. Hirsch, Phys. Rev. B 41, 6435 (1990).
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Phys. Rev. B 103, 205411 (2021).
  * Song  _et al._ [2021] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Phys. Rev. B 103, 205412 (2021).
  * Wong  _et al._ [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 582, 198 (2020).
  * Cea and Guinea [2020] T. Cea and F. Guinea, Phys. Rev. B 102, 045107 (2020).
  * Ochoa [2020] H. Ochoa, Phys. Rev. B 102, 201107 (2020).
  * San-Jose  _et al._ [2012] P. San-Jose, J. González, and F. Guinea, Phys. Rev. Lett. 108, 216802 (2012).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Phys. Rev. Lett. 122, 106405 (2019).
  * Alicea and Fisher [2006] J. Alicea and M. P. A. Fisher, Phys. Rev. B 74, 075422 (2006).
  * Gonzalez-Arraga  _et al._ [2017] L. A. Gonzalez-Arraga, J. L. Lado, F. Guinea, and P. San-Jose, Phys. Rev. Lett. 119, 107201 (2017).
  * Sboychakov  _et al._ [2019] A. O. Sboychakov, A. V. Rozhkov, A. L. Rakhmanov, and F. Nori, Phys. Rev. B 100, 045111 (2019).
  * Klebl  _et al._ [2021] L. Klebl, Z. A. H. Goodwin, A. A. Mostofi, D. M. Kennes, and J. Lischner, Phys. Rev. B 103, 195127 (2021).


# A New Algorithm for Discrete-Time Parameter Estimation

Yingnan Cui  Massachusetts Institute of Technology  Joseph E. Gaudio  Aurora
Flight Sciences  Anuradha M. Annaswamy  Massachusetts Institute of Technology

###### Abstract

We propose a new discrete-time adaptive algorithm for parameter estimation of
a class of time-varying plants. The main contribution is the inclusion of a
time-varying gain matrix in the adjustment of the parameter estimates. We show
that in the presence of time-varying unknown parameters, the parameter
estimation error converges uniformly to a compact set under conditions of
persistent excitation, with the size of the compact set proportional to the
time-variation of the unknown parameters. Under conditions of finite
excitation, the convergence is asymptotic and non-uniform.

##  1 Introduction

A central task in adaptive control is the compensation for online parametric
uncertainties in the processes of adaptive filtering, prediction and
control[1, 2, 3, 4, 5, 6]. Over the past several decades, several discrete-
time online estimation schemes have been designed, including conditions for
parametric convergence[7, 8, 9]. Most of the methods in these works assume
that the unknown parameters are constant. In this paper, we propose a new
online algorithm for estimating parameters that may be time-varying in a
discrete-time plant.

In many of the adaptive filtering, prediction, and control problems, the
underlying structure consists of a linear regression relation between the
prediction error and parameter error[5, 10, 11, 12, 13, 14]. A typical
gradient-based algorithm, for example, minimizes a quadratic criterion in
terms of the prediction error[4, 6]. A necessary and sufficient condition,
denoted as persistent excitation, is needed to ensure uniform asymptotic
convergence of the estimates to their true values is noted as persistent
excitation in adaptive systems[15, 16, 17]. If the excitation holds over a
finite interval, then the convergence of the estimation errors to zero can be
shown to be asymptotic but non-uniform. This leads to the well-known recursive
least squares (RLS) algorithm, where the learning rates are adjusted through
the use of the outer product of the underlying regressor[4, 6].

One of the main inadequacies of the RLS algorithm is its inability to ensure
estimation in the presence of time-varying parameters, which is due to the
fact that the learning rates themselves may converge to zero and as a result,
unable to track the changes in the parameters. While methods such as
covariance resetting[18], fixed or varying forgetting factors [19] can be
used, they introduce constraints and do not lead to efficient trade-offs
between speed and accuracy. Additionally, they can introduce discontinuities
in the estimation algorithm. The algorithm proposed in this paper removes some
of these inadequacies, and articulates conditions under which exponential
convergence to a compact set can occur. In particular, we focus on the
convergence properties of the proposed algorithms in the presence of both
persistent and finite excitation.

As mentioned earlier, there has been a significant amount of papers and books
that have addressed online parameter estimation algorithms in discrete-time[7,
8, 16, 9, 15, 17, 10, 11, 12, 18, 19]. In many of them, the learning rates are
assumed to be constant. Of these, works such as [18, 19, 20, 9] have focused
on a RLS-type algorithm and address variations thereof. In [18], combined
variable-rate and variable-direction forgetting is adopted to address some of
the inadequacies of the RLS algorithm. In [19], a data-dependent forgetting
rate is used in the update of the covariance matrix. A nonlinear weighting
coefficient is devised in the update of the covariance matrix in [20]. Another
interesting approach has been suggested in [9], which consists of a dynamic
regressor extension and mixing (DREM) idea, and it tries to lift the
limitation of the input signal by a linear, BIBO-stable operator. In all of
these works, the unknown parameters are assumed to be constant. In addition,
the focus in most of these papers is more in showing boundedness of the
parameter estimates, rather than the speed of convergence of these estimates
to their true values.

This paper proposes a new algorithm for fast parameter estimation in discrete-
time plants under varied conditions of excitation. A time-varying gain matrix
is included in the algorithm, but no covariance resetting or forgetting factor
is introduced. It is shown that this algorithm leads to convergence of the
parameter error to a compact set that is proportional to the rate of variation
of the unknown parameters. The convergence is asymptotic under finite
excitation, and becomes uniform and asymptotic under persistent excitation.
This compact set is shown to be directly proportional to the rate of variation
of the unknown parameters.

Section 2 presents the definitions of discrete-time persistent and finite
excitation, and some mathematical preliminaries. The underlying problem and
the proposed algorithm are presented in Section 3. Section 4 presents the
algorithm. Section 5.1 and Section 5.2 show the properties of the information
matrix and time-varying gain we adopted in the algorithm. Section 5.3 provides
a proof of stability and parameter convergence. Section 6 provides a summary
and concluding remarks.

##  2 Preliminaries

We denote $\mathbb{N}^{+}$ as the set of positive integers, and $\|\cdot\|$ as
the Euclidean norm.

######  Definition 1 ([21]).

A bounded function $\phi_{k}\mathrel{\mathop{\mathchar
58\relax}}[k_{0},\infty)\rightarrow\mathbb{R}^{N}$ is persistently exciting
(PE) at a level $\alpha$ if there exists $\Delta T\in\mathbb{N}^{+}$ and
$\alpha>0$ such that

$\sum_{i=k+1}^{k+\Delta T}\phi_{i}\phi_{i}^{\top}\succeq\alpha I,\;\;\forall
k\geq k_{0}.$

######  Definition 2.

A bounded function $\phi_{k}\mathrel{\mathop{\mathchar
58\relax}}[k_{0},\infty)\rightarrow\mathbb{R}^{N}$ is finitely exciting (FE)
over an interval $[k_{1},k_{2}]$ at a level $\alpha$ if there exists $\Delta
T\in\mathbb{N}^{+}$ and $\alpha>0$ such that

$\sum_{i=k_{1}}^{k_{2}}\phi_{i}\phi_{i}^{\top}\succeq\alpha I,\;\;\text{for
some }k_{2}\geq k_{1}+\Delta T$

In both Definition 1 and Definition 2, the degree of excitation is given by
$\alpha$. It can be noted that the PE condition in Definition 1 pertains to a
property over a moving window for all $k\geq k_{0}$, whereas the FE condition
in Definition 2 pertains to a single interval $[k_{1},k_{2}]$.

The following matrix inversion lemma is useful in deriving the results that
follow.

######  Lemma 1.

(Kailath Variant of the Woodbury Identity [22].)

$(A+BC)^{-1}=A^{-1}-A^{-1}B(I+CA^{-1}B)^{-1}CA^{-1}$

##  3 Problem Formulation and The Algorithm

The class of discrete-time models we consider in this paper is of the form[4,
6]

$\begin{split}y_{k}=&-\sum_{i=1}^{n}a_{i,k}^{*}y_{k-i}+\sum_{j=1}^{m}b_{j,k}^{*}u_{k-j-d}\\\
&+\sum_{\ell=1}^{p}c^{*}_{\ell,k}f_{\ell}(y_{k-1},\ldots,y_{k-n},u_{k-1-d},\ldots,u_{k-m-d}),\end{split}$
(1)

where $a_{i,k}^{*}$, $b_{j,k}^{*}$ and $c_{\ell,k}^{*}$ are unknown parameters
that may vary with time and need to be identified, and $d$ is a known time-
delay. The function $f_{\ell}$ is an analytic function of its arguments and is
assumed to be such that the system in (1) is bounded-input-bounded-output
(BIBO) stable. Denote $z_{k-1}=[y_{k-1},\ldots,y_{k-n}]^{\top}$ and
$v_{k-d-1}=[u_{k-1-d},\ldots,u_{k-m-d}]^{\top}$. We rewrite (1) in the form of
a linear regression

$y_{k}=\phi_{k}^{\top}\theta_{k}^{*},$ (2)

where
$\phi_{k}=[z_{k-1}^{\top},v_{k-d-1}^{\top},f_{1}(z_{k-1}^{\top},v_{k-d-1}^{\top}),\ldots,\allowbreak
f_{p}(z_{k-1}^{\top},v_{k-d-1}^{\top})]^{\top}$ is a regressor determined by
exogenous signals and
$\theta_{k}^{*}=[a_{1,k}^{*},\ldots,a_{n,k}^{*},b_{1,k}^{*},\ldots,b_{m,k}^{*},c_{1,k}^{*},\ldots,c_{\ell,k}^{*},]^{\top}$
is the underlying unknown parameter vector. We propose to identify the
parameter $\theta^{*}_{k}$ as $\theta_{k}$ using an estimator

$\hat{y}_{k}=\phi_{k}^{\top}\theta_{k}.$ (3)

This leads to a prediction error

$e_{y,k}=\phi_{k}^{\top}\tilde{\theta}_{k},$ (4)

where $e_{y,k}=\hat{y}_{k}-y_{k}$ is the output prediction error and the
parameter error is $\tilde{\theta}_{k}=\theta_{k}-\theta_{k}^{*}$.

######  Assumption 1.

Defining the time variation in $\theta_{k}^{*}$ as
$\Delta\theta_{k}^{*}=\theta_{k}^{*}-\theta_{k-1}^{*}$, in what follows we
will assume that the time variation $\Delta\theta_{k}^{*}$ is bounded, i.e.,
there exists a $\Delta^{*}\geq 0$ such that

$\|\Delta\theta_{k}^{*}\|\leq\Delta^{*}.$ (5)

######  Assumption 2.

We also assume that the unknown parameters are bounded, i.e., there exists a
$\theta^{*}_{\max}\geq 0$ such that

$\|\theta_{k}^{*}\|\leq\theta_{\max}^{*}.$ (6)

##  4 Adaptive Law with a Time-Varying Learning Rate

We propose the following new recursive algorithm for estimating
$\theta_{k}^{*}$ as $\theta_{k}$, for all $k\geq 0$:

$\displaystyle\theta_{k+1}$
$\displaystyle=\theta_{k}-\frac{\lambda_{\Gamma}\kappa\Gamma_{k}\phi_{k}e_{k}}{1+\|\phi_{k}\|^{2}},$
(7) $\displaystyle\Omega_{k}$
$\displaystyle=(1-\lambda_{\Omega})\Omega_{k-1}+\frac{\phi_{k}\phi_{k}^{\top}}{1+\|\phi_{k}\|^{2}},$
(8) $\displaystyle\Gamma_{k}$
$\displaystyle=\Gamma_{k-1}+\lambda_{\Gamma}(\Gamma_{k-1}-\kappa\Gamma_{k-1}\Omega_{k}\Gamma_{k-1}),$
(9)

where $\lambda_{\Gamma}$, $\kappa$ are positive scalars and
$\Gamma_{k}\in\mathbb{R}^{N\times N}$ is the time-varying gain matrix. The
initial value $\Gamma_{0}$ is a symmetric positive definite matrix chosen so
that $\Gamma_{0}\preceq\Gamma_{\max}I$, where $\Gamma_{\max}\in\mathbb{R}^{+}$
is the maximum learning rate that is assumed to be given. The information
matrix $\Omega_{k}\in\mathbb{R}^{N\times N}$ captures any excitation that may
be present in the regressor $\phi_{k}$. The initial value of $\Omega_{k}$ is a
symmetric positive definite matrix with $0\prec\Omega_{0}\preceq I$. The
parameter $0<\lambda_{\Omega}<1$ is used for the adjustment of the evolution
of $\Omega_{k}$ and is directly related to the upper bound of $\Omega_{k}$.
The scalars $\lambda_{\Gamma}$, $\kappa$ and $\lambda_{\Omega}$ represent
various weights of the proposed algorithm.

The main innovation in the algorithm above is in the evolution of
$\Gamma_{k}$. Unlike the recursive least squares algorithm (see [4] for
example) with covariance resetting introduced via an externally chosen period,
the algorithm in (9) automatically ensures that $\Gamma_{k}$ stays within
certain bounds, through a combination of linear and nonlinear components. As
will be shown in the following section, stability and exponential decrease of
parameter error towards a compact set are established under appropriate
excitation conditions. In what follows, we assume that the regressors have
some level of excitation, satisfying either Definition 1 or Definition 2. When
these conditions are not satisfied, addition of a suitable projection can
ensure stable estimation [13], details of which are not addressed in this
paper.

##  5 Stability and Convergence Analysis

Before we address the boundedness of the parameter estimates and their
convergence,, we establish the boundedness of the information matrix
$\Omega_{k}$ in (8) and time-varying gain matrix $\Gamma_{k}$ defined in (9).

###  5.1 Boundedness of $\Omega_{k}$

We first show that for any arbitrary regressor $\phi_{k}$, $\Omega_{k}$ is a
positive semi-definite matrix with an upper bound.

######  Lemma 2.

For the dynamic system in (8), if $0<\lambda_{\Omega}<1$ and $\Omega_{0}$ is
such that $0\preceq\Omega_{0}\preceq I$, then
$0\preceq\Omega_{k}\preceq\Omega_{\max}I$, for all $k\geq 1$, where
$\Omega_{\max}=1/\lambda_{\Omega}$.

###### Proof.

We consider the lower bound and the upper bound of $\Omega_{k}$ separately.

  1. 1.

First we prove $\Omega_{k}\succeq 0$ by induction. Let $v\in\mathbb{R}^{N}$
and $v\neq 0$. Since $0<\lambda_{\Omega}<1$ and $\Omega_{0}\succeq 0$,

$\displaystyle
v^{\top}\Omega_{1}v=(1-\lambda_{\Omega})v^{\top}\Omega_{0}v+\frac{\|\phi_{1}^{\top}v\|^{2}}{1+\|\phi_{1}\|^{2}}\geq
0$

Therefore, $\Omega_{1}\succeq 0$. Similarly if
$0\preceq\Omega_{k-1}\preceq\frac{1}{\lambda_{\Omega}}I$, it follows that

$v^{\top}\Omega_{k}v=(1-\lambda_{\Omega})v^{\top}\Omega_{k-1}v+\frac{\|\phi_{k}^{\top}v\|^{2}}{1+\|\phi_{k}\|^{2}}\geq
0.$

  2. 2.

Next we prove $\Omega_{k}\preceq\Omega_{\max}I$. Using convolution,
$\Omega_{k}$ in (8) can be written as

$\displaystyle\Omega_{k}$
$\displaystyle=(1-\lambda_{\Omega})^{k}\Omega_{0}+\sum_{i=0}^{k-1}(1-\lambda_{\Omega})^{k-1-i}\frac{\phi_{i}\phi_{i}^{\top}}{1+\|\phi_{i}\|^{2}}$
$\displaystyle\preceq(1-\lambda_{\Omega})^{k}\Omega_{0}+I\sum_{i=0}^{k-1}(1-\lambda_{\Omega})^{i}.$

Since $\Omega_{0}\leq I$, using simple algebra, we can then show that
$\Omega_{k}\leq\Omega_{\max}I$, for all $k\geq 0$.

∎

Lemma 2 guarantees that $\Omega_{k}$ in (8) is a BIBO-stable system and
determines an explicit upper bound of $\Omega_{k}$. In the next lemma, we
focus on its lower bound, and show that $\Omega_{k}$ is positive definite,
under appropriate excitation conditions.

######  Lemma 3.

If $0<\lambda_{\Omega}<1$, and $\phi_{k}$ is persistently exciting over a
period $\Delta T$ at a level $\alpha$ for all $k\geq k_{1}$, then
$\Omega_{k}\succeq\Omega_{PE}I$, $\forall k\geq k_{2}$, where
$\Omega_{PE}=(1-\lambda_{\Omega})^{\Delta T-1}\alpha$ and $k_{2}=k_{1}+\Delta
T$.

###### Proof.

Using convolution, $\Omega_{k}$ can be rewritten as

$\displaystyle\Omega_{k}$
$\displaystyle=(1-\lambda_{\Omega})^{k}\Omega_{0}+\sum_{i=0}^{k-1}(1-\lambda_{\Omega})^{k-i-1}\frac{\phi_{i}\phi_{i}^{\top}}{1+\|\phi_{i}\|^{2}}$

Since $\phi_{k}$ is persistently exciting,

$\displaystyle\Omega_{k}$
$\displaystyle\succeq(1-\lambda_{\Omega})^{k-1}\sum_{i=k-\Delta
T}^{k-1}(1-\lambda_{\Omega})^{-i}\frac{\phi_{i}\phi_{i}^{\top}}{1+\|\phi_{i}\|^{2}}$
$\displaystyle\succeq\Omega_{PE}I$

∎

A weaker condition on $\phi_{k}$ is finite excitation. In Lemma 4, we show the
lower bound of $\Omega_{k}$ when $\phi_{k}$ is finitely exciting.

######  Lemma 4.

If $0<\lambda_{\Omega}<1$ in (8), $\phi_{k}$ is finitely exciting at a level
$\alpha$ over $[k_{1}^{\prime},k_{2}^{\prime}]$, and
$0\preceq\Omega_{0}\preceq I$, then $\Omega_{k}\succeq\Omega_{FE}I,\;\forall
k\in[k_{2}^{\prime},k_{3}^{\prime}]$, where $k_{3}^{\prime}$ is a finite
instant with $k_{3}^{\prime}>k_{2}^{\prime}$, and
$\Omega_{FE}=\alpha\lambda_{\Omega}(1-\lambda_{\Omega})^{k_{3}^{\prime}-k_{1}^{\prime}-1}$.

###### Proof.

Using convolution, we can write $\Omega_{k}$ as

$\displaystyle\Omega_{k}$
$\displaystyle=(1-\lambda_{\Omega})^{k}\Omega_{0}+\sum_{i=0}^{k-1}(1-\lambda_{\Omega})^{k-i-1}\frac{\phi_{i}\phi_{i}^{\top}}{1+\|\phi_{i}\|^{2}}$
$\displaystyle\succeq(1-\lambda_{\Omega})^{k-1}\sum_{i=k_{1}^{\prime}}^{k_{2}^{\prime}}(1-\lambda_{\Omega})^{-i}\frac{\phi_{i}\phi_{i}^{\top}}{1+\|\phi_{i}\|^{2}}$
$\displaystyle\succeq\alpha(1-\lambda_{\Omega})^{k-k_{1}^{\prime}-1}I$

Since $(1-\lambda_{\Omega})^{k-k_{1}^{\prime}-1}$ is a decreasing function of
$k$ and $\Omega_{k}\succeq\Omega_{FE}I$ for $k_{2}^{\prime}\leq k\leq
k_{3}^{\prime}$, the lemma holds. ∎

######  Remark 1.

It should be noted that when $\phi_{k}$ is finitely exciting, $\Omega_{k}$ is
positive definite but not uniformly positive definite. This follows since
$\Omega_{FE}$ is a decreasing function of $k_{3}^{\prime}$. The effects of
persistent excitation and finite excitation on the regressor, $\Omega_{k}$ and
$\Gamma_{k}$ are shown in Table 1 and 2, respectively.

$k\in$ | $[k_{1},k_{2}]$ | $k_{2}$ | $[k_{2},\infty)$
---|---|---|---
$\|\phi_{k}\phi_{k}^{\top}\|$ | $\geq 0$ |  |
$\|\sum_{k=k_{1}}^{k_{2}}\phi_{k}\phi_{k}^{\top}\|$ |  | $\geq\alpha$ |
$\|\Omega_{k}\|$ |  |  | $\geq\Omega_{PE}$
$\|\Gamma_{k}\|$ |  |  | $\geq\Gamma_{PE}$
Table 1: Norm of Signals in Different Phases of Persistent Excitation Propagation $k\in$ | $[k_{1}^{\prime},k_{2}^{\prime}]$ | $k_{2}^{\prime}$ | $[k_{2}^{\prime},k_{3}^{\prime}]$
---|---|---|---
$\|\phi_{k}\phi_{k}^{\top}\|$ | $\geq 0$ |  |
$\|\sum_{k=k_{1}^{\prime}}^{k_{2}^{\prime}}\phi_{k}\phi_{k}^{\top}\|$ |  | $\geq\alpha$ |
$\|\Omega_{k}\|$ |  |  | $\geq\Omega_{FE}$
$\|\Gamma_{k}\|$ |  |  | $\geq\Gamma_{FE}$
Table 2: Norm of Signals in Different Phases of Finite Excitation Propagation

###  5.2 Boundedness of $\Gamma_{k}$

One of the major contributions of this paper is the adjustment of the time-
varying learning rate matrix $\Gamma_{k}$ in (9). Unlike covariance resetting,
(9) allows the linear part and the quadratic part to automatically balance
each other so that $\Gamma_{k}$ is always bounded, under conditions of
persistent excitation. In this section we will show this result formally.

We first consider the scalar case of (9), given by

$\gamma_{k}=\gamma_{k-1}+\lambda_{\Gamma}(\gamma_{k-1}-\kappa\omega_{k}\gamma_{k-1}^{2}),$
(10)

where

$0<\omega_{\min}\leq\omega_{k}\leq\omega_{\max}.$ (11)

Let $\gamma_{\max}>0$ and

$\gamma_{\min}=\min\left\\{\frac{1}{\kappa\omega_{\max}},\gamma_{\max}+\lambda_{\Gamma}(\gamma_{\max}-\kappa\omega_{\max}\gamma_{\max}^{2})\right\\}$
(12)

where $\kappa_{\min}\leq\kappa<\kappa_{\max}$, and $\kappa_{\min}$,
$\kappa_{\max}$ are defined as

$\kappa_{\min}=\frac{1}{\gamma_{\max}\omega_{\min}},\;\;\kappa_{\max}=\frac{1+\lambda_{\Gamma}}{\lambda_{\Gamma}\omega_{\max}\gamma_{\max}}.$
(13)

In the update of (10), $\gamma_{\max}$ and $\lambda_{\Gamma}$ are chosen
arbitrarily by the user. The gain $\kappa$ is then chosen using the bounds in
(13). We now state the following lemma that establishes bounds of $\gamma_{k}$
(See Fig. 1 and Fig. 2 for an illustration of the bounds) which in turn depend
on $\lambda_{\Gamma},\gamma_{\max}$ and $\kappa$.

Figure 1: A parabola for illustrating the bounds of $\gamma_{k}$ whose update
is in (10) for the case when
$\gamma_{\max}>\frac{1+\lambda_{\Gamma}}{2\lambda_{\Gamma}\kappa\omega_{k}}$.
For a given $\omega_{k}$, the line $\gamma_{k}=\gamma_{k-1}$ intersects the
parabola at $\gamma_{E}$, giving
$\gamma_{k}=\gamma_{k-1}=\frac{1}{\kappa\omega_{k}}$. When the parabola is
above the straight line, $\gamma_{k}\geq\gamma_{k-1}$; Otherwise
$\gamma_{k}\leq\gamma_{k-1}$. The point $\gamma_{\max}(\omega_{k})$ is the
maximum obtainable value for $\gamma_{k}$ and it will be no larger than
$\gamma_{\max}$. When $\gamma_{k-1}=\gamma_{\max}$, the resulting $\gamma_{k}$
will be no less than $\gamma_{\min}$. Figure 2: Illustration of the upper
bounds of $\gamma_{k}$ as $\omega_{k}$ varies. The blue curve shows how
$\gamma_{\max}(\omega_{k})$ changes. The value of $\omega^{*}$ is defined such
that (14) is satisfied. Under the assumptions for Lemma 5,
$\gamma_{\max}(\omega_{k})$ gets larger as $\omega_{k}$ decreases. The global
maximum can be obtained when $\omega_{k}=\omega_{\min}$ and
$\gamma_{k-1}=\gamma_{\max}$.

######  Lemma 5.

If
$0<\lambda_{\Gamma}<\min\left\\{\frac{1}{\omega_{\max}/\omega_{\min}-1},1\right\\}$,
and $\gamma_{\min}\leq\gamma_{0}\leq\gamma_{\max}$, then the solutions of (10)
are such that $0<\gamma_{\min}\leq\gamma_{k}\leq\gamma_{\max}$, $\forall k\geq
0$.

###### Proof.

Defining the right hand side of Eq. (10) as $f(\gamma_{k-1},\omega_{k})$, it
is easy to see that $f(.,.)$ is quadratic in $\gamma_{k-1}$, with one maximum
and two minima (see Fig. 1 and Fig. 2). Suppose $\omega^{*}$ is such that

$\frac{1+\lambda_{\Gamma}}{2\lambda_{\Gamma}\kappa\omega^{*}}=\gamma_{\max}.$
(14)

It is easy to check that $\omega^{*}\geq\omega_{\min}$. If
$\omega_{k}\geq\omega^{*}$,

$\displaystyle f(\gamma_{k-1},\omega_{k})$
$\displaystyle\leq\gamma_{\max}+\lambda_{\Gamma}(\gamma_{\max}-\kappa\omega^{*}\gamma_{\max}^{2})$
$\displaystyle\leq\gamma_{\max}+\lambda_{\Gamma}(\gamma_{\max}-\kappa\omega_{\min}\gamma_{\max}^{2}).$

If $\omega_{k}\leq\omega^{*}$,

$\displaystyle f(\gamma_{k-1},\omega_{k})$
$\displaystyle\leq\gamma_{\max}+\lambda_{\Gamma}(\gamma_{\max}-\kappa\omega_{k}\gamma_{\max}^{2})$
$\displaystyle\leq\gamma_{\max}+\lambda_{\Gamma}(\gamma_{\max}-\kappa\omega_{\min}\gamma_{\max}^{2}).$

Noticing the bounds of $\lambda_{\Gamma}$ and $\kappa$, simple algebraic
manipulations show that

$\displaystyle\max_{\gamma_{k},\omega_{k}}f(\gamma_{k-1},\omega_{k})$
$\displaystyle=\gamma_{\max}+\lambda_{\Gamma}(\gamma_{\max}-\kappa\omega_{\min}\gamma_{\max}^{2})$
$\displaystyle\leq\gamma_{\max}.$

Similarly, it can be shown by inspection that

$\min_{\omega_{k}}f(\gamma_{k-1},\omega_{k})=\gamma_{k-1}+\lambda_{\Gamma}(\gamma_{k-1}-\kappa\omega_{\max}\gamma_{k-1}^{2}).$

Since the largest decrease of $f$ occurs when $\gamma_{k}$ is at its maximum
value, it follows that

$\min_{\gamma_{k},\omega_{k}}f(\gamma_{k-1},\omega_{k})=\gamma_{\max}+\lambda_{\Gamma}(\gamma_{\max}-\kappa\omega_{\max}\gamma_{\max}^{2}).$

In addition, $\kappa\omega_{\max}>0$. Therefore, the lower bound of
$\gamma_{k}$ and that it is positive are ensured through the choice of
$\gamma_{\min}$ as in (12).∎

######  Remark 2.

It should also be noted that (11) implies that we have implicitly assumed that
$\omega_{k}$ is persistently exciting, i.e., we consider the period
$k\in[k_{2},\infty)$ as shown in Table 1. If instead $\omega_{k}$ is finitely
exciting, the bounds on $\gamma_{k}$ are only assured for a finite period,
i.e., $k\in[k_{2}^{\prime},k_{3}^{\prime}]$ as shown in Table 2, rather than
for all $k\geq 0$.

######  Remark 3.

As we will show next, the properties we obtain in Lemma 5 for the scalar
$\gamma_{k}$ will be directly utilized in showing the properties of the matrix
$\Gamma_{k}$. Specifically, it will be shown that each eigenvalue of
$\Gamma_{k}$ has an update similar to (10).

We now proceed to the matrix case. First we consider the case when $\phi_{k}$
is persistently exciting and thus
$\Omega_{PE}I\preceq\Omega_{k}\preceq\Omega_{\max}I$, $\forall k\geq k_{2}$,
as shown in Lemma 3. Let $\Gamma_{\max}>0$,
$\kappa_{\min}=\frac{1}{\Gamma_{\max}\Omega_{PE}}$,
$\kappa_{\max}=\frac{1+\lambda_{\Gamma}}{\lambda_{\Gamma}\Omega_{\max}\Gamma_{\max}}$,
$\Gamma_{PE}=\min\left\\{\frac{1}{\kappa\Omega_{\max}},\Gamma_{\max}+\lambda_{\Gamma}(\Gamma_{\max}-\kappa\Omega_{\max}\Gamma_{\max}^{2})\right\\}$.

######  Lemma 6.

If $\phi_{k}$ is persistently exciting,
$\Omega_{PE}I\preceq\Omega_{k}\preceq\Omega_{\max}I$, $\forall k\geq k_{2}$,
learning rate gain
$0<\lambda_{\Gamma}<\min\left\\{\frac{1}{\Omega_{\max}/\Omega_{PE}-1},1\right\\}$,
gain $\kappa_{\min}\leq\kappa<\kappa_{\max}$, and initial value
$\Gamma_{PE}I\preceq\Gamma_{k_{2}}\preceq\Gamma_{\max}I$, then
$\Gamma_{PE}I\preceq\Gamma_{k}\preceq\Gamma_{\max}I$, $\forall k>k_{2}$ for
some $k_{2}>k_{1}$.

###### Proof.

We first establish the lower bound $\Gamma_{PE}$ using induction: Let
$\Gamma_{PE}I\preceq\Gamma_{k-1}\preceq\Gamma_{\max}I$. It is easy to check
that $\Gamma_{PE}>0$. Since $\Gamma_{k-1}$ is positive definite, we can
factorize $\Gamma_{k-1}$ as[23]

$\Gamma_{k-1}=Q_{k-1}\Lambda_{k-1}Q_{k-1}^{\top}$ (15)

where $\Lambda_{k-1}$ is diagonal and all of its entries are the eigenvalues
of $\Gamma_{k-1}$, and $Q_{k-1}$ is orthogonal. Using Eq. (9) and the fact
that $Q_{k-1}$ is orthogonal, we obtain that

$\displaystyle\Gamma_{k}$
$\displaystyle=\Gamma_{k-1}+\lambda_{\Gamma}(\Gamma_{k-1}-\kappa\Gamma_{k-1}\Omega_{k}\Gamma_{k-1})$
$\displaystyle\succeq\Gamma_{k-1}+\lambda_{\Gamma}(\Gamma_{k-1}-\kappa\Omega_{\max}\Gamma_{k-1}^{2})$
$\displaystyle=Q_{k-1}\left[\Lambda_{k-1}+\lambda_{\Gamma}(\Lambda_{k-1}-\kappa\Omega_{\max}\Lambda_{k-1}^{2})\right]Q_{k-1}^{\top},$

since

$x^{\top}(\Omega_{\max}I-\Omega_{k})x\geq 0,\;\forall x\neq 0.$

Since $\Gamma_{k-1}$ is positive definite, it follows that

$v^{\top}\Gamma_{k-1}(\Omega_{\max}I-\Omega_{k})\Gamma_{k-1}v\geq 0,\;\forall
v\neq 0.$

Thus

$\Gamma_{k-1}(\Omega_{\max}I-\Omega_{k})\Gamma_{k-1}\succeq 0.$ (16)

Denote the $i$th diagonal element of $\Gamma_{k-1}$ as $\lambda_{k-1,i}$,
$i=1,\cdots,N$. Since $\Gamma_{PE}I\preceq\Gamma_{k-1}\preceq\Gamma_{\max}I$,
we have $\Gamma_{PE}\leq\lambda_{k-1,i}\leq\Gamma_{\max}$, $\forall i$.
Applying Lemma 5, we obtain $\Gamma_{PE}\leq\lambda_{k,i}\leq\Gamma_{\max}$,
$\forall i$ and the lower bound of $\Gamma_{k}$ is proved. Since the lower
bound holds for $k=k_{2}$, proof by induction guarantees that the lower bound
holds for all $k\geq k_{2}$.

To prove the upper bound, similar to Eq. (16), using Lemma 3 we can show that

$\Gamma_{k-1}(\Omega_{PE}I-\Omega_{k})\Gamma_{k-1}\preceq 0.$ (17)

Therefore

$\displaystyle\Gamma_{k}\preceq
Q_{k-1}\left[\Lambda_{k-1}+\lambda_{\Gamma}(\Lambda_{k-1}-\kappa\Omega_{PE}\Lambda_{k-1}^{2})\right]Q_{k-1}^{\top}.$

Using the same procedure as above, it can be shown that the upper bound holds
for $\Gamma_{k}$ if it holds for $\Gamma_{k-1}$, and using induction that the
upper bound holds for all $k$, Lemma 6 is proved. ∎

When $\phi_{k}$ is finitely exciting in a period
$[k_{1}^{\prime},k_{2}^{\prime}]$, we have shown in Lemma 2 and Lemma 4 that
$\Omega_{FE}I\preceq\Omega_{k}\preceq\Omega_{\max}I$, $\forall
k_{2}^{\prime}\leq k\leq k_{3}^{\prime}$. For the following Lemma, we let
$\Gamma_{\max}>0$, and denote
$\kappa_{\min}=\frac{1}{\Gamma_{\max}\Omega_{FE}}$,
$\kappa_{\max}=\frac{1+\lambda_{\Gamma}}{\lambda_{\Gamma}\Omega_{\max}\Gamma_{\max}}$
and
$\Gamma_{FE}=\min\left\\{\frac{1}{\kappa\Omega_{\max}},\Gamma_{\max}+\lambda_{\Gamma}(\Gamma_{\max}-\kappa\Omega_{\max}\Gamma_{\max}^{2})\right\\}$.

######  Lemma 7.

If $\phi_{k}$ is finitely exciting over $[k_{1}^{\prime},k_{2}^{\prime}]$ with
level $\alpha$, and resulting in
$\Omega_{FE}I\preceq\Omega_{k}\preceq\Omega_{\max}I$ over
$[k_{2}^{\prime},k_{3}^{\prime}]$, learning rate gain
$0<\lambda_{\Gamma}<\min\left\\{\frac{1}{\Omega_{\max}/\Omega_{FE}-1},1\right\\}$,
gain $\kappa_{\min}\leq\kappa<\kappa_{\max}$, and initial value
$\Gamma_{FE}I\preceq\Gamma_{k_{2}^{\prime}}\preceq\Gamma_{\max}I$, then
$\Gamma_{FE}I\preceq\Gamma_{k}\preceq\Gamma_{\max}I$, for $k_{2}^{\prime}\leq
k\leq k_{3}^{\prime}$, for some $k_{3}^{\prime}\geq k_{2}^{\prime}+1$.

###### Proof.

The proof is very similar to the proof for Lemma 6. But the guaranteed
boundedness period changes to $[k_{2}^{\prime},k_{3}^{\prime}]$. ∎

######  Remark 4.

Note that Lemma 6 and Lemma 7 are essentially extensions to Lemma 5. The
extension to the matrix $\Gamma_{k}$ from scalar $\gamma_{k}$ has been carried
out by diagonalization and applying Lemma 5. In Lemma 6 and Lemma 7, it should
be noted that $\Omega_{PE}$ and $\Omega_{FE}$ are corresponding to
$\omega_{\min}$ in Lemma 5, respectively.

###  5.3 Stability and Parameter Convergence

We now return to the main algorithm in (7), (8) and (9). Note that in Section
5.1 and Section 5.2 we have shown that both $\Omega_{k}$ and $\Gamma_{k}$ are
bounded under different excitation conditions. In this section we consider the
boundedness of $\theta_{k}$ as well as its convergence.

We first introduce an intermediate variable in the evolution of $\Gamma_{k}$
in (9) in the form of a matrix $\bar{\Gamma}_{k}$,

$\bar{\Gamma}_{k}=\Gamma_{k-1}-\lambda_{\Gamma}\kappa\Gamma_{k-1}\frac{\phi_{k-1}\phi_{k-1}^{\top}}{1+\|\phi_{k-1}\|^{2}}\Gamma_{k-1},$
(18)

and consider a Lyapunov function candidate

$V_{k}=\tilde{\theta}_{k}\bar{\Gamma}_{k}^{-1}\tilde{\theta}_{k}.$ (19)

Under conditions of excitation, $V_{k}$ will be shown to be a Lyapunov
function and converge to a compact set that scales with $\Delta^{*}$.

We first consider the case when $\phi_{k}$ is persistently exciting. In this
case, as proved in Lemma 3, $\Omega_{PE}I\preceq\Omega_{k}$, $\forall k\geq
k_{2}$ for some $k_{2}>0$. Let

$\lambda_{\Omega}>\frac{1}{\Omega_{PE}+1},$ (20)
$0<\lambda_{\Gamma}<\min\left\\{\frac{1}{\Omega_{\max}/\Omega_{PE}-1},1,\Omega_{PE}\right\\}.$
(21)

We define

$\displaystyle\kappa_{\min}$
$\displaystyle=\frac{1}{\Gamma_{\max}\Omega_{PE}},$ (22)
$\displaystyle\kappa_{\max}$
$\displaystyle=\min\Bigg{\\{}\frac{1+\lambda_{\Gamma}}{\lambda_{\Gamma}\Omega_{\max}\Gamma_{\max}},\frac{1}{\lambda_{\Gamma}\Gamma_{\max}},\allowbreak\frac{1}{(1-\lambda_{\Omega})\Omega_{\max}\Gamma_{\max}}\Bigg{\\}},$
(23)

and choose any $\kappa$ that satisfies

$\kappa_{\min}\leq\kappa<\kappa_{\max}.$ (24)

According to Lemma 6, $\Gamma_{PE}I\preceq\Gamma_{k}\preceq\Gamma_{\max}I$,
and $\Gamma_{PE}$ is defined as

$\Gamma_{PE}=\min\left\\{\frac{1}{\kappa\Omega_{\max}},\Gamma_{\max}+\lambda_{\Gamma}(\Gamma_{\max}-\kappa\Omega_{\max}\Gamma_{\max}^{2})\right\\}.$
(25)

Furthermore, let

$\bar{\Gamma}_{\min}=\min\left\\{\Gamma_{\max}-\lambda_{\Gamma}\kappa\Gamma_{\max}^{2},\Gamma_{PE}-\lambda_{\Gamma}\kappa\Gamma_{PE}^{2}\right\\}.$
(26)

We give the expression of the upper-bound for the decay rate as

$\mu_{1}=\Psi_{\min}\bar{\Gamma}_{\min},$ (27)

and consider a parameter $\mu_{2}$ that satisfies $0<\mu_{2}<\mu_{1}$. In
(27),
$\Psi_{\min}=(1+\lambda_{\Gamma}\bar{\Gamma}_{\min}^{-1}\Gamma_{\max})^{-1}\Upsilon_{\min}\Gamma_{\max}^{-2}$,
where
$\Upsilon_{\min}=\min\big{\\{}\lambda_{\Gamma}\Gamma_{\max}\left[1-\kappa(1-\lambda_{\Omega})\Omega_{\max}\Gamma_{\max}\right]$,
$\lambda_{\Gamma}\Gamma_{PE}\left[1-\kappa(1-\lambda_{\Omega})\Omega_{\max}\Gamma_{PE}\right]\big{\\}}$.
Finally, a compact set is defined as

$D_{PE}=\left\\{V|V\leq V_{PE}\right\\},$ (28)

where

$V_{PE}=\left[\frac{\Delta^{*}\left(c_{1}+\sqrt{c_{1}^{2}+4c_{2}(\mu_{1}-\mu_{2})}\right)}{2(\mu_{1}-\mu_{2})}\right]^{2},$
(29) $c_{1}=\bar{\Gamma}_{\max}^{1/2}\bar{\Gamma}_{\min}^{-1},$ (30)
$c_{2}=\bar{\Gamma}_{\min}^{-1}.$ (31)

We now state the stability result with persistent excitation.

######  Theorem 1.

Let $\phi_{k}$ be persistently exciting. For the ARMA model in (1), under the
assumptions in 1 and 2, subject to the bounds in (5), if the hyperparameters
$\lambda_{\Omega}$, $\lambda_{\Gamma}$ and $\kappa$ satisfy (20), (21) and
(24), respectively, the adaptive algorithm in (7), (8) and (9) guarantees that
in $D_{PE}^{c}$

$V_{k}\leq\exp\left[-\mu_{2}(k-k_{2})\right]V_{k_{2}},\;\forall k\geq k_{2}$

###### Proof.

Starting from (7), we get

$\theta_{k}-\theta_{k-1}^{*}=\left(I-\frac{\lambda_{\Gamma}\kappa\Gamma_{k-1}\phi_{k-1}\phi_{k-1}^{\top}}{1+\|\phi_{k-1}\|^{2}}\right)\tilde{\theta}_{k-1}.$
(32)

Using (18), we rewrite $\Gamma_{k}$ as

$\Gamma_{k}=\bar{\Gamma}_{k}+\underbrace{\lambda_{\Gamma}\Gamma_{k-1}-\lambda_{\Gamma}\kappa(1-\lambda_{\Omega})\Gamma_{k-1}\Omega_{k-1}\Gamma_{k-1}}_{\Upsilon_{k}}$
(33)

Now we show that under the conditions in Theorem 1,
$\bar{\Gamma}_{k}\succeq\bar{\Gamma}_{\min}I\succ 0$. With simple
manipulations, it can be shown that

$\bar{\Gamma}_{k}\succeq
Q_{k-1}\left(\Lambda_{k-1}-\lambda_{\Gamma}\kappa\Lambda_{k-1}^{2}\right)Q_{k-1}^{\top}.$
(34)

It is easy to see that $\Gamma_{PE}\leq\lambda_{k-1,i}\leq\Gamma_{\max}$,
$\forall i$, where $\lambda_{k-1,i}$ is the $i$th diagonal element of
$\Gamma_{k-1}$. From the inequality in (34), $\bar{\Gamma}_{k}$ assumes its
lower bound either when $\lambda_{k-1,i}=\Gamma_{PE}$ or when
$\lambda_{k-1,i}=\Gamma_{\max}$. It follows therefore that
$\bar{\Gamma}_{k}\succeq\bar{\Gamma}_{\min}I$, where $\bar{\Gamma}_{\min}>0$
is defined in (26). It is also obvious that
$\bar{\Gamma}_{k}\preceq\bar{\Gamma}_{\max}I$, where
$\bar{\Gamma}_{\max}=\Gamma_{\max}$. Using the same approach, we can also show
that $0\prec\Upsilon_{\min}I\preceq\Upsilon_{k}\preceq\Upsilon_{\max}I$, where
$\Upsilon_{k}$ is defined in (33) and
$\Upsilon_{\max}=\lambda_{\Gamma}\Gamma_{\max}$. With these bounds
established, we consider the Lyapunov function in (19). Since
$\tilde{\theta}_{k}=\theta_{k-1}^{*}-\theta_{k}^{*}+\bar{\Gamma}_{k}\Gamma_{k-1}^{-1}\tilde{\theta}_{k-1}$,
we obtain that

$\displaystyle V_{k}$
$\displaystyle=\left(\theta_{k-1}^{*}-\theta_{k}^{*}+\bar{\Gamma}_{k}\Gamma_{k-1}^{-1}\tilde{\theta}_{k-1}\right)^{\top}\bar{\Gamma}_{k}^{-1}$
$\displaystyle\hskip
14.22636pt\left[\theta_{k-1}^{*}-\theta_{k}^{*}+\left(I-\lambda_{\Gamma}\kappa\frac{\Gamma_{k-1}\phi_{k-1}\phi_{k-1}^{\top}}{1+\|\phi_{k-1}\|^{2}}\right)\tilde{\theta}_{k-1}\right]$
$\displaystyle\leq\tilde{\theta}_{k-1}^{\top}\Gamma_{k-1}^{-1}\tilde{\theta}_{k-1}+2\tilde{\theta}_{k-1}^{\top}\Gamma_{k-1}^{-1}(\theta_{k-1}^{*}-\theta_{k}^{*})$
$\displaystyle\hskip 8.5359pt+(\Delta^{*})^{2}\bar{\Gamma}_{\min}^{-1}.$

Since $\Upsilon_{k}\succ 0$, we can write
$\Upsilon_{k}=\Upsilon_{k}^{1/2}\Upsilon_{k}^{1/2}$. Applying Lemma 1 to (33),
we have

$\Gamma_{k}^{-1}=\bar{\Gamma}_{k}^{-1}-\underbrace{\bar{\Gamma}_{k}^{-1}\Upsilon_{k}^{1/2}\left(I+\Upsilon_{k}^{1/2}\bar{\Gamma}_{k}^{-1}\Upsilon_{k}^{1/2}\right)^{-1}\Upsilon_{k}^{1/2}\bar{\Gamma}_{k}^{-1}}_{\Psi_{k}}.$
(35)

Since $\Upsilon_{k}^{1/2}\succ 0$ and $\bar{\Gamma}_{k}^{-1}\succ 0$, applying
the definition of positive definite matrices we obtain that
$\Upsilon_{k}^{1/2}\bar{\Gamma}_{k}^{-1}\Upsilon_{k}^{1/2}\succ 0$. Also, it
follows that
$\Upsilon_{k}^{1/2}\bar{\Gamma}_{k}^{-1}\Upsilon_{k}^{1/2}\preceq\lambda_{\Gamma}\bar{\Gamma}_{\min}^{-1}\Gamma_{\max}I$.
Therefore
$\left(I+\Upsilon_{k}^{1/2}\bar{\Gamma}_{k}^{-1}\Upsilon_{k}^{1/2}\right)^{-1}\succeq(1+\lambda_{\Gamma}\bar{\Gamma}_{\min}^{-1}\Gamma_{\max})^{-1}I$.
Applying the definition of positive definite matrices again, we get
$\Psi_{k}\succeq\Psi_{\min}I\succ 0$. Substitute $\Gamma_{k}^{-1}$ in (35)
into $V_{k}$, we get

$\displaystyle\begin{split}V_{k}&\leq\tilde{\theta}_{k-1}^{\top}\Gamma_{k-1}^{-1}\tilde{\theta}_{k-1}+\Delta^{*}\bar{\Gamma}_{\max}^{1/2}\bar{\Gamma}_{\min}^{-1}\sqrt{V_{k-1}}+(\Delta^{*})^{2}\bar{\Gamma}_{\min}^{-1}\\\
&\leq(1-\underbrace{\Psi_{\min}\bar{\Gamma}_{\min}}_{\mu_{1}})V_{k-1}+\Delta^{*}\underbrace{\bar{\Gamma}_{\max}^{1/2}\bar{\Gamma}_{\min}^{-1}}_{c_{1}}\sqrt{V_{k-1}}\\\
&\hskip
8.5359pt+(\Delta^{*})^{2}\underbrace{\bar{\Gamma}_{\min}^{-1}}_{c_{2}},\end{split}$
(36)

where we have replaced $\tilde{\theta}_{k-1}$ with an upper bound involving
$V_{k-1}$. From (35) and $\Gamma_{k}^{-1}\succ 0$, it follows that
$\Psi_{\min}\bar{\Gamma}_{\min}<1$. Therefore, we have

$\displaystyle\Delta V_{k}$ $\displaystyle=V_{k}-V_{k-1}$
$\displaystyle\leq-\mu_{1}V_{k-1}+c_{1}\Delta^{*}\sqrt{V_{k-1}}+c_{2}(\Delta^{*})^{2},$

where $\mu_{1}$ is defined in (27). Choosing any $0<\mu_{2}<\mu_{1}$, it can
be noted that $\Delta V_{k}<0$ in $D_{PE}^{c}$, where the compact set $D_{PE}$
is defined in (28). For any $k$ such that $V_{k}>V_{PE}$, it follows that

$V_{k}\leq(1-\mu_{2})V_{k-1}$

We therefore obtain that when $k\geq k_{2}$, for all $V_{k}>V_{PE}$

$\displaystyle V_{k}$ $\displaystyle\leq(1-\mu_{2})^{k-k_{2}}V_{k_{2}}$
$\displaystyle\leq\exp\left[-\mu_{2}(k-k_{2})\right]V_{k_{2}}.$

∎

######  Remark 5.

The key steps that enabled the proof of Theorem 1 was 1) to express
$\Gamma_{k}=\bar{\Gamma}_{k}+\Upsilon_{k}$, and 2) to define $V_{k}$ only
using $\bar{\Gamma}_{k}$ as in (18) rather than $\Gamma_{k}$. This provided
the necessary tractability as it avoided matrix inversions involving
$\Gamma_{k}$. From (18), it can be seen that $\bar{\Gamma}_{k}$ is similar to
the standard gain matrix used in recursive least squares [4] and therefore
retains the desired positive-definiteness property with suitable choice of the
hyperparameters in (25). This in turn allowed the proof of Theorem 1.

We now address the case when $\phi_{k}$ is finitely exciting. As proved in
Lemma 4, it follows that for $k\in[k_{2}^{\prime},k_{3}^{\prime}]$, we have
$\Omega_{FE}I\preceq\Omega_{k}$ and
$\Gamma_{FE}I\preceq\Gamma_{k}\preceq\Gamma_{\max}I$. Let

$\lambda_{\Omega}>\frac{1}{\Omega_{FE}+1},$ (37)
$0<\lambda_{\Gamma}<\min\left\\{\frac{1}{\Omega_{\max}/\Omega_{FE}-1},1,\Omega_{FE}\right\\}.$
(38)

We define

$\displaystyle\kappa_{\min}$
$\displaystyle=\frac{1}{\Gamma_{\max}\Omega_{FE}},$ (39)
$\displaystyle\kappa_{\max}$
$\displaystyle=\min\left\\{\frac{1+\lambda_{\Gamma}}{\lambda_{\Gamma}\Omega_{\max}\Gamma_{\max}},\frac{1}{\lambda_{\Gamma}\Gamma_{\max}},\frac{1}{(1-\lambda_{\Omega})\Omega_{\max}\Gamma_{\max}}\right\\},$
(40)

and choose any $\kappa$ that satisfies

$\kappa_{\min}\leq\kappa<\kappa_{\max}.$ (41)

According to Lemma 7, $\Gamma_{FE}I\preceq\Gamma_{k}\preceq\Gamma_{\max}I$,
and $\Gamma_{FE}$ is defined as

$\Gamma_{FE}=\min\left\\{\frac{1}{\kappa\Omega_{\max}},\Gamma_{\max}+\lambda_{\Gamma}(\Gamma_{\max}-\kappa\Omega_{\max}\Gamma_{\max}^{2})\right\\}.$
(42)

Furthermore, let

$\bar{\Gamma}_{\min}=\min\left\\{\Gamma_{\max}-\lambda_{\Gamma}\kappa\Gamma_{\max}^{2},\Gamma_{FE}-\lambda_{\Gamma}\kappa\Gamma_{FE}^{2}\right\\},$
(43)

The upper-bound for the decay rate will be

$\mu_{1}^{\prime}=\Psi_{\min}\bar{\Gamma}_{\min},$ (44)

and $\mu_{2}^{\prime}$ is any value that satisfies
$0<\mu_{2}^{\prime}<\mu_{1}^{\prime}$. In (44),
$\Psi_{\min}=(1+\lambda_{\Gamma}\bar{\Gamma}_{\min}^{-1}\Gamma_{\max})^{-1}\Upsilon_{\min}\Gamma_{\max}^{-2}$,
where
$\Upsilon_{\min}=\min\big{\\{}\lambda_{\Gamma}\Gamma_{\max}\left[1-\kappa(1-\lambda_{\Omega})\Omega_{\max}\Gamma_{\max}\right]$,
$\lambda_{\Gamma}\Gamma_{FE}\left[1-\kappa(1-\lambda_{\Omega})\Omega_{\max}\Gamma_{FE}\right]\big{\\}}$.
For the purpose of convergence statement, we consider the same update for
$\bar{\Gamma}_{k}$ in (18) and Lyapunov function in (19). Finally, a compact
set is defined as

$D_{FE}=\left\\{V|V\leq V_{FE}\right\\},$ (45)

where

$V_{FE}=\left[\frac{\Delta^{*}\left(c_{1}^{\prime}+\sqrt{(c_{1}^{\prime})^{2}+4c_{2}^{\prime}(\mu_{1}^{\prime}-\mu_{2}^{\prime})}\right)}{2(\mu_{1}^{\prime}-\mu_{2}^{\prime})}\right]^{2},$
(46) $c_{1}^{\prime}=\bar{\Gamma}_{\max}^{1/2}\bar{\Gamma}_{\min}^{-1},$ (47)
$c_{2}^{\prime}=\bar{\Gamma}_{\min}^{-1}.$ (48)

######  Theorem 2.

Let $\phi_{k}$ be finitely exciting for $k\in[k_{1}^{\prime},k_{2}^{\prime}]$.
For the ARMA model in (1), under the assumptions in 1 and 2, subject to the
bounds in (5), if the hyperparameters $\lambda_{\Omega}$, $\lambda_{\Gamma}$
and $\kappa$ satisfy (37), (38) and (41), respectively, the adaptive algorithm
in (7), (8) and (9) guarantees that in $D_{FE}^{c}$

$V_{k}\leq\exp\left[-\mu_{2}^{\prime}(k-k_{2}^{\prime})\right]V_{k_{2}^{\prime}},\;k\in[k_{2}^{\prime},k_{3}^{\prime}]$
(49)

###### Proof.

The proof is similar to that of Theorem 1, with the distinction in the lower
bound of $\Omega_{k}$ and the effective period. In particular, the decay of
$V_{k}$ becomes

$\displaystyle\Delta V_{k}$ $\displaystyle=V_{k}-V_{k-1}$
$\displaystyle\leq-\mu_{1}^{\prime}V_{k-1}+c_{1}^{\prime}\Delta^{*}\sqrt{V_{k-1}}+c_{2}^{\prime}(\Delta^{*})^{2},$

which holds only for $k_{2}^{\prime}\leq k\leq k_{3}^{\prime}$, where
$k_{3}^{\prime}\geq k_{2}^{\prime}+1$. ∎

######  Remark 6.

Unlike (49) in Theorem 2, in Theorem 1, the convergence properties hold for
all $k\geq k_{2}$. In the case of finite excitation, after $k=k_{3}^{\prime}$,
$\Omega_{k}$ reaches its lower bound $\Omega_{FE}I$ and the lower bound of
$\Gamma_{k}$ does not hold anymore. Thus exponential decay of the parameter
error cannot be ensured when $k>k_{3}^{\prime}$.

######  Remark 7.

Theorem 1 guarantees convergence of $V_{k}$ to a compact set $D_{PE}$, while
Theorem 2 guarantees that $V_{k}$ approaches $D_{FE}$. For static parameters,
i.e., $\theta_{k}^{*}\equiv\theta^{*}$ and $\Delta^{*}=0$, from (29) and (46)
it can be shown that under PE, $\tilde{\theta}_{k}$ converges exponentially
fast towards the origin; While under FE, the trajectories of
$\tilde{\theta}_{k}$ tend exponentially fast towards the origin.

######  Remark 8.

Note that under assumptions 1 and 2, (29) can also be written as

$V_{PE}=2\Delta^{*}\theta_{\max}^{*}\left[\frac{c_{1}+\sqrt{c_{1}^{2}+4c_{2}(\mu_{1}-\mu_{2})}}{2(\mu_{1}-\mu_{2})}\right]^{2},$
(50)

and (46) can also be written as

$V_{FE}=2\Delta^{*}\theta_{\max}^{*}\left[\frac{c_{1}^{\prime}+\sqrt{{c_{1}^{\prime}}^{2}+4c_{2}^{\prime}(\mu_{1}^{\prime}-\mu_{2}^{\prime})}}{2(\mu_{1}^{\prime}-\mu_{2}^{\prime})}\right]^{2}.$
(51)

From (50) and (51), note that the compact set scales linearly with both the
bounds of the unknown parameter $\theta_{\max}^{*}$ and the time variation of
the unknown parameter $\Delta^{*}$.

######  Remark 9.

A few comments about the choice of the hyperparameters in (7), (8) and (9) for
both cases of excitation. We assume $\Omega_{PE}$ and $\Omega_{FE}$ are known
and $\Gamma_{\max}$ is given. Parameter $\lambda_{\Omega}$ is chosen according
to (20) and (37), and $\Omega_{\max}$ is available according to Lemma 2. Then
$\lambda_{\Gamma}$ and $\kappa$ will in turn be chosen from the given bounds
in (21), (24), (38) and (41). Thus all the hyperparameters are now available
and $\Gamma_{PE}$ and $\Gamma_{FE}$ can be derived from (25) and (42).

##  6 CONCLUSIONS

We have presented a new adaptive algorithm for discrete-time parameter
estimation of a class of time-varying plants. We showed that in the presence
of time-varying parameters, the parameter estimation error converges uniformly
to a compact set under conditions of persistent excitation, with the size of
the compact set proportional to the time-variation of the unknown parameters.
Under conditions of finite excitation, the convergence is shown to be
asymptotic and non-uniform. The bounds for the compact set are shown to scale
linearly with both the bounds of the unknown parameter and the time variation
of the unknown parameter. All results derived assume conditions of excitation
on exogenous signals. These may be relaxed through the use of a projection
operator and are the subject of future work.

## References

  * [1] K. S. Narendra and A. M. Annaswamy, _Stable Adaptive Systems_. Dover, 2005.
  * [2] ——, _Stable Adaptive Systems_. NJ: Prentice-Hall, Inc., 1989, (out of print).
  * [3] E. Lavretsky and K. A. Wise, _Robust and Adaptive Control with Aerospace Applications_. Springer London, 2013.
  * [4] G. C. Goodwin and K. S. Sin, _Adaptive Filtering Prediction and Control_. Prentice Hall, 1984.
  * [5] P. A. Ioannou and J. Sun, _Robust Adaptive Control_. PTR Prentice-Hall, 1996.
  * [6] I. D. Landau, R. Lozano, M. M’Saad, and A. Karimi, _Adaptive Control: Algorithms, Analysis and Applications_. Springer Science & Business Media, 2011.
  * [7] L. Ljung, _System Identification: Theory for the User_. Prentice-Hall, 1987.
  * [8] B. D. Anderson and C. Johnson, “Exponential convergence of adaptive identification and control algorithms,” _Automatica_ , vol. 18, no. 1, pp. 1–13, 1982.
  * [9] R. Ortega, S. Aranovskiy, A. Pyrkin, A. Astolfi, and A. Bobtsov, “New results on parameter estimation via dynamic regressor extension and mixing: Continuous and discrete-time cases,” _IEEE Transactions on Automatic Control_ , 2020.
  * [10] S. Aranovskiy, A. Belov, R. Ortega, N. Barabanov, and A. Bobtsov, “Parameter identification of linear time-invariant systems using dynamic regressor extension and mixing,” _International Journal of Adaptive Control and Signal Processing_ , vol. 33, no. 6, pp. 1016–1030, 2019.
  * [11] S. Aranovskiy, A. Bobtsov, R. Ortega, and A. Pyrkin, “Parameters estimation via dynamic regressor extension and mixing,” in _2016 American Control Conference (ACC)_. IEEE, 2016.
  * [12] K. S. Narendra and Z. Han, “Discrete-time adaptive control using multiple models,” in _Proceedings of the 2011 American Control Conference_. IEEE, 2011, pp. 2921–2926.
  * [13] J. E. Gaudio, A. M. Annaswamy, E. Lavretsky, and M. A. Bolender, “Parameter estimation in adaptive control of time-varying systems under a range of excitation conditions,” _IEEE Transactions on Automatic Control (submitted)_ , 2020.
  * [14] J. E. Gaudio, “Fast learning and adaptation in control and machine learning,” Ph.D. dissertation, Massachusetts Institute of Technology, 2020.
  * [15] S. Boyd and S. Sastry, “On parameter convergence in adaptive control,” _Systems & Control Letters_, vol. 3, no. 6, pp. 311–319, 1983.
  * [16] S. Boyd and S. S. Sastry, “Necessary and sufficient conditions for parameter convergence in adaptive control,” _Automatica_ , vol. 22, no. 6, pp. 629–639, 1986.
  * [17] K. S. Narendra and A. M. Annaswamy, “Persistent excitation in adaptive systems,” _International Journal of Control_ , vol. 45, no. 1, pp. 127–160, 1987.
  * [18] A. L. Bruce, A. Goel, and D. S. Bernstein, “Recursive least squares with matrix forgetting,” in _2020 American Control Conference (ACC)_. IEEE, 2020, pp. 1406–1410.
  * [19] S. A. U. Islam, T. W. Nguyen, I. V. Kolmanovsky, and D. S. Bernstein, “Adaptive control of discrete-time systems with unknown, unstable zero dynamics,” in _2020 American Control Conference (ACC)_. IEEE, 2020, pp. 1387–1392.
  * [20] M. Yu and D. Huang, “Discrete-time adaptive learning control for parametric uncertainties with unknown periods,” in _52nd IEEE Conference on Decision and Control_. IEEE, 2013, pp. 1786–1791.
  * [21] T.-H. Lee and K. S. Narendra, “Robust adaptive control of discrete-time systems using persistent excitation,” _Automatica_ , vol. 24, no. 6, pp. 781–788, 1988.
  * [22] K. Petersen and M. Pedersen, “The matrix cookbook. technical university of denmark,” _Technical Manual_ , 2008.
  * [23] S. Axler, _Linear algebra done right_. Springer, 2015, vol. 2.